program EventShotModel
  use sep

  implicit none
  real,dimension(:,:),allocatable:: head,dat
  real,dimension(:),allocatable::wave
  integer,dimension(:),allocatable::timearray
  real:: lx,ly,x2,s2,s,x,y,lt,t0,ll,d1,sx,sy,m
  integer:: n1,n2,n3,i,j,k,kk,nh,ntr,a,b,tt,wn1

  call sep_init()

  call to_history('esize',4)
  call from_history("n1",wn1);call from_history('d1',d1,0.004)
  call from_param('n1',n1,nint(4.5/d1));
  call from_aux("geom",'n1',nh)
  call from_aux("geom",'n2',ntr)
  call to_history ('n1',n1);call to_history('n2',ntr); call to_history('n3',ntr)

  allocate(head(nh,ntr));head=0.
  call sep_read(head,"geom")

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


  call from_param("x",x,2000.);x=x+sx
  call from_param("y",y,-40000.);y=y+sy
  call from_param("s",s,1462.);s=1./s
  call sep_close()
  allocate(dat(n1,ntr));dat=0.
  allocate(wave(wn1));wave=0.
  call sep_read(wave)

  do j=1,ntr
     sx=head(2,j)
     sy=head(1,j)
     lx=(x-sx)**2
     ly=(y-sy)**2
     x2=lx+ly    
     s2=s*s
     lt=s2*x2    
     t0=sqrt(lt)
     dat=0.
     do i=1,ntr
        lx=(head(2,i)-x)**2
        ly=(head(1,i)-y)**2
        x2=lx+ly
        ll=1.+(sqrt(s2*x2)-t0)/d1
        tt=nint(ll)
        n3=tt-wn1/2
        a=max(tt-wn1/2,1)
        b=min(wn1-1,wn1+n3-1)
        kk=1
        if( n3.lt.1) kk=1-n3
        b=min(a+b,n1)
        do k=a,b
           dat(k,i)=wave(kk)
           kk=kk+1
        end do
     end do
     call sep_write(dat)
  end do

end program EventShotModel
