Editing Guide to madagascar programs (section)

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