function [sigma_12r, sigma_12theta, sigma_3r, sigma_3theta sigma_incr sigma_inctheta] = sigma(n,r)
global kp_a kp_b ks_a ks_b lambda_a lambda_b mu_a mu_b

sigma_12r = [kp_b**2*((-lambda_b-2*mu_b/(kp_b*r)**2*(n+(kp_b*r)**2-n**2))*besselh(n,1,kp_b*r)+2*mu_b/(kp_b*r)*besselh(n+1,1,kp_b*r)), 2*mu_b*n/r*(-1/r*(1-n)*besselh(n,1,ks_b*r)-ks_b*besselh(n+1,1,ks_b*r)), kp_b**2*((-lambda_b+2*mu_b/(kp_b*r)**2*(n-(kp_b*r)**2+n**2))*besselh(n,2,kp_b*r)-2*mu_b/(kp_b*r)*besselh(n-1,2,kp_b*r)), 2*mu_b*n/r*(-1/r*(1+n)*besselh(n,2,ks_b*r)+ks_b*besselh(n-1,2,ks_b*r))];
sigma_12theta =[-2*mu_b*n/r*(-kp_b*besselh(n+1,1,kp_b*r)+(n-1)/r*besselh(n,1,kp_b*r)), 2*mu_b/r*(1/r*(-n**2+n+0.5*(ks_b*r)**2)*besselh(n,1,ks_b*r)-ks_b*besselh(n+1,1,ks_b*r)), -2*mu_b*n/r*(kp_b*besselh(n-1,2,kp_b*r)-(n+1)/r*besselh(n,2,kp_b*r)), 2*mu_b/r*(-1/r*(n**2+n-0.5*(ks_b*r)**2)*besselh(n,2,ks_b*r)+ks_b*besselh(n-1,2,ks_b*r))];
sigma_3r =[kp_a**2*((-lambda_a-2*mu_a/(kp_a*r)**2*(n+(kp_a*r)**2-n**2))*besselj(n,kp_a*r)+2*mu_a/(kp_a*r)*besselj(n+1,kp_a*r)), 2*mu_a*n/r*(-1/r*(1-n)*besselj(n,ks_a*r)-ks_a*besselj(n+1,ks_a*r))];
sigma_3theta = [-2*mu_a*n/r*(-kp_a*besselj(n+1,kp_a*r)+(n-1)/r*besselj(n,kp_a*r)), 2*mu_a/r*(1/r*(-n**2+n+0.5*(ks_a*r)**2)*besselj(n,ks_a*r)-ks_a*besselj(n+1,ks_a*r))];
sigma_incr =epsilon(n)*1i**n*((-lambda_b*kp_b**2-2*mu_b/r**2*(n-n**2+(kp_b*r)**2))*besselj(n,kp_b*r)+2*mu_b*kp_b/r*besselj(n+1,kp_b*r));
sigma_inctheta =-2*mu_b*epsilon(n)*1i**n*n/r*(-(1-n)/r*besselj(n,kp_b*r)-kp_b*besselj(n+1,kp_b*r));

endfunction


function epsilon = epsilon(n)
  if n==0
    epsilon = 1;
  else
    epsilon = 2;
  end
endfunction  