module OW_ctrillwei

use sep
use OWvelo_areal
use OW_type
use OW_fftw_areal
use OW_parms_areal
use OW_exec_time
use OW_rndseed

implicit none

integer,private					:: iwa, ispa
contains

!---------------------------------------------------------
subroutine owctrlwei_init(iwa_in,ispa_in)

integer						:: iwa_in,ispa_in

iwa=iwa_in
ispa=ispa_in

end subroutine
!---------------------------------------------------------
subroutine owctrlwei(taper,vsal,vel,recset,souset,GG1,GG2,rand1,rand2,wwr,wws)

complex,optional					:: wwr(:,:),wws(:,:)
real,optional					:: rand1,rand2(:)
integer						:: iz,ikx,iref,ixw,nout,seed_1,seed_2
real						:: minv,maxv,vsal,sqa,f,g,rnd_phase
real						:: vel(:,:),recset(:,:),souset(:,:),taper(:)
real, allocatable					:: vref(:),sqa0(:)!,rnd_phase(:)
complex, allocatable				:: ws(:),gs0(:),gs(:,:)
complex, allocatable				:: wr(:),gr0(:),gr(:,:)
complex, allocatable				:: tfld(:),refl(:)
complex						:: GG1(:),GG2(:)

if (allocated(ws)) deallocate(ws);allocate(ws(kxm%n));	ws=0
if (allocated(gs)) deallocate(gs);allocate(gs(kxm%n,nref));	gs=0
if (allocated(gs0)) deallocate(gs0);allocate(gs0(kxm%n));	gs0=0
if (allocated(wr)) deallocate(wr);allocate(wr(kxm%n));	wr=0
if (allocated(gr)) deallocate(gr);allocate(gr(kxm%n,nref));	gr=0
if (allocated(gr0)) deallocate(gr0);allocate(gr0(kxm%n));	gr0=0
if (allocated(sqa0)) deallocate(sqa0);allocate(sqa0(kxm%n));sqa0=0
if (allocated(vref)) deallocate(vref);allocate(vref(nref));	vref=0
if (allocated(tfld)) deallocate(tfld);allocate(tfld(kxm%n));tfld=0
if (allocated(refl)) deallocate(refl);allocate(refl(kxm%n));refl=0

!--------------------------------------------------
! Setting FFTW plans
!
call owfftw(tfld)
do iz=sog%z%n,nwbott,-1
!--------------------------------------------------
! Compute vref
!
   call owvref(nref,vsal,vel(iz,:),vref,minv,maxv)

   tfld=0.
   tfld=cmplx(recset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=0.
   tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
   wr=wr+refl

   tfld=0.
   tfld=cmplx(souset(iz,:),0); call sfftw_execute(planF);
   refl=tfld/sqrt(1.*kxm%n)
   tfld=0.
   tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
   ws=ws+refl

   if (minv .eq. maxv) then
      gr0=wr;
      gs0=ws;

!      Operations in the wavenumber domain
!      phase shift

      do ikx = 1,kxm%n

         sqa = sou%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2

         if (sqa > 0) then
            gr0(ikx)=gr0(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
            gs0(ikx)=gs0(ikx)*exp(cmplx(0,-v%z%d*sqrt(sqa)))
         else
            gr0(ikx)=gr0(ikx)*exp(-sqrt(-sqa)*v%z%d)
            gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
         end if
      end do

      tfld=gr0; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
      tfld=gs0; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)

   else
      do iref=1,nref                 
         gr(:,iref)=wr
         gs(:,iref)=ws
         do ikx = 1,kxm%n
!      phase shift
            sqa = sou%w%all(iwa)**2 / (vref(iref))**2 - kxm%all(ikx)**2;
            if (sqa > 0) then
               gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
               gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa)));
            else
               gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
               gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
            end if   
         end do

! Operations in x            

         tfld=gr(:,iref); call sfftw_execute(planI); gr(:,iref)=tfld/sqrt(1.*kxm%n)
         tfld=gs(:,iref); call sfftw_execute(planI); gs(:,iref)=tfld/sqrt(1.*kxm%n)

 !      split step
         gr(:,iref) = gr(:,iref) * exp(cmplx(0,v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
         gs(:,iref) = gs(:,iref) * exp(cmplx(0,-v%z%d*sou%w%all(iwa))*(1/vel(iz,:)-1/vref(iref)));
      end do

      do iref=1,nref-1
         where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
            wr=gr(:,iref)
            ws=gs(:,iref)
         end where
      end do
      where (vel(iz,:) >= vref(nref))
         wr=gr(:,nref)
         ws=gs(:,nref)
      end where 
   end if
   wr=wr*taper
   ws=ws*taper
   if (encode) then
      if(encz) then
          if (funct.eq.1) then
             rnd_phase=sin(2*pi*rand2(iz))
          else if (funct.eq.2) then
             rnd_phase=cos(2*pi*rand2(iz))
          else if (funct.eq.3) then
             rnd_phase=tan(2*pi*rand2(iz))
          else
             rnd_phase=rand2(iz)    
          end if
          wr=wr*exp(cmplx(0,eps*2*pi*rnd_phase))
          ws=ws*exp(cmplx(0,eps*2*pi*rnd_phase))
      end if
   end if
   if (exist_file("wwr")) then
      wwr(iz,:)=wr
      wws(iz,:)=ws
   end if
end do


if (encode) then
   if (encw) then
      wr=wr*exp(cmplx(0,(2*pi*rand1)))
      ws=ws*exp(cmplx(0,(2*pi*rand1)))
   else if(enclin) then
      wr=wr*exp(cmplx(0,(t0*sou%w%all(iwa))))
      ws=ws*exp(cmplx(0,(t0*sou%w%all(iwa))))
   end if
end if

if (nwbott.ne.1) then
call owvref(1,minval(vel(nwbott-1,:)),vel(iz,:),vref,minv,maxv)
tfld=0.
tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n)
tfld=0.
tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n)
!      Operations in the wavenumber domain
!      phase shift
   do ikx = 1,kxm%n
      sqa = sou%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2
      if (sqa > 0) then
         wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa)))
         ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa)))
      else
         wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott)
         ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott)
      end if
   end do
   tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n)
   tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
   wr=wr*taper
   ws=ws*taper
end if
GG1=wr
GG2=ws

!
deallocate(ws,gs0,gs,wr,gr0,gr,sqa0,vref,tfld)
end subroutine

!-------------------------------------------
subroutine legendre(r,g,s)

integer				:: g
real				:: r,s

if (g==0) then
   s=1
else if (g==1) then
   s=r
else if (g==2) then
   s=0.5*(3*r**2-1)
else if (g==3) then
   s=0.5*(5*r**3-3*r)
else if (g==4) then
   s=1/8*(35*r**4-3*r**2+3)
else if (g==5) then
   s=1/8*(63*r**5-70*r**3+15*r)
else if (g==6) then
   s=1/16*(231*r**6-315*r**4+105*r**2-5)
else if (g==7) then
   s=1/16*(429*r**7-693*r**5+315*r**3-35*r)
else if (g==8) then
   s=1/128*(6435*r**8-12012*r**6+6930*r**4-1260*r**2+35)
else if (g==9) then
   s=1/128*(12155*r**9-25740*r**7+18018*r**5-4620*r**3+315*r)
else if (g==10) then
   s=1/256*(46189*r**10-109395*r**8+90090*r**6-30030*r**4+3465*r**2-63)
end if

end subroutine

end module
