c
      program dynamic2

      IMPLICIT REAL*8(A-H,O-Z)
      complex*16 im,ss,pp,d1,d2,dc1,dc2,dsum1,dsum2
      real*8 mu
c  
      vs=3000.d0
      vp=5000.d0
      a=vs*vs/(vp*vp)
      mu=0.5d0*(1.d0-2.d0*a)/(1.d0-a)
      write(*,1000)vs,vp,mu   
 1000 format('vs,vp,mu   =',3e12.4)   
c
      slip =1.d0
      write(*,*)'Enter Zs'
      read(*,*)Zs
c
      R1=Zs*30.d0
      R2=10000.d0
      m1=1000
      m2=2000
      if(R1.ge.R2)then
        R1=R2
        m2=0
      end if
      write(*,1010)slip,R1,R2
 1010 format('slip,R1,R2 =',3e12.4)
      write(*,1020)m1,m2
 1020 format('m1,m2      =',2i12)
      write(*,*)
c
      pi=dacos(-1.d0)
c      write(*,*)pi
   10 write(*,*)'Enter frequency (Hz)'
      read(*,*)freq
      waves=vs/freq
      wavep=vp/freq
      write(*,*)'Wave Lengths for S and P=',waves,wavep
      om=freq*2.d0*pi

      im=dCMPLX(0.d0,1.d0)      
      ss=+im*om/vs
      pp=+im*om/vp
c
      c3=-1.d0/(8.d0*pi*(1.d0-mu))
      c4=1.d0-2.d0*mu
c
      dR=R1/m1
      write(*,*)'dR1 =',dR
c
      open(9,file='d2_d',status='unknown')
      write(9,2000)m1,m2,R1,R2
 2000 format(2i12,2e12.4)
c
      rd=0.d0
      u1=0.d0
      d1=(0.d0,0.d0)
c      write(*,9999)rd,-pi*c3*Zs*u1,d1
      write(9,2010)rd,-pi*c3*Zs*u1,d1
 2010 format(4e12.4)
      sum1=0.d0
      dsum1=(0.d0,0.d0)
      do j=1,m1
        rd=dR*j
        rr=sqrt(rd*rd+Zs*Zs)
        rr2=rd*rd/(rr*rr)
        call dconst(rr,ss,pp,dc1,dc2)
        r3=rr*rr*rr
        u2=(2.d0*C4 +3.d0*rr2)*rd/r3
        d2=(2.d0*dc1+dc2*rr2)*rd/r3
        sum1=sum1+0.5d0*(u1+u2)
        dsum1=dsum1+0.5d0*(d1+d2)
        u1=u2
        d1=d2
c        write(*,9999)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
 9999   format('rd,u1,d1=',f12.4,3e12.4)
        write(9,2010)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
      end do
      sum1=-pi*c3*Zs*sum1*dR
      dsum1=-Zs*dsum1*dR/4.d0
c
      if(R1.ge.R2)goto 20
      if(m2.eq.0)goto 20
      dR=(R2-R1)/m2
      write(*,*)'m2,dR',m2,dR
c
      rd=R1
      rr=sqrt(rd*rd+Zs*Zs)
      rr2=rd*rd/(rr*rr)
      call dconst(rr,ss,pp,dc1,dc2)
      r3=rr*rr*rr
      u1=(2.d0*C4 +3.d0*rr2)*rd/r3
      d1=(2.d0*dc1+dc2*rr2)*rd/r3
c      write(*,9999)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
      write(9,2010)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
      sum2=0.d0
      dsum2=(0.d0,0.d0)
      do j=1,m2
        rd=R1+dR*j
        rr=sqrt(rd*rd+Zs*Zs)
        rr2=rd*rd/(rr*rr)
        call dconst(rr,ss,pp,dc1,dc2)
        r3=rr*rr*rr
        u2=(2.d0*C4 +3.d0*rr2)*rd/r3
        d2=(2.d0*dc1+dc2*rr2)*rd/r3
        sum2=sum2+0.5d0*(u1+u2)
        dsum2=dsum2+0.5d0*(d1+d2)
        u1=u2
        d1=d2
c        write(*,9999)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
        write(9,2010)rd,-pi*c3*Zs*u1,-Zs*d1/4.d0
      end do
      sum2=-pi*c3*Zs*sum2*dR
      dsum2=-Zs*dsum2*dR/4.d0
c
 20   write(9,2010)Zs,sum1,sum2,sum1+sum2
      write(9,2010)dsum1,dsum2
      close(9)
c
      write(*,*)
      write(*,3000)sum1,sum2,sum1+sum2
      write(*,3010)dsum1,dsum2,dsum1+dsum2
 3000 format('Static : sum1,sum2,sum1+2='/,3(e12.4,12x))
 3010 format('Dynamic: sum1,sum2,sum1+2='/,6e12.4)
c
c      goto 10     
c
      stop
      end
c
c
c
      subroutine dconst(rr,ss,pp,dc1,dc2)

      IMPLICIT REAL*8(A-H,O-Z)
      complex*16 ss,sr,sr2,sr3,pp,pr,pr2,pr3,ps2,esr,epr
      complex*16 kai,kair,fair,dc1,dc2
c
      sr=ss*rr
      sr2=sr*sr
      sr3=sr*sr2
      pr=pp*rr
      pr2=pr*pr
      pr3=pr*pr2
      ps2=pp*pp/(ss*ss)
c
      esr=cdexp(sr)
      epr=cdexp(pr)
c
      kai=(3.d0/sr2-3.d0/sr+1.d0)*esr
      kai=kai-ps2*(3.d0/pr2-3.d0/pr+1.d0)*epr
      kair=(-9.d0/sr3+9.d0/sr2-4.d0/sr+1.d0)*esr*sr
      kair=kair-ps2*(-9.d0/pr3+9.d0/pr2-4.d0/pr+1.d0)*epr*pr
      fair=(-3.d0/sr3+3.d0/sr2-2.d0/sr+1.d0)*esr*sr
      fair=fair-ps2*(-3.d0/pr3+3.d0/pr2-1.d0/pr)*epr*pr
c
      dc1=fair-kai
      dc2=2.d0*(2.d0*kai-kair)
c
c      write(*,*)'r*kai',dreal(ps2-1.d0)*0.5d0,kai
c      write(*,*)'kair ',dreal(1.d0-ps2)*0.5d0,kair
c      write(*,*)'fair ',-dreal(ps2+1.d0)*0.5d0,fair
c      write(*,*)'dc1',-dreal(ps2),dc1
c      write(*,*)'dc2',3.d0*dreal(ps2-1.d0),dc2
c
      return
      end
