#define FIGDIR ./Fig
#define ULTIMOD Ultiall
#define FIGLIST equivalns mlong+fshort
#define NEEDF90 FIGLIST
#include "../Cakerules/SEP.martin.rules"
#define BINDIR ../Bin/MTYPE
#include <SEP.defs>
#define WAVELETARGS n1=128 d1=0.001  tdelay=0.060
#define INPUTDIM nx=208 nz=104 dx=10. dz=10
$#define NONSTAGGERED  ncoef=16  coeff=long_nonstag.H
$#define NONSTAGSHORT  ncoef2=10 coeff2=short_nonstag.H

#define NONSTAGGERED  ncoef=1  coeff=coeffc.H
#define NONSTAGSHORT  ncoef2=1 coeff2=coeffd.H


#defien FATMULT 2
#define GRAPHLABELS label1="Time [sec]" label2="Amplitude" plotfat=6,6,6

equivalns.slidefigures.resize&:  backcol.v
        vppen backcol.v FIGDIR/equivalns.v rotate=90 |vppen align=lb| SLIDETUBE scale=0.75 xshift=3 yshift=2

mlong+fshort.slidefigures.resize&:  backcol.v
        vppen backcol.v FIGDIR/mlong+fshort.v rotate=90 |vppen align=lb| SLIDETUBE scale=0.75 xshift=3 yshift=2

	
	
/* is very close to the composite long, but with an operator reduction of 1:2*/
/* at cost of a little memory */
mlong+fshort.H: layerall.H wavelet.vert.H coeffc.H coeffd.H BINDIR/ULTIMOD if not exist mlong+fshort.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk    split=1 zseis=zseisd.split.H zsnap=zsnapd.split.H xseis=xseisd.split.H xsnap=xsnapd.split.H NOstaggrid=0 ncoef=1  coeff=coeffd.H ncoef2=1 coeff2=coeffc.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk   split=0 NOstaggrid=0 ncoef=1  coeff=coeffd.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk   split=0 zseis=zseisc.split.H zsnap=zsnapc.split.H xseis=xseisc.split.H xsnap=xsnapc.split.H NOstaggrid=0 ncoef=1  coeff=coeffc.H
	Merge xseis xseisd.split.H xseisc.split.H axis=1 space=y >j1
	Merge zseis zseisd.split.H zseisc.split.H axis=1 space=y >j2
	Merge j1 j2 axis=2 space=y >mlong+fshort.H

FIGDIR/mlong+fshort.v: mlong+fshort.H
	Merge xseis xseisd.split.H xseisc.split.H axis=3 space=n >j1
	Merge zseis zseisd.split.H zseisc.split.H axis=3 space=n >j2
	Window n1=1 f1=90 <j1 |Graph title="x long+short" GRAPHLABELS out=jx.v >/dev/null
	Window n1=1 f1=90 <j2 |Graph title="z long+short" GRAPHLABELS out=jz.v >/dev/null
	vp_OverUnderIso jx.v jz.v >FIGDIR/mlong+fshort.v





long+short.H: layerall.H wavelet.vert.H coeffc.H coeffd.H BINDIR/ULTIMOD if not exist long+short.H
$	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk    split=1 zseis=zseisd.split.H zsnap=zsnapd.split xseis=xseisd.split.H xsnap=xsnapd.split.H NOstaggrid=0 NONSTAGGERED NONSTAGSHORT
$	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk   split=0 NOstaggrid=0 NONSTAGGERED
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk    split=1 zseis=zseisd.split.H zsnap=zsnapd.split.H xseis=xseisd.split.H xsnap=xsnapd.split.H NOstaggrid=0 ncoef=1  coeff=coeffd.H ncoef2=1 coeff2=coeffc.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk   split=0 NOstaggrid=0 ncoef=1  coeff=coeffd.H
	Merge xseis xseisd.split.H axis=1 space=y >j1
	Merge zseis zseisd.split.H axis=1 space=y >j2
	Merge j1 j2 axis=2 space=y >long+short.H

FIGDIR/long+short.v: long+short.H
	Merge xseis xseisd.split.H axis=3 space=n >j1
	Merge zseis zseisd.split.H axis=3 space=n >j2
	Window n1=1 f1=90 <j1 |Graph title="x long+short" GRAPHLABELS out=jx.v >/dev/null
	Window n1=1 f1=90 <j2 |Graph title="z long+short" GRAPHLABELS out=jz.v >/dev/null
	vp_OverUnderIso jx.v jz.v >FIGDIR/long+short.v


#define NONSTAGGERED  ncoef=16  coeff=long_nonstag.H
#define NONSTAGSHORT  ncoef2=10 coeff2=short_nonstag.H

/* equivalence for the non staggered grid */
equivalns.H:   layerall.H wavelet.vert.H BINDIR/ULTIMOD long_nonstag.H   if not exist equivalns.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk    split=1 zseis=zseisd.split.H zsnap=zsnapd.split.H xseis=xseisd.split.H xsnap=xsnapd.split.H staggrid=0 ncoef=16  coeff=long_nonstag.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layerall.H par=parfile.transmission <wavelet.vert.H >junk   split=0 staggrid=0 ncoef=16  coeff=long_nonstag.H
	Merge xseis xseisd.split.H axis=1 space=y >j1
	Merge zseis zseisd.split.H axis=1 space=y >j2
	Merge j1 j2 axis=2 space=y >equivalns.H

FIGDIR/equivalns.v: equivalns.H
	Merge xseis xseisd.split.H axis=3 space=n >j1
	Merge zseis zseisd.split.H axis=3 space=n >j2
	Window n1=1 f1=90 <j1 |Graph title="x component" GRAPHLABELS out=jx.v >/dev/null
	Window n1=1 f1=90 <j2 |Graph title="z component" GRAPHLABELS out=jz.v >/dev/null
	vp_OverUnderIso jx.v jz.v >FIGDIR/equivalns.v


FIGDIR/equiwava.v: equivala.H
$	Merge xsnap xsnapd.split.H axis=3 space=n >j1
$	Merge zsnap zsnapd.split.H axis=3 space=n >j2
$	Window n3=1 f3=4 <j1 >jj1
$	Byte <jj1 |Ta2vplot title="x component snapshot" out=jx.v >/dev/null
$	Window n3=1 f3=4 <j2 >jj2
$	Byte <jj2 |Ta2vplot  title="z component snapshot" out=jz.v >/dev/null
$	vp_OverUnderIso jx.v jz.v >FIGDIR/equiwava.v
	Window <xsnap n3=1 f3=6 >jj1
	 Byte <jj1 |Ta2vplot title="x component snapshot" out=jx.v >/dev/null
	cp jx.v FIGDIR/equiwava.v


/* layered model regular , source on top, receivers on bottom  */

transmission.H&: /* layer.h wavelet.vert.H  */
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layer.h par=parfile.transmission <wavelet.vert.H >junk    split=1 zseis=zseis.split.H zsnap=zsnap.split.H xseis=xseis.split.H xsnap=xsnap.split.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layer.h par=parfile.transmission <wavelet.vert.H >junk   split=0 



/* layered model, source on top, receivers on top  */
reflection.H&: /* homoiso.h wavelet.vert.H  */
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layer.h par=parfile.reflection <wavelet.vert.H >junk   split=1 zseis=zseis.split.H zsnap=zsnap.split.H xseis=xseis.split.H xsnap=xsnap.split.H
	time BINDIR/ULTIMOD  type=4 modeling=1 moduli=layer.h par=parfile.reflection <wavelet.vert.H >junk   split=0 


short_nonstag.H:
        echo "-0.2713077307504363E-02 0.9277622829117149E-02 -0.3045387209576234E-01 0.1178572129899180E+00 -0.1247830531633298E+01 0.1247830531633298E+01 -0.1178572129899180E+00  0.3045387209576234E-01 -0.9277622829117149E-02  0.2713077307504363E-02 " >jascii
	atoF <jascii >short_nonstag.H

long_nonstag.H:
$        echo "0.1326224603581672E-04 -0.1260422755482438E-03 0.6766549997685445E-03 -0.2713077307504363E-02 0.9277622829117149E-02 -0.3045387209576234E-01 0.1178572129899180E+00 -0.1247830531633298E+01 0.1247830531633298E+01 -0.1178572129899180E+00  0.3045387209576234E-01 -0.9277622829117149E-02  0.2713077307504363E-02 -0.6766549997685445E-03  0.1260422755482438E-03 -0.1326224603581672E-04" >jascii
        echo "0.1326224603581672E-04 -0.1260422755482438E-03 0.6766549997685445E-03 -0.2713077307504363E-02 0.9277622829117149E-02 -0.3045387209576234E-01 0.1178572129899180E+00 -0.1247830531633298E+01 0 0.1247830531633298E+01 -0.1178572129899180E+00  0.3045387209576234E-01 -0.9277622829117149E-02  0.2713077307504363E-02 -0.6766549997685445E-03  0.1260422755482438E-03 -0.1326224603581672E-04" >jascii
	atoF <jascii >long_nonstag.H

coeffc.H:
        echo "0.1326224603581672E-04 -0.1260422755482438E-03 0.6766549997685445E-03 -0.2713077307504363E-02 0.9277622829117149E-02 -0.3045387209576234E-01 0.1178572129899180E+00 -0.1247830531633298E+01 0.1247830531633298E+01 -0.1178572129899180E+00  0.3045387209576234E-01 -0.9277622829117149E-02  0.2713077307504363E-02 -0.6766549997685445E-03  0.1260422755482438E-03 -0.1326224603581672E-04" >jascii
	atoF <jascii >coeffc.H


coeffd.H :
$         echo "-0.0095703125 0.079752604 -1.1962891 1.1962891 -0.079752604 0.0095703125" >jascii
         echo "0.079752604 -1.1962891 1.1962891 -0.079752604" >jascii
	atoF <jascii >coeffd.H



/* vertical source */
wavelet.vert.H:
	Wavelet WAVELETARGS  domain=time wavelet=ricker2 fund=15 fhigh=30. | Scale dscale=1000000. >junk
	Spike WAVELETARGS n2=1 n3=1 k1=-1 k2=-1 k3=-1 >junk2
	Merge junk junk2 space=n axis=2 >wavelet.vert.H

/* vertical source 1ms */
wavelet.fine.vert.H:
	Wavelet n1=128 d1=0.001 domain=time wavelet=ricker2 tdelay=0.060 fund=15 fhigh=30. | Scale dscale=1000000. >junk
	Spike n1=128 n2=1 n3=1 d1=0.001 k1=-1 k2=-1 k3=-1 >junk2
	Merge junk junk2 space=n axis=2 >wavelet.fine.vert.H


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

layerall.H:  BINDIR/Geninput
	BINDIR/Geninput INPUTDIM  symmetry=2 vp=4.5 vs=2.5  rho=2.5   >j1
	BINDIR/Geninput INPUTDIM  symmetry=2 vp=3.5 vs=1.9  rho=2.5   >j2
	Merge j1 j2 space=n axis=1 >layerall.H


/* isotropic  input   */
homoiso.h: BINDIR/Geninput
	BINDIR/Geninput nx=208 nz=208 dx=40. dz=40.  symmetry=2 vp=4.0 vs=2.3  rho=2.5   >homoiso.h

homoiso.fine.h: BINDIR/Geninput
	BINDIR/Geninput nx=208 nz=208 dx=10. dz=10.  symmetry=2 vp=4.0 vs=2.3  rho=2.5   >homoiso.fine.h

homojustiso.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*



haskell& junk.solid2:   Mult.a
        Haskellpsv nt=512 dt=.004 nx=256 dx=20. layers=Mult.a receiver=uz lsource=2 lreceiver=4 >junk.solid2


data2.H: junk.solid2
        Window <wavelet.vert.H n2=1 >junk.wav
        Scale dscale=1.e14 <junk.solid2 | Conv filtin=junk.wav >temp.H
        Mute vmute=6000 <temp.H >junk.mute
        Bandpass flo=2. fhi=100 <junk.mute  >data2.H

sem&:
	tube yscale=0.8 coeffa.slide coeffb.slide coeffc.slide coeffd.slide coeffn.slide
	tube yscale=0.8 FIGDIR/equivala.v FIGDIR/long+short.v

tcpr&:
	vppen rotate=90 scale=2. <FIGDIR/equivala.v |vppen align=lb|tcprpen
	vppen rotate=90 scale=2. <FIGDIR/long+short.v |vppen align=lb|tcprpen

slides&:  coeffb.H coeffc.H coeffd.H  coeffe.H 
	echo "0,0,0,-1,1,0,0,0" > jco.as
	atoF <jco.as >jco.H
	Ftplot < jco.H    title="-1,1 central difference" out=coeffa.slide 
	Ftplot < coeffb.H title="Francis central limit" out=coeffb.slide
	Ftplot < coeffc.H title="Holberg optimized 16" out=coeffc.slide
        Ftplot < coeffd.H title="Holberg optimized 6"  out=coeffd.slide
        Ftplot < coeffe.H title="Holberg optimized 4"  out=coeffe.slide
	echo "0,0,0,-0.5,0,0.5,0,0" > jco.as
	atoF <jco.as >jco_nonstag.H
	Ftplot < jco_nonstag.H title="-1/2,0,1/2" out=coeffn.slide


BINDIR/Geninput: if not {{ cd ../Bin ; cake -q MTYPE/Geninput }}
	( cd ../Bin ; cake MTYPE/Geninput )
	
BINDIR/ULTIMOD: if not {{ cd ../Bin ; cake -q MTYPE/ULTIMOD }}
	( cd ../Bin ; cake MTYPE/ULTIMOD )
	
clean&: jclean
	-RM_CMD zsnap zseis  xseis xsnap 
	-RM_CMD j* 
	-/bin/rm -f *.slide

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