from rsf.proj import *

xmax = 6.0
zmax = 2.0

layers = ((0.30,0.50,0.20,0.30),
          (0.55,0.75,0.45,0.55),
          (0.65,0.85,0.55,0.65),
          (1.30,1.30,1.60,1.20))

velocities = (1.508,
              1.581,
              1.690,
              1.826,
              2.000)

def arr2str(array,sep=' '):
    return string.join(map(str,array),sep)

vstr = arr2str(velocities,',')

n1 = len(layers[0])
n2 = len(layers)

Flow('layers.asc',None,
     '''
     echo %s
     n1=%d n2=%d o1=0 d1=%g
     data_format=ascii_float in=$TARGET     
     ''' % (string.join(map(arr2str,layers),' '),
            n1,n2,xmax/(n1-1)))
Flow('layers','layers.asc','dd form=native')

d = 0.0101 # non-round for reproducibility

Flow('refs','layers',
     'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d)))
Flow('dips','refs','deriv scale=y')

Flow('mod1','refs',
     '''
     unif2 d1=%g n1=%d v00=%s 
     ''' % (d,int(1.5+zmax/d),vstr))

Result('mod1',
       '''
       grey color=j title="Model 1" 
       screenratio=%g screenht=4
       mean=y titlesz=8 labelsz=6
       label1="Depth (km)"
       label2="Distance (km)"
       ''' % (zmax/xmax))

shots = (2.0,2.6,3.3,3.95)
plots = []
for s in range(4):
    shot = 'shot%d' % s
    Flow(shot,'refs dips',
         '''
         kirmod_newton nt=751 dt=0.004 freq=15
         ns=1 s0=%g ds=0.01 nh=60 dh=0.04 h0=-1.2 verb=y 
         vstatus=0 velocity=%s debug=n fwdxini=y
         xref=0 zref=0 dip=${SOURCES[1]}
         ''' % (shots[s],vstr))
    Plot(shot,
         '''
         wiggle transp=y yreverse=y poly=y
         title="Shot at %g km"
         ''' % shots[s])
    plots.append(shot)

Result('shots',plots,'SideBySideAniso')

End()