module shootit_mod
use sep
use sep3d_struct_mod
use shoot3d_mod
use shot_type_mod
use shoot_mod
use basic_util_mod
implicit none
type(shot_par) :: par
real, private, save, pointer :: vel2d(:,:),vel3d(:,:,:)
real(kind=wp),allocatable,private  :: x_rays(:,:),z_rays(:,:),v_rays(:,:)
real,allocatable,private  ::  posx(:,:),posz(:,:),negx(:,:),negz(:,:)
complex,pointer, save,private :: match2d(:,:)
real,pointer, save,private :: rays3d(:,:),s_3d_p(:,:,:),s_3d_n(:,:,:)
contains

logical function init_shot_pars(vtag,nh,rays_s)
character(len=*) :: vtag
type(sep3d)      :: rays_s,vel_s
integer          :: ndim,nh,n(3),n2
integer,external :: sreed
real             :: maxt,o(3),d(3),zero

init_shot_pars=.false.
call init_sep3d(vtag,vel_s,"INPUT")
par%ndims=sep3d_ndims(vel_s)


call from_param("maxt",maxt,3.)
call from_param("dt",par%dt,.008)
call from_param("max_off",par%max_off,5.)
call from_param("jt",par%jt,4)
par%nt=maxt/par%dt+1

call from_param("nang",par%ntheta,20)
call from_param("oang",par%otheta,0.)
call from_param("dang",par%dtheta,2.0)

!for now we never shoot in the phi direction
par%nphi=1
par%ophi=0.
par%dphi=1.

call init_sep3d(rays_s,"INPUT","FLOAT","GRID")
! s_x,s_y,s_z,dip_xz,dip_yz,theta,phi,time
call axis_allocate(rays_s,4)
call key_allocate(rays_s,11)
rays_s%n=(/6*par%nt,par%ntheta,par%nphi,nh/)
rays_s%o=(/0.,par%otheta,par%ophi,0./)
rays_s%d=(/par%dt,par%dtheta,par%dphi,1./)
rays_s%label(1)="time"
rays_s%label(2)="theta"
rays_s%label(3)="phi"
rays_s%label(4)="pt"
rays_s%keytype="scalar_float"
rays_s%keyfmt="xdr_float"
rays_s%keyname(1)="pos_x"
rays_s%keyname(2)="pos_y"
rays_s%keyname(3)="pos_z"
rays_s%keyname(4)="dip_x"
rays_s%keyname(5)="dip_y"
rays_s%keyname(6)="gamma"
rays_s%keyname(7)="var"
rays_s%keyname(8)="theta"
rays_s%keyname(8)="theta"
rays_s%keyname(9)="phi"
rays_s%keyname(10)="time"
rays_s%keyname(11)="vel"
call sep3d_set_sep3d(rays_s)
call set_n_o_d_3d(vel_s,n,o,d)
par%min_theta=-(par%otheta+par%dtheta*(par%ntheta-1))
par%max_theta= par%otheta+par%dtheta*(par%ntheta-1)
par%nrays=1
par%nsamp=par%nt*par%jt


if(par%ndims==2) then
  par%nrays=(par%max_theta-par%min_theta)/par%dtheta+1
  zero=(par%min_theta)/par%dtheta
  if(abs(zero-nint(zero))>par%dtheta/10.) then
    write(0,*) "tracing range must include the zero ray change dtheta or otheta"
    return
  end if
  allocate(vel2d(n(1),n(2)))
  call  init_shot(n(1:2),o(1:2),d(1:2), par%nt*par%jt, &
    par%dt/par%jt,(par%ntheta-1)*2+1)
  if(4*size(vel2d)/=sreed("vel",vel2d,size(vel2d)*4)) &
    call seperr("trouble reading in velocity")
  allocate(x_rays(par%nsamp,par%nrays))
  allocate(z_rays(par%nsamp,par%nrays))
  allocate(v_rays(par%nsamp,par%nrays))
  allocate(negx(par%nt,par%ntheta))
  allocate(negz(par%nt,par%ntheta))
  allocate(posx(par%nt,par%ntheta))
  allocate(posz(par%nt,par%ntheta))
  allocate(match2d(par%nt*2, par%ntheta))
else
  allocate(vel3d(n(1),n(2),n(3)))
  if(4*size(vel3d)/=sreed("vel",vel3d,size(vel3d)*4)) &
    call seperr("trouble reading in velocity")
  call  init_shot3d(vel3d,n,o,d,par%dt)
  allocate(rays3d(3,par%nsamp))
  allocate(s_3d_n(3,par%nt, par%ntheta))
  allocate(s_3d_p(3,par%nt, par%ntheta));
end if


init_shot_pars=.true.
end function


logical function make_rays(x,y,z,dx,dy,valid,times,rays)
real    :: x,y,z,dx,dy
logical :: valid(:)
real    :: times(:)
real    :: rays(:,:,:)


if(par%ndims==2) then
   make_rays=make_rays_2d(x,y,z,dx,dy,valid,times,rays)
else
   make_rays=make_rays_3d(x,y,z,dx,dy,valid,times,rays)
end if
end function

logical function make_rays_3d(x,y,z,dx,dy,valid,times,rays)
real    :: x,y,z,dx,dy
logical :: valid(:)
real    :: times(:)
real    :: rays(:,:,:)
complex :: range
integer :: iang
real    :: pt2theta,normalx,i1,i2,normaly,ang,p2theta


make_rays_3d=.false.
p2theta=45./atan(1.)
normalx=(4.*atan(1.)-atan(dx))*p2theta
normaly=(4.*atan(1.)-atan(dy))*p2theta
normaly=0

do iang=1,par%ntheta
  ang=par%otheta+par%dtheta*(iang-1)
  if(.not. shoot3d_rays(rays3d,(/z,x,y/),(normalx-ang),&
    normaly)) write(0,*) "trouble shooting",iang,x,z,y
  s_3d_n(:,:,iang)=rays3d(:,::par%jt)
  if(.not. shoot3d_rays(rays3d,(/z,x,y/),(normalx+ang),&
    normaly)) write(0,*) "trouble shooting",iang,x,z,y
  s_3d_p(:,:,iang)=rays3d(:,::par%jt)
end do


if(.not. match_rays_3d(s_3d_n,s_3d_p,rays,valid,times)) then
   write(0,*) "trouble matching rays section"
   return
end if
make_rays_3d=.true.
end function


logical function make_rays_2d(x,y,z,dx,dy,valid,times,rays)
real    :: x,y,z,dx,dy
logical :: valid(:)
real    :: times(:)
real    :: rays(:,:,:)
complex :: range
real    :: pt2theta,normal,p2theta,theta
integer :: i1,i2,ipos,ineg,i


make_rays_2d=.false.
p2theta=45./atan(1.)
normal=(4.*atan(1.)-atan(dx))*p2theta


call shoot_rays_new(x_rays,z_rays,v_rays,vel2d,cmplx(x,z), &
  cmplx(normal+par%min_theta,normal+par%max_theta))

do i=1,par%ntheta
  theta =(i-1)*par%dtheta+par%otheta
  ipos=(theta-par%min_theta)/par%dtheta+1
  ineg=(-theta-par%min_theta)/par%dtheta+1
  negx(:,i)=x_rays(::par%jt,ineg)
  negz(:,i)=z_rays(::par%jt,ineg)
  posx(:,i)=x_rays(::par%jt,ipos)
  posz(:,i)=z_rays(::par%jt,ipos)
end do


match2d=cmplx(0.,0.)
if(.not. match_rays_2d(negx,negz,posx,posz,match2d,valid,times)) then
  write(0,*) "trouble matching ray section"
  return
end if
!call srite("bb.H",match2d,size(match2d)*8)

do i2=1,size(rays,3)
  do i1=1,size(rays,2)
    rays(1,i1,i2)=real(match2d(i1,i2))
    rays(2,i1,i2)=0.
    rays(3,i1,i2)=aimag(match2d(i1,i2))
  end do
end do
make_rays_2d=.true.
end function


logical function match_rays_3d(neg,pos,match,valid,time) 
real           :: neg(:,:,:),pos(:,:,:),match(:,:,:)
logical        :: valid(:)
real           :: time(:)
real           :: rayt(3,par%nt,2)
integer        :: iang,i
real           :: ang
                                                                                
match_rays_3d=.false.
valid=.false.
                                                                                
do iang=1,par%ntheta
  rayt(:,:,1)=neg(:,:,iang)
  rayt(:,:,2)=pos(:,:,iang)
  if(cat_ray_3d(rayt,match(:,:,iang),time(iang))) then
    if(abs(match(1,1,iang)-match(1,size(match,1),iang)) <= par%max_off) &
       valid(iang)=.true.
  end if
end do
                                                                                
time=time*par%dt/par%jt
match_rays_3d=.true.
                                                                                
end function

logical function cat_ray_3d(rays,match_rays,time_ref)
integer ::  nmatches,nmax,ilow,ihi
real :: rays(:,:,:)
real :: match_rays(:,:),ff(3)
integer i1,i2,nhi,nlow,iloc,icross
real    ::  loc,f,time_ref
logical :: crossed,bad1,bad2
                                                                                
  icross=1;iloc=1;
  time_ref=0.
  do while((icross<size(rays,2) .and. rays(3,icross,1)> 0.) .or.  icross<2)  
    icross=1+icross
  end do
  if(rays(3,icross,1) >0.) then
     bad1=.true.
  else
    bad1=.false.
     do i1=2,icross-1
       if(rays(3,i1,1)> rays(3,i1-1,1)) bad1=.true.
     end do 
                                                                                
     time_ref=icross-1+f
                                                                                
     time_ref=icross-1+f
     f=rays(3,icross-1,1)  /(rays(3,icross-1,1)-rays(3,icross,1))
     ff=rays(1:3,icross,1)*f+rays(1:3,icross-1,1)*(1.-f)
     match_rays(1:3,iloc)=ff
     do i1=icross-1,1,-1
       match_rays(1:3,iloc)=rays(1:3,i1,1)
       iloc=iloc+1
     end do
  end if
                                                                                
  icross=1;
  do while((icross<size(rays,2) .and. rays(3,icross,2)> 0. ).or. icross<2)  
    icross=icross+1
  end do
                                                                                
  if(rays(3,icross,2) >0.) then
     bad2=.true.
  else
    bad2=.false.
     do i1=2,icross-1
        if(rays(3,i1,2)> rays(3,i1-1,2)) bad2=.true.
     end do
     do i1=1,icross-1
       match_rays(1:3,iloc)=rays(1:3,i1,2)
       iloc=iloc+1
     end do
     f=rays(3,icross-1,2)  /(rays(3,icross-1,2)-rays(3,icross,2))
     match_rays(1:3,iloc)=rays(1:3,icross,2)*f+rays(1:3,icross-1,2)*(1.-f)
     match_rays(1,iloc:)=match_rays(1,iloc)
     match_rays(2,iloc:)=match_rays(2,iloc)
     match_rays(3,iloc:)=match_rays(3,iloc)
     time_ref=time_ref+icross-1+f
  end if
  if(bad1 .or. bad2)  then
    cat_ray_3d=.false.
  else 
    cat_ray_3d=.true.
  end if
                                                                                
end function
                                                                                

logical function match_rays_2d(negx,negz,posx,posz,rays,valid,time)
real  :: negx(:,:),negz(:,:)
real  :: posx(:,:),posz(:,:)
complex           :: rays(:,:)
logical        :: valid(:)
real           :: time(:)
real(kind=wp)  :: xray(par%nt,2),zray(par%nt,2)
integer        :: iang,i,j
real           :: ang
                                                                                
match_rays_2d=.false.
valid=.false.
                                                                                
do iang=1,par%ntheta
  xray(:,1)=negx(:,iang)
  zray(:,1)=negz(:,iang)
  xray(:,2)=posx(:,iang)
  zray(:,2)=posz(:,iang)
  if(cat_ray(xray,zray,rays(:,iang),time(iang))) then
    if(abs(real(rays(1,iang))-real(rays(size(Rays,1),iang))) <= par%max_off) &
       valid(iang)=.true.
  end if
end do
time=time*par%dt/par%jt
match_rays_2d=.true.
end function



logical function cat_ray(x_rays,z_rays,match_rays,time_ref)
integer ::  nmatches,nmax,ilow,ihi
real(kind=wp) :: x_rays(:,:),z_rays(:,:)
complex :: match_rays(:),ff
integer i1,i2,nhi,nlow,iloc,icross
real    ::  loc,f,time_ref
logical :: crossed,bad1,bad2
                                                                                
match_rays=0.
icross=1;iloc=1;
time_ref=0.
do while((icross<size(z_rays,1) .and. z_rays(icross,1)> 0.) .or.  icross<2) 
   icross=icross+1
end do

if(z_rays(icross,1) >0.) then
  bad1=.true.
else
  bad1=.false.

  do i1=2,icross-1
    if(z_rays(i1,1)> z_rays(i1-1,1)) bad1=.true.
  end do 

  f=z_rays(icross-1,1)    /(z_rays(icross-1,1)-z_rays(icross,1))
  ff=cmplx(x_rays(icross,1)*f+x_rays(icross-1,1)*(1.-f),&
    z_rays(icross,1)*f+z_rays(icross-1,1)*(1.-f))
  match_rays(iloc)=ff
  time_ref=icross-1+f
  do i1=icross-1,1,-1
    match_rays(iloc)=cmplx(x_rays(i1,1),z_rays(i1,1))
    iloc=iloc+1
  end do
end if


icross=1;
do while((icross<size(z_rays,1) .and. z_rays(icross,2)> 0. ).or. icross<2)
  icross=icross+1
end do
if(z_rays(icross,2) >0.) then
  bad2=.true.
else
  bad2=.false.
  do i1=2,icross-1
    if(z_rays(i1,2)> z_rays(i1-1,2)) bad2=.true.
  end do
  do i1=1,icross-1
    match_rays(iloc)=cmplx(x_rays(i1,2),z_rays(i1,2))
    iloc=iloc+1
  end do
  if(abs(z_rays(icross-1,2)-z_rays(icross,2))<.000001) then
    f=0.
  else
    f=z_rays(icross-1,2)    /(z_rays(icross-1,2)-z_rays(icross,2))
  end if
  match_rays(iloc)=cmplx(x_rays(icross,2)*f+x_rays(icross-1,2)*(1.-f), &
   z_rays(icross,2)*f+z_rays(icross-1,2)*(1.-f))

  match_rays(iloc:)=match_rays(iloc)
  iloc=iloc+1
  time_ref=time_ref+icross-1+f
end if

if(bad1 .or. bad2) then
   cat_ray=.false.
else 
   cat_ray=.true.
end if
                                                                                
end function








end module
