program aa
  use sep
    implicit none
  integer :: n(5)
  integer :: i,npts,ipts,ia,ib,i1,i2
  real :: o1,o2,d1(5),d2(5)
  real :: mymax(2)
  real :: dist
  character(len=128) :: label
  integer :: nz,nang,nx,nf,iang,ngs
  
  real :: oz,oang,ox,dx,dz,dang
  real,allocatable :: pts(:,:),pos(:,:)
  real :: locx,p2a,last(2)

  call sep_init()


  call from_history("n2",npts)
  call sep_get_data_axis_par("semb",1,nz,oz,dz,label)
  call sep_get_data_axis_par("semb",4,nx,ox,dx,label)
  call sep_get_data_axis_par("semb",2,nang,oang,dang,label)
  call from_param("dist",dist,.1)
  
  call from_param("nangs",ngs)
  call to_history("n1",ngs)
 
 
  allocate(pts(4,npts),pos(2,npts))
  call from_param("locx",locx) 
  
  p2a=atan(1.)/45.
  nf=0
  iang=0
  call sreed("in",pts,npts*16)
  last=-7
  do ipts=1,npts
     if(abs(pts(2,ipts) - locx) < dist) then
       if(any(abs(last -pts(:2,ipts)) > .001)) then
          last=pts(:2,ipts)
          
          if(iang >0) then
            do i=iang,ngs
               pos(:,i)=pos(:,iang)
            end do
            call srite("out",pos,ngs*8)
            nf=nf+1
            iang=0
          end if
        end if
        iang=iang+1
        pos(2,iang)=pts(1,ipts)+(1.-pts(4,ipts))*pts(1,ipts)*tan(atan(1.)/45.*pts(3,ipts))**2
        pos(1,iang)=pts(3,ipts)
        write(0,*) "MET",iang,nf
     end if
   end do
   call to_history("esize",8)
   call to_history("n2",nf+1)
    do i=iang,ngs
               pos(:,i)=pos(:,iang)
            end do
                        call srite("out",pos,ngs*8)
 
  

end program
