#############################################################################
###################### COMMON DEFINITIONS - DO NOT ALTER ####################
#############################################################################
from rsf.proj import *
#############################################################################
###################### END COMMON DEFINITIONS ###############################
#############################################################################

######################## LOCAL DEFINITIONS ##################################

######### fetch list - non-reproducible data fetched from web archive.
# the format for this dictionary is 
#    filename : [<subdir of "data" directory in archive>, <archive URL>]
# Example: 
# fetches = {
#            'line_fix.su' : ['marmousi', 'http://www.trip.caam.rice.edu'],
#            'velocity.HH' : ['marmousi', 'http://www.trip.caam.rice.edu'],     
#           }
fetches = {
           'vel_4layer.HH': ['obs', 'http://www.trip.caam.rice.edu']
          }
for file in fetches.keys():
    Fetch(file,fetches[file][0],server=fetches[file][1])

######### non-IWAVE flows - include these in standard Madagascar form
Flow('null',None,'spike n1=651 d1=0.008 n2=400 o2=1 d2=1 mag=0')
Flow('head','null','window n1=1')

Flow('tracl','head','math output=x1               | dd type=int')
Flow('selev','head','math output=-1875            | dd type=int')
Flow('gelev','head','math output=-6               | dd type=int')
Flow('sx',   'head','math output=12000            | dd type=int')
Flow('gx',   'head','math output="5000+25*(x1-1)" | dd type=int')
Flow('offset','sx gx','add scale=-1,1 ${SOURCES[1]}')

Flow('tnull','null tracl selev gelev sx gx offset',
     '''
     segyheader 
     tracl=${SOURCES[1]} 
     selev=${SOURCES[2]}  
     gelev=${SOURCES[3]}  
     sx=${SOURCES[4]}  
     gx=${SOURCES[5]}  
     offset=${SOURCES[6]} 
     ''')

Flow('hdr12000.su','null tnull','suwrite tfile=${SOURCES[1]} endian=0')

Flow('null12',None,'spike n1=651 d1=0.008 n2=400 o2=1 d2=1 mag=0 n3=12 o3=0 d3=1')
Flow('head12','null12','window n1=1')

Flow('tracl12','head12','math output="x1+400*x2"      | dd type=int')
Flow('selev12','head12','math output=-1875            | dd type=int')
Flow('gelev12','head12','math output=-6               | dd type=int')
Flow('sx12',   'head12','math output="8000+400*x2"    | dd type=int')
Flow('gx12',   'head12','math output="5000+25*(x1-1)" | dd type=int')
Flow('offset12','sx12 gx12','add scale=-1,1 ${SOURCES[1]}')

Flow('tnull12','null12 tracl12 selev12 gelev12 sx12 gx12 offset12',
     '''
     segyheader 
     tracl=${SOURCES[1]} 
     selev=${SOURCES[2]}  
     gelev=${SOURCES[3]}  
     sx=${SOURCES[4]}  
     gx=${SOURCES[5]}  
     offset=${SOURCES[6]} 
     ''')

Flow('hdr8-12km.su','null12 tnull12','suwrite tfile=${SOURCES[1]} endian=0')

# sfadd writes output to stdout, so it must be left active - 
# however there is not pipe input to this command, so stdin=0
Flow('vel_4layer.rsf', 'vel_4layer.HH', 'dd form=native');
Flow('csq_4layer','vel_4layer.rsf',
     'mul ${SOURCES[0]} | put data_type=csq ')

# smooth csq_4layer to make migration velocity model
Flow('csq_4layersm','csq_4layer','smooth rect1=50 rect2=50 repeat=4')
Flow('csq_4layersm2','csq_4layer','smooth rect1=5 rect2=5 repeat=4')
Flow('dcsq_4layer','csq_4layer csq_4layersm2','add scale=1,-1 ${SOURCES[1]}')

# create base wavelet (just time series, without source position
# information) via suwaveform

import math

fpeak=5
t0 = 1.0/fpeak
s  = 1.0 /(math.sqrt(2)*math.pi*fpeak)

# Gaussian derivative wavelet (type=gaussd in suwaveform)
Flow('wavelet',None,
     '''
     math n1=101 d1=0.004 
     output="%g*(x1-%g)*exp(-%g*(x1-%g)^2)" | put o1=-0.1
     ''' % (-1.e6/(s*s*s*math.sqrt(2*math.pi)),t0,1/(2*s*s),t0))
Flow('twavelet','wavelet','segyheader')
Flow('wavelet_base.su','wavelet twavelet','suwrite endian=0 tfile=${SOURCES[1]}')

# add source coordinates from hdrfile to source AND receiver 
# coordinates from wavelet to create "dressed" wavelet for array
# source option in iwave. Note that iwave reads a source GATHER by
# detecting new source coordinates (sx, sy, selev) but assigns source
# trace GRID POSITIONS in the array by receiver coordinates (gx, gy, 
# gelev). The towed array app sets these coordinates up so that each 
# shot may have an array of sources, with the source traces in the 
# same position relative to the data source coordinates - hence 
# "towed_array".

# use naming convention: time series stored in wavelet_base, 
# headers for experiment foo stored in hdrfoo.su, wavelet in
# waveletfoo.su

for foo in ['12000', '8-12km']:
    Flow('wavelet' + foo + '.su', ['wavelet_base.su', 'hdr' + foo + '.su'],
         '''
         towed_array
         data=${SOURCES[1]} 
         src=${SOURCES[0]} 
         towed=${TARGETS[0]}
         ''',
         stdin=0, stdout=0)

# initial data - near-impulse in u0, u1=0 - very wierd, should be fixed!!
Flow('init',None,'makevel n1=416 n2=800 d1=25 d2=25 dlens=200 tlens=200 v000=0 vlens=100 x1lens=1875 x2lens=12000 | sfput id1=0 id2=1 id3=2 dim=2 gdim=3 n3=1 o3=0 d3=1')

# movie-init output file
Flow('movieinit','init csq_4layer',
     '''
     makevel n1=416 n2=800 n3=33 d1=25 d2=25 d3=160 v000=0  | 
     put id1=0 id2=1 id3=2 dim=2 gdim=3 > $TARGET &&                            
     acd deriv=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} init=${SOURCES[0]} sampord=1 movie=$TARGET
     ''',stdin=0,stdout=-1,workdir='movieinit.work')

# movie-src output file
Flow('moviesrc','wavelet12000.su csq_4layer',
     '''
     makevel n1=416 n2=800 n3=33 d1=25 d2=25 d3=160 v000=0  | 
     put id1=0 id2=1 id3=2 dim=2 gdim=3 > $TARGET &&
     acd deriv=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} source=${SOURCES[0]} sampord=1 movie=$TARGET
     ''',stdin=0,stdout=-1,workdir='moviesrc.work')

# trace-src output file
Flow('shot12000.su','wavelet12000.su csq_4layer hdr12000.su',
     '''
     /bin/cp ${SOURCES[2]} $TARGET &&
     acd deriv=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} source=${SOURCES[0]} sampord=1 data=$TARGET
     ''',stdin=0,stdout=-1,workdir='shot1200.work')

# linearized modeling
Flow('born12000.su','wavelet12000.su csq_4layersm dcsq_4layer hdr12000.su',
     '''
     /bin/cp ${SOURCES[3]} $TARGET &&
     acd deriv=1 adjoint=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} csq_d1=${SOURCES[2]} source=${SOURCES[0]} sampord=1 data=$TARGET
     ''',stdin=0,stdout=-1,workdir='born12000.work')

# rtm
Flow('migr12000','born12000.su csq_4layersm dcsq_4layer wavelet12000.su',
     '''
     scale < ${SOURCES[2]} > $TARGET dscale=0.0 &&
     acd deriv=1 adjoint=1 nsnaps=10 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} csq_b1=$TARGET source=${SOURCES[3]} sampord=1 data=${SOURCES[0]}
     ''',stdin=0,stdout=-1,workdir='migr12000.work')

# 12 shot "line"
Flow('line8-12km.su','wavelet8-12km.su csq_4layer hdr8-12km.su',
     '''
     /bin/cp ${SOURCES[2]} $TARGET &&
     acd deriv=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} source=${SOURCES[0]} sampord=1 data=$TARGET
     ''',stdin=0,stdout=-1,workdir='line8-12km.work')

# 12 shot linearized simulation
Flow('born8-12km.su','wavelet8-12km.su csq_4layersm dcsq_4layer hdr8-12km.su',
     '''
     /bin/cp ${SOURCES[3]} $TARGET &&
     acd deriv=1 adjoint=0 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} csq_d1=${SOURCES[2]} source=${SOURCES[0]} sampord=1 data=$TARGET
     ''',stdin=0,stdout=-1,workdir='born8-12km.work')

# rtm
Flow('migr8-12km','born8-12km.su csq_4layersm dcsq_4layer wavelet8-12km.su',
     '''
     scale < ${SOURCES[2]} > $TARGET dscale=0.0 &&
     acd deriv=1 adjoint=1 nsnaps=10 order=2 cfl=0.5 cmin=1.0 cmax=5.0 
     csq=${SOURCES[1]} csq_b1=$TARGET source=${SOURCES[3]} sampord=1 data=${SOURCES[0]}
     ''',stdin=0,stdout=-1,workdir='migr8-12km.work')

Result('csq-4layer', 'csq_4layer', 'put label1=Depth unit1=m label2=Distance unit2=m label="V\_p\^\^2\_" unit="m\^2\_/ms\^2\_" | grey color=c mean=y title="Atlantis Cartoon" scalebar=y barreverse=y')

Result('csq-4layersm', 'csq_4layersm', 'put label1=Depth unit1=m label2=Distance unit2=m label="V\_p\^\^2\_" unit="m\^2\_/ms\^2\_" | grey color=c mean=y title="MigVel - Atlantis Cartoon" scalebar=y barreverse=y')

Result('dcsq-4layer', 'dcsq_4layer', 'put label1=Depth unit1=m label2=Distance unit2=m label="V\_p\^\^2\_" unit="m\^2\_/ms\^2\_" | grey color=c mean=y title="Delta V\_p\^\^2\_ - Atlantis Cartoon" scalebar=y barreverse=y')

Result('gauss', 'init', 'put label1=Depth unit1=m label2=Distance unit2=m label="Pressure" unit="GPa" | grey color=c mean=y clip=100 title="Gaussian Data t=0.0 s" scalebar=y barreverse=y')

Result('frameinit', 'movieinit', 'window f3=32 n3=1 |put label1=Depth unit1=m label2=Distance unit2=m label="Pressure" unit="GPa" | grey color=c mean=y clip=300 title="Gaussian data t=5.12 s" scalebar=y barreverse=y')

Result('wavelet','wavelet_base.su', 'suread endian=0 read=data | put label1=Time label2=Pressure unit1=s unit2=GPa title="Gaussian deriv fmax=5 Hz" unit="GPa" |sfgraph')

Result('framesrc', 'moviesrc', 'window f3=32 n3=1 |put label1=Depth unit1=m label2=Distance unit2=m label="Pressure" unit="GPa" | grey color=c mean=y clip=1.e+6 title="Point source 5 Hz t=5.12 s" scalebar=y barreverse=y')

Result('shot12000','shot12000.su', 'suread endian=0 read=data | put labeldcsq_4layer.1=Time label2=Distance d2=0.025 o2=5 unit1=s unit2=km title="Node at 12 km" label="Pressure" unit="GPa" | grey scalebar=y barreverse=y')

Result('born12000','born12000.su', 'suread endian=0 read=data | put label1=Time label2=Trace unit1=s title="Born x_s=12000 m" unit="GPa" | grey')

Result('migr12000', 'migr12000', 'put label1=Depth unit1=m label2 = Distance unit2=m label="V\_p\^\^2\_" unit="m\^2\_/ms\^2\_" |window max1=8000 min2=5000 max2=15000 | grey color=c mean=y title="RTM of Born Shot x_s = 12000" clip=1.e+8')

Result('line8-12km','line8-12km.su', 'suread endian=0 read=data | put label1=Time label2=Node d2=0.4km o2=8 unit1=s unit2=km title="Nodes @0.4km" label="Pressure" unit="GPa" | grey scalebar=y barreverse=y')

Result('born8-12km','born8-12km.su', 'suread endian=0  read=data | put label1=Time label2=Trace unit1=s title="Born x_s=8-12km" unit="GPa" | grey')

Result('migr8-12km', 'migr8-12km', 'put label1=Depth unit1=m label2 = Distance unit2=m label="V\_p\^\^2\_" unit="m\^2\_/ms\^2\_" |window max1=8000 min2=5000 max2=15000 | grey color=c mean=y title="RTM of Born Shot x_s = 8-12km" clip=1e+8')

End()