from rsf.proj import *

nt=351
dt=0.004

Flow('source',None,
     '''
     spike n1=%d d1=%g k1=101 | 
     ricker1 frequency=15
     '''%(nt,dt))

#Parameter specification
Flow('velx',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output="1500+30*(x2-1.5)*(x2-1.5)+30*(x3-1)*(x3-1)+40*(x1)*(x1)"
| scale dscale=0.001
     ''')
# Result('velx',
#        '''
#        byte clip=820 bias=860 bar=bar1.rsf barreverse=y |
#        grey3 transp=n pclip=100 screenration=1 color=j
#        label1="X" label2="Y" label3="Z" scalebar=y poly=y bias=1000
#        frame1=200 frame2=550 frame3=600
#        unit1= unit2=  title="V1" 
#        ''')
Flow('vely',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output="1500+40*(x2-1.5)*(x2-1.5)+40*(x3-1)*(x3-1)+60*(x1)*(x1)"
| scale dscale=0.001
     ''')
# Result('vely',
#        '''
#        byte clip=820 bias=860 bar=bar2.rsf barreverse=y |
#        grey3 transp=n pclip=100 screenration=1 color=j
#        frame1=200 frame2=550 frame3=600
#        unit1= unit2=   title="V2" 
#        ''')
Flow('velz',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output="1500+35*(x2-1.5)*(x2-1.5)+40*(x3-1)*(x3-1)+50*(x1)*(x1)"
| scale dscale=0.001
     ''')
#Result('velz',
#        '''
#        byte clip=820 bias=860 bar=bar.rsf barreverse=y |
#        grey3 transp=n pclip=100 screenration=1 color=j
#        frame1=200 frame2=550 frame3=600
#       label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
#       unit1=km unit2=km unit3=km flat=y screenratio=0.8 
#       scalebar=y barlabel="V\_p" barunit="km/s"
#       wanttitle=n
#        ''')

Flow('eta1',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output=0.3
     ''')
Flow('eta2',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output=0.1
     ''')
Flow('gamma',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=3 k2=251 k3=126 |
     math output=1.03
     ''')

# Conversion to cij
Flow('c33','velz',
     '''
     math output='input*input'
     ''')
Flow('c44','vely',
     '''
     math output='0.5*input*input'
     ''')
Flow('c55','velx',
     '''
     math output='0.5*input*input'
     ''')
Flow('c66','velx',
     '''
     math output='0.25*input*input'
     ''')

Flow('c13','velx c33 c55',
     '''
     math c33=${SOURCES[1]} c55=${SOURCES[2]} output='-c55+sqrt(c55*c55+(input*input*(c33-c55)-c33*c55))'
     ''')
Flow('c23','vely c33 c44',
     '''
     math c33=${SOURCES[1]} c44=${SOURCES[2]} output='-c44+sqrt(c44*c44+(input*input*(c33-c44)-c33*c44))'
     ''')


Flow('c11','eta1 c33 c13 c55',
     '''
     math c33=${SOURCES[1]} c13=${SOURCES[2]} c55=${SOURCES[3]} output='(input+0.5)*(2*c13*(c13+2*c55)+2*c33*c55)/(c33-c55)'
     ''')
Flow('c22','eta2 c33 c23 c44',
     '''
     math c33=${SOURCES[1]} c23=${SOURCES[2]} c44=${SOURCES[3]} output='(input+0.5)*(2*c23*(c23+2*c44)+2*c33*c44)/(c33-c44)'
     ''')
Flow('c12','gamma c11 c66',
     '''
     math c11=${SOURCES[1]} c66=${SOURCES[2]} output='sqrt((input*input-1)*(c11*(c11-c66))+(c11-c66)*(c11-c66))-c66'
     ''')

# Calculation of q
Flow('q12','c11 c13 c33 c55',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))'
    ''')
Flow('q32','c11 c13 c33 c55',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))'
    ''')
Flow('q21','c22 c23 c33 c44',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))'
    ''')
Flow('q31','c22 c23 c33 c44',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))'
    ''')
Flow('q23','c22 c12 c11 c66',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(input-c55)+(c55+c13)*(c55+c13))/(c33*(input-c55))'
    ''')
Flow('q13','c22 c12 c11 c66',''' 
    math c13=${SOURCES[1]} c33=${SOURCES[2]} c55=${SOURCES[3]}
    output='(c55*(c33-c55)+(c55+c13)*(c55+c13))/(input*(c33-c55))'
    ''')

Flow('refl',None,
     '''
     spike n1=256 n2=256 n3=256 d1=0.025 d2=0.025 d3=0.025 o1=-3.175 o2=-3.175 o3=-3.175 unit1=km unit2=km unit3=km k1=128 k2=128 k3=128 |
     smooth rect1=2 rect2=2 rect3=3 repeat=3
     ''')
Flow('seta1','c11','math output=45')
Flow('seta2','c11','math output=45')
#Flow('seta','seta1 seta2','cat ${SOURCES[1]}  axis=4') 

Flow('vp','c33','''math output='sqrt(input)' ''')

# Lowrank decomposition #####################################################
Flow('fftc','vp','fft1 pad=1 | fft3 axis=2 pad=1 | fft3 axis=3 pad=1')

# Exact ######################################################################
Flow('right0 left0','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2',
     '''
     ortholr3 seed=2015 npk=10 eps=0.00001 dt=%g approx=0 tilt=y mode=0
     fft=${SOURCES[1]} c12=${SOURCES[2]} c13=${SOURCES[3]} 
     c22=${SOURCES[4]} c23=${SOURCES[5]} c33=${SOURCES[6]} 
     c44=${SOURCES[7]} c55=${SOURCES[8]} c66=${SOURCES[9]}
     seta1=${SOURCES[10]} seta2=${SOURCES[11]}
     left=${TARGETS[1]} right=${TARGETS[0]}
     ''' % dt)

Flow('wave0 snaps0','source refl left0 right0',
     '''
     fftwave3p ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y snap=50 snaps=${TARGETS[1]}
     ''')

Result('wave0',
       '''
       window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all |
       grey3 frame1=92 frame2=92 frame3=92 screenratio=1
       labelfat=6 titlefat=6 labelsz=9 titlesz=9
       wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
       ''')

Plot('snaps0',
     '''
     byte gainpanel=all |
     grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7
     wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n
     ''')

# Zone approximation ######################################################################
Flow('right1 left1','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2',
     '''
     ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=1 relation=3,3,3 tilt=y mode=0 
     fft=${SOURCES[1]} c12=${SOURCES[2]} c13=${SOURCES[3]} 
     c22=${SOURCES[4]} c23=${SOURCES[5]} c33=${SOURCES[6]} 
     c44=${SOURCES[7]} c55=${SOURCES[8]} c66=${SOURCES[9]}
     seta1=${SOURCES[10]} seta2=${SOURCES[11]}
     left=${TARGETS[1]}
     right=${TARGETS[0]}
     ''' % dt)

Flow('wave1 snaps1','source refl left1 right1',
     '''
     fftwave3p ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y snap=50 snaps=${TARGETS[1]}
     ''')

Result('wave1',
       '''
       window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all |
       grey3 frame1=92 frame2=92 frame3=92 screenratio=1
       labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
       ''')

Plot('snaps1',
     '''
     byte gainpanel=all |
     grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7
     wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n
     ''')


# Acoustic approximation ######################################################################
Flow('right2 left2','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2',
     '''
     ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=2 tilt=y mode=0 
     fft=${SOURCES[1]} c12=${SOURCES[2]} c13=${SOURCES[3]} 
     c22=${SOURCES[4]} c23=${SOURCES[5]} c33=${SOURCES[6]} 
     c44=${SOURCES[7]} c55=${SOURCES[8]} c66=${SOURCES[9]}
     seta1=${SOURCES[10]} seta2=${SOURCES[11]}
     left=${TARGETS[1]}
     right=${TARGETS[0]}
     ''' % dt)

Flow('wave2 snaps2','source refl left2 right2',
     '''
     fftwave3p ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y snap=50 snaps=${TARGETS[1]}
     ''')

Result('wave2',
       '''
       window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all |
       grey3 frame1=92 frame2=92 frame3=92 screenratio=1
       labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
       ''')

Plot('snaps2',
     '''
     byte gainpanel=all |
     grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7
     wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n
     ''')


# Weak Anisotropy approximation ######################################################################

Flow('right3 left3','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2',
     '''
     ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=3 tilt=y mode=0 
     fft=${SOURCES[1]} c12=${SOURCES[2]} c13=${SOURCES[3]} 
     c22=${SOURCES[4]} c23=${SOURCES[5]} c33=${SOURCES[6]} 
     c44=${SOURCES[7]} c55=${SOURCES[8]} c66=${SOURCES[9]}
     seta1=${SOURCES[10]} seta2=${SOURCES[11]}
     left=${TARGETS[1]}
     right=${TARGETS[0]}
     ''' % dt)

Flow('wave3 snaps3','source refl left3 right3',
     '''
     fftwave3p ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y snap=50 snaps=${TARGETS[1]}
     ''')

Result('wave3',
       '''
       window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 | byte gainpanel=all |
       grey3 frame1=92 frame2=92 frame3=92 screenratio=1
       labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
       ''')

Plot('snaps3',
     '''
     byte gainpanel=all |
     grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7
     wanttitle=n  flat=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
     ''')

# Zone approximation (nine) ######################################################################
Flow('right4 left4','c11 fftc c12 c13 c22 c23 c33 c44 c55 c66 seta1 seta2',
     '''
     ortholr3 seed=2010 npk=10 eps=0.00001 dt=%g approx=1 nine=y tilt=y mode=0 
     fft=${SOURCES[1]} c12=${SOURCES[2]} c13=${SOURCES[3]} 
     c22=${SOURCES[4]} c23=${SOURCES[5]} c33=${SOURCES[6]} 
     c44=${SOURCES[7]} c55=${SOURCES[8]} c66=${SOURCES[9]}
     seta1=${SOURCES[10]} seta2=${SOURCES[11]}
     left=${TARGETS[1]}
     right=${TARGETS[0]}
     ''' % dt)

Flow('wave4 snaps4','source refl left4 right4',
     '''
     fftwave3p ref=${SOURCES[1]} left=${SOURCES[2]} right=${SOURCES[3]} verb=y snap=50 snaps=${TARGETS[1]}
     ''')

Result('wave4',
       '''
       window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3 |byte gainpanel=all |
       grey3 frame1=92 frame2=92 frame3=92 screenratio=1
       labelfat=6 titlefat=6 labelsz=9 titlesz=9 wanttitle=n flat=y label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
       ''')

Plot('snaps4',
     '''
     byte gainpanel=all |
     grey4 frame1=128 frame2=128 frame3=128 point1=0.4 point2=0.7
     wanttitle=n label1="n\_3\^" label2="n\_1\^" label3="n\_2\^" flat=n
     ''')

clip2=0.0003
clip1=0.00015
for case in range(4):
        clip=clip1
        if case+1 == 1:
                method = 'proposed approximation'
        elif case+1 == 2:
                method = 'acoustic approximation'
        elif case+1 == 3:
                method = 'Weak-anisotropy approximation'
                clip=clip2
        else:
                method = 'proposed approximation (nine)'
        Flow('error%d' % (case+1),'wave0 wave%d' % (case+1),'add scale=-1,1 ${SOURCES[1]}| window max1=2.3 min1=-2.3 max2=2.3 min2=-2.3 max3=2.3 min3=-2.3')
        Result('error%d' % (case+1),
        ''' 
        byte gainpanel=all bar=bar.rsf clip=%g allpos=y maxval=%g minval=0.0|
        grey3 label1="n\_3\^" label2="n\_1\^" label3="n\_2\^"
        unit1=km unit2=km unit3=km screenratio=0.8
        scalebar=y barlabel="Absolute Error" title="Amplitude error of the %s"
        frame1=92 frame2=92 frame3=92
        flat=y wanttitle=n labelfat=6 titlefat=6 labelsz=9 titlesz=9
        ''' % (clip,clip,method) )


End()