from rsf.proj import *
import fdmodds

yshot=4
yshot2=4.025
yshot3=4.05
zshot=0.0
dl=0.025
dll=-0.025

# Fetch Files from repository
raw=['marmvel.hh','marmsmooth.HH']
for file in raw:
    Fetch(file,"marm")
    if file is 'marmvel.hh':
        d=.004
        fileOut='marmvel'
        t='Velocity\ Model'
    if file is 'marmsmooth.HH':
        d=.024
        fileOut='marmsmooth'
        t='Smoothed\ Velocity\ Model'
# Convert Files to RSF and update headers
    Flow(fileOut,file,'''dd form=native |
        scale rscale=.001 | put
        label1=Depth label2=Position unit1=km unit2=km
        d1=%f d2=%f''' % (d,d))
# Plotting Section
    Result(fileOut,'''window $SOURCE  |
        grey color=I gainpanel=a allpos=y scalebar=y
        title=%s barlabel=\(km\/s\) screenratio=.326
        screenht=3 wherettle=t labelsz=4 ttlesz=6 ''' % t)

Flow('marm1','marmsmooth','window min2=1 max2=7 max1=2.9 j1=1 j2=1')
Flow('marm2','marmsmooth','window min2=1.025 max2=7.025 max1=2.9 j1=1 j2=1')
Flow('marm3','marmsmooth','window min2=1.05 max2=7.05 max1=2.9 j1=1 j2=1')
Plot('marm1','grey allpos=y screenratio=.327 screenht=4.5 labelsz=12 wanttitle=n')
Plot('marm2','grey allpos=y screenratio=.327 screenht=4.5 wanttitle=n')
Plot('marm3','grey allpos=y screenratio=.327 screenht=4.5 wanttitle=n')
Flow('marmO','marmvel','window min2=1 max2=7')
Plot('marmO','grey allpos=y screenratio=.327 screenht=4.5 wanttitle=n scalebar=y labelsz=12 barlabel=Velocity\(km\/s\)')
Result('medium','marmO marm1','SideBySideIso')

# ------------------------------------------------------------
par = {
    'nt':10000, 'dt':0.00025,'ot':0,'lt':'t','ut':'s',
    'nx':251, 'ox':1, 'dx':.024, 'lx':'x','ux':'km',
    'nz':122,  'oz':0, 'dz':.024, 'lz':'z','uz':'km',
    'kt':400    # wavelet delay
    }
# add F-D modeling parameters
fdmodds.param(par)
# ------------------------------------------------------------
# wavelet
Flow('wav',None,
         '''spike nsp=1 mag=1 n1=%(nt)d d1=%(dt)g o1=%(ot)g k1=%(kt)d |
         ricker1 frequency=15 | scale axis=123 | 
         put label1=t label2=x label3=y | transp''' % par)
Result('wav','transp | window n1=1000 | graph title="" label1="t" label2=')
# ------------------------------------------------------------
# experiment setup
Flow('r_',None,'math n1=%(nx)d d1=%(dx)g o1=%(ox)g output=0' % par)
Flow('s_',None,'math n1=1      d1=0      o1=0      output=0' % par)
# receiver positions
Flow('zr','r_','math output=.025')
Flow('xr','r_','math output="x1"')
Flow('rr',['xr','zr'],'''cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} | transp
     ''', stdin=0)
Plot('rr',fdmodds.rrplot('',par))
# source positions
Flow('zs','s_','math output=.01')
Flow('xs','s_','math output=5.0')
Flow('rs','s_','math output=1')
Flow('ss',['xs','zs','rs'],'''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} ${SOURCES[2]} | transp
     ''', stdin=0)
Plot('ss',fdmodds.ssplot('',par))
# source positions 2
Flow('zss','s_','math output=.01')
Flow('xss','s_','math output=5.05')
Flow('rss','s_','math output=1')
Flow('sss',['xss','zss','rss'],'''
     cat axis=2 space=n
     ${SOURCES[0]} ${SOURCES[1]} ${SOURCES[2]} | transp
     ''', stdin=0)
Plot('sss',fdmodds.ssplot('',par))
# ------------------------------------------------------------
# density
Flow('vel','marm1',
      '''
      put o1=%(oz)g d1=%(dz)g  o2=%(ox)g d2=%(dx)g
      ''' % par)
Plot('vel',fdmodds.cgrey('''allpos=y bias=1.5 pclip=97 title=Survey\ Design 
                  color=F ttlesz=6 labelsz=4 wherettle=t barrevers=y''',par))
Result('vel',['vel','rr','ss'],'Overlay')
# ------------------------------------------------------------
# density
Flow('den','vel','math output=1')
# ------------------------------------------------------------
# density 2
Flow('vels','marm2',
      '''
      put o1=%(oz)g d1=%(dz)g  o2=%(ox)g d2=%(dx)g
      ''' % par)
Plot('vels',fdmodds.cgrey('''allpos=y bias=1.5 pclip=97 title=Survey\ Design 
                  color=F ttlesz=6 labelsz=4 wherettle=t barrevers=y''',par))
Result('vels',['vels','rr','ss'],'Overlay')
# ------------------------------------------------------------
# density
Flow('dens','vels','math output=1')
# ------------------------------------------------------------
# finite-differences modeling
fdmodds.awefd('dat','wfl','wav','vel','den','ss','rr','free=y dens=y',par)

Plot('wfl',fdmodds.wgrey('pclip=99',par),view=1)
Plot('dat','transp |' + fdmodds.dgrey('''pclip=99 title=Data\ Record label2=Offset clip=.000002
            wherettle=t ttlesz=6 labelsz=12''',par))

times=['.4','0.8','1.2','1.6']
cntr=0
for item in ['16','32','48','64']:
    Plot('time'+item,'wfl','''window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a wantunits=n wanttitle=n unit1= unit2=
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    cntr=cntr+1

# ------------------------------------------------------------
# finite-differences modeling new source
fdmodds.awefd('dats','wfls','wav','vels','dens','ss','rr','free=y dens=y',par)

Plot('wfls',fdmodds.wgrey('pclip=99',par),view=1)
Plot('dats','transp |' + fdmodds.dgrey('''pclip=99 title=Data\ Record label2=Offset clip=.000002
            wherettle=t ttlesz=6 labelsz=12''',par))

times=['.4','0.8','1.2','1.6']
cntr=0
for item in ['16','32','48','64']:
    Plot('times'+item,'wfls','''window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a wanttitle=n unit1= unit2=
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Plot('timed'+item,['wfls','wfl'],'''difference subtracter=${SOURCES[1]} |
              window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 wanttitle=n unit1= unit2=
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Result('timeSide2'+item,['time'+item,'times'+item,'timed'+item],'SideBySideIso')
    cntr=cntr+1
# ------------------------------------------------------------
# finite-differences modeling perturbation
fdmodds.awefdds('datds','wflds','wav','vel','den','ss','rr','free=y dens=y dl=0.025',par)

Plot('wflds',fdmodds.wgrey('pclip=99',par),view=1)
Plot('datds','transp |' + fdmodds.dgrey('''pclip=99 title=Data\ Record label2=Offset clip=.000002 
            wherettle=t ttlesz=6 labelsz=12''',par))

times=['.4','0.8','1.2','1.6']
cntr=0
for item in ['16','32','48','64']:
    Plot('timeds'+item,'wflds','''window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a wanttitle=n unit1= unit2=
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Plot('timedds'+item,['wfls','wflds'],'''difference subtracter=${SOURCES[1]} |
              window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a clip=.000005 wanttitle=n unit1= unit2=
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12  label1='Depth (km)' label2='Lateral (km)' 
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Result('timeSide3'+item,['timeds'+item,'times'+item,'timedds'+item],'SideBySideIso')
    cntr=cntr+1
# ------------------------------------------------------------
# finite-differences modeling perturbation 2nd
fdmodds.awefdds2nd('datds2','wflds2','wav','vel','den','ss','rr','free=y dens=y dl=0.025',par)

Plot('wflds2',fdmodds.wgrey('pclip=99',par),view=1)
Plot('datds2','transp |' + fdmodds.dgrey('''pclip=99 title=Data\ Record label2=Offset clip=.000002
            wherettle=t ttlesz=6 labelsz=12''',par))

times=['.4','0.8','1.2','1.6']
cntr=0
for item in ['16','32','48','64']:
    Plot('timeds2'+item,'wflds2','''window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a wanttitle=n unit1= unit2=
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Plot('timedds2'+item,['wfls','wflds2'],'''difference subtracter=${SOURCES[1]} |
              window f3=%s n3=1 min1=0 min2=2.8 | grey gainpanel=a clip=.000005 label1='Depth (km)' label2='Lateral (km)' 
              pclip=99 wantframenum=y ttle=Wavefield\ at\ %s\ \(s\) labelsz=12 wanttitle=n unit1= unit2=
              ttlesz=9 screenratio=.5 screenht=7.5 wherettle=t''' % (item,times[cntr]))
    Result('timeSide4'+item,['timeds2'+item,'times'+item,'timedds2'+item],'SideBySideIso')
    Result('timeSide5'+item,['timed'+item,'timedds'+item,'timedds2'+item],'SideBySideIso')
    cntr=cntr+1

Plot('datdiff',['dats','dat'],'difference subtracter=${SOURCES[1]} | transp |'+ fdmodds.dgrey('''pclip=99 
            title=Data\ Record label2=Offset clip=.000002 o2=-5 wanttitle=n
            wherettle=t ttlesz=6 labelsz=12''',par))
Plot('datdiff2',['dats','datds'],'difference subtracter=${SOURCES[1]} | transp |'+ fdmodds.dgrey('''pclip=99 
            ttle=Data\ Record label2=Offset clip=.000002 o2=-5 wanttitle=n
            wherettle=t ttlesz=6 labelsz=12''',par))
Plot('datdiff3',['dats','datds2'],'difference subtracter=${SOURCES[1]} | transp |'+ fdmodds.dgrey('''pclip=99 
            ttle=Data\ Record label2=Offset clip=.000002 o2=-5 wanttitle=n
            wherettle=t ttlesz=6 labelsz=12''',par))
Result('datdiff',['dats','datdiff','datdiff2'],'SideBySideIso')
Result('datdiffF',['datdiff','datdiff2','datdiff3'],'SideBySideIso')

End()