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

  use sep
  use mute_mod
  use cgstep_mod
  use adj_mod
  use solver_reg_mod
  use solver_prc_mod
  use solver_tiny_mod
  use avomodel
  use tcai1
  use igrad1

  implicit none
  
  real		  		    		:: theta, time, pi
  integer	  			    	:: n, npe, npi, nsolv
  integer					:: niter
  integer                         		:: i1,i2,i3,j1,j2,nAop
  integer					:: n1,n2,n3
  real                            		:: o1,o2,o3,d1,d2,d3,eps
  real				    		:: emt(10),emx(10),imt(10),imx(10)
  real,dimension(:,:),allocatable 		:: cdp,dado,mask,A,B,C
  real,dimension(:),allocatable   		:: m,mask1,data
  real,dimension(:),pointer			:: aa

  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)

  !! get info from parameter
  call from_param("niter",niter,50)
  call from_param("eps",eps,1.)
  call from_param("solv",nsolv,3)

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

!! allocate working space
  allocate(cdp(n1,n2),mask(n1,n2),data(n1*n2*n3),mask1(n1*n2*n3),aa(3))
!!
A=0.;B=0.;C=0.
nAop=size(m)+size(aa)

mask=1.;m=0.;
do i3=1,n3
  call sreed_window("in",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,cdp)
  do i1=1,n1
	time=o1+(i1-1)*d1
	if ((npe.ne.0).or.(npi.ne.0)) then
		call mute(n2,o2,d2,npe,npi,emx,emt,imx,imt,time,mask(i1,:))
	end if
     do i2=1,n2
	mask1((i1-1)+(i2-1)*n1+(i3-1)*n2+1)=mask(i1,i2)
	data((i1-1)+(i2-1)*n1+(i3-1)*n2+1)=cdp(i1,i2)
     end do
  end do
end do
deallocate(cdp,mask)
allocate(m(3*n1*n3))

call avomodel_init(n1,n2,n3,o2,d2,mask1)

if (nsolv.eq.1) then
  call solver_tiny(m=m, d=data, Fop=avomodel_lop,stepper=cgstep,niter=niter)
elseif (nsolv.eq.2) then
  call solver_reg(m=m, d=data, &
      Fop=avomodel_lop,Aop=igrad1_lop,stepper=cgstep,nAop=nAop,niter=niter,eps=eps,&
      verb=.false.)
else
  call solver_prc(m=m, d=data, &
      Fop=avomodel_lop,Sop=igrad1_lop,stepper=cgstep,nSop=nAop,niter=niter,eps=eps,&
      verb=.false.)
end if

  call avomodel_close()
  !!
  !! Creates modeled data
  !! 
  deallocate(data,mask1)
!!  allocate(A(n1,n3),B(n1,n3),C(n1,n3))
  allocate(dado(n1,n2))
  A(:,i3)=m(1:n1);B(:,i3)=m(n1+1:2*n1);C(:,i3)=m(2*n1+1:3*n1)
  do i1=1,n1
	A(i1,i3)=m(i1)
	B(i1,i3)=m(n1+i1)
	C(i1,i3)=m(2*n1+i1)
      do i2=1,n2
	 theta=o2 + (i2-1)*d2
	 dado(i1,i2)=m(i1)+m(n1+i1)*sind(theta)**2+m(2*n1+i1)*sind(theta)**2*tand(theta)**2
      end do
  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 ("A_out",A,n1*n3*4)
call srite ("B_out",B,n1*n3*4)
call srite ("C_out",C,n1*n3*4)
call cgstep_close()
close (unit=2)
!! clean up
deallocate(cdp,dado,data,mask,mask1,m)
call sep_end_prog()
end program AVO
