module RWE2D_MODEL_wemig
  use sep            !! basic SEPlib routine I/O module
  use util_          !! handle axes module
  use Kernel2D       !! GOGI's 80 degree code
  use fft2D_mod      !! Set up FFTW

  implicit none  
  integer, private :: ne1,ne3,nw,nxtap,nx1,nx3
  integer, private :: ns1,nr1,method,nh1,na1
  real, private :: de1,de3,dw,oe1,oe3, ow,dr1,or1,oh1,oa1,da1
  real, private ::  pi,r0,eps,ox1,ox3,dx1,dx3,os1,dh1
  real, private :: ds1,minang,maxang
  real, allocatable,private :: ww(:),tap(:),Timg(:,:,:,:)
  real, allocatable,private :: vmin(:)
  complex, allocatable, private :: rax(:,:)
  logical, private :: verbose,source_norm,timing
  !!
contains
  !----------------------------------------------------------------
  !! . . Initialization routine - parameter passing and array allocation
  subroutine wemig_axes_init(aE1,aE3,aW,aX1,aX3,aS1,aR1,aH1,aA1)
    type(myaxis)          :: aE1,aE3,aW,aX1,aX3,aS1,aR1,aH1,aA1

    !! . . RWE
    ne1=aE1%n; ne3=aE3%n;
    de1=aE1%d; de3=aE3%d; 
    oe1=aE1%o; oe3=aE3%o; 

    !! . . Cartesian
    nx1=aX1%n; nx3=aX3%n
    ox1=aX1%o; ox3=aX3%o
    dx1=aX1%d; dx3=aX3%d

    !! . . Input data
    nr1=aR1%n
    dr1=aR1%d
    or1=aR1%o

    !! . . Shot axis
    ns1=aS1%n
    os1=aS1%o
    ds1=aS1%d
    
    !! . . Subsurface offset
    nh1=aH1%n
    dh1=aH1%d
    oh1=aH1%o
 
    !! . . Subsurface angle
    na1=aA1%n
    da1=aA1%d
    oa1=aA1%o

    !! . . Frequency
    pi=acos(-1.)
    nw=aW%n; dw=aW%d*2.*pi; ow=aW%o*2.*pi

  end subroutine wemig_axes_init
 
 !----------------------------------------------------------------
  !! . . Initialization routine - parameter passing and array allocation
  subroutine wemig_param_init(verbose_in,r0_in,minang_in,maxang_in,eps_in,method_in)
    integer :: method_in
    real    :: r0_in,rbeg_in,rend_in,eps_in,minang_in,maxang_in
    logical :: verbose_in

    !! For OMP
    integer, external :: omp_get_num_threads,omp_get_thread_num
    integer  :: nth,ith,i1,ie3,ie1,ih1

    pi=acos(-1.)
    eps=eps_in    
    method = method_in
    verbose=verbose_in;    
    r0=r0_in
    minang = minang_in
    maxang = maxang_in
    call from_param("nxtap" ,nxtap, 40) !! Width of cosine taper
    call from_param("source_norm",source_norm,.false.)
    call from_param("timing",timing,.false.)


    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

    !! Velocity and fields
    allocate( vmin(ne3) )
    allocate( tap(ne1),rax(ne1,nth),ww(nth))
    allocate( Timg(ne1,ne3,nh1,nth) )

    !! Initalize
    !$OMP PARALLEL DO PRIVATE(ie3,ih1,ie1) 
    do ih1=1,nh1
       do ie3=1,ne3
          do ie1=1,ne1
             Timg(ie1,ie3,ih1,1)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    if (nth .gt. 1) then
       !$OMP PARALLEL DO PRIVATE(ith,ie3,ih1,ie1) 
       do ith=2,nth
          do ih1=1,nh1
             do ie3=1,ne3
                do ie1=1,ne1
                   Timg(ie1,ie3,ih1,ith)=0.
                end do
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
    call rwetaper_init() !! Initialize the taper program

  end subroutine wemig_param_init
  
  !----------------------------------------------------------------
  !! . . Main wavefield extrapolation and imaging routine
  subroutine wemig(rwf,RImg,ishot,RVel,fields)
    integer :: iw,ie1,ih1,ie3,inn,is,id,ip,ishot,iih
    real    :: vtest,rarg,sarg,ss,w
    real    :: RImg(:,:,:),fields(:,:),RVel(:,:)
    complex :: rwf(:,:)
    logical :: planewave

    !! For OMP
    integer, external :: omp_get_num_threads,omp_get_thread_num
    integer  :: nth,ith
    logical :: dofilt

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

    call from_param("planewave",planewave,.false.)

    sarg =-1.; 
    rarg = 1.

    !! Get Minimum values for the Fourier filtering
    do ie3=1,ne3
       vmin(ie3) = minval(RVel(:,ie3)/fields(:,ie3))
    end do

    if (method .eq. 1) then
       call shot_ker_init(nth, ne1,    de1,   de3, 0., 0., .11)
    else if  (method .eq. 2) then
       call shot_ker_init(nth, ne1,r0* de1,r0*de3, 0., 0., .11)
    end if
 
    !! . . Initialize FFT2D algorithm
    call fftw_init(rax,rax)

    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(iw,id,ie3,ih1,ie1) &
    SHARED(nw,ow,dw,ne3,nh1,ne1,rarg,sarg,Timg,rax,fields,tap)
    Frequency_loop: do iw=1,nw  !! FREQUENCY LOOP         
       id = omp_get_thread_num()+1

       ww(id) = dble(iw-1.)*dw+ow

       !! . . Define Receiver wavefield 
       do ie1=1,ne1
          rax(ie1,id) = rwf(ie1,iw)
       end do

       !! . . Taper wavefields
       call rwetaper(rax,id) 

       Extrap_loop: do ie3=2,ne3  

          !! . . Step Receiver wavefields
          call shot_ker_onestep(id,rax,(Rvel(:,ie3)+Rvel(:,ie3-1))/2.,&
          (fields(:,ie3)+fields(:,ie3-1))/2.,ww(id),rarg)

          !! . . Taper
          call rwetaper(rax,id) 

          !! . . High-angle filter  FFT to kx
          call RAX_Forward(rax,id)

          !! . . Phase correction routine
          call wemig_phs_correction(id,rax,ww(id),vmin(ie3),rarg,&
          (fields(:,ie3)+fields(:,ie3-1))/2.)

          !! . . FFT to x
          call RAX_Backward(rax,id)

          !! . . Taper
          call rwetaper(rax,id) 

          !! . . Imaging condition
          do ie1=1,ne1
               Timg(ie1,ie3,1,id)=Timg(ie1,ie3,1,id)+real(rax(ie1,id))
          end do

       end do Extrap_loop

       if (verbose .and. id .eq. 1 ) &
       write(0,'(a,i5,i5,i5,2f11.6,2f11.6)')  &
       '         * ST: F%,is,ns,w,im: ',100*iw/nw,ishot,ns1,& 
       ww(id)/2./pi,maxval(Timg(:,:,:,id))

    end do Frequency_loop
    !$OMP END PARALLEL DO

    !! SUMMATION OVER OUTSIDE IMAGE CONDITION
    do ith=1,nth
       !$OMP PARALLEL DO PRIVATE(ie3,ih1,ie1) 
       do ih1=1,nh1
          do ie3=1,ne3
             do ie1=1,ne1
                RImg(ie1,ie3,ih1)=RImg(ie1,ie3,ih1) + Timg(ie1,ie3,ih1,ith)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end do

    !$OMP PARALLEL DO PRIVATE(ie3,ih1,ie1,ith) 
    do ith=1,nth
       do ih1=1,nh1
          do ie3=1,ne3
             do ie1=1,ne1
                Timg(ie1,ie3,ih1,ith)=0.
             end do
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !write(0,*) 'WRITING OUT RImage.H'
    !call srite("RImage.H",RImg,4*ne3*ne1*nh1)

    call shot_ker_release()

  end subroutine wemig

  !----------------------------------------------------------------
  !! . . Apply Taper to wavefield
  subroutine rwetaper(dax,id)
    complex :: dax(:,:)
    integer :: ie1,id
    do ie1=1,ne1
       dax(ie1,id) = dax(ie1,id) * tap(ie1)
    end do
  end subroutine rwetaper

  !----------------------------------------------------------------
  !! . . Read in the parameters required for setting up geometry...etc
  subroutine Read_In_Parameters(aE3,verbose,norm,nsx,nsz,eps,minang,&
                                maxang,method,newimage,nh)
    type(myaxis) :: aE3
    integer:: norm,nsx,nsz,method,nh
    real :: eps,minang,maxang,amax,b,c
    logical :: verbose,newimage
    call from_param("norm",norm,1)  
    call from_param("nsx",nsx,2)
    call from_param("nsz",nsz,2)
    call from_param("ne3",ne3,500)
    call from_param("oe3",oe3,0.) 
    call from_param("de3",de3,0.003)
    call from_param("eps",eps,0.001)
    call from_param("minang",minang, 10.)
    call from_param("maxang",maxang,170.)
    call from_param("method",method,3)
    call from_param("verbose",verbose,.false.) 
    call from_param("newimage",newimage,.false.)
    call from_Param("nh",nh,1)
    aE3%n=ne3
    aE3%o=oe3
    aE3%d=de3
  end subroutine Read_In_Parameters

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine RWF_Info(aE1,aW,aS1,aR1)
    type(myaxis) ::   aE1,aW,aS1,aR1
    character(len=128) :: name
    name = 'rwf'
    call ainaxis(aE1,1,name)
    call ainaxis(aW ,2,name)
    call ainaxis(aS1,3,name)
    call ainaxis(aR1,1,name)

  end subroutine RWF_Info

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine IMAGE_Out_Info(aX1,aX3,aA1)
    type(myaxis) ::         aX1,aX3,aA1
    character(len=128) :: name
    name = 'image'
    call aouaxis(aX3,1,name)
    call aouaxis(aX1,2,name)
    call aouaxis(aA1,3,name)
    call auxputch("esize","i",4,name)
  end subroutine IMAGE_Out_Info

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine IMAGE_Info(aX1,aX3,aA1)
    type(myaxis) ::     aX1,aX3,aA1
    character(len=128) :: name
    name = 'image'
    call ainaxis(aX3,1,name)
    call ainaxis(aX1,2,name)
    call ainaxis(aA1,3,name)
  end subroutine IMAGE_Info

 !----------------------------------------------------------------
 !! . . Read in the header information describing Cartesian velocity/image space
  subroutine VEL_Info(aX1,aX3)
    type(myaxis) ::   aX1,aX3
    character(len=128) :: name
    name = 'vel'
    call ainaxis(aX3,1,name)
    call ainaxis(aX1,2,name)

  end subroutine VEL_Info

 !----------------------------------------------------------------
 !! . . Write out header information on the image
  subroutine write_out(aX1,aX3,aE1,aE3,aS1,aR1,aH1,aA1)
    type(myaxis) ::    aX1,aX3,aE1,aE3,aS1,aR1,aH1,aA1
    character (len=128 ) :: name

    !! . . Create output header file associated with RImage
    name= "RImage.H"
    call aouaxis(aE1,1,name)
    call aouaxis(aE3,2,name)
    call aouaxis(aA1,3,name)
    call aouaxis(aS1,4,name)
    call auxputch("esize","i",4,name)

    !! . . Create output header file associated with RVel
    name= "RVel.H"
    call aouaxis(aE1,1,name)
    call aouaxis(aE3,2,name)
    call aouaxis(aS1,3,name)
    call auxputch("esize","i",4,name)

    !! . . Create output header file associated with the rays (QC)
    name= "rays"
    call aouaxis(aE1,1,name)
    call aouaxis(aE3,2,name)
    call auxputch("n3","i",2,name)
    call auxputch("o3","r",1.,name)
    call auxputch("d3","r",1.,name)
    call auxputch("esize","i",4,name)

    !! . . Write(0,*) out taper
    name="taper.H"
    call aouaxis(aE1,1,name)
    call auxputch("esize","i",4,name)

  end subroutine write_out

 !----------------------------------------------------------------
 !! . . Zero a bunch of important arrays
  subroutine Zero_arrays(RVel,RImg,ncz,ncx,nch,aE1,aE3,aH1)
    type(myaxis) :: aE1,aE3,aH1
    integer :: ix,iz,ih,ie1,ie3,ncz,ncx,nch
    real    :: RVel(:,:),RImg(:,:,:)

    !$OMP PARALLEL DO PRIVATE(ie3,ie1)
    do ie3=1,aE3%n
       do ie1=1,aE1%n
          RVel(ie1,ie3)=0.
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(ie3,ih,ie1)
    do ih=1,aH1%n
       do ie3=1,aE3%n
          do ie1=1,aE1%n
             RImg(ie1,ie3,ih)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO


  end subroutine Zero_arrays
  
 !----------------------------------------------------------------
 !! . . Set up the Cartesian geometry for migration
  subroutine CART_Setup(raysin,fields,gnorm,n3,o3,d3)
    integer :: iE1,iE2,iE3,n3
    real:: o3,d3
    real :: raysin(:,:,:),fields(:,:),gnorm(:,:)


    !! Do Cartesian coordinate setup
    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          raysin(iE1,iE3,1) =  dble(iE1-1.)*de1 + oe1 + os1 
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          raysin(iE1,iE3,2) =  dble(iE3-1.)*de3 + oe3 
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          fields(iE1,iE3)= 1.
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          gnorm(iE1,iE3)= 1.
       end do
    end do
    !$OMP END PARALLEL DO

    write(0,*) '        ************************************************'
    write(0,*) '        * '
    write(0,*) '        *      2D  Cartesian Migration'
    write(0,*) '        * '
    write(0,*) '        *            written by Jeff Shragge ' 
    write(0,*) '        *  '
    write(0,*) '        ************************************************'
    write(0,*) ' ' 
    write(0,*) '        ************************************************'
    write(0,*) '        * CARTESIAN GRID CHECK '
    write(0,'(a,2f11.4,2f11.4)') '         * X-dir',minval(raysin(:,:,1)),maxval(raysin(:,:,1))
    write(0,'(a,2f11.4,2f11.4)') '         * Z-dir',minval(raysin(:,:,2)),maxval(raysin(:,:,2))
    write(0,*) '        ************************************************'
    write(0,*) ' '

  end subroutine CART_Setup
  
 !----------------------------------------------------------------
 !! . .Set up the cylindrical elliptic  coordinate system for migration
  subroutine ELLIP_Setup(raysin,fields,gnorm,n1,o1,d1)
    integer              :: iE1,iE3,n1
    real                 :: o1,d1
    real                 :: fields(:,:),raysin(:,:,:),gnorm(:,:)
    real,    allocatable :: rline(:),tline(:)
    allocate( tline(ne3),rline(ne1) )

    !! . . Define the extrapolation axis
    d1 = (maxang-minang)/180.*pi/dble(ne1-1.)
    o1 = pi*minang/180.
    r0 = dble(ne1-1.)*de1/2.

    rline=(/ (iE1,iE1=0,ne1-1) /)*d1+o1
    tline=(/ (iE3,iE3=0,ne3-1) /)*de3+oe3

    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          raysin(iE1,iE3,1)=-abs(r0)*cosh(tline(iE3))*cos(rline(iE1))+oe1+os1+r0
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          raysin(iE1,iE3,2)=abs(r0)*sinh(tline(iE3))*sin(rline(iE1))
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . Slowness stretch factor
    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          fields(iE1,iE3)= sqrt(sinh(tline(iE3))**2+sin(rline(iE1))**2)
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . GNORM
    !$OMP PARALLEL DO PRIVATE(iE3,iE1)
    do iE3=1,ne3
       do iE1=1,ne1
          gnorm(iE1,iE3)= sqrt(sinh(tline(iE3))**2+sin(rline(iE1))**2)
       end do
    end do
    !$OMP END PARALLEL DO


    write(0,*) '        ************************************************'
    write(0,*) '        * 2D RWE migration in ELLIPIC coordinates '
    write(0,*) '        * '
    write(0,*) '        *       written by Jeff Shragge ' 
    write(0,*) '        *  '
    write(0,*) '        ************************************************'
    write(0,*) ' ' 
    write(0,*) '        ************************************************'
    write(0,*) '        * ELLIPTIC GRID CHECK:'
    write(0,'(a,2f11.4,2f11.4)') '         * MIN/MAX Angles  :',minval(rline)*180./pi,maxval(rline)*180./pi
    write(0,'(a,2f11.4,2f11.4)') '         * MIN/MAX Ex step :',minval(tline),maxval(tline)
    write(0,'(a,2f11.4,2f11.4)') '         * Rline min/max   :', minval(rline)*180./pi,maxval(rline)*180./pi
    write(0,*) '        ************************************************'
    write(0,*) ' '
    deallocate( tline, rline )

  end subroutine ELLIP_Setup

 !----------------------------------------------------------------
 !! . . Write out information for visual verification
  subroutine Geometry_QC(raysin,aE1,aE3,aS1,aW,aR1,aH1,aA1)
    type(myaxis)                           :: aE1,aE3,aS1,aW,aR1,aH1,aA1
    real :: raysin(:,:,:)
    !! . . Screen output for QC
    write(0,*) '        ************************************************'
    write(0,*) '        * GEOMETRY CHECK'
    write(0,'(a,2f11.4,2f11.4)')    "         * E1-coord min/max         :", minval(raysin(:,1,1)),maxval(raysin(:,1,1))
    write(0,'(a,2f11.4,2f11.4)')    "         * E3-coord min/max         :", minval(raysin(:,:,2)),maxval(raysin(:,:,2))
    write(0,'(a,i4,2f11.4,2f11.4)') '         * RECEIVERS inline    n,o,d:',ne1,oe1,de1
    write(0,'(a,i4,2f11.4,2f11.4)') '         * EXTRAPOLATION STEPS n,o,d:',ne3,oe3,de3
    write(0,'(a,i4,2f11.4,2f11.4)') '         * SHOTS inline        n,o,d:',aS1%n,aS1%o,aS1%d
    write(0,'(a,i4,2f11.4,2f11.4)') '         * FREQUENCIES         n,o,d:',aW%n,aW%o,aW%d
    write(0,'(a,i4,2f11.4,2f11.4)') '         * Subsurface H        n,o,d:',aH1%n,aH1%o,aH1%d
    write(0,'(a,i4,2f11.4,2f11.4)') '         * Subsurface A        n,o,d:',aA1%n,aA1%o,aA1%d
    write(0,*) '        ************************************************'
    write(0,*) ' '
  end subroutine Geometry_QC

 !----------------------------------------------------------------
 !! . . Write out velocity information for visual verification
  subroutine Velocity_Check(CVel,RVel)
    real :: RVel(:,:),CVel(:,:)
    !! . . Screen output for QC
    write(0,*) '        *************************************************'
    write(0,*) '        * VELOCITY CHECK'
    write(0,'(a,i6,i6,2f11.4,2f11.4)')      '         * Cart VEL n1,n3,min,max:',size(CVel,2),size(CVel,1),minval(CVel),maxval(CVel)
    write(0,'(a,i6,i6,2f11.4,2f11.4)')      '         * RWE  VEL n1,n3,min,max:',size(RVel,1),size(RVel,2),minval(RVel),maxval(RVel)
    write(0,*) '        ************************************************'
    write(0,*) ' '
  end subroutine Velocity_Check

 !----------------------------------------------------------------
 !! . . Reposition the wavefield on elliptic coordinates
  subroutine ELLIP_Position_XWF(xwfin,xwfout,raysin,shotloc)
    integer :: iE1,iX1,jx
    real :: raysin(:,:,:),x,dx1,dx2,dd1,dd2,shotloc,r,s,xx,ddx,dr
    complex  :: xwfin(:,:),xwfout(:,:)

!    do iE1 = 1,ne1
       !! . . Generate interpolation coefficients
       !! . . independent of y-direction
!       x = raysin(iE1,1,1) - shotloc
!       iX1 = floor( ( x - or1 ) / dr1) + 1.5
!       dx1 = x - (or1 + (iX1-1.)*dr1 )
!       dx2 = (or1 + dble(iX1)*dr1 ) - x
!       dd1 = dx2/(dx1+dx2)
!       dd2 = dx1/(dx1+dx2)
!       !! . . Boundary Check
!       if (iX1 .ge. ne1) iX1 = ne1-1
!       if (iX1 .le. 1 )  iX1 = 1
!       xwfout(iE1,:) = dd1*xwfin(iX1,:)+dd2*xwfin(iX1+1,:)
!     end do
    dr = sqrt(dr1**2)
    do iE1=1,ne1
       x = raysin(iE1,1,1) - shotloc
       iX1 = floor( ( x - or1 ) / dr1) + 1.5
	   if (iX1 .lt. 1 .or. iX1 .gt. nr1 ) cycle
       do jx=iX1-2, iX1+2
          if ( jx .lt. 1 .or. jx .gt. nr1) cycle
          xx = or1 + (jx-1.)*dr1
          ddx = xx-x
          r = sqrt( ddx**2) / dr
          s = min( sin(r) / r , 1.0 )
          xwfout(iE1,:) = xwfout(iE1,:) + s*xwfin(jx,:)
       end do
    end do
    
    write(0,*) '        ************************************************'
    write(0,*) '        * WAVEFIELD CHECK'
    write(0,*) '        * MIN/MAX Wavefield: ',minval(real(xwfout)),maxval(real(xwfout))
    write(0,*) '        ************************************************'
    write(0,*) ' '
  end subroutine ELLIP_Position_XWF

  !----------------------------------------------------------------
 !! . . Find the corresponding Cartesian velocity indicies
 subroutine Get_Vel_Indicies(raysin,iix1,oo1)
   integer :: iix1(:)
   real :: oo1(:)
   real :: ee,raysin(:,:,:)

   !! X1 index
   ee = minval(raysin(:,:,1))
   iix1(1) = floor ((ee - ox1)/dx1) + 1.5
   oo1(1) =ee
   ee = maxval(raysin(:,:,1))
   iix1(2) = ceiling ((ee - ox1)/dx1) + 1.5
   oo1(2) = ee
   if (iix1(1) .le. 0) iix1(1)=1
   if (iix1(1) .gt. nx1)  iix1(1)=nx1
   if (iix1(2) .le. 0) iix1(2)=1
   if (iix1(2) .gt. nx1	) iix1(2)=nx1

 end subroutine Get_Vel_Indicies

  !----------------------------------------------------------------
  !! . . Taper Initialization Program
 subroutine rwetaper_init()
   integer :: i1,j1,ie1,i,jtap

   do ie1=1,ne1
      tap(ie1)=1.
   end do

   do i1=1,nxtap
      j1 = abs(nxtap-i1)
      tap(i1) = cos(pi/2.* j1/(nxtap-1.))
   end do

   do i1=1,nxtap
      j1 = abs(nxtap-i1)
      tap(ne1-i1+1) = cos(pi/2.* j1/(nxtap-1.))
   end do

   call srite("taper.H",tap,4*size(tap) )
 end subroutine rwetaper_init
 !----------------------------------------------------------------

end module  RWE2D_MODEL_wemig
