1! surface_waves_concave 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
22subroutine surface_waves_concave(du,dv,dut,dvt,X,Y,t,nmode,c,B)
23  !
24  implicit none
25  integer :: nmode
26  double precision :: x,y,du,dv,dvt,dut,t
27  double precision :: r
28  double complex   :: B,q1,q2,v1,v2,q,v,c,ka,kb,omega
29  double precision , parameter :: pi = acos(-1.d0)
30  double precision :: A,sigma,alpha,beta,lambda,mu,rho
31
32  integer, parameter  :: n = 10
33  integer :: nm,np
34  double precision :: ZR, ZI, FNU
35  double precision :: CYRP(n), CYRM(n),cyr(n),CYI(n)
36  integer  :: KODE, M, NZ, IERR
37
38  r = dsqrt(X**2+Y**2)
39  sigma = atan2(Y,X)
40  rho    = 1.d0
41  lambda = 1.d0
42  mu     = 1.d0
43  omega = dble(nmode)*c
44  alpha = dsqrt((lambda+2.d0*mu)/rho)
45  beta = dsqrt(mu/rho)
46  Ka = omega/alpha
47  Kb = omega/beta
48  A = 1.d0
49
50  FNU = 0.d0
51  KODE = 1
52  M = 2
53  nm = nmode-1
54  np = nmode+1
55
56  ZR = dreal(Ka*r)
57  ZI = dimag(Ka*r)
58
59  call ZBESH(ZR, ZI, FNU, KODE, M, n, CYR, CYI, NZ, IERR)
60  q1 = A*Ka/2.d0*(dcmplx(cyr(nm+1),cyi(nm+1))-dcmplx(cyr(np+1),cyi(np+1)))
61  v1 = A*(0.d0,1.d0)*(dble(nmode)/r)*dcmplx(cyr(nmode+1),cyi(nmode+1))
62  ZR = dreal(Kb*r)
63  ZI = dimag(Kb*r)
64  call ZBESH(ZR, ZI, FNU, KODE, M, n, CYR, CYI, NZ, IERR)
65  q2 = B*(0.d0,1.d0)*(dble(nmode)/r)*dcmplx(cyr(nmode+1),cyi(nmode+1))
66  v2 = -B*Kb/2.d0*(dcmplx(cyr(nm+1),cyi(nm+1))-dcmplx(cyr(np+1),cyi(np+1)))
67
68  Q = (q1+q2)
69  V = (v1+v2)
70  du = dreal((cos(sigma)*Q-sin(sigma)*V)*zexp((0.d0,1.d0)*(omega*t+dble(nmode)*sigma)))
71  dv = dreal((sin(sigma)*Q+cos(sigma)*V)*zexp((0.d0,1.d0)*(omega*t+dble(nmode)*sigma)))
72
73  dut = real(omega*(0.d0,1.d0)*(cos(sigma)*Q-sin(sigma)*V)*&
74       exp((0.d0,1.d0)*(omega*t+dble(nmode)*sigma)))
75  dvt = real(omega*(0.d0,1.d0)*(sin(sigma)*Q+cos(sigma)*V)*&
76       exp((0.d0,1.d0)*(omega*t+dble(nmode)*sigma)))
77
78end subroutine surface_waves_concave
79