!
!                  This file is copyrighted by Colorado
!                  State University. It may be distributed
!                  freely, as long as it is distributed 
!                  completely, including this copyright
!                  message.
!
module kinds
!
  implicit none
  integer, parameter :: kindi = selected_int_kind(9)
  integer, parameter :: kindr = selected_real_kind(p=11)
end module kinds
! -------------------------------------------------------------------
module random
  use kinds
!
  integer(kind=kindi), parameter :: modulus = 2147483647
  integer(kind=kindi) :: seed(17), it1 = 17, it2 = 12
!
  real(kind=kindr), parameter :: odenom = 1.d0/(2.d0**31)
end module random
! -------------------------------------------------------------------
module constants
  use kinds
!
  real(kind=kindr) ::  pi, p2i, pi2, pi4, pi4inv, dtor, rtod
end module constants
! -------------------------------------------------------------------
module cylinders
  use kinds
!
  real(kind=kindr) :: zmax_cyls, zmin_cyls, zbot_cyls
  real(kind=kindr) :: xmin, ymin, zmin, xmax, ymax, zmax
  integer(kind=kindi), allocatable, dimension (:) :: matno_glass, matno_cyli 
  real(kind=kindr), allocatable, dimension(:) :: y0, z0, ro, ri, &
      length_cyl, beta_cyl, betad
  real(kind=kindr), allocatable, dimension(:) :: xmin_cyl, xmax_cyl, &
       ymin_cyl, ymax_cyl, zmin_cyl, zmax_cyl
end module cylinders
! -------------------------------------------------------------------
module surfaces
  use kinds
!
  integer(kind=kindi), allocatable, dimension (:) :: matno_surf
  real(kind=kindr), allocatable, dimension (:,:) :: x, y, z
  real(kind=kindr), allocatable, dimension (:) :: cax, cbx, cgx, &
      cay, cby, cgy, caz, cbz, cgz
end module surfaces
! -------------------------------------------------------------------
module directions
  use kinds
!
  real(kind=kindr), allocatable, dimension(:) :: phid, altituded, phi, & 
    altitude, thetatd, thetald, ta, tp, cp, ca, sp, sa, exa, eya, eza  
  real(kind=kindr), allocatable, dimension(:) :: thetal, thetat, ttt, ttl
  real(kind=kindr) :: ex, ey, ez
  real(kind=kindr) :: ex_plane, ey_plane, ez_plane
end module directions
! -------------------------------------------------------------------
module materials
  use kinds
!
  integer(kind=kindi) :: itheta_i, nc, nc_int, nc_save, matnob
  real(kind=kindr) :: c_theta_i, theta_i, theta_o, phi_i, phi_o, &
      theta_o_90, c_theta_o, s_theta_o, c_phi_o, s_phi_o, e_r, e_theta, &
      delta_theta, eysave, a11, a12, beta, sb, cb, l_glass
!
  real(kind=kindr), allocatable, dimension(:) :: rho_s, rho_ss, rho_d, &
      r_ss, r_d, rp1inv_ss, rp1inv_d
  real(kind=kindr), allocatable, dimension(:,:) ::  pdf_surf, rho_glass, &
      tau_glass, drho_glass, dtau_glass
  real(kind=kindr), allocatable, dimension(:)  :: n_glass, k_glass, t
end module materials
! -------------------------------------------------------------------
module control
  use kinds
!
  character(len=1), dimension(72) :: head
  integer(kind=kindi) :: num_photons, nplost_max, num_ploops
  integer(kind=kindi) :: seed0, seed_from_time, nd, np, nploop
  integer(kind=kindi) :: num_dirs, n_long, n_tran, num_cyls, num_surfs, &
     num_mats_surf, num_mats_glass, endcap, norm
  integer(kind=kindi) :: idebug, idata, iprint, iprint_ran, iglass 
  real(kind=kindr) :: tol
!
end module control
! -------------------------------------------------------------------
module photons
  use kinds
!
  integer(kind=kindi) :: npabs_run, nplost_run, int_no, iabs
  integer(kind=kindi) :: int0_cyl, int0_glass, int0_surf, &
     int1_cyl, int1_glass, int1_surf
  integer(kind=kindi), allocatable, dimension (:,:) :: npabs
  integer(kind=kindi), allocatable, dimension(:) :: npabs_cyls, npabs_tot, &
     npabs_glass, npabs_surf, nplost, npabs_sum, matno_int, npabs_a, &
     np_tot, npabs_caps
  real(kind=kindr), allocatable, dimension (:) :: iam_cyls, iam_glass, & 
     tau_alpha
  real(kind=kindr) :: dzint_old, dzintb, tad, glass_frac 
  real(kind=kindr) :: xi(2), yi(2), zi(2), dzint(2)
  real(kind=kindr) :: xe, ye, ze, a, b, b1, c, c1, ez2, ey2, rad
  real(kind=kindr) :: yehat, zehat, nx, ny, nz, xint, yint, zint
end module photons
! -------------------------------------------------------------------
module planes
  use kinds
!
  real(kind=kindr) :: xminb, xmaxb, yminb, ymaxb, xmina, xmaxa, &
    ymina, ymaxa, zmaxe, ztop, zbot = 0.d0
  real(kind=kindr) :: xep(3,4), areaa, areab, areae
end module planes
! -------------------------------------------------------------------
module dates
  use kinds
!
  character(len=8) :: adate, atime
  integer :: valuesa(8,0:2)
end module dates
!
!*********************************************************************
!                                                                    *
program iam                                                          
!                                                                    *
!*********************************************************************
!*     program to compute IAM's for simple cylindrical geometries    *
!*********************************************************************
!
!     - modules
!
  use kinds
  use random
  use constants
  use cylinders
  use surfaces
  use directions
  use materials
  use control
  use photons
  use planes
  use dates
!
  implicit none
!
!  Variable definitions
!
  character(len=1), dimension(40) :: headin
  character(len=8) :: version
  character(len=8) :: adate_compiled
  integer(kind=kindi) :: i, iu, nlines7, npdo, m, n
  integer(kind=kindi) :: nd_converged, nd_not_converged
  real(kind=kindr) :: tol_check, zz = 1.96
!
!  get initial date and time
!
  call getdt (0)
  call getdt (1) ! Just in case of a fatal error during input
!
!  Version Information
!
  adate_compiled(1:8) = '1-11-98'
  version(1:8)        = 'ver. 2.2'
!
!  initialize constants
!
  pi = 4.d0*atan(1.d0)
  p2i = 2.d0*pi
  pi2 = 0.5d0*pi
  pi4 = 0.25d0*pi
  pi4inv = 1.d0/pi4
  dtor = pi/180.d0
  rtod = 1.d0/dtor
  if (idebug /= 0) then
    write (9, 7777) pi, p2i, pi2, pi4, dtor, rtod
  end if
!
!     open files
!
!      - file 1: input file w/o comments (& in col. 1), containing:
!            + control variables
!            + directions
!            + geometry
!            + material properties
!      - file 2: output file of iam's in tabular format
!      - file 3: output file of geometry
!      - file 4: output file of lost photons
!      - file 7: input file containing:
!            + control variables
!            + directions
!            + geommetry
!            + material properties
!      - file 8: debug info. file 
!      - file 9: output file 
!
  iu = 1
  open (iu, file = 'iam.tmp',  status = 'unknown', err = 4000)
  iu = 2
  open (iu, file = 'iam.tab',  status = 'unknown', err = 4000)
  iu = 3
  open (iu, file = 'iam.plt',  status = 'unknown', err = 4000)
  iu = 4
  open (iu, file = 'iam.lst',  status = 'unknown', err = 4000)
  iu = 7
  open (iu, file = 'iam.in',   status = 'old',     err = 4000)
  iu = 9
  open (iu, file = 'iam.out',  status = 'unknown', err = 4000)
  iu = 10
  open (iu, file = 'iam.chk',  status = 'unknown', err = 4000)
!
!
  write (*,*) ' '
  write (*,*) '                 Copyrighted by Colorado State University'
  write (*,*) ' '
  write (9,*) '                 Copyrighted by Colorado State University'
  write (9,*) ' '
!
!      convert file  (strip comments w/ & in col. 1)
!
  iu = 7
  call files (iu, 1, nlines7)
  close (iu)
  write (*,*) '  ', nlines7,' - lines read from input file '
  write (*,2004) 
!
!      read control
!
  iu = 1
  read (iu, 1000) (headin(i), i=1,40)
  read (iu, *) num_photons, nplost_max, num_ploops, seed0, tol
  read (iu, *) n_long, n_tran, num_cyls, num_surfs, num_mats_surf, &
               num_mats_glass
  read (iu,*) iprint, iprint_ran, idebug, idata, iglass
  num_dirs = n_long*n_tran
!
!  aperture plane
!
  read (iu, *)  xmina, xmaxa, ymina, ymaxa
!
! back plane
!
  read (iu, *)  xminb, xmaxb, yminb, ymaxb, matnob
!
!  Check Input
!
  call check_control
!
!  Select initial seed, if 0, get from time
!
  if (seed0 == 0) then
    seed0 = seed_from_time
  elseif (seed0 < 0) then
    write (*,*) 'Error -- seed0 must be >= 0 -- stopping'
    write (9,*) 'Error -- seed0 must be >= 0 -- stopping'
    call adios (2)
  end if
  if (mod(seed0, 2) /= 1) then
    write (*,*) 'Error -- seed0 must be odd - adding one and continuing'
    write (9,*) 'Error -- seed0 must be odd - adding one and continuing'
    seed0 = seed0 + 1
  end if
!
!  Debug Information
!
  if (idebug /= 0) then
    iu = 8
    open (iu, file = 'iam.dbg',  status = 'unknown', err = 4000)
  end if
!
!  write title page
!
  head = ' '
  head(1:40) = headin(1:40)
  do i = 1,8
    head(54+i) = adate(i:i)
    head(64+i) = atime(i:i)
  end do
  call header (9)
  write (9, 2010) head(1:40), adate, atime, version, adate_compiled
  write (9, 2011) num_photons, num_ploops, tol, seed0
  write (9, 2012) num_dirs, num_cyls, num_surfs, num_mats_surf, num_mats_glass
  write (9, 2013) xmina, xmaxa, ymina, ymaxa, &
     xminb, xmaxb, yminb, ymaxb, matnob
  write (9, 2014) iprint, iprint_ran, idebug, idata, iglass
!
!  Calculation of indices
!
!   - npabs(0, nd): top plane
!   - npabs(1, nd): bottom plane (outside back plane)
!   - npabs(2, nd): back plane (outside aperture plane)
!   - npabs(int0_cyl:int1_cyl, nd): absorption by inside cylinder
!   - npabs(int0_glass:int1_glass,nd): absorption by glass
!   - npabs(int0_surf:int1_surf,nd): absorption by surfaces
!
!     - indices
!
  int0_cyl = 3
  int1_cyl = int0_cyl + num_cyls - 1
  int0_glass = int1_cyl + 1   
  int1_glass = int0_glass + num_cyls - 1   
  int0_surf = int1_glass + 1   
  int1_surf = int0_surf + num_surfs - 1
  int1_surf = max(int1_surf, int0_surf)
  write (9, 2015) int0_cyl, int1_cyl, int0_glass, int1_glass, &
     int0_surf, int1_surf
!
!  Initialize RNG
!
  call raninit (iprint_ran, seed0)
!
!  Allocate Arrays
!
!   - directions
!
  write (*,*) 'Allocating direction arrays of length - ', num_dirs
  if (num_dirs > 0) then
    allocate (thetal(n_long), thetat(n_tran))
    allocate (thetald(n_long+1), thetatd(n_tran+1))
    allocate (phi(num_dirs),  altitude(num_dirs))
    allocate (phid(num_dirs), altituded(num_dirs))
    allocate (cp(num_dirs), sp(num_dirs), tp(num_dirs), ttt(n_tran))
    allocate (ca(num_dirs), sa(num_dirs), ta(num_dirs), ttl(n_long))
    allocate (exa(num_dirs), eya(num_dirs), eza(num_dirs))
  end if
!
  call read_directions
!
!  - cylinders
!
  write (*,*) 'Allocating cylinder arrays of length - ', num_cyls
  allocate (xmax_cyl(num_cyls), y0(num_cyls), z0(num_cyls))
  allocate (ro(num_cyls), ri(num_cyls), length_cyl(num_cyls))
  allocate (beta_cyl(num_cyls), betad(num_cyls))
  allocate (matno_glass(num_cyls), matno_cyli(num_cyls))
  allocate (xmin_cyl(num_cyls), ymax_cyl(num_cyls), ymin_cyl(num_cyls))
  allocate (zmin_cyl(num_cyls), zmax_cyl(num_cyls))
  call read_cylinders
  ztop = zmax_cyls *1.01
!
!  Absorption Arrays
!
  write (*,*) 'Allocating absorption arrays'
!
!     - 1-D arrays
!
  allocate (npabs_cyls(num_dirs), npabs_glass(num_dirs)) 
  allocate (npabs_surf(num_dirs), npabs_caps(num_dirs)) 
  allocate (iam_cyls(num_dirs), iam_glass(num_dirs), tau_alpha(num_dirs)) 
  allocate (nplost(num_dirs), npabs_sum(num_dirs), npabs_tot(num_dirs))
  allocate (npabs_a(num_dirs), np_tot(num_dirs))
!
!     - 2-D arrays
!
  allocate (npabs(0:int1_surf, num_dirs))
!
!   - surfaces
!
  write (*,*) 'Allocating surface arrays of length - ', num_surfs
  if (num_surfs > 0) then
    allocate (matno_surf(num_surfs))
    allocate (x(4,num_surfs), y(4,num_surfs), z(4,num_surfs))
    allocate (cax(num_surfs), cay(num_surfs), caz(num_surfs))
    allocate (cbx(num_surfs), cby(num_surfs), cbz(num_surfs))
    allocate (cgx(num_surfs), cgy(num_surfs), cgz(num_surfs))
    call read_surfaces
  end if
!
  areaa = (xmaxa-xmina)*(ymaxa-ymina)
!
!   - materials
!
  write (*,*) 'Allocating surface material arrays of length - ', num_mats_surf
  allocate (matno_int(0:int1_surf))
  allocate (r_d(num_mats_surf), r_ss(num_mats_surf))
  allocate (rp1inv_d(num_mats_surf), rp1inv_ss(num_mats_surf))
  allocate (rho_s(num_mats_surf), rho_ss(num_mats_surf), rho_d(num_mats_surf))
  allocate (pdf_surf(3,num_mats_surf))
  if (num_mats_glass > 0) then
    write (*,*) 'Allocating glass material arrays of length - ', num_mats_glass
    allocate (n_glass(num_mats_glass),k_glass(num_mats_glass),t(num_mats_glass))
    allocate (rho_glass(0:91,num_mats_glass), drho_glass(0:91,num_mats_glass))
    allocate (tau_glass(0:91,num_mats_glass), dtau_glass(0:91,num_mats_glass))
  end if
  call read_materials
!
!  Set Intersection Material Numbers   
!
  matno_int = 0
!  
!   - top and bottom planes
!
  matno_int(0:1) = 0
!  
!   - back plane
!
  matno_int(2) = matnob
!
!    - cylinders: surfaces
!
  do nc = 1, num_cyls
    matno_int(int0_cyl+nc-1) = matno_cyli(nc)
  end do
!
!    - cylinders: glass
!
  do nc = 1, num_cyls
    matno_int(int0_glass+nc-1) = matno_glass(nc)
  end do
!
!    - surfaces
!
  do nc = 1, num_surfs
    matno_int(int0_surf+nc-1) = matno_surf(nc)
  end do
  if (idebug /= 0) then
    call header (9)
    write (9,*) 'Intersection Array of Material Numbers from 0:',int1_surf
    write (9,*)
    write (9,*)
    write (9,*) matno_int
    write (9,*)
    write (9,*)
  end if
!
!  finished input phase
!
  write (*, *) '         Finished Data Check '
  write (9,*)
  write (9, *) '         Finished Data Check '
  call getdt (1)
!
!  data check termination
!
  if (idata /= 0) then
    call adios(1)
  end if
!
! Photon initializations for run
!
!     - constants
!
  npabs_run = 0
  nd_converged = 0
  nd_not_converged = 0
!
!     - 1-D arrays
!
  nplost = 0
  npabs_glass = 0
  npabs_a = 0
  np_tot = 0
!
!     - 2-D arrays
!
  npabs = 0
!
! Write Convergence Header
!
  write (*,2003)
!
!     photon emission loop
!
  do nd = 1, num_dirs
!
!     Initialize array counting photons absorbed by cylinder end caps:
!
    npabs_caps(nd) = 0
!
!     compute emission direction, ex, ey and ez
!        - here +tive from origin to plane
!        - formulated for quadrant 1, not tested otherwise
!
    ex = exa(nd)
    ey = eya(nd)
    ez = eza(nd)
    call geteplane
!
!     reverse e: now negative (pointing from e-plane back down toward
!             origin)
!
    ex_plane = -ex
    ey_plane = -ey
    ez_plane = -ez
!
!     photon convergence loop - blocked
!
    npdo = 0
    do nploop = 1, num_ploops
      npdo = npdo + 1
!
!    intersection
!
      if (iglass == 0) then
        call intsec_ng
      else
        call intsec
      end if
!
!  photon and iam computations
!     - now, iam is based on aperture plane
!
      npabs_cyls(nd) = sum(npabs(int0_cyl:int1_cyl,nd))
      if (iglass == 0) then
        npabs_glass(nd) = 0
      else
        npabs_glass(nd) = sum(npabs(int0_glass:int1_glass,nd))
      end if
      if (num_surfs == 0) then
        npabs_surf(nd) = 0
      else
        npabs_surf(nd) = sum(npabs(int0_surf:int1_surf,nd))
      end if
      npabs_sum(nd) = npabs_cyls(nd) + npabs_surf(nd) + npabs_glass(nd) &
                    + npabs_caps(nd)
      npabs_tot(nd) = sum(npabs(:,nd)) + npabs_caps(nd)
      np_tot(nd) = npabs_tot(nd) + nplost(nd)
!
! Printing
!
      if (idebug /= 0) then
        iu = 8
        call print_npabs (iu)
      end if
      if (iprint < 0) then
        iu = 9
        call print_npabs (iu)
      end if
!
      if (npabs_sum(nd) > 0) then
        tau_alpha(nd) = real(npabs_cyls(nd))/real(npabs_a(nd))
      else
        write (*,*) 'Error -- no absorbed photons: nd, nploop = ', nd, nploop
        write (9,*) 'Error -- no absorbed photons: nd, nploop = ', nd, nploop
        call adios (2)
      end if
!
!  Check Convergence
!
     if (tau_alpha(nd) > 0.) then
       tol_check = zz*sqrt(abs((1.-tau_alpha(nd))/(tau_alpha(nd)* &
                   npabs_a(nd))))
     else if (tau_alpha(nd) < 0.) then
       write (*,3070) nd, nploop, tau_alpha(nd)
       write (9,3070) nd, nploop, tau_alpha(nd)
       call adios (2)
     else
       tol_check = 1.d10
     end if
     if (tol_check <= tol) then
       write (*, 2000) nd, nploop, npabs_tot(nd), tol_check, tol
       nd_converged = nd_converged + 1
       go to 600
     else
       write (*, 2001) nd, nploop, npabs_tot(nd), tol_check, tol
     end if
   end do
   nd_not_converged = nd_not_converged + 1
 600 write (*,*) 'Finished Run for Direction No. - ', nd
!
!  compute final photon quantities for direction nd
!
      if (iglass == 0) then
        iam_glass(nd) = 0
      else
        iam_glass(nd) = real(sum(npabs(int0_glass:int1_glass,nd)))/ &
                        real(npabs_sum(nd))
      end if
!
!  Output for direction nd
!
   write (9, 2050) nd, npabs_tot(nd), npabs_cyls(nd), npabs_a(nd), &
       npabs_glass(nd), tau_alpha(nd),iam_glass(nd),tol_check
!
!  Check for Lost Photons
!
   if (npabs_tot(nd) /= npdo*num_photons) then
     write (*,5000) nd, npabs_tot(nd), nplost(nd), npdo*num_photons
     write (9,5000) nd, npabs_tot(nd), nplost(nd), npdo*num_photons
   end if
 end do
!
!  Calculate IAM's
!
  Do nd = 1, num_dirs
    iam_cyls(nd) = tau_alpha(nd)/tau_alpha(norm)
  end do
!
!  Final Photon Tally
!
  npabs_run = sum(npabs_tot) 
  nplost_run = sum(nplost)
  glass_frac = real(sum(npabs_glass))/real(npabs_run)
  write (*, 2080) npabs_run, nplost_run, glass_frac
  write (9, 2080) npabs_run, nplost_run, glass_frac
  call tau_alpha_d
!
  write (*,*) 'Direction No.  tau-alpha      __IAM__  npabs_cyls  npabs_glass' 
  do nd = 1, num_dirs
    write (*,2021) nd,tau_alpha(nd),iam_cyls(nd),npabs_cyls(nd),npabs_glass(nd)
  end do
!
!  Write IAM table output
!
  do n = 0, n_long - 1
    write (2,1001) (iam_cyls(m), m = n*n_tran + 1, (n + 1)*n_tran)
  end do
!
  write (2,*) ' '
  write (2,1002) 'tau_alpha normal                    = ', tau_alpha(norm)
  write (2,1002) 'tau_alpha diffuse                   = ', tad
  write (2,1002) 'aperture area [m^2]                 = ', areaa/10000
  write (2,1002) 'backplane area/aperture area        = ', areab/areaa
  write (2,1002) 'backplane-tube center distance [cm] = ', z0(int0_cyl)
  write (2,1002) 'backplane diffuse reflectance       = ', rho_d(matno_int(2)) 
  write (2,1002) 'backplane specular reflectance      = ', rho_s(matno_int(2))
  write (2,1002) 'backplane semi-specular reflectance = ', rho_ss(matno_int(2))
  write (2,1002) 'backplane r_d power                 = ', r_d(matno_int(2))
  write (2,1002) 'backplane r_ss power                = ', r_ss(matno_int(2))
!
!  ERRORS in the indices.  Fixed by Jeff Miller, 11/6/97
!
!  write (2,1002) 'glass thickness*extinction coef.    = ', &
!               t(matno_glass(int0_cyl))*k_glass(matno_glass(int0_cyl))
!  write (2,1002) 'index of refraction                 = ', &
!               n_glass(matno_glass(int0_cyl))
!  write (2,1002) 'outside radius of glass             = ', ro(int0_cyl)
!  write (2,1002) 'radius of collector cylinder        = ', ri(int0_cyl)
!  write (2,1005) 'length of cylinders                 = ', length_cyl(int0_cyl)
  write (2,1002) 'glass thickness*extinction coef.    = ', &
               t(matno_int(int0_glass))*k_glass(matno_int(int0_glass))
  write (2,1002) 'index of refraction                 = ', &
               n_glass(matno_int(int0_glass))
  write (2,1002) 'outside radius of glass             = ', ro(1)
  write (2,1002) 'radius of collector cylinder        = ', ri(1)
  write (2,1005) 'length of cylinders                 = ', length_cyl(1)
  write (2,1002) 'tolerance                           = ', tol
  write (2,1004) '# of directions converged on        = ', nd_converged
  write (2,*) ' '
  write (2,*) 'Tube Centers:'
  write (2,1003) (y0(m), m = 1, num_cyls)
!
!  finished solution phase - end of run
!
  call adios (1)
!
!  trap errors opening files
!
  4000 continue
  write (*,4001) iu
  call adios (2)
!
!  format statements
!
 1000 format (40a1)
 1001 format (18(f5.3,1x))
 1002 format (a38, f7.4)
 1003 format (15(f6.2,1x))
 1004 format (a38, i3)
 1005 format (a38, f7.2)
 2000 format (1x,'    Converged - ', 1x, 2(i8,1x), i9, 4x, f10.6, 2x, f9.3)
 2001 format (1x,'Not Converged - ', 1x, 2(i8,1x), i9, 4x, f10.6, 2x, f9.3)
 2003 format (////18x,'S o l u t i o n   P h a s e'/// &
  / 20x,' direction photon  tot. no.   confidence   tolerance' &
  /  20x,'    no.     loop      of       interval' &
  /  20x,'             no.    photons')
 2004 format (////21x,'I n p u t   P h a s e'/// )
 2010 format (//20x,' I A M   P R O G R A M'// &
  / 5x,'Problem heading       = ',40a1 / &
  / 5x,'Date of run           = ',a8 / &
  / 5x,'Time of Run           = ',a8 / &
  / 5x,'IAM program version   = ',a8 / &
  / 5x,'Compile date          = ',a8 )
 2011 format (////12x,'P h o t o n   C o n t r o l   I n f o .'/ &
  / 5x,'Number of photons per photon loop     = ',i9 &
  / 5x,'Maximum number of photon loops        = ',i9 &
  / 5x,'Convergence tolerance for photons     = ',f9.3 &
  / 5x,'Initial seed for RNG                  = ',i9 &
  / 8x,' if 0 entered, generated ''randomly'' from the time'/)
 2012 format (////13x,'O t h e r   C o n t r o l   I n f o .'/ &
  / 5x,'Number of directions            = ',i5 &
  / 5x,'Number of cylinders             = ',i5 &
  / 5x,'Number of planar surfaces       = ',i5 &
  / 5x,'Number of Surface Materials     = ',i5 &
  / 5x,'Number of Glass Materials       = ',i5/)
 2013 format (////18x,'G e o m e t r y   I n f o .'/ &
  / 5x,'Aperture Plane' &
  / 8x,'xmin           = ',f10.4 &
  / 8x,'xmax           = ',f10.4 &
  / 8x,'ymin           = ',f10.4 &
  / 8x,'ymax           = ',f10.4 &
  / 5x,'Back Plane' &
  / 8x,'xmin           = ',f10.4 &
  / 8x,'xmax           = ',f10.4 &
  / 8x,'ymin           = ',f10.4 &
  / 8x,'ymax           = ',f10.4 &
  / 8x,'material no.   = ',i10/)
 2014 format (////16x,'P r i n t i n g   C o n t r o l' / &
  / 5x,'Print control (iprint)                 = ',i5, &
  / 8x,' =  1,  modest printing of results in output file' &
  / 8x,' =  0,      no printing of results in output file' &
  / 8x,' = -1, massive printing of results in output file' &
  / 5x,'Print control for RNG (iprint_ran)     = ',i5, &
  / 5x,'Debugging flag (idebug)                = ',i5 &
  / 8x,' =  1, modest debug info. written' &
  / 8x,' =  0, no debug info. written' &
  / 8x,' = -1, massive debug info. written' &
  / 5x,'Data check flag (idata)                = ',i5 &
  / 8x,'/=  0, data check only' &
  / 8x,' =  0, normal execution' &
  / 5x,'Glass only flag (iglass)               = ',i5 &
  / 8x,' =  0, no glass is around any cylinder' &
  / 8x,' =  1, glass is around at least one cylinder' //)
 2015 format(//5x,'I n t e r s e c t i o n   I n d i c e s   F o ll o w'/ &
  / 5x,'top plane                                       = ',4x,'0' &
  / 5x,'bottom, outside backplane                       = ',4x,'1' &
  / 5x,'back plane (outside aperture plane)             = ',4x,'2' &
  / 5x,'Cylinders:' &
  / 8x,'start of cylinders (int0_cyl)                = ', i5 &
  / 8x,'  end of cylinders (int1_cyl)                = ', i5 &
  / 5x,'Glass:' &
  / 8x,'start of glass (int0_glass)                  = ', i5 &
  / 8x,'  end of glass (int1_glass)                  = ', i5 &
  / 5x,'Surfaces:' &
  / 8x,'start of surfaces (int0_surf)                = ', i5 &
  / 8x,'  end of surfaces (int1_surf)                = ', i5 //)
 2021 format(i8, 9x, f7.5, 7x, f7.5, i9, i9)
!
 2050 format (//5x,'Output for direction  -',i3/ &
  / 10x,'npabs_tot    = ',i9 &
  / 10x,'npabs_cyls   = ',i9 &
  / 10x,'npabs_a      = ',i9 &
  / 10x,'npabs_glass  = ',i9 &
  / 10x,'tau_alpha    = ',f9.6 &
  / 10x,'iam_glass    = ',f9.6 &
  / 10x,'tol_check    = ',f9.6/)
 2051 format (5x,'Photon absorption numbers: ' &
  / 5x,'npabs_bot          = ',i10 &
  / 5x,'npabs_back         = ',i10 &
  / 5x,'npabs_ap           = ',i10 &
  / 5x,'npabs_cyl          = ',i10 &
  / 5x,'npabs_top          = ',i10/)
 2060 format (5x,'Photon absorption IAM''s: ' &
  / 5x,'iam_cyl           = ',f12.6 &
  / 5x,'iam_glass_o       = ',f12.6 &
  / 5x,'iam_glass_i       = ',f12.6 &
  / 5x,'iam_glass         = ',f12.6/)
 2080 format (5x,'Photon numbers for the entire run : ' &
  / 10x,'Absorbed photons           = ',i10 &
  / 10x,'Lost photons               = ',i10 &
  / 10x,'Fraction absorbed by glass = ',f10.6 //)
 3070 format (//2x,'Error in tau_alpha for nd, nploop = ',2i5, &
  2x,'tau_alpha = ',f10.5)
  4001 format (1x,'Error opening file - ',i5,'  - stopping')
  5000 format (1x,'Warning -- truly lost photons for nd - ', i5 &
  / 5x,'Number absorbed               = ',i10 &
  / 5x,'Number lost                   = ',i10 &
  / 5x,'Number emitted                = ',i10 )
 7777 format (//10x,'Circular Constants'// &
 / 5x,'pi                         = ',f10.6 &
 / 5x,'p2i                        = ',f10.6 &
 / 5x,'pi2                        = ',f10.6 &
 / 5x,'pi4                        = ',f10.6 &
 / 5x,'dtor                       = ',f10.6 &
 / 5x,'rtod                       = ',f10.6 //)
!
end program iam
!
!*********************************************************************
!
subroutine adios (k)
!
!*********************************************************************
!*     subroutine to terminate execution                             *
!*********************************************************************
!
  use kinds
  use control
  use dates
  use photons
!
  implicit none
!
  integer(kind=kindi) :: k, imon, idays, ihrs, imins, isecs
  integer, dimension(12) :: mondays =  &
      (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
  real(kind=kindr)  :: seconds, rate
!
  call getdt (2)
!
  call header (9)
  write (9, 1000)
!
  if (k == 1) then
    write (*,*) 'Normal Termination'
    write (9,*) 'Normal Termination'
  else
    write (*,*) 'Error Termination - see output file'
    write (9,*) 'Error Termination'
  endif
!
!  compute time and execution rate
!   (for now, based on wall clock time - no seconds function yet)
!
  idays = valuesa(3,2) - valuesa(3,1)
  if (idays < 0) then
    imon = valuesa(2,1)
    idays = mondays(imon) - valuesa(3,1) + valuesa(3,2)
  end if
!
  if (idays < 0) then
    write (9, *) 'Error in days --  setting to zero'
    idays = 0
  end if
  ihrs = valuesa(5,2)-valuesa(5,1) + idays*24
  imins = valuesa(6,2)-valuesa(6,1)
  isecs = valuesa(7,2)-valuesa(7,1)
  write (9, 1001) ihrs, imins, isecs
!
!  compute rate, based on wall clock time
!
  seconds = ((idays*24. + ihrs)*60. + imins)*60. + isecs
  if (seconds > 0.) then
    rate = real(sum(npabs_tot(1:num_dirs)))/seconds
  else
    rate = 0.
  end if
  write (9, 1002) seconds, npabs_run, nplost_run, rate
!
  stop
!
 1000 format (////5x,'R u n   T e r m i n a t i o n   I n f o r m a t i o n'//)
 1001 format (////5x,'Elapsed Wall Clock Times' / &
  / 8x,'Hours                = ', i5 &
  / 8x,'Minutes              = ', i5 &
  / 8x,'Seconds              = ', i5 //)
 1002 format (//5x,'Execution Rate' / &
  / 8x,'Total seconds                 = ', f10.2, ' (seconds)' &
  / 8x,'Total No. of Photons Emitted  = ', i10 &
  / 8x,'Total No. of Photons Lost     = ', i10 &
  / 8x,'Execution Rate                = ', f10.2, ' (photons/sec)')
!
 end subroutine adios
!
!***********************************************************************
!
subroutine getdt (k)
!
!***********************************************************************
!*    Obtain ASCII date and time, and 'random' initial seed            *
!***********************************************************************
!
  use kinds
  use random
  use control
  use dates
!
  implicit none
!
  character(len=5) :: zone
  character(len=10) :: time
  character(len=8) :: date
!
  integer(kind=kindi) :: iyr, imon, iday, ihr, imin, isec, imil
  integer :: values(8), k
!
  call date_and_time (date, time, zone, values)
  valuesa(1:8, k) = values
  if (k == 0) then
!
!     get ASCII date and ASCII Time
!
    adate(1:2) =  date(5:6)
    adate(3:3) = '-'
    adate(4:5) =  date(7:8)
    adate(6:6) = '-'
    adate(7:8) =  date(3:4)
!
    atime(1:2) = time(1:2)
    atime(3:3) = ':'
    atime(4:5) = time(3:4)
    atime(6:6) = ':'
    atime(7:8) = time(5:6)
!
!  get integer values
!
    iyr =  values(1)
    if (iyr > 99) then
      iyr = mod(iyr, 100)
    end if
    imon = values(2)
    iday = values(3)
    ihr =  values(5)
    imin = values(6)
    isec = values(7)
    imil = values(8)
!
!     get initial seed, seed_from_time
!
    seed_from_time=isec+60*(imin+60*(ihr+24*(iday+28*(imon+iyr))))
    if (mod(seed_from_time,2) == 0) then
       seed_from_time = ior(seed_from_time, 1)
    end if
!
    if (seed_from_time < 0) then
       write (*,*) 'seed_from_time is negative = ',seed_from_time,' correcting'
       write (9,*) 'seed_from_time is negative = ',seed_from_time,' correcting'
       seed_from_time = seed_from_time + modulus
    end if
  end if
!
  return
!
 1000 format (5x,'  Date and Time Follow for k = ', i5/ &
 / 5x,'iyr                   = ', i5 &
 / 5x,'imon                  = ', i5 &
 / 5x,'iday                  = ', i5 &
 / 5x,'ihr                   = ', i5 &
 / 5x,'imin                  = ', i5 &
 / 5x,'isec                  = ', i5 &
 / 5x,'imil                  = ', i5 &
 / 5x,'seed_from_time        = ', i20//)
!
end subroutine getdt
!
!*********************************************************************
!
subroutine header (iu)
!
!*********************************************************************
!*     subroutine to print header to file iu                         *
!*********************************************************************
!
  use kinds
  use control
!
  implicit none
!
  integer ::  iu, n
  character(len=8) :: stars8='********'
  character(len=7) :: stars7='*******'
!
  write (iu, 1000) (stars8, n=1,9), stars7
  write (iu, 1001) (head(n), n=1,72)
  write (iu, 1002) (stars8, n=1,9),stars7
!
  return
!
 1000 format ('1',9a8,a7)
 1001 format (4x, 72a1)
 1002 format (' ',9a8,a7)
!
end subroutine header
!
!***********************************************************************
!***********************************************************************
!**                                                                   **
!**                   Input Phase Routines Follow                     **
!**                                                                   **
!***********************************************************************
!***********************************************************************
!
!
!*********************************************************************
!
subroutine check_control
!
!*********************************************************************
!*     subroutine to check control information                       *
!*********************************************************************
!
  use kinds
  use control
  use materials
  use planes
!
  implicit none
!
  integer(kind=kindi) :: ierr = 0
!
  if (num_photons == 0) then
    num_photons = 1000
   endif
   if (num_ploops <= 0) then
     num_ploops = 1
   endif
   if (tol < 1.e-6) then
     tol = 0.01
   endif
   if (iglass < 0 .or. iglass > 1) then
     write (*,*) 'Input Error -- iglass = ',iglass,' out of range [0,1]'
     write (9,*) 'Input Error -- iglass = ',iglass,' out of range [0,1]'
     ierr = -1
   endif
   if (num_dirs <= 0) then
     write (*,*) 'Input Error -- number of directions must be > 0'
     write (9,*) 'Input Error -- number of directions must be > 0'
     ierr = -1
   endif
   if (num_cyls <= 0) then
     write (*,*) 'Input Error -- number of cylinders must be >= 0'
     write (9,*) 'Input Error -- number of cylinders must be >= 0'
     ierr = -1
   endif
   if (num_mats_surf <= 0) then
     write (*,*) 'Input Error -- no. of surface materials must be > 0'
     write (9,*) 'Input Error -- no. of surface materials  must be > 0'
     ierr = -1
   endif
!
!  check planar geometries
!
   if (xminb  >= xmaxb) then
     write (*,*) 'Input Error -- xminb must be < xmaxb'
     write (9,*) 'Input Error -- xminb must be < xmaxb'
     ierr = -1
   endif
   if (yminb  >= ymaxb) then
     write (*,*) 'Input Error -- yminb must be < ymaxb'
     write (9,*) 'Input Error -- yminb must be < ymaxb'
     ierr = -1
   endif
!
   if (xmina  >= xmaxa) then
     write (*,*) 'Input Error -- xmina must be < xmaxa'
     write (9,*) 'Input Error -- xmina must be < xmaxa'
     ierr = -1
   endif
   if (ymina  >= ymaxa) then
     write (*,*) 'Input Error -- ymina must be < ymaxa'
     write (9,*) 'Input Error -- ymina must be < ymaxa'
     ierr = -1
   endif
!
  if (xmina < xminb .or. xmaxa > xmaxb) then
     write (*,*) 'Input Error -- x coordinates for aperture plane must be '
     write (*,*) '  within x coordinates for backplane'
     write (9,*) 'Input Error -- x coordinates for aperture plane must be '
     write (9,*) '  within x coordinates for backplane'
     ierr = -1
   endif
  if (ymina < yminb .or. ymaxa > ymaxb) then
     write (*,*) 'Input Error -- y coordinates for aperture plane must be '
     write (*,*) '  within y coordinates for backplane'
     write (9,*) 'Input Error -- y coordinates for aperture plane must be '
     write (9,*) '  within y coordinates for backplane'
     ierr = -1
   endif
!
  if (matnob < 1 .or. matnob > num_mats_surf) then
     write (*,*) 'Input Error -- Material no. for aperture plane must be'
     write (*,*) '  > 0 and <=', num_mats_surf
     write (9,*) 'Input Error -- Material no. for aperture plane must be'
     write (9,*) '  > 0 and <=', num_mats_surf
     ierr = -1
   endif
!
  if (ierr == -1) then
    write (*,*) 'Stopping in check_control -- errors detected'
    write (9,*) 'Stopping in check_control -- errors detected'
    call adios (2)
  end if
!
  return
!
end subroutine check_control
!
!*********************************************************************
!
subroutine files (unit_in, unit_out, nlines)
!
!*********************************************************************
!*     subroutine to strip comment (& in col. 1) lines from files    *
!*********************************************************************
!
  use kinds
!
  implicit none
!
  integer(kind=kindi) :: unit_out, unit_in, ierr, nlines, nchars, n, i
  integer(kind=kindi), parameter :: nlines_max = 10000
  character(len=1) :: aline(80)
!
  do n = 1, nlines_max
    read (unit_in, 1000, end = 100, err = 200) aline
    if ( aline(1) /= '&') then
      do i = 80, 1, -1
         if (aline(i) /= ' ') then
           go to 10
         end if
      end do
 10   nchars = i + 1
      write (unit_out, 1000) aline(1:nchars)
    end if
  end do
!
  write (*,*) 'Error from files - too many lines in file no. - ', unit_in
  write (*,*) ' --- Stopping ---'
  call adios (2)
!
 100 nlines = n - 1
  rewind (unit_out)
  return
!
 200 write (*,*) 'Error reading line ',n,' from file no. - ', unit_in
     write (9,*) 'Error reading line ',n,' from file no. - ', unit_in
  ierr = 1
  call adios (2)
!
 1000 format (80a1)
!
end subroutine files
!
!***********************************************************************
!
subroutine raninit (iprnt, seed0)
!
!***********************************************************************
!*    subroutine to initialize pseudo-random number generator          *
!***********************************************************************
!
  use kinds
  use random
!
  implicit none
!
  integer(kind=kindi) :: n, iprnt, seed0
  integer(kind=kindi), parameter :: ia_mult = 16807
!
!  generate initial seeds using multiplicative, congruential generator
!
  seed(1) = mod(ia_mult*seed0,modulus)
  do n = 2, 17
     seed(n) = mod(ia_mult*seed(n-1),modulus)
     if (seed(n) < 0) then
        seed(n) = seed(n) + modulus
     endif
  end do
!
!     check for an odd value
!
  do n=1, 17
    if (mod(seed(n),2) == 1) then
      go to 250
    endif
  end do
  seed(1) = ior (seed(1), 1)
!
!     output of initial variables
!
 250 continue
  if (iprnt == 0 ) then
     return
  endif
  call header (9)
  write (9,1000) 
  n = 0
  write (9,1001) n,seed0,seed0
  do n=1, 17
     write (9,1001) n,seed(n),seed(n)
  end do
!
  return
!
 1000 format (//10x,'i n i t i a l   s e e d s'/ &
     / 7x,'n',12x,'seed',8x,'seed' &
     / 20x,'(dec)',8x,'(hex)') 
 1001 format (5x,i5,5x,i12,5x,z8)
!
end subroutine raninit
!
!*********************************************************************
!
subroutine read_cylinders
!
!*********************************************************************
!*     subroutine to read cylinder input                             *
!*********************************************************************
!
  use control
  use cylinders
  use constants
  use planes
!
  implicit none
!
  integer(kind=kindi) :: n, nn, ierr, mg, mn, nc
  integer(kind=kindi), allocatable, dimension (:) :: itest
  real(kind=kindr) x0n, y0n, z0n, ln, ron, rin, betadn
!
  write (*,*) '  - reading cylinders'
  call header (9)
  write (9, 1000)
!
  ierr = 0
!
!  allocate itest array
!
  allocate (itest(num_cyls))
  itest = 0
!
  do n = 1, num_cyls
    read (1, *, end = 3000, err = 4000) nn, x0n, y0n, z0n, ln, ron, &
      rin, betadn, mg, mn
!
    if (nn < 1 .or. nn > num_cyls) then
      ierr = -1
      write (*,*) 'Error --  reading cylinders on line - ',n
      write (*,*) 'Cylinder no. -- ',nn,' out of range'
      write (9,*) 'Error --  reading cylinders on line - ',n
      write (9,*) 'Cylinder no. -- ',nn,' out of range'
    else
      xmax_cyl(nn) = x0n
      y0(nn) = y0n
      z0(nn) = z0n
      length_cyl(nn) = ln
      ro(nn) = ron
      ri(nn) = rin
      betad(nn) = betadn
      matno_glass(nn) = mg
      matno_cyli(nn) = mn
    end if
!
    if (itest(nn) /= 0) then
      ierr = 1
      write (*,*) 'Error --  reading cylinders on line - ',n
      write (*,*) 'Cylinder no. -- ',nn,' already entered'
      write (9,*) 'Error --  reading cylinders on line - ',n
      write (9,*) 'Cylinder no. -- ',nn,' already entered'
    else
      itest(nn) = 1
    end if
  end do
!
!  test if all cylinders have been entered
!
  do n = 1, num_cyls
    if (itest(n) == 0) then
      write (*,*) 'Error -- in cylinders entered' 
      write (*,*) 'Cylinder no. -- ',n,' is missing'
      write (9,*) 'Error -- in cylinders entered' 
      write (9,*) 'Cylinder no. -- ',n,' is missing'
      ierr = -1
    end if
  end do
!
!  Compute maximum and minimum extents of cylinders
!  - cylinder (array) quantities
!
  xmin_cyl = xmax_cyl - length_cyl
  ymax_cyl = y0 + ro
  ymin_cyl = y0 - ro
  zmax_cyl = z0 + ro
  zmin_cyl = z0 - ro
!
!  global max, min quantities
!
  xmax = maxval (xmax_cyl)
  ymax = maxval (ymax_cyl)
  zmax_cyls = maxval (zmax_cyl)
  zmax = zmax_cyls
  xmin = minval (xmin_cyl)
  ymin = minval (ymin_cyl)
  zmin = minval (zmin_cyl)
!
! Check extents of cylinders vs. aperture plane
!
  if (xmax > xmaxa .or. xmin < xmina) then
    write (*,*) 'Error -- cylinders entend in x beyond aperture plane'
    write (9,*) 'Error -- cylinders entend in x beyond aperture plane'
    ierr = -1
  end if
  if (ymax > ymaxa .or. ymin < ymina) then
    write (*,*) 'Error -- cylinders entend in y beyond aperture plane'
    write (9,*) 'Error -- cylinders entend in y beyond aperture plane'
    ierr = -1
  end if
!
!  process and check input
!
  do nc = 1, num_cyls
    if (ro(nc) <= ri(nc)) then
      write (*,*) 'Error -- ro(',nc,') < ri(',nc,')'
      write (9,*) 'Error -- ro(',nc,') < ri(',nc,')'
      ierr = -1
    end if
    if (ro(nc) >= z0(nc)) then
      write (*,*) 'Error -- ro(',nc,') >= z0(',nc,')'
      write (9,*) 'Error -- ro(',nc,') >= z0(',nc,')'
      ierr = -1
    end if
    if (length_cyl(nc) <= 0.) then
      write (*,*) 'Error -- length_cyl(',nc,') <= 0.'
      write (9,*) 'Error -- length_cyl(',nc,') <= 0.'
      ierr = -1
    end if
!
    if (betad(nc) /= 0.) then
      write (*,*) 'Error -- betad(',nc,') must = 0 for this version'
      write (9,*) 'Error -- betad(',nc,') must = 0 for this version'
      ierr = -1
    end if
!
    if (iglass /= 0) then
      if (matno_glass(nc) < 1 .or. matno_glass(nc) > num_mats_glass) then
        write (*,*) 'Error -- matno_glass(',nc,') = ',matno_glass(nc), & 
                     ': out of range'
        write (9,*) 'Error -- matno_glass(',nc,') = ',matno_glass(nc), &
                     ': out of range'
        ierr = -1
      end if
    end if
    if (matno_cyli(nc) < 1 .or. matno_cyli(nc) > num_mats_surf) then
      write (*,*) 'Error -- matno_cyli(',nc,') = ',matno_cyli,': out of range'
      write (9,*) 'Error -- matno_cyli(',nc,') = ',matno_cyli,': out of range'
      ierr = -1
    end if
!
!  conversions
!
    beta_cyl(nc) = betad(n)*dtor
  end do
!
  do nc = 1, num_cyls
    write (9, 1011) nc, xmax_cyl(nc), y0(nc), z0(nc), length_cyl(nc), ro(nc)
  end do
  if (iglass /= 0) then
    write (9,1001)
  else
    write (9,1101)
  end if
  do nc = 1, num_cyls
    if (iglass /= 0) then
      write (9, 1012) nc, ri(nc), betad(nc), matno_glass(nc), matno_cyli(nc)
    else
      write (9, 1112) nc, betad(nc), matno_cyli(nc)
    end if
  end do
  write (9,1002)
  do n = 1, num_cyls
    write (9, 1013) n, xmin_cyl(n), xmax_cyl(n), ymin_cyl(n), &
      ymax_cyl(n), zmin_cyl(n), zmax_cyl(n)
  end do
!
  deallocate (itest)
!
  if (ierr == -1) then
    write (*,*) 'Stopping in read_cylinders -- errors detected'
    write (9,*) 'Stopping in read_cylinders -- errors detected'
    call adios (2)
  end if
!
  return
!
 3000 write (*,*) 'Error from read_cylinders -- ', &
   ' reached end of file at line - ',n
  call adios (2)
!
 4000 write (*,*) 'Error from read_cylinders -- ', &
   ' error reading input on line - ',n
  call adios (2)
!
 1000 format (//13x,' O u t p u t   o f   C y l i n d e r s'// &
  4x,' No.',4x,'x0',10x,'y0',10x,'z0',10x,'len',10x,'ro'/)
 1001 format (//4x, ' No.',4x,'ri',7x,'beta(deg)', 6x,'matno_glass',2x, &
  'matno_cyli'/)
 1101 format (//11x,6x,'beta(deg)', 4x,'matno_cyli'/)
 1002 format (//11x,'Mins and Maxs of Cylinders'/ &
  / 5x,'No.',7x,'x-min',7x,'x-max',7x,'y-min',7x,'y-max',7x,'z_min',7x,'z-max'/)
 1011 format (2x,i5,'.',7(1x,g11.4), 2x,g9.2,3x,i5,8x,i5)
 1012 format (2x,i5,'.',2(1x,g11.4),5x,i5,7x,i5)
 1112 format (2x,i5,'.',1(1x,g11.4),1(5x,i5))
 1013 format (2x,i5,'.',6(2x,f10.2))
!
end subroutine read_cylinders
!
!*********************************************************************
!
subroutine read_directions
!
!*********************************************************************
!*     subroutine to read and convert input directions               *
!*********************************************************************
!
  use control
  use directions
  use constants
!
  implicit none
!
  integer(kind=kindi) :: n, nl, nt, ierr, nn
  real(kind=kindr) :: long, tran
  integer(kind=kindi), allocatable, dimension(:) :: ltest, ttest
!
  write (*,*) '  - reading directions'
  call header (9)
  write (9,1000)
!
  ierr = 0
!
!  allocate test arrays
!
  allocate (ltest(n_long), ttest(n_tran))
  ltest = 0
  ttest = 0
!
! Read Longitudinal Incident Angles
!
  do n = 1, n_long
    read (1, *, end = 3000, err = 4000) nn, long
!
    if (nn < 1 .or. nn > n_long) then
      ierr = -1
      write (*,*) 'Error --  reading thetald on line - ',n
      write (*,*) 'thetald no. -- ',nn,' out of range'
      write (9,*) 'Error --  reading thetald on line - ',n
      write (9,*) 'thetald no. -- ',nn,' out of range'
    else
      thetald(nn) = long  
    end if
!
!
    if (ltest(nn) /= 0) then
      ierr = -1
      write (*,*) 'Error --  reading thetald on direction - ',n
      write (*,*) 'thetald no. -- ',nn,' already entered'
      write (9,*) 'Error --  reading thetald on direction - ',n
      write (9,*) 'thetald no. -- ',nn,' already entered'
    else
      ltest(nn) = 1
    end if
  end do
!
  do n = 1, n_tran
    read (1, *, end = 3000, err = 4000) nn, tran
!
    if (nn < 1 .or. nn > n_tran) then
      ierr = -1
      write (*,*) 'Error --  reading thetatd on line - ',n
      write (*,*) 'thetatd no. -- ',nn,' out of range'
      write (9,*) 'Error --  reading thetatd on line - ',n
      write (9,*) 'thetatd no. -- ',nn,' out of range'
    else
      thetatd(nn) = tran
    end if
!
!
    if (ttest(nn) /= 0) then
      ierr = -1
      write (*,*) 'Error --  reading thetatd on direction - ',n
      write (*,*) 'thetatd no. -- ',nn,' already entered'
      write (9,*) 'Error --  reading thetatd on direction - ',n
      write (9,*) 'thetatd no. -- ',nn,' already entered'
    else
      ttest(nn) = 1
    end if
  end do
!
!  test if all directions have been entered
!
  do n = 1, n_long
    if (ltest(n) == 0) then
      write (*,*) 'Error -- in thetald entered' 
      write (*,*) 'thetald no. -- ',n,' is missing'
      write (9,*) 'Error -- in thetald entered' 
      write (9,*) 'thetald no. -- ',n,' is missing'
      ierr = -1
    end if
  end do
!
  do n = 1, n_tran
    if (ttest(n) == 0) then
      write (*,*) 'Error -- in thetatd entered'
      write (*,*) 'thetatd no. -- ',n,' is missing'
      write (9,*) 'Error -- in thetatd entered'
      write (9,*) 'thetatd no. -- ',n,' is missing'
      ierr = -1
    end if
  end do
!
  if (ierr == -1) then
    write (*,*) 'Stopping in read_directions - errors detected'
    write (9,*) 'Stopping in read_directions - errors detected'
    call adios (2)
  end if
!
!  process and check input
!
  norm = 0
  n = 0
  do nl = 1, n_long
    if (thetald(nl) < 0. .or. thetald(nl) >= 90.) then
      write (*,*) 'Error -- thetald(',nl,') out of range: [0, 90) deg.'
      write (9,*) 'Error -- thetald(',nl,') out of range: [0, 90) deg.'
      ierr = -1
    end if
!
    thetal(nl) = thetald(nl)*dtor
    ttl(nl) = tan(thetal(nl))
!
    do nt = 1, n_tran
      n = n + 1
      if (thetatd(nt) <= -90. .or. thetatd(nt) >= 90.) then
        write (*,*) 'Error -- thetatd(',nt,') out of range: (-90, 90) deg.'
        write (9,*) 'Error -- thetatd(',nt,') out of range: (-90, 90) deg.'
        ierr = -1
      end if
!
      if (ierr == -1) then
        write (*,*) 'Stopping in read_directions - errors detected'
        write (9,*) 'Stopping in read_directions - errors detected'
        call adios (2)
      end if
!
! conversions
!
      thetat(nt) = thetatd(nt)*dtor
      ttt(nt) = tan(thetat(nt))
      if (thetal(nl) == 0) then
        if (thetat(nt) == 0) then
          norm = n
        end if
        if (thetatd(nt) > 0) then
          phi(n) = 90*dtor
          altitude(n) = (90 - thetatd(nt))*dtor
        else
          phi(n) = -90*dtor
          altitude(n) = (90 + thetatd(nt))*dtor
        end if
        sp(n) = sin(phi(n))
      elseif (thetat(nt) == 0) then
        phi(n) = 0.
        sp(n) = sin(phi(n))
        altitude(n) = (90 - thetald(nl))*dtor
      else
        phi(n) = atan(ttt(nt)/ttl(nl))
        sp(n) = sin(phi(n))
        altitude(n) = atan(sp(n)/ttt(nt))
      end if
!
      ca(n) = cos(altitude(n))
      sa(n) = sin(altitude(n))
      ta(n) = tan(altitude(n))
      cp(n) = cos(phi(n))
      tp(n) = tan(phi(n))
      altituded(n) = altitude(n)*rtod
      phid(n) = phi(n)*rtod
      exa(n) = ca(n)*cp(n)
      eya(n) = -ca(n)*sp(n)
      eza(n) = sa(n)
    end do
  end do
!
! make sure normal incidence was entered
!
  if (norm == 0) then
    write (*,*) 'Stopping in read_directions - normal incidence needed'
    write (9,*) 'Stopping in read_directions - normal incidence needed'
    call adios (2)
  end if
!
  deallocate (ltest, ttest)
!
!  Output
!
  do n = 1, num_dirs
    write (9,1001) n, altituded(n), phid(n), exa(n), eya(n), eza(n)
  end do
! do n = 1, n_long
!   write (2,1002) thetald(n)
! end do
! do n = 1, n_tran
!   write (2,1003) thetatd(n)
! end do
      write (2, 1002) (thetald(n), n = 1, n_long)
      write (2, 1002) (thetatd(n), n = 1, n_tran)
!
  if (ierr == -1) then
    write (*,*) 'Errors detected in read_directions - stopping '
    write (9,*) 'Errors detected in read_directions - stopping '
    call adios (2)
  end if
!
  return
!
 3000 write (*,*) 'Error in read_directions -- reached end of file -- Stopping'
      write (9,*) 'Error in read_directions -- reached end of file -- Stopping'
  call adios (2)
  return
!
 4000 write (*,*) 'Error from read_directions -- ', &
   ' error reading input on line - ',n,' -- stopping'
      write (9,*) 'Error from read_directions -- ', &
   ' error reading input on line - ',n,' -- call stopping'
  call adios (2)
!
 1000 format (////22x,'D i r e c t i o n s'/ &
  / 5x,'No.',9x,'Altitude',8x,'Azimuth',6x,'  EX',10x,'EY',10x,'EZ' &
  / 18x,'(deg.)',10x,'(deg.)'/)
 1001 format (2x,i5,'.',6x,f10.3,5x,f10.3,1x,3(1x,f11.4))
 1002 format (15(f5.0, 1x))
!
end subroutine read_directions
!
!*********************************************************************
!
subroutine read_materials
!
!*********************************************************************
!*     subroutine to read materials from input file                  *
!*********************************************************************
!
  use kinds
  use constants
  use control
  use materials
!
  implicit none
!
  integer(kind=kindi) :: n, ierr, nn, it
  integer(kind=kindi), allocatable, dimension (:) :: itest
  real(kind=kindr) :: ng, kg, rhos, rhoss, rhod, rss, rd, tn, l
  real(kind=kindr) :: th1, th2, ct1, ct2, st1, st2, tau
  real(kind=kindr) :: rho_par, rho_perp, tau_par, tau_perp, alpha_glass
!
  write (*,*) '  - reading surface materials'
  call header (9)
  write (9,1000)
!
  ierr = 0
!
!  Initialization of PDF array
!
    pdf_surf = 0.
!
!  allocate itest array
!
  allocate (itest(num_mats_surf))
  itest = 0
!
  do n = 1, num_mats_surf
    read (1, *, end = 3000, err = 4000) nn, rhos, rhoss, rhod, rss, rd
    if (nn < 1 .or. nn > num_mats_surf) then
      ierr = -1
      write (*,*) 'Error --  reading surface materials on line - ',n
      write (*,*) 'Surface material no. -- ',nn,' out of range'
      write (9,*) 'Error --  reading surface materials on line - ',n
      write (9,*) 'Surface material no. -- ',nn,' out of range'
    end if
    rho_s(nn) = rhos
    rho_ss(nn) = rhoss
    rho_d(nn) = rhod
    r_ss(nn) = rss
    r_d(nn) = rd
!
!  Test if already entered
!
    if (itest(nn) /= 0) then
      ierr = 1
      write (*,*) 'Error --  reading surface materials on material - ',n
      write (*,*) 'Surface material no. -- ',nn,' already entered'
      write (9,*) 'Error --  reading surface  materials on material - ',n
      write (9,*) 'Surface material no. -- ',nn,' already entered'
    else
      itest(nn) = 1
    end if
  end do
!
!  test if all surfaces have been entered
!
  do n = 1, num_mats_surf
    if (itest(n) == 0) then
      write (*,*) 'Error -- in surface materials entered' 
      write (*,*) 'Surface material no. -- ',n,' is missing'
      write (9,*) 'Error -- in surface materials entered' 
      write (9,*) 'Surface material no. -- ',n,' is missing'
      ierr = -1
    end if
  end do
!
!  Check Values
!
  do n = 1, num_mats_surf
     if (rho_s(n) < 0. .or. rho_s(n) > 1.) then
      write (9, 3001) n
      write (9, 3003) rho_s(n)
     ierr = -1
    end if
    if (rho_ss(n) < 0. .or. rho_ss(n) > 1.) then
      write (9, 3001) n
      write (9, 3004) rho_ss(n)
      ierr = -1
    end if
    if (rho_d(n) < 0. .or. rho_d(n) > 1.) then
      write (9, 3001) n
      write (9, 3005) rho_d(n)
      ierr = -1
    end if
    if (r_ss(n) < 0. .or. r_ss(n) > 10.) then
      write (9, 3001) n
      write (9, 3006) r_ss(n)
      ierr = -1
    end if
    if (r_d(n) < 0. .or. r_d(n) > 10.) then
      write (9, 3001) n
      write (9, 3007) r_d(n)
      ierr = -1
    end if
  end do
!
!  Output of surface materials
!
  do n = 1, num_mats_surf
    pdf_surf(1,n) = rho_ss(n)
    pdf_surf(2,n) = pdf_surf(1,n) + rho_s(n)
    pdf_surf(3,n) = pdf_surf(2,n) + rho_d(n)
!
    rp1inv_ss(n) = 1.d0/(1.d0+r_ss(n))
    rp1inv_d(n) = 1.d0/(1.d0+r_d(n))
!
    write (9, 1001) n, rho_s(n), rho_ss(n), rho_d(n), r_ss(n), r_d(n)
!
!  Test range
!
    if (pdf_surf(3,n) > 1.) then
      write (*,3001) n
      write (*,3010) pdf_surf(3,n)
      ierr = -1
    end if
  end do
!
  deallocate (itest)
!
! Glass Materials
!
  if (num_mats_glass > 0) then
    allocate (itest(num_mats_glass))
    itest = 0
    do n = 1, num_mats_glass
      read (1, *, end = 3000, err = 4000) nn, ng, kg, tn
      if (nn < 1 .or. nn > num_mats_glass) then
        ierr = -1
        write (*,*) 'Error --  reading glass materials on line - ',n
        write (*,*) 'Glass material no. -- ',nn,' out of range'
        write (9,*) 'Error --  reading glass materials on line - ',n
        write (9,*) 'Glass material no. -- ',nn,' out of range'
      end if
      n_glass(nn) = ng
      k_glass(nn) = kg
      t(nn)  = tn
!
!  Test if already entered
!
      if (itest(nn) /= 0) then
        ierr = 1
        write (*,*) 'Error --  reading glass materials on material - ',n
        write (*,*) 'Glass material no. -- ',nn,' already entered'
        write (9,*) 'Error --  reading glass  materials on material - ',n
        write (9,*) 'Glass material no. -- ',nn,' already entered'
      else
        itest(nn) = 1
      end if
    end do
!
!  Check Values
!
    do n = 1, num_mats_glass
      if (n_glass(n) < 1. .or. n_glass(n) > 10.) then
        write (9, 3011) n
        write (9, 3012) n_glass(n)
        ierr = -1
      end if
      if (k_glass(n) <= 0. .or. k_glass(n) > 1.) then
        write (9, 3011) n
        write (9, 3013) k_glass(n)
        ierr = -1
      end if
      if (t(n) < 1.e-6) then
        write (*,*) 'Error -- t(',n,') < 1.e-6'
        write (9,*) 'Error -- t(',n,') < 1.e-6'
        ierr = -1
      end if
    end do
!
!  Output of glass materials
!
    call header (9)
    write (9, 1010)
    do n = 1, num_mats_glass
      do it = 0, 89
        th1 =it*dtor
        ct1 = cos(th1)
        st1 = sin(th1)
        st2 = st1/n_glass(n)
        th2 = asin(st2)
        ct2 = cos(th2)
        l = t(n)/ct2
        tau = exp(-k_glass(n)*l)
        rho_par = ((ct2 - n_glass(n)*ct1)**2)/((ct2 + n_glass(n)*ct1)**2)
        rho_perp = ((ct1 - n_glass(n)*ct2)**2)/((ct1 + n_glass(n)*ct2)**2)
        rho_par = max(rho_par, 0.d0)
        rho_perp = max(rho_perp, 0.d0)
        tau_par = ((1.-rho_par)**2)*tau/(1.-tau**2*rho_par**2)
        tau_perp = ((1.-rho_perp)**2)*tau/(1.-tau**2*rho_perp**2)
        tau_glass(it, n) = 0.5*(tau_par + tau_perp)
        rho_par = rho_par*(1. + (1.-rho_par)**2*tau**2/(1.-rho_par**2*tau**2))
        rho_perp =rho_perp*(1.+(1.-rho_perp)**2*tau**2/(1.-rho_perp**2*tau**2))
        rho_glass(it, n) = 0.5*(rho_par + rho_perp)
      end do
      tau_glass(90, n) = 0.d0
      rho_glass(90, n) = 1.d0
      write (9, 1011) n, n_glass(n), k_glass(n)
!
!  compute increment in glass properties
!
      do it = 0, 89
        drho_glass(it,n) = rho_glass(it+1,n) - rho_glass(it,n)
        dtau_glass(it,n) = tau_glass(it+1,n) - tau_glass(it,n)
      end do
      drho_glass(90,n) = 0.
      dtau_glass(90,n) = 0.
!
      if (idebug /=0) then
        write (9,1012)
        do it = 0, 90
          alpha_glass = 1.d0 - tau_glass(it,n) - rho_glass(it,n)
          write (9,1013) it, tau_glass(it,n), rho_glass(it,n), alpha_glass
        end do
      end if
    end do
!
    deallocate (itest)
  end if
!
!  Error termination
!
  if (ierr == -1) then
    write (*,*) 'Fatal errors detected in read_materials -- stopping'
    write (9,*) 'Fatal errors detected in read_materials -- stopping'
    call adios (2)
  end if
!
  return
!
 3000 write (*,*) 'Error from read_materials --', &
  ' reached end of file at material sequence - ',n    
      write (9,*) 'Error from read_materials --', &
  ' reached end of file at material sequence - ',n
  call adios (2)
!
 4000 write (*,*) 'Error from read_materials -- ', &
   ' error reading input on line - ',n
      write (9,*) 'Error from read_materials-- ', &
   ' error reading input on line - ',n
  call adios (2)
!
 1000 format (//13x,' O u t p u t   o f   M a t e r i a l s'/ &
  / 10x,'S u r f a c e   M a t e r i a l s :' &
  / 5x,'No.     rho_s    rho_ss     rho_d        r_ss     r_d')
 1001 format (2x,i5,'.',3(2x,f8.3),5x,2(2x,f6.3))
!
 1010 format (//13x,' O u t p u t   o f   M a t e r i a l s'/ &
  / 10x,'G l a s s   M a t e r i a l s :' &
  / 5x,'No.   n_glass          k_glass ')
 1011 format (2x,i5,'.',2x,g12.5,2x,g14.6)
 1012 format (//5x,'Glass Properties' &
    / 7x,'theta',9x,'tau',12x,'rho',11x,'alpha')
 1013 format (5x,i5,'.',3(3x, f12.5))
! 
 3001 format (1x,'Error from read_materials, surf. mat. no. - ',i5)
 3003 format (1x,' -- specular reflectance out of range; [0,1], rho_s = ',g10.5)
 3004 format (1x,' -- semi-specular reflectance out of range; [0,1], rho_ss = ',g10.5)
 3005 format (1x,' -- diffuse reflectance out of range; [0,1], rho_d = ',g10.5)
 3006 format (1x,' -- semi-specular power out of range; [0,10], r_ss = ',g10.5)
 3007 format (1x,' -- diffuse power out of range; [0,10], r_d = ',g10.5)
 3010 format (1x,' -- material property sum out of range; [0,1], sum = ',g10.5)
 3011 format (1x,'Error from read_materials, glass . mat. no. - ',i5)
 3012 format (1x,' -- index of refraction out of range [0,5], n_g = ',g10.5)
 3013 format (1x,' -- extinction coefficient out of range [0,1], k_g = ',g10.5)
!
end subroutine read_materials
!
!*********************************************************************
!
subroutine read_surfaces
!
!*********************************************************************
!*     subroutine to read surfaces from input file                   *
!*********************************************************************
!
  use kinds
  use control
  use surfaces
!
  implicit none
!
  integer(kind=kindi) :: n, ierr, nn, mn
  integer(kind=kindi), allocatable, dimension (:) :: itest
  real(kind=kindr) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4
!
  write (*,*) 'Reading Surfaces'
  call header (9)
  write (9,1000)
!
  ierr = 0
!
!  allocate itest array
!
  allocate (itest(num_surfs))
  itest = 0
!
  do n = 1, num_surfs
    read (1, *, end = 3000, err = 4000) nn, mn
    read (1, *, end = 3000, err = 4000) x1, y1, z1
    read (1, *, end = 3000, err = 4000) x2, y2, z2
    read (1, *, end = 3000, err = 4000) x3, y3, z3
    read (1, *, end = 3000, err = 4000) x4, y4, z4
!
    if (nn < 1 .or. nn > num_surfs) then
      ierr = -1
      write (*,*) 'Error --  reading surfaces on line - ',n
      write (*,*) 'Surface no. -- ',nn,' out of range'
      write (9,*) 'Error --  reading surfaces on line - ',n
      write (9,*) 'Surface no. -- ',nn,' out of range'
    else
      matno_surf(nn) = mn
      x(1,nn) = x1
      x(2,nn) = x2
      x(3,nn) = x3
      x(4,nn) = x4
      y(1,nn) = y1
      y(2,nn) = y2
      y(3,nn) = y3
      y(4,nn) = y4
      z(1,nn) = z1
      z(2,nn) = z2
      z(3,nn) = z3
      z(4,nn) = z4
    end if
!
    if (itest(nn) /= 0) then
      ierr = 1
      write (*,*) 'Error --  reading surfaces on surface - ',n
      write (*,*) 'Surface no. -- ',nn,' already entered'
      write (9,*) 'Error --  reading surfaces on surface - ',n
      write (9,*) 'Surface no. -- ',nn,' already entered'
    else
      itest(nn) = 1
    end if
!
    if (matno_surf(nn) <=  0 .or. matno_surf(nn) > &
        num_mats_surf) then
      ierr = -1
      write (*,*) 'Error --  reading surfaces on surface - ',n
      write (*,*) 'Material no. for surface no. -- ',nn,' = ', matno_surf(nn), &
       ' is out of range'
      write (9,*) 'Error --  reading surfaces on surface - ',n
      write (9,*) 'Material no. for surface no. -- ',nn,' = ', matno_surf(nn), &
       ' is out of range'
    end if
  end do
!
!  test if all surfaces have been entered
!
  do n = 1, num_surfs
    if (itest(n) == 0) then
      write (*,*) 'Error -- in surfaces entered' 
      write (*,*) 'Surface no. -- ',n,' is missing'
      write (9,*) 'Error -- in surfaces entered' 
      write (9,*) 'Surface no. -- ',n,' is missing'
      ierr = -1
    end if
  end do
!
  do n = 1, num_surfs
    write (9, 1001) n, matno_surf(n)
    write (9, 1002)  x(1,nn), y(1,nn), z(1,nn), &
     x(2,nn), y(2,nn), z(2,nn), &
     x(3,nn), y(3,nn), z(3,nn), &
     x(4,nn), y(4,nn), z(4,nn)
  end do
!
  deallocate (itest)
!
  if (ierr == -1) then
    write (*,*) 'Errors detected in read_surfaces - stopping '
    write (9,*) 'Errors detected in read_surfaces - stopping '
    call adios (2)
  end if
!
  return
!
 3000 write (*,*) 'Error from read_directions --', &
  ' reached end of file at line - ',n
      write (9,*) 'Error from read_directions --', &
  ' reached end of file at line - ',n
  call adios (2)
  return
!
 4000 write (*,*) 'Error from read_directions -- ', &
   ' error reading input on line - ',n
      write (9,*) 'Error from read_directions -- ', &
   ' error reading input on line - ',n
  call adios (2)
  return
!
 1000 format (//14x,' O u t p u t   o f   S u r f a c e s'/ &
  / 5x,'  No.  ',10x,'Material No.  / 4 Corner Points (x, y, z)'/)
  1001 format (5x,i5,'.',10x,i5)
  1002 format (4( 1x, '(', f11.4, ',', 1x, f11.4, ',', 1x, f11.4, ')' / ))
!
end subroutine read_surfaces
!
!***********************************************************************
!***********************************************************************
!**                                                                   **
!**                  Solution Phase Routines Follow                   **
!**                                                                   **
!***********************************************************************
!***********************************************************************
!***********************************************************************
!
subroutine absref
!
!***********************************************************************
!*    subroutine to determine photon-material interactions             *
!*      output is iabs and new direction (not new X)                 *
!*      iabs =  2 for error                                            *
!*              1 for transmission                                     *
!*              0 for absorption                                       *
!*             -1 for diffuse reflection                               *
!*             -2 for specular reflection                              *
!*             -3 for semi-specular reflection                         *
!* Note: except for absorption, glass is treated as infinitely thin    *
!***********************************************************************
!*
!*  transformations
!*
!*   - forward
!*    {ex    }   | 1   0    0 | {ex}
!*    {er    } = | 0  ny   nz | {ey}
!*    {etheta}   | 0 -nz   ny | {ez}
!*
!*   - back
!*    {ex}   | 1   0    0 | {ex}
!*    {ey} = | 0  ny  -nz | {er}
!*    {ez}   | 0  nz   ny | {etheta}
!*
!*
  use kinds
  use constants
  use control
  use cylinders
  use materials
  use directions
  use photons
  use random
  use surfaces
!
  implicit none
!
  integer (kind=kindi) :: mn
  real(kind=kindr) ::  rhos, taus, frac
  real(kind=kindr) :: rn_in, rn1_out, rn2_out, ranf
  real(kind=kindr) :: exl, eyl, ezl, e_theta_i
!
  iabs = 2
  mn = matno_int(int_no)
  rn_in = ranf()
  if (idebug < 0) then
    write (8, 2001) np, int_no, mn, nx, ny, nz
  end if
!
!  surface material types
!
  if (int_no < int0_glass .or. int_no > int1_glass) then
!
!  Surface **************************************************
!
    if (rn_in > pdf_surf(3,mn)) then
!
!      - absorbed
!
      iabs = 0
!
    elseif (rn_in > pdf_surf(2, mn)) then
!
!   diffuse reflection
!
      iabs = -1
      rn1_out = ranf()
      theta_o = acos(rn1_out**rp1inv_d(mn))
      rn2_out = ranf()
      phi_o = p2i*rn2_out
      s_theta_o = sin(theta_o)
      exl = cos(phi_o)*s_theta_o
      eyl = sin(phi_o)*s_theta_o
      ezl = cos(theta_o)
      if (idebug < -1) then
        write (8, 3065) rn1_out, rn2_out, exl, eyl, ezl
      end if
      if (int_no == 2) then
        ex = exl
        ey = eyl
        ez = ezl
      elseif (int_no >= int0_cyl .and. int_no <= int1_cyl) then
        ex = exl
        ey =  nz*eyl + ny*ezl
        ez = -ny*eyl + nz*ezl
        if (idebug < -1) then
          write (8, 3066) ny, nz, ex, ey, ez
        end if
      else
!
!  error
!
        iabs = 2
      endif
    elseif (rn_in > pdf_surf(1, mn)) then
!
!   specular reflection
!
      iabs = -2
      if (int_no == 2) then
        ez = -ez
      elseif (int_no >= int0_cyl .and. int_no <= int1_cyl) then
        a11 = 1. - 2.*ny**2
        a12 = -2.*nz*ny
!
!   note: a21 = a12, a22 = -a11
!
        eysave = ey
        ey = a11*eysave + a12*ez
        ez = a12*eysave - a11*ez
      end if
    else
!
!      - semi-specular reflection
!
      iabs = -3
      c_theta_i = -(ex*nx+ey*ny+ez*nz)
      theta_i = acos(c_theta_i)
      e_theta_i = -nz*ey + ny*ez
!
! ERROR fix:  Jeff Miller, 11/14/97.
!
      IF (e_theta_i==0.0_kindr .and. ex==0.0_kindr) THEN
        phi_i = 0.0_kindr
      ELSE
        phi_i = atan2(e_theta_i, ex)
      END IF
      rn1_out = ranf()
      theta_o_90 = acos(rn1_out**rp1inv_ss(mn))
      delta_theta = (pi2-theta_i)*theta_o_90/pi2
      rn2_out = ranf()
      beta = p2i*rn2_out
      sb = sin(beta)
      cb = cos(beta)
!jeff      if (theta_i < dtor) then
!jeff        theta_o = delta_theta
!jeff        phi_o = beta
!jeff      else
      theta_o = theta_i + delta_theta*cb
      phi_o = phi_i + pi + delta_theta*sb
!jeff      end if
      if (theta_o < 0.) then
        theta_o = -theta_o
        phi_o = phi_o + pi
      end if
      s_theta_o = sin(theta_o)
      exl = cos(phi_o)*s_theta_o
      eyl = sin(phi_o)*s_theta_o
      ezl = cos(theta_o)
      if (int_no == 2) then
        ex = exl
        ey = eyl
        ez = ezl
      elseif (int_no >= int0_cyl .and. int_no <= int1_cyl) then
        ex = exl
        ey = (nz*eyl + ny*ezl)
        ez = (-ny*eyl + nz*ezl)
      else
        iabs = 2
      end if
    end if
  else
!
!  Glass ****************************************************
!    
    nc = mod(int_no, int0_glass) + 1
    c_theta_i = -(ny*ey + nz*ez)
    if (c_theta_i > 1.) then
      theta_i = 90.d0
      write (*, 6000) np, int_no, nc, ny, nz, ex, ey, ez, c_theta_i
      write (*, 6001) xe, ye, ze, xint, yint, zint
      write (9, 6000) np, int_no, nc, ny, nz, ex, ey, ez, c_theta_i
      write (9, 6001) xe, ye, ze, xint, yint, zint
      if (idebug /= 0) then
        write (8, 6000) np, int_no, nc, ny, nz, ex, ey, ez, c_theta_i
        write (8, 6001) xe, ye, ze, xint, yint, zint
      end if
    else
      theta_i = acos(c_theta_i)*rtod
    end if
    itheta_i = theta_i
    frac = theta_i - itheta_i   
!
!   - test if reflected
!
    rhos = rho_glass(itheta_i,mn) + frac*drho_glass(itheta_i,mn)
    taus = tau_glass(itheta_i,mn) + frac*dtau_glass(itheta_i,mn)
    if (idebug < 0) then
      write (8, 2002) nc, itheta_i, frac, taus, rhos
    end if
    if (rn_in < taus) then
      iabs = 1
    elseif (rn_in < taus + rhos) then
!
!   specular reflection matrix: a12 = a21, a22 = - a11
!
      iabs = -2
      a11 = 1. - 2.*ny**2
      a12 = -2.*ny*nz
      eysave = ey
      ey = a11*eysave + a12*ez
      ez = a12*eysave - a11*ez
    else
!
!  Absorbed
!
      iabs = 0
    end if
  end if
!
  if (idebug < 0) then
    write (8,1000) int_no, iabs, rn_in, ex, ey, ez
  end if
!
  return
!
 1000 format (' absref: int_no, iabs, rn_in, E_o: ',2i3,1x,f6.4,3x,3(1x,f8.4))
 2001 format(' absref: np, int_no, mn, nxyz = ',i5,2i3,3(1x,f7.4))
 2002 format(' -  glass: nc, itheta, frac, taus, rhos = ',2i3,3(1x,f7.4))
 3065 format(3x,'rn1&2_out, exl, eyl, ezl - ',5(1x,e12.5))
 3066 format(3x,'ny, nz, exyz - ',5(1x,e12.5))
 6000 format(2x,'Error in absref - theta > 90 - setting to 90:' &
  / 5x,'np, int_no, nc, ny, nz, exyz, cth - ',i10, 2i5, &
  / 5x, 6(1x,e12.5))
 6001 format(2x,'Xe, Xint(3) = ',6(1x,f10.5))
!
end subroutine absref
!
!***********************************************************************
!
function ranf ()
!
!***********************************************************************
!*    function to generate pseudo-random numbers                       *
!***********************************************************************
!
  use random
!
  implicit none
!
  real(kind=kindr) :: ranf
!
!     get new seed
!
  seed(it1) = seed(it1) + seed(it2)
  seed(it1) = iand(seed(it1), modulus)
!
!     generate normalized random number
!
  ranf = seed(it1)*odenom
!
!     cyclic permutation of taps
!
  it1 = it1 - 1
  it2 = it2 - 1
  if (it1 < 1) then
     it1 = 17
  elseif (it2 < 1) then
     it2 = 17
  endif
!
  return
!
end function ranf
!
!*********************************************************************
!
subroutine geteplane
!
!*********************************************************************
!*     subroutine to compute geometry of emission plane              *
!*********************************************************************
!
  use control
  use directions
  use planes
  use cylinders
!
  implicit none
!
  integer(kind=kindi) :: n, m
  real(kind=kindr) :: dxb, dyb, b, c, r, s, t, ix, iy, px, py, dx, dy
  real(kind=kindr) :: dist1, dist2, slope1, A, E, mx, my, zbase
!
  call header (9)
  write (9,1000) nd, phid(nd), altituded(nd), ztop
  if (idebug /= 0) then
    call header (8)
    write (8,1000) nd, phid(nd), altituded(nd), ztop
  end if
!
! compute size of backplane
!
    dxb = xmaxb - xminb
    dyb = ymaxb - yminb
    areab = dxb*dyb
!
! altitude = 90 degrees
!
  if (altituded(nd) == 90.) then
    xep(3,1:4) = ztop
    xep(1,1:4:3) = xmaxb
    xep(1,2:3)   = xminb
    xep(2,1:2)   = yminb
    xep(2,3:4)   = ymaxb
!
! altitude /= 90
!
!
! Points 1 and 4
!
!***********************************************************************
!
! New method for computing points 1 & 4 on the emission plane
! (px,py) is the point in the x-y plane the emission plane pivots around
! --->Changed by Joe Ryan (6/27/96)
!
!***********************************************************************
!
  else
    b = sp(nd)*zmax*ta(nd)     
    c = cp(nd)*zmax*ta(nd)
    if ((xmax + c) > xmaxb) then
      px  = xmax + c
      dxb = px - xminb 
    else 
      px = xmaxb
    end if
    if (phid(nd) <= 0) then
      if ((ymax - b) > ymaxb) then
        py = ymax - b
        dyb = py - yminb
      else
        py = ymaxb
      end if 
    else
      if ((ymin - b) < yminb) then
        py = ymin - b
        dyb = ymaxb - py 
      else
        py = yminb
      end if
    end if
!
! compute (xmine, ymine), lower left corner - Point 1
!
    if (phid(nd) <= 0) then
      xep(1,1) = px - dyb*cp(nd)*sp(nd)
      xep(2,1) = py - dyb*cp(nd)**2
    else
      xep(1,1) = px - dxb*sp(nd)**2
      xep(2,1) = py - dxb*sp(nd)*cp(nd)
    end if  
    xep(3,1) = 0.
!
! compute (xmaxe, ymaxe), lower right corner - Point 4
!
    if (phid(nd) <= 0) then
      xep(1,4) = px - dxb*sp(nd)**2
      xep(2,4) = py - dxb*sp(nd)*cp(nd)
    else
      xep(1,4) = px + dyb*sp(nd)*cp(nd)
      xep(2,4) = py + dyb*(1 - sp(nd)**2)
    end if
    xep(3,4) = 0.
!
! Points 2 and 3
!
!********************************************************************
!
! New method for computing points 2 & 3 on the emission plane
! (ix,iy) is the point in the x-y plane which must be included in the 
! projected emission plane 
! --->Changed by Joe Ryan (6/28/96)
!
!********************************************************************
!
    r = zmax*cp(nd)/ta(nd)
    s = zmax*sp(nd)/ta(nd)
!  
    ix = xmin - r
    if (ix > xminb) then
      ix = xminb
    end if 
    if (phid(nd) <= 0) then
      iy = ymin + s
      if (iy > yminb) then
        iy = yminb
      end if
    else 
      iy = ymax + s
      if (iy < ymaxb) then
        iy = ymaxb
      end if
    end if
!
    if (phid(nd) == 0) then
      dist1 = xep(1,1) - ix
      dist2 = xep(1,1) - xmaxb
    else
      slope1 = (xep(2,1) - xep(2,4))/(xep(1,1) - xep(1,4))
      A = -slope1
      E = -xep(2,1) + slope1*xep(1,1)
      dist1 = abs((A*ix + iy + E)/sqrt(A**2 + 1))
      if (phid(nd) <= 0) then
        dist2 = abs((A*xmaxb + ymaxb + E)/sqrt(A**2 + 1)) 
      else 
        dist2 = abs((A*xmaxb + yminb + E)/sqrt(A**2 + 1))
      end if
    end if
!
    ztop = dist1*ca(nd)*sa(nd)
    t = ztop*ta(nd)
    dx = t*cp(nd)
    dy = t*sp(nd)
!
!
! compute (xmine, ymine), upper left corner - Point 2
!
    xep(1,2) = xep(1,1) - dx 
    xep(2,2) = xep(2,1) + dy
    xep(3,2) = ztop
!
! compute (xmaxe, ymaxe), uppper right corner - Point 3
!
    xep(1,3) = xep(1,4) - dx
    xep(2,3) = xep(2,4) + dy
    xep(3,3) = ztop
!
! Determine if emission plane is too big
!
    if (dist2 > 0) then
      zbase = sa(nd)*ca(nd)*dist2
      mx = dist2*(sa(nd)**2)*cp(nd)
      my = dist2*(sa(nd)**2)*sp(nd)
!
! New Point 1
!
      xep(1,1) = xep(1,1) - mx
      xep(2,1) = xep(2,1) + my
      xep(3,1) = zbase
!
! New Point 4
!
      xep(1,4) = xep(1,4) - mx
      xep(2,4) = xep(2,4) + my
      xep(3,4) = zbase
!
    end if
  end if
!
  zmaxe = xep(3,3)
!
!  area
!
  areae = sqrt((xep(1,1)-xep(1,4))**2+(xep(2,1)-xep(2,4))**2) * &
    sqrt((xep(1,2)-xep(1,3))**2+(xep(2,2)-xep(2,3))**2)+(xep(3,2)-xep(3,3)**2)
!
!  Print Out Points
!
  do n = 1, 4
    write (9,1001) n, (xep(m,n), m = 1,3)
    if (idebug /= 0) then
      write (8,1001) n, (xep(m,n), m = 1,3)
    end if
  end do
!
  return
!
 1000 format (5x,'E m i s s i o n   P l a n e   F o r   D i r e c t i o n' &
   ,'   -', i5,/ &
 / 10x, 'azimuthal angle, phi (degrees)       = ',f10.2 &
 / 10x, 'altitude angle, altitude (degrees)   = ',f10.2 &
 / 10x, 'z_top                                = ',f10.2 &
 // 22x,'x',12x,'y',12x,'z')
 1001 format (1x,'Point -', i5, '.', 4(1x,f12.5))
!
end subroutine geteplane
!
!***********************************************************************
!
subroutine intsec
!
!***********************************************************************
!*    subroutine to perform intersection calculations                  *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use random
  use constants
  use cylinders
  use directions
  use materials
  use photons
  use planes
!
  implicit none
!
  integer(kind=kindi) :: reflect
  real(kind=kindr) :: ranf
  real(kind=kindr) :: rn1, rn2, rn1c, rn2c, rn12, rn1c2, rn12c, rn1c2c
  real(kind=kindr) :: xinta, yinta, slopex, slopey
!
!  Note: except for absorption, glass is treated as infinitely thin
!
!  Initial Debug Info.
!
  if (idebug < 0) then
    write (8, 3000) nd, nploop, ex_plane, ey_plane, ez_plane
    write (8,*) 'num_photons = ',num_photons
  end if
!
! Initialization
!
  xint = 0.
  yint = 0.
  zint = -1.
  do np = 1, num_photons
!
!    emission point
!
 50 continue
    rn1 = ranf()
    rn2 = ranf()
    rn1c = 1. - rn1
    rn2c = 1. - rn2
    rn12 = rn1*rn2
    rn1c2 = rn1c*rn2
    rn12c = rn1*rn2c
    rn1c2c = rn1c*rn2c
    xe = xep(1,1)*rn1c2c+xep(1,2)*rn1c2+xep(1,3)*rn12+xep(1,4)*rn12c
    ye = xep(2,1)*rn1c2c+xep(2,2)*rn1c2+xep(2,3)*rn12+xep(2,4)*rn12c
    ze = xep(3,1)*rn1c2c+xep(3,2)*rn1c2+xep(3,3)*rn12+xep(3,4)*rn12c
!
!  Emission Plane Directions
!
    ex = ex_plane
    ey = ey_plane
    ez = ez_plane
!
    slopex = -ex/ez
    slopey = -ey/ez
!
!  Intersection w/ aperture plane
!
    xinta = slopex*ze + xe
    yinta = slopey*ze + ye
    if ( (xinta-xmina)*(xinta-xmaxa) <= 0. .and. &
         (yinta-ymina)*(yinta-ymaxa) <= 0. ) then
        npabs_a(nd) = npabs_a(nd) + 1
    end if
!
!  Here for reflection
!
    nc_save = 0
    reflect = -1
 100 reflect = reflect + 1
    ey2 = ey**2
    ez2 = ez**2
    a = ey2 + ez2
    if (idebug < 0) then
      write (8,3002) np, reflect, xe, ye, ze, ex, ey, ez
    end if
!
!  Intersection with Cylinders
!
!  compute possible intersection with all cylinders
!   - int_no 
!         =-1 for lost
!         = 0 for top plane
!         = 1 for bottom plane
!         = 2 for backplane
!         = int0_cyl to int0_cyl+num_cyls for cylinders
!         = int0_glass to int0_glass+num_cyls for glass
!         = int0_surf to int0_surf+num_surfs for surfaces
!
    int_no = -1
    dzint_old = 1.e30
!
!  Intersection with Cylinders
!
    call intsec_cyls
    if (nc_int /= 0) then
      nc_save = nc_int
      int_no = nc_int + int0_glass - 1
      if (idebug < 0) then
        write (8, 3018) int_no, xint, yint, zint, dzint_old
      end if
    end if
!
!  Intersection with Surfaces
!
    if (num_surfs > 0) then
      call intsec_surfs
      if (idebug < 0) then
        write (8, 3019) int_no, xint, yint, zint, dzint_old
      end if
    end if
!
!  Intersection w/ Backplane
!
!    int_no = 0: top plane
!        = 1: bottom plane
!        = 2: back plane
!
    if (int_no < 0) then
!
!     - intersection with bottom plane
!
      if (ez < 0.) then
!       write (8,*) 'Xe, E:',xe,ye,ze,ex,ey,ez
        xint = -ex*ze/ez + xe
        yint = -ey*ze/ez + ye
        zint = 0.
        if ( (xint-xminb)*(xint-xmaxb) <= 0. .and. &
                (yint-yminb)*(yint-ymaxb) <= 0. ) then
          int_no = 2
		  nc_save=0    !this line corrects an error.  Jeff Miller 12/97
        else
          int_no = 1
          npabs(int_no, nd) = npabs(int_no, nd) + 1
          if (idebug < 0) then
            write (8, 3016) int_no
          end if
          cycle
        end if
        if (idebug < 0) then
          write (8, 3012) int_no
        end if
!
!      - intersection with top plane
!
      else
        if (idebug < 0) then
          write (8, 3013) int_no
        end if
      int_no = 0
      npabs(0, nd) = npabs(0, nd) + 1
      cycle
      end if
    end if
!
! Final Accumulation and Test
!
    if (idebug < 0) then
      write (8, 3025) int_no, xint, yint, zint, dzint_old
    end if
!
!  Photon Interaction
!    iabs =  1 for transmission
!            0 for absorption
!           -1 for diffuse reflection
!           -2 for specular reflection
!           -3 for semi-specular reflection
!
    if (int_no < 0) then
!
!  Lost
!
      nplost(nd) = nplost(nd) + 1
      write (*,*) ' - lost photon no. - ',nplost(nd),'" direction - ',nd
      write (4,3040) nd, nplost(nd), xe, ye, ze, xe+ex, ye+ey, ze+ez
      if (nplost(nd) > nplost_max) then
        write (*, 4000) nd, nploop, np, nplost(nd)
        write (9, 4000) nd, nploop, np, nplost(nd)
        call adios (2)
      end if
      go to 50
!
!  Back Plane
!
    elseif (int_no == 2) then
      nx = 0.
      ny = 0.
      nz = 1.d0
      call absref
      if (iabs == 0) then 
        npabs(int_no,nd) = npabs(int_no,nd) +1
      elseif (iabs < 0) then
        xe = xint
        ye = yint
        ze = 0
        if (idebug < -1) then
          write (8,3031) int_no
        end if
        go to 100
      end if
!
!  Outside cylinders - glass
!
    elseif (int_no >= int0_glass .and. int_no <= int1_glass) then
      if (endcap == 1) then
!
! Absorbed by cylinder endcap (cylinder)
!
        iabs = 0
        npabs_caps(nd) = npabs_caps(nd) + 1
      else
        nx = 0.
        ny = (yint - y0(nc_int))/ro(nc_int)
        nz = (zint - z0(nc_int))/ro(nc_int)
        call absref
        if (iabs == 0) then
!
!  Absorbed
!
          npabs (int_no, nd) = npabs (int_no, nd) + 1
          if (idebug < -1) then
            write (8, 3030) int_no
          end if
        elseif (iabs < 0) then
!
!  Reflected
!
          xe = xint
          ye = yint
          ze = zint
          if (idebug < -1) then
            write (8, 3031) int_no
          end if
          go to 100
        else
!
!  Transmitted
!
          xe = xint
          ye = yint
          ze = zint
          if (idebug < -1) then
            write (8, 3032) int_no
          end if
          call intsec_inside
!
!  Error fix.  Jeff Miller, 12/97.
!
		  if (endcap==1) then
		    npabs_caps(nd) = npabs_caps(nd) + 1
		    cycle
	      end if
!
!  Now, either absorbed, or exited from cylinder (outside)
!
          if (iabs == 0) then
            npabs(int_no, nd) = npabs(int_no, nd) + 1
          else
            xe = xint
            ye = yint
            ze = zint
            go to 100
          end if
        end if
      end if
!
! Ordinary Surface
!
!   elseif (int_no >= int0_surf .and. int_no <= int1_surf) then
!     ns = mod(int_no, int0_surf) + 1
!     nx = cgx(ns)
!     ny = cgy(ns)
!     nz = cgz(ns)
!     call absref
    else
      write (9, 3020) nd, np, nploop, int_no
      write (*, 3020) nd, np, nploop, int_no
      go to 50
    end if
  end do
!
  return
!
 3000 format (1x,'*** intsec - ', &
     'nd,nploop= ', 2i3, 'ex,ey,ez=', 3(1x,e11.4))
 3002 format (1x,'np,reflect- ',2i5,2x,'Xe= ', 3(1x,e11.4),', E= ',3(1x,e11.4))
 3012 format (1x,'int_no ',i5,' - int_no. back')
 3013 format (1x,'int_no ',i5,' - int_no. top')
 3014 format (5x,'np,      int_no,       xint,       yint,         zint')
 3015 format (3x,i5,3x,i5,3(1x,e11.4))
 3016 format (1x,'int_no ',i5,' - int_no. bottom')
 3018 format (1x,'cyls:  int_no, xyz-int,dzint_old,=',i5, 4(1x,g12.5))
 3019 format (1x,'surfs: int_no,xyz-int,dzint_old,=', i5, 4(1x,g12.5))
 3020 format (1x,'intsec error - int out of range: nd, np, nploop, int_no= ',4i5)
 3025 format (1x,'final: int_no,xyz-int,dzint_old,=', i5, 4(1x,g12.5))
 3040 format (1x,'intsec lost: nd, #lost,Xe, Xe+E= ',2i5,6(1x,g12.5))
 3030 format (1x,'intsec (absorbed): int_no =', 2i5)
 3031 format (1x,'intsec (reflected): int_no =', 2i5)
 3032 format (1x,'intsec (transmitted): int_no =', 2i5)
 4000 format (1x,'Error - lost photons exceeds specified limit', &
  / 5x,'at direction no. ',i5,'  - nploop no. ',i5,'  - np ',i5, &
    '  - no. lost this direction = ',i10)
!
end subroutine intsec
!
!***********************************************************************
!
subroutine intsec_cyls
!
!***********************************************************************
!*    subroutine to perform intersection calculations for outside cyls *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use random
  use constants
  use materials
  use cylinders
  use directions
  use photons
!
  implicit none
!
  real(kind=kindr) :: zz
!
  nc_int = 0
  endcap = 0
  do nc = 1, num_cyls
    if (nc /= nc_save) then
!
!  skip cylinders "behind" emission point
!
      if (ey >= 0.) then
        if ((ye - ymax_cyl(nc)) > 0.) then
           cycle
        end if
      else
        if ((ymin_cyl(nc) - ye) > 0.) then
           cycle
        end if
      end if
!
!  now for cylinders in front of emission point
!
      yehat = ye - y0(nc)
      zehat = ze - z0(nc)
      b1 = ez*yehat - ey*zehat
      b = 2.d0*ey*b1
      c = b1**2 - ez2*ro(nc)**2
!
      rad = b**2 - 4.d0*a*c
      if (idebug < -1) then
        write (8, 3010) nc, rad, a, b, c
      end if
      if (rad > 0.) then
        zz = -b/a
        zi(1) = 0.5d0*(zz + sqrt(rad)/a) + z0(nc)
        xi(1) = ex*(zi(1) - ze)/ez + xe
        dzint(1) = (zi(1) - ze)*ez
        zi(2) = 0.5d0*(zz - sqrt(rad)/a) + z0(nc)
        xi(2) = ex*(zi(2) - ze)/ez + xe
        dzint(2) = (zi(2) - ze)*ez
!
        if (idebug < -1) then
          write (8, 3080) xi, yi, zi, dzint
        end if
!
        if (dzint(1) > 0. .and. dzint(2) > 0.) then
          if (dzint(1) < dzint(2)) then
!
! See if photon intersects end of cylinder first
!
            if (xi(1) > xmax_cyl(nc) .and. xi(2) < xmax_cyl(nc)) then
              endcap = 1
              nc_int = nc
              xint = xmax_cyl(nc)
              yint = ey*(xint - xe)/ex + ye
              zint = ez*(xint - xe)/ex + ze
              dzint_old = (zint - ze)*ez
              cycle
!
            elseif (xi(1) < xmin_cyl(nc) .and. xi(2) > xmin_cyl(nc)) then
              endcap = 1
              nc_int = nc
              xint = xmin_cyl(nc)
              yint = ey*(xint - xe)/ex + ye
              zint = ez*(xint - xe)/ex + ze
              dzint_old = (zint - ze)*ez
              cycle
!
            elseif (xi(1) < xmin_cyl(nc) .or. xi(1) > xmax_cyl(nc)) then
              cycle
            end if
!
            if (dzint(1) < dzint_old) then
              nc_int = nc
              dzint_old = dzint(1)
              zint = zi(1)
              xint = xi(1)
              yint = ey*(zint - ze)/ez + ye
            end if
          else
!
! See if photon intersects end of cylinder first
!
            if (xi(2) > xmax_cyl(nc) .and. xi(1) < xmax_cyl(nc)) then
              endcap = 1
              nc_int = nc
              xint = xmax_cyl(nc)
              yint = ey*(xint - xe)/ex + ye
              zint = ez*(xint - xe)/ex + ze
              dzint_old = (zint - ze)*ez
              cycle
!
            elseif (xi(2) < xmin_cyl(nc) .and. xi(1) > xmin_cyl(nc)) then
              endcap = 1
              nc_int = nc
              xint = xmin_cyl(nc)
              yint = ey*(xint - xe)/ex + ye
              zint = ez*(xint - xe)/ex + ze
              dzint_old = (zint - ze)*ez
              cycle    
!
            elseif (xi(2) < xmin_cyl(nc) .or. xi(2) > xmax_cyl(nc)) then
              cycle
            end if
!
            if (dzint(2) < dzint_old) then
              nc_int = nc
              dzint_old = dzint(2)
              zint = zi(2)
              xint = xi(2)
              yint = ey*(zint - ze)/ez + ye
            end if
          end if
!
        elseif (dzint(1) > 0. .and. dzint(2) < 0.) then
!
! See if photon intersects end of cylinder first
!
          if (xi(1) < xmax_cyl(nc) .and. xi(2) > xmax_cyl(nc)) then
            endcap = 1
            nc_int = nc
            xint = xmax_cyl(nc)
            yint = ey*(xint - xe)/ex + ye
            zint = ez*(xint - xe)/ex + ze
            dzint_old = (zint - ze)*ez
            cycle
!
          elseif (xi(1) > xmin_cyl(nc) .and. xi(2) < xmin_cyl(nc)) then
            endcap = 1
            nc_int = nc
            xint = xmin_cyl(nc)
            yint = ey*(xint - xe)/ex + ye
            zint = ez*(xint - xe)/ex + ze
            dzint_old = (zint - ze)*ez
            cycle
!
          elseif (xi(1) < xmin_cyl(nc) .or. xi(1) > xmax_cyl(nc)) then
            cycle
          end if
!
          if (dzint(1) < dzint_old) then
            nc_int = nc
            dzint_old = dzint(1)
            zint = zi(1)
            xint = xi(1)
            yint = ey*(zint - ze)/ez + ye
          end if
!
        elseif (dzint(1) < 0. .and. dzint(2) > 0.) then
!
! See if photon intersects end of cylinder first
!
          if (xi(2) < xmax_cyl(nc) .and. xi(1) > xmax_cyl(nc)) then
            endcap = 1
            nc_int = nc
            xint = xmax_cyl(nc)
            yint = ey*(xint - xe)/ex + ye
            zint = ez*(xint - xe)/ex + ze
            dzint_old = (zint - ze)*ez
            cycle
!
          elseif (xi(2) > xmin_cyl(nc) .and. xi(1) < xmin_cyl(nc)) then
            endcap = 1
            nc_int = nc
            xint = xmax_cyl(nc)
            yint = ey*(xint - xe)/ex + ye
            zint = ez*(xint - xe)/ex + ze
            dzint_old = (zint - ze)*ez
            cycle
!
          elseif (xi(2) < xmin_cyl(nc) .or. xi(2) > xmax_cyl(nc)) then
            cycle
          end if
!
          if (dzint(2) < dzint_old) then
            nc_int = nc
            dzint_old = dzint(2)
            zint = zi(2)
            xint = xi(2)
            yint = ey*(zint - ze)/ez + ye
          end if
        end if
      end if
    end if
  end do
!
  return
!
 3010 format (1x, 'intsec_cyls: nc, rad, a, b, c = ',i5, 4(2x, g11.4))
 3080 format (2x,'xyzi(2),dzint(2) = ',8(1x,f7.2))
!
end subroutine intsec_cyls
!
!***********************************************************************
!
subroutine intsec_inside
!
!***********************************************************************
!*    subroutine to perform intersection calculations inside glass     *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use random
  use constants
  use cylinders
  use materials
  use directions
  use photons
  use planes
!
  implicit none
!
  integer(kind=kindi) :: reflect
  real(kind=kindr) :: zz
!
  nx = 0.
  int_no = -1
  reflect = -1
!
  if (idebug < 0) then
    write (8, 3040) nc_int, xe, ye, ze, ex, ey, ez
  end if
!
 100 reflect = reflect + 1
  ey2 = ey**2
  ez2 = ez**2
  a = ey2 + ez2
  yehat = ye - y0(nc_int)
  zehat = ze - z0(nc_int)
  b1 = ez*yehat - ey*zehat
  b = 2.d0*ey*b1
!
!  intersection with inside cylinder
!
  c = b1**2 - ez2*ri(nc_int)**2
!
  rad = b**2 - 4.d0*a*c
  if (idebug < -1) then
    write (8, 3010) nc_int, rad, a, b, c
  end if
  if (rad > 0.) then
    zz = -b/a
    zi(1) = 0.5d0*(zz + sqrt(rad)/a) + z0(nc_int)
    dzint(1) = (zi(1) - ze)*ez
    zi(2) = 0.5d0*(zz - sqrt(rad)/a) + z0(nc_int)
    dzint(2) = (zi(2) - ze)*ez
!
    if (idebug < -1) then
      write (8, 3080) zi, dzint
    end if
!
!  - now, both dzints must be positive
!     Note: if inside, may or may not hit inside cyl
!      - don't check for being within x extents
!
    if (dzint(1) < dzint(2)) then
      zint = zi(1)
      xint = ex*(zi(1) - ze)/ez + xe
      yint = ey*(zi(1) - ze)/ez + ye
    else
      zint = zi(2)
      xint = ex*(zi(2) - ze)/ez + xe
      yint = ey*(zi(2) - ze)/ez + ye
    end if
!
! First see if photon intersects end of cylinder.  Error correction, Jeff Miller 12/97.
!
    if (xint > xmax_cyl(nc_int)) then
      endcap = 1
      int_no = nc_int + int0_glass - 1
      xint = xmax_cyl(nc_int)
      yint = ey*(xint - xe)/ex + ye
      zint = ez*(xint - xe)/ex + ze
      iabs = 0
      goto 200
    end if
!
    if (xint < xmin_cyl(nc_int)) then
      endcap = 1
      int_no = nc_int + int0_glass - 1
      xint = xmax_cyl(nc_int)
      yint = ey*(xint - xe)/ex + ye
      zint = ez*(xint - xe)/ex + ze
	  iabs = 0
      goto 200
    end if
!
!  End of correction.
!
    int_no = nc_int + int0_cyl - 1
    ny = (yint -y0(nc_int))/ri(nc_int)
    nz = (zint -z0(nc_int))/ri(nc_int)
    call absref
    if (iabs == 0) then
!
!    absorbed
!
      go to 200
    else
!
!    reflected from inside cylinder - no transmission
!
      xe = xint
      ye = yint
      ze = zint
!
!  get exit point from cylinder
!
      ey2 = ey**2
      ez2 = ez**2
      a = ey2 + ez2
      yehat = ye - y0(nc_int)
      zehat = ze - z0(nc_int)
      b1 = ez*yehat - ey*zehat
      b = 2.d0*ey*b1
      c = b1**2 - ez2*ro(nc_int)**2
!
      rad = b**2 - 4.d0*a*c
      if (idebug < -1) then
        write (8, 3011) nc_int, rad, a, b, b1, c
        write (8, 3050) xe, ye, ze, ex, ey, ez
        write (8, 3051) nx, ny, nz
        write (8, 3052) y0(nc_int), z0(nc_int), yehat, zehat
      end if
!
      if (rad < 0.) then
        write (9, 4011) np, reflect, nc_int, rad, a, b, c
        write (*, 4011) np, reflect, nc_int, rad, a, b, c
        iabs = 1
        return
      end if
!
      zz = -b/a
      zi(1) = 0.5d0*(zz + sqrt(rad)/a) + z0(nc_int)
      dzint(1) = (zi(1) - ze)*ez
      zi(2) = 0.5d0*(zz - sqrt(rad)/a) + z0(nc_int)
      dzint(2) = (zi(2) - ze)*ez
!
      if (idebug < -1) then
        write (8, 3081) zi, dzint
      end if
!
!  - now, only one dzint is positive, the other is zero
!     Note: if inside, end check taken out
!
      if (dzint(1) > dzint(2)) then
        zint = zi(1)
        xint = ex*(zint - ze)/ez + xe
        yint = ey*(zint - ze)/ez + ye
      else
        zint = zi(2)
        xint = ex*(zint - ze)/ez + xe
        yint = ey*(zint - ze)/ez + ye
      end if
!
! First see if photon intersects end of cylinder.  Error correction, Jeff Miller 12/97.
!
      if (xint > xmax_cyl(nc_int)) then
        endcap = 1
        int_no = nc_int + int0_glass - 1
        xint = xmax_cyl(nc_int)
        yint = ey*(xint - xe)/ex + ye
        zint = ez*(xint - xe)/ex + ze
    	iabs = 0
        goto 200
      end if
!
      if (xint < xmin_cyl(nc_int)) then
        endcap = 1
        int_no = nc_int + int0_glass - 1
        xint = xmax_cyl(nc_int)
        yint = ey*(xint - xe)/ex + ye
        zint = ez*(xint - xe)/ex + ze
	    iabs = 0
        goto 200
      end if
!
!  End of correction.
!
    end if
  else
!
!  missed inside cylinder
!    - intersect with outside cylinder
!
    c = b1**2 - ez2*ro(nc_int)**2
    rad = b**2 - 4.d0*a*c
    if ( rad < 0.) then
      write (9, 4050) np, nc_int, xe, ye, ze, ex, ey, ez
      write (9, 4051) rad, a, b, b1, c
      write (*, 4050) np, nc_int, xe, ye, ze, ex, ey, ez
      write (*, 4051) rad, a, b, b1, c
    else
      zz = -b/a
      zi(1) = 0.5d0*(zz + sqrt(rad)/a) + z0(nc_int)
      dzint(1) = (zi(1) - ze)*ez
      zi(2) = 0.5d0*(zz - sqrt(rad)/a) + z0(nc_int)
      dzint(2) = (zi(2) - ze)*ez
      if (idebug < -1) then
        write (8, 3082) zi, dzint
      end if
      if (dzint(1) > dzint(2)) then
        zint = zi(1)
        xint = ex*(zint - ze)/ez + xe
        yint = ey*(zint - ze)/ez + ye
      else
        zint = zi(2)
        xint = ex*(zint - ze)/ez + xe
        yint = ey*(zint - ze)/ez + ye
      end if
      if (idebug < -1) then
       write (8, 4060) xint, yint, zint
      end if
    end if
  end if
!
! compute photon interaction on outside cylinder
!
! First see if photon intersects end of cylinder.  Error correction, Jeff Miller 12/97.
!
  if (xint > xmax_cyl(nc_int)) then
    endcap = 1
    int_no = nc_int + int0_glass - 1
    xint = xmax_cyl(nc_int)
    yint = ey*(xint - xe)/ex + ye
    zint = ez*(xint - xe)/ex + ze
	iabs = 0
    goto 200
  end if
!
  if (xint < xmin_cyl(nc_int)) then
    endcap = 1
    int_no = nc_int + int0_glass - 1
    xint = xmax_cyl(nc_int)
    yint = ey*(xint - xe)/ex + ye
    zint = ez*(xint - xe)/ex + ze
	iabs = 0
    goto 200
  end if
!
!  End of correction.
!
  int_no = nc_int + int0_glass - 1
  ny = -(yint-y0(nc_int))/ro(nc_int)
  nz = -(zint-z0(nc_int))/ro(nc_int)
  if (idebug < -1) then
    write (8, 3060) nc_int, xint, yint, zint, ex, ey, ez
    write (8, 3061) nx, ny, nz
  end if
  call absref
  if (iabs >= 0) then
!
!  0 - absorbed, 1 - transmitted: return
!
    go to 200
  else
!
!  reflected internally
!
    xe = xint
    ye = yint
    ze = zint
    if (idebug < -1) then
      write (8, 3070) xe, ye, ze
    end if
    go to 100
  end if
!
 200 continue
  if (idebug < -1) then
    write (9,3020) int_no, iabs
  end if
!
  return
!
 3010 format (1x, 'going  in: nc_int, rad, a, b, c = ',i5, 4(2x, g11.4))
 3011 format (1x, 'going out: nc_int,rad,a,b,b1,c = ',i5, 5(2x, f10.4))
 3020 format (1x,' exiting intsec_inside, int_no, iabs = ', 2i5)
 3030 format (1x,' inside: int_no, xyzint = ', i5,3(1x,f12.5))
 3040 format (1x,' intsec_inside: nc_int, X, E = ', i5,6(1x,f12.5))
 3050 format(' pre-int w/ out cyl: Xe, E = ',6(1x,f7.3))
 3051 format('  N_cyl = ',6(1x,f7.3))
 3052 format('  y0,z0,yehat,zehat = ',4(1x,g12.5))
 3060 format(' pre-absref for out cyl: nc, Xint, E = ',i3,6(1x,f7.3))
 3061 format('  N_cyl = ',6(1x,f7.3))
 3070 format (2x,' reflected internally: Xe = ',3(1x,f10.5))
 3080 format (2x,' inside: zi1&2, dzint(1)&2 = ',4(1x,f7.2))
 3081 format (2x,'outside: zi1&2, dzint(1)&2 = ',4(1x,f7.2))
 3082 format (2x,'Missed inside - outside: zi1&2, dzint(1)&2 = ',4(1x,f7.2))
 4011 format (1x, 'intsec_inside error - rad < 0 for outside cylinder:' &
  / 5x,' np,reflect,nc_int,rad,a,b,c = ',i7, 2i5, 4(2x, g11.4))
 4050 format (1x, 'Error from intsec_inside - photon no. - ',i5 &
  / 5x,'inside cylinder no - ',i5,'  no intersection found with' &
  / 5x,'outside cylinder; Xe, E = ',3(1x,f11.5) / 5x,3(1x,f10.4))
 4051 format (2x,'a, b, b1, c = ',4(1x,e12.5))
 4060 format (1x, 'missed inside cylinder, Xi = ',3(1x,f11.5))
!
end subroutine intsec_inside
!
!***********************************************************************
!
subroutine intsec_ng
!
!***********************************************************************
!*    subroutine to perform intersection calculations for no glass     *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use random
  use constants
  use cylinders
  use materials
  use directions
  use photons
  use planes
!
  implicit none
!
  integer(kind=kindi) :: reflect
  real(kind=kindr) :: ranf
  real(kind=kindr) :: rn1, rn2, rn1c, rn2c, rn12, rn1c2, rn12c, rn1c2c
  real(kind=kindr) :: xinta, yinta, slopex, slopey
!
!  Initial Debug Info.
!
  if (idebug < 0) then
    call header (8)
    write (8, 3000) nd, nploop, ex_plane, ey_plane, ez_plane
  end if
!
!  Initialization
!
!
!  Error correction.  Jeff Miller, 12/97
!  These lines are not in the correct position in this subroutine, but
!    are in the subroutine intsec.  See comments below.
!
!  slopex = -ex/ez
!  slopey = -ey/ez
  xint = 0.
  yint = 0.
  zint = -1.
!
!  Photon Emission Loop.
!
  do np = 1, num_photons
!
!    Emission point
!
 50 continue
    rn1 = ranf()
    rn2 = ranf()
	rn1 = 0.5_kindr
	rn2 = 0.5_kindr
    rn1c = 1. - rn1
    rn2c = 1. - rn2
    rn12 = rn1*rn2
    rn1c2 = rn1c*rn2
    rn12c = rn1*rn2c
    rn1c2c = rn1c*rn2c
    xe = xep(1,1)*rn1c2c+xep(1,2)*rn1c2+xep(1,3)*rn12+xep(1,4)*rn12c
    ye = xep(2,1)*rn1c2c+xep(2,2)*rn1c2+xep(2,3)*rn12+xep(2,4)*rn12c
    ze = xep(3,1)*rn1c2c+xep(3,2)*rn1c2+xep(3,3)*rn12+xep(3,4)*rn12c
!
!  Emission Plane Directions
!
    ex = ex_plane
    ey = ey_plane
    ez = ez_plane
!
!  Error correction.  Jeff Miller, 12/97
!  These lines were in the correct position in the subroutine intsec, but
!    out of place in this subroutine.  See comments above.
!
    slopex = -ex/ez
    slopey = -ey/ez
!
!  Intersection w/ aperture plane
!
!	IF (np==41) THEN
!	  WRITE(8,*) ' '
!    END IF
    xinta = slopex*ze + xe
    yinta = slopey*ze + ye
    if ( (xinta-xmina)*(xinta-xmaxa) <= 0. .and. &
         (yinta-ymina)*(yinta-ymaxa) <= 0. ) then
        npabs_a(nd) = npabs_a(nd) + 1
     end if
     if (idebug < 0) then
       write (8, 3002) np, xe, ye, ze, ex, ey, ez
    end if
!
!  Here for reflection

    nc_save = 0
    reflect = -1
 100 reflect= reflect + 1!
!    787 FORMAT(' ',2i5,3f10.5)
! 	IF (nploop==3) WRITE(10,787) np,int_no-2,ex,ey,ez
!
!  constants with emission direction
!
    ey2 = ey**2
    ez2 = ez**2
    a = ey2 + ez2
!
!  Intersection with Cylinders
!
!  compute possible intersection with all cylinders
!   - int_no 
!         =-1 for lost
!         = 0 for top plane
!         = 1 for bottom plane
!         = 2 for backplane
!         = int0_cyl to int0_cyl+num_cyls for cylinders
!         = int0_glass to int0_glass+num_cyls for glass
!         = int0_surf to int0_surf+num_surfs for surfaces
!
    int_no = -1
    dzint_old = 1.e30
    xint = 0.
    yint = 0.
    zint = -1.
!
!  Intersection with Cylinders
!
    call intsec_cyls
    if (nc_int /= 0) then
      nc_save = nc_int
      int_no = nc_int + int0_cyl - 1
      if (idebug < 0) then
        write (8, 3018) int_no, xint, yint, zint, dzint_old
      end if
    end if
!
!  Intersection with Surfaces
!
    if (num_surfs > 0) then
      call intsec_surfs
      if (idebug < 0) then
        write (8, 3019) int_no, xint, yint, zint, dzint_old
      end if
    end if
!
!  Intersection w/ Backplane
!
!    int_no = 0: top plane
!        = 1: bottom plane
!        = 2: back plane
!
    if (ez < 0.) then
!
!      - intersection with bottom plane
!
      if (int_no < 0) then
        xint = -ex*ze/ez + xe
        yint = -ey*ze/ez + ye
        zint = 0.
        if ( (xint-xminb)*(xint-xmaxb) <= 0. .and. &
                 (yint-yminb)*(yint-ymaxb) <= 0. ) then
          int_no = 2
		  nc_save=0    !this line corrects an error.  Jeff Miller 12/97
        else
          int_no = 1
          npabs (int_no, nd) = npabs(int_no, nd) + 1
          if (idebug < 0) then
            write (8, 3012) int_no
          end if
          cycle
        end if
        if (idebug < 0) then
          write (8, 3012) int_no
        end if
      end if
!
!      - intersection with top plane
!
    else
      if (int_no < 0) then
        int_no = 0
        if (idebug < 0) then
          write (8, 3013) int_no
        end if
        npabs (int_no, nd) = npabs(int_no, nd) + 1
        cycle
      end if
    end if
!
! Final Accumulation and Test
!
    if (idebug < 0) then
      write (8, 3030) int_no, xint, yint, zint, dzint_old
    end if
!
!  Photon Interaction ***************************************************
!
    if (int_no < 0) then
!
!  Lost
!
      nplost(nd) = nplost(nd) + 1
      write (*, 3090) nplost(nd), nd, nploop
      write (*,*) '   - XE(3), E(3) = ', xe, ye, ze, ex, ey, ez
      write (4, 3091) nd, np, nploop, nplost(nd),xe,ye,ze,xe+ex,ye+ey,ze+ez
      if (nplost(nd) > nplost_max) then
        write (*, 4000) nd, nploop, np, nplost(nd)
        write (9, 4000) nd, nploop, np, nplost(nd)
        call adios (2)
      end if
      go to 50
!
    elseif (endcap == 1) then
!
!  Absorbed by cylinder endcap (cylinder!)
!
      iabs = 0
!     
    elseif (int_no == 2) then
!
!  Back Plane
!
      nx = 0.
      ny = 0.
      nz = 1.
      call absref
!
    elseif (int_no >= int0_cyl .or. int_no <= int1_cyl) then
!
!  Cylinder
!
      nx = 0.
      ny = (yint - y0(nc_int))/ro(nc_int)
      nz = (zint - z0(nc_int))/ro(nc_int)
      if (idebug < -1) then
        write (8,*) 'nc_int, y0, z0, ro = ',nc_int, y0(nc_int), z0(nc_int), ro(nc_int)
      end if
      call absref
    else
!
!  Anything else - error
!
      write (*,*) ' -- error from intsec_ng: bad int_no =', int_no,' --- stopping'
      write (*,*) '    np, nploop, nd = ',np, nploop, nd
      write (*,*) '    XINT(3) = ', xe, ye, ze, ex, ey, ez
      write (*,*) '    XE(3), E(3) = ', xe, ye, ze, ex, ey, ez
      write (9,*) ' -- error from intsec_ng: bad int_no =', int_no,' --- stopping'
      write (9,*) '    np, nploop, nd = ',np, nploop, nd
      write (9,*) '    XINT(3) = ', xe, ye, ze, ex, ey, ez
      write (9,*) '    XE(3), E(3) = ', xe, ye, ze, ex, ey, ez
      stop
    end if
!
!  Photon disposition
!
    select case (iabs)
    case (0)
!
!  Absorbed
!
      npabs(int_no, nd) = npabs(int_no, nd) + 1
!	  WRITE(8,*) np,int_no-1
      cycle
    case default
!
!  Reflected
!
      xe = xint
      ye = yint
      ze = zint
      go to 100
    end select
!
  end do
!
  return
!
 3000 format (1x,'intsec_ng - ', &
     'nd,nploop= ', 2i5, 'ex,ey,ez= ', 3(1x,e11.4))
 3002 format (1x,'np - ',i5,2x,'Em pt = ', 3(1x,e11.4),', Dir = ',3(1x,e11.4))
 3012 format (1x,'int_no ',i5,' - int_no. back')
 3013 format (1x,'int_no =',i5,' - top plane')
 3014 format (5x,'np,      int_no,       xint,       yint,         zint')
 3015 format (3x,i5,3x,i5,3(1x,e11.4))
 3016 format (1x,'int_no =',i5,' - bot plane')
 3018 format (1x,'cyls:  int_no,xyz-int,dzint_old,=', i5, 4(1x,g12.5))
 3019 format (1x,'surfs: int_no,xyz-int,dzint_old,=', i5, 4(1x,g12.5))
 3020 format (1x,'int_no error: int_no,xyz-int,dzint_old,=', i5, 4(1x,g12.5))
 3030 format (1x,'intsec: int_no,xyz-int,dzint_old= ',i5,4(1x,g12.5))
 3090 format (1x,'Lost photon # ',i6,5x,'direction - ',i3,5x,'nploop - ',i3)
 3091 format (3i2,i5,6(1x,e11.4))
 4000 format (1x,'Error - lost photons exceeds specified limit', &
  / 5x,'at direction no. ',i5,'  - nploop no. ',i5,'  - np ',i5, &
    '  - no. lost this direction = ',i10)
!
end subroutine intsec_ng
!
!***********************************************************************
!
subroutine intsec_surfs
!
!***********************************************************************
!*    subroutine to perform intersection calculations with surfaces    *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use random
  use constants
  use directions
  use photons
  use planes
!
  implicit none
!
! integer(kind=kindi) :: ns
!
  write (*,*) 'Error - surfaces not yet implemented - stopping!'
  write (4,*) 'Error - surfaces not yet implemented - stopping!'
  write (9,*) 'Error - surfaces not yet implemented - stopping!'
  stop
! do ns = 1, num_surfs
! end do
!
end subroutine intsec_surfs
!
!***********************************************************************
!
subroutine print_npabs (iu)
!
!***********************************************************************
!*    subroutine to perform intersection calculations                  *
!***********************************************************************
!
!     - modules
!
  use kinds
  use control
  use photons
  use planes
!
  implicit none
!
  integer :: iu
  integer(kind = kindi) :: n
!
  call header (iu)
!
!  print oput npabs
!
  write (iu, 1000) nd, nploop, npabs_cyls(nd), npabs_glass(nd), &
     npabs_surf(nd), npabs_sum(nd), nplost(nd), npabs_tot(nd) 
!
!  back & top plane absorptances
!
  write (iu, 1012) (npabs(n, nd), n = 0, int0_cyl-1), npabs_a(nd)
!
!  write cylinder absorptions
!
  write (iu, 1001) 
  write (iu, 1011) (npabs(int0_cyl+n-1, nd), n = 1, num_cyls)
!
!  write glass absorptions
!
  if (iglass /= 0) then
    write (iu, 1002) 
    write (iu, 1011) (npabs(int0_glass+n-1, nd), n = 1, num_cyls)
  end if
!
!  Surfaces
!
  if (num_surfs > 0) then
    write (iu, 1004) 
    write (iu, 1011) (npabs(int0_surf+n-1, nd), n = 1, num_surfs)
  end if
!
  return
!
 1000 format (//5x,'P h o t o n   A b s o r p t i o n   A r r a y' &
  / 5x,'For direction  = ', i5,', nploop  = ',i5, &
  / 5x,'Photons absorbed for this direction by:' &
  / 7x,'Cylinders                        = ', i10,' photons' &
  / 7x,'Glass                            = ', i10,' photons' &
  / 7x,'Surfaces                         = ', i10,' photons' &
  / 7x,'Cylinders, glass and surfaces    = ', i10,' photons' &
  / 7x,'Lost                             = ', i10,' photons' &
  / 7x,'Total (including top and bottom) = ', i10,' photons'/) 
 1001 format (//5x,' C y l i n d e r   A b s o r p t a n c e s')
 1002 format (//5x,' G l a s s   A b s o r p t a n c e s')
 1004 format (//5x,' S u r f a c e   A b s o r p t a n c e s')
 1011 format (8i10)
 1012 format (//5x,' P l a n e   A b s o r p t a n c e s' &
  / 13x,'Top plane                  = ',i10,' photons' &
  / 13x,'Bottom plane               = ',i10,' photons' &
  / 13x,'Back plane                 = ',i10,' photons' &
  / 13x,'Aperture plane (virtual)   = ',i10,' photons' //)
!
end subroutine print_npabs


!
!*********************************************************************
!
subroutine tau_alpha_d
!
!*********************************************************************
!*     subroutine to calculate tau-alpha diffuse                     *
!*********************************************************************
!
  use control
  use photons
  use directions
  use constants
!
  implicit none
!
  integer(kind=kindi) :: i, j, k, l, m, q, t, ndelt
  real(kind=kindr) :: tad1, tad2, talp, ddelt, dthet, dthett, dthetl, sintt, &
                      costt, tantt, azm
  real(kind=kindr), allocatable, dimension(:,:) :: taul
!
  allocate (taul(n_long+1,n_tran+1))
!
!
! Make a 2-D array of tau-alpha, with a 90 degree buffer value:
!
  q = 0
  do i = 1, n_long
    do j = 1, n_tran
      q = q + 1
      taul(i,j) = tau_alpha(q)
    end do
    taul(i, n_tran+1) = taul(i, n_tran)
  end do
  do j = 1, n_tran
    taul(n_long+1,j) = taul(n_long,j)
  end do
  taul(n_long+1,n_tran+1) = taul(n_long,n_tran)
  thetatd(n_tran+1) = 90
  thetald(n_long+1) = 90
!
  tad2 = 0
  ndelt = 30
!
! Integrate tau_alpha over the azimuth (azm)
!
  do i = 1, ndelt
    tad1 = 0
    ddelt = pi2/ndelt
    azm = (i - 1)*ddelt + ddelt/2
!
! Integrate tau_alpha over the zenith (dthet)
!
    do j = 1, ndelt
      dthet = (j-1)*ddelt + ddelt/2
      sintt = sin(dthet)
      costt = cos(dthet)
      tantt = tan(dthet)
!
! Convert zenith and azimuth to transverse and longitudinal angles
!
      dthett = atan(tantt*sin(azm))*rtod
      dthetl = atan(tantt*cos(azm))*rtod
!
! Determine the interval dthett & dtehtl fall into
!
      t = 0
      l = 0
      do k = 1, n_tran
        if (dthett >= thetatd(k) .and. dthett < thetatd(k+1)) then
          t = k
        end if
      end do
      if (t == 0) then 
        t = n_tran
      end if
      do m = 1, n_long
        if (dthetl >= thetald(m) .and. dthetl < thetald(m+1)) then
          l = m
        end if
      end do
      if (l == 0) then 
        l = n_long
      end if
!
! Use bilinear interpolation to determine tau-alpha(dthetl,dthett)
!
      a = (dthetl - thetald(l))/(thetald(l+1) - thetald(l))
      b = (dthett - thetatd(t))/(thetatd(t+1) - thetatd(t))
      talp = (1-a)*(1-b)*taul(l,t) + a*(1-b)*taul(l+1,t) + a*b*taul(l+1,t+1) &
           + (1-a)*b*taul(l,t+1)
      tad1 = tad1 + talp*costt*sintt*ddelt
    end do
    tad2 = tad2 + tad1*ddelt
  end do
  tad = 4*tad2/pi
!
return
!
end subroutine tau_alpha_d