program reflection_pts
  use sep3d_struct_mod
  use sep
  use pts_select_mod
  implicit none
  type(sep3d) :: image_s,dipx_s,dipc_s,gamma_s,var_s,dipy_s,vel_s,pts_s
  integer     :: ndim
  integer,external :: sreed
  real, allocatable,dimension(:,:,:) :: dipc,image,var
  integer :: nenv(3),ncoh(3),nvar(3),nvel(3),n(3)
  real    :: oenv(3),denv(3),ocoh(3),dcoh(3),ovar(3),dvar(3),ovel(3),dvel(3)
  real    :: o(3),d(3)
  real, pointer :: pts(:,:)
  real, pointer :: temp1(:),temp2(:,:),temp3(:,:,:)
  integer :: i2


  call sep_init(SOURCE)

!read in parameters
  call init_sep3d("vel",vel_s,"INPUT")
  ndim=sep3d_ndims(vel_s)
  if(ndim < 2 .or. ndim >3) &
    call seperr("expecting velocity to be 2 or 3-D")
  call set_n_o_d_3d(vel_s,nvel,ovel,dvel)
  if(.not. grab_pts_select_pars().or..not.grab_parse_pts_pars(nvel,ovel,dvel))&
    call seperr("trouble grabbing parameters")


!read in files needed for pt selection

  !read in the files needed for pt sellection
  call init_sep3d("envelope",image_s,"INPUT")
  call set_n_o_d_3d(image_s,nenv,oenv,denv)
  allocate(image(nenv(1),nenv(2),nenv(3)))
  if(4*size(image)/=sreed("envelope",image,size(image)*4)) &
    call seperr("trouble reading in envelope")


  call init_sep3d("dip_coher",dipc_s,"INPUT")
  call set_n_o_d_3d(dipc_s,ncoh,ocoh,dcoh)
  allocate(dipc(ncoh(1),ncoh(2),ncoh(3)))
  if(4*size(dipc)/=sreed("dip_coher",dipc,size(dipc)*4)) &
    call seperr("trouble reading in dip coherence")

  call init_sep3d("var",var_s,"INPUT")
  call set_n_o_d_3d(var_s,nvar,ovar,dvar)
  allocate(var(nvar(1),nvar(2),nvar(3)))
  if(4*size(var)/=sreed("var",var,size(var)*4) )&
    call seperr("trouble reading in variance")

  !check to see if the files we will need later are here
  call init_sep3d("dip_xz",dipx_s,"INPUT")
  if(ndim==3) call init_sep3d("dip_yz",dipy_s,"INPUT")
  call init_sep3d("in",gamma_s,"INPUT")

!pick pts
  if(.not. select_pts(pts,vel_s,var_s,var,image_s,image,dipc_s,dipc))&
    call seperr("trouble sellecting pts")

!do i2=1,size(pts,2)
!  write(0,*) "CHECK ME",i2,pts(:,i2)
!end do

!parse pts 
  if(.not. parse_pts(pts)) call seperr("trouble parsing pts")
do i2=1,size(pts,2)
  write(0,*) "CHECK ME",i2,pts(:,i2)
end do

  deallocate(image,dipc)
!create the sep3d struture with the pts
 call init_sep3d(pts_s,"OUTPUT","FLOAT","HEADER",n=(/1,size(pts,2)/))
 call axis_allocate(pts_s,2)
 call key_allocate(pts_s,7)
 pts_s%n=(/1,size(pts,2)/)
 pts_s%keyname(1)="pos_x"
 pts_s%keyname(2)="pos_y"
 pts_s%keyname(3)="pos_z"
 pts_s%keyname(4)="var"
 pts_s%keyname(5)="dip_x"
 pts_s%keyname(6)="dip_y"
 pts_s%keyname(7)="vel"
 pts_s%keytype="scalar_float"
 pts_s%keyfmt="xdr_float"
 allocate(temp2(1,size(pts,2)))
 allocate(temp1(size(pts,2)))
 call sep3d_set_sep3d(pts_s)
 call sep3d_set_number_headers(pts_s,size(pts,2))
 call sep3d_set_key_vals(pts_s,"pos_z",pts(1,:))
 call sep3d_set_key_vals(pts_s,"pos_x",pts(2,:))
 call sep3d_set_key_vals(pts_s,"pos_y",pts(3,:))
 call sseek("var",0 ,0)
 call extract_vals("var",pts,temp1)
 call sep3d_set_key_vals(pts_s,"var",temp1)
 deallocate(var)
 call extract_vals("vel",pts,temp1)
 call sep3d_set_key_vals(pts_s,"vel",temp1)
 call extract_vals("dip_xz",pts,temp1)
 call sep3d_set_key_vals(pts_s,"dip_x",temp1)
 if(nvel(3)/=1) then
   call extract_vals("dip_yz",pts,temp1)
 else 
   temp1=0.
 end if
 call sep3d_set_key_vals(pts_s,"dip_y",temp1)
 call extract_vals("in",pts,temp1)
 call sep3d_inorder(pts_s)
 call sep3d_write_description("out",pts_s)
 temp2(1,:)=temp1 
 if(.not. sep3d_write_data("out",pts_s,temp2)) &
   call seperr("trouble writing out data")
 
  




end program
