1function  [A1n,A2n, B1n,B2n] = cylindrical_wallout_abc(n,kp,ks,a,b,lambda,mu,rho,omega)
2
3  %Amplitude of incomming wave
4  phi_0 = 1;
5  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%%
6  C1 = lambda+2*mu;
7  C2 = ((lambda/b) - i*kp*(lambda+2*mu));
8  C3 = ((mu/b) + i*ks*mu);
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) = -phi_0*epsilon * ((-1i)**n) * dr_j_p;
34  F(2) = (n/a) * phi_0 * epsilon* ((-1i)**n) *besselj(n,kp*a);
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  %end
85
86end
87