1function  [A1n,A2n, B1n,B2n] = cylindrical_wall_s_out_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  psi_0 = 1;
7  D1= lambda+2*mu;
8  [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] =...
9 bessel_family(n,kp,ks,a,b);
10  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%%
11  [Aj,Bj] = ComplexPadeISqrt2(L,theta);
12  zp=-(n/(kpeps*b))^2;
13  zs=-(n/(kseps*b))^2;
14  xip = 0;
15  xis = 0;
16  for ind=0:L-1
17    xip = xip + Aj(ind+1)./(Bj(ind+1)+zp);
18    xis = xis + Aj(ind+1)./(Bj(ind+1)+zs);
19  end
20  xip=xip/kp;
21  xis=xis/ks;
22
23  det = 1.d0 + ((n/b)**2) * xip * xis;
24
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    %    disp(['n:', num2str(n),'M:', num2str(M(1,1))])
50
51  F(1) = psi_0*epsilon*((-1i)**n)*(n/a)*besselj(n,ks*a);
52  F(2) = psi_0*epsilon*((-1i)**n)*dr_j_s;
53  F(3) = 0;
54  F(4) = 0;
55
56    % Calculate the inverse determinant of the matrix
57    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))...
58    +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)...
59      *(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)))+...
60      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)))...
61      - 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))));
62
63    % Calculate the inverse of the matrix B=M^{-1}
64    Minv(1,1) = detinv*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))...
65     +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)));
66    Minv(2,1) = detinv*(M(2,1)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))...
67     +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)));
68    Minv(3,1) = detinv*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))...
69     +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)));
70    Minv(4,1) = detinv*(M(2,1)*(M(3,3)*M(4,2)-M(3,2)*M(4,3))...
71     +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)));
72    Minv(1,2) = detinv*(M(1,2)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))...
73     +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)));
74    Minv(2,2) = detinv*(M(1,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))...
75     +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)));
76    Minv(3,2) = detinv*(M(1,1)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))...
77     +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)));
78    Minv(4,2) = detinv*(M(1,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))...
79     +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)));
80    Minv(1,3) = detinv*(M(1,2)*(M(2,3)*M(4,4)-M(2,4)*M(4,3))...
81     +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)));
82    Minv(2,3) = detinv*(M(1,1)*(M(2,4)*M(4,3)-M(2,3)*M(4,4))...
83     +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)));
84    Minv(3,3) = detinv*(M(1,1)*(M(2,2)*M(4,4)-M(2,4)*M(4,2))...
85     +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)));
86    Minv(4,3) = detinv*(M(1,1)*(M(2,3)*M(4,2)-M(2,2)*M(4,3))...
87     +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)));
88    Minv(1,4) = detinv*(M(1,2)*(M(2,4)*M(3,3)-M(2,3)*M(3,4))...
89     +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)));
90    Minv(2,4) = detinv*(M(1,1)*(M(2,3)*M(3,4)-M(2,4)*M(3,3))...
91     +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)));
92    Minv(3,4) = detinv*(M(1,1)*(M(2,4)*M(3,2)-M(2,2)*M(3,4))...
93     +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)));
94    Minv(4,4) = detinv*(M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))...
95     +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)));
96
97    % Compute AB = M^{-1}F
98    A1n = Minv(1,1)*F(1) + Minv(1,2)*F(2) + Minv(1,3)*F(3) + Minv(1,4)*F(4);
99    A2n = Minv(2,1)*F(1) + Minv(2,2)*F(2) + Minv(2,3)*F(3) + Minv(2,4)*F(4);
100    B1n = Minv(3,1)*F(1) + Minv(3,2)*F(2) + Minv(3,3)*F(3) + Minv(3,4)*F(4);
101    B2n = Minv(4,1)*F(1) + Minv(4,2)*F(2) + Minv(4,3)*F(3) + Minv(4,4)*F(4);
102