!
! program written by Dimitri Komatitsch (University of Pau, France)
! in 2000 while working at Harvard University (Cambridge, USA)
! based on information provided by Jose Carcione (OGS, Trieste, Italy)
!
! Analytical solution for a 3D TI medium along the symmetry axis
! taken from Carcione et al., Geophysics, vol. 57, p. 1606 (1992)
!
! Computes the component of the displacement vector recorded along
! the symmetry axis caused by a point force acting on that same axis,
! convolved with a Ricker time signal.
!
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
! THE FOLLOWING PARAMETERS CAN BE CHANGED BY THE USER:
! - The material properties: line 32
! - The offset source-receiver: line 35
! - The seismograms time length and number of timesteps: line 39
! - The source type:  line 42 (if changed, the user will have to provide a new function similarly as the one in line 162)
!
! In case of using the default Ricker time signal, one can change the center frequency and time delay in lines: 171-2)
! Also the amplitudes can be scaled by changing line: 175
!
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  program TI_analytical_3D

  implicit none

! parameters of the TI medium (Mesaverde Clay shale in this example)
! Any c_ij values can be entered as long as the material remains transversely isotropic
  double precision, parameter :: c11 = 66.6d9, c33 = 39.9d9, c55 = 10.9d9, c12 = 19.7d9, c13 = 39.4d9, c44 = c55, rho = 2590.d0

! source-receiver distance
  double precision, parameter :: z = 728.856213417182d0

! total duration of the seismograms and total number of time steps
  double precision, parameter :: total_duration = 1.5d0
  integer, parameter :: nstep = 250000

! derivative of a Ricker wavelet
  double precision, external :: derivricker

! local variables
  integer it,i,j,k,k0,k1,i1,i2,ival
  double precision alpha,beta,gamma,zbar,V_s,t,ts,tp,t1,Dz
  double precision zbar1,hz,uz,vr,vp,t2,t3,t4,sumconv,deltat
  double precision, dimension(-nstep:nstep-1) :: s
  double precision, dimension(0:nstep-1) :: r,time_v,storesumconv

! beginning of the main program
  V_s = dsqrt(c55/rho)
  ts = z / V_s
  tp = z / dsqrt(c33/rho)
  alpha=c33/c44
  beta=c11/c44
  gamma=1.+alpha*beta-(c13/c44+1.)**2
  vr=dsqrt(c44/rho)
  vp=vr*dsqrt(alpha)
  t1=gamma*(beta+1.)
  t2=2.*beta*(alpha+1.)
  t3=1.+alpha*beta-gamma
  t4=alpha+beta-gamma
  zbar1=dsqrt(t1-t2+2.*dsqrt(beta*t3*t4))/(beta-1.)
  t1 = ts / zbar1

! this analytical solution is valid only for a subset of materials
  if(gamma >= beta + 1.d0) stop 'invalid solution 1'
  if(gamma**2 - 4.d0*alpha*beta >= 0.d0) stop 'invalid solution 2'

! duration of a time step
  deltat = total_duration / dble(nstep-1)

! loop on all the time steps
  do it = 0,nstep-1

! compute real time
  t = deltat * dble(it)

! if we are outside of the time interval, do not compute anything
  if(t <= tp) then
    uz = 0.d0
    goto 20
  else if(t > t1) then
    uz = 1.d0
    goto 20
  endif

  zbar = z / (V_s*t)

  Dz = (gamma-(beta+1.d0)*zbar**2)**2 - 4.d0*beta*(alpha-zbar**2)*(1.d0-zbar**2)

  hz = 0.5d0 - (2.d0*(1.d0-zbar**2) - gamma + (beta+1.d0)*zbar**2)/(2.d0*dsqrt(Dz))

  if(t <= tp) then
    uz = 0.d0
  else if(t <= ts) then
    uz = hz
  else if(t <= t1) then
    uz = 2.d0 * hz
  else
    uz = 1.d0
  endif

! store the result for convolution
 20 time_v(it) = t
  r(it) = uz

  enddo

! fill up the source array
  do i=0,nstep-1
    s(i) = derivricker(time_v(i))
  enddo

! bounds between which the Ricker is non-zero
  t1 = 0.d0
  t2 = 0.17d0
  i1 = nint(t1/deltat)
  i2 = nint(t2/deltat)

! zero padding
  do i=-nstep,-1
    s(i) = 0.d0
  enddo

! perform the convolution in time
  do j=0,nstep-1

  sumconv = 0.d0

! bounds for the sum over k
  k0 = 0
  k1 = nstep-1

! limit the sum to the non-zero part of the source time function
  ival = j-i1
  if(ival < k1 .and. ival >= 0) k1 = ival
  ival = j-i2
  if(ival > k0 .and. ival <= nstep-1) k0 = ival

  do k=k0,k1
    sumconv = sumconv + s(j-k) * r(k)
  enddo

  storesumconv(j) = sumconv

  enddo

! save the convolved solution
  open(unit=10,file='output.dat',status='unknown')
  do j=0,nstep-1
    write(10,*) time_v(j),storesumconv(j)
  enddo
  close(10)

  end program TI_analytical_3D

!
! Derivative of Ricker for the convolution
!
  double precision function derivricker(t)

  implicit none

  double precision, parameter :: pi = 3.141592653589793d0

  double precision t,f0,t0,factor,a

! parameters of the source
  f0 = 16.d0
  t0 = 0.07d0

! constant scaling factor (can be changed to anything because the equation is linear)
  factor = 10.d0

! derivative of a Ricker
  a = pi*pi*f0*f0
  derivricker = - factor * 2.d0 * a * (t-t0) * (-3.d0+2.d0*a*(t-t0)**2)*exp(-a*(t-t0)**2)

  end function derivricker

