program EventSum
  use sep

  implicit none
  real,dimension(:,:),allocatable:: head,dat
  integer,dimension(:),allocatable::timearray
  real,dimension(:),allocatable::power
  real:: lx,ly,x2,s2,s,x,y,lt,t0,ll,d1,sx,sy,m
  integer:: n1,n2,n3,i,j,k,gtime,nh,ntr,gloc,mm
  logical:: event
  call sep_init()
  call from_history('n1',n1);call from_history('d1',d1);call from_history('n2',n2);call from_history('n3',n3)
  call to_history ('n1',n3);call to_history('n2',1); call to_history('n3',1);call to_history('d1',1.)
  call to_history('esize',4);call from_param('event',event,.true.)

  call from_aux("geom",'n1',nh)
  call from_aux("geom",'n2',ntr)
  if(n2.ne.ntr) call seperr('n2/=ntr')
  allocate(head(nh,ntr));head=0.
  call sep_read(head,"geom")
  
  call from_param('gloc',gloc,15)
  sx=head(2,gloc)
  sy=head(1,gloc)
  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,n2));dat=0.
  allocate(power(n3));power=0.
  allocate(timearray(n2));timearray=1

  !! distances
  lx=(x-sx)**2
  ly=(y-sy)**2
  x2=lx+ly    
  s2=s*s
  lt=s2*x2    
  t0=sqrt(lt)
  
  if(event) then
     trc:do i=1,ntr
        lx=(head(2,i)-x)**2
        ly=(head(1,i)-y)**2
        x2=lx+ly
        ll=1.5+(sqrt(s2*x2)-t0)/d1
        timearray(i)=int(ll)
     end do trc
     gtime=maxval(timearray)-1
     timearray=timearray-gtime
  end if

  panel:do i=1,n3
     call sep_read(dat)
     timeshift:do k=1,n1
        lx=0.
        trace:do j=1,ntr
           if(timearray(j).ge.1.and.timearray(j).le.n1) then
              lx=lx+dat(timearray(j),j)
           end if
        end do trace
        power(i)=power(i)+abs(lx)
        timearray=timearray+1
     end do timeshift    
     timearray=timearray-n1
     if(mod(i,n3/5)==0) then
        write(0,'(2x,a,i5,a,i5,a,i3,a)') "   Finished",i,' of',n3,' panels:',nint(100.*real(i)/real(n3)),'%'     
     end if
  end do panel
  call sep_write(power)
end program EventSum
