#SEPINC = /net/koko/bob/SEP_BOB/include

include $(SEPINC)/SEP.top


F90LIBS =-lsupersetf90 -lsepgeef90 -lsuperset -lsepauxf90 -lsepmathf90 -lsep2df90 -lsep3df90 -lsep3d -lsepf90 -lsep


RESULTSER =

RESULTCR = 
RESULTSNR = 

#UF90FLAGS = -O0 -g -traceback -fpe:0 -check all

BINDIR = ./bin
SRCDIR = ./src
RESDIR = ./fig
OBJDIR = ./obj

B = ${BINDIR}
S = ${SRCDIR}
R = ${RESDIR}
O = ${OBJDIR}

HTMLDIR = .

default1:

default:
	make -k -i default1; make -k -i default1

dn = >/dev/null
grey = Grey $(dn)

#clip = pclip=99 wantscalebar=1
clip = newclip=1 crowd=0.80 # wantscalebar=1
FAT = 3

$B/Simpvelmodel.x: ${OBJDIR}/simpvelmod.o

#################### Finite Difference Modeling 4th order in time and 10th order in space
seis.H: $B/acoustic2d_4t8x.x
	$B/acoustic2d_4t8x.x par=par.P vel=vel-homo.H snap=snap.H

#### Window and subsample the marmousi model to dx=10m dz=4m ######
vel/vel-smarm.H: ~xukai/huo/fortran/Vmodel/marmousi/marmvel.H
	< $< Window3d min2=3000 max2=7000 > vel\vel1.H
	< vel\vel1.H Interp d1out=8 d2out=10 >$@

#### Generate the simple two layer velocity model with a dip velocity boundary using c
n1 = 200
n2 = 250

vel/simpmod-flat.H: $B/SimpVelModel2d.x
	$< d1=5 d2=5 n1=${n1} n2=${n2}  > $@

vel/simpmod-dipping.H: $B/SimpVelModel2d.x
	$< d1=5 d2=5 n1=${n1} n2=${n2} slope=0.3  > $@

vel/simpmod-steep.H: $B/SimpVelModel2d.x
	$< d1=5 d2=5 n1=${n1} n2=${n2} slope=-2.0  > $@

#### Create the data file from modeling program

### Use model Marmousi
#seis-smarm.H: vel/vel-smarm.H  $B/acoustic2d_4t10x.x
#	$B/acoustic2d_4t10x.x par=par/marm1.P < $< Snap=snap-smarm.H >$@


### Use model simple
#seis-simp-flat.H: vel/simpmod-flat.H  $B/acoustic2d_4t10x.x
#	$B/acoustic2d_4t10x.x par=par/simp.P < $< Snap=snap-simp-flat.H srcwav=wavlet.H >$@

centerFreq = 15
maxFreq = 25
#generate the wavelet
#since peak freq=30Hz. the length should be at least 1/30*2/dt = 16pts

#ricker.H:
#	Wavelet n1=256 d1=0.001 domain=time wavelet=ricker2 fund=${centerFreq}\
#		 tdelay=0.128 phase=none \
#	  >tmp.H
#		Math file1=tmp.H exp='file1+0.0' >$@
#	< tmp.H Add scale=-1 | Bandpass flo=2 nplo=5 nphi=5 fhi=${maxFreq} > $@
#	< tmp2.H Bandpass flo=2 fhi={maxFreq} > $@  || this will make the zero-freq worse
#	echo peakFreq=${maxFreq} >> $@

#	< $< Bandpass flo=10 fhi=60 phase=0 >$@
delta.H: 
	Spike nsp=1 n1=64 k1=1 mag=1.0 d1=0.004 >$@
	echo freq=40 >> $@

#Impulse response of wave modeling
dseis-simp-flat.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x par/modelingSimpVel.P
	$B/Modeling2d.x par=par/modelingSimpVel.P vel=$< \
	src_snap=snap-simp-flat.H wavelet=delta.H \
	wavelet_interped=wavInterped.H recv_snap=dummy.H \
	seis_record=$@ 

ricker.H:
	cp ~/Research/matlab/ricker.20.H $@
	echo '' >> $@
	echo  'peakFreq=40' >> $@

#bench mark Fdmod by SEPlib
seis-std.H: vel/simpmod-flat.H
	Fdmod intag=$< outtag=snap-std.H hsfile=$@ mt=28 no_stdout=0 tmax=2.0 fpeak=30 xs=625 zs=500

#Do the modeling using new RTM code
seis.simp.flat.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x par/modelingSimpVel.P
	$B/Modeling2d.x par=par/RTMSimpVel.P vel=$< \
	src_snap=snap.s.f.H wavelet=ricker.H \
	recv_snap=dummy.H \
	seis_record=$@ dump=dump.H dump2=dump2.H

#Do the RTM
vel/simpmod-flat-smooth.H: vel/simpmod-flat.H
	< $< Smooth rect1=20 rect2=10 > $@

#Adj RTM
refl.s.f.H: vel/simpmod-flat.H
	<$< $B/Grad2d.x | Scale scale_to=1.0 >$@


# RTM
img.s.H imgs.s.H: vel/simpmod-flat-smooth.H $B/RTM2d.x par/RTMSimpVel.P
	$B/RTM2d.x par=par/RTMSimpVel.P vel=$< \
	src_snap=snap.s.H wavelet=ricker.H \
	recv_snap=recv.snap.s.H \
	seis_record=seis.lin.s.H RTM_image=imgs.s.H \
  	dump=dump.H final_RTM_image=img.s.H 

#Make simple uniform vel

vel/vel1.H:
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=10 d2=10 > $@

## %0.5 vel perturbation
vel/vel2.H: vel/vel1.H
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3015 d1=10 d2=10 > $@

# %33 perturb
vel/vel3.H:vel/vel1.H
	Vel n1=140 n2=300 n3=1 o3=0 d3=1 o1=0 o2=0 vc=4000 d1=10 d2=10 > $@

# -%10 perturb
vel/vel4.H:vel/vel1.H
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=2700 d1=10 d2=10 > $@


#True vel, with the anomaly
vel/vela1.H:
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 xan=2500 zan=500 \
	ranx=400 ranz=200 rankx=200 rankz=100 exank=1.0 exan=1.0 dvank=-300 \
	dvan=-300 d1=10 d2=10 > $@
#without the anamoly
vel/vela2.H: 
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=10 d2=10 > $@
vel/vela3.H:
	Vel n1=140 n2=500 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=10 d2=10 > $@

img1.H imgs1.H : seis-simp.H vel/vel1.H ricker.H $B/RTM2d.x par/beamSimp.P
	$B/RTM2d.x par=par/beamSimp.P vel=vel/vel1.H \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$< dump=wav2.H final_RTM-image=$@ RTM_image=imgs1.H

#Calculate the gradient of the local objective function

seis.lin.s.H: vel/simpmod-flat-smooth.H $B/AdjRTM2d.x par/RTMSimpVel.P 
	$B/AdjRTM2d.x par=par/RTMSimpVel.P vel=$< \
	src_snap=snap.s.H wavelet=ricker.H \
	recv_snap=recv.snap.s.H refl_model=refl.s.f.H \
	seis_record=$@ dump=dump.H 

#Single Shots
seis1.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp.P
	$B/Modeling2d.x par=par/beamSimp.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

grad1.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#Yi Luo's Method
grad1-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#classical FWI gradient
grad1-3.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=1 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#slow vel
seis3.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp.P
	$B/Modeling2d.x par=par/beamSimp.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

grad3.H: vel/vel4.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=3 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#impluse response Yi Lu
seisImp.H: vel/vel1.H ricker.H $B/Modeling2d.x 
	$B/Modeling2d.x par=par/beamSimpImpulse.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

gradImp.H: vel/vel3.H ricker.H $B/beamWET.x 
	$B/beamWET_grad.x par=par/beamSimpImpulse.P vel=$< data_obs=seisImp.H \
	wavelet=ricker.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
#Multiple Shots
seis2.H: vel/vel1.H ricker.H $B/Modeling2d.x par/beamSimp2.P
	$B/Modeling2d.x par=par/beamSimp2.P vel=$< \
	src_snap=snap1.H wavelet=ricker.H \
	seis_record=$@ dump=wav2.H

#Yi Luo's gradient
grad2-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp2.P seis2.H
	$B/beamWET_grad.x par=par/beamSimp2.P vel=$< data_obs=seis2.H \
	wavelet=ricker.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	

#Do the FWI for multiple iterations
#YL's method multiple shots
velMod2-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp2.P seis2.H
	$B/beamWET.x par=par/beamSimp2.P vel=$< data_obs=seis2.H \
	wavelet=ricker.H recv_snap=snap1B.H grad=grad2-2.H \
	objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H dump9=dump9.H> finalData=dat_@ $@

#single shot
velMod1-2.H: vel/vel3.H ricker.H $B/beamWET.x par/beamSimp.P seis1.H
	$B/beamWET.x par=par/beamSimp.P vel=$< data_obs=seis1.H \
	wavelet=ricker.H recv_snap=snap1B.H grad=grad2-2.H \
	objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H dump9=dump9.H >$@ grad=grads1-2.H finalData=dat_@

#regular 1D DaveHale's boundary
#There might be coding bugs. It is less effect than boundary_type=1
seis-simp-flat-regular.H: vel/simpmod-flat.H ricker.H $B/Modeling2d.x par/modelingSimpVel.P
	$B/Modeling2d.x par=par/modelingSimpVel.P vel=$< \
	src_snap=snap-simp-flat-regular.H wavelet=ricker.H \
	wavelet_interped=wavInterped.H recv_snap=dummy.H \
	seis_record=$@ boundary_type=3

#Do the RTM
image-simp-flat.H: vel/simpmod-flat-smooth.H $B/RTM2d.x par/RTMSimpVel.P
	$B/RTM2d.x par=par/RTMSimpVel.P vel=$< \
	src_snap=src-snap-simp.H wavelet=ricker2.H \
	wavelet_interped=wavInterped.H recv_snap=recv-snap-simp.H \
	seis_record=seis-simp-flat.H boundary_type=1 \
	RTM_image=$@ dump1=dump1.H dump2=dump2.H

#build the refl model, maybe calculate the gradient is a good option
refl-simpmod-flat.H: vel/simpmod-flat.H
	<$< $B/Grad2d.x | Scale scale_to=1.0 >$@

#RTM modeling
seis-lin-simp-flat.H: vel/simpmod-flat-smooth.H $B/adjRTM2d.x par/modelingSimpVel.P 
	$B/adjRTM2d.x par=par/RTM_modelingSimpVel.P vel=$< \
	src_snap=src-snap-simp.H wavelet=ricker2.H \
	wavelet_interped=wavInterped.H recv_snap=recv-snap-simp.H \
	seis_record=seis-lin-simp-flat.H boundary_type=1 \
	refl_model=refl-simpmod-flat.H dump1=dump1.H dump2=dump2.H

cropped-image-simp-flat.H: image-simp-flat.H
	< $< Window3d min1=300.0 > $@

imag-stack.H: cropped-image-simp-flat.H
	< $< Stack axis=3 >$@

###  compare Ali
vel/v1.H:
	Vel n1=481 n2=601 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3000 d1=25 d2=25 > $@

vel/v2.H: vel/v1.H
	Vel n1=481 n2=601 n3=1 o3=0 d3=1 o1=0 o2=0 vc=2900 d1=25 d2=25 > $@


vel/v3.H: vel/v1.H
   Vel n1=481 n2=601 n3=1 o3=0 d3=1 o1=0 o2=0 vc=3100 d1=25 d2=25 > $@


wav.H:
   Wavelet n1=550 d1=0.0038 wavelet=ricker2 phase=0 fund=15 delay=0.057 > $@
   
seisa1.H: vel/v1.H ricker.H $B/Modeling2d.x par/beamSimp.P
	$B/Modeling2d.x par=par/ vel=$< \
	src_snap=snap1.H wavelet=wav.H \
	seis_record=$@ dump=wav2.H

grada3.H: vel/v3.H ricker.H $B/beamWET.x par/beamSimp.P
	$B/beamWET_grad.x par=par/beamSimp.P vel=$< data_obs=seisa1.H \
	wavelet=wav.H src_snap=snap1A.H recv_snap=snap1B.H \
	grad=$@ objFunc=2 dump1=dump1.H dump2=dump2.H dump3=dump3.H dump4=dump4.H \
	dump5=dump5.H dump6=dump6.H dump7=dump7.H dump8=dump8.H	
	

gcn = Grey pclip=100 color=j wantscalebar=1 gainpanel=e

snap:
	< snap-simp-flat.H Window n3=1 f3=10 > temp.H;
	${gcn} < temp.H | Tube

%.ps: %.v
	pstexpen $< $@ color=n fat=${FAT} fatmult=1.5 invras=n
%.pdf: %.ps
	epstopdf $<

clean: jclean
	-rm ${BINDIR}/*.* -f
	-rm ${OBJDIR}/*.* -f
	-rm *.mod *.o *.H -f
	-rm $R/*.* -f
#	@(cd $(srcdir);${MAKE} clean)

ss:
	rm ${BINDIR}/* -f

.PHONY: clean ss default1

#Type 'make' to build all the figures
#And then type 'scons' to build the report. 

include $(SEPINC)/SEP.bottom


