!!
!! 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		  		    		:: cmpx,ax
  integer                         		:: i3,i2,ia,ib
  integer					:: n1,n2,n3,n4
  real                            		:: o1,o2,o3,o4,d1,d2,d3,d4
!  complex,dimension(:,:,:),allocatable 		:: cmp
  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("cmp_x",cmpx,2000.)

  !! 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))
  allocate(cmp(n1,n2))
open(unit=2,file="diag.txt")
!do i3=1,n4
do i3=1,n3
   cmp=0.
!   call sreed("in",cmp,n1*n2*n3*8)
   call sreed("in",cmp,n1*n2*4)
!   do i2=1,n1
   do i2=1,n2
!      ax = ((i2-1)*d1+o1)/2+(i3-1)*d4+o4
      ax = ((i2-1)*d2+o2)/2+(i3-1)*d3+o3
      if (ax.eq.cmpx) then
	write(2,*) i3,i2,ax
      else
!	cmp(i2,:,:)=0
	cmp(:,i2)=0
      end if
   end do
   call srite("out",cmp,n1*n2*4)
end do
!! clean up
close(2)
deallocate(cmp)
call sep_end_prog()
end program Drop
