!--------------------------------------------------------------------
! DESCRPTION:
!	program to apply a top mute to a set of 2D CMP gathers.
! 	The mute is applied as a function of offset and time and
!	is intended to remove NMO-stretch.
!--------------------------------------------------------------------	
!	Input:
!		picks_file : name of ASCII file with the mute picks information.
!			The file contains three columns: cmp, time and offset.
!			The cmp value is repeated as many times as time-offset
!			picks are given for that cmp. The program interpolates
!			linearly between offsets in a given cmp and linearly between
!			cmps.
!		verb : =1 for verbose output
!
!	EXAMPLES:
!		cmp_mute.x < cmps_nmo.H picks_file="mute_pick" > cmps_nmo_mute.H
!---------------------------------------------------------------------
program cmp_mute

	use sep
	use constants

	implicit none

	! program arguments
	character(len=20)				:: picks_file
	logical							:: mask
	integer							:: verb

	!! internal variables
	real,allocatable,dimension(:) 	:: trace,temp
	real,allocatable,dimension(:) 	:: mute_cmps
	real,allocatable,dimension(:,:)	:: mute_times,mute_offs
	real,allocatable,dimension(:,:)	:: cmp_offs_tmute
	integer							:: i,icmps,ihx,ihh,imx,it,ihx1,icmp1
	integer							:: ioff,j 
	integer							:: nt,nhx,nmx,esize
	real							:: o1,ohx,omx,dhx,dmx,dt
	real							:: cmp,off,time
	real							:: old_cmp,old_off,old_time
	real							:: hx1,hxn,hh1,hh2,hhn
	real							:: a,b,t1,t2,tmin,t
	real							:: cmp0,cmp1,cmp2,cmpn
	real							:: hx,mx
	integer							:: nmute_cmps,nmute_off,nmute_noffmax
	integer							:: infp,stat
	logical							:: first_rec
	character(len=6)				:: label1,label2,label3
	
	!! initialize SEP stuff
	call initpar()

	!! read from history file
    call from_history("n1",nt)
    call from_history("n2",nhx)
    call from_history("n3",nmx)
    call from_history("o1",o1)
    call from_history("o2",ohx)
    call from_history("o3",omx)
    call from_history("d1",dt)
    call from_history("d2",dhx)
    call from_history("d3",dmx)
    call from_history("label1",label1,"time")
    call from_history("label2",label2,"offset")
    call from_history("label3",label3,"cmp")
    call from_history("esize",esize)

	!! write to history file
    call to_history("n1",nt)
    call to_history("n2",nhx)
    call to_history("n3",nmx)
    call to_history("o1",o1)
    call to_history("o2",ohx)
    call to_history("o3",omx)
    call to_history("d1",dt)
    call to_history("d2",dhx)
    call to_history("d3",dmx)
    call to_history("label1",label1)
    call to_history("label2",label2)
    call to_history("label3",label3)
    call to_history("esize",esize)

	!! get name of file with the mute picks
	call from_param("picks_file",picks_file)
	call from_param("verb",verb,0)
	call from_param("mask",mask,.false.)

	!! open file to read in the mute information
	infp=20
	open(unit=infp,file=picks_file,status='old',action='read',iostat=stat)
	if (stat/=0) then
		write(0,*)'ERROR: could not open file',picks_file,'for input'
		stop
	end if

	!! count cmps at which picks are given and the largest 
	!! number of time-offset pairs
	nmute_cmps=1	!! number of cmps at which mute pairs were given
	nmute_off=0		!! number of offsets at which mute pairs were given
	nmute_noffmax=0	!! maximum number of offset-time pairs in a cmp
	first_rec=.true.
	i=0
	do
		!read(infp,'(F10.2)',iostat=stat) cmp
		read(infp,*,iostat=stat) cmp
		if (stat<0) exit			!! found end of file
		i=i+1
		write(0,*)'i=',i,'cmp=',cmp
		if (first_rec .eq. .true.) then				!! first file record
			old_cmp=cmp
			nmute_off=nmute_off+1
			first_rec=.false.
		else						!! other records
			if (cmp==old_cmp) then	!! same cmp	
				nmute_off=nmute_off+1
			else					!! new cmp
				nmute_cmps=nmute_cmps+1		!! increase cmp count
				old_cmp=cmp			!! reset cmp
				if(nmute_off>nmute_noffmax) nmute_noffmax=nmute_off
				nmute_off=0			!! reset mute offset count
			end if	
		end if
	end do
	if (i==0) then
		write(0,*)'ERROR: NOTHING WAS READ FROM INPUT FILE'
		stop
	else
		write(0,*)'READ',i,'RECORDS','nmute_cmps=',nmute_cmps,'nmute_noffmax=',nmute_noffmax
	end if
	rewind(infp)
			
	!! allocate working space
	allocate(mute_cmps(nmute_cmps))
	allocate(mute_times(nmute_cmps,nmute_noffmax))
	allocate(mute_offs(nmute_cmps,nmute_noffmax))
	allocate(cmp_offs_tmute(nmute_cmps,nhx))

	!! read in the mute information
	first_rec=.true.
	ioff=0
	icmps=1
	i=0
	do
		!read(infp,'(3(F10.2))',iostat=stat) cmp,time,off
		read(infp,*,iostat=stat) cmp,time,off
		if (stat < 0) exit
		i=i+1
		!write(0,*)'i=',i,'ioff=',ioff,'cmp=',cmp,'off=',off,'time=',time
	
		!! check to see if this is the first mute pair read
		if (first_rec .eq. .true.) then
			ioff=ioff+1
			mute_cmps(icmps)=cmp
			mute_offs(icmps,ioff)=off
			mute_times(icmps,ioff)=time
			old_cmp=cmp
			old_off=off
			old_time=time
			first_rec=.false.
!write(0,*)'i=',i,'icmps=',icmps,'cmp=',cmp,'mute_cmp=',mute_cmps(icmps)
		else
	
			!! check to see if this is still the same cmp
			if (cmp==old_cmp) then
				ioff=ioff+1
				mute_offs(icmps,ioff)=off
				mute_times(icmps,ioff)=time
				old_off=off
				old_time=time
			else
				!! found new cmp. 
				old_cmp=cmp

				!! make sure all cmps have the same number of mute
				!! time-offset pairs
				do j=ioff,nmute_noffmax
					mute_offs(icmps,j)=old_off
		            mute_times(icmps,j)=old_time
				end do

				!! save this mute time-offset pair
				icmps=icmps+1
				ioff=1
				mute_cmps(icmps)=cmp
				mute_offs(icmps,ioff)=off
				mute_times(icmps,ioff)=time
				old_off=off
				old_time=time
!write(0,*)'i=',i,'icmps=',icmps,'cmp=',cmp,'mute_cmp=',mute_cmps(icmps)
			end	if
		end if	
	end do	
	close(infp)
		
!!do icmps=1,nmute_cmps
!!	write(0,*)'cmp=',mute_cmps(icmps)
!!	do ioff=1,nmute_noffmax
!!		write(0,*)'time=',mute_times(icmps,ioff),'off=',mute_offs(icmps,ioff) 
!!	end do
!!end do

	!! interpolate the mute offset-time pairs for all offsets
	!! at the cmps at which they were given
	do icmps=1,nmute_cmps
		do ihx=1,nhx
			hx=ohx+(ihx-1)*dhx
			hx1=mute_offs(icmps,1)
			hxn=mute_offs(icmps,nmute_noffmax)
			if (hx<hx1) then 
				cmp_offs_tmute(icmps,ihx)=0.
			else if (hx>=hxn) then
				cmp_offs_tmute(icmps,ihx)=mute_times(icmps,nmute_noffmax)
			else
				!! do linear interpolation
				do ihh=1,nmute_noffmax-1
					hh1=mute_offs(icmps,ihh)
					hh2=mute_offs(icmps,ihh+1)
					if (hx==hh1) then
						cmp_offs_tmute(icmps,ihx)=mute_times(icmps,ihh)
						exit
					else 
						if (hx>hh1.AND.hx<hh2) then
							a=(hx-hh1)/(hh2-hh1)
							b=1-a
							t1=mute_times(icmps,ihh)
							t2=mute_times(icmps,ihh+1)
							cmp_offs_tmute(icmps,ihx)=b*t1+a*t2
						end if
					end if
				end do
			end if
		end do
	end do
				
!!do icmps=1,nmute_cmps
!!	write(0,*)'cmp=',mute_cmps(icmps)
!!	do ihx=1,nhx
!!		hx=ohx+(ihx-1)*dhx
!!		write(0,*)'hx=',hx,'cmp_offs_tmute=',cmp_offs_tmute(icmps,ihx)
!!	end do
!!end do
	!! allocate space for the traces
	deallocate(mute_offs,mute_times)
	allocate(trace(nt),temp(nt))
	
	!! interpolate mute to all cmps
	do imx=1,nmx

		mx=omx+(imx-1)*dmx

		!! compute interpolation distances for this cmp with
		!! respect to the two cmps that bracket it and for which
		!! mute time-offset pairs were actually supplied 
		!
		do icmps=1,nmute_cmps-1
			cmp0=mute_cmps(1)
			cmp1=mute_cmps(icmps)
			cmp2=mute_cmps(icmps+1)
			cmpn=mute_cmps(nmute_cmps)
			if (mx<=cmp0) then 	!! extrapolate to the left
				icmp1=-1
				exit
			else if	(mx>=cmpn) then		!! extrapolate to the right
				icmp1=-2
				exit
			else	
				!! interpolate
				if (mx>cmp1 .and. mx<cmp2) then
					a=(mx-cmp1)/(cmp2-cmp1)
					b=1-a
				end if
				icmp1=icmps
			end if	
		end do

	
		!! loop over offsets to apply the mute
		do ihx=1,nhx

			hx=ohx+(ihx-1)*dhx

			!! read trace
			call sreed("in",trace,4*nt)

			!! compute the mute time for this offset according
			!! to its interpolation distances to the corresponding
			!! offset trace in the two cmps that bracket this cmp
			!! and for which mute time-offset pairs were given
			!
			if (icmp1==-1) then 	!! extrapolate to the left
				tmin=cmp_offs_tmute(1,ihx)
			else if (icmp1==-2) then	!! extrapolate to the right
				tmin=cmp_offs_tmute(nmute_cmps,ihx)
			else
				t1=cmp_offs_tmute(icmp1,ihx)
				t2=cmp_offs_tmute(icmp1+1,ihx)
				tmin=b*t1+a*t2
			end if

if (mx==2300) then
	write(0,*)'hx=',hx,'tmin=',tmin
	if (hx==1045) then
		write(0,*)'o1=',o1,'dt=',dt
	end if
end if
			!! mute the input trace
			do it=1,nt
				t=o1+(it-1)*dt
	
				if (t<tmin) then
					temp(it)=0.
				else
					if (mask) then
						temp(it)=1.
					else 
						temp(it)=trace(it)
					end if
				end if
			end do

			!! output muted trace
			call srite("out",temp,4*nt)
		end do
	end do

	!! deallocate work space
	deallocate(trace,temp)
	deallocate(cmp_offs_tmute,mute_cmps)

end program cmp_mute
