from rsf.proj import *

nz = 257
dz = 0.004
nx = 256
dx = 0.008


# Dome model
############
z0 = (nz/3.0 + 1.0) * 0.25
z1 = nz/3.0 + 0.5 + 3.*z0
x0 = dx*(nx/2.0 + 1.0) * 0.5

Flow('bot',None,
     'spike n1=115 d1=%g o1=-0.456 mag=%g' % (dx,dz*z1))
Flow('top','bot',
     'math output="%g - %g*sqrt(1-x1*x1*%g)" ' % \
     (dz*z1,dz*2.*z0,4./(3.*x0*x0)))
Flow('two','top bot',
     'cat axis=2 ${SOURCES[1]} | unif2 d1=%g n1=%d v00=1,2,3' % (dz,nz))
Flow('dome','two',
     '''
     igrad |
     bandpass flo=10 fhi=50 |
     pad beg2=71 end2=71 |
     put n3=1 
     ''')

# Lines model
#############
nl = 5
lines = []
for n in range(1,nl+1):
    line = 'line%d' % n
    Flow(line,None,
         'math n1=257 d1=%g output="%g + x1*%g" ' % (dx,0.05*n,0.2*n))
    lines.append(line)
Flow('mod',lines,
     '''
     cat axis=2 ${SOURCES[1:%d]} |
     unif2 d1=0.004 n1=257 v00=%s
     ''' % (nl,string.join(map(str,range(1,nl+2)),',')))
Flow('lines','mod','igrad | bandpass flo=10 fhi=50 | put n3=1')

# General rules
###############

def plots(one,two,name,axis,f3,forw):
    Plot(one,
         'grey label1=Time label2=Space unit1=s unit2=km title=%s' % name)
    Plot(two+'0',two,
         '''
         window n3=1 f3=0 |
         grey label1=Time label2=Midpoint unit1=s unit2=km
         title="Zero %s"
         ''' % axis)
    Plot(two+'1',two,
         '''
         window n3=1 f3=%d |
         grey label1=Time label2=Midpoint unit1=s unit2=km
         title="Far %s"
         ''' % (f3,axis))
    if forw:
        show = [one,two+'0',two+'1']
    else:
        show = [two+'0',two+'1',one]
    Result(two,show,'SideBySideIso')

def agmig(name,modl,vel):
    # Modeling
    data = 'data-'+name
    Flow(data,modl,
         '''
         halfint inv=1 |
         preconstkirch zero=y inv=y h0=0 dh=%g nh=48 vel=%g |
         window
         ''' % (dx,vel))
    plots(modl,data,'Model','Offset',47,1)
    # Common-offset migration
    offset = 'offset-'+name
    Flow(offset,data,
         '''
         transp plane=43 | preconstkirch vel=1.5 |
         halfint inv=y adj=y | put n3=48 n4=1
         ''')
    gather = 'gather-'+name
    Plot(gather,offset,
         '''
         window f2=125 n2=1 n3=25 |
         grey label1=Time unit1=s label2=Half-Offset unit2=km
         title="Offset Gather"
         ''')
    stack = 'stack-'+name
    Flow(stack,offset,'stack axis=3')
    plots(stack,offset,'Stack','Offset',47,0)
    # Angle-gather migration
    angle = 'angle-'+name
    Flow(angle,data,
         '''
         agmig g0=0 ng=48 dg=1 vel=1.5 amax=70 |
         halfint inv=y adj=y
         ''')
    astack = 'astack-'+name
    Flow(astack,angle,'stack axis=3')
    plots(astack,angle,'Stack','Angle',25,0)
    for n3 in (25,48):
        Plot(angle+str(n3),angle,
             '''
             window f2=125 n2=1 n3=%d |
             grey label1=Time unit1=s label2=Angle unit2=degree
             title="Reflection Angle Gather"
             ''' % n3)
    agather = 'agather-'+name
    Result(agather,[gather,angle+'25',angle+'48'],'SideBySideIso')

agmig('dome','dome',1.5)
agmig('lines','lines',1.5)
agmig('dome-fast','dome',3)

End()