#define FIGDIR ../Fig
#define BINDIR ../Bin/MTYPE

/*this is in (10e9)N/m^2*/
#define TISHALE c11=29.81 c13=10.35 c33=20.7 c44=5.18 c66=9.73 rho=2.3
#define MINTP min2=0 max2=4.096 min1=-0.675 max1=0.675 yreverse=1 dash=0,3 plotfat=0,6
#define MINTT min2=0 max2=4.000 min1=-4.12 max1=14.44  yreverse=1 dash=0,3 plotfat=0,6
#define MINGRAPH min2=0 max2=1.000
#define PSARGS fat=6
#define NOHEAD >/dev/null

#include <SEP.defs>

figures&: FIGDIR/ppcalibrate.v  FIGDIR/extraiso.v

annotate&:
	vp_annotate batch=y <FIGDIR/ppcalibrate.v                                       text=ppcalibrate.annotate >FIGDIR/ppcalibrate.va

/* overlay isotropic tau - p  curves  */
overlay.iso.tp&: pz.iso.tp.V px.iso.tp.V

iso.ttpp %.iso.ext.pp: ../Psource/CMP/%seis.iso.cmp.SUB.SL
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=3.6 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > iso.ttpp extract=%.iso.ext.pp

iso.ttss %.iso.ext.ss: ../Psource/CMP/%seis.iso.cmp.SUB.SL
        TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL> iso.ttss extract=%.iso.ext.ss

iso.ttps %.iso.ext.ps: ../Psource/CMP/%seis.iso.cmp.SUB.SL
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel1=3.6 vel2=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > iso.ttps extract=%.iso.ext.ps

iso.ttsp %.iso.ext.sp: ../Psource/CMP/%seis.iso.cmp.SUB.SL
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel2=3.6 vel3=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > iso.ttsp extract=%.iso.ext.sp

%.iso.tp.V : BINDIR/TXTPcurve ../Psource/CMP/%seis.iso.cmp.SUB.SL iso.ttpp iso.ttss iso.ttps iso.ttsp
	Taplot <../Psource/CMP/%seis.iso.cmp.SUB.SL |Ta2vplot >j1.V
	Graph <iso.ttpp   MINTP >j2.V
	Graph <iso.ttss   MINTP >j3.V
	Graph <iso.ttps   MINTP >j4.V
	Graph <iso.ttsp   MINTP >j5.V
	Vppen <j1.V j2.V j4.V  erase=once vpstyle=n >%.iso.tp.V

/*
%.iso.tp.V : BINDIR/TXTPcurve ../Psource/CMP/%seis.iso.cmp.SUB.SL
	Taplot <../Psource/CMP/%seis.iso.cmp.SUB.SL |Ta2vplot >j1.V
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=3.6 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > ttpp extract=%.iso.ext.pp
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL> ttss extract=%.iso.ext.ss
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel1=3.6 vel2=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > ttps extract=%.iso.ext.ps
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel2=3.6 vel3=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.iso.cmp.SUB.SL > ttsp extract=%.iso.ext.sp
	Graph <ttpp   MINTP >j2.V
	Graph <ttss   MINTP >j3.V
	Graph <ttps   MINTP >j4.V
	Graph <ttsp   MINTP >j5.V
	Vppen <j1.V j2.V j4.V  erase=once vpstyle=n >%.iso.tp.V
	Vppen <j1.V j3.V j5.V  erase=once vpstyle=n >%.iso.tp2.V
*/

/* overlay  homogeneous isotropic values for calibration */

$ do same as abov for iso and reproduce

%.homoiso.tp.V %.homoiso.ext.pp: BINDIR/TXTPcurve ../Psource/CMP/%seis.homoiso.cmp.SUB.SL
	Taplot <../Psource/CMP/%seis.homoiso.cmp.SUB.SL |Ta2vplot >j1.V
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=3.6 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.homoiso.cmp.SUB.SL > ttpp extract=%.homoiso.ext.pp
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.homoiso.cmp.SUB.SL> ttss extract=%.homoiso.ext.ss
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel1=3.6 vel2=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.homoiso.cmp.SUB.SL > ttps extract=%.homoiso.ext.ps
	TXTPcurve taup=1 wantextract=1 extrwin=10 velocity=-1 vel2=3.6 vel3=2.06 z0=1.0 tshift=0.08 <../Psource/CMP/%seis.homoiso.cmp.SUB.SL > ttsp extract=%.homoiso.ext.sp
	Graph <ttpp   MINTP >j2.V
	Graph <ttss   MINTP >j3.V
	Graph <ttps   MINTP >j4.V
	Graph <ttsp   MINTP >j5.V
	Vppen <j1.V j2.V j4.V  erase=once vpstyle=n >%.homoiso.tp.V
	Vppen <j1.V j3.V j5.V  erase=once vpstyle=n >%.homoiso.tp2.V


FIGDIR/ppcalibrate.v:  pz.iso.tp.V px.iso.tp.V pz.homoiso.tp.V px.homoiso.tp.V 
	Window f1=10 n1=1 f2=73 n2=56 < pz.iso.ext.pp > jjz1
	Window f1=10 n1=1 f2=73 n2=56 < px.iso.ext.pp > jjx1
	Window f1=10 n1=1 f2=73 n2=56 < pz.homoiso.ext.pp > jjz0
	Window f1=10 n1=1 f2=73 n2=56 < px.homoiso.ext.pp > jjx0
	Pad < pz.iso.ext.pp  n1=0 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2 >jj
	Pad < pz.homoiso.ext.pp  n1=0 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2 >jjj
	Add mode=divide jj jjj | Window f2=74 n2=55  n1=1 f1=3 >jref 
	Window n2=1 n1=72 <../Zoeppritz/iso.reflectivity >jj1
	echo "d1=0.002777" >> jj1
	Graph <jj1  pad=1 min2=0 max2=1.0 min1=0 max1=0.2  title="P-P Reflectivity Calibration" label1="Horizontal Slowness [s/km]" label2="Energy Coefficient" plotcol=6 >jzoepp.ref.V 
	Graph <jref pad=1 min2=0 max2=1.0 min1=0 max1=0.2  title=" " label1=" "         label2=" "   plotcol=5  dash=3 plotfat=6 >jiso.ref.V
	Vppen <jzoepp.ref.V jiso.ref.V erase=once vpstyle=n out=FIGDIR/ppcalibrate.v  NOHEAD


FIGDIR/extraiso.v: pz.iso.ext.pp pz.homoiso.ext.pp 
	Taplot <pz.iso.ext.pp |Ta2vplot title="Extracted Reflected Amplitudes" label2="Horiz.Slowness" label1="Differential travel time curve from centroid" out=FIGDIR/pz.iso.ext.pp.v
	Taplot <pz.homoiso.ext.pp |Ta2vplot title="Extracted Dual Amplitudes" label2="Horiz.Slowness" label1="Differential travel time curve from centroid" out=FIGDIR/pz.homoiso.ext.pp.v  NOHEAD
	vp_OverUnderIso FIGDIR/pz.iso.ext.pp.v FIGDIR/pz.homoiso.ext.pp.v >FIGDIR/extraiso.v


pseis.%.tp.pp:   pxseis.%.tp pzseis.%.tp  pxseis.homoti.tp pzseis.homoti.tp
	Pad < px.%.ext.pp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxnumerator
	Pad < pz.%.ext.pp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjznumerator
	Pad < pz.homoti.ext.pp  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjzdivisor
	Pad < px.homoti.ext.pp  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxdivisor	
	Add mode=sum jjxdivisor jjzdivisor >jj
	Add mode=sum jjxnumerator jjznumerator >jjj
	Add mode=divide  jjj jj > pseis.%.tp.pp 
$	Window <pseis.%.tp.pp f1=10 n1=1 f2=85  n2=59 |Graph >pseis.%.tp.pp.V

sseis.%.tp.ss:   sxseis.%.tp szseis.%.tp /* sx.homoti.tp sz.homoti.tp  */
	Pad < px.%.ext.ss      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxnumerator
	Pad < pz.%.ext.ss      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjznumerator
	Pad < pz.homoti.ext.ss  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjzdivisor
	Pad < px.homoti.ext.ss  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxdivisor	
	Add mode=sum jjxdivisor jjzdivisor >jj
	Add mode=sum jjxnumerator jjznumerator >jjj
	Add mode=divide  jjj jj > pseis.%.tp.ss 
$	Window <pseis.%.tp.ss f1=10 n1=1 f2=85  n2=59 |Graph >pseis.%.tp.ss.V

sseis.%.tp.sp:   sxseis.%.tp szseis.%.tp /* sxseis.homoti.tp szseis.homoti.tp  */
	Pad < px.%.ext.sp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxnumerator
	Pad < pz.%.ext.sp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjznumerator
	Pad < pz.homoti.ext.sp  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjzdivisor
	Pad < px.homoti.ext.sp  n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxdivisor	
	Add mode=sum jjxdivisor jjzdivisor >jj
	Add mode=sum jjxnumerator jjznumerator >jjj
	Add mode=divide  jjj jj > pseis.%.tp.sp 
$	Window <pseis.%.tp.sp f1=10 n1=1 f2=85  n2=59 |Graph >pseis.%.tp.sp.V



/* no clever calibration is  possible for converted waves, 
                        so we only give reflrected amplitude */
/* that is ok, if we compare azimuthally  */
pseis.%.tp.ps:  pxseis.%.tp pzseis.%.tp  
	 Pad < px.%.ext.ps      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxnumerator
	 Pad < pz.%.ext.ps      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjznumerator
$	Add mode=sum jjxdivisor jjzdivisor >jj
	Add mode=sum jjxnumerator jjznumerator >pseis.%.tp.ps
$        Window < pseis.%.tp.ps n1=1 f1=10 f2=88 n2=73 |Graph title="P-S Reflected Amplitudes (%)" label1="Horizontal Slowness"   label="Magnitude" out=FIGDIR/pseis.%.tp.ps.v        NOHEAD

pseis.%.tp.sp: pxseis.%.tp pzseis.%.tp 
	 Pad < px.%.ext.sp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjxnumerator
	 Pad < pz.%.ext.sp      n1=128 | Rtoc| Ft3d sign1=1 sign2=0 sign3=0 |Cabs2         > jjznumerator
$	Add mode=sum jjxdivisor jjzdivisor >jj
	Add mode=sum jjxnumerator jjznumerator >pseis.%.tp.sp
$        Window < pseis.%.tp.sp n1=1 f1=10 f2=88 n2=73 |Graph title="S-P Reflected Amplitudes (%)" label1="Horizontal Slowness"   label="Magnitude" out=FIGDIR/pseis.%.tp.sp.v        NOHEAD

/* compare the azimuthal amplitudes   in the converted P S  case */



/* The final comarison figures   */

/* some wimpy tests */
/* isotropic moveout curves in tau p  */
doitiso&: 
	TPisocurve vp=3.6 vs=2.06 z0=1.0 n1=512 n2=360 d1=0.008 o2=-0.675 d2=0.00375 > taupcurve

doit&: 
	Window n1 <ptigroup.h >ptigroup2.h
	TXTPcurve vel=ptigroup2.h z0=1.0 n1=2000 n2=464 d1=0.002 o2=-2.06 d2=0.020 velocity=0 > curves

doit2&: 
	Window n1 <ptigroup.h >ptigroup2.h
	Extract vel=ptigroup2.h z0=1.0 n1=2000 n2=464 d1=0.002 center=104 d2=0.020 velocity=0 > junk


/* Here I am generating the velocity curves using Ellipse which lives in  */
/* Model/Aniso/2D/Stolt                                                   */
/* what = 1->p  2->sv  3->sh  */

ptigroup.h pellip.h:
	Ellipse TISHALE n1=1000 what=1 >junk
	Window <junk f2=1 n2=1>ptigroup.h            /*TI group*/
	Window <junk f2=3 n2=1>pellip.h            /* ELLIP group */

svtigroup.h svellip.h:
	Ellipse TISHALE n1=1000 what=2 >junk
	Window <junk f2=1 n2=1>svtigroup.h            /*TI group*/
	Window <junk f2=3 n2=1>svellip.h            /* ELLIP group */

shtigroup.h shellip.h:
	Ellipse TISHALE n1=1000 what=3 >junk
	Window <junk f2=1 n2=1>shtigroup.h            /*TI group*/
	Window <junk f2=3 n2=1>shellip.h            /* ELLIP group */

ptiphase.h:
	Ellipse TISHALE  n1=1000 what=1 >junk
	Window <junk f2=0 n2=1>ptiphase.h
	Window <junk f2=2 n2=1>ptiphase.ell.h

svtiphase.h:
	Ellipse TISHALE  n1=1000 what=2 >junk
	Window <junk f2=0 n2=1>svtiphase.h
	Window <junk f2=2 n2=1>svtiphase.ell.h

shtiphase.h:
	Ellipse TISHALE  n1=1000 what=3 >junk
	Window <junk f2=0 n2=1>shtiphase.h
	Window <junk f2=2 n2=1>shtiphase.ell.h

clean&:
	-RM_CMD *.x *.o *.f  a.out core *.T *.V j* ttpp ttps ttss ttsp

Ellipse: Ellipse.o ti.o ti2ell.o ellip.o phase2group.o 
	LD LDOPTS Ellipse.o ti.o ti2ell.o ellip.o phase2group.o SEPLIB FLIB2 SYSLIB -o BINDIR/Ellipse

Ti2ell: Ti2ell.o 
	LD LDOPTS Ti2ell.o SEPLIB FLIB2 SYSLIB -o BINDIR/Ti2ell
	
%.ps : %.v
	pstexpen %.v %.ps

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