from rsf.proj import *

y1=0.18
y2=0.30
x1=25
x2=150

a=(y1-y2)/(x1-x2)
b=(y2*x1-y1*x2)/(x1-x2)

Flow('layer1',None,'spike n1=500 d1=1 mag=0.18 label1=Trace unit1=')

Flow('layer2','layer1','math output="%g*x1+%g" | clip2 lower=0.18 upper=0.3' % (a,b))

Flow('layers','layer1 layer2','cat axis=2 ${SOURCES[1]}')

Result('layers','graph max2=0.6 min2=0 yreverse=y label2=Time unit2=s title=Layers')

Flow('grid','layers','unif2 v00=1,2,3 n1=601 d1=0.001 label1=Time unit1=s')

Result('grid','grey color=J mean=y wanttitle=n')

vp=(3300,3200,3300)
rho=(2600,2550,2650)

ai = [str(vp[k]*rho[k]/1.0e6) for k in range(3)]

Flow('ai','layers','unif2 v00=%s n1=601 d1=0.001 label1=Time unit1=s' % ','.join(ai))

Result('ai','grey color=lb mean=y title="Acoustic Impedance" barlabel=AI scalebar=y')

Plot('sand1',None,'box x0=7.5 y0=7.5 label="shale 1" xt=0 yt=0')
Plot('shale',None,'box x0=7.5 y0=5.9 label="sand"    xt=0 yt=0')
Plot('sand2',None,'box x0=7.5 y0=4.3 label="shale 2" xt=0 yt=0')

for freq in (5,10,15):
    seismic = 'seismic%d' % freq
    Flow(seismic,'ai','ai2refl | ricker1 frequency=%d' % freq)
    Plot(seismic,'grey color=lb title="%d Hz wavelet" ' % freq)
    Result(seismic,seismic+' sand1 shale sand2','Overlay')

End()