from rsf.proj import *

h1=1.
h2=0.75
t1=1.

Flow('x0',None,
     '''
     math
     n1=401 n2=401 n3=2
     o1=-1.8 o2=-0.45 o3=0.05235987756
     d1=0.009 d2=0.00225 d3=0.47123889804
     output="x1-x2*cos(x3)/sin(x3)"
     ''')
Flow('vk','x0','math output="(%g*sin(x3)^2-x2*x2)/(%g*sin(x3)^2-(x1*sin(x3)-x2*cos(x3))^2)" ' % (h2*h2,h1*h1))
Flow('vmask','vk','mask min=0 | dd type=float')
Flow('time','vk','math output="%g*sqrt(abs(input))" ' % (t1*h1/h2))

for case in range(2):
    v1 = (2.0,0.01)[case]
    tmean = (0.8,0.4)[case]

    r=v1*t1/2
    dt=2*h1/v1
    b=t1*t1/(t1*t1+dt*dt)
    r=r*r

    rx = 'rx%d' % case
    ry = 'ry%d' % case

    Flow(rx,'x0','scale dscale=%g' % (1/(1-b)))
    Flow(ry,['x0',rx],
         '''
         math x0=${SOURCES[0]} rx=${SOURCES[1]}
         output="(x0-rx)*cos(x3)/sin(x3)-x2*((x0-rx)^2-%g*rx*rx+%g)/(%g*sin(x3)^2-x2*x2)"
         ''' % (b,r,h2*h2))

    time = 'time-%d' % case

    Flow(time,[rx,ry,'time','vmask'],
         '''
         math rx=${SOURCES[0]} ry=${SOURCES[1]} output="%g-%g*rx*rx-ry*ry" |
         mask min=0 | dd type=float | 
         math t=${SOURCES[2]} m=${SOURCES[3]} output="input*m*abs(t-%g)"
         ''' % (r,b,tmean))

Flow('time4','time-0','window n3=1')
Flow('time2','time-0','window f3=1')
Flow('time3','time-1','window n3=1')
Flow('time1','time-1','window f3=1')

title = (
    "phi=30, v=10 m/s",
    "phi=30, v=2000 m/s",
    "phi=3, v=10 m/s",
    "phi=3, v=2000 m/s")

for n in range(4):
    time = 'time%d' % (n+1)
    Plot(time,
         '''
         contour title="AMO impulse response: %s"
         nc=%d c0=0.01 dc=0.02 transp=n yreverse=n 
         ''' % (title[n],(200,50,200,50)[n]))

Result('amoapp','time1 time3 time2 time4','TwoRows')

End()