#define BINDIR ../../Bin/MTYPE
$#define FIGLIST elvisduff diffvisduff3 diffvisduff2 
#define FIGLIST elvis  diffvisduff3 
$#define FIGDIR ./Fig
#define FIGDIR ../Fig
#include <SEP.defs>
#define ULTIMOD Ultiall
#define NULL head=/dev/null

#define CUBEFLAT flat=y
#define CUBEFLAT flat=n
#define CUBEFRAMES frame3=15 frame2=140
#define CUBEFRAMES frame3=150 frame2=140
#define CUBELABELS label1="Depth [km]" label2="Distance [km]" label3="Time [sec]"
#define THLABELS label3="Amplitude" label1="Depth [km]" label2="Distance [km]"

FIGDIR/diffvisduff2.v:
	Window < zsnap.duff-0.001-0.00001-0.00001.H f3=150 n3=1 j2=4 j1=4 >jjd
	Window < zsnap.vis-0.001.H f3=150 n3=1 j2=4 j1=4 >jjv
	Add scale=1,-1 jjv jjd >jjsub
	Clip <jjv >jjvc
	Thplot <jjd clip=`Get clip <jjvc parform=n` title="High-Order" THLABELS out=jjd.v NULL
	Thplot <jjsub clip=`Get clip <jjvc parform=n` title="Difference from Visco-Elastic" THLABELS out=jjsub.v NULL
	vp_SideBySideAniso jjd.v jjsub.v >FIGDIR/diffvisduff2.v
	
FIGDIR/diffvisduff3.v: zsnap.diffduff.v zsnap.vis-0.001.H  zsnap.duff-0.001-0.00001-0.00001.H
	Window < zsnap.vis-0.001.H f3=150 n3=1 >jjv
	Clip <jjv >jjvc
	Add scale=1,-1 zsnap.vis-0.001.H zsnap.duff-0.001-0.00001-0.00001.H >jjsub
	Taplot < zsnap.vis-0.001.H clip=`Get clip <jjvc parform=n` > zsnap.vis-0.001.T
	Taplot < jjsub  clip=`Get clip <jjvc parform=n` > jjsub.T
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < jjsub.T out=zsnap.dvd.v NULL title="Higher-Order - Visco Difference"
$
$	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.vis-0.001.T out=zsnap.vis.0.001.v NULL title="Visco Elastic"
$	vppen scale=0.5 align=lb zsnap.vis-0.001.v >j1.v
$
	vppen scale=0.5 align=lb zsnap.diffduff.v  >j1.v
	vppen scale=0.5 align=lb zsnap.dvd.v       >j2.v
$	vp_OverUnderIso j1.v j2.v >FIGDIR/diffvisduff3.v
	vp_OverUnderIso j2.v j1.v >FIGDIR/diffvisduff3.v
	vp_Movie zsnap.vis.v zsnap.dvd.v zsnap.diffduff.v >FIGDIR/diffvisduff3.v3
$	vp_Movie zsnap.vis-0.001.v  zsnap.dvd.v >FIGDIR/diffvisduff3.v3

zsnap.diffduff.v: zsnap.duff-0.001-0.00001-0.H zsnap.duff-0.001-0-0.00001.H
	Add scale=1,-1 zsnap.duff-0.001-0.00001-0.H zsnap.duff-0.001-0-0.00001.H >jjsub.H
	Taplot gainpanel=150 <jjsub.H >jjsub.T
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < jjsub.T out=zsnap.diffduff.v NULL title="Higher Order Term Difference"

/* subtract the snap shot cubes */



FIGDIR/elvis.v elvis.v: zsnap.vis.v zsnap.diff.el.vis.v
	vppen scale=0.5 align=lb zsnap.vis.v >j1.v
	vppen scale=0.5 align=lb zsnap.diff.el.vis.v       >j2.v
	vp_OverUnderIso j1.v j2.v >FIGDIR/elvis.v
	vp_Movie zsnap.vis.v  zsnap.diff.el.vis.v >FIGDIR/elvis.v3

diffvisduff.v: zsnap.diff.duff.v zsnap.diff.vis.duff.v
	vp_OverUnderAniso zsnap.diff.duff.v zsnap.diff.vis.duff.v >diffvisduff.v

elvisduff.v: zsnap.el.v zsnap.diff.el.vis.v zsnap.diff.duff.v
	vp_OverUnderAniso zsnap.el.v zsnap.diff.el.vis.v zsnap.diff.duff.v >elvisduff.v

/* the individual Cube plots */
zsnap.el.v:  zsnap.el.S.T    /* just elastic */
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.el.S.T out=zsnap.el.v NULL        title="Elastic"

zsnap.diff.el.vis.v:  zsnap.diff.el.vis.T  /* el - visco */
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.diff.el.vis.T out=zsnap.diff.el.vis.v NULL title="Visco - Elastic Difference"

zsnap.diff.duff.v: zsnap.diff.duff.T   /* el - duff */
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.diff.duff.T out=zsnap.diff.duff.v NULL title="Higher-Order - Elastic Difference"

zsnap.diff.vis.duff.v: zsnap.diff.vis.duff.T   /* visco - duff */
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.diff.vis.duff.T out=zsnap.diff.vis.duff.v NULL title="Higher-Order - Visco Difference"

zsnap.vis.v: zsnap.vis.S.T    /* visco */
	Cubeplot CUBEFRAMES CUBELABELS CUBEFLAT < zsnap.vis.S.T out=zsnap.vis.v NULL title="Visco Elastic"


zsnap.diff.vis.duff.T:  zsnap.vis.S.T zsnap.diff.vis.duff.H
	Taplot gainpanel=150 clip=`Get < zsnap.vis.S.T clip parform=n` < zsnap.diff.vis.duff.H > zsnap.diff.vis.duff.T

zsnap.diff.el.vis.T:  zsnap.vis.S.T zsnap.diff.el.vis.H
	Taplot gainpanel=150 clip=`Get < zsnap.vis.S.T clip parform=n` < zsnap.diff.el.vis.H > zsnap.diff.el.vis.T


%.T: %.H
	Taplot gainpanel=150 <%.H >%.T

/* take the z component cube of  el vis and duff subtract */
zsnap.diff.el.vis.H: zsnap.el.S.H zsnap.vis.S.H
	Add scale=1,-1  zsnap.el.S.H zsnap.vis.S.H > zsnap.diff.el.vis.H

zsnap.diff.duff.H: zsnap.el.S.H zsnap.duff.S.H
	Add scale=1,-1  zsnap.el.S.H zsnap.duff.S.H >zsnap.diff.duff.H

zsnap.diff.vis.duff.H: zsnap.vis.S.H zsnap.duff.S.H
	Add scale=1,-1  zsnap.el.S.H zsnap.duff.S.H >zsnap.diff.vis.duff.H

/* scale them to a common maximum to conentrate on differences in wavelet and
   not bias in the overall amplitude  */

zsnap.el.S.H: zsnap.el.H
	Scale axis=123 <zsnap.el.H >zsnap.el.S.H

zsnap.duff.S.H: zsnap.duff-0.001-0.00001-0.00001.H
	Scale axis=123 <zsnap.duff-0.001-0.00001-0.00001.H >zsnap.duff.S.H

zsnap.vis.S.H:  zsnap.vis-0.001.H 
	Scale axis=123 <zsnap.vis-0.001.H >zsnap.vis.S.H



run&: el duff-0.001-0.00001-0.00001 vis-0.001
	cake only

only&: duff-0.001-0.00001-0  duff-0.001-0-0.00001


/* just elastic */
el& zsnap.el.H: wavelet.vertfine.H homojustiso.H
	time BINDIR/ULTIMOD type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.el.H ysnap=ysnap.el.H zsnap=zsnap.el.H xseis=xseis.el.H yseis=yseis.el.H zseis=zseis.el.H <wavelet.vertfine.H >junk nt=2000

/* duffing equation */
duff-%1-%2-%3& zsnap.duff-%1-%2-%3.H: homojustiso.H wavelet.vertfine.H
	 Scale dscale=%1 <homojustiso.H >homojustiso.visco.H
	time BINDIR/ULTIMOD visco=1 viscomoduli=homojustiso.visco.H duffing=1 hold=%2 spring=%3 type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.duff-%1-%2-%3.H ysnap=ysnap.duff-%1-%2-%3.H zsnap=zsnap.duff-%1-%2-%3.H xseis=xseis.duff-%1-%2-%3.H yseis=yseis.duff-%1-%2-%3.H zseis=zseis.duff-%1-%2-%3.H <wavelet.vertfine.H >junk nt=2000

/* just viscous */
vis-%1& zsnap.vis-%1.H: homojustiso.H wavelet.vertfine.H
	 Scale dscale=%1 <homojustiso.H >homojustiso.visco.H
	time BINDIR/ULTIMOD visco=1 viscomoduli=homojustiso.visco.H type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.vis-%1.H ysnap=ysnap.vis-%1.H zsnap=zsnap.vis-%1.H xseis=xseis.vis-%1.H yseis=yseis.vis-%1.H zseis=zseis.vis-%1.H <wavelet.vertfine.H >junk nt=2000

compare:
	Scale <xseis.VISCO.H axis=123  >j1
	Scale <xseis.ISO.duff.H axis=123 >j2
	Merge j1 j2 axis=3 space=n >j3
	Window f1=80 n1=1 <j3 |Graph |Tube

compare:
	Scale <zseis.VISCO.H axis=123  >j1
	Scale <zseis.ISO.duff.H axis=123 >j2

/* 
visco 0.005 and hold 0.1 and spring 0.1   fails nt=200
visco 0.005 and hold 0.01 and spring 0.01   fails
visco 0.005 and hold 0.001 and spring 0.001   fails

visco 0.005 and hold 0.0001 and spring 0.0001 works on nt=200, fails on nt=2000

visco 0.001 and hold 0.0001 and spring 0.0001   works on nt=2000
visco 0.0005 and hold 0.0001 and spring 0.0001   works on nt=2000

*/

testdufffine&: homojustiso.H wavelet.vertfine.H
	 Scale dscale=0.001 <homojustiso.H >homojustiso.visco.H
        time BINDIR/ULTIMOD  visco=1 viscomoduli=homojustiso.visco.H duffing=1 hold=0.00001 spring=0.00001 type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.ISO.duff.H ysnap=ysnap.ISO.duff.H zsnap=zsnap.ISO.duff.H xseis=xseis.ISO.duff.H yseis=yseis.ISO.duff.H zseis=zseis.ISO.duff.H <wavelet.vertfine.H >junk nt=2000
        time BINDIR/ULTIMOD visco=1 viscomoduli=homojustiso.visco.H type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.VISCO.H ysnap=ysnap.VISCO.H zsnap=zsnap.VISCO.H xseis=xseis.VISCO.H yseis=yseis.VISCO.H zseis=zseis.VISCO.H <wavelet.vertfine.H >junk nt=2000


testduff&: homojustiso.H wavelet.vert.H
        time BINDIR/ULTIMOD duffing=1 hold=0 spring=0 type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.ISO.duff.H ysnap=ysnap.ISO.duff.H zsnap=zsnap.ISO.duff.H xseis=xseis.ISO.duff.H yseis=yseis.ISO.duff.H zseis=zseis.ISO.duff.H <wavelet.vert.H >junk

testnorm&: homojustiso.H wavelet.vert.H
        time BINDIR/ULTIMOD type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.NORM.H ysnap=ysnap.NORM.H zsnap=zsnap.NORM.H xseis=xseis.NORM.H yseis=yseis.NORM.H zseis=zseis.NORM.H <wavelet.vert.H >junk

testnormfine&: homojustiso.H wavelet.vertfine.H
        time BINDIR/ULTIMOD type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.NORM.H ysnap=ysnap.NORM.H zsnap=zsnap.NORM.H xseis=xseis.NORM.H yseis=yseis.NORM.H zseis=zseis.NORM.H <wavelet.vertfine.H >junk  nt=2000

/* at 2ms
    ok: 0.00005
    	0.0001      in peaks visible

not ok:
	+-0.001       blows up

with fine: 0.1 ms
	0.001	works nicely
	-0.001 blows up as expected
	0.005   works well,  scale the NORM and VISCO data 1 each and substract
*/

testvisco&: homojustiso.H wavelet.vertfine.H 
	Scale dscale=0.005 <homojustiso.H >homojustiso.visco.H
        time BINDIR/ULTIMOD visco=1 viscomoduli=homojustiso.visco.H type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.VISCO.H ysnap=ysnap.VISCO.H zsnap=zsnap.VISCO.H xseis=xseis.VISCO.H yseis=yseis.VISCO.H zseis=zseis.VISCO.H <wavelet.vertfine.H >junk nt=2000

testfin&: homojustiso.H wavelet.vert.H
        time BINDIR/ULTIMOD finite=1 type=2 src_type=0 ncomp=2 modeling=1 moduli=homojustiso.H par=parfile  xsnap=xsnap.ISO.fin ysnap=ysnap.ISO.fin zsnap=zsnap.ISO.fin xseis=xseis.ISO.fin yseis=yseis.ISO.fin zseis=zseis.ISO.fin <wavelet.vert.H >junk


BINDIR/Sinusoid: Sinusoid.o
	FLINK(Sinusoid.o,BINDIR/Sinusoid)

wavelet.vertfine.H:
	Wavelet n1=256 d1=0.0005 domain=time wavelet=ricker1 tdelay=0.060 fund=15 fhigh=30. | Scale dscale=10000. >junk
	Spike n1=256 n2=1 n3=1 d1=0.0005 k1=-1 k2=-1 k3=-1 >junk2
	Merge junk junk2 space=n axis=2 >wavelet.vertfine.H

wavelet.vert.H:
	Wavelet n1=128 d1=0.002 domain=time wavelet=ricker1 tdelay=0.060 fund=15 fhigh=30. | Scale dscale=10000. >junk
	Spike n1=128 n2=1 n3=1 d1=0.002 k1=-1 k2=-1 k3=-1 >junk2
	Merge junk junk2 space=n axis=2 >wavelet.vert.H

/* make a two component source */
homoiso.H: BINDIR/Geninput
	BINDIR/Geninput nx=208 nz=208 dx=20. dz=20.  symmetry=2 vp=4.0 vs=2.3  rho=2.5   >homoiso.H

homojustiso.H: homoiso.H
	Window n3=1 <homoiso.H >junk.c11
	Window f3=3 n3=1  <homoiso.H >junk.c55
	Window f3=4 <homoiso.H >junk.rho
	Merge junk.c11 junk.c55 junk.rho axis=3 space=n >homojustiso.H

homoac.h:
	Window n3=1 <homoiso.H >junk.c11
	Window f3=4 <homoiso.H >junk.rho
	Merge junk.c11 junk.rho axis=3 space=n >homoac.h

homoiso3c.h: homoiso.H
	Window n3=4 <homoiso.H >junk      /* get sitffnesses only */
	Window f3=4 <homoiso.H >junk2     /* get rho */
	Window f3=3 n3=1 <junk >junk3    /* get c55*/
	Merge junk junk3 junk2 axis=3 space=n >homoiso3c.h
	RM_CMD junk*
	
isolayer.h:
	BINDIR/Geninput nx=208 nz=70 dx=20. dz=20.  symmetry=2 vp=4.0 vs=2.3  rho=2.5   >junk0
	BINDIR/Geninput nx=208 nz=138 dx=20. dz=20.  symmetry=2 vp=4.5 vs=3.2  rho=2.0   >junk1
	Merge junk0 junk1  space=n axis=2 |Transp plane=12 >isolayer.h


clean&: jclean
	-RM_CMD *.x *.o  a.out core junk*  j* *.T

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