program clean_arrival
use sep
use sorting_mod
use sep3d_struct_mod
type(sep3d) :: err_s,res_s,ang_s,baz_s
real,allocatable :: good(:,:),err(:,:),res(:,:),ang(:,:),small(:),baz(:,:)
integer :: i1,i2,nwind,res_wind,ib,ie
real :: ang_cut,err_cut,res_cut,avg,a,max_az,dir

call sep_init(SOURCE)
call init_sep3d("in",err_s,"INPUT")
call init_sep3d("res",res_s,"INPUT")
call init_sep3d("ang",ang_s,"INPUT")
call init_sep3d("back",baz_s,"INPUT")

allocate(err(err_s%n(1),err_s%n(2)))
allocate(good(err_s%n(1),err_s%n(2)))
allocate(baz(err_s%n(1),err_s%n(2)))
allocate(res(err_s%n(1),err_s%n(2)))
allocate(ang(err_s%n(1),err_s%n(2)))

call sreed("in",err,size(err)*4)
call sreed("res",res,size(err)*4)
call sreed("ang",ang,size(err)*4)
call sreed("back",baz,size(err)*4)
good=1.
call from_param("ang_cut",ang_cut,15.)
call from_param("err_cut",err_cut,.3)
call from_param("res_cut",res_cut,15.)
call from_param("dir",dir,150.)
call from_param("max_az",max_az,30.)
call from_param("res_wind",res_wind,2)
ang_cut=ang_cut/100.
res_cut=res_cut/100.
allocate(small(err_s%n(2)))
!where(abs(err)>err_cut) good=0.

do i1=1,err_s%n(1)
    i2=nint((baz(i1,1)-dir)/180.)
    if(abs(i2*180-(baz(i1,1)-dir)) >max_az) good(i1,:)=0.
    if(abs( err(i1,i2)) < .00001) good(i1,i2)=0.
    
    if(err(i1,i2)<.0001) good(i1,i2)=0.
!    if(err(i1,i2)/err(i1,i2)/=1.) write(0,*) "this test",err(i1,i2)
!    if(err(i1,i2)/err(i1,i2)/=1.) write(0,*) "this test",err(i1,i2)
  do i2=1,err_s%n(2)
    if(abs( err(i1,i2)) < .00001) good(i1,i2)=0.
    if(abs( err(i1,i2)) > err_cut) good(i1,i2)=0.
    if(abs( ang(i1,i2)) < .0001)  ang(i1,i2)=0.
  end do
  
  small=ang(i1,:)
  call quick_sort(small,err_s%n(2))
  avg=small(nint(size(small)/2.))
  do i2=1,err_s%n(2)
    if(abs((ang(i1,i2) -avg)/avg ) >ang_cut) good(i1,i2)=0.
    ib=max(i2-res_wind,1)
    ie=min(i2+res_wind,err_s%n(2))
    small(1:ie-ib+1)=res(i1,ib:ie)
    a=small(nint((ie-ib+1)/2.))
    if(abs(a-res(i1,i2))/max(a,res(i1,i2)) > res_cut .and.abs(a-res(i1,i2)) >0.5 ) then
       write(0,*) "kill3",i1,i2,res(i1,i2),a
       good(i1,i2)=0.
    end if
    if(abs(res(i1,i2)) >.3) then
       write(0,*) "kill4",i1,i2,res(i1,i2)
       good(i1,i2)=0.
    end if
  end do
end do

call srite("out",good,size(good)*4)

end program
