# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!!$
!!$=head1 NAME
!!$
!!$ ShotMig - 2/3D Shot profile migration (Green's function modeling)
!!$ with nearest-neighbor interpolated PSPI. WEIlib based.
!!$ Written by Brad Artman.
!!$
!!$ MAKE SURE YOUR SHOT LOCATION IS NOT IN THE TAPER AREA!!
!!$
!!$ Performs converted wave imaging, direct passive migration,
!!$ and forward scattering imaging.
!!$ Makes Image-space SRME model and illumination normalization fields.
!!$ Imaging conditions sum over all frequencies and shot locations.
!!$
!!$ Some output required. What that might be is flexible.
!!$ OMP parallelized. Ready for python (man Parallel) parallel run.
!!$
!!$ Prepare data (as described below) with:
!!$ Transf shotmig=y < time.H > freq.H !!$ f_min=minHz f_min1=polelow f_max1=polehi f_max=maxHz

!!$ If 2D, use param: is2d=y (see Transf self-doc).
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ShotMig D=shotcube.H D_d=sourcefunc.H par=parfile.P !!$ S(V)=slo(vel).H R=image.H nref=[1]

!!$
!!$=head1 Files
!!$
!!$=over 4
!!$
!!$=item D=shotcube.H (upward propogated wavefield U)
!!$
!!$ sep77 hypercube where n1=gx n2=gy n3=freq n4=sx n5=sy
!!$ Do not use negative deltas. Instead, use large negative
!!$ o's and positive d's. Origin for sx(sy) very important
!!$ even if only migrating one shot.
!!$
!!$ If D is not a sep77, can include a header file.
!!$ Dimensions of D are then
!!$ n1=frequency, n2=ntraces, n3=shots
!!$
!!$=item HD=xy.header.H (optional)
!!$
!!$ sep77 hypercube where n1=2, n2=ntraces, n3=shots.
!!$ o3,d3 of D must=1. Must have source headers (file HD_d).
!!$ See parameter xyreverse.
!!$
!!$=item D_d=sourcefunc.H (downward propogated wavefield D)
!!$
!!$ sep77 hypercube where n1=1 n2=1 n3=freq (n4=1 n5=1)
!!$ (o3,n3,d3 must be identical to those for shotcube.H)
!!$ n1,n2 can be size of a shotgather to make areal sources.
!!$
!!$=item HD_d=xy.header.shots.H (optional)
!!$
!!$ sep77 hypercube where n1=2, and n2=shots from D.
!!$ see parameter xyreverse.
!!$
!!$=item S(V)=slo(vel).H
!!$
!!$ sep77 hypercube where n1=x n2=y n3=z
!!$ S=slo.H automatically tells program units are time/distance
!!$ V=vel.H automatically tells program units are distance/time
!!$ Velocity file must be sampled as gx,gy axes, and be at least
!!$ as large as the entire output domain (will not pad for you).
!!$
!!$=item S(V)_R=rcvr_slowness(vel).H (optional, must have S(V)_S)
!!$
!!$ Receiver side slowness(velocity) model for converted wave imaging.
!!$
!!$=item S(V)_S=src_slowness(vel).H (optional, must have S(V)_R)
!!$
!!$ Source side slowness(velocity) model for converted wave imaging.
!!$
!!$=item R=image.H migration result with imaging condition UD*
!!$
!!$ sep77 hypercube where n1=x n2=y n3=halfoffset_x n4=halfoffset_y n5=z
!!$ add=y sums previous result and dictates size of image domain.
!!$
!!$=item M=mults.H image space multiple model with SRME imaging condition UU
!!$
!!$ sep77 hypercube where n1=x n2=y n3=halfoffset_x n4=halfoffset_y n5=z
!!$ add=y sums previous result and dictates size of image domain.
!!$
!!$=item N_d=normd.H image space illumination model DD*
!!$
!!$ sep77 hypercube where n1=x n2=y n3=halfoffset_x n4=halfoffset_y n5=z
!!$ add=y sums previous result and dictates size of image domain.
!!$
!!$=item N_u=normu.H image space illumination model UU*
!!$
!!$ sep77 hypercube where n1=x n2=y n3=halfoffset_x n4=halfoffset_y n5=z
!!$ add=y sums previous result and dictates size of image domain.
!!$
!!$=item W=wavfield_up.H Greens function from continuation of D=
!!$
!!$ sep77 hypercube where n1=gx n2=gy n3=freq n4=z
!!$ add=y does not effect wavefields.
!!$
!!$=item W_d=wavfield_down.H Greens function from continuation of D_d=
!!$
!!$ sep77 hypercube where n1=gx n2=gy n3=freq n4=z
!!$ add=y does not effect wavefields.
!!$
!!$=head1 Green's function modeling
!!$
!!$ ShotMig D_d=source.H S=slow.H W_d=wfld_d.H nref=[1] !!$ nx= ny= ox= oy= dx= dy= nsx= nsy= osx= osy=

!!$ produce only Greens functions w/ size of domain specified
!!$ by n,o,d parameters included in command
!!$
!!$=head1 Source location search
!!$
!!$ ShotMig D=data.H V=vel.H W=wfld.H nref=[1] !!$ nx= ny= ox= oy= dx= dy= nsx= nsy= osx= osy=

!!$ produce only acausally propagated data wavefields
!!$ to locate a subsurface source function.
!!$ Size of domain specified by n,o,d parameters
!!$ included in command
!!$
!!$=head1 Parallel implimentation
!!$
!!$ Parallelized for python (and OMP) distribution.
!!$ When splitting shot axes, important parameters
!!$ are total shots from survey nx(y)_shots= and
!!$ origin of total shots from survey ox(y)_shots= .
!!$
!!$ IMPORTANT: environmental variable controlling
!!$ stacksize for OMP operations needs to be huge
!!$ for some reason. setenv KMP_STACKSIZE 600000000
!!$ in your ~/.cshrc should do the trick. Regardless,
!!$ try to keep memory consumption below
!!$ MEMORYUSAGE*OMP_NUM_THREADS=1 (don't know why)
!!$ by using the nws and nzs parameters.
!!$
!!$ nohup Parallel command=ShotMig par=par.P !!$ nx_shots=`Get n4 parform=n < shotcube.H` !!$ ox_shots=`Get o4 parform=n < shotcube.H` !!$ ny_shots=`Get n5 parform=n < shotcube.H` !!$ ox_shots=`Get o5 parform=n < shotcube.H` !!$ files="shotcube.H,source.H,vel.H,image.H" !!$ tags="D,D_d,V,R" axis=4,3,4,6 !!$ usage="INPUT,INPUT,INPUT,OUTPUT" !!$ file_type="DISTRIBUTE,COPY,COPY,COPY" !!$ pverb=1 nblock=your_problem_not_mine
# 147 "<stdin>"
!!$
!!$=head1 Parameters
!!$ WARNING: SOME DEFAULTS CHANGE WITH COMBINATIONS OF PARAMETERS
!!$ RUN THE PROGRAM AND READ THE OUTPUT CAREFULLY!!!
!!$
!!$=item io_file=[SCREEN] NONE (character)
!!$
!!$ Optional file to write all verbosity.
!!$
!!$=item badsign=[n]/y Reverses propagation causality
!!$
!!$ So, you FTd a huge data set with the wrong sign. Use this parameter
!!$ instead of remaking. Always be sure D and D_d used same transform sign.
!!$
!!$=item add=[n]/y adds previous output
!!$
!!$ Does not work for first run in "Parallel",
!!$ but anything goes for serial runs. Parameters
!!$ taken from image(multiple model) file.
!!$
!!$=item adj=[y]/n migration vs. modeling
!!$
!!$ Does not perform migration. Only Greens function modeling
!!$ (Wavefield). If datum=y, uses D instead of D_d.
!!$
!!$=item nref=[1] (integer)
!!$
!!$ Number of Lloyds algorithm selected reference velocities for PSPI
!!$
!!$=item watercolumn=[-99999.99] (real)
!!$
!!$ Depth of water column for datuming with wtr_steps.
!!$ Defines firstz for image at new datum. Velocity still
!!$ extracted from V, so can be any homogeneous first layer.
!!$
!!$=item wtr_steps=[2] (integer)
!!$
!!$ Number of phaseshift steps through water column
!!$
!!$=item tshift=[0.] (real)
!!$
!!$ Time shift applied to source function to correspond to early
!!$ time windowed from the data with no events (i.e. water column).
!!$
!!$=item minx=[calculate] (real)
!!$
!!$ Minimum x value for image output (R,M,N)
!!$
!!$=item miny=[calculated] (real)
!!$
!!$ Minimum y value for image output (R,M,N)
!!$
!!$=item maxx=[calculate] (real)
!!$
!!$ Maximum x value for image output (R,M,N)
!!$
!!$=item maxy=[calculate] (real)
!!$
!!$ Maximum y value for image output (R,M,N)
!!$
!!$=item nx(y)_shots=[n4(n5)] (integer)
!!$
!!$ Number of shots in total survey- only important
!!$ if splitting shot axis for parallelization.
!!$
!!$=item ox(y)_shots=[o4(o5)] (real)
!!$
!!$ Origin of shots in total survey- only important
!!$ if splitting shot axis for parallelization.
!!$
!!$=item nhx(y)=[1] (integer)
!!$
!!$ Number subsurface offsets to calculate
!!$ I(x,h)=U(x-h/2)D*(x+h/2)
!!$
!!$=item jhx(y)=[1] (integer)
!!$
!!$ Set delta h = delta x * jhx
!!$
!!$=item item dhx(y)=[dx*jhx] (real)
!!$
!!$ Alternative to jhx. mod(dhx,dx)==0
!!$
!!$=item item hneg(y)=[false] (logical)
!!$
!!$ Make all offsets negative (oh=-nh/dh)
!!$
!!$=item item hpos(y)=[false] (logical)
!!$
!!$ Make all offsets posative (oh=0.)
!!$
!!$=item item ohx(y)=[calculated] (real)
!!$
!!$ Specify origin for offset axes
!!$
!!$=item marine=[y]/n OUTDATED
!!$
!!$ Geometry flag indicating y sets : offset=y and rollalong=y
!!$ n sets : offset=n and rollalong=n
!!$
!!$=item offset=[y]/n (logical)
!!$
!!$ Receiver coordinates defined by distance from source location.
!!$
!!$=item rollalong=[y]/n (logical)
!!$
!!$ Rolling acquisition for each shot. False means working space
!!$ is the same for every shot.
!!$
!!$=item decon=[n]/y (logical)
!!$
!!$ Divide R,M,N_d,N_u imaging conditions by DD*
!!$
!!$=item eps=[0.001] (real)
!!$
!!$ epsilon value check for division by DD*
!!$
!!$=item passive=[n]/y (logical)
!!$
!!$ Direct passive migration- no D_d file. D_d=conj(D)
!!$ Also true when called with only D= and R= files.
!!$ Assumes n4(*n5) of D is "event(time-window)"-axis.
!!$
!!$=item datum=[n]/y (logical)
!!$
!!$ Uses only D(ata) and makes only W(avefield) to back-propagate
!!$ events (locating earthquakes or datuming).
!!$
!!$=item autonzs=[y]/n (logical)
!!$
!!$ Holds all depths in memory for each shot panel migration/modeling.
!!$ False if nws not equal to depth axis of slowness(velocity) model.
!!$
!!$=item autonws=[y]/n (logical)
!!$
!!$ Holds all frequencies in memory for each shot panel migration/modeling.
!!$ False if nws not equal to frequency axis of input.
!!$
!!$=item nzs=[n3 of vel(slo).H] (integer)
!!$
!!$ Number of depths per memory block
!!$
!!$=item nws=[n3 of D_d= (D= if passive)] (integer)
!!$
!!$ Number of frequencies per memory block
!!$
!!$=item fwdscat=[n]/y (logical)
!!$
!!$ Forward scatterd imaging with converted modes.
!!$ Changes propagation direction for source.H
!!$ See SEP113 page 157 "Imaging with buried sources"
!!$
!!$=item scale=[1.] (real)
!!$
!!$ data=scale*data
!!$
!!$=item taper=[y]/n (logical)
!!$
!!$ Cosine taper edges of wavefield domain
!!$ MAKE SURE YOUR SHOT LOCATION IS NOT IN THE TAPER AREA!!
!!$
!!$=item taper_pcent=[7/naxis] (real)
!!$
!!$ Percent (on each side) of the wavefield to apply cosine taper
!!$
!!$=item minx(y)1=[7] (integer)
!!$
!!$ Trace to end taper (up) from front edge
!!$
!!$=item maxx(y)1=[7] (integer)
!!$
!!$ Trace to begin taper (down) to far edge
!!$
!!$=item smooth=[n]/y (logical)
!!$
!!$ First smooth DD* (over x-y plane) before using for division.
!!$ Uses triangle smoothing (2*rect+1 long)
!!$
!!$=item rectx=[2] (integer)
!!$
!!$ Specifies x-length of box for smoothing
!!$
!!$=item recty=[rectx] (integer)
!!$
!!$ Specifies y-length of box for smoothing
!!$
!!$=item antialias=[n]/y (logical)
!!$
!!$ Cosine taper wavenumber components from wavefields
!!$ before extrapolation. Defined, as a function of
!!$ frequency from ray-parameter limit parameters.
!!$ Only filters before first raw data. Remove
!!$ conditional line 116 in wem.F90 for all z levels.
!!$
!!$=item pmin_x(y) (real)
!!$
!!$ Maximum negative ray-parameter to allow.
!!$ Must be negative real number.
!!$
!!$=item pmax_x(y) (real)
!!$
!!$ Maximum posative ray-parameter to allow.
!!$ Must be posative real number.
!!$
!!$=item image_real=[y]/n (logical)
!!$
!!$ Output real (complex) values for R=
!!$
!!$=item firstz=[1] (integer)
!!$
!!$ First depth level to write into R,M,N_d,N_u
!!$
!!$=item firstzW=[1] (integer)
!!$
!!$ First depth level to write into W_d,W
!!$
!!$=item xyreverse=[n] (logical)
!!$
!!$ If 1-axis of headerfile is (y,x) instead of (x,y)
!!$
!!$=item debug=[n] (logical)
!!$
!!$ Set for memory allocation sums
!!$
!!$=item maxang=[50.] (real)
!!$
!!$ Maximum angle through which to propogate. SSR is only good to about 45.
!!$

!!$=back
!!$
!!$=head1 KEYWORDS
!!$
!!$ subsurface offset,common image gather, ODCIG, wave-equation
!!$ migration, Greens function, IS-SRME, multiples, passive,
!!$ direct migration
!!$
!!$=head1 CATEGORY
!!$
!!$B<seis/image>
!!$
!!$=cut
!!$
!!$ Brad Artman (brad@sep.stanford.edu)
!!$ 06/6/06
!!$
!!$
!>
!%%

program Shotprofile
  use sep
  use SPmig_mod
  use wem_mod
  use image_mod
  use ssf_mod
  use sll_mod


  logical::verb
  integer::stat,ierr,iverb,n,i
  integer:: time0(8), time1(8)
  real :: start, finish

  !----------------------------------------------------------------

  start=0.;finish=0.
  call date_and_time(values=time0)
  start_time = time0 (5) * 3600 + time0 (6) * 60 &
  + time0 (7) + 0.001 * time0 (8)

  call sep_init("/homes/sep/claudio/research.ln/Others/ShotMig/ShotMig.F90")
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  !----------------------------------------------------------------
  call SPmig_init(SLin=sll_init,&
     WCin=kxmigN_init,&
     FKin=wem_init,&
     FXin=ssf_init,&
     IGin=ig_init)

  call date_and_time(values=time1)
  finish_time = time1 (5) * 3600 + time1 (6) * 60 &
  + time1 (7) + 0.001 * time1 (8)
  if(run%verb)write(run%io,*) '        * '
  if(run%verb)write(run%io,'(8x,a,f12.3,a)')' * Main initiallization: ', &
  finish_time - start_time,' (s)'


  stat=SPmig(down%adj,down%add,&
   SLop=sll, &
   WCop=kxmigN_SSR, &
   FKop=wem,&
   FXop=ssf,&
   IGop=offIgop,&
   MIGop=multIGop,&
   NDIGop=normIGop)


  call date_and_time(values=time1)
  finish_time = time1 (5) * 3600 + time1 (6) * 60 &
  + time1 (7) + 0.001 * time1 (8)

  if(run%verb)write(run%io,'(8x,a,f10.2,a)')' Total elapsed time:', &
  (finish_time - start_time)/60.,' (m)'

  call sep_end_prog()
  write(run%io,*)'Gracias a Dios-  TERMINE !'
  call exit(0)

end program Shotprofile
