```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()```