!!
!! 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[tan()]^2
!!
!!                   
!! It looks like LIFT, Paradigm Geophysical.
!!
!!
program LSQpr

  use sep

  implicit none
  
  real theta, pi, x(3)
  integer                         :: i1,i2,i3,n1,n2,n3
  integer			  :: n
  real                            :: o1,o2,o3,d1,d2,d3
  real				  :: s2,s4,t2,t4,s2t2,det
  real                            :: fang,lang	 
  real,dimension(:,:,:),allocatable :: cdp0
  real,dimension(:,:),allocatable :: dado, Q, cdp
  real,dimension(:), allocatable  :: m, d

  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)
  call from_history("d1",d1)
  call from_history("d2",d2)
  call from_history("d3",d3)
  call from_history("o1",o1)
  call from_history("o2",o2)
  call from_history("o3",o3)

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

  call sep_close()

!! get processing parameters
!!  call from_param("fang",fang)
!!  call from_param("lang",lang)


  pi=atan(1.)*4.
!!  do i2=1,n2
!!	  write(0,*)'Processing CDP',i2,'out of',n2  
!!	  call sreed_window("in",3,(/n1,n2,n3/),(/n1,1,n3/),(/0,i2-1,0/),(/1,1,1/),4,cdp)
!!	  allocate(cdp(n1,n2,n3))
	allocate(cdp0(n1,n2,n3))
!!	  call sreed("in",cdp0,n1*n2*n3*4)
	  call sep_read(cdp0,'in')
	allocate(cdp(n1,n2)) 
	do i1=1,n1
	do i2=1,n2
	cdp(i1,i2)=cdp0(i1,i2,1)
	end do 
	end do 
	do i1=1,n1
!!
!! Creates the inverted matrix
!!
	     s2=0;s4=0;t2=0;t4=0;s2t2=0;d(1)=0;d(2)=0;d(3)=0;
	     Q(1,2)=0;Q(1,3)=0;Q(2,1)=0;Q(2,2)=0;Q(2,3)=0;Q(3,1)=0;Q(3,2)=0;Q(3,3)=0;
	     do i3=1,n3
		theta=o3 + (i3-1)*d3
!!	        d(1)=d(1)+cdp(i1,i3)
!!	        d(2)=d(2)+cdp(i1,i3)*sin(theta*pi/180)**2
!!	        d(3)=d(3)+cdp(i1,i3)*tan(theta*pi/180)**2
		s2=s2+sin(theta*pi/180)**2
		t2=t2+tan(theta*pi/180)**2
		s4=s4+sin(theta*pi/180)**4
		t4=t4+tan(theta*pi/180)**4
		s2t2=s2t2+sin(theta*pi/180)**2*tan(theta*pi/180)**2
	     end do
	     det=n*s4*t4+2*s2*t2*s2t2-s4*t2**2-n*s2t2**2-t4*s2**2
             Q(1,1)=(s4*t4-(s2*t2)**2)/det
	     Q(1,2)=-((s2*t4-t2*s2t2))/det
	     Q(1,3)=(s2*s2t2-s4*t2)/det
	     Q(2,1)=(Q(1,2))/det
	     Q(2,2)=(n*t4-t2**2)/det
	     Q(2,3)=-((n*s2t2-s2*t2))/det
	     Q(3,1)=(Q(1,3))/det
	     Q(3,2)=(Q(2,3))/det
	     Q(3,3)=(n*s4-s2**2)/det
!!
!! Obtains A, B, C
!!
	     m=matmul(Q,d)
!!
!! Creates modeled data
!!
	     do i3=1,n3
		theta=o3 + (i3-1)*d3
		dado(i1,i3)=m(1)+m(2)*sin(theta*pi/180)**2+m(3)*tan(theta*pi/180)**2
	     end do
	  end do
!!          call srite_window("out",3,(/n1,n2,n3/),(/n1,1,n3/),(/0,i2-1,0/),(/1,1,1/),4,dado)
	  call srite("out",dado,n1*n3*4)
!!  end do
!! clean up
  deallocate(cdp,Q,d,m,dado)

  call sep_end_prog()
  end program LSQpr
