!APPLY A COSINE TAPER
module taper_mod
  use down_mod
  use kmap_mod

  implicit none
  real,private,save,allocatable,dimension(:,:) ::cos_taper_filtery,cos_taper_filterx
  !! space domain taper for U, D
  real,dimension(:,:),private,save,allocatable::taper
  real,dimension(:,:,:),private,save,allocatable::top,bott,right,left
  integer:: x1,x2,y1,y2

  contains
!----------------------------------------------------------------
!! cosine^2 taper for the wavefields
  subroutine xytaper_init()
    real:: pcentdist,a,t
    integer::er,maxd,dist
#ifdef DEBUG
    character(len=128), parameter :: unit='xytaper_init'
    if(run%depar) call in(unit)
#endif

    !! default is 1 trace zeroed, 7 traces damped at each side 
    allocate(taper(down%amx%n,down%amy%n),stat=er)

    if(er/=0) call seperr('cant allocate taper array (image.F90)')
    call from_param('taper_pcent',pcentdist,100.*8./real(down%amx%n))

    if(pcentdist==0.)  then
       down%taperme=.false.
       write(run%io,*)'        *       (no taper applied)'
       return
    end if
    pcentdist=pcentdist/100.
    t=pcentdist
    er=nint(pcentdist*down%amx%n)

    call from_param("minx1",down%minx1,er)
    call from_param("maxx1",down%maxx1,down%minx1)

    if (down%amy%n.gt.1) then
       pcentdist=max(t,8./real(down%amy%n))
       er=nint(pcentdist*down%amy%n)
       call from_param("miny1",down%miny1,er)
       call from_param("maxy1",down%maxy1,down%miny1)
    end if

    call init_cos_taper()
    call taper_report() 

    a=4*size(taper)
    down%mem=down%mem+a
#ifdef MEMCHECK
       if(run%debug) write(run%io,*) 'Taper.F90 taper: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
#ifdef DEBUG
    if(run%depar) call out(unit)
#endif
     end subroutine xytaper_init
!----------------------------------------------------------------
! cos taper function (cos^2 by changing 4 lines after 118, and 2 after 187)
                                             !! 0.                     (n-1)*d  
                                             !! |                         |
     subroutine init_cos_taper()              !!_______/----------\_________
       real,allocatable::cos_taper_filter(:)  !! 0n   1n          1x       0x
       integer :: i1_x,i0_x,i1_n,i0_n         !!1x,1n are pole location indices
       integer :: i1,i2,i3,n1,ix,iy
       real    :: pi2,d,d1,min1,max1,a
#ifdef DEBUG
       character(len=128), parameter :: unit='init_cos_taper'
       if(run%depar) call in(unit)
#endif
       pi2=atan(1.)*2. !pi/2.

       n1=down%amx%n   

       allocate(cos_taper_filter(n1))
       a=4*size(cos_taper_filter)
       down%mem=down%mem+a
#ifdef MEMCHECK
       if(run%debug) write(run%io,*) 'Taper.F90 tmptaper: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
       i0_n=max(1,nint(0.1*real(down%minx1)))
       i1_n=max(i0_n+1,down%minx1)
       i1_x=min(n1-1,max(i1_n+1,n1-down%maxx1+1))
       i0_x=max(i1_x+1,min(n1,1+n1-nint(0.1*down%maxx1)))

       !! indices to stop and begin tapering wflds
       !! really dont need to allocate all of taper if only using such a small part
       x1=i1_n-1;x2=i1_x+1

       cos_taper_filter(1:i0_n)=0.
       cos_taper_filter(i0_x:)=0.
       cos_taper_filter(i1_n:i1_x)=1.
       if(down%amy%n.gt.1) then  !! tapery kills all xline if xline<5
          do i1=i0_n,i1_n 
             !            cos_taper_filter(i1)=(1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) )
             cos_taper_filter(i1)=sqrt(1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) )
          end do
          do i1=i1_x,i0_x 
             !            cos_taper_filter(i1)=(1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) )
             cos_taper_filter(i1)=sqrt(1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) )
          end do
       else
          do i1=i0_n,i1_n 
             !            cos_taper_filter(i1)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))**2
             cos_taper_filter(i1)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))
          end do
          do i1=i1_x,i0_x 
             !            cos_taper_filter(i1)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))**2
             cos_taper_filter(i1)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))
          end do

       end if
       do i2=1,down%amy%n
          do i1=1,down%amx%n
             taper(i1,i2)=cos_taper_filter(i1)
          end do
       end do

       !----------------------------------------------------------------
       !!Y direction
       if (down%amy%n.gt.1) then !! tapery kills all xline if xline<5
          a=4*size(cos_taper_filter)
          down%mem=down%mem-a
#ifdef MEMCHECK
          if(run%debug) write(run%io,*) 'Taper.F90 tmptaper: deallocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
          deallocate(cos_taper_filter)
          n1=down%amy%n
          allocate(cos_taper_filter(n1))
          a=4*size(cos_taper_filter)
          down%mem=down%mem+a
#ifdef MEMCHECK
          if(run%debug) write(run%io,*) 'Taper.F90 tmptaper: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
          i0_n=max(1,nint(0.1*real(down%miny1)))
          i1_n=max(i0_n+1,down%miny1)
          i1_x=min(n1-1,max(i1_n+1,n1-down%maxy1+1))
          i0_x=max(i1_x+1,min(n1,1+n1-nint(0.1*down%maxy1)))

          y1=i1_n-1;y2=i1_x+1
          cos_taper_filter(1:i0_n)=0.
          cos_taper_filter(i0_x:)=0.
          cos_taper_filter(i1_n:i1_x)=1.
          do i1=i0_n,i1_n 
             !            cos_taper_filter(i1)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))
             cos_taper_filter(i1)=sqrt((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))
          end do
          do i1=i1_x,i0_x 
             !            cos_taper_filter(i1)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))
             cos_taper_filter(i1)=sqrt((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))
          end do
         
          do i2=1,down%amy%n
             do i1=1,down%amx%n
                taper(i1,i2)=cos_taper_filter(i2)*taper(i1,i2)
             end do
          end do

          allocate(left(down%amx%n,y1,run%o_nth),right(down%amx%n,down%amy%n-y2+1,run%o_nth))
          allocate(top(x1,y2-y1-1,run%o_nth),bott(down%amx%n-x2+1,y2-y1-1,run%o_nth))

          do i3=1,run%o_nth
             do i2=1,y1
                do i1=1,down%amx%n
                   left(i1,i2,i3)=taper(i1,i2)
                end do
             end do
             
             iy=0
             do i2=y2,down%amy%n
                iy=iy+1
                do i1=1,down%amx%n
                   right(i1,iy,i3)=taper(i1,i2)
                end do
             end do
          end do
       else
          y1=0;y2=2
          allocate(top(x1,1,run%o_nth),bott(down%amx%n-x2+1,1,run%o_nth))
       end if

       !! top stripe
       do i3=1,run%o_nth
          iy=0
          do i2=y1+1,y2-1
             iy=iy+1
             do i1=1,x1
                top(i1,iy,i3)=taper(i1,i2)
             end do
          end do
          !!bottom stripe
          iy=0
          do i2=y1+1,y2-1
             iy=iy+1
             ix=0
             do i1=x2,down%amx%n
                ix=ix+1
                bott(ix,iy,i3)=taper(i1,i2)
             end do
          end do
       end do

       a=4*(size(left)+size(right)+size(bott)+size(top))
       down%mem=down%mem+a
#ifdef MEMCHECK
       if(run%debug) write(run%io,*) 'Taper.F90 liltapers: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
       a=4*(size(taper)+size(cos_taper_filter))
       down%mem=down%mem-a
#ifdef MEMCHECK
       if(run%debug) write(run%io,*) 'Taper.F90 tmptaper: deallocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
       deallocate(cos_taper_filter,taper)
#ifdef DEBUG
       if(run%depar) call out(unit)
#endif

  end subroutine init_cos_taper
!----------------------------------------------------------------
  subroutine taperwfld(X,iw,iith)
    !    complex,dimension(:,:,:),allocatable::wfld
    type(dat)::X
    integer,intent(in):: iw,iith
    integer::i1,i2,i3,iy,ix,iws,ith=1
#ifdef DEBUG
    character(len=128), parameter :: unit='taperwfld'
    if(run%depar) call in(unit)
#endif

    if(trim(X%tag).eq."D_d".and..not.down%passive) then
       if (down%amx%n.gt.1) then
          if(X%indx.lt.size(top,1))then
             write(run%io,*)"Source X location within small X taper zone."
             write(run%io,*)" All zero image for current shot"
             down%validshot=.false.
             if(down%shotnum.lt.1) call seperr("Migration terminated due to bad padding.")
          end if

          if(X%indx.gt.down%amx%n-size(bott,1)) then
             write(run%io,*)"Source X location within big X taper zone."
             write(run%io,*)" All zero image for current shot"
             down%validshot=.false.
             if(down%shotnum.lt.1) call seperr("Migration terminated due to bad padding.")
          end if
       end if

       if (down%amy%n.gt.1) then
          if(X%indy.lt.size(left,2)) then
             write(run%io,*)"Source Y location within small Y taper zone."
             write(run%io,*)" All zero image for current shot"
             down%validshot=.false.
             if(down%shotnum.lt.1) call seperr("Migration terminated due to bad padding.")
          end if
          if(X%indy.gt.down%amy%n-size(right,2))then
             write(run%io,*)"Source Y location within big Y taper zone."
             write(run%io,*)" All zero image for current shot"
             down%validshot=.false.
             if(down%shotnum.lt.1) call seperr("Migration terminated due to bad padding.")
          end if
       end if
    end if
    if(.not.down%validshot) return

#ifdef SEP_OMP
!$OMP PARALLEL PRIVATE(ith)
    ith=1+omp_get_thread_num()
!$OMP DO PRIVATE(iws,ix,iy,i1,i2)
#endif
    do iws=1,down%nws
       if(down%amy%n.gt.1) then
          !! left stripe
          do i2=1,y1
             do i1=1,down%amx%n
                !             wfld(i1,i2,iws)=wfld(i1,i2,iws)*left(i1,i2,ith)
                X%bin(i1,i2,iws)=X%bin(i1,i2,iws)*left(i1,i2,ith)
             end do
          end do
          !! right stripe
          iy=0
          do i2=y2,down%amy%n
             iy=iy+1
             do i1=1,down%amx%n
                X%bin(i1,i2,iws)=X%bin(i1,i2,iws)*right(i1,iy,ith)
                !             wfld(i1,i2,iws)=wfld(i1,i2,iws)*right(i1,iy,ith)
             end do
          end do
       else
          y1=0;y2=2
       end if
       !! top stripe
       iy=0
       do i2=y1+1,y2-1
          iy=iy+1
          do i1=1,x1
             X%bin(i1,i2,iws)=X%bin(i1,i2,iws)*top(i1,iy,ith)
          end do
       end do
       !!bottom stripe
       iy=0
       do i2=y1+1,y2-1
          iy=iy+1
          ix=0
          do i1=x2,down%amx%n
             ix=ix+1
             X%bin(i1,i2,iws)=X%bin(i1,i2,iws)*bott(ix,iy,ith)
          end do
       end do
    end do
#ifdef SEP_OMP
!$OMP END DO
!$OMP END PARALLEL 
#endif
#ifdef DEBUG
    if(run%depar) call out(unit)
#endif

  end subroutine taperwfld
!----------------------------------------------------------------
  subroutine cos_taper_cmplx_allocate()
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='cos_taper_cmplx_allocate'
    if(run%depar) call in(unit)
#endif

    allocate(cos_taper_filterx(down%amx%n,run%o_nth))
    if(down%amy%n.gt.1) &
    allocate(cos_taper_filtery(down%amy%n,run%o_nth))
    cos_taper_filterx=1.
    cos_taper_filtery=1.
    a=4*size(cos_taper_filterx)+4*size(cos_taper_filtery)
    down%mem=down%mem+a

#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'Taper.F90 cmplxtaper: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif

#ifdef DEBUG
    if(run%depar) call out(unit)
#endif
  end subroutine cos_taper_cmplx_allocate
!----------------------------------------------------------------
  subroutine cos_taper_cmplx_deallocate()
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='cos_taper_cmplx_deallocate'
    if(run%depar) call in(unit)
#endif
    a=4*size(cos_taper_filterx)+4*size(cos_taper_filtery)
    down%mem=down%mem-a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'Taper.F90 cmplxtaper: deallocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
    deallocate(cos_taper_filterx)

    if(down%amy%n.gt.1) &
    deallocate(cos_taper_filtery)
#ifdef DEBUG
    if(run%depar) call out(unit)
#endif
  end subroutine cos_taper_cmplx_deallocate
!----------------------------------------------------------------
! cos^2 taper
                                                !! 0       Ny         -dk   
                                                !! |        |          |
  subroutine cos_taper_cmplx(wfld,iws,ith,w)    !!-----\________/------      
    complex,dimension(:,:,:), allocatable::wfld !!   1x,0x    0n,1n    
    integer :: i1_x,i0_x,i1_n,i0_n              !!1x,1n are pole location indices
    integer :: i1,i2,n1
    integer,intent(in)::ith,iws
    real,intent(in)::w
    real    :: pi2,d1,scale,max0,max1,min0,min1,firstk,lastk
#ifdef DEBUG
    character(len=128), parameter :: unit='cos_taper_cmplx'
    if(run%depar) call in(unit)
#endif

    pi2=atan(1.)*2. !pi/2.
    n1=down%amx%n

    max0=down%pmax_x*w
    min0=down%pmin_x*w

    scale=0.03*(max0-min0)

    max1=max0-scale
    min1=min0+scale

    lastk=kmap%kx(n1)
    firstk=kmap%kx(1)
    d1=(kmap%kx(2)-kmap%kx(1))

    !!roll down 
    i1_x = min(n1/2-1,max(1,nint((max1-firstk)/d1+1.)))
    i0_x = min(n1/2,max(i1_x+1,nint((max0-firstk)/d1+1.)))

    !!roll up
    i0_n = min(n1-1,max(n1/2+1,n1+1-nint( abs(min0-lastk)/d1+1.) ))
    i1_n = min(n1,max(i0_n+1,n1+1-nint( abs(min1-lastk)/d1+1.) ))

    if (i0_n-i0_x .gt. 1 )then 

       !         if (i1_x .ge. i0_x ) then 
       !            call erexit("Gonna have NANs on the cmplx_x_filter roll-down side (bitch). ")
       !         end if
       !         if (i0_n .ge. i1_n ) then 
       !            call erexit("Gonna have NANs on the cmplx_x_filter roll-up side (bitch). ")
       !         end if
       !
       !         if ( i0_n<n1/2 .or. i0_n>n1-1 ) then
       !            call erexit("invalid pmin_x to init_cos_taper_cmplx")
       !         end if
       !         if ( i1_n<n1/2+1 .or. i1_n>n1 ) then
       !            call erexit("invalid pmin_x1 in init_cos_taper_cmplx")
       !         end if
       !         if ( i0_x<2 .or. i0_x>n1/2 ) then
       !            call erexit("invalid pmax_x to init_cos_taper_cmplx")
       !         end if
       !         if ( i1_x<1 .or. i1_x>n1/2-1 ) then
       !            call erexit("invalid pmax_x1 in init_cos_taper_cmplx")
       !         end if

       cos_taper_filterx(1:i0_x,ith)=1.
       cos_taper_filterx(i0_x:i0_n,ith)=0.
       cos_taper_filterx(i0_n:,ith)=1.

       if(down%amy%n.gt.5) then
          do i1=i0_n,i1_n 
             cos_taper_filterx(i1,ith)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))
          end do
          do i1=i1_x,i0_x 
             cos_taper_filterx(i1,ith)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))
          end do
       else
          do i1=i0_n,i1_n 
             cos_taper_filterx(i1,ith)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))**2
          end do
          do i1=i1_x,i0_x 
             cos_taper_filterx(i1,ith)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))**2
          end do
       end if
       do i2=1,down%amy%n
          do i1=1,down%amx%n
             wfld(i1,i2,iws)=cos_taper_filterx(i1,ith)*wfld(i1,i2,iws)
          end do
       end do

!       write(run%io,*)' ith, w',ith,w
!       write(run%io,*)' size(wfld,1),size(filterx,1),2',size(wfld,1),size(cos_taper_filterx,1),size(cos_taper_filterx,2)
!
!       do i1=1,size(wfld,1)
!          write(run%io,*) i1,'X:',ith,':',cos_taper_filterx(i1,ith)
!       end do
!!       if (slo%source) then
!          if(iws.eq.1) then  
!             write(run%io,*)'X DIMENSION'
!             write(run%io,*) '   i1_x, i0_x', i1_x, i0_x
!             write(run%io,*) '   i0_n,i1_n', i0_n,i1_n
!             write(run%io,*) '   w',w
!             write(run%io,*) '   max0,max1',max0,max1
!             write(run%io,*) '   min0,min1',min0,min1
!             write(run%io,*) ' filt',ith,':',sum( cos_taper_filterx(:,ith))
!          end if
!!       end if
!call seperr("asdf")
    end if
    !----------------------------------------------------------------
    !!Y direction
    if (down%amy%n.gt.5) then

       n1=down%amy%n

       max0=down%pmax_y*w
       min0=down%pmin_y*w

       scale=0.075*(max0-min0)
       max1=max0-scale
       min1=min0+scale

       lastk=kmap%ky(n1)
       firstk=kmap%ky(1)
       d1=(kmap%ky(2)-kmap%ky(1))

       !!roll up
       i0_n = min(n1-1,max(n1/2+1,n1+1-nint( abs(min0-lastk)/d1+1.) ))
       i1_n = min(n1,max(i0_n+1,n1+1-nint( abs(min1-lastk)/d1+1.) ))

       !!roll down 
       i1_x = min(n1/2-1,max(1,nint((max1-firstk)/d1+1.)))
       i0_x = min(n1/2,max(i1_x+1,nint((max0-firstk)/d1+1.)))

       if (i0_n-i0_x .gt. 1 )then 
!          if (i1_x .ge. i0_x ) then 
!             call erexit("Gonna have NANs on the cmplx_y_filter roll-down side (bitch). ")
!          end if
!          if (i0_n .ge. i1_n ) then 
!             call erexit("Gonna have NANs on the cmplx_y_filter roll-up side (bitch). ")
!          end if
!          
!          if ( i0_n<n1/2 .or. i0_n>n1-1 ) then
!             call erexit("invalid pmin_y to init_cos_taper_cmplx")
!          end if
!          if ( i1_n<n1/2+1 .or. i1_n>n1 ) then
!             call erexit("invalid pmin_y1 to init_cos_taper_cmplx")
!          end if
!          if ( i0_x<2 .or. i0_x>n1/2 ) then
!             call erexit("invalid pmax_y to init_cos_taper_cmplx")
!          end if
!          if ( i1_x<1 .or. i1_x>n1-1 ) then
!             call erexit("invalid pmax_y1 to init_cos_taper_cmplx")
!          end if

          cos_taper_filtery(1:i0_x,ith)=1.
          cos_taper_filtery(i0_x:i0_n,ith)=0.
          cos_taper_filtery(i0_n:,ith)=1.

          do i1=i0_n,i1_n 
             cos_taper_filtery(i1,ith)=((1.-cos( (real(i1-i0_n)/real(i1_n-i0_n))*pi2) ))
          end do
          do i1=i1_x,i0_x 
             cos_taper_filtery(i1,ith)=((1.-cos( pi2-(real(i1-i1_x)/real(i0_x-i1_x))*pi2) ))
          end do

          do i2=1,down%amy%n
             do i1=1,down%amx%n
                wfld(i1,i2,iws)=cos_taper_filtery(i2,ith)*wfld(i1,i2,iws)
             end do
          end do

!          write(run%io,*)'Y DIMENSION'
!          write(run%io,*) '   i1_x, i0_x', i1_x, i0_x
!          write(run%io,*) '   i0_n,i1_n', i0_n,i1_n
!          write(run%io,*) '   w',w
!          write(run%io,*) '   max0,max1',max0,max1
!          write(run%io,*) '   min0,min1',min0,min1
!          write(run%io,*) ' filt',ith,':',sum( cos_taper_filtery(:,ith))
!          write(run%io,*) ' wfld',cabs(sum(wfld))
!          do i1=1,size(wfld,2)
!             write(run%io,*) i1,'Y:',ith,':',cos_taper_filtery(i1,ith)
!          end do
!!$OMP BARRIER
!          call seperr("fuckoff!")
       end if
    end if

#ifdef DEBUG
    if(run%depar) call out(unit)
#endif

  end subroutine cos_taper_cmplx
!----------------------------------------------------------------
  subroutine taper_report()
    integer:: a,b
    if(run%verb.le.2) return
    write(run%io,*)'        ************************************************'
    write(run%io,*)'        *   Taper params (index in workspace wavefield)'
    if (down%amx%n.gt.1) then
       a=down%minx1
       b=1+down%amx%n-down%maxx1
       write(run%io,'(a,i2,a,i5)')'         *          x poles at: ',a,' and',b
    end if
    if (down%amy%n.gt.1) then
       a=down%miny1
       b=1+down%amy%n-down%maxy1
       write(run%io,'(a,i2,a,i5)')'         *          y poles at: ',a,' and',b
    end if
  end subroutine taper_report
!----------------------------------------------------------------
end module taper_mod
