!!
!! The routine calculates synthetic primaries from a dataset
!! previously submitted to any noise attenuation process. It
!! is based on the Aki-Richards AVO equation
!!
!!               r()=A + B[sin()]^2 + C[sin()*tan()]^2
!!
!!          
!!
!!

  program Drop

  use sep

  implicit none
  
  real		  		    		:: theta,amin, amax
  integer                         		:: i3,i2,ia,ib
  integer					:: n1,n2,n3,n4
  real                            		:: o1,o2,o3,o4,d1,d2,d3,d4
  complex,dimension(:,:,:),allocatable 		:: cmp
  
  call sep_init()
  call sep_begin_prog()

  !! get info from history file
    call from_history("n1",n1);   call from_history("n2",n2);   call from_history("n3",n3,1)
    call from_history("d1",d1);   call from_history("d2",d2);   call from_history("d3",d3,1.)
    call from_history("o1",o1);   call from_history("o2",o2);   call from_history("o3",o3,1.)
    call from_history("n4",n4);   call from_history("o4",o4);   call from_history("d4",d4,1.)

  !! get info from parameter
  call from_param("amin",amin,10.);   call from_param("amax",amax,30.);

  !! write info out to history file
  call to_history("n1",n1);call to_history("n2",n2);call to_history("n3",n3);call to_history("n4",n4);
  call to_history("d1",d1);call to_history("d2",d2);call to_history("d3",d3);call to_history("d4",d4);
  call to_history("o1",o1);call to_history("o2",o2);call to_history("o3",o3);call to_history("o4",o4);

!! allocate working space
  allocate(cmp(n1,n2,n3))
do i3=1,n4
   cmp=0.
   call sreed("in",cmp,n1*n2*n3*8)
   ia = int((amin-o1)/d1 + 1); ib = int((amax-o1)/d1 + 1);write (0,*) n4, i3
   cmp(ia:ib,:,:)=0.
!   write(0,*) maxval(cmp(:,int((ia+ib)/2)))
   call srite("out",cmp,n1*n2*n3*8)
end do
!! clean up
deallocate(cmp)
call sep_end_prog()
end program Drop
