program OW_semblance
  use sep
  use cgstep_mod1
  use laplac2
  use solver_reg_mod1
  use owsemb

  implicit none
  logical                                         :: adj
  real, allocatable		  		:: in(:),out(:),curve(:,:)
  integer					:: n1,n2,n3,nrho,i1,i2,i3,irho,izd,stat,niters
  real                                            :: d1,d2,d3,o1,o2,o3,orho,drho,alfa1,alfa2,rho
  real                                            :: gama,sgama,cgama,zm,zd,eps,z1,z2,rhov
  call sep_init()
  call sep_begin_prog()
  call from_history("n1",n1);call from_history("d1",d1);call from_history("o1",o1)
  call from_history("n2",n2);call from_history("d2",d2);call from_history("o2",o2)
  call from_history("n3",n3);call from_history("d3",d3);call from_history("o3",o3)  
  call from_param("orho",orho,.8);call from_param("drho",drho,.016);call from_param("nrho",nrho,101)
  call from_param("adj",adj,.false.);call from_param("alfa1",alfa1,18.);call from_param("alfa2",alfa2,0.)
  call from_param("eps",eps,0.1);call from_param("niter",niters,50)
  call to_history("n1",n1);call to_history("d1",d1);call to_history("o1",o1)
  call to_history("n2",nrho);call to_history("d2",drho);call to_history("o2",orho)
  call to_history("n3",n3);call to_history("d3",d3);call to_history("o3",o3)  
  if (exist_file("curve")) then
     call to_history("n1",2,"curve");call to_history("d1",1.,"curve");call to_history("o1",1.,"curve")
     call to_history("n2",n2,"curve");call to_history("d2",d2,"curve");call to_history("o2",o2,"curve")
     call to_history("n3",1,"curve");call to_history("d3",1.,"curve");call to_history("o3",1.,"curve")  
     allocate(curve(2,n2))
     call from_param("z1",z1);  call from_param("z2",z2);  call from_param("rhov",rhov);
  end if  
  alfa1=cosd(alfa1);alfa2=cosd(alfa2)
  allocate(in(n1*n2*n3),out(n1*nrho*n3))
  call sreed("in",in,4*n1*n2*n3)
  stat=owsem_ini(n1,o1,d1,n2,o2,d2,n3,nrho,orho,drho,alfa1,alfa2)
  if (adj) then
     stat=owsemb_lop(.true.,.false.,out,in)
  else
     call laplac2_init(n1,nrho)
     call solver_reg(m=out, d=in,Fop=owsemb_lop,Aop=laplac2_lop,stepper=cgstep,nAop=n1*nrho*n3,niter=niters,eps=eps,verb=.true.)
  end if
  call srite("out",out,4*n1*nrho*n3)
  if (exist_file("curve")) then
     do i2=1,n2
        gama=(i2-1)*d2+o2;sgama=sind(gama)**2;cgama=cosd(gama)**2
        do i1=1,2
           if (i1 == 1) then
              zd=(1+(rhov-1)/alfa1*sgama/(alfa1**2-sgama))*z1
           else
              zd=(1+(rhov-1)/alfa2*sgama/(alfa2**2-sgama))*z2
           end if
           curve(i1,i2)=zd
        end do
     end do
     call srite("curve",curve,8*n2)
     deallocate(curve)
  end if
  deallocate(in,out)
  call sep_end_prog()
end program 
