# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program SourceLocate
  use sep
  implicit none

  integer:: n1,n2,n3,nx,ny,nz,ns,nt,nh,ntr,verb
  integer:: it,itr,ix,iy,is,iz,i3,its,update,norm
  real::o1,d1,t,z,x,y,s,ot,oz,ox,oy,os,dt,dz,dx,dy,ds
  real:: lz,lx,ly,lt,gz,mint,jj,t2,d2,v
  real:: tout,zout,xout,yout,vout,eps
  real,allocatable,dimension(:,:):: head
  real,allocatable,dimension(:) :: time

  call sep_init("/homes/sep/brad/Research/Valhall/Src/SourceLocate.F90")

  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_history("n2",n2)
  call from_history("n3",n3)

  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,6500.);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("ot",ot,2.500);call from_param("dt",dt,0.01);call from_param("nt",nt,50)
  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,10)

  call to_history('n1',n1,"OUT");call to_history('o1',o1,"OUT");call to_history('d1',d1,"OUT")
  if(ntr.ne.n2*n3*n1) call seperr("number headers not matched to input")

  call sep_close()

  allocate(time(n1),head(nh,ntr))
  head=0.;time=0.
  if(nh.eq.2) gz=68.

  do i3=1,n3
     mint=100000.
     call sep_read(time)
     call sep_read(head,"geom")
     ox=head(2,1)+ox
     oy=head(1,1)+oy
     oz=gz +oz
     write(0,*) 'Panel ',i3,'eps',eps
     write(0,'(a,4f10.2,f10.7,f10.2)')' ox,oy,oz,ot,os,ov',ox,oy,oz,ot,os,1/os
     write(0,'(a,3f10.2,2f10.7)')' dx,dy,dz,dt,ds',dx,dy,dz,dt,ds
     write(0,'(a,5i4)')' nx,ny,nz,nt,ns',nx,ny,nz,nt,ns
     write(0,*) '**************************************************'
     write(0,*) ' '
     write(0,*)' Min/max in distance from source trace (x,y):',head(2,1),head(1,1)
     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 (t,s,v)'
     write(0,*) ot,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 (t,s,v)'
     write(0,*) (nt-1)*dt+ot,((ns-1)*ds+os),1/((ns-1)*ds+os)
     write(0,*) ' '
     write(0,*) '**************************************************'
     its=0
     update=-1
       slo:do is=1,ns
        s=(os+ (is-1)*ds)**2
        v=1/(os+ (is-1)*ds)
        write(0,*) 'velocity',v
        tau:do it=1,nt
           t=ot+ (it-1)*dt
           t2=t*t
           dep:do iz=1,nz
              z=oz+ (iz-1)*dz
              xpo:do ix=1,nx
                 x=ox+ (ix-1)*dx
                 ypo:do iy=1,ny
                    y=oy+ (iy-1)*dy
                    lt=0.
                    its=its+1
                    trc:do itr=1,ntr
                       lx=(x-head(2,itr))**2
                       ly=(y-head(1,itr))**2
                       lz=(z-gz)**2
                       d2=lx + ly + lz
                       jj=sqrt(t2+d2*s) -t

                       jj=time(itr) - jj
                       if(norm==2) lt=lt+ jj*jj
                       if(norm==1) lt=lt+ abs(jj)
                    end do trc
                    if(norm==2) lt=sqrt(lt)
                    if(lt.lt.mint) then
                       xout=x-head(2,1)
                       yout=y-head(1,1)
                       zout=z-gz
                       tout=t
                       vout=v
                       if(verb.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,t,v',xout,yout,zout,tout,vout
                          write(0,*) ' -        -        -        -        -        -'
                       end if
                       mint=lt
                       update=update+1
                    end if
                 end do ypo
              end do xpo
           end do dep
        end do tau
     end do slo
     write(0,*)' Min/max in distance from source trace (x,y):',head(2,1),head(1,1)
     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 (t,s,v)'
     write(0,*) ot,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 (t,s,v)'
     write(0,*) (nt-1)*dt+ot,((ns-1)*ds+os),1/((ns-1)*ds+os)
     write(0,*) ' '
     write(0,*) '**************************************************'

     if(update.gt.0.or.mint.le.eps) then
        write(0,*) 'Final residual :',mint
        write(0,*) 'Updated :',update,'times'
        write(0,*)' x,y,z,t,v',xout,yout,zout,tout,vout
        write(0,*)' model space loops',its
     else
        write(0,*)' failed to update'
     end if
     do itr=1,ntr
        lx=xout-head(2,itr)+head(2,1)
        ly=yout-head(1,itr)+head(1,1)
        lz=zout-gz
        time(itr)=-tout+sqrt(lx*lx + ly*ly + lz*lz)/vout
     end do
  end do
  call sep_write(time)
end program SourceLocate
