# . .  Coordinate System and 
#
# Test1.rays.HH
#
# . . and velocity
#
# Test1.vel.HH
#
nt=2048
dt=0.004
ot=0.
rweparx = nref=3 verbose=1 forward=1 kinematic=1. min_region_pct=.1\
		 del_dist=.05 verb=0 niter_lloyd=20

# Cartesian Output
nx1=881
ox1=-11000
dx1=25

ny1=1
oy1=0
dy1=50

nz1=441
dz1=25
oz1=0
XBIN = ./Bin/
CCoutx  =azn=$(nz1) azo=$(oz1) azd=$(dz1)\
	 axn=$(nx1) axo=$(ox1) axd=$(dx1)\
	 ayn=$(ny1) ayo=$(oy1) ayd=$(dy1)\
	 nsx=3 nsz=3 nsy=0 adj=1

HWTgraph = Graph crowd=${crowd} wantaxis=n title=" " label1=" " labe
CPGraph = min1=-8000 max1=5000 min2=0  max2=8000 crowd=0.8 labelsz=12
CPGrey  = min2=-8000 max2=8000 min1=0  max1=8000 crowd=0.8 pclip=99.5 labelsz=12
CPRGraph = min1=2000 max1=12000 min2=0 max2=400 crowd=0.8 labelsz=12
CPRGrey  = min2=2000 max2=12000 min1=0 max1=400 crowd=0.8 pclip=99.75 labelsz=12

CPnzx=400
CPozx=0 
CPdzx=1

CPnxx=500
CPoxx=200
CPdxx=20

CPcartparxxxxx= a1n=$(CPnzx) a1d=$(CPdzx) a1o=$(CPozx) \
	    a2n=$(CPnxx) a2d=$(CPdxx) a2o=$(CPoxx) 

Matlab = /net/koko/share/local/matlab_old/bin/matlab

Test1.Vel.HH Test1.rays.HH:  
	$(MATLAB) < ./matlab/Test1.m &


Tmp.H:
	Vel vc=1498 n1=441 n2=881 n3=1 d3=1 o3=0 o1=0 o2=-11000 d1=25 d2=25 > $@
TVel.H:   Tmp.H $(XBIN)/Vxz.x
	< Tmp.H $(XBIN)/Vxz.x vo=1500 gx=0 gz=0.35 >$@

Rtest.H: TVel.H
	<TVel.H Hwt2d xsou=1000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t1.H
	<TVel.H Hwt2d xsou=1500 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t2.H
	<TVel.H Hwt2d xsou=2000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t3.H
	<TVel.H Hwt2d xsou=2500 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t4.H
	<TVel.H Hwt2d xsou=3000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t5.H
	<TVel.H Hwt2d xsou=3500 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t6.H
	<TVel.H Hwt2d xsou=4000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t7.H
	<TVel.H Hwt2d xsou=4500 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t8.H
	<TVel.H Hwt2d xsou=5000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t9.H
	Cat t[123456789].H axis=2 space=n > $@
Rtest2.H: TVel.H
	<TVel.H Hwt2d xsou=1000 zsou=0 oT=0 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t1.H
	<TVel.H Hwt2d xsou=1500 zsou=0 oT=1 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t2.H
	<TVel.H Hwt2d xsou=2000 zsou=0 oT=2 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t3.H
	<TVel.H Hwt2d xsou=2500 zsou=0 oT=3 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t4.H
	<TVel.H Hwt2d xsou=3000 zsou=0 oT=4 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t5.H
	<TVel.H Hwt2d xsou=3500 zsou=0 oT=5 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t6.H
	<TVel.H Hwt2d xsou=4000 zsou=0 oT=6 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t7.H
	<TVel.H Hwt2d xsou=4500 zsou=0 oT=7 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t8.H
	<TVel.H Hwt2d xsou=5000 zsou=0 oT=8 nT=350 dT=0.01 oG=-42 nG=1 dG=1> t9.H
	Cat t[123456789].H axis=2 space=n > $@

T1RAYS.H T1rays.H: Test1.rays.HH
	cp Test1.rays.HH $@
	echo " d1=10" >> $@
	cp Test1.rays.HH tt.H
	echo "d1=10" >> Test1.rays.HH
	< Test1.rays.HH Window3d n3=1 f3=0 > xx.H
	< Test1.rays.HH Window3d n3=1 f3=1 > zz.H
	Cmplx xx.H zz.H > T1RAYS.H
	Cat axis=3 space=n xx.H zz.H > T1rays.H

NVel.%.H tNVel.%.H:  TVel.HH $(XBIN)/CC2RC.x $(XBIN)/Vxz.x
	< Test1.Vel.HH $(XBIN)/Vxz.x vo=1500 gx=0 gz=0.$* | Pad beg2=1000 extend=1 > t$@ 
	< t$@ $(XBIN)/CC2RC.x rays=T1RAYS.H > $@

Test1.data.H:                   TVel.HH
	Get n2 n3 o2 o3 d2 d3 < TVel.HH >par.P
	Spike n1=2048 o1=0. d1=$(dt) par=par.P \
		nsp=4 k1=100,200,300,400 k2=0,0,0,0 k3=0,0,0,0 mag=2,1.9,1.8,1.7|\
		Smooth rect1=2 repeat=2 tridiag=y maxmem=1000000000  |\
		Bandpass flo=2.5 fhi=40. > $@

offset = n1=2048 n2=1024 n3=1 d2=10 d3=10 b=0.0005 c=0.0
Test2.data.H:    $(XBIN)/plane.x
		 $(XBIN)/plane.x d=1.5       $(offset) > t$@
	< t$@    $(XBIN)/plane.x d=2.2 add=y $(offset) > tt$@
	< tt$@   $(XBIN)/plane.x d=2.9 add=y $(offset) > ttt$@
	< ttt$@  $(XBIN)/plane.x d=3.6 add=y $(offset) > tttt$@
	< tttt$@ $(XBIN)/plane.x d=0.8 add=y $(offset) > ttttt$@
	< ttttt$@ Bandpass flo=1.5 fhi=50 > $@
	\rm t$@ tt$@ ttt$@ tttt$@


offset = n1=4096 n2=924 n3=1 o1=-1. d2=5 d3=10 b=0.0005 c=0.0 a=-1
Test.data.H:
		 $(XBIN)/plane.x d=-0.2       $(offset) b=0.00045 > t$@
	< t$@    $(XBIN)/plane.x d=-0.6 add=y $(offset) b=0.00045 > tt$@
	< tt$@   $(XBIN)/plane.x d=-1.0 add=y $(offset) b=0.00045 > ttt$@
	< ttt$@  $(XBIN)/plane.x d=-1.4 add=y $(offset) b=0.00045 > tttt$@
	< tttt$@ $(XBIN)/plane.x d=-1.8 add=y $(offset) b=0.00045 > ttttt$@
	< ttttt$@ Bandpass flo=1.5 fhi=50 > q1$@
		 $(XBIN)/plane.x d=-2.2       $(offset) b=0.00045 > t$@
	< t$@    $(XBIN)/plane.x d=-2.6 add=y $(offset) b=0.00045 > tt$@
	< tt$@   $(XBIN)/plane.x d=-3.0 add=y $(offset) b=0.00045 > ttt$@
	< ttt$@  $(XBIN)/plane.x d=-3.4 add=y $(offset) b=0.00045 > tttt$@
	< tttt$@ $(XBIN)/plane.x d=-3.8 add=y $(offset) b=0.00045 > ttttt$@
	< ttttt$@ Bandpass flo=1.5 fhi=50 > q2$@
	Add q1$@ q2$@ > p$@
	Pad < p$@ beg2=100  |\
		 Smooth rect2=3 rect1=3 repeat=5 tridiag=y maxmem=3000000000000 > $@ 	
	\rm t$@ tt$@ ttt$@ tttt$@
	echo "o2=2000 d2=10" >> $@

weiparxxx= min5=1.5 max5=27.5 j5=1 squeeze=n

ww.Test1.H: Test.data.H
	<   Test.data.H  Transf wei=y verb=y >t$@
	<  t$@ Window3d $(weiparxxx) squeeze=n| Transp plane=35 > $@
	\rm t$@ 
rweparx = nref=3 verbose=1 forward=1 kinematic=1. min_region_pct=.1\
		 del_dist=.05 verb=0 niter_lloyd=20
CCoutx  = azn=$(nz1) azo=$(oz1) azd=$(dz1)\
         axn=$(nx1) axo=$(ox1) axd=$(dx1)\
         ayn=$(ny1) ayo=$(oy1) ayd=$(dy1)\
         nsx=2 nsz=2 nsy=0 adj=1
HWTgraph = Graph crowd=${crowd} wantaxis=n title=" " label1=" " labe
CPGraph = min1=-8000 max1=8000 min2=0 max2=8000 crowd=0.8
CPGrey  = min2=-8000 max2=8000 min1=0 max1=8000 crowd=0.8
CPnzx=1600
CPozx=-8000
CPdzx=10
CPnxx=600
CPoxx=0
CPdxx=10 

CCtest1.H: Test1.mig.H $(B)/C2Rsinc3D.x                Test1.rays.HH
	<  Test1.mig.H $(B)/C2Rsinc3D.x $(CCoutx) rays=Test1.rays.HH > $@ 

Test1.mig.H: ww.Test1.H  Test1.Vel.HH  $(XBIN)/RWE3D_source.x
	<   ww.Test1.H $(B)/RWE3D_source.x rays=Test1.rays.HH \
		vel=Test1.Vel.HH swf=ww.Test1.H $(rweparx)> $@

CCout2  = azn=1100 azd=10 azo=0 axn=2201 axd=10 axo=-11000 adj=1 norm=1 nsx=2 nsz=2

Polar.v: Test1.mig.H Test1.rays.HH TVel.HH $(XBIN)/Lin2.x
	< Test1.rays.HH Window3d n4=2 f4=0 j4=2 > tt.H
	< tt.H Window3d n3=1 f3=0 f2=10 > xx.H
	< tt.H Window3d n3=1 f3=1 f2=10 > zz.H
	Cmplx xx.H zz.H > r.H	
	<  r.H Window3d j1=1 j2=50  > kk1.H
	< kk1.H ${HWTgraph} $(dn) out=jr.v plotcol=7 $(CPGraph) wantaxis=0 
	<  r.H Window3d j1=20 j2=1 | Transp > kk2.H
	< kk2.H ${HWTgraph} label2="Depth (m)" $(CPGraph) \
	label1="Distance (m)" labelrot=n wantaxis=n $(dn) out=jw.v plotcol=7
	vp_Overlay jr.v jw.v > ul.v
	< TVel.HH Grey newclip=1 bpclip=0 epclip=100 out=uu.v $(dn) title=" "\
		label1="Depth (m)" label2="Distance (m)" $(CPGrey) wantaxis=y \
		allpos=y newclip=1 labelsz=11 title="Physical Domain" titlesz=14 \
		color= 
	vp_Overlay uu.v ul.v > UL.v

	< Test1.Vel.HH Scale rscale=0.001 > t1.H
	echo "d1=1. o1=0. o2=0 d2=15" >> t1.H
	< t1.H Grey newclip=1 out=ur.v $(dn)  \
		label1="Extrapolation Step" label2="Distance (m) " title="Riemannian Domain " \
		$(CPRGrey) min2=2000 wantaxis=y color= allpos=y labelsz=11 titlesz=14
	$(XBIN)/Lin2.x axis=2 a2n=11 a2d=1000 a2o=2000 a1n=5 a1d=100 a1o=0 >xx.H
	$(XBIN)/Lin2.x axis=1 a2n=11 a2d=1000 a2o=2000 a1n=5 a1d=100 a1o=0 >zz.H 
	Cmplx xx.H zz.H > r.H
	< r.H ${HWTgraph} $(dn) out=jr.v plotcol=7 $(CPRGraph) wantaxis=0 
	< r.H Transp | ${HWTgraph} label2="Depth (m)" $(CPRGraph) \
	label1="Distance (m)" labelrot=n wantaxis=n $(dn) out=jw.v plotcol=7
	vp_Overlay jr.v jw.v > ul.v
	\rm jr.v jw.v 
	vp_Overlay ur.v ul.v > UR.v

	< CCtest1.H Clip clip=-500 chop=less |\
		Grey max1=8000 min2=-8000 max2=8000 out=BL.v $(dn)  \
		label1="Depth (m)" label2="Distance (m) " \
		title=" " $(CPGrey) \
		wantaxis=y pclip=99.75 labelsz=11

	echo "d1=1" >> Test1.mig.H
	< Test1.mig.H Clip clip=-500 chop=less |Grey out=BR.v $(dn)  \
		label1="Extrapolation Step" label2="Distance (m)" \
		title=" " crowd=0.8 wantaxis=y \
		pclip=100 labelsz=11 $(CPRGrey) min2=2000

	vp_SideBySideAniso UL.v UR.v > top.v
	vp_SideBySideAniso BL.v BR.v > bottom.v
	vp_OverUnderAniso top.v bottom.v > $@
#	\rm -f tt.H xx.H zz.H r.H kk1.H kk2.H ul.v jr.v jw.v V1.H V2.H V3.H \
##		uu.v ul.v UL.v lb.v tt.v BL.v t1.H t1.H ur.v\
#		ul.v UR.v top.v bottom.v

############################################################
#
# . . Gaussian Slowness example
#
Gvel.H: $(XBIN)/Vxz.x
	Gauss n1=501 o1=0. d1=10 n2=1000 o2=-5000 d2=10 ngauss=1 \
		maxvel=1000 velback=300 max1=500 max2=0     var=20000 > g1.H
	Gauss n1=501 o1=0. d1=10 n2=1000 o2=-5000 d2=10 ngauss=1 \
		maxvel=-600 velback=300 max1=1000 max2=-1000 var=20000 > g2.H
	Gauss n1=501 o1=0. d1=10 n2=1000 o2=-5000 d2=10 ngauss=1 \
		maxvel=-1200 velback=400 max1=1500 max2=1500 var=20000 > g3.H
	Math file1="g1.H" file2="g2.H" file3="g3.H" exp="file1+file2+file3" > t$@
	< t$@ $(XBIN)/Vxz.x vo=1000 gx=0.0 gz=0.2 > tt$@
	Math file1="t$@" file2="tt$@" exp="file1+file2" > $@

Grays.H:  Gvel.H
	< Gvel.H Hwt2d xsou=0 zsou=0 oT=0. nT=321 dT=0.0046875\
		 oG=-80 nG=321 dG=0.5 > $@ 

CPGraph2 = min2=0 max2=4000 min1=-4000 max1=4000 crowd=$(GEc)
CPGrey2 = min1=0 max1=4000 min2=-4000 max2=4000 crowd=$(GEc)



Geikon.H: Gvel.H
	FMeikonal < Gvel.H  order=3 zshot=0 xshot=0 yshot=0 |\
		 Smooth rect1=2 rect2=2 repeat=2> $@


GSrays.HH GSrays2.HH: Geikon.H
	matlab < ./matlab/Cnt.m

GSmesh.H: GSrays.HH
	< GSrays.HH $(XBIN)/PreRegRays.x > $@

GEc=.75
GSrays.v: GSmesh.H
	< GSmesh.H Window3d j1=1 j2=4 > kk1.H
	< kk1.H ${HWTgraph} $(dn) out=jr.v plotcol=7 $(CPGraph2) wantaxis=0 
	< GSmesh.H Window3d j1=25 j2=1 | Transp > kk2.H
	< kk2.H ${HWTgraph} label2="Depth (m)" $(CPGraph2) \
	label1="Distance (m)" labelrot=n wantaxis=y $(dn) wantaxis=n out=jw.v plotcol=7
	vp_Overlay jr.v jw.v > t1.v
	< Gvel.H Grey color= newclip=1 bpclip=0 epclip=100 out=t2.v $(dn) \
		$(CPGrey2) crowd=$(GEc) title=" " label1="Depth (m)" \
		 label2="Distance (m)" wantaxis=y labelsz=10 polarity=-1	
	vp_Overlay t2.v t1.v > $@

GGeikon.v: Geikon.H Gvel.H
	< Geikon.H Contour nc=16 dc=0.09750 out=t1.v $(dn) $(CPGrey2) \
		plotcol=7 crowd=$(GEc) wantaxis=n title=" " label1=" " label2=" "
	< Gvel.H Grey color= newclip=1 bpclip=0 epclip=100 out=t2.v $(dn) \
		$(CPGrey2) crowd=$(GEc) title=" " label1="Depth (m)"\
		labelsz=10  label2="Distance (m)" polarity=-1
	vp_Overlay t2.v t1.v > $@

GGrays.v:  Grays.H Gvel.H
	< Grays.H Window3d j1=1 j2=4 > kk1.H
	< kk1.H ${HWTgraph} $(dn) out=jr.v plotcol=7 $(CPGraph2) wantaxis=0 
	< Grays.H Window3d j1=25 j2=1 | Transp > kk2.H
	< kk2.H ${HWTgraph} label2="Depth (m)" $(CPGraph2) \
	label1="Distance (m)" labelrot=n wantaxis=n $(dn) out=jw.v plotcol=7
	vp_Overlay jr.v jw.v > t1.v
	< Gvel.H Grey color= newclip=1 bpclip=0 epclip=100 out=t2.v $(dn) \
		$(CPGrey2) crowd=$(GEc) title=" " label1="Depth (m)" \
		labelsz=10 label2="Distance (m)" polarity=-1
	vp_Overlay t2.v t1.v > $@
$(R)/GaussExample.v: GSrays.v GGrays.v GGeikon.v
	vp_OverUnderAniso GGrays.v GGeikon.v GSrays.v > $@