#define FIGDIR ../Fig
#define NEEDF90 zsnap.demo seismo before novar xquizsyntB windowxx

#include "../../Cakerules/SEP.martin.rules"
/* beware of cgstep, large dynamic range data have to be rescaled ! */

#define BINDIR ../../Bin/MTYPE
#define LIBDIR ../../Lib/MTYPE
#define FIGLIST wavelet2N filt2N diffspec zsnap.demo seismo before novar  xquizsyntB windowxx

#include <SEP.defs>
#define TRACEAXES label1="Time [sec]" label2="Displacement"

#define SCALEDATA |Scale dscale=0.00001

#define SHOTLIST 100 120 140 160 180 200 220 240 260 280
#define ALLZ [[ sub -i X Zshot.X  SHOTLIST ]]
#define ALLX [[ sub -i X Hshot.X  SHOTLIST ]]

#define SEISMO useindex=1

/* Shot variations */

#define VAR1 | Bandpass flo=5 fhi=30
#define VAR2 | Scale dscale=2
#define VAR3 | Bandpass flo=5 fhi=15 | Scale dscale=0.7
#define VAR4 | Conv  filtin=static10.H
#define VAR5 | Bandpass flo=10 fhi=25
#define VAR6 | Conv  filtin=static17.H | Scale dscale=0.1
#define VAR7 | Bandpass flo=20 fhi=50
#define VAR8 | Conv  filtin=static5.H | Bandpass flo=20 fhi=50
#define VAR9 | Noise var=0.333 


midpoint-offset&:
	Cubify axis3=CDP-X axis2=SOFFSET <Zzall.H attrfile=Zzall.hdr.H >jcpdzz.H


/* tests:
	2 spikes shifted belongin to 2 shots: testsimple
				add Noise   : testsimpleN
	2 wavelets shifted                  : testwavelet2
	 			add Noise   : testwavelet2N
*/

/* this is the 2 wavelet filter example */
#define THPLOTLABEL label1="Time Lag" label2="Amplitude"
#define WIGGLELABEL label1="Time [sec]" label2="Trace Number" o2num=1 d2num=1

FIGDIR/wavelet2N.v : wavelet2N.H wavelet2N.eq.H 
	cp wavelet2N.H jw.H
	echo "d2=1 o2=0" >> jw.H
	cp wavelet2N.eq.H jw.eq.H
	echo "d2=1 o2=0" >> jw.eq.H
	Wiggle pclip=100 <jw.H title="Input" WIGGLELABEL out=j1.v NULL
	Wiggle pclip=100 <jw.eq.H title="Equalized" WIGGLELABEL out=j2.v NULL
	vp_SideBySideIso j1.v j2.v >FIGDIR/wavelet2N.v	

#define CLIPVAL `Window f3=39 n3=1  <wavelet2N.filt.H |Clip pclip=100|Get clip parform=n`

FIGDIR/filt2N.v: wavelet2N.filt.H  
	Window f3=39 n3=1  <wavelet2N.filt.H | Dots pclip=100 title="Filter After 40 Iterations" THPLOTLABEL out=FIGDIR/filt2N.v NULL
	Dotmovie.pl wavelet2N.filt.H j3.v 40 CLIPVAL  out=j3.v NULL
	vppen scale=0.5 <j3.v >FIGDIR/filt2N.v3

#define FILTLABEL label1="Time Lag" label2="Filter No."
FIGDIR/filt%.v: filter.%.H
	Window f3=39 n3=1 <filter.%.H >jfil.H
	Thplot < jfil.H title="% Filter" FILTLABEL out=j1.v NULL alpha=45
	Taplot < jfil.H pclip=90 | Ta2vplot title=" " FILTLABEL out=j2.v NULL
	vppen big=n  <j1.v |vppen align=lb |vppen xsize=4 ysize=3 >j3.v
	vppen <j2.v  xsize=4 ysize=3 >j4.v
	vppen grid=-1 gridnum=2,1  gridsize=4.4,3.4 j3.v j4.v | vppen scale=0.5 |vppen align=lb >FIGDIR/filt%.v
	Thplot < filter.%.H title="% Filter" FILTLABEL out=FIGDIR/filt%.v3 NULL
	
FIGDIR/filtras.v: filter.Xx.H filter.Xz.H filter.Zx.H filter.Zz.H
	Window f3=39 n3=1 <filter.Xx.H |Taplot pclip=90 >jfilXx.H
	Window f3=39 n3=1 <filter.Xz.H |Taplot pclip=90 >jfilXz.H
	Window f3=39 n3=1 <filter.Zx.H |Taplot pclip=90 >jfilZx.H
	Window f3=39 n3=1 <filter.Zz.H |Taplot pclip=90 >jfilZz.H
	Merge axis=2 jfilXx.H jfilXz.H >j1.H
	Merge axis=2 jfilZx.H jfilZz.H >j2.H
	Merge axis=1 j1.H j2.H  >j3.H
	Ta2vplot <j3.H title="Equalization Filters" wanttitle=n wantlabel=n wantaxis=n  out=FIGDIR/filtras.v NULL

FIGDIR/filtwig.v: filter.Xx.H filter.Xz.H filter.Zx.H filter.Zz.H
	Window f3=39 n3=1 <filter.Xx.H |Thplot  out=jfilXx.v NULL
	Window f3=39 n3=1 <filter.Xz.H |Thplot  out=jfilXz.v NULL
	Window f3=39 n3=1 <filter.Zx.H |Thplot  out=jfilZx.v NULL
	Window f3=39 n3=1 <filter.Zz.H |Thplot  out=jfilZz.v NULL
	vppen big=n  <jfilXx.v |vppen align=lb |vppen xsize=4 ysize=4 >j1.v
	vppen big=n  <jfilXz.v |vppen align=lb |vppen xsize=4 ysize=4 >j2.v
	vppen big=n  <jfilZx.v |vppen align=lb |vppen xsize=4 ysize=4 >j3.v
	vppen big=n  <jfilZz.v |vppen align=lb |vppen xsize=4 ysize=4 >j4.v
	vppen grid=-1 gridnum=2,2  gridsize=4.4,4.4 j1.v j2.v j3.v j4.v | vppen scale=0.5 |vppen align=lb >FIGDIR/filtwig.v
	rm jfilXx.v jfilXz.v jfilZx.v jfilZz.v j1.v j2.v  j3.v j4.v

#define SEISMOLABEL label1="Time [sec]" label2="Distance [km]"
FIGDIR/seismo.v: Zzseis.200.t
	Taplot <Zzseis.200.t pclip=95 |Ta2vplot title="Zz" SEISMOLABEL out=Zzseis.200.v NULL
	vp_annotate <Zzseis.200.v text=seismo.ann  batch=y >FIGDIR/seismo.v



FIGDIR/zsnap.demo.v3:  zsnap.demo.T 
	Cubeplot <zsnap.demo.T n1pix=400 n2pix=400 flat=n movie=3 nframe=15 dframe=10 frame2=100 title="Z Comp.  Marmousi TI" label1="Depth [km]" label2="Distance [km]" label3="Time [sec]" out=FIGDIR/zsnap.demo.v3 NULL 

FIGDIR/zsnap.demo.v: zsnap.demo.T
	Cubeplot <zsnap.demo.T n1pix=400 n2pix=400 flat=n movie=0  frame3=170 frame2=100 title="Zz  Marmousi TI" label1="Depth [km]" label2="Distance [km]" label3="Time [sec]" out=FIGDIR/zsnap.demo.v NULL 


static%.H:
	Spike n1=20 k1=% nsp=1 >static%.H


$ INTERACTIVE
dot.action&: BINDIR/Dot  /* run Dot for a variety of filter length,and seed */
	xtpanel -file dot.panel -var DOT BINDIR/Dot

/* Dot product test of the equalizaton operator, error is ca 1e-8 */
dot&: BINDIR/Dot
	BINDIR/Dot --

/* THIS IS A kludge, since CAKE, chokes on too complicated rules */

norm%.H: if exist %all.H and exist %all.hdr.H
        BINDIR/Count <%all.H attrfile=%all.hdr.H >norm%.H recifile=reci%.H scdfile=scd%.H

norm%.H: if not exist %all.H and not exist %all.hdr.H
	cake %all.H %all.hdr.H
	cake norm%.H
	
/* END of kludge */

FIGDIR/before.v: BINDIR/Kill   normXx.H normXz.H normZx.H normZz.H 
	BINDIR/Kill < normXx.H >normXx.K.H attrin=Xxall.hdr.H attrout=Xxall.hdr.K.H
	BINDIR/Kill < normXz.H >normXz.K.H attrin=Xzall.hdr.H attrout=Xzall.hdr.K.H
	BINDIR/Kill < normZx.H >normZx.K.H attrin=Zxall.hdr.H attrout=Zxall.hdr.K.H
	BINDIR/Kill < normZz.H >normZz.K.H attrin=Zzall.hdr.H attrout=Zzall.hdr.K.H
	Window < normXx.K.H n2=90 |Taplot pclip=90 >j1.H
	Window < normXz.K.H n2=90 |Taplot pclip=90 >j2.H
	Window < normZx.K.H n2=90 |Taplot pclip=90 >j3.H 
	Window < normZz.K.H n2=90 |Taplot pclip=90 >j4.H
	Merge j1.H j2.H axis=2 >j11.H
	Merge j3.H j4.H axis=2 >j22.H
	Merge j11.H j22.H axis=1 > jj.H
	<jj.H Ta2vplot title="Before Equalization" wanttitle=n wantaxis=n wantlabel=n NULL out=FIGDIR/before.v

	
FIGDIR/after.v: BINDIR/Kill  normXx.eq.H normXz.eq.H normZx.eq.H normZz.eq.H  
	BINDIR/Kill < normXx.eq.H >normXx.eq.K.H attrin=Xxall.hdr.H attrout=Xxall.hdr.K.H
	BINDIR/Kill < normXz.eq.H >normXz.eq.K.H attrin=Xzall.hdr.H attrout=Xzall.hdr.K.H
	BINDIR/Kill < normZx.eq.H >normZx.eq.K.H attrin=Zxall.hdr.H attrout=Zxall.hdr.K.H
	BINDIR/Kill < normZz.eq.H >normZz.eq.K.H attrin=Zzall.hdr.H attrout=Zzall.hdr.K.H
	Window < normXx.eq.K.H n2=90 |Taplot pclip=90 >j1.H
	Window < normXz.eq.K.H n2=90 |Taplot pclip=90 >j2.H
	Window < normZx.eq.K.H n2=90 |Taplot pclip=90 >j3.H 
	Window < normZz.eq.K.H n2=90 |Taplot pclip=90 >j4.H
	Merge j1.H j2.H axis=2 >j11.H
	Merge j3.H j4.H axis=2 >j22.H
	Merge j11.H j22.H axis=1 > jj.H
	Ta2vplot <jj.H title="After Equalization" wanttitle=n wantaxis=n wantlabel=n NULL out=FIGDIR/after.v


FIGDIR/Zzbefore+after.v:  normZz.H normZz.eq.H  BINDIR/Kill
	BINDIR/Kill < normZz.H >normZz.K.H attrin=Zzall.hdr.H attrout=Zzall.hdr.K.H
	BINDIR/Kill < normZz.eq.H >normZz.eq.K.H attrin=Zzall.hdr.H attrout=Zzall.hdr.K.H
	Merge axis=2 normZz.K.H normZz.eq.K.H >jcomp.H
	Taplot <jcomp.H |Ta2vplot out=FIGDIR/Zzbefore+after.v NULL

FIGDIR/Xxbefore+after.v: normXx.H normXx.eq.H BINDIR/Kill
	BINDIR/Kill < normXx.H >normXx.K.H attrin=Xxall.hdr.H attrout=Xxall.hdr.K.H
	BINDIR/Kill < normXx.eq.H >normXx.eq.K.H attrin=Zzall.hdr.H attrout=Xxall.hdr.K.H
	Merge axis=2 normXx.K.H normXx.eq.K.H >jcomp.H
	Taplot <jcomp.H |Ta2vplot out=FIGDIR/Xxbefore+after.v NULL

FIGDIR/Zxbefore+after.v:  normZx.H normZx.eq.H BINDIR/Kill
	BINDIR/Kill < normZx.H >normZx.K.H attrin=Zxall.hdr.H attrout=Zxall.hdr.K.H
	BINDIR/Kill < normZx.eq.H >normZx.eq.K.H attrin=Zxall.hdr.H attrout=Zxall.hdr.K.H
	Merge axis=2 normZx.K.H normZx.eq.K.H >jcomp.H
	Taplot <jcomp.H |Ta2vplot out=FIGDIR/Zxbefore+after.v NULL

FIGDIR/Xzbefore+after.v:  normXz.H normXz.eq.H BINDIR/Kill
	BINDIR/Kill < normXz.H >normXz.K.H attrin=Xzall.hdr.H attrout=Xzall.hdr.K.H
	BINDIR/Kill < normXz.eq.H >normXz.eq.K.H attrin=Xzall.hdr.H attrout=Xzall.hdr.K.H
	Merge axis=2 normXz.K.H normXz.eq.K.H >jcomp.H
	Taplot <jcomp.H |Ta2vplot out=FIGDIR/Xzbefore+after.v NULL

$ Xquizit tests
FIGDIR/xquizsynt.v: BINDIR/Xquizitav BINDIR/Kill
	< normZz.H BINDIR/Kill  attrin=Zzall.hdr.H  >normZz.K.H 
	< reciZz.H BINDIR/Kill  attrin=Zzall.hdr.H  >reciZz.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.zz.H filt=xqfilt.zz.H <normZz.K.H in2=reciZz.K.H >xq.normZz.eq.K.H out2=xq.reciZz.eq.K.H
	< normXx.H BINDIR/Kill  attrin=Xxall.hdr.H  >normXx.K.H 
	< reciXx.H BINDIR/Kill  attrin=Xxall.hdr.H  >reciXx.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.xx.H filt=xqfilt.xx.H <normXx.B.K.H in2=reciXx.K.H >xq.normXx.eq.K.H out2=xq.reciXx.eq.K.H
	< normXz.H BINDIR/Kill  attrin=Xzall.hdr.H  >normXz.K.H 
	< reciXz.H BINDIR/Kill  attrin=Xzall.hdr.H  >reciXz.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.xz.H filt=xqfilt.xz.H <normXz.K.H in2=reciXz.K.H >xq.normXz.eq.K.H out2=xq.reciXz.eq.K.H
	< normZx.H BINDIR/Kill  attrin=Zxall.hdr.H  >normZx.K.H 
	< reciZx.H BINDIR/Kill  attrin=Zxall.hdr.H  >reciZx.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.zx.H filt=xqfilt.zx.H <normZx.K.H in2=reciZx.K.H >xq.normZx.eq.K.H out2=xq.reciZx.eq.K.H
	Window < xq.normXx.eq.K.H n2=90 |Taplot pclip=90 >j1.H
	Window < xq.normXz.eq.K.H n2=90 |Taplot pclip=90 >j2.H
	Window < xq.normZx.eq.K.H n2=90 |Taplot pclip=90 >j3.H
	Window < xq.normZz.eq.K.H n2=90 |Taplot pclip=90 >j4.H
	Merge j1.H j2.H axis=2 >j11.H
	Merge j3.H j4.H axis=2 >j22.H
	Merge j11.H j22.H axis=1 > jj.H
	Ta2vplot <jj.H title="Tracepair Equal." wanttitle=n wantaxis=n wantlabel=n NULL out=FIGDIR/xquizsynt.v

FIGDIR/xquizsyntB.v: BINDIR/Xquizitav BINDIR/Kill all10shots Zzall.hdr.H normZz.H normXx.H Xxall.hdr.H normXz.H reciXz.H Xzall.hdr.H normZx.H reciZx.H Zxall.hdr.H
	< normZz.H BINDIR/Kill  attrin=Zzall.hdr.H | Balance  >normZz.B.K.H 
	< reciZz.H BINDIR/Kill  attrin=Zzall.hdr.H | Balance  >reciZz.B.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.zz.H filt=xqfilt.zz.H <normZz.B.K.H in2=reciZz.B.K.H >xq.normZz.eq.B.K.H out2=xq.reciZz.eq.B.K.H
	< normXx.H BINDIR/Kill  attrin=Xxall.hdr.H | Balance  >normXx.B.K.H 
	< reciXx.H BINDIR/Kill  attrin=Xxall.hdr.H | Balance  >reciXx.B.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.xx.H filt=xqfilt.xx.H <normXx.B.K.H in2=reciXx.B.K.H >xq.normXx.eq.B.K.H out2=xq.reciXx.eq.B.K.H
	< normXz.H BINDIR/Kill  attrin=Xzall.hdr.H | Balance  >normXz.B.K.H 
	< reciXz.H BINDIR/Kill  attrin=Xzall.hdr.H | Balance  >reciXz.B.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.xz.H filt=xqfilt.xz.H <normXz.B.K.H in2=reciXz.B.K.H >xq.normXz.eq.B.K.H out2=xq.reciXz.eq.B.K.H
	< normZx.H BINDIR/Kill  attrin=Zxall.hdr.H | Balance  >normZx.B.K.H 
	< reciZx.H BINDIR/Kill  attrin=Zxall.hdr.H | Balance  >reciZx.B.K.H
	BINDIR/Xquizitav what=1  niter=5 bestfilt=xqbestfilt.zx.H filt=xqfilt.zx.H <normZx.B.K.H in2=reciZx.B.K.H >xq.normZx.eq.B.K.H out2=xq.reciZx.eq.B.K.H
	Window < xq.normXx.eq.B.K.H n2=90 |Taplot pclip=90 >j1.H
	Window < xq.normXz.eq.B.K.H n2=90 |Taplot pclip=90 >j2.H
	Window < xq.normZx.eq.B.K.H n2=90 |Taplot pclip=90 >j3.H
	Window < xq.normZz.eq.B.K.H n2=90 |Taplot pclip=90 >j4.H
	Merge j1.H j2.H axis=2 >j11.H
	Merge j3.H j4.H axis=2 >j22.H
	Merge j11.H j22.H axis=1 > jj.H
	Ta2vplot <jj.H title="Tracepair Equal." wanttitle=n wantaxis=n wantlabel=n NULL out=FIGDIR/xquizsyntB.v

xq2:      normZz.H reciZz.H
	< normZz.H Balance >normZz.B.H
	< reciZz.H Balance >reciZz.B.H
	BINDIR/Xquizits what=1  bestfilt=xqbestfilt.zz.H filt=xqfilt.zz.H <normZz.B.H in2=reciZz.B.H >xq.normZz.eq.B.H out2=xq.reciZz.eq.B.H
	BINDIR/Kill < normZz.B.H >normZz.B.K.H attrin=Zzall.hdr.H 
	BINDIR/Kill < xq.normZz.eq.B.H >xq.normZz.eq.B.K.H attrin=Zzall.hdr.H 

xq.zz.H:
	BINDIR/Philtrei < Zzall.H bestfilt=xqbestfilt.xx.H > xq.zz.H
	
xqfilter.Zz.K.H:
        BINDIR/Kill < normZz.H >normZz.K.H attrin=Zzall.hdr.H attrout=Zzall.hdr.K.H
        BINDIR/Kill < reciZz.H >reciZz.K.H attrin=Zzall.hdr.H 
	BINDIR/Xquiziti what=1 bestfilt=xqbest.xx.K.H filt=xqfilt.xx.K.H <normZz.K.H in2=reciZz.K.H >xq.normZz.K.eq.H out2=xq.reciZz.K.eq.H


$ EQUALIZE

filter.Zz.H equalZz&: Zzall.hdr.H Zzall.H 
        Cp < Zzall.hdr.H >jZzall.hdr.H
        BINDIR/Count <Zzall.H attrfile=jZzall.hdr.H >normZz.H recifile=reciZz.H scdfile=scdZz.H
	time BINDIR/EqualizePre5a <normZz.H attrfile=jZzall.hdr.H recifile=reciZz.H scdfile=scdZz.H restart=d >normZz.eq.H filtermovie=filter.Zz.H par=parfile.equal niter=40 damp=0.01

filter.Xx.H equalXx&:    Zzall.hdr.H Xxall.H 
        Cp < Xxall.hdr.H >jXxall.hdr.H
        BINDIR/Count <Xxall.H attrfile=jXxall.hdr.H >normXx.H recifile=reciXx.H scdfile=scdXx.H
	time BINDIR/EqualizePre5a <normXx.H attrfile=jXxall.hdr.H recifile=reciXx.H scdfile=scdXx.H restart=d >normXx.eq.H filtermovie=filter.Xx.H par=parfile.equal niter=40 damp=0.01
	

filter.Zx.H equalZxXz&:     Zxall.H Xzall.H  Zxall.hdr.H Xzall.hdr.H  
	Merge space=n axis=3 Zxall.H Xzall.H >ZxXzall.H
	Merge space=n axis=3 Zxall.hdr.H Xzall.hdr.H >ZxXzall.hdr.H
        Cp < ZxXzall.hdr.H >jZxXzall.hdr.H
        BINDIR/Count6 <ZxXzall.H attrfile=jZxXzall.hdr.H >normZxXz.H recifile=reciZxXz.H scdfile=scdZxXz.H
	Window f3=0  n3=10 < normZxXz.H >normZx.H
	Window f3=10 n3=10 < normZxXz.H >normXz.H
	Window f3=0  n3=10 < reciZxXz.H >reciZx.H
	Window f3=10 n3=10 < reciZxXz.H >reciXz.H
	Window f3=0  n3=10 < scdZxXz.H >scdZx.H
	Window f3=10 n3=10 < scdZxXz.H >scdXz.H
	Window f3=0  n3=10 < jZxXzall.hdr.H >jZxall.hdr.H
	Window f3=10 n3=10 < jZxXzall.hdr.H >jXzall.hdr.H
$ now we have everything sorted and split up again, proceed with estimation	
	time BINDIR/EqualizePre5a <normZx.H attrfile=jZxall.hdr.H recifile=reciZx.H scdfile=scdZx.H restart=d >normZx.eq.H filtermovie=filter.Zx.H par=parfile.equal niter=40 damp=0.01
$	time BINDIR/EqualizePre5a <normXz.H attrfile=jXzall.hdr.H recifile=reciXz.H scdfile=scdXz.H restart=d >normXz.eq.H filtermovie=filter.Xz.H par=parfile.equal niter=40 damp=0.01

filter.Xz.H equalXzZx&:     Xzall.H Zxall.H  Xzall.hdr.H Zxall.hdr.H 
	Merge space=n axis=3 Xzall.H Zxall.H >XzZxall.H
	Merge space=n axis=3 Xzall.hdr.H Zxall.hdr.H  >XzZxall.hdr.H
        Cp < XzZxall.hdr.H >jXzZxall.hdr.H
        BINDIR/Count6 <XzZxall.H attrfile=jXzZxall.hdr.H >normXzZx.H recifile=reciXzZx.H scdfile=scdXzZx.H
	Window f3=0  n3=10 < normXzZx.H >normXz.H
	Window f3=10 n3=10 < normXzZx.H >normZx.H
	Window f3=0  n3=10 < reciXzZx.H >reciXz.H
	Window f3=10 n3=10 < reciXzZx.H >reciZx.H
	Window f3=0  n3=10 < scdXzZx.H >scdXz.H
	Window f3=10 n3=10 < scdXzZx.H >scdZx.H
	Window f3=0  n3=10 < jXzZxall.hdr.H >jXzall.hdr.H
	Window f3=10 n3=10 < jXzZxall.hdr.H >jZxall.hdr.H
$ now we have everything sorted and split up again, proceed with estimation	
	time BINDIR/EqualizePre5a <normXz.H attrfile=jXzall.hdr.H recifile=reciXz.H scdfile=scdXz.H restart=d >normXz.eq.H filtermovie=filter.Xz.H par=parfile.equal niter=40 damp=0.01
$	time BINDIR/EqualizePre5a <normZx.H attrfile=jZxall.hdr.H recifile=reciZx.H scdfile=scdZx.H restart=d >normZx.eq.H filtermovie=filter.Zx.H par=parfile.equal niter=40 damp=0.01
	
	
filter.ZxXz.H equalZxXz&:     Zxall.H Xzall.H  Zxall.hdr.H Xzall.hdr.H 
	Merge space=n axis=3 Zxall.H Xzall.H >ZxXzall.H
	Merge space=n axis=3 Zxall.hdr.H Xzall.hdr.H >ZxXzall.hdr.H
        BINDIR/Count6 <ZxXzall.H attrfile=jZxXzall.hdr.H >normZxXz.H recifile=reciZxXz.H scdfile=scdZxXz.H
	time BINDIR/EqualizePre5a <normZxXz.H attrfile=Zxall.hdr.H recifile=reciZxXz.H scdfile=scdZxXz.H restart=d >normZxXz.eq.H filtermovie=filter.ZxXz.H par=parfile.equal niter=40 damp=0.01



filter.Xz.H equalXzZx&:     Xzall.H Zxall.H  Xzall.hdr.H Zxall.hdr.H 

equalize&:  equalZz  equalXx equalZxXz equalXzZx


/* test if matching two spikes with noise is possible */

$ works fine,  finds the shift in the first iteration  
testsimple&: jsimple.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.simple niter=10 <jsimple.H >jequal.H

$ works with 30 percent noise, gets the shift ok, with the random wavelet
testsimpleN&: jsimple.H BINDIR/EqualizePre5 
	Noise <jsimple.H var=0.0333 >jnoise.H
        BINDIR/EqualizePre5 restart=n par=parfile.simple niter=10 <jnoise.H >jequal.H


testspike&: jtwospikes.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.spike niter=10 <jtwospikes.H >jequal.H

testsimple4&: jfourspikes.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.spike niter=10 <jfourspikes.H >jequal.H

testwavelet2&: wavelet2.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.spike niter=40 damp=0.01 <wavelet2.H >jequal.H


$ 2 wavelets with 30 percent noise
testwavelet2N& wavelet2N.filt.H wavelet2N.eq.H: wavelet2.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.spike niter=40 damp=0.01 <wavelet2N.H >wavelet2N.eq.H filtermovie=wavelet2N.filt.H constrain=1


testwavelet4&: jwavelet4.H BINDIR/EqualizePre5 
        BINDIR/EqualizePre5 restart=n par=parfile.spike niter=100 <jwavelet4.H >jequal.H


testspike6&: jtwospikes.H BINDIR/EqualizePre6 
        BINDIR/EqualizePre6 restart=n par=parfile.spike niter=10 <jtwospikes.H >jequal.H
        Graph <filter1  |Tube

jsimple.H:
        echo " 1 1 0 0 10 0 1 0 2 0 0 0 1 0 0 1 0 0 " >junk
        echo " 2 1 0 0 20 0 2 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        atoF <junk >jhdr.H
        echo "n1=18 n2=2 n3=1" >>jhdr.H
        echo "hdrkey1=SHOT hdrkey2=SEQNO" >>jhdr.H
        echo "hdrkey3=junk hdrkey4=junk hdrkey5=OFFSET hdrkey6=junk">>jhdr.H
        echo "hdrkey7=SHT-X hdrkey8=SHT-Y hdrkey9=REC-X hdrkey10=REC-Y" >>jhdr.H
        echo "hdrkey11=junk hdrkey12=junk">>jhdr.H
        echo "hdrkey13=SCOMPX hdrkey14=SCOMPY hdrkey15=SCOMPZ" >>jhdr.H
        echo "hdrkey16=RCOMPX hdrkey17=RCOMPY hdrkey18=RCOMPZ" >>jhdr.H
        Spike  o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=25 n2=2 n3=1 k1=10,15 k2=1,2 k3=1,1 nsp=2 >jsimple.H
        echo "attrfile=jhdr.H" >>jsimple.H
        echo "o2=1.5 d2=0.5 o3=0.01 d3=0.1" >>jhdr.H

jtwospikes.H:
        echo " 1 1 0 0 10 0 1 0 2 0 0 0 1 0 0 1 0 0 " >junk
        echo " 2 1 0 0 20 0 2 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        atoF <junk >jhdr.H
        echo "n1=18 n2=2 n3=1" >>jhdr.H
        echo "hdrkey1=SHOT hdrkey2=SEQNO" >>jhdr.H
        echo "hdrkey3=junk hdrkey4=junk hdrkey5=OFFSET hdrkey6=junk">>jhdr.H
        echo "hdrkey7=SHT-X hdrkey8=SHT-Y hdrkey9=REC-X hdrkey10=REC-Y" >>jhdr.H
        echo "hdrkey11=junk hdrkey12=junk">>jhdr.H
        echo "hdrkey13=SCOMPX hdrkey14=SCOMPY hdrkey15=SCOMPZ" >>jhdr.H
        echo "hdrkey16=RCOMPX hdrkey17=RCOMPY hdrkey18=RCOMPZ" >>jhdr.H
        Spike  o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250 n2=2 n3=1 k1=37,39,38 k2=1,2,1 k3=1,1,1 nsp=3 | Noise var=0.0333  >jtwospikes.H
        echo "attrfile=jhdr.H" >>jtwospikes.H
        echo "o2=1.5 d2=0.5 o3=0.01 d3=0.1" >>jhdr.H

$ these are two independent shots, so one is constraint the other is not.
$ EqualizePre should check if on trace pair a particular shot has only one
$ trace pair, then itshould constrain, otherwise all are linked.
jfourspikes.H:
        echo " 1 1 0 0 10 0 1 0 2 0 0 0 1 0 0 1 0 0 " >junk
        echo " 2 1 0 0 20 0 2 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        echo " 3 1 0 0 10 0 3 0 4 0 0 0 1 0 0 1 0 0 " >>junk
        echo " 4 1 0 0 20 0 4 0 3 0 0 0 1 0 0 1 0 0 " >>junk
        atoF <junk >jhdr.H
        echo "n1=18 n2=4 n3=1" >>jhdr.H
        echo "hdrkey1=SHOT hdrkey2=SEQNO" >>jhdr.H
        echo "hdrkey3=junk hdrkey4=junk hdrkey5=OFFSET hdrkey6=junk">>jhdr.H
        echo "hdrkey7=SHT-X hdrkey8=SHT-Y hdrkey9=REC-X hdrkey10=REC-Y" >>jhdr.H
        echo "hdrkey11=junk hdrkey12=junk">>jhdr.H
        echo "hdrkey13=SCOMPX hdrkey14=SCOMPY hdrkey15=SCOMPZ" >>jhdr.H
        echo "hdrkey16=RCOMPX hdrkey17=RCOMPY hdrkey18=RCOMPZ" >>jhdr.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=4  k1=37,47,100,115 k2=1,2,3,4 nsp=4      >jfourspikes.H
        echo "attrfile=jhdr.H" >>jfourspikes.H
        echo "o2=1.5 d2=0.5 o3=0.01 d3=0.1" >>jhdr.H

wavelet2N.H: wavelet2.H
	Noise <wavelet2.H var=0.00633 >wavelet2N.H

wavelet2.H wavelet2.H:
        echo " 1 1 0 0 10 0 1 0 2 0 0 0 1 0 0 1 0 0 " >junk
        echo " 2 1 0 0 20 0 2 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        atoF <junk >jhdr.H
        echo "n1=18 n2=2 n3=1" >>jhdr.H
        echo "hdrkey1=SHOT hdrkey2=SEQNO" >>jhdr.H
        echo "hdrkey3=junk hdrkey4=junk hdrkey5=OFFSET hdrkey6=junk">>jhdr.H
        echo "hdrkey7=SHT-X hdrkey8=SHT-Y hdrkey9=REC-X hdrkey10=REC-Y" >>jhdr.H
        echo "hdrkey11=junk hdrkey12=junk">>jhdr.H
        echo "hdrkey13=SCOMPX hdrkey14=SCOMPY hdrkey15=SCOMPZ" >>jhdr.H
        echo "hdrkey16=RCOMPX hdrkey17=RCOMPY hdrkey18=RCOMPZ" >>jhdr.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=37 k2=1 nsp=1  VAR1>j1.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=47 k2=1 nsp=1  VAR8>j2.H
	Merge j1.H j2.H axis=2 space=n>wavelet2.H
        echo "attrfile=jhdr.H" >>wavelet2.H
        echo "o2=1.5 d2=0.5 o3=0.01 d3=0.1" >>jhdr.H

jwavelet4.H:
        echo " 1 1 0 0 10 0 1 0 2 0 0 0 1 0 0 1 0 0 " >junk
        echo " 2 1 0 0 20 0 2 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        echo " 1 1 0 0 10 0 1 0 3 0 0 0 1 0 0 1 0 0 " >>junk
        echo " 3 1 0 0 20 0 3 0 1 0 0 0 1 0 0 1 0 0 " >>junk
        atoF <junk >jhdr.H
        echo "n1=18 n2=4 n3=1" >>jhdr.H
        echo "hdrkey1=SHOT hdrkey2=SEQNO" >>jhdr.H
        echo "hdrkey3=junk hdrkey4=junk hdrkey5=OFFSET hdrkey6=junk">>jhdr.H
        echo "hdrkey7=SHT-X hdrkey8=SHT-Y hdrkey9=REC-X hdrkey10=REC-Y" >>jhdr.H
        echo "hdrkey11=junk hdrkey12=junk">>jhdr.H
        echo "hdrkey13=SCOMPX hdrkey14=SCOMPY hdrkey15=SCOMPZ" >>jhdr.H
        echo "hdrkey16=RCOMPX hdrkey17=RCOMPY hdrkey18=RCOMPZ" >>jhdr.H
$	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=4  k1=37,47,100,115 k2=1,2,3,4 nsp=4      >jwavelet4.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=37 k2=1 nsp=1  VAR1>j1.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=47 k2=1 nsp=1  VAR2>j2.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=100 k2=1 nsp=1 VAR1>j3.H
	Spike o2=1.5 d2=0.5 o3=0.01 d3=0.1 n1=250  n2=1  k1=115 k2=1 nsp=1 VAR3>j4.H
	Merge j1.H j2.H j3.H j4.H axis=2 space=n>jwavelet4.H
        echo "attrfile=jhdr.H" >>jwavelet4.H
        echo "o2=1.5 d2=0.5 o3=0.01 d3=0.1" >>jhdr.H


/* HERE IS THE MARMOUSI SYNTHETIC generation */

progs&: BINDIR/Makeelastic BINDIR/Makeaniso


/* create all the 10 shots with  variations if VAR1-9 are defined  */
all10shots&:  /*  ALLZ ALLX  */
	cake allmerge
	cake allheaders

allheaders&: [[ sub -i X Zheader.X SHOTLIST ]] [[ sub -i X Hheader.X SHOTLIST ]]

allmerge&: mergedata mergehdr


$ LOOK AT SHOTVARIATIONS

#define SPIKEARGS n1=256 k1=128 mag=1
#define SPECTRUM | Spectrum  
/* first shot has not variations, is the reference shot */
shotvar0.H:
	Spike n1=129 mag=0 d1=0.976562 >shotvar0.H

shotvar1.H: 
	Spike SPIKEARGS VAR1 SPECTRUM >shotvar1.H

shotvar2.H:  /* a spike in time = 1 in omega */
	Spike n1=129  d1=0.976562 mag=`Clip <shotvar1.H|Get clip parform=n` |Scale dscale=2 >shotvar2.H
/*	Spike SPIKEARGS VAR2 SPECTRUM >shotvar2.H*/

shotvar3.H: 
	Spike SPIKEARGS VAR3 SPECTRUM >shotvar3.H

shotvar4.H: 
	Spike SPIKEARGS VAR4 SPECTRUM >shotvar4.H

shotvar5.H: 
	Spike SPIKEARGS VAR5 SPECTRUM >shotvar5.H

shotvar6.H:  /* a spike in time = 1 in omega */
	Spike n1=129 d1=0.976562 mag=`Clip <shotvar1.H|Get clip parform=n` |Scale dscale=0.1 >shotvar6.H
/*	Spike SPIKEARGS VAR6 SPECTRUM >shotvar6.H*/

shotvar7.H: 
	Spike SPIKEARGS VAR7 SPECTRUM >shotvar7.H

shotvar8.H: 
	Spike SPIKEARGS VAR8 SPECTRUM >shotvar8.H

shotvar9.H: /* 1+ random in omega */
	Spike SPIKEARGS VAR9 SPECTRUM >shotvar9.H


shotvar.H: [[ sub -i X shotvarX.H 0 1 2 3 4 5 6 7 8 9 ]]
	Merge space=n axis=2 shotvar0.H shotvar1.H shotvar2.H shotvar3.H shotvar4.H shotvar5.H shotvar6.H shotvar7.H shotvar8.H shotvar9.H | Scale dscale=1 axis=1 >shotvar.H

#define THRESHHOLD 0.1

FIGDIR/diffspec.v: shotvar.H 
	echo "o2=1 d2=1 o2num=1 d2num=1 " >>shotvar.H
	Window n1=64 <shotvar.H |Taplot allpos=y pclip=100 |Ta2vplot title="Spectrum of Variations" label1="Frequency [Hz]" label2="Shot" out=FIGDIR/diffspec.v NULL

/* apply here the variations */

%.100.t: %.100 if not exist %.100
	Transp plane=12 <  %.100  >%.100.t

%.120.t: %.120 if not exist %.120
	Transp plane=12 <  %.120  VAR1 >%.120.t

%.140.t: %.140 if not exist %.140
	Transp plane=12 <  %.140  VAR2 >%.140.t

%.160.t: %.160 if not exist %.160
	Transp plane=12 <  %.160  VAR3 >%.160.t

%.180.t: %.180 static10.H if not exist %.180
	Transp plane=12 <  %.180  VAR4 >%.180.t

%.200.t: %.200 if not exist %.200
	Transp plane=12 <  %.200  VAR5 >%.200.t

%.220.t: %.220 static17.H if not exist %.220 
	Transp plane=12 <  %.220  VAR6 >%.220.t

%.240.t: %.240 if not exist %.240
	Transp plane=12 <  %.240  VAR7 >%.240.t

%.260.t: %.260 static5.H if not exist %.260
	Transp plane=12 <  %.260  VAR8 >%.260.t

%.280.t: %.280 if not exist %.280
	Transp plane=12 <  %.280  VAR9 >%.280.t


transpZdata&: Zxseis.100.t Zxseis.120.t Zxseis.140.t Zxseis.160.t Zxseis.180.t Zxseis.200.t Zxseis.220.t Zxseis.240.t Zxseis.260.t Zxseis.280.t Zzseis.100.t Zzseis.120.t Zzseis.140.t Zzseis.160.t Zzseis.180.t Zzseis.200.t Zzseis.220.t Zzseis.240.t Zzseis.260.t Zzseis.280.t Zzseis.100.t Zzseis.120.t Zzseis.140.t Zzseis.160.t Zzseis.180.t Zzseis.200.t Zzseis.220.t Zzseis.240.t Zzseis.260.t Zzseis.280.t Zzseis.100.t Zzseis.120.t Zzseis.140.t Zzseis.160.t Zzseis.180.t Zzseis.200.t Zzseis.220.t Zzseis.240.t Zzseis.260.t Zzseis.280.t

transpXdata&: Xxseis.100.t Xxseis.120.t Xxseis.140.t Xxseis.160.t Xxseis.180.t Xxseis.200.t Xxseis.220.t Xxseis.240.t Xxseis.260.t Xxseis.280.t Xzseis.100.t Xzseis.120.t Xzseis.140.t Xzseis.160.t Xzseis.180.t Xzseis.200.t Xzseis.220.t Xzseis.240.t Xzseis.260.t Xzseis.280.t Xzseis.100.t Xzseis.120.t Xzseis.140.t Xzseis.160.t Xzseis.180.t Xzseis.200.t Xzseis.220.t Xzseis.240.t Xzseis.260.t Xzseis.280.t Xzseis.100.t Xzseis.120.t Xzseis.140.t Xzseis.160.t Xzseis.180.t Xzseis.200.t Xzseis.220.t Xzseis.240.t Xzseis.260.t Xzseis.280.t


ZxXzall.H:
	Merge space=n axis=3 Zxall.H Xzall.H >ZxXzall.H

ZxXzall.hdr.H:
	Merge space=n axis=3 Zxall.hdr.H Xzall.hdr.H >ZxXzall.hdr.H


#define ZXDEPT Zxseis.100.t Zxseis.120.t Zxseis.140.t Zxseis.160.t Zxseis.180.t Zxseis.200.t Zxseis.220.t Zxseis.240.t Zxseis.260.t Zxseis.280.t

Zxall.H:  ZXDEPT
	Merge axis=3 space=n ZXDEPT  SCALEDATA >Zxall.H

#define ZZDEPT Zzseis.100.t Zzseis.120.t Zzseis.140.t Zzseis.160.t Zzseis.180.t Zzseis.200.t Zzseis.220.t Zzseis.240.t Zzseis.260.t Zzseis.280.t

Zzall.H:  ZZDEPT
	Merge axis=3 space=n ZZDEPT  SCALEDATA >Zzall.H

#define XXDEPT Xxseis.100.t Xxseis.120.t Xxseis.140.t Xxseis.160.t Xxseis.180.t Xxseis.200.t Xxseis.220.t Xxseis.240.t Xxseis.260.t Xxseis.280.t

Xxall.H:  XXDEPT
	Merge axis=3 space=n XXDEPT SCALEDATA >Xxall.H

#define XZDEPT Xzseis.100.t Xzseis.120.t Xzseis.140.t Xzseis.160.t Xzseis.180.t Xzseis.200.t Xzseis.220.t Xzseis.240.t Xzseis.260.t Xzseis.280.t

Xzall.H:  XZDEPT
	Merge axis=3 space=n  SCALEDATA >Xzall.H

/* merge and scale data to something close to 1 */
mergedata&: transpZdata transpXdata
	Merge axis=3 space=n Zxseis.100.t Zxseis.120.t Zxseis.140.t Zxseis.160.t Zxseis.180.t Zxseis.200.t Zxseis.220.t Zxseis.240.t Zxseis.260.t Zxseis.280.t  SCALEDATA >Zxall.H
	Merge axis=3 space=n Zzseis.100.t Zzseis.120.t Zzseis.140.t Zzseis.160.t Zzseis.180.t Zzseis.200.t Zzseis.220.t Zzseis.240.t Zzseis.260.t Zzseis.280.t  SCALEDATA >Zzall.H
	Merge axis=3 space=n Xxseis.100.t Xxseis.120.t Xxseis.140.t Xxseis.160.t Xxseis.180.t Xxseis.200.t Xxseis.220.t Xxseis.240.t Xxseis.260.t Xxseis.280.t  SCALEDATA >Xxall.H
	Merge axis=3 space=n Xzseis.100.t Xzseis.120.t Xzseis.140.t Xzseis.160.t Xzseis.180.t Xzseis.200.t Xzseis.220.t Xzseis.240.t Xzseis.260.t Xzseis.280.t  SCALEDATA >Xzall.H

#define ZXDEPHT Zxseis.100.hdr.H Zxseis.120.hdr.H Zxseis.140.hdr.H Zxseis.160.hdr.H Zxseis.180.hdr.H Zxseis.200.hdr.H Zxseis.220.hdr.H Zxseis.240.hdr.H Zxseis.260.hdr.H Zxseis.280.hdr.H
Zxall.hdr.H: ZXDEPHT
	Merge axis=3 space=n  ZXDEPHT >Zxall.hdr.H

#define ZZDEPHT Zzseis.100.hdr.H Zzseis.120.hdr.H Zzseis.140.hdr.H Zzseis.160.hdr.H Zzseis.180.hdr.H Zzseis.200.hdr.H Zzseis.220.hdr.H Zzseis.240.hdr.H Zzseis.260.hdr.H Zzseis.280.hdr.H
Zzall.hdr.H: ZZDEPHT
	Merge axis=3 space=n  ZZDEPHT >Zzall.hdr.H

#define XZDEPHT Xzseis.100.hdr.H Xzseis.120.hdr.H Xzseis.140.hdr.H Xzseis.160.hdr.H Xzseis.180.hdr.H Xzseis.200.hdr.H Xzseis.220.hdr.H Xzseis.240.hdr.H Xzseis.260.hdr.H Xzseis.280.hdr.H
Xzall.hdr.H: XZDEPHT
	Merge axis=3 space=n  XZDEPHT >Xzall.hdr.H

#define XXDEPHT Xxseis.100.hdr.H Xxseis.120.hdr.H Xxseis.140.hdr.H Xxseis.160.hdr.H Xxseis.180.hdr.H Xxseis.200.hdr.H Xxseis.220.hdr.H Xxseis.240.hdr.H Xxseis.260.hdr.H Xxseis.280.hdr.H
Xxall.hdr.H: XXDEPHT
	Merge axis=3 space=n  XXDEPHT >Xxall.hdr.H

mergehdr&:
	Merge axis=3 space=n Zxseis.100.hdr.H Zxseis.120.hdr.H Zxseis.140.hdr.H Zxseis.160.hdr.H Zxseis.180.hdr.H Zxseis.200.hdr.H Zxseis.220.hdr.H Zxseis.240.hdr.H Zxseis.260.hdr.H Zxseis.280.hdr.H >Zxall.hdr.H
	Merge axis=3 space=n Zzseis.100.hdr.H Zzseis.120.hdr.H Zzseis.140.hdr.H Zzseis.160.hdr.H Zzseis.180.hdr.H Zzseis.200.hdr.H Zzseis.220.hdr.H Zzseis.240.hdr.H Zzseis.260.hdr.H Zzseis.280.hdr.H >Zzall.hdr.H
	Merge axis=3 space=n Xxseis.100.hdr.H Xxseis.120.hdr.H Xxseis.140.hdr.H Xxseis.160.hdr.H Xxseis.180.hdr.H Xxseis.200.hdr.H Xxseis.220.hdr.H Xxseis.240.hdr.H Xxseis.260.hdr.H Xxseis.280.hdr.H >Xxall.hdr.H
	Merge axis=3 space=n Xzseis.100.hdr.H Xzseis.120.hdr.H Xzseis.140.hdr.H Xzseis.160.hdr.H Xzseis.180.hdr.H Xzseis.200.hdr.H Xzseis.220.hdr.H Xzseis.240.hdr.H Xzseis.260.hdr.H Xzseis.280.hdr.H >Xzall.hdr.H

show260&:
	Byte <Xxseis.260 |Ta2vplot title="Xx" out=Xxseis.260.v head=/dev/null
	Byte <Xzseis.260 |Ta2vplot title="Xz" out=Xzseis.260.v head=/dev/null
	Byte <Zxseis.260 |Ta2vplot title="Zx" out=Zxseis.260.v head=/dev/null
	Byte <Zzseis.260 |Ta2vplot title="Zz" out=Zzseis.260.v head=/dev/null
	vp_SideBySideAniso Xxseis.260.v Xzseis.260.v >jj1.v
	vp_SideBySideAniso Zxseis.260.v Zzseis.260.v >jj2.v
	vp_OverUnderAniso jj1.v jj2.v >panel.260.v



Zheader.%&: BINDIR/CreateHeader
	BINDIR/CreateHeader SEISMO < Zxseis.%.t  > Zxseis.%.hdr.H shot=%  shtx=% recx=1 scompx=0 scompy=0 scompz=1 rcompx=1 rcompy=0 rcompz=0 
$	BINDIR/CreateHeader SEISMO < Zyseis.%.t  > Zyseis.%.hdr.H shot=%  shtx=% recx=1 scompx=0 scompy=0 scompz=1 rcompx=0 rcompy=1 rcompz=0 
	BINDIR/CreateHeader SEISMO < Zzseis.%.t  > Zzseis.%.hdr.H shot=%  shtx=% recx=1 scompx=0 scompy=0 scompz=1 rcompx=0 rcompy=0 rcompz=1 

Hheader.%& Xheader.%&: BINDIR/CreateHeader
	BINDIR/CreateHeader SEISMO < Xxseis.%.t  > Xxseis.%.hdr.H shot=%  shtx=% recx=1 scompx=1 scompy=0 scompz=0 rcompx=1 rcompy=0 rcompz=0 
$	BINDIR/CreateHeader SEISMO < Xyseis.%.t  > Xyseis.%.hdr.H shot=%  shtx=% recx=1 scompx=1 scompy=0 scompz=0 rcompx=0 rcompy=1 rcompz=0 
	BINDIR/CreateHeader SEISMO < Xzseis.%.t  > Xzseis.%.hdr.H shot=%  shtx=% recx=1 scompx=1 scompy=0 scompz=0 rcompx=0 rcompy=0 rcompz=1 


Ztest.%&: tismall.H wavelet.Z.H  BINDIR/Ultiall
	time Ultiall  par=parfile type=4 ncomp=3 modeling=1 src_type=0 nsrc=1 src_1=% moduli=tismall.H <wavelet.Z.H >junk snap_i=100 surf_type=0 tmax=1


Zshot.%& Zxseis.% Zzseis.%: tismall.H wavelet.Z.H  BINDIR/Ultiall
        time Ultiall  par=parfile type=4 ncomp=3 modeling=1 src_type=0 nsrc=1 src_1=% moduli=tismall.H xseis=Zxseis.% yseis=Zyseis.% zseis=Zzseis.% <wavelet.Z.H >junk snap_i=100


Hshot.%& Xxseis.% Zzseis.%: tismall.H wavelet.X.H  BINDIR/Ultiall
        time Ultiall  par=parfile type=4 ncomp=3 modeling=1 src_type=0 nsrc=1 src_1=% moduli=tismall.H xseis=Xxseis.% yseis=Xyseis.% zseis=Xzseis.% <wavelet.X.H >junk


#define DEMO 200
Zcube.demo zsnap.demo&: tismall.H wavelet.Z.H BINDIR/Ultiall
        time Ultiall  par=parfile type=4 ncomp=3 modeling=1 src_type=0 nsrc=1 src_1=DEMO moduli=tismall.H xseis=Zxseis.demo.DEMO yseis=Zyseis.demo.DEMO zseis=Zzseis.demo.DEMO <wavelet.Z.H >junk snap_i=10 zsnap=zsnap.demo xsnap=xsnap.demo ysnap=ysnap.demo


elastic_vertical1&:  /* tismall.H wavelet.vert.H */ BINDIR/Ultiall
        time Ultiall  type=4 ncomp=3 modeling=1 src_type=0 nsrc=1 src_1=104  moduli=tismall.H par=parfile  xsnap=xsnap1 zsnap=zsnap1 ysnap=ysnap1 xseis=xseis1 yseis1 zseis=zseis1 <wavelet.vert.H >junk


zsnap.demo.T: zsnap.demo
	Taplot gainpanel=all < zsnap.demo >zsnap.demo.T


/* no source behavior variation at all */
novar&:
	Merge  Xxseis.100 Xxseis.120 Xxseis.140 Xxseis.160 Xxseis.180 Xxseis.200 Xxseis.220 Xxseis.240 Xxseis.260 Xxseis.280 axis=3 space=n >jXx
	Transp plane=12 <jXx >jXx.t
	Merge  Xzseis.100 Xzseis.120 Xzseis.140 Xzseis.160 Xzseis.180 Xzseis.200 Xzseis.220 Xzseis.240 Xzseis.260 Xzseis.280 axis=3 space=n >jXz
	Transp plane=12 <jXz >jXz.t
	Merge  Zxseis.100 Zxseis.120 Zxseis.140 Zxseis.160 Zxseis.180 Zxseis.200 Zxseis.220 Zxseis.240 Zxseis.260 Zxseis.280 axis=3 space=n >jZx
	Transp plane=12 <jZx >jZx.t
	Merge  Zzseis.100 Zzseis.120 Zzseis.140 Zzseis.160 Zzseis.180 Zzseis.200 Zzseis.220 Zzseis.240 Zzseis.260 Zzseis.280 axis=3 space=n >jZz
	Transp plane=12 <jZz >jZz.t

novar2&: novar 
	BINDIR/Count <jXx.t attrfile=Xxall.hdr.H >jnormxx.H recifile=jrecixx.H scdfile=jscdxx.H
	BINDIR/Count <jXz.t attrfile=Xzall.hdr.H >jnormxz.H recifile=jrecixz.H scdfile=jscdxz.H
	BINDIR/Count <jZx.t attrfile=Zxall.hdr.H >jnormzx.H recifile=jrecizx.H scdfile=jscdzx.H
	BINDIR/Count <jZz.t attrfile=Zzall.hdr.H >jnormzz.H recifile=jrecizz.H scdfile=jscdzz.H
	BINDIR/Kill < jnormxx.H >jnormxx.K.H attrin=Xxall.hdr.H attrout=jxxall.hdr.K.H
	BINDIR/Kill < jnormxz.H >jnormxz.K.H attrin=Xzall.hdr.H attrout=jxzall.hdr.K.H
	BINDIR/Kill < jnormzx.H >jnormzx.K.H attrin=Zxall.hdr.H attrout=jzxall.hdr.K.H
	BINDIR/Kill < jnormzz.H >jnormzz.K.H attrin=Zzall.hdr.H attrout=jzzall.hdr.K.H
	echo jnorm .. K.H contains the shots with no variations
	
FIGDIR/novar.v novar3&: novar2
	Merge jnormxx.K.H jnormxz.K.H axis=2 >j1.H
	Merge jnormzx.K.H jnormzz.K.H axis=2 >j2.H
	Merge j1.H j2.H axis=1 >j3.H
	Taplot pclip=90 <j3.H |Ta2vplot title="No Variations" out=FIGDIR/novar.v NULL wanttitle=n wantaxis=n wantlabel=n


/* this is a kludge again, since cake can digest too complicated rules */

#define WINDOWSIZE f1=316 n1=60 axis=123
FIGDIR/windowxx.v window0&: 
	cake FIGDIR/novar.v FIGDIR/xquizsyntB.v
	Window WINDOWSIZE  <jnormxx.K.H >jnovarxx.H
	Window WINDOWSIZE  <xq.normXx.eq.B.K.H >jeqxx.H
	Window WINDOWSIZE  <normXx.K.H >jxx.H
	Taplot < jnovarxx.H gainpanel=every pclip=98 |Ta2vplot title="Xx Blow Up With No Variations" wanttitle=n wantlabel=n wantaxis=n out=junk1.v NULL xsize=6 ysize=1
	Taplot < jeqxx.H gainpanel=every pclip=98 |Ta2vplot title="Xx Blow Up of Equalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk2.v NULL xsize=6 ysize=1
	Taplot < jxx.H gainpanel=every pclip=98 |Ta2vplot title="Xx Blow Up of Unequalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk3.v NULL xsize=6 ysize=1
	vp_OverUnderIso junk1.v junk2.v junk3.v >junk4.v
	vppen junk4.v xsize=6 ysize=3 >FIGDIR/windowxx.v


FIGDIR/windowxz.v :
	Window WINDOWSIZE <jnormxz.K.H >jnovarxz.H
	Window WINDOWSIZE <xq.normXz.eq.B.K.H >jeqxz.H
	Window WINDOWSIZE <normXz.K.H >jxz.H
	Taplot < jnovarxz.H gainpanel=every pclip=98 |Ta2vplot title="Xz Blow Up With No Variations" wanttitle=n wantlabel=n wantaxis=n out=junk1.v NULL xsize=6 ysize=1
	Taplot < jeqxz.H gainpanel=every pclip=98 |Ta2vplot title="Xz Blow Up of Equalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk2.v NULL xsize=6 ysize=1
	Taplot < jxz.H gainpanel=every pclip=98 |Ta2vplot title="Xz Blow Up of Unequalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk3.v NULL xsize=6 ysize=1
	vp_OverUnderIso junk1.v junk2.v junk3.v >junk4.v
	vppen junk4.v xsize=6 ysize=3 >FIGDIR/windowxz.v

FIGDIR/windowzx.v :
	Window WINDOWSIZE  <jnormzx.K.H >jnovarzx.H
	Window WINDOWSIZE  <xq.normZx.eq.B.K.H >jeqzx.H
	Window WINDOWSIZE  <normZx.K.H >jzx.H
	Taplot < jnovarzx.H gainpanel=every pclip=98 |Ta2vplot title="Zx Blow Up With No Variations" wanttitle=n wantlabel=n wantaxis=n out=junk1.v NULL xsize=6 ysize=1
	Taplot < jeqzx.H gainpanel=every pclip=98 |Ta2vplot title="Zx Blow Up of Equalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk2.v NULL xsize=6 ysize=1
	Taplot < jzx.H gainpanel=every pclip=98 |Ta2vplot title="Zx Blow Up of Unequalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk3.v NULL xsize=6 ysize=1
	vp_OverUnderIso junk1.v junk2.v junk3.v >junk4.v
	vppen junk4.v xsize=6 ysize=3 >FIGDIR/windowzx.v

FIGDIR/windowzz.v :
	Window f1=316 n1=60 axis=123  <jnormzz.K.H >jnovarzz.H
	Window f1=316 n1=60 axis=123  <xq.normZz.eq.B.K.H >jeqzz.H
	Window f1=316 n1=60 axis=123  <normZz.K.H >jzz.H
	Taplot < jnovarzz.H gainpanel=every pclip=98 |Ta2vplot title="Zz Blow Up With No Variations" wanttitle=n wantlabel=n wantaxis=n out=junk1.v NULL xsize=6 ysize=1
	Taplot < jeqzz.H gainpanel=every pclip=98 |Ta2vplot title="Zz Blow Up of Equalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk2.v NULL xsize=6 ysize=1
	Taplot < jzz.H gainpanel=every pclip=98 |Ta2vplot title="Zz Blow Up of Unequalized Region" wanttitle=n wantlabel=n wantaxis=n out=junk3.v NULL xsize=6 ysize=1
	vp_OverUnderIso junk1.v junk2.v junk3.v >junk4.v
	vppen junk4.v xsize=6 ysize=3 >FIGDIR/windowzz.v




/* 3 sources */
$elastic_vertical%&:  /* tismall.H wavelet.vert.H */ BINDIR/Ultiall
$        time Ultiall  type=4 ncomp=3 modeling=1 src_type=0 nsrc=3 src_1=52 src_inc=52 moduli=tismall.H par=parfile  xsnap=xsnap% zsnap=zsnap% ysnap=ysnap% xseis=xseis% yseis=yseis% zseis=zseis% <wavelet.vert.H >junk


wavelet.H:
        Wavelet n1=256 d1=0.001 domain=time wavelet=ricker1 tdelay=0.060 fund=15 fhigh=22. >junk.H
	 Scale dscale=1000000. <junk.H >wavelet.H
$        Wavelet n1=256 d1=0.001 domain=time wavelet=ricker1 tdelay=0.040 fund=25 fhigh=40. >junk.H


/* vertical source   2components */
wavelet.Z.H: wavelet.H if not exist wavelet.Z.H
	Spike n1=256 n2=2 n3=1 nsp=0 >junk
	Merge wavelet.H junk axis=2 space=n >wavelet.Z.H
	echo "d1=0.0005" >>wavelet.Z.H

wavelet.X.H: wavelet.H if not exist wavelet.X.H
	Spike n1=256 n2=1 n3=1 nsp=0 >junk1
	Spike n1=256 n2=1 n3=1 nsp=0 >junk3
	Merge junk1 wavelet.H junk3 axis=2 space=n >wavelet.X.H
	echo "d1=0.0005" >>wavelet.X.H

wavelet.Y.H: wavelet.H if not exist wavelet.Y.H
	Spike n1=256 n2=1 n3=1 nsp=0 >junk1
	Spike n1=256 n2=1 n3=1 nsp=0 >junk3
	Merge junk1 junk3 wavelet.H  axis=2 space=n >wavelet.Y.H
	echo "d1=0.0005" >>wavelet.Y.H


/* horizontal source   2components */
wavelet.hori.H: wavelet.H
	Spike n1=256 n2=2 n3=1 nsp=0 >junk
	Merge junk wavelet.H axis=2 space=n >wavelet.hori.H


tismall.H: tifull.H  if not exist tismall.H
	Window  j1=3 f2=1 j2=3          <tifull.H >junk.H
	Window  n1=208  f2=206 n2=464   <junk.H   >tismall.H

tifull.H: vpvs.H rho.H BINDIR/Makeaniso
	Merge vpvs.H rho.H axis=3 space=n >junk.H
	echo "d1=0.004 d2=0.004" >> junk.H
	BINDIR/Makeaniso < junk.H >tifull.H

timiddle.H: tifull.H
	 Window f2=500 n2=1500 n1=750 <tifull.H >timiddle.H



vpvs.H: /q2/marmousi/marmvel.H BINDIR/Makeelastic
	Scale dscale=0.001 </q2/marmousi/marmvel.H >jj.H
	BINDIR/Makeelastic <jj.H >vpvs.H

rho.H: /q2/marmousi/marmdens.H 
        Cp </q2/marmousi/marmdens.H |Scale dscale=.25 >jj.H   /* because it is a convex number */
	Scale dscale=0.001 <jj.H >rho.H

/* copy it so that it will reside on the scalable disk array in case of CM5 */
marmvel.H: if not exist marmvel.H
	Cp <../../Data/Marm/marmvel.HH >marmvel.H

marmdens.H: if not exist marmdens.H
	Cp <../../Data/Marm/marmdens.HH >marmdens.H

acoustic.H: marmvel.H marmdens.H
	Cp <marmvel.H >jj
	Window <jj j1=3 f2=1 j2=3 >junk
	Cp <marmdens.H  >jj 
	Window <jj j1=3 f2=1 j2=3 >jrho
	Cat junk >junk2
	Add mode=product junk junk2 >junk3
	Add mode=product junk3 jrho >junk4
	Cat junk4 jrho >junk6
	Window n1=208  f2=206 n2=464 <junk6 >acoustic.H
$
$	Scale dscale=1e-9 <junk4 >junk5
$	Scale dscale=0.001 <jrho >jrho2


clean&: jclean
	-RM_CMD *.x *.o  a.out core *.fcm *.f

veryclean&:
$	-RM_CMD *.x *.o  a.out core *.H 
	-rm *.x *.o  a.out core *.H Xxseis.* Xyseis.* Xzseis.* Zzseis.* Zyseis.* Zxseis.* *.v *.v3 *snap *snap.*


#ifndef SPPOPTS
#define SPPOPTS
#endif
%.fcm:  SRCDIR/%.rcm            if exist SRCDIR/%.rcm
        SPP SPPOPTS < SRCDIR/%.rcm | RATFOR -6"&" -C | FPP -DSOURCE="'"FULLSRC/%.rcm"'" -Dsource="'"FULLSRC/%.rcm"'"   > %.fcm

%.f: SRCDIR/%.rs		if exist SRCDIR/%.rs
	SPP SPPOPTS <SRCDIR/%.rs >%.spp
	SAW %.spp | RATFOR | FPP -Dsource="'"FULLSRC/%.rs"'" MACH ICAT(SRCDIR)  FIXNULL > %.f
	-rm %.spp

BINDIR/%: if not {{ cd ../../Bin ; cake -q %}}
	( cd ../../Bin ; cake %)

#include <SEP.obj.rules>
#include <SEP.prog.rules>
#include <SEP.idoc.rules>
