1function [A1n,A2n, B1n,B2n] = cylindrical_wall_s_out_abc2(n,kp,ks,eps_p,eps_s,a,b,lambda,mu,rho,omega) 2 3 kpeps = kp +i* eps_p; 4 kseps = ks +i* eps_s; 5 %Amplitude of incomming wave 6 psi_0 = 1; 7 xip = (1/kp) * ( 1/(sqrt(1-(n/(kpeps*b))^2 ) )); 8 xis = (1/ks) * ( 1/(sqrt(1-(n/(kseps*b))^2 ) )); 9 10 det = 1 + ((n/b)**2) * xip * xis ; 11 12 D1 = lambda+2*mu; 13 D2 = (rho*(omega**2)*xip * xis) / det; 14 D3 = lambda/b - (i*rho*(omega**2)*xip)/det; 15 D4 = -mu/b - (i*rho*(omega**2)*xis)/det; 16 17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%% 18[epsilon,dr_j_p,dr_j_s,dr_h1_p_a,dr_h2_p_a,dr_h1_s_a,dr_h2_s_a,dr_h1_p,dr_h2_p,dr_h1_s,dr_h2_s,dr2_h1_p,dr2_h2_p,dr2_h1_s,dr2_h2_s] =... 19 bessel_family(n,kp,ks,a,b); 20 21 M(1,1) = dr_h1_p_a; 22 M(1,2) = dr_h2_p_a; 23 M(1,3) = (-n/a) * besselh(n,1,ks*a); 24 M(1,4) = (-n/a) * besselh(n,2,ks*a); 25 26 M(2,1) = (n/a) * besselh(n,1,kp*a); 27 M(2,2) = (n/a) * besselh(n,2,kp*a); 28 M(2,3) = -dr_h1_s_a; 29 M(2,4) = -dr_h2_s_a; 30 31 M(3,1) = D1*dr2_h1_p + (D1-D2)*(-(n/b)**2)*besselh(n,1,kp*b) + D3*dr_h1_p; 32 M(3,2) = D1*dr2_h2_p + (D1-D2)*(-(n/b)**2)*besselh(n,2,kp*b) + D3*dr_h2_p; 33 M(3,3) = -D1*( -(n/(b**2)) *besselh(n,1,ks*b) + (n/b)*dr_h1_s ) - (D1-D2)*(-n/b)*dr_h1_s - D3*(n/b)*besselh(n,1,ks*b); 34 M(3,4) = -D1*( -(n/(b**2)) *besselh(n,2,ks*b) + (n/b)*dr_h2_s ) - (D1-D2)*(-n/b)*dr_h2_s - D3*(n/b)*besselh(n,2,ks*b); 35 36 M(4,1) = mu* (-(n/(b**2))*besselh(n,1,kp*b) + (n/b)*dr_h1_p) + (D2-mu)* (n/b)*dr_h1_p + D4*(n/b)*besselh(n,1,kp*b); 37 M(4,2) = mu* (-(n/(b**2))*besselh(n,2,kp*b) + (n/b)*dr_h2_p) + (D2-mu)* (n/b)*dr_h2_p + D4*(n/b)*besselh(n,2,kp*b); 38 M(4,3) = -mu*dr2_h1_s + (D2-mu)*(-(n/b)**2)*besselh(n,1,ks*b) - D4*dr_h1_s; 39 M(4,4) = -mu*dr2_h2_s + (D2-mu)*(-(n/b)**2)*besselh(n,2,ks*b) - D4*dr_h2_s; 40 41 F(1) = psi_0*epsilon*((-1i)**n)*(n/a)*besselj(n,ks*a); 42 F(2) = psi_0*epsilon*((-1i)**n)*dr_j_s; 43 F(3) = 0; 44 F(4) = 0; 45 46 % Calculate the inverse determinant of the matrix 47 detinv = 1/(M(1,1)*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))... 48 +M(2,4)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))) - M(1,2)*(M(2,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)... 49 *(M(3,4)*M(4,1)-M(3,1)*M(4,4))+M(2,4)*(M(3,1)*M(4,3)-M(3,3)*M(4,1)))+... 50 M(1,3)*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))+M(2,2)*(M(3,4)*M(4,1)-M(3,1)*M(4,4))+M(2,4)*(M(3,1)*M(4,2)-M(3,2)*M(4,1)))... 51 - M(1,4)*(M(2,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))+M(2,2)*(M(3,3)*M(4,1)-M(3,1)*M(4,3))+M(2,3)*(M(3,1)*M(4,2)-M(3,2)*M(4,1)))); 52 53 % Calculate the inverse of the matrix B=M^{-1} 54 Minv(1,1) = detinv*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))... 55 +M(2,3)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))+M(2,4)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))); 56 Minv(2,1) = detinv*(M(2,1)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))... 57 +M(2,3)*(M(3,1)*M(4,4)-M(3,4)*M(4,1))+M(2,4)*(M(3,3)*M(4,1)-M(3,1)*M(4,3))); 58 Minv(3,1) = detinv*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))... 59 +M(2,2)*(M(3,4)*M(4,1)-M(3,1)*M(4,4))+M(2,4)*(M(3,1)*M(4,2)-M(3,2)*M(4,1))); 60 Minv(4,1) = detinv*(M(2,1)*(M(3,3)*M(4,2)-M(3,2)*M(4,3))... 61 +M(2,2)*(M(3,1)*M(4,3)-M(3,3)*M(4,1))+M(2,3)*(M(3,2)*M(4,1)-M(3,1)*M(4,2))); 62 Minv(1,2) = detinv*(M(1,2)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))... 63 +M(1,3)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))+M(1,4)*(M(3,3)*M(4,2)-M(3,2)*M(4,3))); 64 Minv(2,2) = detinv*(M(1,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))... 65 +M(1,3)*(M(3,4)*M(4,1)-M(3,1)*M(4,4))+M(1,4)*(M(3,1)*M(4,3)-M(3,3)*M(4,1))); 66 Minv(3,2) = detinv*(M(1,1)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))... 67 +M(1,2)*(M(3,1)*M(4,4)-M(3,4)*M(4,1))+M(1,4)*(M(3,2)*M(4,1)-M(3,1)*M(4,2))); 68 Minv(4,2) = detinv*(M(1,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))... 69 +M(1,2)*(M(3,3)*M(4,1)-M(3,1)*M(4,3))+M(1,3)*(M(3,1)*M(4,2)-M(3,2)*M(4,1))); 70 Minv(1,3) = detinv*(M(1,2)*(M(2,3)*M(4,4)-M(2,4)*M(4,3))... 71 +M(1,3)*(M(2,4)*M(4,2)-M(2,2)*M(4,4))+M(1,4)*(M(2,2)*M(4,3)-M(2,3)*M(4,2))); 72 Minv(2,3) = detinv*(M(1,1)*(M(2,4)*M(4,3)-M(2,3)*M(4,4))... 73 +M(1,3)*(M(2,1)*M(4,4)-M(2,4)*M(4,1))+M(1,4)*(M(2,3)*M(4,1)-M(2,1)*M(4,3))); 74 Minv(3,3) = detinv*(M(1,1)*(M(2,2)*M(4,4)-M(2,4)*M(4,2))... 75 +M(1,2)*(M(2,4)*M(4,1)-M(2,1)*M(4,4))+M(1,4)*(M(2,1)*M(4,2)-M(2,2)*M(4,1))); 76 Minv(4,3) = detinv*(M(1,1)*(M(2,3)*M(4,2)-M(2,2)*M(4,3))... 77 +M(1,2)*(M(2,1)*M(4,3)-M(2,3)*M(4,1))+M(1,3)*(M(2,2)*M(4,1)-M(2,1)*M(4,2))); 78 Minv(1,4) = detinv*(M(1,2)*(M(2,4)*M(3,3)-M(2,3)*M(3,4))... 79 +M(1,3)*(M(2,2)*M(3,4)-M(2,4)*M(3,2))+M(1,4)*(M(2,3)*M(3,2)-M(2,2)*M(3,3))); 80 Minv(2,4) = detinv*(M(1,1)*(M(2,3)*M(3,4)-M(2,4)*M(3,3))... 81 +M(1,3)*(M(2,4)*M(3,1)-M(2,1)*M(3,4))+M(1,4)*(M(2,1)*M(3,3)-M(2,3)*M(3,1))); 82 Minv(3,4) = detinv*(M(1,1)*(M(2,4)*M(3,2)-M(2,2)*M(3,4))... 83 +M(1,2)*(M(2,1)*M(3,4)-M(2,4)*M(3,1))+M(1,4)*(M(2,2)*M(3,1)-M(2,1)*M(3,2))); 84 Minv(4,4) = detinv*(M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))... 85 +M(1,2)*(M(2,3)*M(3,1)-M(2,1)*M(3,3))+M(1,3)*(M(2,1)*M(3,2)-M(2,2)*M(3,1))); 86 87 % Compute AB = M^{-1}F 88 A1n = Minv(1,1)*F(1) + Minv(1,2)*F(2) + Minv(1,3)*F(3) + Minv(1,4)*F(4); 89 A2n = Minv(2,1)*F(1) + Minv(2,2)*F(2) + Minv(2,3)*F(3) + Minv(2,4)*F(4); 90 B1n = Minv(3,1)*F(1) + Minv(3,2)*F(2) + Minv(3,3)*F(3) + Minv(3,4)*F(4); 91 B2n = Minv(4,1)*F(1) + Minv(4,2)*F(2) + Minv(4,3)*F(3) + Minv(4,4)*F(4); 92 93