from rsf.proj import *

# Download from http://www.freeusp.org/2007_BP_Ani_Vel_Benchmark/
tgz = 'ModelParams.tar.gz'

Fetch(tgz,'BP',top=os.environ.get('DATAPATH'),server='local')

pars = Split('epsilon delta vp theta')

sgy = {}
for par in pars:
    sgy[par] = os.path.join('ModelParams',par.capitalize() + '_Model.sgy')

Flow(sgy.values(),tgz,'zcat $SOURCE | tar -xvf -',stdin=0,stdout=-1)

for par in pars:
    Flow([par,par+'.asc',par+'.bin'],sgy[par],
         '''
         segyread hfile=${TARGETS[1]} bfile=${TARGETS[2]} read=d |
         put
         o2=0 d2=0.00625 label2=Distance unit2=km
         o1=0 d1=0.00625 label1=Depth unit1=km %s |
          transp plane=12 
         ''' % ('','| scale dscale=0.001')[par=='vp'])
#    Flow([par,par+'.asc',par+'.bin'],sgy[par],
#         '''
#         segyread hfile=${TARGETS[1]} bfile=${TARGETS[2]} |
#         put
#         o2=0 d2=0.00625 label2=Distance unit2=km
#         o1=0 d1=0.00625 label1=Depth unit1=km %s
#         ''' % ('','| scale dscale=0.001')[par=='vp'])
    Result(par,
           '''
           window j1=8 j2=2 |
           grey color=j barlabel=%s scalebar=y
           screenwd=12.595 screenht=1.8
           labelsz=4 titlesz=5 barreverse=y
           wanttitle=n allpos=%d bias=%g barunit=%s
           parallel2=n transp=n
           ''' % (par.capitalize(),
                  par != 'theta',
                  (0,1.5)[par=='vp'],
                  ('','km/s')[par=='vp']))
Flow('vx','vp epsilon',
     '''
     math e=${SOURCES[1]} output="input*sqrt(2.0*e+1.0)"
     ''')
Flow('yita','epsilon delta',
     '''
     math e=${SOURCES[0]} d=${SOURCES[1]} output="((1.0+2.0*e)/(1.0+2.0*d)-1.0)/2.0"
     ''')
for par in ('vx','yita'):
    Result(par,
           '''
           window j1=8 j2=2 |
           grey color=j barlabel=%s scalebar=y
           screenwd=12.595 screenht=1.8
           labelsz=4 titlesz=6 barreverse=y
           wanttitle=y allpos=%d bias=%g barunit=%s
           parallel2=n transp=n title=%s
           ''' % (par.capitalize(),
                  par != 'theta',
                  (0,1.5)[par=='vx'],
                  ('','km/s')[par=='vx'],
                  par.capitalize()))

nt=4300
dt=0.001
Flow('source',None,
     '''
     spike n1=%d d1=%g k1=200 | 
     ricker1 frequency=15 
     '''%(nt,dt))
Result('source','graph  title="Source Wavelet" ')
for par in ('vp','vx','yita','theta'):
    Flow(par+'end2',par,'window j1=2 j2=2 | window  n1=2048 f1=2761 n2=900')
    Plot(par+'end2',
           '''
           grey color=j barlabel=%s scalebar=y
           screenwd=12.595 screenht=1.8
           labelsz=4 titlesz=5 barreverse=y
           wanttitle=n allpos=%d bias=%g barunit=%s
           parallel2=n transp=n title=%s
           ''' % (par.capitalize(),
                  par != 'theta',
                  (0,1.5)[par=='vp'],
                  ('','km/s')[par=='vp'],
                  par.capitalize()))
Flow('wavffdend2',['source','vxend2','vpend2','yitaend2','thetaend2'],
     '''
     ffdantti4b_smsr velx=${SOURCES[1]}  velz=${SOURCES[2]} yita=${SOURCES[3]} seta=${SOURCES[4]} dt=%g nt=%d isx=1024 isz=1 nbt=80 nbb=44 nbl=0 nbr=0 
     ct=0.002 cb=0.002 cl=0.0 cr=0.0 err=0.0001
     '''%(dt,nt) )
Plot('wavffdend2',
     ''' 
     window j3=300 |
     grey poly=y label2="Z" label1="X" title="TTI" transp=n min1=38 max1=56
     yreverse=y transp=n gainpanel=each pclip=100
     ''',view=1)
Result('ttisnapshot','wavffdend2',
       ''' 
       window n3=1 f3=4200 |
       grey poly=y label2="Depth" label1="Distance" title="Snapshot in the BP 2D TTI model" transp=n min1=38 max1=56
       yreverse=y transp=n gainpanel=each pclip=100 unit1=km unit2=km screenwd=14 screenht=10 
       ''' )

Plot('vp2','vp',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel=%s scalebar=y
     labelsz=6 titlesz=8 barreverse=y
     wanttitle=y allpos=%d bias=%g barunit=%s
     parallel2=n transp=n title=a
     screenwd=15 screenht=11
     '''
     % ('Vp',
     1,
     1.5,
     'km/s'
     ))
Plot('epsilon2','epsilon',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel=%s scalebar=y
     labelsz=6 titlesz=8 barreverse=y
     wanttitle=y allpos=%d bias=%g barunit=%s
     parallel2=n transp=n title=b
     screenwd=15 screenht=11
     '''
     % ('Epsilon',
     1,
     0,
     ''
     ))
Plot('delta2','delta',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel=%s scalebar=y
     labelsz=6 titlesz=8 barreverse=y
     wanttitle=y allpos=%d bias=%g barunit=%s
     parallel2=n transp=n title=c
     screenwd=15 screenht=11
     '''
     % ('Delta',
     1,
     0,
     ''
     ))
Plot('theta2','theta',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel=%s scalebar=y
     labelsz=6 titlesz=8 barreverse=y
     wanttitle=y allpos=%d bias=%g barunit=%s
     parallel2=n transp=n title=d
     screenwd=15 screenht=11
     '''
     % ('Theta',
     0,
     0,
     ''
     ))
Result('ttimodel','vp2 epsilon2 delta2 theta2','TwoRows')
Result('vp2','vp',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0
     '''
     % (1,
     1.5,
     ))
Result('epsilon2','epsilon',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0
     '''
     % (1,
     0,
     ))
Result('delta2','delta',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0
     '''
     % (1,
     0,
     ))
Result('theta2','theta',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0

     '''
     % (0,
     0,
     ))
Result('vx2','vx',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0
     '''
     % (1,
     1.5,
     ))
Result('yita2','yita',
     '''
     window j1=4 j2=2 | window  n1=1024 f1=1381 n2=900 |
     grey color=j barlabel= scalebar=y
     labelsz=9 titlesz=10 barreverse=y
     wanttitle=n allpos=%d bias=%g 
     parallel2=y transp=n title= screenratio=1.0
     '''
     % (0,
     0,
     ))
Plot('wavffdend2',
     ''' 
     window j3=300 |
     grey poly=y label2="Z" label1="X" title="TTI" transp=n min1=38 max1=56
     yreverse=y transp=n gainpanel=each pclip=100 unit1="km" unit2="km"
     ''',view=1)
End()