program test
use sep
use reflect2d_mod
use taper_mod
use ricker_mod
use oneway_mod
implicit none

real    , dimension(:,:,:)      , pointer :: data
real    , dimension(:,:)        , pointer :: v, refl
real    , dimension(:)          , pointer :: sou
complex , dimension(:)          , pointer :: sou1,sou2
complex , dimension(:)          , pointer :: source,source1
real    , dimension(:)          , pointer :: x, z, xs, xr, t

integer         :: ix, iz, it, ixr, ixs, iv 
integer         :: nx, ny, nz
integer         :: nxr, nxs
integer         :: nt
integer         :: nref
integer         :: ixpad, imax_z, imin_z
integer         :: jz, jx, jxr, jxs
integer         :: spike
integer         :: esize
real            :: dx, dy, dz, ox, oy, oz  
real            :: dxr, dxs, oxr, oxs  
real            :: dt, ot  
real            :: pi  
real            :: tape 
integer         :: tape_opt 
real            :: freq
integer         :: restart
logical         :: flip, greens

character(len=256) :: vel_="vel", ref_="ref", sou_="sou", data_="dat",sou1_="sou1", vel1_="vel1_", wfld_d="wfld_d", wfld_u="wfld_u" 
tape=.1;
tape_opt=2; !1=triang. 2==hanning taper
pi=3.1415937; 
spike=0;
freq=20
nref=2;
restart=0

call initpar()
call sep_begin_prog()
call from_param("greens",greens,.false.) ;

write(*,*)greens

if (exist_file("vel"))then
                call from_aux("vel","n1",nz); call from_aux("vel","d1",dz); call from_aux("vel","o1",oz)
                call from_aux("vel","n2",nx); call from_aux("vel","d2",dx); call from_aux("vel","o2",ox)
                call from_aux("vel","n3",ny); call from_aux("vel","d3",dy); call from_aux("vel","o3",oy)
else
                call seperr('velocity file does not exist!!!!!!')
endif


allocate(v(nz,nx),refl(nz,nx),x(nx),z(nz))
call sreed("vel", v, 4*nx*ny*nz)

call from_param("nxs",nxs,10) ; if (mod(nxs,2) ==0.) nxs=nxs+1
call from_param("nxr",nxr,nx) ; if (mod(nxr,2) ==0.) nxr=nxr+1

if (nxr==1) then; dxr=1; else; dxr=floor(((nx-1)*dx)/(nxr-1)); endif;  oxr=ox; jxr=1;
if (nxs==1) then; dxs=1; else; dxs=floor(((nx-1)*dx)/(nxs-1)); endif;  oxs=ox; jxs=1;
        allocate(xr(nxr),xs(nxs))
nt=1501; dt=0.004; ot=0.;
        allocate(t(nt))
ixpad=ceiling(4.*nx); 
        do  ix=1, nx;  x( ix) =  ox + ( ix-1)* dx; enddo        
        do  iz=1, nz;  z( iz) =  oz + ( iz-1)* dz; enddo        
        do ixr=1,nxr; xr(ixr) = oxr + (ixr-1)*dxr; enddo        
        do ixs=1,nxs; xs(ixs) = oxs + (ixs-1)*dxs; enddo        
        do  it=1, nt;  t( it) =  ot + ( it-1)* dt; enddo        
imax_z=min(nz,floor(nz/2.)+floor(nz*.05));
imin_z=max(1,floor(nz/3.)-floor(nz*.05));

write(*,*)'imin_z,imax_z', imin_z,imax_z

!v(1:ceiling(nz/4.),:)=2000
!v(ceiling(nz/4.):nz,:)=2300
v(ceiling(nz*.72):nz,:)=3000
!v=3000
!v(:,1:floor(nx/2.))=2000
!v(:,ceiling(nx/2.):nx)=2400
!v=2000
call reflec2d(refl, v, nz, nx)
call taper(refl,refl, dx,nx,    0, tape, tape_opt);
call auxinout(ref_)
call auxputch("n1", "i", nz , ref_); call auxputch("n2", "i", nx , ref_); 
call auxputch("d1", "r", dz , ref_); call auxputch("d2", "r", dx , ref_)
call auxputch("o1", "r", oz , ref_); call auxputch("o2", "r", ox , ref_)
call auxputch("label1", "s", "z" , ref_)
call auxputch("label2", "s", "x" , ref_)
call auxputch("esize" , "i",  4  , ref_)
call srite("ref", refl, 4*nz*nx*ny)

!get the source wavelet
allocate(source1(nt),source(nt),sou(nt))
call ricker(source1,nt,freq,t(2)-t(1),floor(nt*dt),.true.);
source=0;
if(ceiling(nt/2.)-1==nt-ceiling(nt/2.))then; source(1:ceiling(nt/2.)   )=source1(ceiling(nt/2.)-10:nt);
else;                                        source(1:ceiling(nt/2.)+1.)=source1(ceiling(nt/2.)-10:nt);
endif

sou=source
write(*,*)'writing source wavelet.....'
call auxinout(sou_); 
call auxputch("n1    ","i", nt,sou_);call auxputch("d1   ","r",dt,sou_); call auxputch("o1", "r", ot , sou_); 
call auxputch("label1","s","t",sou_);call auxputch("esize","i", 4,sou_); call srite(sou_,sou , 4*nt)

allocate(data(nt,nxr,nxs))
call oneway(data,x, xs, xr, z, t, v, nref, source, ixpad , imin_z, imax_z, jz, jx, jxr, jxs, refl,restart,tape,tape_opt,greens)
write(*,*)'writing data.....'
call auxinout(data_)
call auxputch("n1", "i", nt, data_); call auxputch("n2", "i", nxr, data_); call auxputch("n3", "i", nxs, data_); 
call auxputch("d1", "r", dt, data_); call auxputch("d2", "r", dxr, data_); call auxputch("d3", "r", dxs, data_); 
call auxputch("o1", "r", ot, data_); call auxputch("o2", "r", oxr, data_); call auxputch("o3", "r", oxs, data_); 
call auxputch("label1", "s", "t" , data_)
call auxputch("label2", "s", "xr" , data_)
call auxputch("label3", "s", "xs" , data_)
call auxputch("esize" , "i",  4  , data_)
call srite("dat", data, 4*nt*nxr*nxs)

deallocate(refl, v, x, z, xs, xr, t, source, source1, sou, data)
call sep_end_prog()
end program
