# 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
#               stew 8/6/89 --- Pull tridiagonal code inline so I can use
#                               CMF$LAYOUT directives.
#		stew 8/11/89 --- Riemann sheeted input/output section
#                                to greatly reduce communication cost.
#		stew 9/19/89 --- added modpar() and reinitialization
#		stew 9/20/89 --- added frame buffer display of migration
#
# ----------
#  Keyword:  15-degree time-domain migration
# ----------
integer nt,nx,stepdown,npad
integer nxpow2,nxlog2,ntaupow2
real eps,dt,dx,beta,vc,clip
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, clip=5.0
from par: character taper='no'
from par: integer npad=1
        npad=1
        nx=nx+npad*2
        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(ntau*nx+nt-ntau)
allocate:  real pfe(nt,nx)
allocate:  integer factor1(nxlog2+1),factor2(nxlog2+1)

%      call cm_init
subroutine parmig(eps,dt,dx,beta,nt,nx,vc,taper,stepdown,ntau,npad,nxpow2,nxlog2,ntaupow2,pfe,vfe,factor1,factor2,clip)
       integer nt,nx,ix,it,itau,infd,outfd,velfd
       integer nxpow2,nxlog2,ntaupow2
       integer stepdown,ntau,ittop,min0,npad,i,j
       real eps,clip,ovc
#       real efac(nxpow2)
       real sc,beta,dt,dx,vc,factor,exp,float
       integer factor1(*),factor2(*)
       integer dispid
       real vfe(ntau,nx),pfe(nt,nx)
       integer itemp,input,output
       integer jt,irem
       integer iada,iadb,iadc
       integer iadislot,iadoslot,iadxedges
       integer iadp1,iadp2,iadp3,iadp4
       integer iadr1,iadr2,iadalp
       integer iadn1,iadn2
       integer iadpsec,iadpin,iadpout
       integer iadpsave, iadtemp
       integer iadnonzero
       integer iadnonzero1
       integer iadnonzero2
       integer iadnonzero3 
       integer iadplot
       integer loopcnt,nsheet
       integer lbound,p1addr,vp,geom
       integer modpar
       integer n3
       real tepsilon
       logical rpt,dodata
       character*(*) taper
       integer dscarr(7,2)
%       include '/usr/include/cm/paris-configuration-fort.h'

       infd=input()
       outfd=output()
%       nsheet=1+(nt-1)/ntaupow2
%       if(stepdown.eq.1) nsheet=1
%      dispid=CMFB_make_display_id()
%      call CMFB_attach_display('Mobil Video Room',dispid)
%      call CMFB_initialize_display(dispid,8,1)
%      call CMFB_overlay_disable(dispid)
%      call CMFB_set_colors_linear(dispid,0,127,
%     1          128,128,128,255,255,255)
%      call CMFB_set_colors_linear(dispid,128,255,
%     1          0,0,0,127,127,127)
%      i=max(1,1024/nxpow2)
%      j=max(1,1024/(nsheet*ntaupow2))
%      call CMFB_set_zoom(dispid,i-1,j-1,0)
%      call CMFB_clear(dispid,128,CMFB_green)
       n3=0
%      dscarr(1,1)=nxpow2
%      dscarr(2,1)=4
%      i=getch('xweight','i',dscarr(2,1))
%      dscarr(3,1)=cm_news_order
%      dscarr(4,1)=0
%      dscarr(6,1)=0
%      dscarr(1,2)=ntaupow2
%      dscarr(2,2)=2
%      i=getch('tweight','i',dscarr(2,2))
%      dscarr(3,2)=cm_news_order
%      dscarr(4,2)=0
%      dscarr(6,2)=0
%       geom=CM_create_detailed_geometry(dscarr,2)
%       vp=CM_allocate_vp_set(geom)
%       call CM_set_vp_set(vp)
%       iada=CM_allocate_stack_field(32)
%       iadb=CM_allocate_stack_field(32)
%       iadc=CM_allocate_stack_field(32)
%       iadp1=CM_allocate_stack_field(32)
%       iadp2=CM_allocate_stack_field(32)
%       iadp3=CM_allocate_stack_field(32)
%       iadp4=CM_allocate_stack_field(32)
%       iadr1=CM_allocate_stack_field(32)
%       iadn1=CM_allocate_stack_field(32)
%       iadr2=CM_allocate_stack_field(32)
%       iadn2=CM_allocate_stack_field(32)
%       iadalp=CM_allocate_stack_field(32)
%       iadpsec=CM_allocate_stack_field(32*nsheet)
%       iadpsave=CM_allocate_stack_field(32*nsheet)
%       iadplot=CM_allocate_stack_field(8)
%       iadislot=CM_allocate_stack_field(1)
%       iadoslot=CM_allocate_stack_field(1)
%       iadxedges=CM_allocate_stack_field(1)
%       iadnonzero=CM_allocate_stack_field(1)
%       iadnonzero1=CM_allocate_stack_field(1)
%       iadnonzero2=CM_allocate_stack_field(1)
%       iadnonzero3=CM_allocate_stack_field(1)
#
#     xedges=.false.
#     xedges(:,1)=.true.
#     xedges(:,nx)=.true.
#
%       call CM_set_context
%       call CM_my_news_coordinate_1L(iadn2,0,32)
%       call CM_u_eq_zero_1L(iadn2,32)
%       call CM_store_test_always(iadxedges)
%       call CM_u_eq_constant_1L(iadn2,nx-1,32)
%       call CM_logior_test_always(iadxedges)
%       call CM_store_test_always(iadxedges)
%       call CM_u_ge_constant_1L(iadn2,nx,32)
%       call CM_store_test_always(iadnonzero1)
%       call CM_u_lt_constant_1L(iadn2,nxpow2,32)
%       call CM_logior_test_always(iadnonzero1)
%       call CM_store_test_always(iadnonzero1)
%       call CM_u_eq_constant_1L(iadn2,nx-1,32)
%       call CM_store_test_always(iadnonzero2)
%       call CM_u_eq_zero_1L(iadn2,32)
%       call CM_store_test_always(iadnonzero3)
%       call CM_my_news_coordinate_1L(iadn1,1,32)
%       ix=CM_set_system_leds_mode(CM_leds_off)
#
       print *,'reading input section'
# read input section
       iadpin=iadpsave
       dodata = .true.
%      call CM_set_context
       call getsec(infd,pfe,npad,nt,nx,iadpin,ntaupow2,dodata)
       dodata = .false.
# interactive loop
       while(0.eq.modpar()) {
       n3=n3+1
to history: integer n3
#
# set input and output slot pointers
#
#      islot=.false.
#      islot(1+mod(nt,ntaupow2),:)=.true.
#      if(stepdown.gt.1){
#         oslot=islot
#             } else {
#          oslot=.true.
#          oslot(1,:)=.false.
#     }
%       call CM_u_eq_constant_1L(iadn1,mod(nt,ntaupow2),32)
%       call CM_store_test_always(iadislot)
      if(stepdown.gt.1) {
%       call CM_u_move_always_1L(iadoslot,iadislot,1)
      } else {
%       call CM_u_ne_zero_1L(iadn1,32)
%       call CM_store_test_always(iadoslot)
      }
# copy input section to computational input/output grid
%      call CM_set_context
       iadtemp=iadpsave
       iadpout=iadpsec
       do i=1,nsheet {
%          call cm_f_move_always_1L(iadpout,iadtemp,23,8)
%          iadtemp=cm_add_offset_to_field_id(iadtemp,32)
%          iadpout=cm_add_offset_to_field_id(iadpout,32)
	}
       iadpin=iadpsec
       dodata=.false.
       call getsec(infd,pfe,npad,nt,nx,iadpin,ntaupow2,dodata)
       iadpout=iadpin
#
# STILL NEED TO CONVERT THIS
#
# 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)
#	    r1=spread(efac,1,ntaupow2)
#           iadpout=iadpsec
#	    do it=1,nt,ntaupow2 {
#%      if(it.gt.1)iadpout=CM_add_offset_to_field_id(iadpout,32)
#%      call CM_f_multiply_always_2_1L(iadpout,iadr1,23,8)
#	    }
#        }


# read migration velocities
       print *,'setting up velocity field'
from par: real vc=vc
from par: real beta=beta
from par: real eps=eps
from par: real clip=clip
       ovc=127.0/clip
       sc = - (1./32.)*(dt/dx)**2
#       alp = 0.0
%      call CM_f_move_const_always_1L(iadalp,0.0,23,8)
       if(vc == 0.0) {
	   call getvel(nt,nx,npad,ntau,vfe,stepdown,sc,iadalp)
      } else {  # constant velocity
           sc = sc * stepdown * vc**2
#           alp=sc
%      call CM_f_move_const_always_1L(iadalp,sc,23,8)
           }
#
# 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
#       c(:,1)=-1.0
#       a(:,nx)=-1.0
#
%       call CM_f_add_const_always_3_1L(iada,iadalp,beta,23,8)
%       call CM_load_context(iadnonzero1)
%       call CM_f_move_zero_1L(iada,23,8)
%       call CM_set_context
%       call CM_f_mult_const_add_const_1L(iadb,iada,-2.0,1.0,23,8)
%       call CM_f_move_always_1L(iadc,iada,23,8)
%       call CM_load_context(iadnonzero2)
%       call CM_f_move_constant_1L(iada,-1.0,23,8)
%       call CM_load_context(iadnonzero3)
%       call CM_f_move_constant_1L(iadc,-1.0,23,8)
#
#       where(xedges)b=1.0
%       call CM_load_context(iadxedges)
%       call CM_f_move_constant_1L(iadb,1.0,23,8)
%       call CM_set_context
#
       call trifact(iada,iadb,iadc,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
%       call CM_set_context
%       call CM_f_subfrom_const_always_2_1L(iadalp,beta,23,8)
# initialize parallel subdiagonals
       print *,'starting migration nt=',nt,' nx=',nx,' ntau=',ntau,' ntaupow2=',ntaupow2
#       p1=0
#       p2=0
#       p3=0
#       p4=0
%      call CM_f_move_zero_always_1L(iadp1,23,8)
%      call CM_f_move_zero_always_1L(iadp2,23,8)
%      call CM_f_move_zero_always_1L(iadp3,23,8)
%      call CM_f_move_zero_always_1L(iadp4,23,8)
#
# now for the "parallel" migration
#
#      P1 P4
#      P2 P3
#
#      (1+(beta+alpha)T) (P4+P2) = (1+(beta-alpha)T) (P1+P3)
#
# At top of loop, P1, P2, and P3 are in correct position
#    P1 runs from it to it+ntaupow2-1
#    P2 runs from it+1 to it+ntaupow2
#    P3 runs from it+1 to it+ntaupow2
# Logic:
#    shift input pointer, slot psec(it) into p1(it)
#    compute P4
#    P4 runs from it to it+ntaupow2-1
#    shift output pointer, slot p4(jt) into psec(jt) as needed
#    shift P1 to P2
#    copy  P4 to P1
#    shift P1 to P3
#
       do it=nt,1+stepdown,-1 {  # work from bottom up
#
#
#         islot=cshift(islot,1,1)
%       call CM_get_from_news_always_1L(iadislot,iadislot,1,
%     1      CM_upward,1)
	  if(mod(it,ntaupow2).eq.0) {
%       iadpin=CM_add_offset_to_field_id(iadpin,-32)
         }
%       call CM_load_context(iadislot)
%       call CM_f_move_1L(iadp1,iadpin,23,8)
%       call CM_set_context
#         R1=P1+P3
%       call CM_f_add_always_3_1L(iadr1,iadp1,iadp3,23,8)
#         r2=r1-CSHIFT(r1,2,1)
%           call CM_f_news_sub_always_3_1L(iadr2,iadr1,iadr1,0,
%     1          CM_upward,23,8)
#         r2=r2-CSHIFT(r2,2,-1)
%           call cm_f_news_sub_always_2_1L(iadr2,iadr2,0,
%     1          cm_downward,23,8)
#         r2=r1-alp*r2
%           call cm_f_mult_subf_1L(iadr2,iadalp,iadr2,iadr1,23,8)
#         where(xedges) r2=0.0
%       call CM_load_context(iadxedges)
%       call CM_f_move_zero_1L(iadr2,23,8)
%       call CM_set_context
#
#####    call factapp(r2,ntaupow2,nxpow2,factor1,factor2,nxlog2,loopcnt)
          lbound=1
           do i=1,loopcnt {
#         p3=r2
%       call CM_f_move_always_1L(iadp3,iadr2,23,8)
#          r1=cshift(p3,dim=2,shift=-lbound)
%       call CM_get_from_power_two_always_1L(iadr1,iadp3,0,i-1,
%     1      CM_downward,32)
%          call CM_f_mult_add_1L(iadr2,factor1(i),iadr1,iadr2,23,8)
#           r1=cshift(p3,dim=2,shift=lbound)
%       call CM_get_from_power_two_always_1L(iadr1,iadp3,0,i-1,
%     1      CM_upward,32)
%          call CM_f_mult_add_1L(iadr2,factor2(i),iadr1,iadr2,23,8)
          lbound=lbound*2
%       call CM_f_ne_zero_1L(iadr2,23,8)
%       call CM_store_test(iadnonzero)
%       call CM_latch_leds_always(iadnonzero)
          }
%          call cm_f_multiply_2_1L(iadr2,factor1(loopcnt+1),23,8)
#          p4=r2-p2
%       call CM_f_mult_sub_1L(iadp4,iadr2,factor1(loopcnt+1),
%     1      iadp2,23,8)
#           if(it>325 && it<365) {
#%         call CM_f_move_always_1L(iada,iadp4,23,8)
#	       call putcomp(nxpow2,ntaupow2,iada,pfe)
#	       }
           if(stepdown == 1) {
               if(it == 2) {
%       call CM_load_context(iadoslot)
%       call CM_f_move_1L(iadpsec,iadp4,23,8)
                   }
           } else {
               itau=(it-2)/(stepdown-1)
               irem=it-2-itau*(stepdown-1)
	       jt=it+itau-1
               if(itau>0 && itau <= ntau && jt <= nt) {
#      print *,' it=',it,' itau=',itau,' irem=',irem,' jt=',jt
#		    oslot=cshift(oslot,1,1)
%       call CM_get_from_news_always_1L(iadoslot,iadoslot,1,
%     1      CM_upward,1)
                    if(mod(jt,ntaupow2).eq.0){
%       iadpout=CM_add_offset_to_field_id(iadpout,-32)
          }
%       call CM_load_context(iadoslot)
%       call CM_f_move_1L(iadpout,iadp4,23,8)
               }
               if(irem.eq.0) {
                    itau=itau-1
		    jt=jt-1
                    if(itau>0 && itau <= ntau && jt <= nt) {
#      print *,' it=',it,' itau=',itau,' irem=',irem,' jt=',jt
#		    oslot=cshift(oslot,1,1)
%       call CM_get_from_news_always_1L(iadoslot,iadoslot,1,
%     1      CM_upward,1)
                    if(mod(jt,ntaupow2).eq.0){
%       iadpout=CM_add_offset_to_field_id(iadpout,-32)
          }
%       call CM_load_context(iadoslot)
%       call CM_f_move_1L(iadpout,iadp4,23,8)
                    }
               }
           }
%       call CM_set_context
#          p2=cshift(p1,1,1)
%       call CM_get_from_news_always_1L(iadp2,iadp1,1,
%     1      CM_upward,32)
#	   p1=p4
%       call CM_f_move_always_1L(iadp1,iadp4,23,8)
#	   p3=cshift(p1,1,1)
%       call CM_get_from_news_always_1L(iadp3,iadp1,1,
%     1      CM_upward,32)
#
      }

# all done migrating; write output section

%      call cm_stop_timer(1)
      print *,'writing migrated section'
       iadpout=iadpsec
       call putsec(outfd,nt,nx,npad,ntaupow2,iadpout,pfe)
# deallocate tridiagonal factors
       do i=0,loopcnt {
%       call cm_deallocate_heap_field(factor1(i+1))
	}
# put migrated image on display
       iadpin=iadpsec
       do i=1,nsheet {
%      call cm_f_max_constant_3_1L(iadr1,iadpin,-clip,23,8)
%      call cm_f_min_constant_2_1L(iadr1,clip,23,8)
%      call cm_f_multiply_const_always_2_1L(iadr1,ovc,23,8)
%      call cm_s_f_floor_2_2L(iadplot,iadr1,8,23,8)
%      call CMFB_write_always(dispid,CMFB_green,iadplot,
%     1     0,(i-1)*ntaupow2)
%      iadpin=cm_add_offset_to_field_id(iadpin,32)
	}
       }
#
      print *,'done'
      return
      end