Editing Guide to madagascar programs (section)

=user/psava programs=

==sfawefd2d==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | acoustic time-domain FD modeling 
|-
! colspan="4" | sfawefd < Fwav.rsf vel=Fvel.rsf sou=Fsou.rsf rec=Frec.rsf wfl=Fwfl.rsf > Fdat.rsf den=Fden.rsf ompchunk=1 ompnth=0 verb=n snap=n free=n expl=n jdata=1 jsnap=nt nq1=sf_n(a1) nq2=sf_n(a2) oq1=sf_o(a1) oq2=sf_o(a2)
|-
| ''file   '' || '''den=''' ||   || 	auxiliary input file name
|-
| ''bool   '' || '''expl=n''' ||  [y/n] || 	"exploding reflector"
|-
| ''bool   '' || '''free=n''' ||  [y/n] || 	free surface flag
|-
| ''int    '' || '''jdata=1''' ||   || 
|-
| ''int    '' || '''jsnap=nt''' ||   || 	save wavefield every *jsnap* time steps
|-
| ''int    '' || '''nq1=sf_n(a1)''' ||   || 
|-
| ''int    '' || '''nq2=sf_n(a2)''' ||   || 
|-
| ''int    '' || '''ompchunk=1''' ||   || 	OpenMP data chunk size
|-
| ''int    '' || '''ompnth=0''' ||   || 	OpenMP available threads
|-
| ''float  '' || '''oq1=sf_o(a1)''' ||   || 
|-
| ''float  '' || '''oq2=sf_o(a2)''' ||   || 
|-
| ''file   '' || '''rec=''' ||   || 	auxiliary input file name
|-
| ''bool   '' || '''snap=n''' ||  [y/n] || 	wavefield snapshots flag
|-
| ''file   '' || '''sou=''' ||   || 	auxiliary input file name
|-
| ''file   '' || '''vel=''' ||   || 	auxiliary input file name
|-
| ''bool   '' || '''verb=n''' ||  [y/n] || 	verbosity flag
|-
| ''file   '' || '''wfl=''' ||   || 	auxiliary output file name
|}

Below, we will demonstrate an example of a model with nx=nz=200 and dx=dz=4m (size: 800x800m). There are two layers: the first is 100x200 samples in (z,x), and the velocity is 1500m/s; the second layer has the same dimension, and the velocity is 3000m/s. The density is set to 1 for the whole grid. A source and a receiver are co-located at x=400 and z=100. The full wavefield for the entire model aperture will be saved every 10th time step.

<ol>
  <li>Velocity model:
<syntaxhighlight lang="bash">
sfspike > Fvel.rsf mag=1500,3000 nsp=2 k1=1,101 l1=100,200 d1=4 d2=4 \
label1=z label2=x n1=200 n2=200 o1=2 o2=2 unit1=m unit2=m
</syntaxhighlight>
</li>

<li>Density model:
<syntaxhighlight lang="bash">
sfspike > Fden.rsf mag=1 nsp=1 k1=1 l1=200 d1=4 d2=4 label1=z \
label2=x n1=200 n2=200 o1=2 o2=2 unit1=m unit2=m
</syntaxhighlight>
</li>

<li>Source position (x,z):
<syntaxhighlight lang="bash">
sfspike n1=2 nsp=2 k1=1,2 mag=400,100 o1=0 o2=0 > Fsou.rsf
</syntaxhighlight>
</li>

<li>Receiver position (x,z):
<syntaxhighlight lang="bash">
sfspike n1=2 nsp=2 k1=1,2 mag=400,100 o1=0 o2=0 > Frec.rsf
</syntaxhighlight>
</li>

<li>Source wavelet:
<syntaxhighlight lang="bash">
sfspike nsp=1 n1=2000 d1=0.0005 k1=200 | sfricker1 frequency=20 |\
sftransp > Fwav.rsf
</syntaxhighlight>
</li>

<li>Creating data at specified receiver + saving full wavefield every 10th step:
<syntaxhighlight lang="bash">
sfawefd2d < Fwav.rsf vel=Fvel.rsf sou=Fsou.rsf rec=Frec.rsf wfl=Fwfl.rsf \
den=Fden.rsf > Fdat.rsf verb=y free=y expl=y snap=y dabc=y jdata=1 jsnap=10
echo 'label1=z unit1=m label2=x unit2=m' >> Fwfl.rsf
</syntaxhighlight>
</li>

<li>View the wavefield movie:
<syntaxhighlight lang="bash">
< Fwfl.rsf sfgrey gainpanel=a pclip=99 color=j scalebar=y | sfpen
</syntaxhighlight>
</li>

<li>View a wavefield snapshot:
<syntaxhighlight lang="bash">
< Fwfl.rsf sfwindow f3=80 n3=1 |\
sfgrey pclip=99 color=j title='snapshot at t=0.4s' |\
sfpen
</syntaxhighlight>
</li>

<li>View the data recorded at the receiver:
<syntaxhighlight lang="bash">
< Fdat.rsf sfwindow |\
sfgraph title='Data recorded at receiver' unit2='' label2=amplitude |\
sfpen
</syntaxhighlight>
</li>
</ol>

[[Image:sfawefd_wfld.jpg|frame|center|sfawefd wavefield screenshot]]
[[Image:sfawefd_dat.png|frame|center|sfawefd data screenshot]]

Attention: time steps that are too large can result in numerical instability.

==sfsrmig3==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | 3-D S/R migration with extended SSF
|-
! colspan="4" | sfsrmig3 slo=Fs_s.rsf sls=Fs_r.rsf < Fw_s.rsf rwf=Fw_r.rsf > Fi.rsf cig=Fc.rsf ompchunk=1 ompnth=0 verb=y eps=0.01 twoway=n nrmax=1 dtmax=0.004 pmx=0 pmy=0 tmx=0 tmy=0 vpvs=1. hsym=n nht=1 oht=0 dht=0.1 nht=1 oht=0 dht=0.1 hsym=n nhh=1 ohh=0 dhh=0.1 nha=180 oha=0 dha=2.0 nhb=180 ohb=0 dhb=2.0 itype=
|-
| ''file   '' || '''cig=''' ||   ||     auxiliary output file name
|-
| ''float  '' || '''dha=2.0''' ||   ||
|-
| ''float  '' || '''dhb=2.0''' ||   ||
|-
| ''float  '' || '''dhh=0.1''' ||   ||
|-
| ''float  '' || '''dht=0.1''' ||   ||
|-
| ''float  '' || '''dtmax=0.004''' ||   ||      max time error
|-
| ''float  '' || '''eps=0.01''' ||   ||         stability parameter
|-
| ''bool   '' || '''hsym=n''' ||  [y/n] ||
|-
| ''string '' || '''itype=''' ||   ||   imaging condition type
:o = zero lag (default)
:e = extended
:x = space-lags
:h = space-lags magnitude
:t = time-lag
|-
| ''int    '' || '''nha=180''' ||   ||
|-
| ''int    '' || '''nhb=180''' ||   ||
|-
| ''int    '' || '''nhh=1''' ||   ||
|-
| ''int    '' || '''nht=1''' ||   ||
|-
| ''int    '' || '''nrmax=1''' ||   ||  max number of refs
|-
| ''float  '' || '''oha=0''' ||   ||
|-
| ''float  '' || '''ohb=0''' ||   ||
|-
| ''float  '' || '''ohh=0''' ||   ||
|-
| ''float  '' || '''oht=0''' ||   ||
|-
| ''int    '' || '''ompchunk=1''' ||   ||       OpenMP data chunk size
|-
| ''int    '' || '''ompnth=0''' ||   ||         OpenMP available threads
|-
| ''int    '' || '''pmx=0''' ||   ||    padding on x
|-
| ''int    '' || '''pmy=0''' ||   ||    padding on y
|-
| ''file   '' || '''rwf=''' ||   ||     auxiliary input file name
|-
| ''file   '' || '''slo=''' ||   ||     auxiliary input file name
|-
| ''string '' || '''sls=''' ||   ||     auxiliary input file name
|-
| ''int    '' || '''tmx=0''' ||   ||    taper on x
|-
| ''int    '' || '''tmy=0''' ||   ||    taper on y
|-
| ''bool   '' || '''twoway=n''' ||  [y/n] ||    two-way traveltime
|-
| ''bool   '' || '''verb=y''' ||  [y/n] ||      verbosity flag
|-
| ''float  '' || '''vpvs=1.''' ||   ||  Vp/Vs ratio
|}

This program performs 3-D and 2-D shot-record (a.k.a. shot-profile) migration with an extended Split-Step Fourier (SSF) extrapolator with multiple reference velocities (hence "extended"). It takes as input a shot wavefield (<tt>stdin</tt>), receiver wavefield (<tt>rwf=</tt>) and slowness model (<tt>slo=</tt>). Outputs are an image (<tt>stdout</tt>) and a cube of Common Image Gathers (<tt>cig=</tt>). An important parameter is <tt>nrmax</tt>, the number of reference velocities. Its default value is 1, but it should be 5 or so for reasonable results. It is also good to specify nonzero taper values (<tt>tmx</tt> and, for 3-D, <tt>tmy</tt> as well). The values of padding parameters <tt>pmx</tt> and <tt>pmy</tt> are split in two by the program, i.e., if your data x-axis is 501-points long, specify pmx=11 to get a value of 512 that will result in fast Fourier Transforms. 

The program will also migrate converted-wave data if a file with the S-wave slowness model (<tt>sls=</tt>) is provided.

The <tt>vpvs</tt> parameter is only used when <tt>itype=h</tt>. Do not specify a <tt>vpvs</tt> value unless you know really well what you are doing.

===Usage example===
The commands below, slightly modified from [http://rsf.svn.sourceforge.net/viewvc/rsf/trunk/book/data/sigsbee/ptest/SConstruct?revision=3993&view=markup RSFSRC/book/data/sigsbee/ptest], show how to prepare the [https://ahay.org/RSF/book/data/sigsbee Sigsbee 2A] data and velocity for migration.

Convert input data (shots) from SEG-Y to RSF:
<syntaxhighlight lang="bash">
sfsegyread tape=sigsbee2a_nfs.segy tfile=tdata.rsf hfile=/dev/null bfile=/dev/null > ddata.rsf
</syntaxhighlight>
Convert trace headers to float (required by <tt>sfheadermath</tt>):
<syntaxhighlight lang="bash">
< tdata.rsf sfdd type=float > trchdr.rsf
</syntaxhighlight>
Shot positions:
<syntaxhighlight lang="bash">
< trchdr.rsf sfheadermath output="fldr + 10925/150" | sfwindow squeeze=y > tsi.rsf
</syntaxhighlight>
Extract offset positions from the trace header files, eliminate length-1 axis, scale, and create a header for binning (required by <tt>sfintbin</tt>):
<syntaxhighlight lang="bash">
< trchdr.rsf sfheadermath output="offset" |\
sfwindow squeeze=y |\
sfmath output="input/75" |\
sfcat axis=2 space=n tsi.rsf |\
sftransp |\
sfdd type=int > tos.rsf
</syntaxhighlight>
Binning and muting:
<syntaxhighlight lang="bash">
< ddata.rsf sfintbin head=tos.rsf xkey=0 ykey=1 |\
sfput label1=Time unit1=s d2=0.075 o2=0.0 label2=hx d3=0.150 o3=10.925 label3=sx |\
sfmutter half=false t0=1.0 v0=6.0 |\
sfput d2=0.02286 o2=0 unit2=km d3=0.04572 o3=3.32994 unit3=km > shots.rsf
</syntaxhighlight>
Keeping only 20 shots so that this 1-node job will not take forever, FFT-ing, decimating frequency slices (same as shortening the time axis), and creating y and hy axes of length 1:
<syntaxhighlight lang="bash">
< shots.rsf sfwindow n3=20 f3=10 j3=20 |\
sffft1 |\
sfwindow n1=200 min1=1 j1=3 |\
sfspray axis=3 n=1 o=0 d=1 label=hy |\
sfspray axis=5 n=1 o=0 d=1 label=sy > rfft.rsf
</syntaxhighlight>
The dimensions of the cube thus created are:
<pre>
$ sfin rfft.rsf trail=n
rfft.rsf:
    in="/var/tmp/rfft.rsf@"
    esize=8 type=complex form=native
    n1=200         d1=0.25        o1=1          label1="Frequency" unit1="Hz"
    n2=348         d2=0.02286     o2=0          label2="hx" unit2="km"
    n3=1           d3=1           o3=0          label3="hy" unit3="km"
    n4=20          d4=0.9144      o4=3.78714    label4="sx" unit4="km"
        1392000 elements 11136000 bytes
</pre>
Create the source wavelet (limited to the same frequency band as the data) and Fourier transform it:
<syntaxhighlight lang="bash">
sfspike k1=1 n1=1500 d1=0.008 |\
sfbandpass flo=15 fhi=25 |\
sffft1 |\
sfwindow n1=200 min1=1 j1=3 |\
sfput label1=freq > sfft.rsf
</syntaxhighlight>
This creates a frequency-domain wavelet:
<pre>
$ sfin sfft.rsf
sfft.rsf:
    in="/var/tmp/sfft.rsf@"
    esize=8 type=complex form=native
    n1=200         d1=0.25        o1=1          label1="freq" unit1="Hz"
        200 elements 1600 bytes
</pre>
Create "synched" source and receiver wavefields with <tt>srsyn</tt> from wavelet and data frequency slices. The receiver and shot frequency slices are "placed" at the right location and padded with zeros up to the x-axis dimension specified below.
<syntaxhighlight lang="bash">
< rfft.rsf sfsrsyn nx=1067 dx=0.02286 ox=3.05562 wav=sfft.rsf swf=swav.rsf > rwav.rsf
</syntaxhighlight>
This creates frequency slices ready for migration for both source and receiver, only axis 1 (frequency) must become axis 3, for both datasets:
<syntaxhighlight lang="bash">
< swav.rsf sftransp plane=12 | sftransp plane=23 > stra.rsf
</syntaxhighlight>
<syntaxhighlight lang="bash">
< rwav.rsf sftransp plane=12 | sftransp plane=23 > rtra.rsf
</syntaxhighlight>
This creates a surface receiver wavefield ready for input to migration. Axis 4 is the shot number. The values of axis 4 are arbitrary because each shot has been padded with zeros so that it covers the entire velocity model. Therefore, the aperture of the downward continuation for each shot will be as large as that of the whole survey.
<pre>
sfin trail=n rtra.rsf
rtra.rsf:
    in="/var/tmp/rtra.rsf@"
    esize=8 type=complex form=native
    n1=1067        d1=0.02286     o1=3.05562    label1="x" unit1="km"
    n2=1           d2=1           o2=0          label2="y" unit2="km"
    n3=200         d3=0.25        o3=1          label3="w" unit3="Hz"
    n4=20          d4=1           o4=0          label4="e" unit4="km"
        4268000 elements 34144000 bytes
</pre>
Convert the velocity model from SEG-Y to RSF, decimate, convert from feet to km, transpose, convert to slowness, and insert an additional axis:
<syntaxhighlight lang="bash">
sfsegyread tape=sigsbee2a_migvel.sgy tfile=/dev/null hfile=/dev/null bfile=/dev/null |\
sfput o1=0 d1=0.00762 label1=z unit1=km o2=3.05562 d2=0.01143 label2=x unit2=km |\
sfwindow j1=4 j2=2 |\
sfscale rscale=0.0003048 |\
sftransp |\
sfmath output="1/input" |\
sfspray axis=2 n=1 d=1 o=0 |\
sfput label2=y > slow.rsf
</syntaxhighlight>
This creates a slowness file ready for input to migration, with an x-axis identical to the x-axis of the wavefield files:
<pre>
$ sfin slow.rsf
slow.rsf:
    in="/var/tmp/slow.rsf@"
    esize=4 type=float form=native
    n1=1067        d1=0.02286     o1=3.05562    label1="x" unit1="km"
    n2=1           d2=1           o2=0          label2="y" unit2="km"
    n3=301         d3=0.03048     o3=0          label3="z" unit3="km"
        321167 elements 1284668 bytes
</pre>
Finally, the migration command (for a 4-processor machine, hence the <tt>ompnth</tt> value). We choose not to compute any image gathers (<tt>itype=o</tt>), but due to the construction of the program we still have to explicitly assign the <tt>cig</tt> tag, or else an RSF file with the name of the tag and no rsf extension will be created:
<syntaxhighlight lang="bash">
< stra.rsf sfsrmig3 nrmax=20 dtmax=5e-05 eps=0.01 verb=y ompnth=4 \
tmx=16 rwf=rtra.rsf slo=slow.rsf itype=o cig=/dev/null > img.rsf
</syntaxhighlight>
The migration of 20 shots takes approximately three hours on a 4-processor machine (1 shot=9 minutes). Without the frequency slice decimation by a factor of 3 and the depth axis decimation by a factor of 4, it would have taken twelve times as much. The resulting image has a y-axis of length 1:
<pre>
$ sfin img.rsf trail=n
img.rsf:
    in="/var/tmp/img.rsf@"
    esize=4 type=float form=native
    n1=1067        d1=0.02286     o1=3.05562    label1="x" unit1="km"
    n2=1           d2=1           o2=0          label2="y" unit2="km"
    n3=301         d3=0.03048     o3=0          label3="z" unit3="km"
        321167 elements 1284668 bytes
</pre>
To properly visualize the image, we need to eliminate the axis of length 1, then transpose the x and z axes to their natural position:
<syntaxhighlight lang="bash">
<img.rsf sfwindow squeeze=y | sftransp | sfgrey > img.vpl
</syntaxhighlight>