"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/contrib/pewe/fortran/cylindrical_walloutabc2.f90" (31 Jul 2018, 18125 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Fortran 90 source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "cylindrical_walloutabc2.f90" see the Fossies "Dox" file reference documentation.

    1 ! cylindrical_cavity This program computes particular solutions to
    2 ! the elastic wave equation in cylindrical geometries,
    3 ! see: https://bitbucket.org/appelo/pewe
    4 !
    5 ! Copyright (C) 2015 Kristoffer Virta & Daniel Appelo
    6 !
    7 !
    8 !  This program is free software: you can redistribute it and/or modify
    9 !  it under the terms of the GNU General Public License as published by
   10 !  the Free Software Foundation, either version 3 of the License, or
   11 !  (at your option) any later version.
   12 !
   13 !  This program is distributed in the hope that it will be useful,
   14 !  but WITHOUT ANY WARRANTY; without even the implied warranty of
   15 !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   16 !  GNU General Public License for more details.
   17 !
   18 !  You should have received a copy of the GNU General Public License
   19 !  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   20 !
   21 
   22 
   23 !  Modified by Vanessa Mattesi for inclusion in GetDP:
   24 !
   25 !  Exact solution for a cylindrical wall (zero displacement on the
   26 !  boundary).
   27 
   28 subroutine cylindrical_walloutabc2(du,dv,dut,dvt,X,Y,t,omega,lambda,mu,rho,a,b)
   29 
   30   implicit none
   31   double precision :: x,y,du,dv,dvt,dut,t
   32   double precision :: r,theta
   33   double precision :: lambda,mu,rho
   34   double precision :: cp,cs,omega,kp,ks,a,b
   35   double precision :: phi_0,epsilon_0,epsilon_1
   36 
   37   double precision :: a0_n,a1_n
   38   double complex M(4,4)
   39   double complex Minv(4,4)
   40   double complex AB(4)
   41   double complex F(4)
   42   double complex   :: i,p,q,detinv,kpeps,kseps,xip,xis,det
   43 
   44   integer :: n
   45   double precision , parameter :: pi = acos(-1.d0)
   46 
   47   ! for GetDP
   48   integer :: ns
   49  
   50   double complex, external :: besselh
   51   t=0
   52   i = (0.d0,1.d0)
   53   ! Compute radius r and angle theta
   54   r = sqrt(X**2+Y**2)
   55   theta = atan2(Y,X)
   56 
   57   ! P and S wave speeds
   58   cp = sqrt((lambda+2.d0*mu)/rho)
   59   cs = sqrt(mu/rho)
   60   ! To satisfy elastic wave equation
   61   kp = omega/cp
   62   ks = omega/cs
   63 
   64   kpeps = kp + i* 0.4*(1.0d0/b)**(2.0d0/3) * kp**(1.0d0/3);
   65   kseps = ks + i* 0.4*(1.0d0/b)**(2.0d0/3) * ks**(1.0d0/3);
   66 
   67   ! Amplitude of incomming wave
   68   phi_0 = 1
   69   ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Series expansion of solution%%%%%%%%%%%
   70   n = 0
   71   epsilon_0 = 1
   72 
   73   F(1) = phi_0*epsilon_0*kp*besjn(1,kp*a)
   74   F(2) = 0
   75 
   76   M(1,1) = -kp*besselh(1,1,kp*a)
   77   M(1,2) = -kp*besselh(1,2,kp*a)
   78   M(2,1) = -(lambda+2*mu) * ((kp**2)/2.0d0) * (besselh(0,1,kp*b) - besselh(2,1,kp*b) )&
   79          - ((lambda/dble(b)) - i*kp*(lambda+2*mu))*kp*besselh(1,1,kp*b)
   80   M(2,2) = -(lambda+2*mu) * ((kp**2)/2.0d0) * (besselh(0,2,kp*b) - besselh(2,2,kp*b) )&
   81          - ( (lambda/dble(b)) - i*kp*(lambda+2*mu))*kp*besselh(1,2,kp*b) 
   82 
   83   AB(1) = 1.0d0/(M(1,1)*M(2,2) - M(1,2)*M(2,1))*( M(2,2)*F(1)-M(1,2)*F(2))
   84   AB(2) = 1.0d0/(M(1,1)*M(2,2) - M(1,2)*M(2,1))*(-M(2,1)*F(1)+M(1,1)*F(2))
   85 
   86   p = -kp*AB(1)*besselh(1,1,kp*r) - kp*AB(2)*besselh(1,2,kp*r)
   87 
   88   q=0
   89   epsilon_1 = 2.0d0
   90   ! for GetDP
   91   ns = 55;
   92 
   93   do n = 1,1
   94 
   95     a0_n = 0
   96     a1_n = 3
   97 
   98     xip = (1.0d0/kp) * ( 1.0d0/(sqrt(1-(dble(n)**2/((kpeps**2)* b**2 )) ) ));
   99     xis = (1.0d0/ks) * ( 1.0d0/(sqrt(1-(dble(n)**2/((kseps**2)* b**2 )) ) ));
  100 
  101     det = 1 + ((dble(n)**2)/(b**2)) * xip * xis ;
  102 
  103     M(1,1) = (kp/2.d0) * ( besselh(n-1,1,kp*a)-besselh(n+1,1,kp*a) )
  104     M(1,2) = (kp/2.d0) * ( besselh(n-1,2,kp*a)-besselh(n+1,2,kp*a) )
  105     M(1,3) = ( dble(n)/a) * besselh(int(n),1,ks*a)
  106     M(1,4) = ( dble(n)/a) * besselh(int(n),2,ks*a)
  107 
  108     M(2,1) = (-dble(n)/a)  * besselh(int(n),1,kp*a)
  109     M(2,2) = (-dble(n)/a)  * besselh(int(n),2,kp*a)
  110     M(2,3) = (-ks/2.d0) * ( besselh(n-1,1,ks*a)-besselh(n+1,1,ks*a) )
  111     M(2,4) = (-ks/2.d0) * ( besselh(n-1,2,ks*a)-besselh(n+1,2,ks*a) )
  112 
  113     M(3,1) = (lambda + 2*mu) * ((kp**2)/4.d0) * (- a1_n*besselh(int(n),1,kp*b)+besselh(n+2,1,kp*b))&
  114            + (lambda-((rho*xip*xis*omega**2)/ det) + 2*mu)&
  115            *(-dble(n)**2)/(b**2) * besselh(int(n),1,kp*b)&
  116            + (lambda/dble(b))-( (i*rho*xip*omega**2)/det )&
  117            *(kp/2.d0)*(besselh(n-1,1,kp*b)-besselh(n+1,1,kp*b))
  118     M(3,2) = (lambda + 2*mu) * ((kp**2)/4.d0) * (- a1_n*besselh(int(n),2,kp*b)+besselh(n+2,2,kp*b))&
  119            + (lambda-((rho*xip*xis*omega**2)/(det)) + 2*mu)&
  120            *(-dble(n)**2)/(b**2) * besselh(int(n),2,kp*b)&
  121            + (lambda/dble(b))-( (i*rho*xip*omega**2)/det )&
  122            *(kp/2.d0)*(besselh(n-1,2,kp*b)-besselh(n+1,2,kp*b))
  123     M(3,3) = (rho*xip*xis*omega**2/det)*(dble(n)*ks/(2*b))&
  124              *( besselh(n-1, 1 ,ks*b) - besselh(n+1,1,ks*b) )&
  125            + (-2*mu/dble(b) -(i*rho*xip*omega**2)/(det) )*(dble(n)/b)*besselh(int(n),1,ks*b)
  126 
  127     M(3,4) = (rho*xip*xis*omega**2/det)*(dble(n)*ks/(2*b))*(besselh(n-1,2,ks*b)-besselh(n+1,2,ks*b))&
  128             + (-2*mu/dble(b) -(i*rho*xip*omega**2)/(det) )*(dble(n)/b)*besselh(int(n),2,ks*b)
  129 
  130 
  131 
  132     M(4,1) = (-(rho*xip*xis*omega**2)/det ) * (dble(n)*kp)/(2*b) * (besselh(n-1,1,kp*b) - besselh(n+1,1,kp*b))&
  133            + ((2*mu/dble(b)) + (i*rho*xis*omega**2)/det ) * (dble(n)/b) * besselh(int(n),1,kp*b)
  134 
  135     M(4,2) = (-(rho*xip*xis*omega**2)/det ) * (dble(n)*kp)/(2*b) * (besselh(n-1,2,kp*b) - besselh(n+1,2,kp*b))&
  136            + ((2*mu/dble(b)) + (i*rho*xis*omega**2)/det ) * (dble(n)/b) * besselh(int(n),2,kp*b)
  137 
  138     M(4,3) = -mu*(ks**2/4.0d0) * (-a1_n*besselh(int(n),1,ks*b)+besselh(n+2,1,ks*b))&
  139            + (mu - (rho*xip*xis*omega**2)/det) * (dble(n)**2/b**2) * besselh(int(n),1,ks*b)&
  140            + ((mu/dble(b)) + (i*rho*xis*omega**2)/det) * (ks/2.d0) * (besselh(n-1,1,ks*b)-besselh(n+1,1,ks*b))
  141 
  142     M(4,4) = -mu*(ks**2/4.0d0) * (-a1_n*besselh(int(n),2,ks*b)+besselh(n+2,2,ks*b))&
  143            + (mu - (rho*xip*xis*omega**2)/det) * (dble(n)**2/b**2) * besselh(int(n),2,ks*b)&
  144            + ((mu/dble(b)) + (i*rho*xis*omega**2)/det) * (ks/2.d0) * (besselh(n-1,2,ks*b)-besselh(n+1,2,ks*b))
  145 
  146     F(1) = -phi_0*epsilon_1* (0.d0,-1.d0)**n * (kp/2.d0) * (besjn(n-1,kp*a)-besjn(n+1,kp*a))
  147     F(2) = (dble(n)/a) * phi_0 * epsilon_1* (0.d0,-1.d0)**n *besjn(n,kp*a)
  148 
  149     F(3) = 0
  150     F(4) = 0
  151     !F(3) = (phi_0*epsilon_1* (-i)**n) * (-(lambda+2*mu) * (kp**2/4.d0) * (-a1_n*besjn(n,ks*b)+besjn(n+2,ks*b))&
  152     !     + (lambda*dble(n)**2/b**2) * besjn(n,kp*b)&
  153     !     - ((lambda/b)-i*kp*(lambda+2*mu))* (kp/2.d0) * (besjn(n-1,kp*b)-besjn(n+1,kp*b)) )
  154     !F(4) = (phi_0*epsilon_1* (-i)**n) * ((mu*dble(n)/b) * kp * (besjn(n-1,kp*b)-besjn(n+1,kp*b))&
  155     !     - ((mu/b)-i*ks*mu)* (dble(n)/b)* besjn(n,kp*b) )
  156 
  157     ! Calculate the inverse determinant of the matrix
  158     detinv = &
  159       1/(M(1,1)*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)&
  160       *(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)))&
  161        - M(1,2)*(M(2,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)&
  162       *(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)))&
  163        + M(1,3)*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))+M(2,2)&
  164       *(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)))&
  165        - M(1,4)*(M(2,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))+M(2,2)&
  166       *(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))))
  167 
  168     ! Calculate the inverse of the matrix B=M^{-1}
  169     Minv(1,1) = detinv*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))&
  170      +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)))
  171     Minv(2,1) = detinv*(M(2,1)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))&
  172      +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)))
  173     Minv(3,1) = detinv*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))&
  174      +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)))
  175     Minv(4,1) = detinv*(M(2,1)*(M(3,3)*M(4,2)-M(3,2)*M(4,3))&
  176      +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)))
  177     Minv(1,2) = detinv*(M(1,2)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))&
  178      +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)))
  179     Minv(2,2) = detinv*(M(1,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))&
  180      +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)))
  181     Minv(3,2) = detinv*(M(1,1)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))&
  182      +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)))
  183     Minv(4,2) = detinv*(M(1,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))&
  184      +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)))
  185     
  186     Minv(1,3) = detinv*(M(1,2)*(M(2,3)*M(4,4)-M(2,4)*M(4,3))&
  187      +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)))
  188     Minv(2,3) = detinv*(M(1,1)*(M(2,4)*M(4,3)-M(2,3)*M(4,4))&
  189      +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)))
  190     Minv(3,3) = detinv*(M(1,1)*(M(2,2)*M(4,4)-M(2,4)*M(4,2))&
  191     +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)))
  192     Minv(4,3) = detinv*(M(1,1)*(M(2,3)*M(4,2)-M(2,2)*M(4,3))&
  193      +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)))
  194     Minv(1,4) = detinv*(M(1,2)*(M(2,4)*M(3,3)-M(2,3)*M(3,4))&
  195      +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)))
  196     Minv(2,4) = detinv*(M(1,1)*(M(2,3)*M(3,4)-M(2,4)*M(3,3))&
  197      +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)))
  198     Minv(3,4) = detinv*(M(1,1)*(M(2,4)*M(3,2)-M(2,2)*M(3,4))&
  199      +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)))
  200     Minv(4,4) = detinv*(M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))&
  201      +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)))
  202 
  203     ! Compute AB = M^{-1}F
  204     AB(1) = Minv(1,1)*F(1) + Minv(1,2)*F(2) + Minv(1,3)*F(3) + Minv(1,4)*F(4)
  205     AB(2) = Minv(2,1)*F(1) + Minv(2,2)*F(2) + Minv(2,3)*F(3) + Minv(2,4)*F(4)
  206     AB(3) = Minv(3,1)*F(1) + Minv(3,2)*F(2) + Minv(3,3)*F(3) + Minv(3,4)*F(4)
  207     AB(4) = Minv(4,1)*F(1) + Minv(4,2)*F(2) + Minv(4,3)*F(3) + Minv(4,4)*F(4)  
  208 
  209     ! Compute p and q
  210     p = p + ( AB(1)*(kp/2.d0)*(besselh(n-1,1,kp*r)-besselh(n+1,1,kp*r))&
  211       + AB(2)*(kp/2.d0)*(besselh(n-1,2,kp*r)-besselh(n+1,2,kp*r))&
  212       + (dble(n)/r)*AB(3)*besselh(int(n),1,ks*r) + (dble(n)/r)*AB(4)*besselh(int(n),2,ks*r) )*cos(dble(n)*theta)
  213 
  214     q = q + ( - (dble(n)/r)*AB(1)*besselh(int(n),1,kp*r) &
  215        - (dble(n)/r)*AB(2)*besselh(int(n),2,kp*r) - AB(3)*(ks/2.d0)*(besselh(n-1,1,ks*r)&
  216        -besselh(n+1,1,ks*r)) - AB(4)*(ks/2.d0)*(besselh(n-1,2,ks*r)-besselh(n+1,2,ks*r)) )*sin(dble(n)*theta)
  217   end do
  218   a0_n = 1
  219   a1_n = 2
  220   ! for GetDP: instead of 2:24
  221 !$OMP PARALLEL DO PRIVATE(n,xip,xis,det,detinv,M,Minv,F,AB) REDUCTION(+:p,q)
  222  do n = 2,ns
  223     xip = (1.0d0/kp) * ( 1.0d0/(sqrt(1-(dble(n)**2/((kpeps**2)* b**2 )) ) ));
  224     xis = (1.0d0/ks) * ( 1.0d0/(sqrt(1-(dble(n)**2/((kseps**2)* b**2 )) ) ));
  225 
  226     det = 1 + ((dble(n)**2)/(b**2)) * xip * xis ;
  227 
  228     M(1,1) = (kp/2.d0) * ( besselh(n-1,1,kp*a)-besselh(n+1,1,kp*a) )
  229     M(1,2) = (kp/2.d0) * ( besselh(n-1,2,kp*a)-besselh(n+1,2,kp*a) )
  230     M(1,3) = ( dble(n)/a) * besselh(int(n),1,ks*a)
  231     M(1,4) = ( dble(n)/a) * besselh(int(n),2,ks*a)
  232 
  233     M(2,1) = (-dble(n)/a)  * besselh(int(n),1,kp*a)
  234     M(2,2) = (-dble(n)/a)  * besselh(int(n),2,kp*a)
  235     M(2,3) = (-ks/2.d0) * ( besselh(n-1,1,ks*a)-besselh(n+1,1,ks*a) )
  236     M(2,4) = (-ks/2.d0) * ( besselh(n-1,2,ks*a)-besselh(n+1,2,ks*a) )
  237 
  238     M(3,1) = (lambda + 2*mu) * ((kp**2)/4.d0) * (a0_n*besselh(n-2,1,kp*b)&
  239            - a1_n*besselh(int(n),1,kp*b)+besselh(n+2,1,kp*b))&
  240            + (lambda-((rho*xip*xis*omega**2)/(det) + 2*mu)&
  241            *(-dble(n)**2)/(b**2) * besselh(int(n),1,kp*b)&
  242            + (lambda/dble(b))-( (i*rho*xip*omega**2)/ det) )&
  243            *(kp/2.d0)*(besselh(n-1,1,kp*b)-besselh(n+1,1,kp*b))
  244     M(3,2) = (lambda + 2*mu) * ((kp**2)/4.d0) * (a0_n*besselh(n-2,2,kp*b)&
  245            - a1_n*besselh(int(n),2,kp*b)+besselh(n+2,2,kp*b))&
  246            + (lambda-((rho*xip*xis*omega**2)/(det)) + 2*mu)&
  247            *(-dble(n)**2)/(b**2) * besselh(int(n),2,kp*b)&
  248            + ((lambda/dble(b))-( (i*rho*xip*omega**2)/(det) ) )&
  249            *(kp/2.d0)*(besselh(n-1,2,kp*b)-besselh(n+1,2,kp*b))
  250 
  251     M(3,3) = (rho*xip*xis*omega**2/det)*(dble(n)*ks/(2*b))*(besselh(n-1,1,ks*b)-besselh(n+1,1,ks*b))&
  252            + (-2*mu/dble(b) -(i*rho*xip*omega**2)/(det) )*(dble(n)/b)*besselh(int(n),1,ks*b)
  253 
  254     M(3,4) = (rho*xip*xis*omega**2/det)*(dble(n)*ks/(2*b))*(besselh(n-1,2,ks*b)-besselh(n+1,2,ks*b))&
  255            + (-2*mu/dble(b) -(i*rho*xip*omega**2)/(det) )*(dble(n)/b)*besselh(int(n),2,ks*b)
  256 
  257     M(4,1) = (-(rho*xip*xis*omega**2)/det ) * (dble(n)*kp)/(2*b) * (besselh(n-1,1,kp*b) - besselh(n+1,1,kp*b))&
  258            + ((2*mu/dble(b)) + (i*rho*xis*omega**2)/det ) * (dble(n)/b) * besselh(int(n),1,kp*b)
  259 
  260     M(4,2) = (-(rho*xip*xis*omega**2)/det ) * (dble(n)*kp)/(2*b) * (besselh(n-1,2,kp*b) - besselh(n+1,2,kp*b))&
  261            + ((2*mu/dble(b)) + (i*rho*xis*omega**2)/det ) * (dble(n)/b) * besselh(int(n),2,kp*b)
  262 
  263     M(4,3) = -mu*(ks**2/4.0d0) * (a0_n*besselh(n-2,1,ks*b) -a1_n*besselh(int(n),1,ks*b)+besselh(n+2,1,ks*b))&
  264            + (mu - (rho*xip*xis*omega**2)/det) * (dble(n)**2/b**2) * besselh(int(n),1,ks*b)&
  265            + ((mu/dble(b)) + (i*rho*xis*omega**2)/det) * (ks/2.d0) * (besselh(n-1,1,ks*b)-besselh(n+1,1,ks*b))
  266 
  267     M(4,4) = -mu*(ks**2/4.0d0) * (a0_n*besselh(n-2,2,ks*b) -a1_n*besselh(int(n),2,ks*b)+besselh(n+2,2,ks*b))&
  268            + (mu - (rho*xip*xis*omega**2)/det) * (dble(n)**2/b**2) * besselh(int(n),2,ks*b)&
  269            + ((mu/dble(b)) + (i*rho*xis*omega**2)/det) * (ks/2.d0) * (besselh(n-1,2,ks*b)-besselh(n+1,2,ks*b))
  270 
  271     F(1) = -phi_0*epsilon_1* (0.d0,-1.d0)**n * (kp/2.d0) * (besjn(n-1,kp*a)-besjn(n+1,kp*a))
  272     F(2) = (dble(n)/a) * phi_0 * epsilon_1* (0.d0,-1.d0)**n *besjn(n,kp*a)
  273 
  274     F(3) = 0
  275     F(4) = 0
  276     !F(3) = (phi_0*epsilon_1* (-i)**n) * (-(lambda+2*mu) * (kp**2/4.d0) * (a0_n*besjn(n-2,ks*b)-a1_n*besjn(n,ks*b)+besjn(n+2,ks*b))&
  277     !     + (lambda*dble(n)**2/b**2) * besjn(n,kp*b)&
  278     !     - ((lambda/b)-i*kp*(lambda+2*mu))* (kp/2.d0) * (besjn(n-1,kp*b)-besjn(n+1,kp*b)) )
  279     !F(4) = (phi_0*epsilon_1* (-i)**n) * ((mu*dble(n)/b) * kp * (besjn(n-1,kp*b)-besjn(n+1,kp*b))&
  280     !     - ((mu/b)-i*ks*mu)* (dble(n)/b)* besjn(n,kp*b) )
  281 
  282     ! Calculate the inverse determinant of the matrix
  283     detinv = &
  284       1/(M(1,1)*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)&
  285       *(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)))&
  286        - M(1,2)*(M(2,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))+M(2,3)&
  287       *(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)))&
  288        + M(1,3)*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))+M(2,2)&
  289       *(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)))&
  290        - M(1,4)*(M(2,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))+M(2,2)&
  291       *(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))))
  292 
  293     ! Calculate the inverse of the matrix B=M^{-1}
  294     Minv(1,1) = detinv*(M(2,2)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))&
  295      +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)))
  296     Minv(2,1) = detinv*(M(2,1)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))&
  297      +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)))
  298     Minv(3,1) = detinv*(M(2,1)*(M(3,2)*M(4,4)-M(3,4)*M(4,2))&
  299      +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)))
  300     Minv(4,1) = detinv*(M(2,1)*(M(3,3)*M(4,2)-M(3,2)*M(4,3))&
  301      +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)))
  302     Minv(1,2) = detinv*(M(1,2)*(M(3,4)*M(4,3)-M(3,3)*M(4,4))&
  303      +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)))
  304     Minv(2,2) = detinv*(M(1,1)*(M(3,3)*M(4,4)-M(3,4)*M(4,3))&
  305      +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)))
  306     Minv(3,2) = detinv*(M(1,1)*(M(3,4)*M(4,2)-M(3,2)*M(4,4))&
  307      +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)))
  308     Minv(4,2) = detinv*(M(1,1)*(M(3,2)*M(4,3)-M(3,3)*M(4,2))&
  309      +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)))
  310     
  311     Minv(1,3) = detinv*(M(1,2)*(M(2,3)*M(4,4)-M(2,4)*M(4,3))&
  312      +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)))
  313     Minv(2,3) = detinv*(M(1,1)*(M(2,4)*M(4,3)-M(2,3)*M(4,4))&
  314      +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)))
  315     Minv(3,3) = detinv*(M(1,1)*(M(2,2)*M(4,4)-M(2,4)*M(4,2))&
  316      +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)))
  317     Minv(4,3) = detinv*(M(1,1)*(M(2,3)*M(4,2)-M(2,2)*M(4,3))&
  318      +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)))
  319     Minv(1,4) = detinv*(M(1,2)*(M(2,4)*M(3,3)-M(2,3)*M(3,4))&
  320      +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)))
  321     Minv(2,4) = detinv*(M(1,1)*(M(2,3)*M(3,4)-M(2,4)*M(3,3))&
  322      +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)))
  323     Minv(3,4) = detinv*(M(1,1)*(M(2,4)*M(3,2)-M(2,2)*M(3,4))&
  324      +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)))
  325     Minv(4,4) = detinv*(M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))&
  326      +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)))
  327 
  328     ! Compute AB = M^{-1}F
  329     AB(1) = Minv(1,1)*F(1) + Minv(1,2)*F(2) + Minv(1,3)*F(3) + Minv(1,4)*F(4)
  330     AB(2) = Minv(2,1)*F(1) + Minv(2,2)*F(2) + Minv(2,3)*F(3) + Minv(2,4)*F(4)
  331     AB(3) = Minv(3,1)*F(1) + Minv(3,2)*F(2) + Minv(3,3)*F(3) + Minv(3,4)*F(4)
  332     AB(4) = Minv(4,1)*F(1) + Minv(4,2)*F(2) + Minv(4,3)*F(3) + Minv(4,4)*F(4)  
  333 
  334     ! Compute p and q
  335     p = p + (AB(1)*(kp/2.d0)*(besselh(n-1,1,kp*r)-besselh(n+1,1,kp*r))&
  336       + AB(2)*(kp/2.d0)*(besselh(n-1,2,kp*r)-besselh(n+1,2,kp*r))&
  337       + (dble(n)/r)*AB(3)*besselh(int(n),1,ks*r) + (dble(n)/r)*AB(4)*besselh(int(n),2,ks*r) )*cos(dble(n)*theta)
  338 
  339     q = q + (- (dble(n)/r)*AB(1)*besselh(int(n),1,kp*r) &
  340        - (dble(n)/r)*AB(2)*besselh(int(n),2,kp*r) - AB(3)*(ks/2.d0)*(besselh(n-1,1,ks*r)&
  341        -besselh(n+1,1,ks*r)) - AB(4)*(ks/2.d0)*(besselh(n-1,2,ks*r)-besselh(n+1,2,ks*r)) )*sin(dble(n)*theta)
  342   end do
  343 !$OMP END PARALLEL DO
  344 
  345   ! disp('DONE COMPUTING COEFICIENTS')
  346   du = dreal(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
  347   dv = dreal(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
  348 
  349   ! for GetDP: return imaginary part in dut, dvt
  350   dut = dimag(exp(omega*(0.d0,1.d0)*t)*(cos(theta)*p-sin(theta)*q))
  351   dvt = dimag(exp(omega*(0.d0,1.d0)*t)*(sin(theta)*p+cos(theta)*q))
  352 
  353 end subroutine cylindrical_walloutabc2