! ! 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