# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!Description:
! Module to trace rays
!
module handle_rays_2d_mod
  use shoot_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,max_off
  real, private,save, pointer :: vel2d(:,:)
  integer,private,save :: n_ang
  integer,private,save :: nt,jt
  real,private,save :: dt,oshot
  type(sep3d),private,save :: vel_s
  contains
  logical function init_2d_ray_tracer(v_s,vel,nang,oang,dang)
    real, target :: vel(:,:)
    type(sep3d) :: v_s
    integer :: nang,n_t,j_t
    real :: oang,dang,d_t
    init_2d_ray_tracer=.false.
    vel_s=v_s
    vel2d=>vel
    n_ang=nang
    o_ang=oang*atan(1.)/45.
    d_ang=dang*atan(1.)/45.
!nt=n_t
!dt=d_t
!jt=j_t
    call from_param("verb",verb,.false.)
    init_2d_ray_tracer=.true.
  end function
  logical function shoot_rays_2d(tag,loc,dips,extra,match)
    character(len=*) :: tag
    type(sep3d) :: struct
    real :: loc(:,:),dips(:,:),extra(:)
    real, allocatable :: normals(:,:),time(:)
    complex,allocatable :: rays(:,:)
    real(kind=wp),allocatable :: x_rays(:,:),z_rays(:,:),v_rays(:,:)
    integer :: ierr,n2,nsz,ndone,iloc
    logical,pointer :: valid(:)
    real :: max_theta,min_theta
    logical :: match
    shoot_rays_2d=.false.
    if (.not. create_ray_file_2d(tag,size(dips,2),struct,match)) then
      write(0,*) "trouble creating ray file"
      return
    end if
    allocate(normals(size(dips,1),size(dips,2)))
    dips=0.
    if (.not. calc_normals_2d(dips,normals)) then
      write(0,*) "trouble converting dips to surface normals"
      return
    end if
    if (match) 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
    allocate(x_rays(nt*jt,n2),z_rays(nt*jt,n2),v_rays(nt*jt,n2),valid&
      &(n2),rays(nt*jt*2,n_ang),time(n_ang),stat=ierr)
    if (ierr.ne.0) then
      write(0,*) "trouble allocating rays"
      return
    end if
    call init_shot(vel_s%n(1:2),vel_s%o(1:2),vel_s%d(1:2), nt*jt,dt/jt&
      &,n2)
    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
      call shoot_rays_new(x_rays,z_rays,v_rays,vel2d,cmplx(loc(1,iloc)&
        &,loc(3,iloc)),cmplx(normals(1,iloc)+min_theta*45./atan(1.)&
        &,normals(1,iloc)+45./atan(1.)*max_theta))
      if (match) then
        if (.not. match_rays_2d(x_rays,z_rays,rays,valid,time)) then
          write(0,*) "trouble matching rays section",iloc
          return
        end if
      else
        if (.not. transfer_ray(x_rays,z_rays,rays,valid,time)) then
          write(0,*) "trouble transfering ray section",iloc
          return
        end if
      end if
      if (.not. write_ray_block(tag,struct,rays,valid,loc(:,iloc),extra&
        &(iloc),time,ndone,iloc-1,normals(1,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)
    deallocate(normals,x_rays,z_rays,v_rays,rays,valid)
    shoot_rays_2d=.true.
  end function
  logical function write_ray_block(tag,struct,rays,valid,loc,extra,time&
    &,ndone,iblock,dipxz)
    character(len=*) :: tag
    type(sep3d) :: struct
    complex :: rays(:,:)
    logical :: valid(:)
    real :: loc(:),extra,time(:)
    integer :: ndone,i,nc,iblock,fwind(3),nwind(3)
    integer :: grid(n_ang),coord(2,n_ang)
    real :: sx(n_ang),sy(n_ang),sz(n_ang),ex(n_ang),ang(n_ang&
      &),dip_xz(n_ang),dipxz,dip_yz(n_ang),theta(n_ang)
    write_ray_block=.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/)
    dip_xz=dipxz
    dip_yz=0.
    call sep3d_set_window(struct,nwind=nwind,fwind=fwind)
!MOVE EVERYTHING TO THE BEGINING OF THE ARRAYS
    grid=-1
    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)
        coord(:,nc)=(/i,fwind(3)+1/)
      else
        grid(i)=-1
      end if
    end do
    if (nc>0) then
      call sep3d_store_grid_values(struct,grid)
! if(.not. sep3d_set_coords(struct,coord))
! call seperr("trouble setting coordinates \n")
      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))
! call sep3d_inorder(struct)
      if (.not. sep3d_write_data(tag,struct,rays(::jt,: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=.true.
  end function
  logical function transfer_ray(x_rays,z_rays,rays,valid,time)
    real(kind=wp) :: x_rays(:,:),z_rays(:,:)
    complex :: rays(:,:)
    logical :: valid(:)
    real :: time(:)
    integer :: iang,i
    real :: f
    transfer_ray=.false.
    valid=.false.
    do iang=1,size(x_rays,2)
      i=1
      do while (i < size(x_rays,1) .and. z_rays(i,iang) > 0.)
        i= i+1
      end do
      if (z_rays(i,iang) < 0.) then
        f=(0.-z_rays(i,iang))/(z_rays(i-1,iang)-z_rays(i,iang))
        z_rays(i,iang)=(1.-f)*z_rays(i-1,iang)+f*z_rays(i,iang)
        x_rays(i,iang)=(1.-f)*x_rays(i-1,iang)+f*x_rays(i,iang)
        rays(1:i,iang)=cmplx(real(x_rays(1:i,iang)),real(z_rays(1:i,iang&
          &)))
        valid(iang)=.true.
        time(iang)=(i+f)*dt/jt
        rays(i+1:,iang)=rays(i,iang)
      end if
    end do
    transfer_ray=.true.
  end function
  logical function match_rays_2d(x_rays,z_rays,rays,valid,time)
    real(kind=wp) :: x_rays(:,:),z_rays(:,:)
    complex :: rays(:,:)
    logical :: valid(:)
    real :: time(:)
    real(kind=wp) :: xray(nt*jt,2),zray(nt*jt,2)
    integer :: iang,i,j
    real :: ang
    match_rays_2d=.false.
    valid=.false.
    do iang=1,n_ang
      ang=o_ang+d_ang*(iang-1)
      i=nint((-ang-oshot)/d_ang+1.)
      xray(:,1)=x_rays(:,i)
      zray(:,1)=z_rays(:,i)
      j=nint((ang-oshot)/d_ang+1.)
      xray(:,2)=x_rays(:,j)
      zray(:,2)=z_rays(:,j)
      if (cat_ray(xray,zray,rays(:,iang),time(iang))) then
        if (abs(real(rays(1,iang))-real(rays(size(Rays,1),iang))) <=&
          & max_off) then
          valid(iang)=.true.
        end if
      end if
    end do
    time=time*dt/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)) then
          bad1=.true.
        end if
      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)) then
          bad2=.true.
        end if
      end do
      do i1=1,icross-1
        match_rays(iloc)=cmplx(x_rays(i1,2),z_rays(i1,2))
        iloc= iloc+1
      end do
      f=z_rays(icross-1,2) /(z_rays(icross-1,2)-z_rays(icross,2))
      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
  logical function calc_normals_2d(dips,normals)
    real :: dips(:,:),normals(:,:)
    real :: p2theta
    calc_normals_2d=.false.
    p2theta=45./atan(1.)
    normals=180.-p2theta*atan(dips)
    calc_normals_2d=.true.
  end function
  logical function create_ray_file_2d(tag,nsz,struct,matched)
    character(len=*) :: tag
    type(sep3d) :: struct
    integer :: nsz
    logical :: matched
    create_ray_file_2d=.false.
    call init_sep3d(struct,"OUTPUT","COMPLEX","GRID",ctag=tag)
    call axis_allocate(struct,3)
    if (matched) then
      struct%n(1:3)=(/nt*2,n_ang,nsz/)
    else
      struct%n(1:3)=(/nt,n_ang,nsz/)
    end if
    struct%d(1:3)=(/dt,d_ang,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"
    struct%keytype="scalar_float"
    struct%keyfmt="scalar_float"
    call sep3d_set_sep3d(struct)
    call sep3d_inorder(struct)
    call sep3d_write_description(tag,struct)
    create_ray_file_2d=.true.
  end function
  logical function read_check_ray_params_2d(tag,ok)
    character(len=*) :: tag
    logical :: stat_ok(20),stat_rerun(20),ok
    read_check_ray_params_2d=.false.
    stat_ok=.true.
    stat_rerun=.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_2d=.true.
    end if
  end function
end module
