program encode

use sep 
use OW_type
use OW_rndseed

implicit none

type (dat_grid)					:: rec, sou
type (mod_grid)					:: img
real						:: pi,xmin,xmax
integer						:: ixs,isp,i,stat,igap,n1,rxs,sxs,n_enc_sht,jgap,ngap
real, allocatable					:: phase(:,:)
complex, allocatable				:: rec1(:,:),sou1(:)
complex, allocatable				:: encrec(:,:,:),encsou(:,:,:),soutemp(:)

 call sep_init()
 call sep_begin_prog()

!---------------------------------------
! IO
!
 call from_history("n1",rec%h%n);	call from_history("d1",rec%h%d);	call from_history("o1",rec%h%o)
 call from_history("n4",rec%x%n);	call from_history("d4",rec%x%d);	call from_history("o4",rec%x%o)  
 call from_aux("sou","n3",sou%w%n);	call from_aux("sou","d3",sou%w%d);	call from_aux("sou","o3",sou%w%o)
 call from_param("xmin",xmin) 
 call from_param("xmax",xmax)
 call from_param("ngap",ngap,1)
 xmin=int(xmin/rec%h%d)*rec%h%d
 xmax=int(xmax/rec%h%d)*rec%h%d
 img%xm%n=(xmax-xmin)/rec%h%d+1.5
 img%xm%d=rec%h%d
 img%xm%o=xmin

 call to_history("n4",ngap,"encrec"); call to_history("d4",1.,"encrec"); call to_history("o4",1.,"encrec")  
 call to_history("n1",img%xm%n,"encrec"); call to_history("d1",img%xm%d,"encrec"); call to_history("o1",img%xm%o,"encrec")
 call to_history("n2",1,"encrec"); call to_history("d2",0.,"encrec"); call to_history("o2",0.,"encrec")
 call to_history("n3",sou%w%n,"encrec") ; call to_history("d3",sou%w%d,"encrec") ; call to_history("o3",sou%w%o,"encrec")
 call to_history("esize",8,"encrec")

 call to_history("n4",ngap,"encsou"); call to_history("d4",1.,"encsou"); call to_history("o4",1.,"encsou")  
 call to_history("n1",img%xm%n,"encsou"); call to_history("d1",img%xm%d,"encsou"); call to_history("o1",img%xm%o,"encsou")
 call to_history("n2",1,"encsou"); call to_history("d2",0.,"encsou"); call to_history("o2",0.,"encsou")
 call to_history("n3",sou%w%n,"encsou") ; call to_history("d3",sou%w%d,"encsou") ; call to_history("o3",sou%w%o,"encsou")
 call to_history("esize",8,"encsou")
!----------------------------------------------------
! Memory Allocation
!
 allocate(rec1(rec%h%n,sou%w%n),sou1(sou%w%n))

!--------------------------------------------------
! Reading source fctn 
 call sreed("sou",sou1,sou%w%n*8)

!--------------------------------------------------
!

!--------------------------------------------------
! Encoding shots
!
!
n_enc_sht=int(rec%x%n/ngap*1.)
pi=4*atan(1.)

allocate(phase(rec%x%n,sou%w%n))
call init_random_seed(rec%x%n,sou%w%n,phase)

allocate(encrec(img%xm%n,sou%w%n,ngap),encsou(img%xm%n,sou%w%n,ngap))
allocate(soutemp(sou%w%n))
encrec=0;encsou=0

do ixs=1,rec%x%n
   soutemp=0;rec1=0
   igap=ixs-int((ixs-1)/ngap)*ngap
   jgap=int((ixs-1)/ngap*1.)+1

   rxs = ceiling(((ixs-1)*rec%x%d+rec%x%o+rec%h%o-img%xm%o)/img%xm%d)
   sxs = ceiling(((ixs-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)

   call sreed("in",rec1,8*rec%h%n*sou%w%n)
   do i=1,rec%h%n
      rec1(i,:)=rec1(i,:)*exp(cmplx(0,2*pi*phase(ixs,:)))
      soutemp=sou1*exp(cmplx(0,2*pi*phase(ixs,:)))
   end do
   encrec(rxs:rxs+rec%h%n-1,:,igap)=encrec(rxs:rxs+rec%h%n-1,:,igap)+rec1
   encsou(sxs,:,igap)=encsou(sxs,:,igap)+soutemp
end do

call srite("encrec",encrec,8*img%xm%n*sou%w%n*ngap)
call srite("encsou",encsou,8*img%xm%n*sou%w%n*ngap)
deallocate(rec1,soutemp,sou1)

call sep_end_prog()
end program encode
