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