!!
!! 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.
!!
!!
  subroutine lsqpr(n1,n2,o2,d2,cdp,dado)

  implicit none
  
  real     theta, pi
  integer  i1,i2,n1,n2
  integer  n
  real     o2,d2
  real	   s2,s4,t2,t4,s2t2,det
  real     fang,lang	 
  real	   cdp(n1,n2), dado(n1,n2), Q(3,3)
  real     m(3), d(3)

  pi=atan(1.)*4.
  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 i2=1,n2
		theta=o2 + (i2-1)*d2
		d(1)=d(1)+cdp(i1,i2)	
		d(2)=d(2)+cdp(i1,i2)*sin(theta*pi/180)**2
		d(3)=d(3)+cdp(i1,i2)*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=n2*s4*t4+2*s2*t2*s2t2-s4*t2**2-n2*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)=(n2*t4-t2**2)/det
	Q(2,3)=-((n2*s2t2-s2*t2))/det
	Q(3,1)=(Q(1,3))/det
	Q(3,2)=(Q(2,3))/det
	Q(3,3)=(n2*s4-s2**2)/det
!!
!! Obtains A, B, C
!!
	m=matmul(Q,d)
!!
!! Creates modeled data
!!
	do i2=1,n2
		theta=o2 + (i2-1)*d2
		dado(i1,i2)=m(1)+m(2)*sin(theta*pi/180)**2+m(3)*tan(theta*pi/180)**2
	end do
  end do
  return; end
!!	  do i1=1,n1
!!	  write(0,*)'Processing CDP',i3,'out of',n3  
!!	  call sreed_window("in",3,(/n1,n2,n3/),(/n1,1,n3/),(/0,i2-1,0/),(/1,1,1/),4,cdp)
!!	  call sreed("in",cdp,n1*n2*4)
!!
!!
