module sourcemod
use sep
implicit none


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


contains
!----------------------------------------------------------------
  subroutine sourcemod_init(parmin,parmax,nposibl)
    real,dimension(:),intent(inout)::parmin
    real,dimension(:),intent(inout)   ::parmax
    integer,dimension(:),intent(in)::nposibl
    integer:: n1,nx,ny,nz,ns,nh
    real::o1,d1,t0
    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 from_param("norm",norm,1)

    call to_history('n1',n1);call to_history('o1',o1);call to_history('d1',d1)


    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.
    if(nh.eq.2) gz=68.

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

    sx=head(2,1)
    sy=head(1,1)

    ox=sx + parmin(1)
    oy=sy + parmin(2)
    oz=gz + parmin(3)
    os=parmin(4)
    parmin(1)=ox
    parmin(2)=oy
    parmin(3)=oz

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

    nx=nposibl(1)
    ny=nposibl(2)
    nz=nposibl(3)
    ns=nposibl(4)

  end subroutine sourcemod_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=sourcemod_op(parent,j)
    call normmeasure(picks,picksout,funcval)    
    count=count+1

  end subroutine func
!----------------------------------------------------------------
  integer function sourcemod_op(modl,j)
    real,dimension(:,:),intent(in)::modl
    integer,intent(in):: j
    real:: s,x,y,z,lx,ly,lz,dx2,t0
    integer:: i
    
    sourcemod_op=1

    picksout=0.

    x=modl(1,j)
    y=modl(2,j)
    z=modl(3,j)
    s=modl(4,j)

    lx=(x-sx)
    lx=lx**2
    ly=(y-sy)
    ly=ly**2
    lz=(z-gz)
    lz=lz**2

    t0=sqrt(lx+ly+lz)*s
    trc:do i=1,ntr
       lx=(head(2,i)-x)**2
       ly=(head(1,i)-y)**2
       lz=(gz       -z)**2
       dx2=lx+ly+lz
       picksout(i)=sqrt(dx2*s*s)-t0
    end do trc
    sourcemod_op=0
  end function sourcemod_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,shz,n1,ox,dx,oy,dy,oz,dz,os,ds,nx,ny,nz,ns,n,v)
  integer:: n1,nx,ny,nz,ns,nt,nh,verb,n,v
  integer:: it,itr,ix,iy,is,iz,i3,its,update
  real::o1,d1,t,z,x,y,s,ot,oz,ox,oy,os,dt,dz,dx,dy,ds,eps,shx,shy,shz
  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)
  call from_param("nkeep",n,5); nretain=n ; write(0,*) 'nretain',nretain
  call from_param("eps",eps,0.05);call from_param("norm",norm,1)
  call from_param("ox",ox,1500.);call from_param("dx",dx,100.);call from_param("nx",nx,10)
  call from_param("oy",oy,-1500.);call from_param("dy",dy,100.);call from_param("ny",ny,10)
  call from_param("oz",oz,9500.);call from_param("dz",dz,100.);call from_param("nz",nz,10)
  call from_param("os",os,0.0002);call from_param("ds",ds,0.0000233);call from_param("ns",ns,20)
  call from_param('verb',verb,4)

  call to_history('n1',n1,"OUT");call to_history('o1',o1,"OUT");call to_history('d1',d1,"OUT")
  if(ntr.ne.n1) call seperr("number headers not matched to input")
  allocate(tmp1(4),tmp2(n))
  call sep_close()
  verb=v
  allocate(picksout(n1),head(nh,ntr))
  head=0.
  call sep_read(head,"geom")
  if(nh.eq.2) gz=68.
  sx=head(2,1)
  sy=head(1,1)
  shx=sx
  shy=sy
  shz=gz
  ox=sx+ox
  oy=sy+oy
  oz=gz+oz
!  write(0,'(a,3f10.2,f10.7,f10.2)')' ox,oy,oz,os,ov',ox,oy,oz,os,1/os
!  write(0,'(a,3f10.2,f10.7)')' dx,dy,dz,ds',dx,dy,dz,ds
!  write(0,'(a,4i4)')' nx,ny,nz,ns',nx,ny,nz,ns
!  write(0,*) '**************************************************' 
!  write(0,*) ' '
!  write(0,*)' Min/max in distance from source trace (x,y):',sx,sy
!  write(0,*) 'min model parameters  (x,y,z)'
!  write(0,*) ox-head(2,1),oy-head(1,1),oz-gz
!  write(0,*) 'min model parameters (s,v)'
!  write(0,*) os,1/os
!  write(0,*) '  '
!  write(0,*) 'max model parameters  (x,y,z)'
!  write(0,*) (nx-1)*dx+ox-head(2,1),(ny+1)*dy+oy-head(1,1),(nz-1)*dz+oz-gz
!  write(0,*) 'max model parameters (s,v)'
!  write(0,*) ((ns-1)*ds+os),1/((ns-1)*ds+os)
!  write(0,*) ' ' 
!  write(0,*) '**************************************************' 
end subroutine exh_search_init
!----------------------------------------------------------------
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

xout=mods(1)
yout=mods(2)
zout=mods(3)
sout=mods(4)

  write(0,*) '**************************************************' 
  write(0,*) ' Model space input parameter ranges '
  write(0,*)' X :',ox,' to ',maxx,' answer:',xout
  write(0,*)'   X :',ox-sx,' to ',maxx-sx,' answer:',xout-sx
  write(0,*)' Y :',oy,' to ',maxy,' answer:',yout
  write(0,*)'   Y :',oy-sy,' to ',maxy-sy,' answer:',yout-sy
  write(0,*)' Z :',oz,' to ',maxz,' answer:',zout
  write(0,*)'   Z :',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
end subroutine reportdone
!----------------------------------------------------------------
end module sourcemod
