from rsf.proj import *

# Reflector model

Flow('dome',None,
     '''
     math d1=0.01 n1=1201 o1=-1 unit1=km label1=Distance
     output="5-3*exp(-(x1-5)^2/9)"
     ''')
Result('dome','graph min2=0 max2=5.5 yreverse=y title=Dome')

Flow('ref','dome',
     '''
     spray axis=2 n=5 o=0 d=0.25 |
     math output="(input-1)*x2*x2+1"
     ''')
Result('ref','graph min2=0 max2=5.5 yreverse=y label2=Depth unit2=km title=')

Flow('dip','ref',
     'math output="2*(x1-5)/3*exp(-(x1-5)^2/9)*x2*x2" ')

# Isotropic Data

Flow('data','ref dip',
     '''
     kirmod cmp=y dip=${SOURCES[1]} 
     nh=129  dh=0.05  h0=0
     ns=513 ds=0.025 s0=0
     freq=20 dt=0.004 nt=1501
     vel=2 
     ''',split=[1,1201],reduce='add')

Result('data',
       '''
       window min3=0 max3=10 |
       byte |
       transp plane=23 |
       grey3 flat=n frame1=500 frame3=20 frame2=320
       label1=Time unit1=s 
       label3=Offset unit3=km 
       label2=Midpoint unit2=km
       title=Data point1=0.8 point2=0.8
       ''')

# Anisotropic Data

Flow('adata','ref dip',
     '''
     kirmod cmp=y dip=${SOURCES[1]} 
     nh=129 dh=0.05  h0=0
     ns=513 ds=0.025 s0=0
     freq=20 dt=0.004 nt=1501 type=aniso
     vel=2 velz=2 eta=0.2 
     ''',split=[1,1201],reduce='add')

Result('adata',
       '''
       window min3=0 max3=10 |
       byte |
       transp plane=23 |
       grey3 flat=n frame1=500 frame3=20 frame2=320
       label1=Time unit1=s 
       label3=Offset unit3=km 
       label2=Midpoint unit2=km
       title=Data point1=0.8 point2=0.8
       ''')

# Make impulse

Flow('impuls',None,
     '''
     spike
     n1=1501 d1=0.004 k1=500
     n2=513  d2=0.025 o2=0 k2=250
     n3=129  d3=0.05  k3=60 label3=Offset |
     ricker1 frequency=5 |
     fft1 | window max1=15 |
     transp plane=12 memsize=1000 |
     transp plane=23 memsize=1000 |
     spray axis=2 n=1
     ''')

# Make slowness model

Flow('slow','dome','math output=0.5 | spray axis=1 d=0.01 n=501 | transp | spray axis=2 n=1')

# DSR migration

Flow('image','impuls slow',
     'camig3 --readwrite=y inv=n eps=0.1 tmx=32 nr=5 pmx=32 pmy=0 phx=0 slo=${SOURCES[1]} ompnth=1',split=[4,93,[0]],reduce='add')

# Make impulse for sfdsr

Flow('impuls2',None,
     '''
     spike
     n1=1501 d1=0.004 k1=500
     n3=513  d3=0.025 o3=0 k3=250
     n2=129  d2=0.025 k2=60 label2=Offset |
     ricker1 frequency=5 |
     cosft sign2=1 sign3=1
     ''')

# Make velocity for sfdsr

Flow('velo','slow','window n1=1 | math output=1/input')
Flow('velz','slow','window n1=1 | math output=1/input')
Flow('eta','slow','window n1=1 | math output=0.2')

m = -1+249*0.025
h = 59*0.05*0.5
s = m-h
r = m+h
v = 2
t = 499*0.004

nw = 400 # maximum number of frequencies

def dsr(data,image):
    Flow(image,[data,'velo'],
         'dsr inv=0 velocity=${SOURCES[1]} eps=1 depth=y nw=%d' % nw,
         split=[3,513,[0]])

def adsr(data,image):
    Flow(image,[data,'velo','velz','eta'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} 
         velz=${SOURCES[2]} eta=${SOURCES[3]} eps=1 depth=y 
         rule=a arule=a nw=%d
         ''' % nw, split=[3,513,[0]])

# DSR migration - impulse response

Flow('image','impuls slow',
     'camig3 --readwrite=y inv=n eps=0.1 tmx=32 nr=5 pmx=32 pmy=0 phx=0 slo=${SOURCES[1]} ompnth=1',split=[4,93,[0]],reduce='add')

dsr('impuls2','dsr')
Flow('image2','dsr','window | cosft sign2=-1')

Plot('image2','grey wanttitle=n')

Plot('time2','image2',
     '''
     math output="(sqrt((%g-x2)^2+x1^2)+sqrt((%g-x2)^2+x1^2))/%g" |
     contour nc=1 c0=%g wanttitle=n wantaxis=n
     ''' % (s,r,v,t))

Result('image','image2 time2','Overlay')

# DSR migration - data

Flow('cdata','data',
     '''
     tpow tpow=1 |
     put d2=0.025 |
     cosft sign2=1 sign3=1
     ''')
dsr('cdata','ddata')
Flow('mig','ddata','window | cosft sign2=-1')
Result('mig','window max2=10 | grey title=Migration')

# DSR migration - anisotropic data

Flow('cadata','adata',
     '''
     tpow tpow=1 |
     put d2=0.025 |
     cosft sign2=1 sign3=1
     ''')
adsr('cadata','dadata')
Flow('amig','dadata','window | cosft sign2=-1')
Result('amig','window max2=10 | grey title=Migration')

def gath(data,image):
    Flow(image,[data,'velo'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} eps=0.1 depth=y na=60 da=1 nw=%d
         ''' % 400,split=[3,513,[0]])
def agath(data,image):
    Flow(image,[data,'velo','velz','eta'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} 
         velz=${SOURCES[2]} eta=${SOURCES[3]} 
         rule=a eps=0.1 depth=y na=60 da=1 nw=%d
         ''' % 400, split=[3,513,[0]])

def aagath(data,image):
    Flow(image,[data,'velo','velz','eta'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} 
         velz=${SOURCES[2]} eta=${SOURCES[3]} 
         rule=a arule=a eps=0.1 depth=y na=30 da=2 nw=%d
         ''' % 400, split=[3,513,[0]])

def aagathtest(data,image):
    Flow(image,[data,'velo','velz','eta'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} 
         velz=${SOURCES[2]} eta=${SOURCES[3]} 
         rule=a arule=i eps=0.1 depth=y na=30 da=2 nw=%d
         ''' % 400, split=[3,513,[0]])

def aigath(data,image):
    Flow(image,[data,'velo','velz','eta'],
         '''
         dsr inv=0 velocity=${SOURCES[1]} 
         velz=${SOURCES[2]} eta=${SOURCES[3]} 
         rule=i arule=a eps=0.1 depth=y na=30 da=2 nw=%d
         ''' % 400, split=[3,513,[0]])

# Angle Gathers - isotropic

gath('cdata','gath')
Flow('agath','gath','cosft sign3=-1 | transp plane=13 memsize=1000 | transp memsize=1000')

Result('agath',
       '''
       window max2=10 max1=4 min2=2 min3=5 max3=45   | 
       byte | 
       grey3 frame1=200 frame2=240 frame3=10 wanttitle=n point1=0.8 point2=0.8
       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_"
       ''')

# Angle Gathers - isotropic construction with anisotropic migration

gath('cadata','gath2')
Flow('agath2','gath2','cosft sign3=-1 | transp plane=13 memsize=1000 | transp memsize=1000')

Result('agath2',
       '''
       window max2=10 max1=4 min2=2 min3=5 max3=45 |
       byte |
       grey3 frame1=200 frame2=240 frame3=10 wanttitle=n point1=0.8 point2=0.8
       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_"
       ''')

#Result('agath2',
#       '''
#       window max2=10 max1=4 min2=2 | 
#       byte |
#       grey3 frame1=200 frame2=160 frame3=10 title="Isotropic Angle Gathers"
#       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_" 
#       ''')
# Angle Gathers - anisotropic migration isotropic angle gathers

agath('cadata','gath3')
Flow('agath3','gath3','cosft sign3=-1 | transp plane=13 memsize=1000 | transp memsize=1000')
.8
Result('agath3',
       '''
       window max2=10 max1=4 min2=2 min3=5 max3=45  |
       byte |
       grey3 frame1=200 frame2=240 frame3=10 wanttitle=n point1=0.8 point2=0.8
       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_"
       ''')


# Angle Gathers - anisotropic migration anisotropic angle gathers

aagath('cadata','gath4')
Flow('agath4','gath4','cosft sign3=-1 | transp plane=13 memsize=1000 | transp memsize=1000')

Result('agath4',
       '''
       window max2=10 max1=4 min2=2 min3=5 max3=45 |
       byte |
       grey3 frame1=200 frame2=240 frame3=10 wanttitle=n point1=0.8 point2=0.8
       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_"
       ''')

# Angle Gathers - anisotropic migration anisotropic angle gathers test

aigath('cadata','gath55')
Flow('agath55','gath55','cosft sign3=-1 | transp plane=13 memsize=1000 | transp memsize=1000')

Result('agath55',
       '''
       window max2=10 max1=4 min2=2 min3=5 max3=45  |
       byte |
       grey3 frame1=200 frame2=240 frame3=10 wanttitle=n point1=0.8 point2=0.8
       label1=Depth unit1=km label2=Distance unit2=km label3=Angle unit3="\^o\_"
       ''')

# Angle Gathers - anisotropic migration anisotropic angle gathers test



End()