# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!Description:
! Module to trace rays
!
module handle_rays_3d_mod
  use shoot3d_mod
  use sep3d_struct_mod
  use param_check_mod
  use sep
  implicit none
  logical, private, save :: verb
  logical, parameter :: T=.true.,F=.false.
  real, private, save :: o_ang,d_ang,oshot,max_off,dt
  integer, private, save :: n_ang,jt,nt
  real, pointer, private, save :: vel3d(:,:,:)
  type(sep3d),private,save :: vel_s
  contains
  logical function init_3d_ray_tracer(v_s,vel,nang,oang,dang)
    real, target :: vel(:,:,:)
    type(sep3d) :: v_s
    integer :: nang,n_t,j_t
    real :: oang,dang,maxoff,d_t
    init_3d_ray_tracer=.false.
    vel_s=v_s
    vel3d=>vel
    n_ang=nang
    o_ang=oang*atan(1.)/45.
    d_ang=dang*atan(1.)/45.
    call from_param("verb",verb,.false.)
    init_3d_ray_tracer=.true.
  end function
  logical function shoot_rays_3d(tag,loc,dips,extra,matched)
    character(len=*) :: tag
    type(sep3d) :: struct
    real :: loc(:,:),dips(:,:),extra(:)
    real, allocatable :: normals(:,:),time(:)
    real, allocatable :: rays(:,:,:),match(:,:,:)
    integer :: ierr,n2,nsz,ndone,iang,iloc
    logical,pointer :: valid(:)
    real :: max_theta,min_theta,ang
    logical :: matched
    shoot_rays_3d=.false.
    write(0,*) "I DIE 1"
    if (.not. create_ray_file_3d(tag,size(dips,2),struct)) then
      write(0,*) "trouble creating ray file"
      return
    end if
    write(0,*) "I DIE a",size(dips,1),size(dips,2)
    allocate(normals(size(dips,1),size(dips,2)))
    write(0,*) "I DIE b"
    if (.not. calc_normals_3d(dips,normals)) then
      write(0,*) "trouble converting dips to surface normals"
      return
    end if
    write(0,*) "I DIE 2",o_ang,d_ang,n_ang
    if (matched) then
      max_theta=(o_ang+d_ang*(n_ang-1))
      n2=1+2*(max_theta/d_ang)
      min_theta=-max_theta
    else
      n2=n_ang
      min_theta=o_ang
      max_theta=o_ang+d_ang*(n_ang-1)
    end if
    oshot=min_theta
    oshot=min_theta
    allocate(valid(n2),rays(6,nt*jt,n2),time(n_ang), match(3,nt*jt*2&
      &,n_ang), stat=ierr)
    write(0,*) "I DIE 3"
    if (ierr.ne.0) then
      write(0,*) "trouble allocating rays"
      return
    end if
    call init_shot3d(vel3d,vel_s%n(1:3),vel_s%o(1:3),vel_s%d(1:3), dt&
      &)
    ndone=0
    do iloc=1,size(dips,2)
      if (verb .and. 1.eq.mod(iloc, nint(size(dips,2)/10.))) then
        write(0,*) "working on ",iloc," of ",size(dips,2)
      end if
      do iang=1,n2
        ang=min_theta+d_ang*(iang-1)
        if (.not. shoot3d_rays(rays(:,:,iang),loc(:,iloc),normals(1,iloc&
          &)+ang*45./atan(1.),normals(2,iloc))) then
          write(0,*) "trouble shooting ",iang,iloc
          return
        end if
! call srite("jb.H",rays(1:3:2,:,iang),size(rays(1,:,1))*8)
      end do
      if (matched) then
        if (.not. match_rays_3d(rays,match,valid,time)) then
          write(0,*) "trouble matching rays section",iloc
          return
        end if
      else
        if (.not. transfer_ray_3d(rays,match,valid,time)) then
          write(0,*) "trouble transfering ray section",iloc
          return
        end if
        write(0,*) "MATFCH",iloc
      end if
      write(0,*) "NO MATFCH",iloc
      if (.not. write_ray_block_3d(tag,struct,rays,valid,loc(:,iloc)&
        &,extra(iloc),time,ndone,iloc-1,normals(1,iloc),normals(2,iloc))) then
        write(0,*) "trouble writing ray block",iloc
        return
      end if
    end do
    call sep3d_set_ntraces(struct%tag,ndone)
    call sep3d_rite_ntraces(tag,struct%tag)
    call sep3df_delete(struct,tag)
!call sep3df_delete(struct)
    deallocate(normals)
    deallocate(match)
    deallocate(rays)
    deallocate(valid)
    shoot_rays_3d=.true.
  end function
  logical function write_ray_block_3d(tag,struct,rays,valid,loc,extra&
    &,time,ndone,iblock,dipxz,dipyz)
    character(len=*) :: tag
    type(sep3d) :: struct
    real :: rays(:,:,:)
    logical :: valid(:)
    real :: loc(:),extra,time(:),dipxz,dipyz
    integer :: ndone,i,nc,iblock,fwind(3),nwind(3)
    integer :: grid(n_ang)
    real :: sx(n_ang),sy(n_ang),sz(n_ang),ex(n_ang),ang(n_ang&
      &),dip_xz(n_ang),dip_yz(n_ang),theta(n_ang)
    write_ray_block_3d=.false.
!SET UP CONSTANTS FOR THIS BLOCK
    nc=0
    sx=loc(1)
    sy=loc(2)
    sz=loc(3)
    ex=extra
    fwind=(/0,0,iblock/)
    nwind=(/nt,n_ang,1/)
!MOVE EVERYTHING TO THE BEGINING OF THE ARRAYS
    do i=1,n_ang
      if (valid(i)) then
        nc= nc+1
        ndone= ndone+1
        grid(i)=ndone
        rays(:,:,nc)=rays(:,:,i)
        time(nc)=time(i)
        theta(nc)=o_ang+d_ang*(i-1)
        dip_xz(n_ang)=dipxz
        dip_yz(n_ang)=dipyz
      else
        grid(i)=-1
      end if
    end do
    if (nc>0) then
      call sep3d_store_grid_values(struct,grid)
      call sep3d_set_number_headers(struct,nc)
      call sep3d_set_key_vals(struct,"s_x",sx(1:nc))
      call sep3d_set_key_vals(struct,"s_y",sy(1:nc))
      call sep3d_set_key_vals(struct,"s_z",sz(1:nc))
      call sep3d_set_key_vals(struct,"dip_xz",dip_xz(1:nc))
      call sep3d_set_key_vals(struct,"dip_yz",dip_yz(1:nc))
      call sep3d_set_key_vals(struct,"theta",theta(1:nc))
      call sep3d_set_key_vals(struct,"description",ex(1:nc))
      call sep3d_set_key_vals(struct,"time",time(1:nc))
      if (.not. sep3d_write_data(tag,struct,reshape(rays(::jt,:,:nc),(&
        &/3*nt,nc/)), nwind=nwind,fwind=fwind)) then
        write(0,*) "trouble writing block",iblock
        return
      end if
    else
      call sep_put_grid_window(tag,struct%ndims ,struct%n(2:), nwind(2:&
        &),fwind(2:),(/1,1,1/),grid)
    end if
    write_ray_block_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= icross+1
    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)) then
          bad1=.true.
        end if
      end do
      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
! match_rays(:,1:iloc)=match_rays(:,iloc:1:-1)
    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)) then
          bad2=.true.
        end if
      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_3d(rays,match,valid,time)
    real :: rays(:,:,:),match(:,:,:)
    logical :: valid(:)
    real :: time(:)
    real :: rayt(3,nt*jt,2)
    integer :: iang,i
    real :: ang
    match_rays_3d=.false.
    valid=.false.
    do iang=1,n_ang
      ang=o_ang+d_ang*(iang-1)
      i=nint((-ang-oshot)/d_ang+1.)
      rayt(:,:,1)=rays(:,:,i)
      i=nint((ang-oshot)/d_ang+1.)
      rayt(:,:,2)=rays(:,:,i)
      if (cat_ray_3d(rayt,match(:,:,iang),time(iang))) then
        if (abs(match(1,1,iang)-match(1,size(match,1),iang)) <= max_off&
          &) then
          valid(iang)=.true.
        end if
      end if
    end do
    time=time*dt/jt
    match_rays_3d=.true.
  end function
  logical function transfer_ray_3d(rays,match,valid,time)
    real :: rays(:,:,:),match(:,:,:)
    logical :: valid(:)
    real :: time(:)
    integer :: iang,i
    real :: f
    transfer_ray_3d=.false.
    valid=.false.
    do iang=1,size(rays,3)
      i=1
      do while (i < size(rays,2) .and. rays(3,i,iang) > 0.)
        i= i+1
      end do
      if (rays(3,i,iang) < 0.) then
        f=(0.-rays(3,i,iang))/(rays(3,i-1,iang)-rays(3,i,iang))
        match(1:3,i,iang)=(1.-f)*rays(1:3,i-1,iang)+f*rays(1:3,i,iang)
        match(1:3,1:i,iang)=rays(:3,:i,iang)
        valid(iang)=.true.
        time(iang)=(i+f)*dt/jt
        match(1,i+1:,iang)=match(1,i,iang)
        match(2,i+1:,iang)=match(2,i,iang)
        match(3,i+1:,iang)=match(3,i,iang)
      end if
    end do
    transfer_ray_3d=.true.
  end function
  logical function create_ray_file_3d(tag,nsz,struct)
    character(len=*) :: tag
    type(sep3d) :: struct
    integer :: nsz
    create_ray_file_3d=.false.
    write(0,*) "CReATING AWI TH TAG"
    call init_sep3d(struct,"OUTPUT","FLOAT","GRID",ctag=trim(tag))
    call axis_allocate(struct,3)
    struct%n(1:3)=(/nt*3,n_ang,nsz/)
    struct%d(1:3)=(/dt,1.,1./)
    struct%label(1)="time"
    struct%unit(1)="s"
    struct%unit(2)="theta"
    struct%label(2)="ang_x"
    struct%label(3)="location"
    call key_allocate(struct,8)
    struct%keyname(1)="s_x"
    struct%keyname(2)="s_y"
    struct%keyname(3)="s_z"
    struct%keyname(4)="dip_xz"
    struct%keyname(5)="dip_yz"
    struct%keyname(6)="theta"
    struct%keyname(7)="description"
    struct%keyname(8)="time"
    write(0,*) "WEH WRONG"
    call sep3d_set_sep3d(struct)
    write(0,*) "WEH2WRONG"
    call sep3d_inorder(struct)
    write(0,*) "WEH3WRONG"
    call sep3d_write_description(tag,struct)
    write(0,*) "WEH4WRONG"
    create_ray_file_3d=.true.
  end function
  logical function calc_normals_3d(dips,normals)
    real :: dips(:,:),normals(:,:)
    real :: p2theta
    calc_normals_3d=.false.
    p2theta=45./atan(1.)
    normals=p2theta*atan(dips)
    calc_normals_3d=.true.
  end function
  logical function read_check_ray_params_3d(tag,ok)
    character(len=*) :: tag
    logical :: ok,stat_ok(20),stat_rerun(20)
    read_check_ray_params_3d=.false.
    stat_rerun=.true.
    stat_ok=.true.
    stat_rerun(1)=check_par_file("dt",dt,tag,.016)
    stat_rerun(2)=check_par_file("nt",nt,tag,nint(3.5/dt))
    stat_rerun(3)=check_par_file("jt",jt,tag,4)
    stat_rerun(4)=check_par_file("max_off",max_off,tag)
    if (all(stat_rerun)) then
      ok=.true.
    else
      ok=.false.
    end if
    if (all(stat_ok)) then
      read_check_ray_params_3d=.true.
    end if
  end function
!#############################
! #
! OLD #
! #
!#############################
end module
