# Time domain 15 degree (t',x',tau') migration of zero-offset section 
#
#Usage: <Headin T15mig n1=nt= d1=dt=  n2=nx=  d2=dx= vc=/vel=
#        trick=.14 taper=no stepdown=1 npad=1 >Headout
#%
#  the retarded coordinate transform:
#	t' 	= 	t + tau
#	x' 	= 	x
#	tau' 	= 	tau       ( tau = 2*z/v -- psudo-depth )
#
#n1=nt,n2=nx	trace length, number of traces
#d1=dt,d2=dx	sampling intervals
#vc=		constant velocity(v=vc); must specify either vc or vel
#vel=		headfile specifying velocity profile(s) of length nt
#		  when v is v(t) or v(t,x)
#trick=.14	1/6. trick in more accurate finite difference approximation 
#		  of second x derivative;  trick<0 when dipfilter needed
#taper=yes	taper side traces to reduce boundary reflections
#stepdown=1    	increment of downward step
#  
# author:       Zhiming Li, 8/17/84
#
#	Modified:	
#		stew        ---	transfered to convex and took out ap code 
#		li 12/6/85  ---	transpose data and velocity array	
#		stew 4/18/86 -- vectorized rewrite based on SEP-38 article
#		stew 4/25/86 -- convert interval vels to tau step internally
#		stew 5/6/86  -- hand roll rhs and triply outer loops for speed
#		stew 5/26/86 -- somewhat faster recoding
#		marta 6/17/86---changed rite to npad+1 instead of npad
#		stew&marta 9/30/86---fixed alp(ntau,nx) alloc and declaration
#				  ---changed min0(ittop,(it-1)/(stepdown-1))
#				     to min0(ittop,(it-1)/(stepdown-1)-1)
#	        marta 10/1/86--- fixed bug that left beta undefined
#               stew 7/27/89 --- For CM2 Fortran 8x
#
# ----------
#  Keyword:  15-degree time-domain migration
# ----------
integer nt,nx,ntp1,stepdown,npad
integer nxpow2,nxlog2,ntaupow2
real eps,dt,dx,beta,vc
character *10 taper
from either: integer n1:nt, n2:nx, stepdown=1
from either: real d1:dt, d2:dx
from par: real beta=.14, vc=0.0, eps=.0001
from par: character taper='no'
from par: integer npad=1
        npad=1
        nx=nx+npad*2
	ntp1=nt+1
	ntau=(nt-1)/stepdown
	nxpow2=1
	nxlog2=0
	while(nxpow2< nx) {
	   nxlog2=nxlog2+1
	   nxpow2=nxpow2*2
	}
	ntaupow2=1
	while(ntaupow2<ntau) {
	   ntaupow2=ntaupow2*2
	}

allocate:  real vfe(ntaupow2*nxpow2+nt-ntau)
allocate:  real pfe(ntp1,nxpow2)
allocate:  integer factor1(nxlog2+1),factor2(nxlog2+1)

subroutine parmig(eps,dt,dx,beta,nt,nx,vc,ntp1,taper,stepdown,ntau,npad,nxpow2,nxlog2,ntaupow2,pfe,vfe,factor1,factor2)
       integer nt,nx,ntp1,ix,it,itau,infd,outfd,velfd
       integer nxpow2,nxlog2,ntaupow2
       integer stepdown,ntau,ittop,min0,npad
       real eps
       real psec(ntp1,nxpow2)
       real p1(ntaupow2,nxpow2),p2(ntaupow2,nxpow2)
       real p3(ntaupow2,nxpow2),p4(ntaupow2,nxpow2)
       real alp(ntaupow2,nxpow2),a(ntaupow2,nxpow2)
       real b(ntaupow2,nxpow2),c(ntaupow2,nxpow2)
       real r1(ntaupow2,nxpow2),r2(ntaupow2,nxpow2),efac(nxpow2)
       real sc,beta,dt,dx,vc,factor,exp,float
       integer factor1(*),factor2(*)
%cmf$layout factor1(:serial)
%cmf$layout factor2(:serial)
       real vfe(ntaupow2,nxpow2),pfe(ntp1,nxpow2)
%cmf$layout pfe(:serial,:serial)
%cmf$layout vfe(:serial,:serial)
       integer reed,rite,icnt,itemp,iprev,inow,auxin,input,output
       integer n1vel,n2vel,lseek
       integer endvec(2),jt,irem
       integer iadr1,iadr2,iadalp
       integer loopcnt
       character*(*) taper
%       include '/usr/include/cm/paris-configuration-fort.h'
%       include '/usr/include/cm/CMF_defs.h'
#%       include 'T15mig.ifc'
define inext iprev

       infd=input()
       outfd=output()
 
       print *,'reading input section'
# read input section
       do ix=1,npad {
	  do it=1,nt {
	     pfe(it,ix)=0.0
       }  }
       do ix=nx-npad+1,nx {
	  do it=1,nt {
	     pfe(it,ix)=0.0
       }  }
       do ix=1+npad,nx-npad {
	   icnt=reed(infd,pfe(1,ix),nt*4)
	   pfe(ntp1,ix)=0.0
       }
       endvec(1)=ntp1
       endvec(2)=nx
       psec=0.0
%      call cmf_to_cm(psec,pfe,endvec)
 
# taper sides to reduce boundary artifacts

       if(taper(:1) == 'y') {
	    efac=0.0
	    do ix=1+npad,nx/32 {
		factor=float(ix-1-nx/32)/float(nx/32)
		efac(ix)=factor
		efac(nx+1-ix)=factor
		}
	    efac=exp(-efac**2)
	    p1=p1*spread(efac,1,ntp1)
	}


# read migration velocities
       print *,'setting up velocity field'
       sc = - (1./32.)*(dt/dx)**2
       alp = 0.0
       if(vc == 0.0) {
	   velfd = auxin('vel')
	   n1vel = nt
	   n2vel = nx
	   icnt= auxpar('n1','d',n1vel,'vel')
	   icnt= auxpar('n2','d',n2vel,'vel')
	   if(n1vel < nt) call erexit('trace and velocity profile length mismatch')
	   if(n2vel != 1 && n2vel < nx) call erexit('too few velocity profiles')
	   do ix=npad,nx-npad {
	       icnt=lseek(velfd,(ix-1)*n1vel*4,0) # goto beginning of vel. profile
	       icnt=reed(velfd,vfe(1,ix),nt*4) # read first nt velocities
	       if(icnt == 0) {  # assume v(tau) and propagate
		    do itau=1,ntau
			vfe(itau,ix) = vfe(itau,ix-1)
	       } else {         # process v(tau,x)
		   do it=1,nt
		       pfe(it,ix)=sc*pfe(it,ix)**2
		   if(stepdown > 1) call vpack(vfe(1,ix),nt,vfe(1,ix),ntau,stepdown)
	       }
	   }
           endvec(1)=ntau
           endvec(2)=nx
%          call cmf_to_cm(alp,vfe,endvec)
      } else {  # constant velocity
	   sc = sc * stepdown * vc**2
	   alp=sc
	   }
#
# set up padding: replicate first and last traces,velocities for zero slope
#

       p1(:,1:npad)=spread(p1(:,npad+1),2,npad)
       alp(:,1:npad)=spread(alp(:,npad+1),2,npad)
       p1(:,nx+1-npad:nx)=spread(p1(:,nx-npad),2,npad)
       alp(:,nx+1-npad:nx)=spread(alp(:,nx-npad),2,npad)
#
# set up and prefactor tridiagonal systems
#
       print *,'factoring tridiagonal system'
       a=beta+alp
       a(:,nx+1:nxpow2)=0.0
       b=1.0-2.0*a
       c=a
#
       a(:,nx)=-1.0
       b(:,nx)=+1.0
#
       b(:,1)=+1.0
       c(:,1)=-1.0
#
       call trifact(a,b,c,ntaupow2,nxpow2,factor1,factor2,nxlog2,eps,loopcnt)
       print *,nxpow2,ntaupow2,nxlog2,loopcnt
%       call cm_reset_timer(0)
%       call cm_start_timer(0)
#
# simplify rhs tridiagonal coefficients
#
       alp=beta-alp
		   
# initialize parallel subdiagonals
       print *,'starting migration'
       p1=0
       p2=0
       p3=0
       p4=0
       inow=1
       iprev=2
%       iadr1=CMF_get_field_id(r1)
%       iadr2=CMF_get_field_id(r2)
%       iadalp=CMF_get_field_id(alp)
#
# now for the "parallel" migration
#
#      P1 P4
#      P2 P3
#
#      (1+(beta+alpha)T) (P4+P2) = (1+(beta-alpha)T) (P1+P3)
#
       do it=nt,1+stepdown,-1 {  # work from bottom up
#
           ittop = ntp1-it
           if(stepdown>1)
	       ittop = min0(ittop,(it-1)/(stepdown-1)-1)
#
           
	   p1=cshift(p3,1,-1)
	   p1(1,:)=psec(it,:)
	   R1=P1+P3
#	   r2=r1-CSHIFT(r1,2,1)
%           call cm_f_news_sub_always_3_1L(iadr2,iadr1,iadr1,2,
%     1          cm_upward,23,8)
#	   r2=r2-CSHIFT(r2,2,-1)
%           call cm_f_news_sub_always_2_1L(iadr2,iadr2,2,
%     1          cm_downward,23,8)
#	   r2=r1-alp*r2
%           call cm_f_mult_subf_1L(iadr2,iadalp,iadr2,iadr1,23,8)
	   r2(:,1)=0.0
	   r2(:,nx)=0.0
 	   call factapp(r2,ntaupow2,nxpow2,factor1,factor2,nxlog2,loopcnt)
	   p4=r2-p2
	   if(stepdown == 1) {
	       if(it == 2) {
		   psec(it:nt,:)=p3(1:ntau,:)
		   }
	   } else {
	       itau=(it-2)/(stepdown-1)
	       irem=it-2-itau*(stepdown-1)
	       if(itau>0 && itau <= ntau) {
		    jt=it+itau-1
		    psec(jt,:)=p3(itau,:)
	       }
	       if(irem.eq.0) {
		    itau=itau-1
		    if(itau>0 && itau <= ntau) {
			jt=it+itau-1
			psec(jt,:)=p3(itau,:)
	            }
	       }
	   }
	   p3=p4
	   p2=p1
#
	   itemp=iprev
	   iprev=inow
	   inow=itemp
      }

# all done migrating; write output section

%      call cm_stop_timer(1)
      print *,'writing migrated section'
      endvec(1)=ntp1
      endvec(2)=nx
%      call cmf_from_cm(pfe,psec,endvec)
      do ix=npad+1,nx-npad {
	   icnt=rite(outfd,pfe(1,ix),nt*4)
      }
#
      print *,'done'
      return
      end

subroutine vpack(vin,nt,vout,ntau,stepdown)
integer nt,ntau,stepdown
real vin(stepdown,ntau),vout(ntau)
integer itau,it

%C$DIR SCALAR
do it=2,stepdown {
%C$DIR NO_RECURRENCE
    do itau=1,ntau {
	vin(1,itau) = vin(1,itau)+vin(it,itau)
    }
}

do itau=1,ntau
    vout(itau) = vin(1,itau)
return
end