#define FIGDIR ../Fig
#define BINDIR ../../Bin/MTYPE
#define SEPLIB_NEW -lsepf -lsepCM5

#define FIGLIST rrot.-60.60.i drot.-60.60.i rrot.-60.60.b drot.-60.60.b medium seismo.i

#define DEBUG
#define CMLIBS
#define ULTIMOD Ultiall
#include <SEP.defs>

$#define MAXMIN min2=-2100 max2=1400 dash=0,1
#define MAXMIN dash=0,1
#define ANNOTSEIS label2="distance [km]" label1="time [sec]" NULL
#define LABELSZ titlesz=20 labelsz=12 titlefat=4 labelfat=4

#define NULL head=/dev/null
#define PSARGS color=y 


movies&: [[sub -i X FIGDIR/X drot.i.v3 rrot.i.v3 drot.b.v3 rrot.b.v3 ]]

seismo.%.action&:
	TUBE TUBEARGS scale=0.5 FIGDIR/seismo.%.v

FIGDIR/drot.0.0.%.v:  panel.0.0.% 
	vp_OverUnderAniso panel.0.0.% >FIGDIR/drot.0.0.%.v
	cake FIGDIR/drot.%.v3


FIGDIR/drot.-60.60.%.v: panel.-60.60.% 
	vp_OverUnderAniso panel.-60.60.% >FIGDIR/drot.-60.60.%.v
	cake FIGDIR/drot.%.v3

FIGDIR/rrot.0.0.%.v: rpanel.0.0.%
	vp_OverUnderAniso rpanel.0.0.% >FIGDIR/rrot.0.0.%.v
	cake FIGDIR/rrot.%.v3

FIGDIR/rrot.-60.60.%.v: rpanel.-60.60.%
	vp_OverUnderAniso rpanel.-60.60.% >FIGDIR/rrot.-60.60.%.v
	cake FIGDIR/rrot.%.v3

FIGDIR/rrot.%.v3: rpanel.0.0.%  rpanel.-30.30.% rpanel.-60.60.%
	vppen vpstyle=n rpanel.0.0.% rpanel.-30.30.% rpanel.-60.60.% >FIGDIR/rrot.%.v3

FIGDIR/drot.%.v3: panel.0.0.% panel.-30.30.% panel.-60.60.%
	vppen vpstyle=n panel.0.0.% panel.-30.30.% panel.-60.60.% >FIGDIR/drot.%.v3

/* 
n2 n3  (n4)
   Xx  Xz
   Zx  Zz  */

FIGDIR/seismo.%.v:  Hxseis1.% Vxseis1.% Hzseis1.% Vzseis1.%
	Byte <Hxseis1.% transp=yes|Ta2vplot out=Hxseis1.%.v LABELSZ titlesz=20 ANNOTSEIS title="Xx"
	Byte <Vxseis1.% transp=yes|Ta2vplot out=Vxseis1.%.v LABELSZ ANNOTSEIS title="Zx"
	Byte <Hzseis1.% transp=yes|Ta2vplot out=Hzseis1.%.v LABELSZ ANNOTSEIS title="Xz"
	Byte <Vzseis1.% transp=yes|Ta2vplot out=Vzseis1.%.v LABELSZ ANNOTSEIS title="Zz"
	vp_SideBySideAniso Hxseis1.%.v Hzseis1.%.v >j1.v
	vp_SideBySideAniso Vxseis1.%.v Vzseis1.%.v >j2.v
	vp_OverUnderAniso  j1.v j2.v |vppen txscale=1.5 >FIGDIR/seismo.%.v
	vppen Hxseis1.%.v Vxseis1.%.v Hzseis1.%.v Vzseis1.%.v >FIGDIR/seismo.%.v3
	

drotate.%1.%2.%3.H:  /* slope.%BINDIR */ BINDIR/Rot  Hxseis1.%3 Vxseis1.%3 Hzseis1.%3 Vzseis1.%3 Hxseis2.%3 Vxseis2.%3 Hzseis2.%3 Vzseis2.%3
	Merge  Hxseis1.%3 Vxseis1.%3 Hzseis1.%3 Vzseis1.%3  >jj.H space=n axis=3
	Window f1=155 n1=1 <jj.H >jj2.H
	echo "n2=1 n3=4" >> jj2.H
	BINDIR/Rot src=%1 geo=%2 < jj2.H >drotate.%1.%2.%3.1.H
	Merge  Hxseis2.%3 Vxseis2.%3 Hzseis2.%3 Vzseis2.%3 >jj.H space=n axis=3
	Window f1=51 n1=1 <jj.H >jj3.H
	echo "n2=1 n3=4" >> jj3.H
	BINDIR/Rot src=%2 geo=%1 < jj3.H >drotate.%1.%2.%3.2.H
	Merge drotate.%1.%2.%3.1.H drotate.%1.%2.%3.2.H >drotate.%1.%2.%3.H axis=2 space=n
	Rm jj2.H jj3.H jj.H

rrotate.%1.%2.%3.H:  drotate.%1.%2.%3.H
	Window <drotate.%1.%2.%3.2.H  f3=0 n3=1 >jXxseis2.H
	Window <drotate.%1.%2.%3.2.H  f3=1 n3=1 >jZxseis2.H
	Window <drotate.%1.%2.%3.2.H  f3=2 n3=1 >jXzseis2.H
	Window <drotate.%1.%2.%3.2.H  f3=3 n3=1 >jZzseis2.H
	Merge jXxseis2.H jXzseis2.H  jZxseis2.H jZzseis2.H  space=n axis=3 >jdrotate.%1.%2.%3.2.H
	Merge drotate.%1.%2.%3.1.H jdrotate.%1.%2.%3.2.H >rrotate.%1.%2.%3.H axis=2 space=n
	Rm j*seis2.H

panel.%1.%2.%3: drotate.%1.%2.%3.H
	Window f3=0 n3=1 <drotate.%1.%2.%3.H >jXx.%1.%2.H
	Window f3=1 n3=1 <drotate.%1.%2.%3.H >jZx.%1.%2.H
	Window f3=2 n3=1 <drotate.%1.%2.%3.H >jXz.%1.%2.H
	Window f3=3 n3=1 <drotate.%1.%2.%3.H >jZz.%1.%2.H
	Graph <jXx.%1.%2.H LABELSZ title="Xx  (%1 %2 degree)" out=jXx.%1.%2.v head=/dev/null MAXMIN
	Graph <jZx.%1.%2.H LABELSZ title="Zx  (%1 %2 degree)" out=jZx.%1.%2.v head=/dev/null MAXMIN
	Graph <jXz.%1.%2.H LABELSZ title="Xz  (%1 %2 degree)" out=jXz.%1.%2.v head=/dev/null MAXMIN
	Graph <jZz.%1.%2.H LABELSZ title="Zz  (%1 %2 degree)" out=jZz.%1.%2.v head=/dev/null MAXMIN
	vp_SideBySideAniso jXx.%1.%2.v jXz.%1.%2.v >jone.v 
	vp_SideBySideAniso jZx.%1.%2.v jZz.%1.%2.v >jtwo.v 
	vp_OverUnderAniso jone.v jtwo.v |vppen txscale=1.5 > panel.%1.%2.%3

rpanel.%1.%2.%3: rrotate.%1.%2.%3.H
	Window f3=0 n3=1 <rrotate.%1.%2.%3.H >jXx.%1.%2.H
	Window f3=1 n3=1 <rrotate.%1.%2.%3.H >jZx.%1.%2.H
	Window f3=2 n3=1 <rrotate.%1.%2.%3.H >jXz.%1.%2.H
	Window f3=3 n3=1 <rrotate.%1.%2.%3.H >jZz.%1.%2.H
	Graph <jXx.%1.%2.H LABELSZ title="Xx recip. (%1 %2 degree)" out=jXx.%1.%2.v head=/dev/null MAXMIN
	Graph <jZx.%1.%2.H LABELSZ title="Zx recip.  (%1 %2 degree)" out=jZx.%1.%2.v head=/dev/null MAXMIN
	Graph <jXz.%1.%2.H LABELSZ title="Xz recip.  (%1 %2 degree)" out=jXz.%1.%2.v head=/dev/null MAXMIN
	Graph <jZz.%1.%2.H LABELSZ title="Zz recip.  (%1 %2 degree)" out=jZz.%1.%2.v head=/dev/null MAXMIN
	vp_SideBySideAniso jXx.%1.%2.v jXz.%1.%2.v >jone.v 
	vp_SideBySideAniso jZx.%1.%2.v jZz.%1.%2.v >jtwo.v 
	vp_OverUnderAniso jone.v jtwo.v |vppen txscale=1.5> rpanel.%1.%2.%3



FIGDIR/medium.v: homoisoslope.H
	Window <homoisoslope.H n3=1 >jj
	Byte <jj allpos=y |Ta2vplot out=jj.v LABELSZ title="Reciprocal Setup" label2="distance [km]" label1="depth [km]" head=/dev/null
$             151*0.04   103*0.04   155*0.04 103*0.04
	echo "2.04 4.12 6.2 4.12" >jjecho.H
	atoF <jjecho.H >jj2.H
	echo esize=8 n1=1 n2=2 >>jj2.H
	Graph min1=0 min2=0 max1=8.28 max2=8.28 yreverse=1 symbolsz=10 symbol="*" wantaxis=n title=" " <jj2.H out=jj2.v head=/dev/null plotcol=0,0 plotfat=10,10
	Graph min1=0 min2=0 max1=8.28 max2=8.28 yreverse=1 symbolsz=10 symbol="*" wantaxis=n title=" " <jj2.H out=jj3.v head=/dev/null 
$	vp_Overlay jj.v jj2.v jj3.v >FIGDIR/medium.v
	vp_Overlay jj.v jj2.v jj3.v >medium.v
	vp_annotate <medium.v text=medium.ann batch=y >FIGDIR/medium.v
	
FIGDIR/medium.ps: FIGDIR/medium.v
	pstexpen FIGDIR/medium.v FIGDIR/medium.ps color=COLOR fat=FAT fatmult=FATMULT invras=y PSTEXOPTS

#define SLOPEDEPS homoisoslope.H wavelet.vert wavelet.horiz
slope.i&:  Vxseis1.i Vzseis1.i Vxseis2.i Vzseis2.i Hxseis1.i Hzseis1.i Hxseis2.i Hzseis2.i

Vxseis1.i Vzseis1.i:  SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.vert >junk  dtseis=0.008 src_1=52 src_depth=104 geo_depth=104 xseis=Vxseis1.i zseis=Vzseis1.i srcextend=0 centered=0

Vxseis2.i Vzseis2.i: SLOPEDEPS
	time  ULTIMOD type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.vert >junk  dtseis=0.008 src_1=156 src_depth=104 geo_depth=104 xseis=Vxseis2.i zseis=Vzseis2.i srcextend=0 centered=0

Hxseis1.i Hzseis1.i: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.horiz >junk  dtseis=0.008 src_1=52 src_depth=104 geo_depth=104 xseis=Hxseis1.i zseis=Hzseis1.i srcextend=0 centered=0

Hxseis2.i Hzseis2.i: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.horiz >junk  dtseis=0.008 src_1=156 src_depth=104 geo_depth=104 xseis=Hxseis2.i zseis=Hzseis2.i srcextend=0 centered=0





slope.b&:  Vxseis1.b Vzseis1.b Vxseis2.b Vzseis2.b Hxseis1.b Hzseis1.b Hxseis2.b Hzseis2.b 

Vxseis1.b Vzseis1.b: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.vert >junk  dtseis=0.008 src_1=52 src_depth=104 geo_depth=104 xseis=Vxseis1.b zseis=Vzseis1.b #srcextend=0 centered=0

Vxseis2.b Vzseis2.b: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.vert >junk  dtseis=0.008 src_1=156 src_depth=104 geo_depth=104 xseis=Vxseis2.b zseis=Vzseis2.b #srcextend=0 centered=0

Hxseis1.b Hzseis1.b: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.horiz >junk  dtseis=0.008 src_1=52 src_depth=104 geo_depth=104 xseis=Hxseis1.b zseis=Hzseis1.b #srcextend=0 centered=0

Hxseis2.b Hzseis2.b: SLOPEDEPS
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoisoslope.H par=parfile.test  <wavelet.horiz >junk  dtseis=0.008 src_1=156 src_depth=104 geo_depth=104 xseis=Hxseis2.b zseis=Hzseis2.b #srcextend=0 centered=0


test&:  homoiso.H wavelet.vert 
	time  ULTIMOD  type=4 src_type=0 ncomp=2 modeling=1 moduli=homoiso.H par=parfile.test  <wavelet.vert >junk  ncoef=0 coeff=ceoff.H dtseis=0.008 nsrc=2 src_inc=10 


/* make a two component source */
wavelet.vert: if not exist wavelet.vert
        Wavelet n1=128 d1=0.004 domain=time wavelet=ricker1 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 >wavelet.vert

wavelet.horiz: if not exist wavelet.horiz
        Wavelet n1=128 d1=0.004 domain=time wavelet=ricker1 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 junk2 junk space=n axis=2 >wavelet.horiz

wavelet.45:
        Wavelet n1=128 d1=0.004 domain=time wavelet=ricker1 tdelay=0.080 fund=15 fhigh=30. | Scale dscale=1000000. >junk
	Cp junk junk2
        Merge junk junk2 space=n axis=2 >wavelet.45

promax-velocity.H:
	Window <homoisoslope.H n3=1 >jj
	Scale <jj dscale=0.4 >jj2
	Add sqrt=y jj2 |Scale dscale=1000 >promax-velocity.H
	echo "d1=40. d2=40." >> promax-velocity.H
	CreatePHeader <promax-velocity.H >hpromax-velocity.H
	echo "attrfile=/homes/sep/martin/Model/Aniso/Recip/hpromax-velocity.H" >> promax-velocity.H

homoiso.H: BINDIR/Geninput if not exist homoiso.H
	BINDIR/Geninput nx=208 nz=208 dx=40. dz=40.  symmetry=2 vp=4.0 vs=2.3  rho=2.5   >homoiso.H

homoisorho.H: homoiso.H if not exist homoisorho.H
	Window f3=4 n3=1 <homoiso.H >rho.H
	Window f3=0 n3=4 <homoiso.H >rest.H
	Window f2=0   n2=104 <rho.H >junk1
	Window f2=104 n2=104 <rho.H >junk2
	Scale <junk2 dscale=0.9 >junk3
	Merge junk1 junk3 axis=2 space=n >rhonew.H
	Merge rest.H rhonew.H axis=3 space=n > homoisorho.H

homoisoslope.H:homoisobump.H BINDIR/Slope if not exist homoisoslope.H
	BINDIR/Slope all=1 scale=0.75 < homoisobump.H > jj.H slope=-4 yinter=416
	Cp jj.H homoisoslope.H
$	Where i1=104 w1=4 i2=156 w2=4 scale=5. all=1 <jj.H >homoisoslope.H

homoisobump.H:homoiso.H BINDIR/Where if not exist homoisobump.H
	BINDIR/Where i1=104 w1=4 i2=156 w2=4 scale=5. all=1 <homoiso.H >homoisobump.H

homoisobump2.H:homoiso.H BINDIR/Where
	BINDIR/Where all=1 i1=104 w1=4 i2=156 w2=4 scale=5. all=1 <homoiso.H >homoisobump2.H


homoisoc.H: homoiso.H
	Window f3=4 n3=1 <homoiso.H >rho.H
	Window f3=0 n3=4 <homoiso.H >rest.H
	Window f2=0   n2=104 <rest.H >junk1
	Window f2=104 n2=104 <rest.H >junk2
	Window f1=0 n1=52 <junk2 >junk2a
	Window f1=52 n1=156 <junk2 >junk2b
	Scale <junk2b dscale=0.8 >junk3
	Merge junk2a junk3 axis=1 space=n >junk3a
	Merge junk1 junk3a axis=2 space=n >new.H
	Merge new.H rho.H axis=3 space=n > homoisoc.H
$ if noise
$	Merge new.H rho.H axis=3 space=n > homoisocn.H
$	Spike n1=208 n2=208 k1=-1 k2=-1 mag=1.0 >scale.H
$	Noise gauss=0 mean=0.05 var=0.0166 rep=0 <scale.H >scalen.H
$	Merge axis=3 space=n scalen.H scalen.H scalen.H scalen.H scalen.H >scalen5.H
$	Add add=0.0,0.8 mode=product homoisocn.H  scalen5.H >homoisoc.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*
	
ortho.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 */
	Cp <junk3 >junk4
	Merge junk junk3 junk4 junk2 axis=3 space=n >ortho.H
	RM_CMD junk*

mono.H: homoiso.H
	Window f3=0 n3=1 <homoiso.H >junk.c11
	Window f3=1 n3=1 <homoiso.H >junk.c13
	Window f3=2 n3=1 <homoiso.H >junk.c33
	Window f3=3 n3=1 <homoiso.H >junk.c55
	Window f3=4 n3=1 <homoiso.H >junk.rho
	Spike n1=208 n2=208 d1=0.040 d2=0.040 nsp=0 >junk.c14
	Cp <junk.c13 >junk.c12
	Cp <junk.c14 >junk.c15
	Cp <junk.c11 >junk.c22
	Cp <junk.c13 >junk.c23
	Cp <junk.c14 >junk.c25
	Cp <junk.c14 >junk.c26
	Cp <junk.c14 >junk.c35
	Cp <junk.c55 >junk.c44
	Cp <junk.c14 >junk.c46
	Cp <junk.c55 >junk.c66
	Merge junk.c11 junk.c12 junk.c13 junk.c15 junk.c22 junk.c23 junk.c25 junk.c33 junk.c35 junk.c44 junk.c46 junk.c55 junk.c66 junk.rho axis=3 space=n >mono.H
	RM_CMD junk*

triclinic.H: homoiso.H
	Window f3=0 n3=1 <homoiso.H >junk.c11
	Window f3=1 n3=1 <homoiso.H >junk.c13
	Window f3=2 n3=1 <homoiso.H >junk.c33
	Window f3=3 n3=1 <homoiso.H >junk.c55
	Window f3=4 n3=1 <homoiso.H >junk.rho
	Cp <junk.c13 >junk.c12
	Spike n1=208 n2=208 d1=0.040 d2=0.040 nsp=0 >junk.c14
	Cp <junk.c14 >junk.c15
	Cp <junk.c14 >junk.c16
	Cp <junk.c11 >junk.c22
	Cp <junk.c13 >junk.c23
	Cp <junk.c14 >junk.c24
	Cp <junk.c14 >junk.c25
	Cp <junk.c14 >junk.c26
	Cp <junk.c14 >junk.c34
	Cp <junk.c14 >junk.c35
	Cp <junk.c14 >junk.c36
	Cp <junk.c55 >junk.c44
	Cp <junk.c14 >junk.c45
	Cp <junk.c14 >junk.c46
	Cp <junk.c14 >junk.c56
	Cp <junk.c55 >junk.c66
	Merge junk.c11 junk.c12 junk.c13 junk.c14 junk.c15 junk.c16 junk.c22 junk.c23 junk.c24 junk.c25 junk.c26 junk.c33  axis=3 space=n >jj
	Merge jj junk.c34 junk.c35 junk.c36 junk.c44 junk.c45 junk.c46 junk.c55 junk.c56 junk.c66 junk.rho axis=3 space=n >triclinic.H
	RM_CMD junk* jj*
	

isolayer.H: BINDIR/Geninput 
	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



layer.H:  BINDIR/Geninput 
	BINDIR/Geninput nx=208 nz=50 dx=40. dz=15. symmetry=2 vp=4.0 vs=3.3  rho=2.5   >junk1
	BINDIR/Geninput nx=208 nz=50 dx=40. dz=15. symmetry=2 c11=34.3 c13=10.7 c33=22.7 c44=5.4 c66=10.6 rho=2.5 >junk0
	BINDIR/Geninput nx=208 nz=50 dx=40. dz=15. symmetry=2 vp=3.0 vs=2.0  rho=2.5     >junk3
	BINDIR/Geninput nx=208 nz=50 dx=40. dz=15. symmetry=2 c11=34.3 c13=10.7 c33=22.7 c44=5.4 c66=10.6 rho=2.5 >junk2
	BINDIR/Geninput nx=208 nz=8 dx=40. dz=15.  symmetry=2 vp=4.0 vs=2.3  rho=2.5     >junk4
	Merge junk0 junk1 junk2 junk3 junk4 space=n axis=2 |Transp plane=12 >layer.H


rmjunk&:
	-Rm j*.H
	-RM_CMD junk* j*

clean&: jclean
	-RM_CMD *.x *.o  a.out core junk* j*  panel* rpanel*  *.fcm *seis{1,2}.*
	-RM_CMD zsnap xsnap wavelet.vert wavelet.horiz

veryclean&: clean
	-RM_CMD *.H wavelet* j* *seis* *snap* *.v  

%.T:	
	Taplot <% gainpanel=every  >%.T

CreatePHeader: CreatePHeader.o
	FC LDOPTS CreatePHeader.o -o CreatePHeader -lsepf -lsepattr -lsep

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

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