program reflection_pts
  use sep3d_struct_mod
  use sep
  use pts_select_mod
  implicit none
  integer     :: i,i1,i2
  integer,external :: sreed
  integer, pointer :: coords(:,:)
  type(sep3d) :: err_s,dip_s,good_s,baz_s,var_s,x_s,y_s,pts_s
  real :: shift
  real, pointer :: dip(:),good(:,:),err(:,:),baz(:),x(:,:),x2(:)
  real, pointer :: dip2(:),baz2(:),good2(:)
  integer ::itot


  call sep_init(SOURCE)
  call init_sep3d("in",err_s,"INPUT")
  call init_sep3d("dip",dip_s,"INPUT")
  call init_sep3d("good",good_s,"INPUT")
  call init_sep3d("baz",baz_s,"INPUT")
  call init_sep3d("x",x_s,"INPUT")


  allocate(good2(good_s%n(1)*good_s%n(2)))
  allocate(good(good_s%n(1),good_s%n(2)))
  allocate(dip(dip_s%n(1)))
  allocate(baz(baz_s%n(1)))
  allocate(err(1,err_s%n(1)*err_s%n(2)))
  allocate(x(x_s%n(1),x_s%n(2)))

  allocate(dip2(x_s%n(1)*err_s%n(2)))
  allocate(baz2(x_s%n(1)*err_s%n(2)))
  allocate(x2(x_s%n(1)*err_s%n(2)))
  allocate(coords(2,err_s%n(1)*err_s%n(2)))



  if(4*size(err)/=sreed("in",err,size(err)*4)) &
    call seperr("trouble reading in errors")

  if(4*size(dip)/=sreed("dip",dip,size(dip)*4)) &
    call seperr("trouble reading in dip")

  if(4*size(good)/=sreed("good",good,size(good)*4)) &
    call seperr("trouble reading in good")

  if(4*size(baz)/=sreed("baz",baz,size(baz)*4)) &
    call seperr("trouble reading in baz")


  if(4*size(x)/=sreed("x",x,size(x)*4)) &
    call seperr("trouble reading in x")

  call from_param("shift",shift,100.)


  call  to_history("n1",err_s%n(2),"lens")

  do i1=1,err_s%n(1)
    do i2=1,err_s%n(2)
      if (good(i1,i2)>0.001) then
        i=i+1
        x(i1,:)=x(i1,i2)
      end if 
    end do
  end do

  itot=0


  i=0
  do i2=1,err_s%n(2)
    do i1=1,err_s%n(1)
      i=i+1
      coords(1,i)=i1
      coords(2,i)=i2
      x2(i)=x(i1,i2)
      dip2(i)=dip(i2)
      if(baz(i2) >30. .and. baz(i2)  < 65) dip2(i)=-dip2(i)
      baz2(i)=baz(i2)
      if(good(i1,i2) > .001) then
        good2(i)=good(i1,i2)
        if(abs(cos(baz2(i)*atan(1.)/45.)) < .3) then
          good(i1,i2)=0.
          good2(i)=0.
!          write(0,*) "BAD",i,baz2(i)
        else 
!          write(0,*) "GOO",i,baz2(i)
        end if
      else 
!        write(0,*) "BD2",i,baz(i2)
      end if
      if(baz(i2) > 210. .and. baz(i2)  < 245) then
!        write(0,*) "CHANGING DIP",dip2
         dip2(i)=-dip2(i)
         good2(i)=0.
      end if
    end do
  end do
  x2=x2-shift

  
  do i2=1,err_s%n(2)
    i=0
    do i1=1,err_s%n(1)
      if (good(i1,i2)>0.001) then
        i=i+1
        x(i1,:)=x(i1,i2)
      end if 
    end do
    call srite("lens",i,4)
  end do


  call init_sep3d(pts_s,"OUTPUT","FLOAT","GRID",n=(/1,good_s%n(1),good_s%n(2)/))
  call key_allocate(pts_s,5)
  pts_s%keyname(1)="pos_x"
  pts_s%keyname(2)="dip"
  pts_s%keyname(3)="baz"
  pts_s%keyname(4)="good"
  pts_s%keyname(5)="var"
  pts_s%keytype="scalar_float"
  pts_s%keyfmt="xdr_float"
  call sep3d_set_sep3d(pts_s)
  call sep3d_set_number_headers(pts_s,size(err))
  if(.not. sep3d_set_coords(pts_s,coords))&
   call seperr("trouble setting coords")
  x2(1)=0
  call sep3d_set_key_vals(pts_s,"pos_x",x2)
  call sep3d_inorder(pts_s)
  dip2=dip2*atan(1.)/45
  call sep3d_set_key_vals(pts_s,"dip",dip2)
  call sep3d_set_key_vals(pts_s,"baz",baz2)
  call sep3d_set_key_vals(pts_s,"good",good2)
  good2=1. 
  call sep3d_set_key_vals(pts_s,"var",good2)

  call sep3d_write_description("out",pts_s)
  if(.not. sep3d_write_data("out",pts_s,err)) &
   call seperr("trouble writing out data")
  if(.not. sep3d_update_ntraces(pts_s))&
   call seperr("trouble updating number of traces")
  call sep3d_rite_num_traces("out",pts_s)

end program
