!!
!! 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 mute_mod
  use cgstep_mod
  use adj_mod
!!  use solver_reg_mod
  use solver_tiny_mod
  use avomodel
!!  use igrad1

  implicit none
  
  real		  		    		:: theta, time, pi
  integer	  			    	:: n, npe, npi
  integer					:: niter, mube, mubi, mufe, mufi
  integer                         		:: i1,i2,i3,j1,j2
  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
  real,dimension(:),allocatable   		:: m,tripa

  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)

!  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)
  call from_param("niter",niter,3)
  call from_param("eps",eps,1.)

  !! 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))
  open( unit=2, file="report.txt", status="unknown")
!! allocate working space
  allocate(cdp(n1,n2),dado(n1,n2),m(3),tripa(n2))
!!
pi = 4*atan(1.)
do i3=1,n3
  cdp=0.;tripa=0;m=0.;dado=0.
  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
	if ((npe.ne.0).or.(npi.ne.0)) then
		call mute(n2,o2,d2,npe,npi,emx,emt,imx,imt,time,mube,mufe,mubi,mufi,tripa,cdp(i1,:))
		if (mod(i1,100).eq.1) write(2,*) 'dado prof =',time,mube,mufe,mubi,mufi,npe,npi
		call avomodel_init(n2,o2,d2,mube,mufe,mubi,mufi)
	else 
		tripa=cdp(i1,:)
		call avomodel_init(n2,o2,d2,mube,mufe,mubi,mufi)
	end if
!!	call solver_reg(m=m, d=tripa, &
!!	    Fop=avomodel_lop,Aop=igrad1,stepper=cgstep,nAop=n2,niter=niter,eps=eps,&
!!	    verb=.true.)
	call solver_tiny(m=m, d=tripa, Fop=avomodel_lop,stepper=cgstep,niter=niter)
!!  if (mod(i3,5).eq.1) write(2,*) i3,'    prof =',time,'m1 = ',m(1),'m2 = ',m(2),'m3 = ',m(3)
  !!
  !! Creates modeled data
  !!
  	do i2=1,n2
		theta=pi*(o2 + (i2-1)*d2)/180
		dado(i1,i2)=m(1)+m(2)*sin(theta)**2+m(3)*sin(theta)**2*tan(theta)**2
!!  if (mod(i1,100).eq.1) write(2,*) 'dado prof =',time,dado(i1,30),dado(i1,130),dado(i1,230)
	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 cgstep_close()
close (unit=2)
!! clean up
deallocate(cdp,dado,tripa,m)
call sep_end_prog()
end program AVO
