module planemod
use sep
implicit none


  real,allocatable,save,public:: head(:,:),tmp1(:),tmp2(:)
  real,allocatable,save,public:: picks(:),picksout(:)
  real,save,public:: sx,sy,gz,pi,ox,oy,os,oa,od,oz,maxx,maxy,maxs,maxa,maxd,maxz,mintau,tau0,t0out
  integer,save,public:: ntr,norm,nretain


contains
!----------------------------------------------------------------
  subroutine planemod_init(parmin,parmax,nposibl)
    real,dimension(:),intent(inout)::parmin
    real,dimension(:),intent(inout)   ::parmax
    integer,dimension(:),intent(in)::nposibl
    integer:: n1,nh
    real::o1,d1,gzi
    call sep_init(SOURCE)
    call from_aux("geom",'n1',nh)
    call from_aux("geom",'n2',ntr)
    call from_history("n1",n1);call from_history('d1',d1);call from_history("o1",o1)
    if(ntr.ne.n1) call seperr("number headers not matched to input")
    call to_history('n1',n1);call to_history('o1',o1);call to_history('d1',d1)
    if(nh==2) then
       call from_param('gz',gzi,68.)
       gz=gzi
    end if
    call from_param("norm",norm,2)
    call to_history("esize",4,"model")
    call to_history("n1",size(parmin)+1,"model")
    call to_history("label1","generations","model")
    call to_history("label2","xyzvadF","model")

    pi=acos(-1.)
    allocate(picks(n1),picksout(n1),head(nh,ntr))
    head=0.;picks=0.;picksout=0.

    call sep_read(picks)
    call sep_read(head,"geom")

    sx=head(2,1)
    sy=head(1,1)
    write(0,*)'  Location of correlated shotpoint (sx,sy):',sx,sy
    ox=parmin(1)+sx 
    oy=parmin(2)+sy
    oz=parmin(3)+gz  !! if gz is commented, z=tau_0 
    os=parmin(4)
    oa=parmin(5)
    od=parmin(6)

    parmin(1)=ox
    parmin(2)=oy
    parmin(3)=oz
    parmin(5)=parmin(5)*pi/180.
    parmin(6)=parmin(6)*pi/180.


    maxx=parmax(1)+sx
    maxy=parmax(2)+sy
    maxz=parmax(3)+gz
    maxs=parmax(4)
    maxa=parmax(5)
    maxd=parmax(6)


    parmax(1)=maxx
    parmax(2)=maxy
    parmax(3)=maxz
    parmax(5)=parmax(5)*pi/180.
    parmax(6)=parmax(6)*pi/180.
    
  end subroutine planemod_init
!----------------------------------------------------------------  
  subroutine func(npopsiz,nparam,j,funcval,parent,flag)
    integer, intent(in) :: npopsiz
    integer, intent(in) :: nparam
    integer, intent(in) :: j
    real, intent(out)   :: funcval
    integer, intent(in), optional :: flag
    real, dimension(nparam,npopsiz), intent(in) :: parent
    integer::st,count
    save count
    
    st=planemod_op(parent,j)
    call normmeasure(picks,picksout,funcval)    
    count=count+1

  end subroutine func
!----------------------------------------------------------------
!!page 144 "lessons in seismic computing",Slotnick
  integer function planemod_op(modl,j)
    real,dimension(:,:),intent(in)::modl
    integer,intent(in):: j
    real:: s,x,y,lx,ly,lt,la,ll,x2,s2,z2,t0,a,d,z,t,az
    integer:: i

    planemod_op=1
    picksout=0.
    !! current model vector
    x=modl(1,j)
    y=modl(2,j)
    z=modl(3,j)
    s=modl(4,j)
    az=modl(5,j)
    d=modl(6,j)

    !! remove azimuth from trial source loc to receiver 

    lx=x-sx
    ly=y-sy
    la=atan2(y-sy,x-sx)
    a=la-az
    if(a.lt.0.)a=pi+a
    a=cos(a)
    d=sin(d)
    
    !! distances
    lx=(x-sx)**2
    ly=(y-sy)**2
    x2=lx+ly    
    lx=sqrt(x2)
    z2=z*z
    !! vel squared
    s2=s*s
    !! term 1
!    lt=4*s2*z2+s2*x2
    lt=s2*z2+s2*x2
    !! term 2
!    ll=4*z*lx*d*a*s2
    !! event at source trace (correlated gathers) always at t=0
!    t0=sqrt(lt+ll)
    t0=sqrt(lt)

    trc:do i=1,ntr
       la=atan2(y-head(1,i),x-head(2,i))
       a=la-az
       if(a.lt.0.)a=pi+a
       a=cos(a)
       lx=(head(2,i)-x)**2
       ly=(head(1,i)-y)**2
       x2=lx+ly
       lx=sqrt(x2)
       ll=4*z*lx*d*a*s2
       lt=4*s2*z2+s2*x2
!       lt=s2*x2
       picksout(i)=sqrt(lt+ll)-t0
!       picksout(i)=sqrt(lt)!-t0
    end do trc

    tau0=2*z*s
    mintau=tau0*cos(asin(a*d))
    t0out=t0
    planemod_op=0

  end function planemod_op

!----------------------------------------------------------------
  subroutine normmeasure(dat,moddat,dist)
    real,dimension(:),intent(in):: dat,moddat
    real,intent(out)::dist
    integer::i1
    
    dist=0.
    if (norm==2) then

       do i1=1,ntr
          dist=dist+(dat(i1)-moddat(i1))**2
       end do
       dist=sqrt(dist)

    elseif(norm==1) then

       do i1=1,ntr
          dist=dist+abs(dat(i1)-moddat(i1))
       end do
    end if

  end subroutine normmeasure
!!----------------------------------------------------------------
subroutine exh_search_init(shx,shy,n1,n)
  integer:: n1,nh,n
  real::o1,d1,shx,shy


  call from_aux("geom",'n1',nh)
  call from_aux("geom",'n2',ntr)
  call from_history("n1",n1);call from_history('d1',d1);call from_history("o1",o1)
  
  call from_param("norm",norm,2)
  call from_param("nkeep",n,5); nretain=n ; write(0,*) 'nretain',nretain
  call to_history('n1',n1);call to_history('o1',o1);call to_history('d1',d1)
  if(ntr.ne.n1) call seperr("number headers not matched to input")
  allocate(tmp2(n))
  call sep_close()
  allocate(picksout(n1),head(nh,ntr))
  head=0.
  call sep_read(head,"geom")

  sx=head(2,1)
  sy=head(1,1)
  shx=sx
  shy=sy
end subroutine exh_search_init
subroutine fwdmod_init(shx,shy,n)
  integer:: n1,nh,n
  real::o1,d1,shx,shy,o2
  call from_history('n1',nh);call from_history('n2',n)
  call from_history('o2',o2)
  call to_history('n1',n);call to_history('o1',o2);call to_history('d1',1.)
  call to_history("n2",1);call to_history('o2',0.);call to_history('d1',1.)
  call sep_close()
  ntr=n
  allocate(picksout(ntr),head(nh,ntr))
  head=0.
  call sep_read(head)
  sx=head(2,1)
  sy=head(1,1)
  shx=sx
  shy=sy
end subroutine fwdmod_init
subroutine fwdmod_write()
  call sep_write(picksout)
end subroutine fwdmod_write

!----------------------------------------------------------------
subroutine getpicks(t)
real::t(:)
t=picksout
end subroutine getpicks
!----------------------------------------------------------------
subroutine rankupdates(modl,mods,res,updates)
  real::modl(:,:),mods(:,:),res
  integer::i,updates,m
  if (updates.le.nretain) then
     m=nretain +1 - updates
     mods(size(modl,1)+1,m)=res
     modl(1:4,1)=mods(1:4,m)
  else
     tmp2=mods(size(modl,1)+1,:)
     do i=1,nretain
        if(res.le.tmp2(i)) exit
     end do
     do m=nretain,i+1,-1
        mods(:,m)=mods(:,m-1)
     end do
     mods(1:size(modl,1),i)=modl(:,1)
     mods(size(modl,1)+1,i)=res
  end if
end subroutine rankupdates
!----------------------------------------------------------------
subroutine reportdone(mods)
  real:: mods(:),xout,yout,zout,sout,aout,dout,t0

  xout=mods(1)
  yout=mods(2)
  zout=mods(3)
  sout=mods(4)
  aout=mods(5)*180./pi
  dout=mods(6)*180./pi
  write(0,*) '**************************************************' 
  write(0,*) ' Model space input parameter ranges '
  write(0,*)' X  :',ox,' to ',maxx,' answer:',xout
  write(0,*)'   relX :',ox-sx,' to ',maxx-sx,' answer:',xout-sx
  write(0,*)' Y  :',oy,' to ',maxy,' answer:',yout
  write(0,*)'   relY :',oy-sy,' to ',maxy-sy,' answer:',yout-sy
  write(0,*)' Horizontal distance from shot trace',sqrt((yout-sy)**2+(xout-sx)**2)
  write(0,*)' Z  :',oz,' to ',maxz,' answer:',zout
  write(0,*)'   relZ  :',oz-gz,' to ',maxz-gz,' answer:',zout-gz
  write(0,*)' S  :',os,' to ',maxs,' answer:',sout
  write(0,*)'       V :',1/os,' to ',1/maxs,' answer:',1/sout
  write(0,'(a,i3,a,i3,a,i3,a,f5.3,a)')' Az :',nint(oa),' to ',nint(maxa),' answer:',nint(aout),' dgs (',mods(5),' radians)'
  write(0,'(a,i3,a,i3,a,i3,a,f5.3,a)')' Dip:',nint(od),' to ',nint(maxd),' answer:',nint(dout),' dgs (',mods(6),' radians)'
  write(0,*)' Minimum time ',mintau
  write(0,*)' Time at source location',tau0
  write(0,*)' Time at shot trace ',t0out
  write(0,*)'  '
  write(0,*)'  '
end subroutine reportdone
!----------------------------------------------------------------
  integer function cmp_planemod_op(modl,j)
    real,dimension(:,:),intent(in)::modl
    integer,intent(in):: j
    real:: s,x,y,lx,ly,lt,la,x2,s2,t0,a,d,z,t,az
    integer:: i


    cmp_planemod_op=1
    picksout=0.
    !! current model vector
    x=modl(1,j)
    y=modl(2,j)
    z=modl(3,j)
    s=modl(4,j)
    az=modl(5,j)
    d=modl(6,j)

    !! remove azimuth from trial source loc to receiver 
    la=atan2(y-sy,x-sx)
    a=az-la  
    !! distances
    lx=(x-sx)**2
    ly=(y-sy)**2
    x2=lx+ly
    !! tau_0
!    lt=(s*z)**2
    lt=z**2
    !! azimuth/dip dependant slowness
    s2=s*s*(1-sin(d)*sin(d)*cos(a)*cos(a))
    !! event at source trace (correlated gathers) always at t=0
    t0=sqrt(lt+x2*s2)

    trc:do i=1,ntr
       la=atan2(y-head(1,i),x-head(2,i))
       a=az-la

       lx=(head(2,i)-x)**2
       ly=(head(1,i)-y)**2
       x2=lx+ly

       s2=s*s*(1-sin(d)*sin(d)*cos(a)*cos(a))
!       lt=(s*z)**2
       lt=(z)**2

       picksout(i)=sqrt(lt+x2*s2)-t0
    end do trc
    cmp_planemod_op=0
  end function cmp_planemod_op
!----------------------------------------------------------------
end module planemod
