!--------------------------------------------------------------------
! DESCRPTION:
!	program to create a 3D velocity model with n plane layers
!	dipping at arbitrary angles with respect to x and y
!--------------------------------------------------------------------	
!	Input:
!		nl   : number of velocity layers (reflectors +1)
!		v    : array[nl] of velocities (m/s)
!		dipx : array[nl-1] of reflectors dips in x-direction (deg)
!		dipy : array[nl-1] of reflectors dips in y-direction (deg)
!		zref0: array[nl-1] of reflectors depth at first cmp (m)
!		nz   : number of depth steps
!		dz   : depth sampling interval (m)
!		oz   : minimum depth (m)
!		ncmpx: number of cmps in x-direction
!		dcmpx: sampling interval of cmps in x-direction (m)
!		ocmpx: first cmp in x-direction (m)
!		ncmpy: number of cmps in y-direction
!		dcmpy: sampling interval of cmps in y-direction (m)
!		ocmpy: first cmp in y-direction (m)
!		verb : =1 for verbose output
!
!	EXAMPLES:
!		make_vels.x nl=2 v=1500,2500 dipx=0,0 dipy=0,0 \
!			zref0=500 nz=1200 dz=5 oz=0 ncmpx=100 dcmpx=20 ocmpx=2000 \
!			ncmpy=1 dcmpy=20 ocmpy=0 verb=1 > cmp_vels_flat.H
!---------------------------------------------------------------------
program make_vels

	use sep
	use constants

	implicit none

	! program arguments
	integer 						:: nl
	real,allocatable,dimension(:) 	:: v
	real,allocatable,dimension(:) 	:: dipx
	real,allocatable,dimension(:) 	:: dipy
	real,allocatable,dimension(:) 	:: zref0
	integer							:: nz
	real							:: dz
	real							:: oz
	integer							:: ncmpx
	real							:: dcmpx
	real							:: ocmpx
	integer							:: ncmpy
	real							:: dcmpy
	real							:: ocmpy
	integer							:: verb

	!! internal variables
	real,allocatable,dimension(:) 	:: zrefl ! reflector depths at
											 ! a given cmpx,cmpy
	real,allocatable,dimension(:) 	:: trace
	real,allocatable,dimension(:)	:: ztempx, ztempy
	real							:: z, cmpx, cmpy, dz2
	integer							:: icmpx, icmpy, iz, il 
	
	!! initialize SEP stuff
	call initpar()

	!! get number of layers
	call from_param("nl",nl)

	!! allocate space for input layer information
	allocate(v(nl),dipx(nl-1),dipy(nl-1),zref0(nl-1))

	!! get the other parameters and set their defaults
	call from_param("v",v)
	call from_param("dipx",dipx)
	call from_param("dipy",dipy)
	call from_param("zref0",zref0)
	call from_param("nz",nz)
	call from_param("dz",dz)
	call from_param("oz",oz)
	call from_param("ncmpx",ncmpx,1)
	call from_param("ocmpx",ocmpx,0.)
	call from_param("dcmpx",dcmpx,20.)
	call from_param("ncmpy",ncmpy,1)
	call from_param("ocmpy",ocmpy,0.)
	call from_param("dcmpy",dcmpy,20.)
	call from_param("verb",verb,0)

	!! write to history file
    call to_history("n1",nz)
    call to_history("n2",ncmpx)
    call to_history("n3",ncmpy)
    call to_history("o1",oz)
    call to_history("o2",ocmpx)
    call to_history("o3",ocmpy)
    call to_history("d1",dz)
    call to_history("d2",dcmpx)
    call to_history("d3",dcmpy)
    call to_history("label1","depth")
    call to_history("label2","cmpx")
    call to_history("label3","cmpy")
    call to_history("esize",4)

	if (verb) then
		write(0,*)'v=',v
		write(0,*)'zref0=',zref0
		write(0,*)'dipx=',dipx
		write(0,*)'dipy=',dipy
	end if

	! convert dip to radians
	dipx=dipx*pi/180.
	dipy=dipy*pi/180.

	!! allocate work space
	allocate(trace(nz),zrefl(nl-1),ztempx(nl-1),ztempy(nl-1))

	! set auxiliary variables
	ztempx=tan(dipx)
	ztempy=tan(dipy)
	dz2=dz/2.

	! deallocate space not needed anymore
	deallocate(dipx,dipy)

	! loop over cmpys
	do icmpy=1,ncmpy
		
		cmpy=ocmpy+(icmpy-1)*dcmpy

		! loop over cmpxs
		do icmpx=1,ncmpx

			cmpx=ocmpx+(icmpx-1)*dcmpx
			if (verb) then
				write(0,*)'icmpy=',icmpy,'icmpx=',icmpx
				write(0,*)'cmpy=',cmpy,'cmpx=',cmpx
			end if

			! create array of reflector depths at this location
			do il=1,nl-1
				zrefl(il)=zref0(il)+cmpy*ztempy(il)+cmpx*ztempx(il)
			end do

			! loop over depths
			do iz=1,nz
				z=(iz-1)*dz;

				! check what layer this depth belongs to and assign
				! the corresponding velocity
				if (z<=zrefl(1)) then
					trace(iz)=v(1)
				else if (z>zrefl(nl-1)) then
					trace(iz)=v(nl)
				else
					do il=2,nl-1
						if (z>zrefl(il-1).and.z<=zrefl(il)) then
							trace(iz)=v(il)
						end if
					end do
				end if
			end do

			! output this trace
			call srite("out",trace,4*nz)
		end do
	end do

	deallocate(v,zref0,ztempx,ztempy)
	deallocate(trace,zrefl)
	call exit(0)
end program make_vels
