PROGRAM Lamb IMPLICIT NONE real*8 G,E,v,d,c1,c2,Err,start,fin,nod,step,F,Pi real*8 zj(0:1,0:20),gp(0:31),gw(0:31),gpp(0:31),gww(0:31) real*8 nstart,per,nnod,rdef,zdef real*8 tstart,T,w,EDAC,rpbase,icr real*8 ur(1000),uz(1000),a_trig,t0_trig INTEGER NJZeros,NF,ngp complex*16 pdef,Vu_r,Vu_z,wc Dimension Vu_R(1000),Vu_Z(1000) Dimension nnod(3,2) Dimension rdef(3),zdef(3),pdef(5) common /com0/ G,E,v,d common /com1/ c1,c2 common /com2/ start,fin,nod,step common /com3/ Err common /com4/ zj,NJZeros common /com5/ gp,gw,gpp,gww,ngp common /com6/ a_trig,t0_trig common /com7/ VU_R,VU_Z common /com8/ nstart,per,nnod common /com9/ rdef,zdef,pdef common /com10/ tstart,T,wc,w,NF,EDAC common /com11/ Pi common /com12/ F common /com13/ rpbase common /com14/ ur,uz C C main program C open (6, file='data.mak') open (7, file='gra.mak') open (8, file='bug.mak') call initialize write(6,*) 'c1=',c1,' c2=',c2,'cr=',icr(),' E=',E,' v=',v,' d=',d write(6,*) 'c1/c2=3**.5 if v=0.25 :',c1/c2,' = ',dsqrt(3) write(6,*) 'rpb ',rpbase,' start:',start,' fin:',fin,' nod:',nod Call Dmenu close(8) close(7) close(6) end c_____________________________________________________________________________ c subroutine Initialize c Initializes most common (global) variables c_____________________________________________________________________________ subroutine INITIALIZE implicit none real*8 G,E,v,d,c1,c2,Err,start,fin,nod,step real*8 nstart,per,nnod,rdef,zdef,pi,F,rpbase real*8 w,EDAC,T,tstart,nr_ray,a_trig,t0_trig complex*16 pdef,wc INTEGER nf,ch Dimension nnod(3,2) Dimension rdef(3),zdef(3),pdef(5) common /com0/ G,E,v,d common /com1/ c1,c2 common /com2/ start,fin,nod,step common /com3/ Err common /com6/ a_trig,t0_trig common /com8/ nstart,per,nnod common /com9/ rdef,zdef,pdef common /com10/ tstart,T,wc,w,NF,EDAC common /com11/ Pi common /com12/ F common /com13/ rpbase C C Material Properties and wave speeds C E=4d6 v=0.25d0 d=2.0d0 c c sandstone constants c E=14.4d6 c v=0.2d0 c d=2.8d0 call setwaves(E,v,d) C Define Pi constant Pi=3.1415926536d0 c C G material constant c G=E/(2*(1+v)) c C Force acting vertically towards the half-space (compression) c F=1 C C Bromwitz Path start and finish imaginary for Plot_Urz C start=0 fin=15000 c C for every 100 nodes it's 20sec time c nod=901d0 step=(fin-start)/(nod-1) C C Plot_Nrz limits definition C nstart=0 c make per=0.99 (better than 0.999 since it tends to cluster points around rp per=0.99d0 nnod(1,1)=600 nnod(2,1)=400 nnod(3,1)=600 nnod(1,2)=20 nnod(2,2)=20 nnod(3,2)=960 C C Plot_Nrz default values to plot C rdef(1)=1 rdef(2)=10 rdef(3)=100 zdef(1)=0 zdef(2)=1 zdef(3)=10 pdef(1)=dcmplx(0,1) pdef(2)=dcmplx(0.5,1) pdef(3)=dcmplx(1,1) pdef(4)=dcmplx(2,1) pdef(5)=dcmplx(3,2) C C Acceptable Error for Newton-Raphson C Err=1.0d-12 c c compute basic rp for p=0+i to use in IRP, after Initializing Err c rpbase=nr_ray(dcmplx(0d0,1d0),1d0) c C Initialize DAC variables c tstart=0.0d0 EDAC=1E-8 c c Initialize triangular Ltrig: t0=rise time, a*t0=total time c a_trig=2 t0_trig=1 c C Initialize array with bessel zeros c Call BZeros c C Initialize Gauss 20p or 32p c write(*,*) 'Select Gauss points (Ngp=20,Ngp=32)' read(*,1000) ch Call SetGauss(ch) 1000 format(I2) return end c_____________________________________________________________________________ c subroutine DMenu C Displays and operates Menu Selection Screen c_____________________________________________________________________________ subroutine DMenu implicit none Integer Choice 300 Continue Write(*,*) 'Gas Program Control Menu' Write(*,*) '------------------------' Write(*,*) '1: Plot Ray, Nr, Nz functions' Write(*,*) '2: Plot Bessel J0, J1 functions' Write(*,*) '3: Test Longmann rule on J0, J1' Write(*,*) '4: Test N-R root finder for Rayleigh' Write(*,*) '5: Plot U_r, U_z along Bromwitz Path' Write(*,*) '6: Perform DAC Inversion on U_r,U_z' Write(*,*) '0: Exit to DOS' Write(*,*) ' ' Write(*,*) 'Your Choice ? : ' Read(*,1400) Choice Select case(choice) case(0) goto 310 case(1) Call PLOT_Nrz case(2) Call PLOT_J case(3) Call Comp_J01 case(4) Call Check_RP case(5) Call PlotU_rz case(6) Call Laplinv case default Goto 300 end select write(*,*) 'Press ' read (*,*) goto 300 310 write(*,*) 'Have a happy dos' 1400 FORMAT ( I1 ) return end c_____________________________________________________________________________ c subroutine werror : write error subroutine. c takes error codes from other routines and displayes the appropriate c error messages c_____________________________________________________________________________ subroutine werror(code) implicit none integer code character*50 emsg select case(code) case (1) emsg='K2A2: Re(k2a1)=0, probably p=0+0i' case (2) emsg='K2A1: Re(k2a2)=0, probably p=0+0i' case (3) emsg='INTGJ: Selections other than n=0,1,2 not supported' case (4) emsg='MZJ: No Jn zeros available for n<0' case (5) emsg='INTGJ: Maximum N is 100' case (6) emsg='subroutine Euler: M>=N' case (7) emsg='subroutine Euler: Sign=0' case (8) emsg='subroutine ibp: i is 1,2 only' case (9) emsg='subroutine Gauss: fun is r,z only' case (10) emsg='subroutine U_rz1 or 2: fun is r,z only' case (11) emsg='subroutine Setgauss: only 20p or 32p supported' case (12) emsg='subroutine U_rz1 or 2: fun is h,t only' case (13) emsg='subroutine cel: qqc must be <>0' case default emsg='Unknown Error' end select write(*,1600) code,emsg 1600 format(1x,'ECode ',I2,' ',a50) write(*,*) 'Program Halted' stop return end c_____________________________________________________________________________ c Calls Test routine Intgj for EPAL, ET. c_____________________________________________________________________________ subroutine Comp_J01 implicit none real*8 Intgj,intj0,intj1 Integer M,N,G c write(*,*) 'Select number of terms to use (max=100)' write(*,*) 'take care : for N>30 1/2**n kills ET and 1/(ei-ej)' write(*,*) 'kills EPAL. But nothing kills EPAL2' read(*,*) N write(*,*) 'Select terms to exclude from Euler (0..20), best=1:' read(*,*) M write(*,*) 'Select Gauss integrations per interval (1..10)' read(*,*) G write(*,*) intj0=IntgJ(0,0,N,G,3) intj1=IntgJ(1,0,N,G,3) write(*,800) 0,intJ0,1-intJ0 write(*,800) 1,intJ1,1-intJ1 write(*,*) intj0=IntgJ(0,0,N,G,2) intj1=IntgJ(1,0,N,G,2) write(*,760) 0,intJ0,1-intJ0 write(*,760) 1,intJ1,1-intJ1 write(*,*) intj0=IntgJ(0,M,N,G,0) intj1=IntgJ(1,M,N,G,0) Write(*,700) 0,intJ0,1-intJ0 Write(*,700) 1,intJ1,1-intJ1 write(*,*) intj0=IntgJ(0,0,N,G,1) intj1=IntgJ(1,0,N,G,1) write(*,750) 0,intJ0,1-intJ0 write(*,750) 1,intJ1,1-intJ1 700 Format(1X,'ET : I(0,+oo) J',I1,'=',F18.12,' =1?, Err:',E14.6) 750 Format(1X,'EPAL : I(0,+oo) J',I1,'=',F18.12,' =1?, Err:',E14.6) 760 Format(1X,'EPAL2: I(0,+oo) J',I1,'=',F18.12,' =1?, Err:',E14.6) 800 Format(1X,'No accel I(0,+oo) J',I1,'=',F18.12,' =1?, Err:',E14.6) return end c_____________________________________________________________________________ C Checks Rayleigh pole accuracy c_____________________________________________________________________________ subroutine Check_Rp implicit none complex*16 p,ray real*8 irp,ab write(*,1060) write(*,*) 'Give an abscissa (ab>0)' Read(*,*) ab P=dcmplx(0,ab) write(*,1061) p,irp(p),ray(irp(p),p) P=dcmplx(1,1) write(*,1061) p,irp(p),ray(irp(p),p) P=dcmplx(0.999,0.999) write(*,1061) p,irp(p),ray(irp(p),p) P=dcmplx(0.8,-0.81) write(*,1061) p,irp(p),ray(irp(p),p) P=dcmplx(0.1,1) write(*,1061) p,irp(p),ray(irp(p),p) write(*,*) 'In general Ray(x) has a true pole only for p=0+g*i' write(*,*) 'And Ray(0) is never 0 so a pole of 0 means no pole' write(*,*) 'the poles found here for p=a+bi(a<>0,b<>0) are for' write(*,*) 'p=0+bi. Stll these are points of anomaly since ray' write(*,*) 'reaches very close to zero' 1060 format(11X,'p',14X,'rpole',14X,'Ray(rpole)=0 ?') 1061 format (1X,2F10.6,3E15.8) return end c_____________________________________________________________________________ C Plot Points For Bessel functions c_____________________________________________________________________________ subroutine PLOT_J implicit none Integer i,j,pin,choice,n real*8 mzj,bj1,bj0,x0,s0,x1,s1,r c call openf(13,'bessel.dat ') write(*,*) 'Default input (0) or user input (1)' read(*,*) choice if (choice.eq.0) then write(*,*) '8-point plot between consecutive zeros,' write(*,*) '121 zeros plotted (21 tabulated + 100 MacMahon)' pin=8 n=121 r=1d0 else Write(*,*) 'How many points plotted between consecutive zeros ?' Read(*,*) pin write(*,*) 'How many zeros plotted?' Read(*,*) n write(*,*) 'Magnifying Parameter r (r>0)' Read(*,*) r endif write(13,1053) pin,N,N+(N-1)*pin,r write(13,1050) c c plot N zeros c do 104 j=0,N x0=mzj(0,j)/r s0=( mzj(0,j+1)-mzj(0,j) ) / (r*(pin+1)) x1=mzj(1,j)/r s1=( mzj(1,j+1)-mzj(1,j) ) / (r*(pin+1)) write(13,1052) x0,bj0(r*x0),x1,bj1(r*x1) c C if j=final zero of Bessel then j+1 zero isn't known, so C why plot more points ? c if (j.lt.N) then do 103 i=1,pin x0=x0 + s0 x1=x1 + s1 write(13,1051) x0,bj0(r*x0),x1,bj1(r*x1) 103 continue endif 104 continue close(13) 1050 format(9X,'x0',11X,'J0(x0)',12X,'x1',11X,'J1(x1)') 1051 format(1X,4F16.10) 1052 format(1X,4F16.10,' j0(x0)=j1(x1)=0') 1053 format(1X,'PIN: ',I5,'Nzeros: ',I5,' Total Points: ',I5, &'Parameter r',f10.4 ) return end c ____________________________________________________________________________ c Opens file for output using either user input or a default filename c_____________________________________________________________________________ subroutine openf(handle,fname) implicit none integer handle,choice character*12 fname write(*,*) 'Output file name, 0=',fname,' ,1=new' read(*,*) choice if (choice.eq.1) then write(*,*) 'Enter something like xxxxx.dat' read(*,2200) fname endif open(handle,file=fname,Status='UNKNOWN') 2200 format(a12) return end c***************************************************************************** c* PLOT Integrands. * c* A collection of useful routines that can plot accurately Nr,Nz * c* or their alternative forms. * c* Including: * c* Plot_nrz : Master input routine * c* Plotone : Plot single Nrz graph * c* Plotalt : Plot alternative Nrz graph * c***************************************************************************** c_____________________________________________________________________________ c Plot Points For Integrands Nr,Nj for various p,r,z c This is a master routine that takes user input and then calls upon c Plotone, or Plotalt to do the actual plotting. c_____________________________________________________________________________ subroutine PLOT_Nrz implicit none complex*16 p,pdef real*8 r,z,a,b real*8 rdef,zdef Dimension rdef(3),zdef(3),pdef(5) Integer i,j,k,ch1,ch2,it common /com9/ rdef,zdef,pdef Write (*,*) 'User Input (0) or Default Input (1) ? ' Read (*,2300) Ch1 Write (*,*) 'Plot around pole (1) or General Plot (2) ? ' Read (*,2300) Ch2 write (*,*) 'function: (0) Nrz, (1) Nrz/Ji(rx)' write(*,*) '(2) Nrz~=Nrz((1+t)/(1-t))*2/(1-t)**2, (3)Nrz~/Ji(x*r)' read(*,2300) it if (ch1.eq.0) then 143 Continue Write(*,*) ' User Input Selected. Provide Data:' Write(*,*) 'Radial Distance ( r>0 ): ' Read(*,*) r Write(*,*) 'Depth ( z>0 ): ' Read(*,*) z Write(*,*) 'Re(P) of Laplace: ' Read(*,*) a Write(*,*) 'Im(P) of Laplace: ' Read(*,*) b p=dcmplx(a,b) call openf(10,'nrz0.dat ' ) if ((it.eq.0).or.(it.eq.1)) then Call Plotone(r,z,p,ch2,it) else Call Plotalt(r,z,p,it) endif close(10) Write (*,*) 'Plot another? ( 1=yes, 0=no ) ' Read (*,2300) Ch1 if (ch1.eq.1) Goto 143 else do 140 i=1,5 do 141 j=1,5 do 142 k=1,5 Call PlotOne(rdef(i),zdef(j),pdef(k),ch2,it) 142 Continue 141 Continue 140 Continue endif 2300 format(i1) return end c_____________________________________________________________________________ c Plots one graph for Nr ,Nz, Ray, detailed or general c c_____________________________________________________________________________ subroutine PlotOne(r,z,p,type,it) implicit none complex*16 ray,nr,nz,p real*8 x,r,z,nstep,nfin,rp,irp,ibp real*8 nstart,per,nnod,mzj,bj0,bj1 Dimension nnod(3,2) common /com8/ nstart,per,nnod Integer i,type,it Write(10,*) 'R: ',r Write(10,*) 'Z: ',z Write(10,*) 'P: ',dreal(p),dimag(p) rp=irp(p) write(10,*) 'bp1: ',ibp(1,p),'bp2: ',ibp(2,p),'rp: ',rp c c if rp=0 no pole was found c IF (rp.ne.0) then if (type.eq.1) then c c Plot around rpole and just a little further c Careful Nr,Nz are zero at zj(n,k)/r not zj(n,k) c because they contain J(x*r) c nstart=0.5*rp nfin=max(8*rp,mzj(0,1)/r,mzj(1,1)/r) else c c else Plot at least five zeros c nstart=0 nfin=100*max(200*rp,mzj(0,20)/r,mzj(1,20)/r) endif else write(10,*) 'Rp for p=0+bi: ',rp if (type.eq.1) then nstart=0.5*rp nfin=max(8*rp,mzj(0,1)/r,mzj(1,1)/r) else nstart=0.5*rp nfin=100*max(200*rp,mzj(0,20)/r,mzj(1,20)/r) endif endif write(10,1002) c c override to force printing from x=0 always. Useful sometimes c nstart=0 C C nnod1,2,3 = number of nodes before, around and after pole C per < 1 = percentage of rp that will be graphed in detail C before and after rp C C Plot points from nstart to rp*per C always lower by one : 81-1=80 c x=nstart nstep=(rp*per)/(nnod(1,type)+1) do 130 i=1,nnod(1,type) if (it.eq.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z),nz(x,p,r,z) elseif(it.eq.1) then if (x.ne.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z)/bj1(x*r), & nz(x,p,r,z)/bj0(x*r) else write (10,1001) x,ray(x,p),0d0,0d0 endif endif x=x+nstep 130 continue c C Plot points from rp*per to (2-per)*rp C Reset x so we avoid errors of x+nstep > rp*per c x=per*rp nstep=rp*(2-2*per)/(nnod(2,type)+1) do 131 i=1,nnod(2,type) if (it.eq.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z),nz(x,p,r,z) elseif(it.eq.1) then if (x.ne.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z)/bj1(x*r), & nz(x,p,r,z)/bj0(x*r) else write (10,1001) x,ray(x,p),0,0 endif endif x=x+nstep 131 continue c C Plot points from (2-per)*rp to nfin c x=(2-per)*rp nstep=(nfin-x)/(nnod(3,type)+1) do 132 i=1,nnod(3,type)+1 if (it.eq.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z),nz(x,p,r,z) elseif(it.eq.1) then if (x.ne.0) then write (10,1001) x,ray(x,p),nr(x,p,r,z)/bj1(x*r), & nz(x,p,r,z)/bj0(x*r) else write (10,1001) x,ray(x,p),0,0 endif endif x=x+nstep 132 continue 1001 format(1X,' ',E18.11,1X,2E18.11,1X,2E18.11,1X,2E18.11) 1002 format(1X,9X,'X',16X,'Real(R)',10X,'Imag(R)',10X,' Real(Nr)' & ,12X,'Imag(Nr)',12X,'Real(Nz)',12X,'Imag(Nz)') return end c_____________________________________________________________________________ c Plots one graph for Nr ,Nz, Ray transformed by the c Stieljes Trasformation c_____________________________________________________________________________ subroutine Plotalt(r,z,p,it) implicit none complex*16 ray,nr,nz,p real*8 start,fin,x,nod,r,z,bj0,bj1,t,dx integer i,it c Write(10,*) 'R: ',r Write(10,*) 'Z: ',z Write(10,*) 'P: ',dreal(p),dimag(p) write(10,1002) c start=-0.999999 fin=+0.999999 nod=2001 c do 100 i=1,nod t=start+i*(fin-start)/(nod-1) x=(1+t)/(1-t) dx=2/(1-t)**2 if (it.eq.2) then write (10,1001) t,ray(x,p),nr(x,p,r,z)*dx,nz(x,p,r,z)*dx else write (10,1001) t,ray(x,p),nr(x,p,r,z)/bj1(x*r), & nz(x,p,r,z)/bj0(x*r) endif 100 continue 1001 format(1X,' ',E18.11,1X,2E18.11,1X,2E18.11,1X,2E18.11) 1002 format(1X,9X,'X',16X,'Real(R)',10X,'Imag(R)',10X,' Real(Nr)' & ,12X,'Imag(Nr)',12X,'Real(Nz)',12X,'Imag(Nz)') return end c***************************************************************************** c* GAUSS Integration Rule * c* Useful routines for initializing and implementing * c* a 20, or 32-point Gauss rule * c***************************************************************************** c_____________________________________________________________________________ c Initialize Gauss points and Weights c_____________________________________________________________________________ subroutine SETGAUSS(N) implicit none real*8 gp(0:31),gw(0:31),gpp(0:31),gww(0:31) real*8 sumgp,sumgw integer ngp,i,N common /com5/ gp,gw,gpp,gww,ngp c ngp=N select case(ngp) case (20) c C Gauss-20 points c gp( 0) =-0.9931285991850949 gp( 1) =-0.9639719272779137 gp( 2) =-0.9122344282513259 gp( 3) =-0.8391169718222188 gp( 4) =-0.7463319064601507 gp( 5) =-0.6360536807265150 gp( 6) =-0.5108670019508270 gp( 7) =-0.3737060887154195 gp( 8) =-0.2277858511416450 gp( 9) =-0.0765265211334973 c C Gauss-20 Weights c gw( 0) = 0.0176140071391521 gw( 1) = 0.0406014298003869 gw( 2) = 0.0626720483341090 gw( 3) = 0.0832767415767047 gw( 4) = 0.1019301198172404 gw( 5) = 0.1181945319615184 gw( 6) = 0.1316886384491766 gw( 7) = 0.1420961093183820 gw( 8) = 0.1491729864726037 gw( 9) = 0.1527533871307258 case(32) c C Gauss-32 points c gp( 0) =-0.9972638618494815 gp( 1) =-0.9856115115452683 gp( 2) =-0.9647622555875064 gp( 3) =-0.9349060759377396 gp( 4) =-0.8963211557660521 gp( 5) =-0.8493676137325699 gp( 6) =-0.7944837959679424 gp( 7) =-0.7321821187402896 gp( 8) =-0.6630442669302152 gp( 9) =-0.5877157572407623 gp(10) =-0.5068999089322293 gp(11) =-0.4213512761306353 gp(12) =-0.3318686022821276 gp(13) =-0.2392873622521370 gp(14) =-0.1444719615827964 gp(15) =-0.0483076656877383 c C Gauss-32 Weights c gw( 0) = 0.0070186100094701 gw( 1) = 0.0162743947309057 gw( 2) = 0.0253920653092621 gw( 3) = 0.0342738629130214 gw( 4) = 0.0428358980222267 gw( 5) = 0.0509980592623761 gw( 6) = 0.0586840934785355 gw( 7) = 0.0658222227763618 gw( 8) = 0.0723457941088485 gw( 9) = 0.0781938957870703 gw(10) = 0.0833119242269468 gw(11) = 0.0876520930044038 gw(12) = 0.0911738786957639 gw(13) = 0.0938443990808046 gw(14) = 0.0956387200792749 gw(15) = 0.0965400885147278 case default call werror(11) end select c c Compute the remaining half gw, gp arrays c do 100 i=0,ngp/2-1 gp(i+ngp/2)=-gp(ngp/2-1-i) gw(i+ngp/2)= gw(ngp/2-1-i) 100 continue sumgp=0 sumgw=0 do 200 i=0,ngp-1 sumgp=sumgp+gp(i) sumgw=sumgw+gw(i) 200 continue write(8,*) 'Gauss-',ngp,' selected' write(8,*) 'SetGauss SumGp:',sumgp,' = 0?' write(8,*) 'SetGauss SumGw:',sumgw,' = 2?' return end c_____________________________________________________________________________ c Adjust Gauss points and Weights from [-1,1] to [a,b] c_____________________________________________________________________________ subroutine ADJGAUSS( A,B ) implicit none DIMENSION gp(0:31),gw(0:31),gpp(0:31),gww(0:31) real*8 gp,gw,gpp,gww,a,b,c,d integer ngp,i common /com5/ gp,gw,gpp,gww,ngp C c=(b-a)/2 d=(a+b)/2 do 100 i=0,ngp-1 gpp(i)= c * gp(i) + d gww(i)= c * gw(i) 100 continue return end c_____________________________________________________________________________ c subroutine Composite GAUSS c uses a gauss 12-point rule to integrate between a,b by c dividing [a,b] in N subdivisions. Totally we get a 12*N gauss rule c_____________________________________________________________________________ complex*16 function CG(A,B,N,fun,p,r,z) implicit none real*8 A,B,r,z,step,bj0,bj1,cn,dn,gpa real*8 gp(0:31),gw(0:31),gpp(0:31),gww(0:31) complex*16 sum,nr,nz,p integer ngp,N,i,j character*1 fun common /com5/ gp,gw,gpp,gww,ngp c c careful A=1 c sum=dcmplx(0,0) step=(B-A)/N cn=step/2 do 100 i=1,n dn=A + (i-0.5d0)*step do 200 j=0,ngp-1 gpa=cn * gp(j) + dn select case(fun) case('r') sum=sum+gw(j)*Nr(gpa,p,r,z) case('z') sum=sum+gw(j)*Nz(gpa,p,r,z) case('0') sum=sum+gw(j)*bj0(gpa) case('1') sum=sum+gw(j)*bj1(gpa) case default call werror(9) end select 200 continue 100 continue cg=cn*sum return end c***************************************************************************** c* Rayleigh poles, branch points and speed * c* a set of functions for easily resolving and computing all these values * c* functions contained * c* ibp(i,p) : branch points * c* irp(p) : pole * c* icr : wave speed cr * c* ispole(p): which p gives a pole * c* nr_ray(p): N-R pole finder * c***************************************************************************** c_____________________________________________________________________________ c Rayleigh wave speed c c2/cr>1 always c_____________________________________________________________________________ real*8 function icr() implicit none real*8 c1,c2,rpbase common /com1/ c1,c2 common /com13/ rpbase icr=1/rpbase return end c_____________________________________________________________________________ c Branch points for Ray,Nr,Nz c 00 the pole of p=0+ib will be given. Check with ispole to be sure c if the value returned is the true pole (p=0+ib) or an anomaly (p=a+ib) c or just check for dreal(p)=0. If p=0+0i then irp=0 automatically. Very c conviniently denoting a no pole (?) condition. c_____________________________________________________________________________ real*8 function Irp( P ) implicit none complex*16 p real*8 rpbase common /com13/ rpbase irp=rpbase*dimag(p) return end c_____________________________________________________________________________ C easy way to accomodate your experience for the existence of C a Ray pole. Change this and every check follows C popular rule : abs(dreal(p)).le.abs(Dimag(p)) ==> ispole c_____________________________________________________________________________ LOGICAL function IsPole( P ) implicit none complex*16 P If ((Dreal(p).eq.0).and.(dimag(p).ne.0)) then ispole=.true. else ispole=.false. endif return end c_____________________________________________________________________________ c RAYLEIGH POLE (with use of Newton-Raphson) c a pole exists only for p=0+bi, b<>0. c this routine needs a Seed to start seeking for Rp and the c complex (not just the imaginary part) of p. A maximum of c 100 iterations will be completed before a "no pole" condition is c declared by indicating rp=0 (rp is always <>0). This might happen even c for p=0+bi if b<1e-8 approximately. Maybe a good seed might solve that. c c Careful. If you build using "Improve Floating point consistency" the compiler c evaluates everything using double accuracy. So the NR_ray((0,1),1) for rbase c will not converge for an err=1e-12, but for err=1e-10. IF you disable this c optimizations the program will be executed on 80487 32bit registers employing c extended accuracy. This permits the best possible accuracy obtained!!!! c So this build enable option doesn't improve anything at all!! c_____________________________________________________________________________ real*8 function NR_Ray( P,Se ) implicit none Integer I logical ispole real*8 Err,x,xk,xs,Se complex*16 p,ray,dray common /com3/ Err c write(8,1102) p c c In general a Se=dimag(p)/100 is very good c x=Se i=0 c c If no pole exists send rp=0 c if (.not.ispole(p)) then NR_Ray=0 goto 210 endif c 200 continue xs=x xk=x - dreal(ray(x,p)/dray(x,p)) x=xk i=i+1 write(8,1100) i,x,ray(x,p),(x-xs),err If(I.gt.100) then Write(8,1103) p c c Pole always <>0 so this is to show that none was found. Exit routine c NR_Ray=0 Goto 210 endIF if ((ABS(x-xs).gt.Err).OR.(abs(ray(x,p)).gt.Err)) goto 200 NR_RAY=x 210 Continue 1100 format(1X,I4,5E20.12) 1101 format(1X,'No rayleigh pole for P= ',2E12.5) 1102 format(1X,'NR_RAY, seeking rayleigh pole for P= ',2E12.5) 1103 format(1X,'Error : NR_RAY not converging. Wrong Seed ? ',2E12.5) return end c***************************************************************************** c* DAC Laplace transform inversion * c* A complete set of routines to perform the inversion * c* * c***************************************************************************** c_____________________________________________________________________________ c subroutine LAPLINV c master routine: Interacts with the user and coodinates the rest of the c routines for optimal results c_____________________________________________________________________________ subroutine LAPLINV implicit none real*8 r,z,urw,uzw,q,tnorm,c1,c2 real*8 T,Tmax,tnods,tstart,w,EDAC complex*16 wc integer nf character*1 Lfun common /com1/ c1,c2 common /com10/ tstart,T,wc,w,NF,EDAC Write(*,*) 'R=' Read(*,*) r Write(*,*) 'Z=' Read(*,*) z Write(*,*) 'Time dependance of load (H:Heaviside, T:Triangular):' Read(*,1000) Lfun Write(*,*) '2*NF+1 terms for Fourier series (NF=2,3,..99)=' Read(*,*) NF Write(*,*) 'Tnorm Maximum: (normalized, Tnorm=2 suggested )' Read(*,*) tnorm Write(*,*) 'Time nodes to be used (tnods<1000):' Read(*,*) tnods C ** just a precaution IF (tnorm.eq.0d0) tnorm=2d0 If (tnods.eq.0) tnods=300 If ((NF.le.0).or.(Nf.gt.99)) nf=10 write(*,*) 'Select displacements output' call openf(12,'urz0.dat ') write(*,*) 'Select Laplace transformed displacements output' call openf(14,'urz_0.dat ') c C De-Normalize time c use rd=(r^2+z^2)^0.5 for proper normalization. Using r works only for z=0 c Tmax=dsqrt(r**2+z**2)*tnorm/c2 c T=0.8*tax or 2*tmax?? T=0.8*tmax c c output primary variables c write(12,*) 'r:',r,' z:',z write(12,*) 'tnorm:',tnorm,' tnods:',tnods,' T:',T,' tmax:',tmax write(12,*) 'NF:',nf,' time dependance ',lfun c c Careful, T is very small so w skyrockets!!! c q=0 w= q - Dlog(EDAC)/(2*T) wc=dcmplx(w,0) Call SetVU_rz(r,z,urw,uzw,Lfun) Call Dac(r,z,urw,uzw,tmax,tnods) Close(12) close(14) 1000 format(a1) return end c_____________________________________________________________________________ c subroutine SetVU_RZ c sets and computes values of u_r, u_z necessary for DAC inversion c returns the global arrays Vu_r, Vu_z, and variables urw, uzw c_____________________________________________________________________________ subroutine SetVU_RZ(r,z,urw,uzw,lfun) implicit none real*8 r,z,start,fin,nod,step,irp,pi real*8 T,tstart,w,urw,uzw,EDAC,rp complex*16 vu_r(1000),vu_z(1000),u_rz2,u_rz,p,wc Integer k,nf character*1 lfun common /com2/ start,fin,nod,step common /com7/ VU_R,VU_Z common /com10/ tstart,T,wc,w,NF,EDAC common /com11/ Pi C C urw=0.5d0*dreal(U_rz2('r',lfun,wc,r,z)) uzw=0.5d0*dreal(U_rz2('z',lfun,wc,r,z)) write(14,*) 'w=',w,' urw=',urw,' uzw=',uzw,' T=',t write(14,*) 'Bromwitz abscissa :',wc write(*,*) 'All poles displayed are for p''=0+dimag(p)*i' write(14,1013) Do 100 K=1, 2*NF+1 p= dcmplx( w, K*Pi/T ) VU_R(k)=U_rz2('r',lfun,p,r,z) VU_Z(k)=U_rz2('z',lfun,p,r,z) rp=irp(p) write(14,1012) p,VU_R(k),VU_z(k),rp write(*,1004) p,VU_R(k),VU_z(k),rp 100 Continue 1004 format(1X,2F10.3,2E12.5,2E12.5,E12.5) 1012 format(1X,2F10.3,2E12.5,2E12.5,E12.5,' Pole for p''=0+bi') 1013 format(5X,'Re(p)',5X,'Im(p)',5X,'Re(U_r)',5X,'Im(U_r)', & 5X,'Re(U_z)',5X,'Im(U_z)',5X,'Pole') return end c_____________________________________________________________________________ c standard EPSILON: series summing accelaration algorithm c_____________________________________________________________________________ subroutine EPSILON(NF,S,Efin) implicit none real*8 E(0:200,1:200),S(200),EFIN INTEGER nf,ip,m DO 10 M=1, 2*NF+1 E(0,M)=0.0D0 E(1,M)=S(M) 10 CONTINUE DO 20 IP=1, 2*NF DO 30 M=1,((2*NF+1)-IP) E(IP+1,M)=E(IP-1,M+1)+(E(IP,M+1)-E(IP,M))**(-1.0D0) c c careful here, for large NF a 0**(-1) might occur c 30 CONTINUE 20 CONTINUE Efin=E(2*NF+1,1) return end c_____________________________________________________________________________ c modified EPSILON2: series summing accelaration algorithm c incorporates dynamic techniques that improve robustness c_____________________________________________________________________________ SUBROUTINE EPSILON2(NF,S,Efin) implicit none DOUBLE PRECISION E(0:200,1:200),S(200),Efin INTEGER NF,IP,I,IC,m C C *** Initialization C * Set E(P=0,M=1..2*N+1)=0, E(P=1,M=1..2*N+1)=S(M) C DO 10 M=1, 2*NF+1 E(0,M)=0.0d0 E(1,M)=S(M) 10 CONTINUE C C Main loop :IP is the P in the symbol e(p,m) that we use in our thesis, c I is used to cut the epsilon tree together C with IC ( makes the inner loop shorter ). This has the effect c of dynamically reducing NF C I=0 DO 20 IP=1, 2*NF IC=I if ( ((2*NF +1)-IP-IC).eq.0) then I=(ip-1)+(mod((ip-1),2)-1) EFIN = E(i,1) write(8,*) 'Epsilon cut at ip:',I goto 40 endif DO 30 M=1,((2*NF+1)-IP)-IC IF (E(IP,M+1).eq.E(IP,M)) then c C Now (e(p,m+1)-e(p,m)=0, so C define how many more E(IP,M) to disregard for next loop C I=(2*NF+1)-IP-M+1 goto 20 ELSE E(IP+1,M)=E(IP-1,M+1)+ 1d0 /(E(IP,M+1)-E(IP,M)) ENDIF 30 CONTINUE 20 CONTINUE c c no cut was needed c EFIN = E(2*NF+1,1) c write(8,*) 'Epsilon full ip:',2*NF+1,' no cut was needed' 40 Continue c c there was a dynamic cut performed to avoid system collapse c RETURN END c_____________________________________________________________________________ c subroutine DAC c Performs DAC inversion by sumation of the Fourier series using EPAL c_____________________________________________________________________________ subroutine DAC(r,z,urw,uzw,tmax,tnods) implicit none real*8 Ur(1000),Uz(1000),Sr(200),Sz(200) real*8 sum1,sum2,Er,Ez,urw,uzw,r,z,rd real*8 T,tmax,tnods,tstart,w,pi,kt,tt real*8 c1,c2,F,G,E,v,d,EDAC,peur,peuz,tstep real*8 peur1,peuz1,LEr,LEz Integer i,k,nf complex*16 VU_R(1000),VU_Z(1000) complex*16 wc common /com0/ G,E,v,d common /com1/ c1,c2 common /com7/ VU_R,VU_Z common /com10/ tstart,T,wc,w,NF,EDAC common /com11/ Pi common /com12/ F common /com14/ ur,uz C C LEr=0d0 LEz=0d0 rd=dsqrt(r**2+z**2) tstep=(tmax-tstart)/(tnods-1) DO 100 i=1, tnods tt=tstart+(i-1)*tstep write(*,*) 'ok here tt=',tt sum1=0d0 sum2=0d0 do 200 k=1, 2*NF+1 kt=k*Pi*tt/T sum1=sum1+dreal(VU_r(k))*dcos(kt) & -dimag(VU_r(k))*dsin(kt) sum2=sum2+dreal(VU_z(k))*dcos(kt) & -dimag(VU_z(k))*dsin(kt) Sr(k)=sum1 Sz(k)=sum2 200 continue C ** \\\\ sum the series using epsilon algorithm Call EPSILON2(NF,Sr,Er) Call EPSILON2(NF,Sz,Ez) ur(i)=( DEXP(w*tt)/T )*( Urw+Er ) uz(i)=( DEXP(w*tt)/T )*( Uzw+Ez ) c c Compute Integral-type error estimator.Use compound trapezoid rule c the final result is LE*tstep for unnormalized graph or LE*tstep*c2/r c for normalized graph. Since analytical results are for z=0 and Heaviside c ignore error for any other z, or other time dependance. c peur1=peur(r,tt) peuz1=peuz(r,tt) if ((i.eq.1).or.(i.eq.tnods))then LEr=LEr + 0.5*abs(peur1-ur(i)) LEz=LEz + 0.5*abs(peuz1-uz(i)) else LEr=LEr + abs(peur1-ur(i)) LEz=LEz + abs(peuz1-uz(i)) endif C ** output normalized results write(12,1600) tt*c2/rd,ur(i)*(Pi*G*rd)/F,uz(i)*(-Pi*G*rd)/F, & peur1*(Pi*G*rd)/F,peuz1*(-Pi*G*rd)/F 100 continue if (z.eq.0) then write(*,*) 'Integral of Error.' write(*,*) 'Valid only for Heaviside time dependance' write(*,*) 'Trapezoid rule Estimation for normalized time' write(*,*) 'LEr:',LEr*tstep*c2/rd write(*,*) 'LEz:',LEz*tstep*c2/rd endif 1600 format(1X,5E12.5) return end c_____________________________________________________________________________ c subroutine PlotU_RZ c Plots results from u_rz, u_rz2 c_____________________________________________________________________________ subroutine PlotU_RZ implicit none real*8 r,z,start,fin,nod,step real*8 ab,pi,irp complex*16 vu_r(1000),vu_z(1000),u_rz,p,bi Integer i Logical ispole character*1 Lfun common /com2/ start,fin,nod,step common /com7/ VU_R,VU_Z common /com11/ Pi C bi=dcmplx(0,1) Write(*,*) 'R=' Read(*,*) r Write(*,*) 'Z=' Read(*,*) z Write(*,*) 'abscissa(ab>0):' Read(*,*) ab Write(*,*) 'Time dependance of load (H:Heaviside, T:Triangular):' Read(*,*) Lfun call openf(11,'u_rz0.dat ') write (11,*) 'Bromwitz abscissa :',ab write (11,*) 'Time dependance :',LFun write (11,1013) P=dcmplx(ab,start) i=1 119 Continue VU_R(i)=U_rz('r',Lfun,p,r,z) VU_Z(i)=U_rz('z',Lfun,p,r,z) IF (ispole(p)) then write(11,1004) p,VU_R(i),VU_z(i),irp(p) write(*,1004) p,VU_R(i),VU_z(i),irp(p) ELSE write(11,1012) p,VU_R(i),VU_z(i) write(*,1012) p,VU_R(i),VU_z(i) endIF i=i+1 P=P+dcmplx(0,step) 1004 format(1X,2F10.3,2E12.5,2E12.5,E12.5) 1012 format(1X,2F10.3,2E12.5,2E12.4,' No pole') 1013 format(5X,'Re(p)',5X,'Im(p)',5X,'Re(U_r)',5X,'Im(U_r)', & 5X,'Re(U_z)',5X,'Im(U_z)',5X,'Pole') If (dimag(p).lt.fin) GOTO 119 close(11) return end c_____________________________________________________________________________ C Euler Transformation C V : Oscillatory series, DV(r,n): Tables of finite differences C NE : number of terms in series(20 here), m : series term where transform starts C sign : sign of first term of series (V0) c c use tiny,precision,huge to determine where to cut off the algorithm. C Needs NE+1 terms ( 0..Ne ) c_____________________________________________________________________________ subroutine EULER(V,NE,M,esign,sum) implicit none real*8 DV(0:100,0:100),V(0:100),sum INTEGER ne,m,n,r,esign,i if (M.ge.NE) Call werror(6) if (esign.eq.0) Call werror(7) DO 10 N=0, NE DV(0,N)=V(N) 10 CONTINUE DO 20 R=1, NE-M-1 DO 30 N=M,NE-R-1 DV(R,N) = DV(R-1,N+1) - DV(R-1,N) 30 CONTINUE 20 CONTINUE sum=dcmplx(0,0) DO 40 N=0, M-1 sum=sum+((-1)**N)*V(N) 40 CONTINUE DO 50 i=m, NE-M-1 sum=sum+((-1)**i)*DV(i-m,m)/(2**(i-m+1)) 50 CONTINUE sum=sum*esign return end c_____________________________________________________________________________ C Laplace transform of Heaviside c_____________________________________________________________________________ complex*16 function LHeavi(p) complex*16 p LHeavi= 1 / p return end c_____________________________________________________________________________ c Laplace transform of Triangular c Triamgular has a t0 rise time and a total time a*t0 c For surface explosion suggested t0=1msec, a*t0=5msec c_____________________________________________________________________________ complex*16 function LTrig(p,a,t0) complex*16 p real*8 a,t0 LTrig= ((a-1) - a*cdexp(-t0*p) + cdexp(-a*t0*p))/ & ((a-1)*t0*(p**2)) return end c_____________________________________________________________________________ C Standard U_rz Function c employes an adaptive Longmann rule to compute the Laplace tranformed c displacements C u_r=S(nr*dx) C u_z=S(nz*dx) c fun='r','z' chooses u_r,u_z c Lfun='h','t' chooses heaviside or triangular time dependance c_____________________________________________________________________________ complex*16 function U_rz( fun,Lfun,P,R,Z ) implicit none complex*16 p,sum,LHeavi,Ltrig,c0,cxV,cg INTEGER i,k,m,imsign,resign,nte,bj real*8 imV(0:100),reV(0:100),resum,imsum real*8 R,Z,irp,mzj,mzj1,mzj2,rp,off,bp1,bp2,ibp real*8 G,E,v,d,Pi,F,a_trig,t0_trig logical found character*1 fun,Lfun common /com0/ G,E,v,d common /com6/ a_trig,t0_trig common /com11/ Pi common /com12/ F c c c0 = complex zero. Just a speed optimizer c Compute Integrals between zeros of bessel j1(x*r) (= zeros of Nr) c which are the zeros of J1(x) divided by r c Compute the following Integrals c c0=dcmplx(0,0) rp=irp(p) bp1=ibp(1,p) bp2=ibp(2,p) write(*,*) 'Ok here rp=',rp i=0 sum=c0 found=.false. c c compute the next two zeros that define the next integration interval c Avoid rp=0=zj(1,0) this is not a pole c if Rp inside i,i+1 bessel zeros then integrate using cauchy principal value c Careful. If you wish place an if-clause to choose the best integration c method, especially when the anomaly is small ( 1 is best then) c c also careful with bp1> rp found at {',i,',',i+1,'} zeros j1' off=1.0E-1*min( rp-mzj1 , mzj2-rp ) sum=sum + cg(mzj1,rp-off,4,fun,p,r,z) sum=sum + cg(rp-off,rp,2,fun,p,r,z) sum=sum + cg(rp,rp+off,2,fun,p,r,z) sum=sum + cg(rp+off,mzj2,4,fun,p,r,z) found=.true. elseif ((mzj1.le.bp1).and.(bp1.le.mzj2)) then sum=sum + cg(mzj1,bp1,2,fun,p,r,z) sum=sum + cg(bp1,mzj2,2,fun,p,r,z) elseif ((mzj1.le.bp2).and.(bp2.le.mzj2)) then sum=sum + cg(mzj1,bp2,2,fun,p,r,z) sum=sum + cg(bp2,mzj2,2,fun,p,r,z) else sum=sum + cg(mzj1,mzj2,4,fun,p,r,z) endif i=i+1 end do c c Start computing terms to Eulerize c number of terms for ET. Total is (nte+1). Exclude m of them c if (i.lt.50*4) then m=10+i/4 else m=60 endif nte=m+30 do 250 k=0,nte mzj1=mzj(bj,i)/r mzj2=mzj(bj,i+1)/r cxV=cg(mzj1,mzj2,2,fun,p,r,z) if (k.eq.0) then c c Re and Im have different signs. Make sure to store them for Euler c resign=sign( 1, dreal(cxV) ) imsign=sign( 1, dimag(cxV) ) write(8,*) 'Signs of:', cxV,' re:',resign,' im:',imsign endif c c make both parts positive for longman rule with use of ET c ReV(k)=Abs(dreal(cxV) ) ImV(k)=Abs(dimag(cxV) ) i=i+1 250 continue Call Euler(ReV,nte,M,resign,resum) Call Euler(ImV,nte,M,imsign,imsum) select case(Lfun) case('h','H') u_rz=F/(2*Pi*G) * Lheavi(p) & *(sum + dcmplx(resum,imsum) ) case('t','T') u_rz=F/(2*Pi*G) * LTrig(p,a_trig,t0_trig) & *(sum + dcmplx(resum,imsum) ) case default call werror(12) end select return end c_____________________________________________________________________________ c alternative function U_rz: c uses EPAL for acceleration c_____________________________________________________________________________ complex*16 function U_rz2( fun,Lfun,P,R,Z ) implicit none complex*16 p,sum,LHeavi,Ltrig,c0,cxV,cg INTEGER i,k,nte,bj,m real*8 imV(200),reV(200),rsum,isum real*8 R,Z,irp,mzj,mzj1,mzj2,rp,off,bp1,bp2,ibp real*8 G,E,v,d,Pi,F,a_trig,t0_trig logical found character*1 fun,Lfun common /com0/ G,E,v,d common /com6/ a_trig,t0_trig common /com11/ Pi common /com12/ F c cxV:complex, reV, imV real and imag parts of cxV c c0 = complex zero. Just a speed optimizer c Compute Integrals between zeros of bessel j1(x*r) (= zeros of Nr) c which are the zeros of J1(x) divided by r c Compute the following Integrals c c0=dcmplx(0,0) rp=irp(p) bp1=ibp(1,p) bp2=ibp(2,p) write(*,*) 'Ok here rp=',rp i=0 sum=c0 found=.false. c c compute the next two zeros that define the next integration interval c Avoid rp=0=zj(1,0) this is not a pole c also careful with bp1> rp found at {',i,',',i+1,'} zeros j1' off=1.0E-1*min( abs(rp-mzj1) , abs(mzj2-rp) ) sum=sum + cg(mzj1,rp-off,4,fun,p,r,z) sum=sum + cg(rp-off,rp,2,fun,p,r,z) sum=sum + cg(rp,rp+off,2,fun,p,r,z) sum=sum + cg(rp+off,mzj2,4,fun,p,r,z) found=.true. elseif ((mzj1.le.bp1).and.(bp1.le.mzj2)) then sum=sum + cg(mzj1,bp1,2,fun,p,r,z) sum=sum + cg(bp1,mzj2,2,fun,p,r,z) elseif ((mzj1.le.bp2).and.(bp2.le.mzj2)) then sum=sum + cg(mzj1,bp2,2,fun,p,r,z) sum=sum + cg(bp2,mzj2,2,fun,p,r,z) else sum=sum + cg(mzj1,mzj2,4,fun,p,r,z) endif i=i+1 end do c c Start computing terms for EPAL c number of terms for EPAL. c cxV=c0 nte=20 if (i.lt.50*4) then m=10+i/4 else m=60 endif c c exclude some terms from epal to avoid pole influence c do 250 k=1,m mzj1=mzj(bj,i)/r mzj2=mzj(bj,i+1)/r sum=sum+cg(mzj1,mzj2,2,fun,p,r,z) i=i+1 250 continue do 260 k=1,nte mzj1=mzj(bj,i)/r mzj2=mzj(bj,i+1)/r cxV=cxV+cg(mzj1,mzj2,2,fun,p,r,z) ReV(k)=dreal(cxV) ImV(k)=dimag(cxV) i=i+1 260 continue c c remember : nte/2 is integer division : 10/2=5, 11/2=5 !! c Call Epsilon2(nte/2,ReV,rsum) Call Epsilon2(nte/2,ImV,isum) select case(Lfun) case('h','H') u_rz2=F/(2*Pi*G) * Lheavi(p) & *(sum + dcmplx(rsum,isum) ) case('t','T') u_rz2=F/(2*Pi*G) * LTrig(p,a_trig,t0_trig) & *(sum + dcmplx(rsum,isum) ) case default call werror(12) end select return end c_____________________________________________________________________________ c function IntgJ c A test function that computes the integrals of Bessel J0, J1 [0,+oo) c using either c al=0 : ET with M excluded terms c al=1 : EPAL with zero excluded terms c al=2 : EPAL 2 c al=3 : no acceleration c bj=0,1 : bessel j0, j1 c the advantage of ET over EPAL here is being more robust for many terms c but EPAL is more accurate for fewer terms. If the evaluation of terms is c complicated and time consuming, then EPAL is the best. But for easily c computed terms like J0,J1 here, ET might prove faster since it involves c a simpler routine. Might not. Dare you check it ?? c_____________________________________________________________________________ real*8 function IntgJ( ebj,M,N,G,al ) implicit none real*8 mzj,mzj1,mzj2 INTEGER i,N,m,al,ebj,G real*8 sum,V(0:100) complex*16 cg c if (n.gt.100) then call werror(5) write(*,*) stop endif c c Compute Integrals between zeros of bessel J0,J1 (Gauss) c use composite gauss rule with G*12 points (G>=1) c sum=0 Do 110 i=0, N mzj1=mzj(ebj,i) mzj2=mzj(ebj,i+1) if (ebj.eq.0) then V(i)=dreal(cg(mzj1,mzj2,G,'0',dcmplx(0,0),0d0,0d0)) else V(i)=dreal(cg(mzj1,mzj2,G,'1',dcmplx(0,0),0d0,0d0)) endif 109 continue sum=sum+v(i) c write(8,*) 'method ',al,' i ',i,' V(i)=',V(i),' sum=',sum c c Euler starts from V(0) but Epsilon from V(1). c still the array is passed as V(0:100) and automatically for Epsilon c V(0) becomes V'(1) (epsilon reads the first term of the array and calls it V'(1) !!) c so no information is lost. c if (al.eq.0) then V(i)=Abs(V(i)) endif 110 continue c c call either EPAL or ET or neither for series sum c if (al.eq.0) then c c sign of v(0) is always "+" for J0, J1 c call Euler(V,N,M,+1,sum) elseif ((al.eq.1).or.(al.eq.2)) then do 120 i=1, N c c V(0)=V(0),V(1)=V(1)+V(0),... etc c v(i)=v(i)+v(i-1) 120 continue if (al.eq.1) then call epsilon(N/2,V,sum) else call epsilon2(N/2,V,sum) endif elseif (al.ne.3) then call werror(3) endif Intgj=sum return end c_____________________________________________________________________________ c function NR c integrant of radial displacement ---> U_R=S(NR*dx) c_____________________________________________________________________________ complex*16 function NR( X,P,R,Z ) implicit none complex*16 p,k2,k2a1,k2a2,ray,ka1,ka2 real*8 X,R,Z,c1,c2,BJ1 common /com1/ c1,c2 C k2= p / dcmplx(c2,0) ka1=k2a1(x,p) ka2=k2a2(x,p) NR=((x**2) & * (-(2*(x**2) + (k2**2)) * CDexp(-ka1*z) & + 2*ka1*ka2 * CDexp(-ka2*z) ) * BJ1(x*r))/ ray(x,p) return end c_____________________________________________________________________________ c function NZ c integrant of vertical displacement ---> U_Z=S(NZ*dx) c_____________________________________________________________________________ complex*16 function NZ( X,P,R,Z ) implicit none complex*16 p,k2,k2a1,k2a2,ray,ka1,ka2 real*8 X,R,Z,c1,c2,BJ0 common /com1/ c1,c2 C C be careful. This expression differs from the Suhubi's original by C having k2a1 in the begining instead of a1 only C Take care when integrating. C k2= p / dcmplx(c2,0) ka1=k2a1(x,p) ka2=k2a2(x,p) NZ=(( x * ka1 ) & * (-(2*(x**2) + (k2**2)) * CDexp(-ka1*z) & + 2*(x**2) * CDexp(-ka2*z) ) * BJ0(x*r))/ ray(x,p) return end c_____________________________________________________________________________ C Rayleigh wave function c_____________________________________________________________________________ complex*16 function RAY( X,P ) implicit none complex*16 p,k2,k2a1,k2a2 real*8 X,c1,c2 common /com1/ c1,c2 k2= p / dcmplx(c2,0) ray=(2*(x**2)+k2**2)**2 - 4*k2a1( X,p )*k2a2( X,p )*(X**2) return end c_____________________________________________________________________________ C Derivative of Rayleigh wave function c_____________________________________________________________________________ complex*16 function DRAY( X,P ) implicit none complex*16 p,k2,k2a1,k2a2,ka1,ka2 real*8 X,c1,c2 common /com1/ c1,c2 k2= p / dcmplx(c2,0) ka1=k2a1(x,p) ka2=k2a2(x,p) dray=8*x*(2*(x**2)+(k2**2)) - ( 4*(x**3)*(ka2**2+ka1**2) & /(ka1*ka2) + 8*x*ka2*ka1 ) return end c_____________________________________________________________________________ C //// k2a1 wave sub-function c_____________________________________________________________________________ complex*16 function K2A1( X,P ) implicit none real*8 X,c1,c2 complex*16 p,k1,ka1 common /com1/ c1,c2 C k1=p/dcmplx(c1,0) ka1=CDsqrt(X**2+(k1)**2) c C Re(k2ai)>0 always, so make sure you get C the right square root c IF (dreal(ka1).gt.0) then k2a1=ka1 ELSEIF (dreal(ka1).lt.0) then k2a1= -ka1 ELSE write(*,*) 'p:',p call werror(1) write(*,*) 'p:',p stop endIF return end c_____________________________________________________________________________ C k2a2 wave sub-function c_____________________________________________________________________________ complex*16 function K2A2( X,P ) implicit none real*8 X,c1,c2 complex*16 p,k2,ka2 common /com1/ c1,c2 C k2=p/DCMPLX(c2,0) ka2=CDsqrt(X**2+(k2)**2) c C Re(k2ai)>0 always, so make sure you get C the right square root c IF (dreal(ka2).gt.0) then k2a2=ka2 ELSEIF (dreal(ka2).lt.0) then k2a2= -ka2 ELSE write(*,*) 'p:',p call werror(2) write(*,*) 'p:',p stop endIF return end c_____________________________________________________________________________ c subroutine Setwaves speed c Computes c1,c2 wave speeds (primary=1, secondary=2) c_____________________________________________________________________________ subroutine SetWaves(E,v,d) implicit none real*8 E,v,d,c1,c2 common /com1/ c1,c2 C c1=dsqrt( E*(1-v) / (d*(1+v)*(1-2*v)) ) c2=dsqrt( E / (2*d*(1+v)) ) return end c***************************************************************************** c* Pekeris's Analytical solution for displacements ur,uz * c* Valid only for Poisson's v=0.25 * c* Uses un-normalized time t * c***************************************************************************** c_____________________________________________________________________________ c Pekeris' solution for ur c_____________________________________________________________________________ real*8 function Peur(r,t) implicit none real*8 F,G,E,v,d,g2,pi,rtn,r,t,tn,c1,c2,k,ks,eik,eip common /com0/ G,E,v,d common /com1/ c1,c2 common /com11/ Pi common /com12/ F tn=c2*t/r g2=(3+dsqrt(3))/4 if (tn.le.(1d0/dsqrt(3d0))) then rtn=0 elseif (tn.lt.1d0) then k=dsqrt((3*tn**2-1)/2) ks=1d0/k rtn=tn/(16*dsqrt(6))*(6*eik(k)-18*eip(8*k**2,k) & +(6-4*dsqrt(3))*eip((20d0-12*dsqrt(3))*k**2,k) & +(6+4*dsqrt(3))*eip((20d0+12*dsqrt(3))*k**2,k)) elseif (tn.lt.dsqrt(g2)) then k=dsqrt((3*tn**2-1)/2) ks=1d0/k rtn=tn*ks/(16*dsqrt(6))*(6*eik(ks)-18*eip(8d0,ks) & +(6-4*dsqrt(3))*eip(20-12*dsqrt(3),ks) & +(6+4*dsqrt(3))*eip(20+12*dsqrt(3),ks)) else k=dsqrt((3*tn**2-1)/2) ks=1d0/k rtn=tn*ks/(16*dsqrt(6))*(6*eik(ks)-18*eip(8d0,ks) & +(6-4*dsqrt(3))*eip(20-12*dsqrt(3),ks) & +(6+4*dsqrt(3))*eip(20+12*dsqrt(3),ks)) & +pi*tn/(24*(tn**2-g2)**0.5) endif peur=3*F*Rtn/(pi**2*G*r) return end c_____________________________________________________________________________ c Pekeris' solution for uz c_____________________________________________________________________________ real*8 function Peuz(r,t) implicit none real*8 F,G,E,v,d,g2,pi,gtn,r,t,tn,c1,c2 common /com0/ G,E,v,d common /com1/ c1,c2 common /com11/ Pi common /com12/ F tn=c2*t/r g2=(3+dsqrt(3))/4 if (tn.le.(1d0/dsqrt(3d0))) then Gtn=0 elseif (tn.lt.1d0) then Gtn=-(Pi/96)*(6-dsqrt(3*dsqrt(3d0)+5)/dsqrt(g2-tn**2)+ & dsqrt(3*dsqrt(3d0)-5)/dsqrt(tn**2+dsqrt(3)/4-0.75)- & dsqrt(3)/dsqrt(tn**2-0.25) ) elseif (tn.lt.dsqrt(g2)) then Gtn=-(Pi/48)*(6-dsqrt(3*dsqrt(3)+5)/dsqrt(g2-tn**2)) else Gtn=-pi/8d0 endif peuz=3*F*gtn/(pi**2*G*r) return end c***************************************************************************** c* Complete Elliptic Integrals of 1st & 3rd kind * c* Approximations. All routines use Burlisch's algorithm * c* for cel( ) general elliptic integral * c***************************************************************************** c_____________________________________________________________________________ c function EIK(m) : Compute Elliptic Integral K(k) c uses cel c_____________________________________________________________________________ real*8 function eik(k) implicit none real*8 k,cel eik=cel(1-k**2,1d0,1d0,1d0) return end c_____________________________________________________________________________ c function EIP(m) : Compute Elliptic Integral P(n,k) c uses cel. Careful, Abramowitz-Stegun use "-n" c_____________________________________________________________________________ real*8 function eip(n,k) implicit none real*8 k,n,cel eip=cel(1-k**2,n+1d0,1d0,1d0) return end c_____________________________________________________________________________ c function cel( ) : Compute general complete Elliptic Integral c Burlisch's algorithm c_____________________________________________________________________________ real*8 FUNCTION CEL(QQC,PP,AA,BB) implicit real*8 (a-h,o-z) PARAMETER (CA=.0003d0, PIO2=1.5707963268d0) IF(QQC.EQ.0.) call werror(13) QC=ABS(QQC) A=AA B=BB P=PP E=QC EM=1. IF(P.GT.0.)THEN P=SQRT(P) B=B/P ELSE F=QC*QC Q=1.-F G=1.-P F=F-P Q=Q*(B-A*P) P=SQRT(F/G) A=(A-B)/G B=-Q/(G*G*P)+A*P ENDIF 1 F=A A=A+B/P G=E/P B=B+F*G B=B+B P=G+P G=EM EM=QC+EM IF(ABS(G-QC).GT.G*CA)THEN QC=SQRT(E) QC=QC+QC E=QC*EM GO TO 1 ENDIF CEL=PIO2*(B+A*EM)/(EM*(EM+P)) RETURN END c***************************************************************************** c* Mathematical Recipes for Bessel function of First Kind * c* All formulas and constants from Abramowitz-Stegun * c***************************************************************************** c_____________________________________________________________________________ c subroutine BZEROS c Create Array of tabulated Zeros of Bessel J0(x) and J1(x) c_____________________________________________________________________________ subroutine BZEROS implicit none DIMENSION zj(0:1,0:20) real*8 zj Integer NJZeros common /com4/ zj,NJZeros c *** j0(s) c J0(0)=1 but we list it here as a zero to aid the algorithms later. c anyway Nz(x=0) is zero so no harm is done. c zj(0, 0)= 0.0 zj(0, 1)= 2.4048255577 zj(0, 2)= 5.5200781103 zj(0, 3)= 8.6537279129 zj(0, 4)=11.7915344391 zj(0, 5)=14.9309177086 zj(0, 6)=18.0710639679 zj(0, 7)=21.2116366299 zj(0, 8)=24.3524715308 zj(0, 9)=27.4934791320 zj(0,10)=30.6346064684 zj(0,11)=33.7758202136 zj(0,12)=36.9170983537 zj(0,13)=40.0584257646 zj(0,14)=43.1997917132 zj(0,15)=46.3411883717 zj(0,16)=49.4826098974 zj(0,17)=52.6240518411 zj(0,18)=55.7655107550 zj(0,19)=58.9069839261 zj(0,20)=62.0484691902 C *** j1(s) zj(1,0)=0.0 zj(1, 1)= 3.8317059702 zj(1, 2)= 7.0155866698 zj(1, 3)=10.1734681351 zj(1, 4)=13.3236919363 zj(1, 5)=16.4706300509 zj(1, 6)=19.6158585105 zj(1, 7)=22.7600843806 zj(1, 8)=25.9036720876 zj(1, 9)=29.0468285349 zj(1,10)=32.1896799110 zj(1,11)=35.3323075501 zj(1,12)=38.4747662348 zj(1,13)=41.6170942128 zj(1,14)=44.7593189977 zj(1,15)=47.9014608872 zj(1,16)=51.0435351836 zj(1,17)=54.1855536411 zj(1,18)=57.3275254379 zj(1,19)=60.4694578453 zj(1,20)=63.61136 C ** Number of Bessel J zeros NJZeros=21 return end c_____________________________________________________________________________ c McMahon's Expansions for Large Zeros of Jn(x) c n : order of Bessel c s : consecutive number of zero c c for n=0,1 and s>=3 a minimum accuracy of 1E-6 c specifically, s:0..3 0-3 decimals c s:4..30 7-8 decimals c s:31..200 6 decimals and slowly less c the accuracy is in decimal digits, not significant. If we consider c signficant digits a constant accuracy of 7 sign.digits is expected. c The accuracy might be better but it can't be observed since the polynomial c approximations used for J0,J1 yield at least such an error. c Careful using zj(0,0)=0 which isn't a zero of J0, just something of a c computational aid for u_rz c We are using a form of the original McMahon's Expansion that is optimized c for high speed computation. Otherwise it is identical to the original c Horner's method of polynomial evaluation is also employed c_____________________________________________________________________________ real*8 function mzj(n,s) implicit none real*8 b,m,zj(0:1,0:20),Pi integer n,s,NJZeros common /com4/ zj,NJZeros common /com11/ Pi if (s.lt.0) call werror(4) if ((s.le.Njzeros-1).and.((n.eq.0).or.(n.eq.1))) then mzj=zj(n,s) else m=4*(n**2) b=(s+0.5d0*n-0.25d0)*Pi mzj=b - (m-1)*( + 1/(8*b) - (7*m-31)/(384*(b**3)) & - ((83*m-982)*m+3779)/(15360*(b**5)) & - (((6949*m-153855)*m+1585743)*m-6277237) & /(3440640*(b**7))) endif return end c_____________________________________________________________________________ C function BJ0 c Approximation for Bessel J0(x) c Horner's method of polynomial evaluation is extensively employed c_____________________________________________________________________________ real*8 function BJ0( X ) implicit none real*8 x3, f0, t0 ,x,x32 If (X.le.3d0) THEN x32=(x/3d0)**2 BJ0=1 + x32*(-2.2499997d0 + x32*(+1.2656208d0 + & x32*(-0.3163866d0 + x32*(+0.0444479d0 + & x32*(-0.0039444d0 + x32*(+0.0002100d0)))))) else x3=3d0/x f0=0.79788456D0 + x3*(-0.00000077D0 + x3*(-0.00552740D0 + & x3*(-0.00009512D0 + x3*(+0.00137237D0 + & x3*(-0.00072805D0 + x3*(+0.00014476D0)))))) t0=X-0.78539816 + x3*(-0.04166397D0 + x3*(-0.00003954D0 + & x3*(+0.00262573D0 + x3*(-0.00054125D0 + & x3*(-0.00029333D0 + x3*(+0.00013558D0)))))) BJ0=f0*DCOS(t0)/DSQRT(X) endIF return end c_____________________________________________________________________________ C function BJ1 c Approximation for Bessel J1(x) c Horner's method of polynomial evaluation is extensively employed c_____________________________________________________________________________ real*8 function BJ1( X ) implicit none real*8 x, x3, f1, t1, x32 If (X.le.3) then x32=(x/3d0)**2 BJ1=x*( 0.5d0 + x32*(-0.56249985d0 + x32*(+0.21093573d0 + & x32*(-0.03954289d0 + x32*(+0.00443319d0 + & x32*(-0.00031761d0 + x32*(+0.00001109d0))))))) else x3=3d0/x f1= 0.79788456D0 + x3*(+0.00000156D0 + x3*(+0.01659667D0 + & x3*(+0.00017105D0 + x3*(-0.00249511D0 + & x3*(+0.00113653D0 + x3*(-0.00020033D0)))))) t1=X -2.35619449D0+ x3*(+0.12499612D0 + x3*(+0.00005650D0 + & x3*(-0.00637879D0 + x3*(+0.00074348D0 + & x3*(+0.00079824D0 + x3*(-0.00029166D0)))))) BJ1=f1*DCOS(t1)/DSQRT(X) endIF return end