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