!!
!! 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
  
  character*20			    :: date, tempo
  real		  		    :: theta, pi, time, rnc
  integer	  		    :: n, npe, npi
  integer	  		    :: if2, if3, k
  real		  		    :: s2,s4,t2,t4,s2t2,det
  integer                           :: i1,i2,i3,n1,n2,n3
  real                              :: o1,o2,o3,d1,d2,d3
  real				    :: emt(10),emx(10),imt(10),imx(10)
  real,dimension(:,:),allocatable :: cdp,dado,A,B,C
  real,dimension(:,:),allocatable   :: Q
  real,dimension(:),allocatable     :: m,d
!  integer,dimension(:), allocatable :: iloc
  integer			    :: iloc(2),il
  integer, external :: sreed

  call sep_init()
  call sep_begin_prog()

  call date_and_time (date,tempo)
  write(0,*) "Begin"
  write(0,*) "Date = ",date," Time = ",tempo


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

!  call from_param("n1",n1)
!  call from_param("n2",n2)
!  call from_param("n3",n3)
!  call from_param("d1",d1)
!  call from_param("d2",d2)
!  call from_param("d3",d3)
!  call from_param("o1",o1)
!  call from_param("o2",o2)
!  call from_param("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)

!!  get mute parameters
  call from_param("npe",npe) 
  call from_param("emt1",emt(1))
  call from_param("emt2",emt(2))
  call from_param("emt3",emt(3))
  call from_param("emt4",emt(4))
  call from_param("emt5",emt(5))
  call from_param("emt6",emt(6))
  call from_param("emt7",emt(7))
  call from_param("emt8",emt(8))
  call from_param("emt9",emt(9))
  call from_param("emt10",emt(10))
  call from_param("emx1",emx(1))
  call from_param("emx2",emx(2))
  call from_param("emx3",emx(3))
  call from_param("emx4",emx(4))
  call from_param("emx5",emx(5))
  call from_param("emx6",emx(6))
  call from_param("emx7",emx(7))
  call from_param("emx8",emx(8))
  call from_param("emx9",emx(9))
  call from_param("emx10",emx(10))
  call from_param("npi",npi) 
  call from_param("imt1",imt(1))
  call from_param("imt2",imt(2))
  call from_param("imt3",imt(3))
  call from_param("imt4",imt(4))
  call from_param("imt5",imt(5))
  call from_param("imt6",imt(6))
  call from_param("imt7",imt(7))
  call from_param("imt8",imt(8))
  call from_param("imt9",imt(9))
  call from_param("imt10",imt(10))
  call from_param("imx1",imx(1))
  call from_param("imx2",imx(2))
  call from_param("imx3",imx(3))
  call from_param("imx4",imx(4))
  call from_param("imx5",imx(5))
  call from_param("imx6",imx(6))
  call from_param("imx7",imx(7))
  call from_param("imx8",imx(8))
  call from_param("imx9",imx(9))
  call from_param("imx10",imx(10))

  if (exist_file("A_out")) then
        call to_history('n1',n1,"A_out")
        call to_history('d1',d1,"A_out")
        call to_history('o1',o1,"A_out")
        call to_history('n2',n3,"A_out")
        call to_history('d2',d3,"A_out")
        call to_history('o2',o3,"A_out")
  end if
  if (exist_file("B_out")) then
        call to_history('n1',n1,"B_out")
        call to_history('d1',d1,"B_out")
        call to_history('o1',o1,"B_out")
        call to_history('n2',n3,"B_out")
        call to_history('d2',d3,"B_out")
        call to_history('o2',o3,"B_out")
  end if
  if (exist_file("C_out")) then
        call to_history('n1',n1,"C_out")
        call to_history('d1',d1,"C_out")
        call to_history('o1',o1,"C_out")
        call to_history('n2',n3,"C_out")
        call to_history('d2',d3,"C_out")
        call to_history('o2',o3,"C_out")
  end if

  pi=atan(1.)*4.

!! allocate working space
  allocate(cdp(n1,n2),dado(n1,n2),Q(3,3),d(3),m(3),A(n1,n3),B(n1,n3),C(n1,n3))

!  call sep_read(cdp,"in")
!  if(size(cdp)*4/=sreed("in",cdp,size(cdp)*4)) call seperr("trouble reading in data")
!  do i1=1,n1
!    write(0,*) "CHECK THE SUM",sum(cdp(i1,128))
!    write(0,*) "CHECK THE VALUE",cdp(i1,128)
!  end do

  pi=atan(1.)*4.
  
!!
!! 
!!
do i3=1,n3
  call sreed_window("in",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,cdp)
  if (mod(i3,1).eq.0) write(0,*) "CMP  ", i3, " of ", n3
  do i1=1,n1
	time=o1+(i1-1)*d1
	if2=0;if3=0
	!!
	!! External mute
	!!
	if (npe.ne.0) then
		if (time <= emt(1)) then
			if2=n2-int((emx(1)-o2)/d2)+1
		elseif(time >= emt(npe)) then 
			if2=n2-int((emx(npe)-o2)/d2)+1
		else
			iloc=minloc(emt, dim=1,MASK=emt>time);il=iloc(1)
			rnc=(emx(il)-emx(il-1))/(emt(il)-emt(il-1))
			if2=n2-int(((rnc*(time-emt(il))+emx(il))-o2)/d2)+1
		end if
	else
		if2=1
	end if

        
	!!
	!! Internal mute
	!!
	if (npi.ne.0) then
		if (time >= imt(npi)) then
			if3=int(imx(npi)/d2)
		elseif((time >= imt(1)).and.(time < imt(npi))) then 
			iloc=minloc(imt, dim=1,MASK=imt<time);il=iloc(1)
			rnc=(imx(il)-imx(il-1))/(imt(il)-imt(il-1))
			if3=int((rnc*(time-imt(il))+imx(il))/d2)
		end if
	end if
!   if (mod(i1,3).eq.1)   write(0,*) "t ", i1, " if2 = ", if2, " if3 = ", if3
!  if (mod(i1,10).eq.1) write(0,*) "t ", i1, " if2 = ", if2, " if3 = ", if3
!  if (mod(i1,5).eq.1) write(0,*) "t ", i1, " if2 = ", if2
!  if (mod(i1,10).eq.1) write(0,*) "t ", i1, " if3 = ", if3

!	mube=int(n2/2)-if2+int(5/d2)+1; mufe=int(n2/2)+if2-int(5/d2)+1
!	mubi=int(n2/2)-if3+1; mufi=int(n2/2)+if3+1
!
!	do i2=1,n2
!		if ((i2<mube).or.(i2>mufe)) then
!			dado(i1,i2)=0.0
!		elseif ((i2>=mubi).and.(i2<=mufi).and.(time>imt(1))) then 
!			dado(i1,i2)=0.0
!		elseif ((i2>=mubi).and.(i2<=mufi).and.(time<=imt(1))) then 
!			dado(i1,i2)=cdp(i1,i2)
!		else
!			dado(i1,i2)=cdp(i1,i2)
!		end if
!	end do
!  end do

	!! 
	!! Creates the inverted matrix
	!!
	s2=0;s4=0;t2=0;t4=0;s2t2=0;d(1)=0;d(2)=0;d(3)=0;k=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;
	
	if ((npi.ne.0).and.(imt(1).ne.0.).and.(time > imt(1))) then
		do i2=if2,int(n2/2)-if3
			theta=o2 + (i2-1)*d2
			d(1)=d(1)+cdp(i1,i2)	
			d(2)=d(2)+cdp(i1,i2)*sind(theta)**2
			d(3)=d(3)+cdp(i1,i2)*tand(theta)**2
			s2=s2+sind(theta)**2
			t2=t2+tand(theta)**2
			s4=s4+sind(theta)**4
			t4=t4+tand(theta)**4
		 	s2t2=s2t2+sind(theta)**2*tand(theta)**2
			k=k+1
		end do
		do i2=int(n2/2)+if3,n2-if2+1
			theta=o2 + (i2-1)*d2
			d(1)=d(1)+cdp(i1,i2)	
			d(2)=d(2)+cdp(i1,i2)*sind(theta)**2
			d(3)=d(3)+cdp(i1,i2)*tand(theta)**2
			s2=s2+sind(theta)**2
			t2=t2+tand(theta)**2
			s4=s4+sind(theta)**4
			t4=t4+tand(theta)**4
		 	s2t2=s2t2+sind(theta)**2*tand(theta)**2
			k=k+1
		end do
	else
		do i2=if2,n2-if2+1
			theta=o2 + (i2-1)*d2
			d(1)=d(1)+cdp(i1,i2)	
			d(2)=d(2)+cdp(i1,i2)*sind(theta)**2
			d(3)=d(3)+cdp(i1,i2)*tand(theta)**2
			s2=s2+sind(theta)**2
			t2=t2+tand(theta)**2
			s4=s4+sind(theta)**4
			t4=t4+tand(theta)**4
		 	s2t2=s2t2+sind(theta)**2*tand(theta)**2
			k=k+1
		end do
	end if
	det=k*s4*t4+2*s2*t2*s2t2-s4*t2**2-k*s2t2**2-t4*s2**2
	Q(1,1)=(s4*t4-(s2t2)**2)/det
	Q(1,2)=-(s2*t4-t2*s2t2)/det
	Q(1,3)=(s2*s2t2-s4*t2)/det
	Q(2,1)=Q(1,2)
	Q(2,2)=(k*t4-t2**2)/det
	Q(2,3)=-(k*s2t2-s2*t2)/det
	Q(3,1)=Q(1,3)
	Q(3,2)=Q(2,3)
	Q(3,3)=(k*s4-s2**2)/det

	!!
	!! Obtains A, B, C
	!!

	m=matmul(Q,d)
!	if ((mod(i1,10).eq.1).and.(mod(i3,50).eq.1)) write(0,*) "p ", i1*d1, "m1= ",m(1),"m2=",m(2),"m3=",m(3)
	!!
	!! Creates modeled data
	!!

	do i2=if2,n2-if2+1
		theta=o2 + (i2-1)*d2
		dado(i1,i2)=(m(1)+m(2)*sind(theta)**2+m(3)*tand(theta)**2)
	end do
!  	if (i1.eq.451) write(0,*) dado(i1,:)
	A(i1,i3)=m(1)
	B(i1,i3)=m(2)
	C(i1,i3)=m(3)
  end do
  call srite_window("out",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,dado)
end do

!  call srite("out",dado,n1*n2*n3*4)
  call srite("A_out",A,n1*n3*4)
  call srite("B_out",B,n1*n3*4)
  call srite("C_out",C,n1*n3*4)

  call date_and_time (date,tempo)
  write(0,*) "End"
  write(0,*) "Date = ",date," Time = ",tempo

!! clean up
  deallocate(cdp,dado,Q,m,d,A,B,C)

  call sep_end_prog()
  end program LSQpr
