!!
!! 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 AVO

  use sep
  use tcai1_mod

  implicit none
  
  real		  		    		:: theta, pi
  integer                         		:: i2,i3,iz,stat
  integer					:: n1,n2,n3,A,B,C,z
  real                            		:: o1,o2,o3,d1,d2,d3
  real,dimension(:,:),allocatable 		:: rc,dado
  real,dimension(:),allocatable			:: w
  
  allocate(w(31))
  w(1)=105;w(2)=544;w(3)=2156;w(4)=6888;w(5)=18292;w(6)=41120;w(7)=79002;w(8)=129948;w(9)=181506;w(10)=209664;w(11)=185640
  w(12)=90624;w(13)=-68184;w(14)=-251328;w(15)=-398517;w(16)=-454950;w(17)=-398517;w(18)=-251328;w(19)=-68184;w(20)=90624;w(21)=185640
  w(22)=209664;w(23)=181506;w(24)=129948;w(25)=79002;w(26)=41120;w(27)=18292;w(28)=6888;w(29)=2156;w(30)=544;w(31)=105

  call sep_init()
  call sep_begin_prog()

  if(exist_file("indata")) then
  !! 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.)
  else
  !! get info from parameters
    call from_param("n1",n1);   call from_param("n2",n2);   call from_param("n3",n3,1)
    call from_param("d1",d1);   call from_param("d2",d2);   call from_param("d3",d3,1.)
    call from_param("o1",o1);   call from_param("o2",o2);   call from_param("o3",o3,1.)
  end if

  !! get info from parameter
  call from_param("A",A);   call from_param("B",B);   call from_param("C",C)
  call from_param("depth",z)

  !! write info out to history file
  call to_history("n1",n1);call to_history("n2",n2);call to_history("n3",1);call to_history("n4",n3);call to_history("n5",1);call to_history("n6",1)
  call to_history("d1",d1);call to_history("d2",d2);call to_history("d3",1.);call to_history("d4",d3);call to_history("d5",1.);call to_history("d6",1.)
  call to_history("o1",o1);call to_history("o2",o2);call to_history("o3",1.);call to_history("o4",o3);call to_history("o5",1.);call to_history("o6",1.)
  call to_history ("label1","z");call to_history ("label2","ax");call to_history ("label3","ay");call to_history ("label4","cmpx");
  call to_history ("label5","cmpy");call to_history ("label6","ro");

!! allocate working space
  allocate(rc(n1,n2),dado(size(w)+n1,n2))

pi = 4*atan(1.)
iz = z/d1 + 1.5
if (iz .gt. n1) call seperr("Reflector is deeper than the maximum depth")
call tcai1_init (size(w),w)
do i3=1,n3
  rc=0.;dado=0.
  if (mod(i3,10).eq.0) write(0,*) "CMP  ", i3, " of ", n3
  do i2=1,n2
     theta=(i2-1)*d2+o2 
     rc(iz,i2)=A+B*sind(theta)**2+C*tand(theta)**2
     stat = tcai1_lop ( .false., .true., rc(:,i2), dado(:,i2))
  end do
end do
rc=0.;rc=dado(int(size(w)/2+1.5):int(size(w)/2+1.5)+n1,:)
call sep_write(rc)
!! clean up
deallocate(rc,dado)
call sep_end_prog()
end program AVO
