program OW3D_windowing

use sep
use OW_type
use OW_exec_time
implicit none

logical					:: taperon,flt,dip
character*80				:: hor
real, allocatable				:: intmp(:,:,:),temp(:,:,:)
real, allocatable				:: pick(:,:,:,:),taper(:),a(:)
real					:: finish,start,pi,b
integer					:: wavz,n,i,ihor,ikx,iz,ih,it,iky,m,nhorz,ikh
integer, allocatable                              :: pick_wavzb(:,:,:),pick_wavzt(:,:,:)
real, allocatable                                 :: pick_smoo(:,:,:,:),pick_dip(:,:,:),pick_cos(:,:,:),pick_sin(:,:,:)
type(mod_grid)                                    :: sog
type(axis)                                        :: shx,xm,ym,z

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%shx%n);call from_history("d2",sog%shx%d);call from_history("o2",sog%shx%o)
 call from_history("n3",sog%ym%n);call from_history("d3",sog%ym%d);call from_history("o3",sog%ym%o)  
 call from_history("n4",sog%z%n);call from_history("d4",sog%z%d);call from_history("o4",sog%z%o)

 write(0,*) sog%z%n,sog%z%d,sog%z%o
 write(0,*) sog%shx%n,sog%shx%d,sog%shx%o
 write(0,*) sog%xm%n,sog%xm%d,sog%xm%o
 write(0,*) sog%ym%n,sog%ym%d,sog%ym%o
 call from_param("dip",dip,.true.)
 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%shx%n);call to_history("d2",sog%shx%d);call to_history("o2",sog%shx%o)
 call to_history("n3",sog%ym%n);call to_history("d3",sog%ym%d);call to_history("o3",sog%ym%o)  
 call to_history("n4",sog%z%n);call to_history("d4",sog%z%d);call to_history("o4",sog%z%o)
 call to_history("n5",nhorz);call to_history("d5",1.);call to_history("o5",1.)  
 call to_history("eout",4);call to_history("label3","")

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

allocate(taper(2*wavz+1))
pi=4*atan(1.)
taper=1.
if (dip) then
allocate(pick_smoo(sog%xm%n,sog%ym%n,sog%shx%n,nhorz),pick_dip(sog%xm%n,sog%ym%n,nhorz))
allocate(pick_cos(sog%xm%n,sog%ym%n,nhorz),pick_sin(sog%xm%n,sog%ym%n,nhorz))
allocate(pick_wavzb(sog%xm%n,sog%ym%n,nhorz),pick_wavzt(sog%xm%n,sog%ym%n,nhorz))
pick_wavzt=wavz
pick_wavzb=wavz
allocate(a(sog%xm%n)); a=0
b=0
do i=1,nhorz
   do ikx=1,sog%xm%n
      a(ikx) = maxval(pick(:,ikx,int(sog%shx%n/2.),i), mask = pick(:,ikx,int(sog%shx%n/2.),i).lt.9999.)      
   end do
   b=maxval(a)
   do ikx=1,sog%xm%n
      do iky=1,sog%ym%n
         do ikh=1,sog%shx%n
            pick_smoo(ikx,iky,ikh,i)=sum(pick(max(1,iky-10):min(iky+10,sog%ym%n),max(1,ikx-10):min(ikx+10,sog%xm%n),ikh,i))/441
         end do
      end do
   end do
write(0,*) "passei 00"
   do ikx=1,sog%xm%n
      do iky=1,sog%ym%n
         do ikh=1,sog%shx%n
            if (pick_smoo(ikx,iky,ikh,i) > b) then
               pick_smoo(ikx,iky,ikh,i)=9999.
            end if
         end do
      end do
   end do
write(0,*) "passei 01"
   do ikx=1,sog%xm%n
      do iky=1,sog%ym%n
         if (ikx==1) then
            pick_dip(ikx,iky,i)=(pick_smoo(ikx+1,iky,int(sog%shx%n/2.),i)-pick_smoo(ikx,iky,int(sog%shx%n/2.),i))/sog%xm%d*sog%z%d
         else if (ikx==sog%xm%n) then
            pick_dip(ikx,iky,i)=(pick_smoo(ikx,iky,int(sog%shx%n/2.),i)-pick_smoo(ikx-1,iky,int(sog%shx%n/2.),i))/sog%xm%d*sog%z%d
         else
            pick_dip(ikx,iky,i)=(pick_smoo(ikx+1,iky,int(sog%shx%n/2.),i)-pick_smoo(ikx-1,iky,int(sog%shx%n/2.),i))/(2*sog%xm%d)*sog%z%d
         end if
         if (mod(ikx,40).eq.1.and.mod(iky,40).eq.1) write(0,*) ikx,iky,pick_dip(ikx,iky,i)
      end do
   end do
   do iky=1,sog%ym%n
      where(pick_smoo(:,iky,int(sog%shx%n/2.),i) == 9999.) pick_dip(:,iky,i)=0.
!      where(pick_dip(:,iky,i) > 50.) pick_dip(:,iky,i)=0.
      pick_cos(:,iky,i)=sqrt(1/(1+pick_dip(:,iky,i)**2))
      pick_sin(:,iky,i)=sqrt(1-pick_cos(:,iky,i)**2)
      where(pick_dip(:,iky,i) < 0) pick_sin(:,iky,i)=-pick_sin(:,iky,i)
      pick_sin(:,iky,i)=pick_sin(:,iky,i)*wavz
      pick_cos(:,iky,i)=pick_cos(:,iky,i)*wavz
   end do
write(0,*) "passei 02"
   do ikx=1,sog%xm%n
      do iky=1,sog%ym%n
         if (pick_sin(ikx,iky,i)<0) then
            if (ikx+int(pick_sin(ikx,iky,i))>=1.and.ikx+int(pick_sin(ikx,iky,i))<=sog%xm%n) then
               if (pick_wavzb(ikx+int(pick_sin(ikx,iky,i)),iky,i).lt.&
                   pick_smoo(ikx+int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)-&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)-int(pick_cos(ikx,iky,i))) then
                   pick_wavzb(ikx+int(pick_sin(ikx,iky,i)),iky,i)=pick_smoo(ikx+int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)-&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)-int(pick_cos(ikx,iky,i))
               end if
            end if
            if (ikx-int(pick_sin(ikx,iky,i))>=1.and.ikx-int(pick_sin(ikx,iky,i))<=sog%xm%n) then
               if (pick_wavzt(ikx-int(pick_sin(ikx,iky,i)),iky,i).lt.&
                   -pick_smoo(ikx-int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)+&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)+int(pick_cos(ikx,iky,i))) then
                   pick_wavzt(ikx-int(pick_sin(ikx,iky,i)),iky,i)=-pick_smoo(ikx-int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)+&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)+int(pick_cos(ikx,iky,i))
               end if
            end if
write(0,*) "passei 03",ikx,iky
         else
            if (ikx+int(pick_sin(ikx,iky,i))>=1.and.ikx+int(pick_sin(ikx,iky,i))<=sog%xm%n) then
               if (pick_wavzb(ikx+int(pick_sin(ikx,iky,i)),iky,i).lt.&
                   pick_smoo(ikx+int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)-&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)-int(pick_cos(ikx,iky,i))) then
                   pick_wavzb(ikx+int(pick_sin(ikx,iky,i)),iky,i)=pick_smoo(ikx+int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)-&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)-int(pick_cos(ikx,iky,i))
               end if
            end if
            if (ikx-int(pick_sin(ikx,iky,i))>=1.and.ikx-int(pick_sin(ikx,iky,i))<=sog%xm%n) then
               if (pick_wavzt(ikx-int(pick_sin(ikx,iky,i)),iky,i).lt.&
                   -pick_smoo(ikx-int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)+&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)+int(pick_cos(ikx,iky,i))) then
                   pick_wavzt(ikx-int(pick_sin(ikx,iky,i)),iky,i)=-pick_smoo(ikx-int(pick_sin(ikx,iky,i)),iky,int(sog%shx%n/2.),i)+&
                   pick_smoo(ikx,iky,int(sog%shx%n/2.),i)+int(pick_cos(ikx,iky,i))
               end if
            end if
write(0,*) "passei 04",ikx,iky

         end if
      end do
   end do
write(0,*) "*****************************************************************"
   do ikx=1,sog%xm%n
      where(pick_smoo(ikx,:,int(sog%shx%n/2.),i).eq.9999.) 
          pick_wavzt(ikx,:,i)=wavz
          pick_wavzb(ikx,:,i)=wavz
      end where
   end do

end do
else
   if (taperon) then
      do i=1,wavz
         taper(i)=(1+cos(pi*(wavz-i+1)/wavz))/2
         taper(size(taper)-i+1)=taper(i)
      end do
   else
      taper=1
   end if
end if
write(0,*) "passei 05"
allocate(temp(sog%xm%n,sog%shx%n,sog%ym%n));
do i=1,nhorz
   do iz=1,sog%z%n
      temp=0.
      call sreed("in",intmp,4*sog%xm%n*sog%shx%n*sog%ym%n)
      do iky=1,sog%ym%n
         do ih=1,sog%shx%n
            do ikx=1,sog%xm%n
               if (dip) then
                  if (iz .ge. int(pick_smoo(ikx,iky,ih,i))-pick_wavzb(ikx,iky,i) .and. &
                      iz .le. int(pick_smoo(ikx,iky,ih,i))+pick_wavzt(ikx,iky,i)) &
                      temp(ikx,ih,iky)=intmp(ikx,ih,iky)
               else
                  if (iz .ge. int(pick_smoo(ikx,iky,ih,i))-wavz .and. iz .le. int(pick_smoo(ikx,iky,ih,i))+wavz) &
                     temp(ikx,ih,iky)=intmp(ikx,ih,iky)*taper(iz-int(pick(iky,ikx,ih,i))+wavz+1)**2
               end if
            end do
         end do
      end do
      call srite("out",temp,4*sog%xm%n*sog%shx%n*sog%ym%n)
   end do
end do

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

call sep_end_prog()
end program
