program OW_windowing

use sep
use OW_type
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
implicit none

logical					:: taperon,dip
character*80				:: hor
real, allocatable				:: intmp(:,:,:),temp(:,:,:)
real, allocatable				:: pick(:,:,:),taper(:)
real					:: finish,start,a
integer					:: wavz,n,i,ihor,ikx,iz,ih,it,par1(3)
integer, allocatable                              :: pick_wavzb(:,:),pick_wavzt(:,:)
real, allocatable                                 :: pick_smoo(:,:),pick_dip(:,:),pick_cos(:,:),pick_sin(:,:)

call sep_init()
call sep_begin_prog()

 call from_param("nhorz",nhorz)
 call from_history("n1",sog%xm%n);call from_history("d1",sog%xm%d);call from_history("o1",sog%xm%o)  
 call from_history("n2",sog%sh%n);call from_history("d2",sog%sh%d);call from_history("o2",sog%sh%o)
 call from_history("n3",sog%z%n);call from_history("d3",sog%z%d);call from_history("o3",sog%z%o)

 write(0,*) sog%z%n,sog%z%d,sog%z%o
 write(0,*) sog%sh%n,sog%sh%d,sog%sh%o
 write(0,*) sog%xm%n,sog%xm%d,sog%xm%o
 write(0,*) nhorz
 call to_history("n1",sog%xm%n);call to_history("d1",sog%xm%d);call to_history("o1",sog%xm%o)  
 call to_history("n2",sog%sh%n);call to_history("d2",sog%sh%d);call to_history("o2",sog%sh%o)
 call to_history("n3",sog%z%n);call to_history("d3",sog%z%d);call to_history("o3",sog%z%o)
 call to_history("n4",nhorz);call to_history("d4",1.);call to_history("o4",1.)  
 call to_history("eout",4);call to_history("label3","")

!--------------------------------------------------
!
! Read rotated CIGs
!
allocate(intmp(sog%xm%n,sog%sh%n,sog%z%n));intmp=0.
call from_param("wavz",wavz,13)
call from_param("taperon",taperon,.true.)
call from_param("dip",dip,.true.)
call sseek("in",0,0)
call sreed("in",intmp,4*sog%z%n*sog%sh%n*sog%xm%n)
call owexec_time(start)
allocate(pick(sog%xm%n,sog%sh%n,nhorz));pick=9999.
do i=1,nhorz
   if (i.eq.1) then
      call from_param('hor1',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
write(0,*) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.2) then
      call from_param('hor2',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.3) then
      call from_param('hor3',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.4) then
      call from_param('hor4',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.5) then
      call from_param('hor5',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.6) then
      call from_param('hor6',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.7) then
      call from_param('hor7',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.8) then
      call from_param('hor8',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.9) then
      call from_param('hor9',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.10) then
      call from_param('hor10',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.11) then
      call from_param('hor11',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   else if (i.eq.12) then
      call from_param('hor12',hor)
      open(51,file=hor,status='old')
      read(51,*,end=1) par1
      do n=1,par1(2)
         read(51,*,end=1) pick(par1(1)+n-1,:,i)
         pick(par1(1)+n-1,:,i)=pick(par1(1)+n-1,:,i)+par1(3)
      end do
   end if
1 close(51)
end do

if (dip) then
allocate(pick_smoo(sog%xm%n,nhorz),pick_dip(sog%xm%n,nhorz))
allocate(pick_cos(sog%xm%n,nhorz),pick_sin(sog%xm%n,nhorz))
allocate(pick_wavzb(sog%xm%n,nhorz),pick_wavzt(sog%xm%n,nhorz))
pick_wavzt=wavz
pick_wavzb=wavz
do i=1,nhorz
   a = maxval(pick(:,int(sog%sh%n/2.),i), mask = pick(:,int(sog%sh%n/2.),i).ne.9999.)
   do ikx=1,sog%xm%n
      pick_smoo(ikx,i)=sum(pick(max(1,ikx-10):min(ikx+10,sog%xm%n),int(sog%sh%n/2.),i))/21
   end do
   where(pick_smoo(:,i) > a) pick_smoo(:,i)=9999.
   do ikx=1,sog%xm%n
      if (ikx==1) then
         pick_dip(ikx,i)=(pick_smoo(ikx+1,i)-pick_smoo(ikx,i))/sog%xm%d*sog%z%d
      else if (ikx==sog%xm%n) then
         pick_dip(ikx,i)=(pick_smoo(ikx,i)-pick_smoo(ikx-1,i))/sog%xm%d*sog%z%d
      else
         pick_dip(ikx,i)=(pick_smoo(ikx+1,i)-pick_smoo(ikx-1,i))/(2*sog%xm%d)*sog%z%d
!         pick_dip(ikx,i)=(pick_smoo(ikx-2,i)-8*pick_smoo(ikx-1,i)+8*pick_smoo(ikx+1,i)-pick_smoo(ikx+2,i))*sog%z%d/(12*sog%xm%d)
      end if
   end do
   where(pick_smoo(:,i) == 9999.) pick_dip(:,i)=0.
   where(pick_dip(:,i) > 50.) pick_dip(:,i)=0.
   pick_cos(:,i)=sqrt(1/(1+pick_dip(:,i)**2))
   pick_sin(:,i)=sqrt(1-pick_cos(:,i)**2)
   where(pick_dip(:,i) < 0) pick_sin(:,i)=-pick_sin(:,i)
   pick_sin(:,i)=pick_sin(:,i)*wavz
   pick_cos(:,i)=pick_cos(:,i)*wavz

   do ikx=1,sog%xm%n
      if (pick_sin(ikx,i)<0) then
         if (ikx+int(pick_sin(ikx,i))>=1) then
            if (pick_wavzb(ikx+int(pick_sin(ikx,i)),i).lt.&
                pick(ikx+int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)-&
                pick(ikx,int(sog%sh%n/2.),i)-int(pick_cos(ikx,i))) then
                pick_wavzb(ikx+int(pick_sin(ikx,i)),i)=pick(ikx+int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)-&
                pick(ikx,int(sog%sh%n/2.),i)-int(pick_cos(ikx,i))
            end if
            if (pick_wavzt(ikx-int(pick_sin(ikx,i)),i).lt.&
                -pick(ikx-int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)+&
                pick(ikx,int(sog%sh%n/2.),i)+int(pick_cos(ikx,i))) then
                pick_wavzt(ikx-int(pick_sin(ikx,i)),i)=-pick(ikx-int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)+&
                pick(ikx,int(sog%sh%n/2.),i)+int(pick_cos(ikx,i))
            end if
         end if
      else
         if (ikx+int(pick_sin(ikx,i))<=sog%xm%n) then
            if (pick_wavzb(ikx+int(pick_sin(ikx,i)),i).lt.&
                pick(ikx+int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)-&
                pick(ikx,int(sog%sh%n/2.),i)-int(pick_cos(ikx,i))) then
                pick_wavzb(ikx+int(pick_sin(ikx,i)),i)=pick(ikx+int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)-&
                pick(ikx,int(sog%sh%n/2.),i)-int(pick_cos(ikx,i))
            end if
            if (pick_wavzt(ikx-int(pick_sin(ikx,i)),i).lt.&
                -pick(ikx-int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)+&
                pick(ikx,int(sog%sh%n/2.),i)+int(pick_cos(ikx,i))) then
                pick_wavzt(ikx-int(pick_sin(ikx,i)),i)=-pick(ikx-int(pick_sin(ikx,i)),int(sog%sh%n/2.),i)+&
                pick(ikx,int(sog%sh%n/2.),i)+int(pick_cos(ikx,i))
            end if
         end if
      end if
   end do
   where(pick_smoo(:,i).eq.9999.) 
       pick_wavzt(:,i)=wavz
       pick_wavzb(:,i)=wavz
   end where

end do
else
allocate(taper(2*wavz+1))
if (taperon) then
   call owtaper_init(wavz)
   call owtaper(taper)
else
   taper=1
end if
end if
allocate(temp(sog%xm%n,sog%sh%n,sog%z%n));
do i=1,nhorz
   temp=0.
   do iz=1,sog%z%n
      do ih=1,sog%sh%n
         do ikx=1,sog%xm%n
            if (dip) then
               if (iz .ge. int(pick(ikx,ih,i))-pick_wavzb(ikx,i) .and. iz .le. int(pick(ikx,ih,i))+pick_wavzt(ikx,i)) &
               temp(ikx,ih,iz)=intmp(ikx,ih,iz)
            else
               if (iz .ge. int(pick(ikx,ih,i))-wavz .and. iz .le. int(pick(ikx,ih,i))+wavz) &
               temp(ikx,ih,iz)=intmp(ikx,ih,iz)*taper(iz-int(pick(ikx,ih,i))+wavz+1)**2
            end if
         end do
      end do
   end do
   call srite("out",temp,4*sog%z%n*sog%xm%n*sog%sh%n)
end do

call owexec_time(finish)
write(0,*)  "FINISHED WINDOWING ",int(finish-start),'(s)'
deallocate(temp,pick,intmp)
if (dip) deallocate(pick_smoo,pick_dip,pick_wavzt)

call sep_end_prog()
end program
