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