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, pointer :: dip(:),good(:),err(:,:),baz(:),var(:),x(:),y(:),x2(:),y2(:)


  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("var",var_s,"INPUT")
  call init_sep3d("x",x_s,"INPUT")
  call init_sep3d("y",y_s,"INPUT")

  allocate(err(1,err_s%n(1)*err_s%n(2)))
  allocate(dip(dip_s%n(1)*dip_s%n(2)))
  allocate(good(good_s%n(1)*good_s%n(2)))
  allocate(baz(baz_s%n(1)*baz_s%n(2)))
  allocate(var(var_s%n(1)*var_s%n(2)))
  allocate(x2(x_s%n(1)*baz_s%n(1)))
  allocate(y2(y_s%n(1)*baz_s%n(1)))
  allocate(x(x_s%n(1)))
  allocate(y(y_s%n(1)))
  allocate(coords(2,baz_s%n(1)*baz_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(var)/=sreed("var",var,size(var)*4)) &
    call seperr("trouble reading in variance")

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

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

  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(i2)
      y2(i)=y(i2)
    end do
  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)
  do i=2,size(x2)
   x2(i)=sqrt((x2(i)-x2(1))**2+(y2(i)-y2(1))**2)
  end do
  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_inorder(pts_s)
  call sep3d_set_key_vals(pts_s,"pos_x",x2)
  call sep3d_set_key_vals(pts_s,"dip",dip)
  call sep3d_set_key_vals(pts_s,"baz",baz)
  call sep3d_set_key_vals(pts_s,"good",good)
  call sep3d_set_key_vals(pts_s,"var",var)

  call sep3d_set_inorder(pts_s)
  if(.not. sep3d_update_ntraces(pts_s)) &
     call seperr("trouble updating ntraces")
  call sep3d_inorder(pts_s)
  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 ntraces")
  call sep3d_rite_num_traces("out",pts_s)

end program
