program PlaneLocate
  use sep
  use planemod
  implicit none

  integer:: n1,n2,n3,nx,ny,nz,ns,nd,na,verbose,n
  integer:: ix,iy,is,iz,its,update,st
  real::z,x,y,s,v,shx,shy,shz,oxm,dx,oym,dy,ozm,dz,osm,ds,oam,da,odm,dd,mx,my,mz,ms,ma,md
  real:: zout,xout,yout,vout,sout,eps,lt,mint
  real,allocatable,dimension(:,:):: head,model,models
  real,allocatable,dimension(:)  :: time,timeout
  logical:: vel
  call sep_init(SOURCE)  
  call from_param("ox",oxm,0.);call from_param("mx",mx,10.);call from_param("nx",nx,400)
  call from_param("oy",oym,-45.);call from_param("my",my,-32.);call from_param("ny",ny,400)
  call from_param("oz",ozm,0.);call from_param("mz",mz,.1);call from_param("nz",nz,1)
  call from_param("os",osm,0.68259);call from_param("ms",ms,osm);call from_param("ns",ns,1)
  call from_param("vel",vel,.false.)
  if (vel) then
     md=ms
     ms=1./osm
     osm=1./md
  end if
  call from_param("od",odm,0.);call from_param("md",md,odm);call from_param("nd",nd,1)
  call from_param("oa",oam,0.);call from_param("md",ma,oam);call from_param("na",na,1)
  call from_param('verb',verbose,4)

  dx=(mx-oxm)/(nx-1)
  dy=(my-oym)/(ny-1)
  dz=(mz-ozm)/max((nz-1),1)
  ds=(ms-osm)/max((ns-1),1)
  da=(ma-oam)/max((na-1),1)
  dd=(md-odm)/max((nd-1),1)
write(0,*)' ds',ds
  call exh_search_init(shx,shy,n1,n)
  if(exist_file('xyplane')) then
     write(0,*)' xyplane write'
     call to_history('n1',nx,'xyplane')
     call to_history('d1',dx,'xyplane')
     call to_history('o1',oxm,'xyplane')
     call to_history('n2',ny,'xyplane')
     call to_history('d2',dy,'xyplane')
     call to_history('o2',oym,'xyplane')
     call to_history("esize",4,'xyplane')
     call to_history("sx",shx,'xyplane')
     call to_history("sy",shy,'xyplane')
  end if
  oxm=oxm*1000.+shx
  oym=oym*1000.+shy
  ozm=ozm*1000.
  osm=osm/1000.
write(0,*)' 1000. osm',osm
  dx=dx*1000.
  dy=dy*1000.
  dz=dz*1000.
  ds=ds/1000.

  allocate(time(n1),timeout(n1),model(7,1),models(size(model,1)+1,n))
  time=0.;timeout=0.;model=0.

  mint=100000.
  models(5,:)=mint
  call sep_read(time)

  its=0
  update=0
  slo:do is=1,ns
     s=(osm+ (is-1)*ds)
     v=1/(osm+ (is-1)*ds)
     write(0,*) 'velocity',v
     write(0,*) 'slowness',s
     dep:do iz=1,nz
        z=ozm+ (iz-1)*dz
        write(0,*)' depth',z
        ypo:do iy=1,ny
           y=oym+ (iy-1)*dy
!           write(0,*)' y',y
           xpo:do ix=1,nx
              x=oxm+ (ix-1)*dx
!              write(0,*)' x',x
              its=its+1
              model(1,1)=x
              model(2,1)=y
              model(3,1)=z
              model(4,1)=s
              model(5,1)=0.
              model(6,1)=0.
              st=planemod_op(model,1)
              call getpicks(timeout)
              call normmeasure(time,timeout,lt)
              if(exist_file('xyplane')) then
                 call srite('xyplane',10./lt,4)
              else
                 if(lt.lt.mint) then
                    if(verbose.gt.5)write(0,*)' Residual',lt,' (previous residual:',mint,')'
                    xout=x
                    yout=y
                    zout=z
                    vout=v
                    sout=s
                    mint=lt
                    update=update+1
                    call rankupdates(model,models,lt,update)
                    if(verbose.gt.9.and.lt.lt.eps) then
                       write(0,*)'Possible winner:'
                       write(0,*)'  error',lt,'update number:',update
                       write(0,'(a,5f10.2)')' x,y,z,v',xout,yout,zout,vout
                       write(0,*) ' -        -        -        -        -        -'
                    end if
                 end if
              end if
           end do xpo
        end do ypo
     end do dep
  end do slo

  if(update.gt.1.or.mint.le.eps) then
     write(0,*) 'Final residual :',mint
     write(0,*) 'Updated :',update,'times'
     write(0,*)' model space loops',its
  else
     write(0,*)' failed to update'
  end if
  
  model(1,1)=xout
  model(2,1)=yout
  model(3,1)=zout
  model(4,1)=sout
  st=planemod_op(model,1)
  if(.not.exist_file('xyplane')) then
     call getpicks(timeout)
     call sep_write(timeout)
     write(0,*)' * ** * ** '
     write(0,*)' Best',n, 'models'
     do ix=1,n
        write(0,*)ix,':',models(1:size(model,1),ix)
     end do
     write(0,*) '**************************************************' 
     write(0,*) ' Model space input parameter ranges '
     write(0,*)' X :',oxm,' to ',(nx-1)*dx+oxm,' answer:',xout
     write(0,*)'   X :',oxm-shx,' to ',(nx-1)*dx+oxm-shx,' answer:',xout-shx
     write(0,*)' Y :',oym,' to ',(ny-1)*dy+oym,' answer:',yout
     write(0,*)'   Y :',oym-shy,' to ',(ny-1)*dy+oym-shy,' answer:',yout-shy
     write(0,*)' Z :',ozm,' to ',(nz-1)*dz+ozm,' answer:',zout
     write(0,*)'   Z :',ozm-shz,' to ',(nz-1)*dz+ozm-shz,' answer:',zout-shz
     write(0,*)' S :',osm,' to ',(ns-1)*ds+osm,' answer:',sout
     write(0,*)'   V :',1/osm,' to ',1/((ns-1)*ds+osm),' answer:',1/sout
  else
     call aux_unlink('out')
  end if
     
end program PlaneLocate
