module photon_reflection_mod

use bmad_struct
use spline_mod

! For a custom reflection calc: 
!   See the photon_reflection_init routine for an example of how to set up the reflection tables.
! Note: It is assumed that for each table that the energy(:) values are equally spaced.
!       Also must have angle(1) = 0 and the last angle(n) = 90.

integer, parameter, private :: n_cheb_term$ = 30

type :: diffuse_param_struct
  real(rp) x, y
  real(rp) lambda, c_norm, chx_norm
  type (spline_struct) prob_spline(0:50)
  integer n_pt_spline
end type

type cheb_diffuse_struct
  real(rp) cch_int(n_cheb_term$), cch(n_cheb_term$)
end type

real(rp), private, parameter :: converge = 1d-9, gmin = 0.001, gmax = 100
integer(rp), private, parameter :: maxsum = 1000, ismax = 20, bmax = 100
logical, private :: appsw = .true.

type diffuse_common_struct
  ! Eventually will eliminate Chebyshev code and this logical. Keep for now for cross-checking.
  logical :: use_spline_fit = .true. 
  real(rp) :: area_err_tol = 4d-3
end type

type (diffuse_common_struct), save :: diffuse_com

private output_specular_reflection_input_params
private zmmax, zzfi, zzfp, hzz, zbessi
private zbessi1, zbessi0, zzexp

contains

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine photon_reflection_std_surface_init (surface)
!
! Routine to initialize the standard proton reflection probability tables.
! The standard tables are for 10 nm C film on Al substrate.
! The surface roughness for diffuse scattering is 200 nm and the
! the surface roughness correlation length is 5.5 um.
!
! Output:
!   surface -- photon_reflect_surface_struct
!-

Subroutine photon_reflection_std_surface_init (surface)

implicit none

type (photon_reflect_surface_struct), target :: surface
type (photon_reflect_table_struct), pointer :: prt

integer :: n_energy, n_angles
integer i, j, k

! Parameters used in generating the reflection values

surface%surface_roughness_rms = 200d-9      ! sigma in Dugan's notation
surface%roughness_correlation_len = 5500d-9 ! T in Dugan's notation
surface%name = 'DEFAULT'
surface%description = '10 nm C film on Al'
surface%reflectivity_file = '<none>'

! There are four tables.

if (allocated (surface%table)) deallocate (surface%table)
allocate(surface%table(4))

! Table 1: 30 - 600 eV, 10eV steps.

prt => surface%table(1)

n_energy = 58;  n_angles = 16
allocate(prt%angle(n_angles), prt%energy(n_energy), prt%p_reflect_scratch(n_angles), prt%int1(n_energy))
allocate(prt%p_reflect(n_angles,n_energy))

prt%angle = [0, 1, 2, 3, 4, 5, 6, 7, 10, 15, 20, 30, 45, 60, 75, 90]

prt%energy = [(i, i = 30, 600, 10)]
prt%max_energy = prt%energy(n_energy)

! Angle:                        0         1         2         3         4         5         6         7        10        15        20        30        45        60        75        90
prt%p_reflect(:,  1) = [1.00, 0.961540, 0.924484, 0.888710, 0.854106, 0.820574, 0.788022, 0.756365, 0.666028, 0.527177, 0.397515, 0.190951, 0.048179, 0.021745, 0.056309, 0.076222]  !   30 eV
prt%p_reflect(:,  2) = [1.00, 0.967596, 0.936125, 0.905445, 0.875429, 0.845955, 0.816905, 0.788167, 0.702620, 0.552681, 0.369345, 0.150822, 0.010838, 0.008902, 0.027377, 0.033684]  !   40 eV
prt%p_reflect(:,  3) = [1.00, 0.969652, 0.940022, 0.910893, 0.882055, 0.853291, 0.824365, 0.795000, 0.699997, 0.477343, 0.309192, 0.094149, 0.002007, 0.004530, 0.008834, 0.009397]  !   50 eV
prt%p_reflect(:,  4) = [1.00, 0.966959, 0.934490, 0.902022, 0.868867, 0.833999, 0.795557, 0.750409, 0.608135, 0.406969, 0.244594, 0.056553, 0.000396, 0.001509, 0.001420, 0.001059]  !   60 eV
prt%p_reflect(:,  5) = [1.00, 0.954062, 0.909654, 0.866230, 0.823382, 0.780835, 0.738427, 0.696093, 0.569940, 0.371238, 0.208157, 0.037080, 0.000076, 0.000191, 0.000024, 0.000209]  !   70 eV
prt%p_reflect(:,  6) = [1.00, 0.960378, 0.921731, 0.883486, 0.845155, 0.806323, 0.766660, 0.725927, 0.596588, 0.370688, 0.182100, 0.020181, 0.000015, 0.000175, 0.001388, 0.002084]  !   80 eV
prt%p_reflect(:,  7) = [1.00, 0.961656, 0.924068, 0.886554, 0.848505, 0.809383, 0.768721, 0.726155, 0.585619, 0.331102, 0.135269, 0.008402, 0.000001, 0.000522, 0.001618, 0.001953]  !   90 eV
prt%p_reflect(:,  8) = [1.00, 0.961891, 0.924373, 0.886625, 0.847897, 0.807499, 0.764815, 0.719341, 0.564311, 0.282178, 0.092338, 0.002615, 0.000001, 0.000532, 0.000927, 0.000897]  !  100 eV
prt%p_reflect(:,  9) = [1.00, 0.962031, 0.924489, 0.886406, 0.846874, 0.805030, 0.760074, 0.711328, 0.539974, 0.232046, 0.057840, 0.000391, 0.000002, 0.000315, 0.000239, 0.000138]  !  110 eV
prt%p_reflect(:, 10) = [1.00, 0.962240, 0.924737, 0.886371, 0.846065, 0.802755, 0.755416, 0.703148, 0.513756, 0.183143, 0.032335, 0.000045, 0.000002, 0.000095, 0.000005, 0.000030]  !  120 eV
prt%p_reflect(:, 11) = [1.00, 0.963095, 0.926233, 0.888135, 0.847520, 0.803059, 0.753389, 0.697236, 0.485294, 0.133094, 0.014273, 0.000439, 0.000002, 0.000003, 0.000114, 0.000230]  !  130 eV
prt%p_reflect(:, 12) = [1.00, 0.963500, 0.926814, 0.888446, 0.846843, 0.800322, 0.747078, 0.685365, 0.444830, 0.087407, 0.004775, 0.000805, 0.000001, 0.000023, 0.000222, 0.000304]  !  140 eV
prt%p_reflect(:, 13) = [1.00, 0.963722, 0.927026, 0.888173, 0.845302, 0.796304, 0.738829, 0.670559, 0.398442, 0.052252, 0.000902, 0.000941, 0.000000, 0.000064, 0.000198, 0.000202]  !  150 eV
prt%p_reflect(:, 14) = [1.00, 0.963887, 0.927106, 0.887652, 0.843290, 0.791383, 0.728878, 0.652720, 0.345612, 0.027296, 0.000121, 0.000840, 0.000000, 0.000074, 0.000093, 0.000061]  !  160 eV
prt%p_reflect(:, 15) = [1.00, 0.963890, 0.926856, 0.886590, 0.840429, 0.785104, 0.716701, 0.631315, 0.289083, 0.012118, 0.000731, 0.000618, 0.000000, 0.000052, 0.000015, 0.000001]  !  170 eV
prt%p_reflect(:, 16) = [1.00, 0.964105, 0.927033, 0.886177, 0.838417, 0.779752, 0.705205, 0.609770, 0.233327, 0.004005, 0.001741, 0.000368, 0.000000, 0.000020, 0.000004, 0.000026]  !  180 eV
prt%p_reflect(:, 17) = [1.00, 0.964604, 0.927729, 0.886412, 0.836931, 0.774249, 0.691786, 0.582910, 0.171951, 0.000731, 0.002397, 0.000158, 0.000000, 0.000002, 0.000033, 0.000063]  !  190 eV
prt%p_reflect(:, 18) = [1.00, 0.964802, 0.927779, 0.885518, 0.833529, 0.765433, 0.672603, 0.546788, 0.115472, 0.000455, 0.002583, 0.000043, 0.000000, 0.000002, 0.000049, 0.000062]  !  200 eV
prt%p_reflect(:, 19) = [1.00, 0.964624, 0.927039, 0.883291, 0.827938, 0.752894, 0.647000, 0.500881, 0.069982, 0.001304, 0.002407, 0.000003, 0.000000, 0.000010, 0.000039, 0.000032]  !  210 eV
prt%p_reflect(:, 20) = [1.00, 0.963996, 0.925350, 0.879436, 0.819627, 0.735677, 0.613436, 0.444127, 0.037954, 0.002322, 0.002010, 0.000011, 0.000000, 0.000015, 0.000016, 0.000006]  !  220 eV
prt%p_reflect(:, 21) = [1.00, 0.962598, 0.922042, 0.872863, 0.806937, 0.711386, 0.569124, 0.376355, 0.018609, 0.003091, 0.001534, 0.000037, 0.000000, 0.000013, 0.000002, 0.000001]  !  230 eV
prt%p_reflect(:, 22) = [1.00, 0.960518, 0.917296, 0.863849, 0.790274, 0.680716, 0.516035, 0.304136, 0.008575, 0.003554, 0.001075, 0.000063, 0.000000, 0.000007, 0.000002, 0.000011]  !  240 eV
prt%p_reflect(:, 23) = [1.00, 0.957459, 0.910473, 0.851283, 0.767857, 0.641203, 0.452606, 0.231057, 0.004581, 0.003670, 0.000678, 0.000076, 0.000000, 0.000002, 0.000010, 0.000019]  !  250 eV
prt%p_reflect(:, 24) = [1.00, 0.952377, 0.899427, 0.831752, 0.734909, 0.587625, 0.377629, 0.163663, 0.004571, 0.003478, 0.000386, 0.000075, 0.000000, 0.000001, 0.000014, 0.000017]  !  260 eV
prt%p_reflect(:, 25) = [1.00, 0.940349, 0.874513, 0.791076, 0.674136, 0.505364, 0.290636, 0.111277, 0.007669, 0.003026, 0.000236, 0.000061, 0.000000, 0.000003, 0.000011, 0.000008]  !  270 eV
prt%p_reflect(:, 26) = [1.00, 0.869775, 0.747058, 0.623730, 0.494226, 0.356060, 0.215413, 0.104113, 0.015620, 0.002326, 0.000417, 0.000045, 0.000000, 0.000005, 0.000008, 0.000009]  !  280 eV
prt%p_reflect(:, 27) = [1.00, 0.614161, 0.360604, 0.189716, 0.079002, 0.020432, 0.009485, 0.020201, 0.012908, 0.000136, 0.000458, 0.000003, 0.000000, 0.000003, 0.000008, 0.000016]  !  290 eV
prt%p_reflect(:, 28) = [1.00, 0.698805, 0.472911, 0.296085, 0.157328, 0.058601, 0.014309, 0.014095, 0.010992, 0.000261, 0.000366, 0.000006, 0.000000, 0.000001, 0.000009, 0.000013]  !  300 eV
prt%p_reflect(:, 29) = [1.00, 0.731910, 0.520344, 0.344498, 0.195659, 0.078838, 0.020168, 0.015091, 0.009853, 0.000264, 0.000356, 0.000012, 0.000000, 0.000000, 0.000009, 0.000008]  !  310 eV
prt%p_reflect(:, 30) = [1.00, 0.751760, 0.549521, 0.374655, 0.218903, 0.089142, 0.023345, 0.016451, 0.008874, 0.000258, 0.000344, 0.000018, 0.000000, 0.000001, 0.000006, 0.000003]  !  320 eV
prt%p_reflect(:, 31) = [1.00, 0.765442, 0.569822, 0.395401, 0.233564, 0.092894, 0.024798, 0.017471, 0.007947, 0.000263, 0.000321, 0.000023, 0.000000, 0.000002, 0.000003, 0.000001]  !  330 eV
prt%p_reflect(:, 32) = [1.00, 0.775695, 0.585046, 0.410532, 0.242650, 0.092061, 0.025400, 0.018079, 0.007056, 0.000283, 0.000288, 0.000024, 0.000000, 0.000003, 0.000001, 0.000002]  !  340 eV
prt%p_reflect(:, 33) = [1.00, 0.783794, 0.597013, 0.421916, 0.247670, 0.087927, 0.025632, 0.018320, 0.006205, 0.000314, 0.000248, 0.000024, 0.000000, 0.000003, 0.000002, 0.000004]  !  350 eV
prt%p_reflect(:, 34) = [1.00, 0.790447, 0.606749, 0.430628, 0.249478, 0.081640, 0.025727, 0.018251, 0.005403, 0.000351, 0.000206, 0.000021, 0.000000, 0.000002, 0.000003, 0.000005]  !  360 eV
prt%p_reflect(:, 35) = [1.00, 0.796128, 0.614976, 0.437472, 0.248806, 0.074561, 0.025673, 0.017884, 0.004670, 0.000386, 0.000166, 0.000017, 0.000000, 0.000001, 0.000004, 0.000004]  !  370 eV
prt%p_reflect(:, 36) = [1.00, 0.801476, 0.622713, 0.443634, 0.246644, 0.067698, 0.025416, 0.017298, 0.004020, 0.000414, 0.000131, 0.000013, 0.000000, 0.000001, 0.000003, 0.000002]  !  380 eV
prt%p_reflect(:, 37) = [1.00, 0.806586, 0.630078, 0.449144, 0.242681, 0.061168, 0.024772, 0.016459, 0.003445, 0.000429, 0.000102, 0.000009, 0.000000, 0.000001, 0.000002, 0.000001]  !  390 eV
prt%p_reflect(:, 38) = [1.00, 0.810943, 0.636204, 0.452929, 0.235850, 0.055084, 0.023860, 0.015499, 0.002944, 0.000433, 0.000080, 0.000007, 0.000000, 0.000001, 0.000001, 0.000001]  !  400 eV
prt%p_reflect(:, 39) = [1.00, 0.814742, 0.641399, 0.455334, 0.226400, 0.049780, 0.022805, 0.014489, 0.002514, 0.000427, 0.000065, 0.000005, 0.000000, 0.000001, 0.000001, 0.000002]  !  410 eV
prt%p_reflect(:, 40) = [1.00, 0.818150, 0.645935, 0.456663, 0.214569, 0.045302, 0.021678, 0.013470, 0.002151, 0.000413, 0.000055, 0.000005, 0.000000, 0.000001, 0.000001, 0.000003]  !  420 eV
prt%p_reflect(:, 41) = [1.00, 0.821200, 0.649852, 0.456919, 0.200383, 0.041459, 0.020495, 0.012463, 0.001846, 0.000392, 0.000050, 0.000005, 0.000000, 0.000001, 0.000002, 0.000002]  !  430 eV
prt%p_reflect(:, 42) = [1.00, 0.823997, 0.653327, 0.456299, 0.184160, 0.038126, 0.019313, 0.011497, 0.001591, 0.000367, 0.000048, 0.000006, 0.000000, 0.000001, 0.000002, 0.000002]  !  440 eV
prt%p_reflect(:, 43) = [1.00, 0.826548, 0.656357, 0.454705, 0.166109, 0.035247, 0.018165, 0.010578, 0.001379, 0.000337, 0.000048, 0.000006, 0.000000, 0.000000, 0.000001, 0.000001]  !  450 eV
prt%p_reflect(:, 44) = [1.00, 0.828891, 0.658997, 0.452143, 0.147172, 0.032754, 0.017043, 0.009700, 0.001205, 0.000305, 0.000050, 0.000006, 0.000000, 0.000000, 0.000001, 0.000001]  !  460 eV
prt%p_reflect(:, 45) = [1.00, 0.831069, 0.661317, 0.448655, 0.128592, 0.030548, 0.015951, 0.008869, 0.001065, 0.000273, 0.000051, 0.000006, 0.000000, 0.000001, 0.000001, 0.000001]  !  470 eV
prt%p_reflect(:, 46) = [1.00, 0.833095, 0.663334, 0.444176, 0.111497, 0.028559, 0.014894, 0.008088, 0.000953, 0.000241, 0.000052, 0.000006, 0.000000, 0.000001, 0.000001, 0.000001]  !  480 eV
prt%p_reflect(:, 47) = [1.00, 0.835065, 0.665217, 0.438839, 0.096707, 0.026742, 0.013874, 0.007358, 0.000864, 0.000211, 0.000052, 0.000005, 0.000000, 0.000001, 0.000001, 0.000001]  !  490 eV
prt%p_reflect(:, 48) = [1.00, 0.836962, 0.666926, 0.432474, 0.084299, 0.025041, 0.012892, 0.006681, 0.000794, 0.000183, 0.000051, 0.000004, 0.000000, 0.000000, 0.000001, 0.000001]  !  500 eV
prt%p_reflect(:, 49) = [1.00, 0.838726, 0.668352, 0.424812, 0.074036, 0.023454, 0.011959, 0.006057, 0.000738, 0.000159, 0.000049, 0.000004, 0.000000, 0.000000, 0.000001, 0.000001]  !  510 eV
prt%p_reflect(:, 50) = [1.00, 0.840371, 0.669505, 0.415681, 0.065628, 0.021974, 0.011075, 0.005482, 0.000694, 0.000137, 0.000045, 0.000003, 0.000000, 0.000000, 0.000001, 0.000001]  !  520 eV
prt%p_reflect(:, 51) = [1.00, 0.841918, 0.670426, 0.404956, 0.058737, 0.020588, 0.010240, 0.004956, 0.000658, 0.000120, 0.000041, 0.000002, 0.000000, 0.000000, 0.000001, 0.000001]  !  530 eV
prt%p_reflect(:, 52) = [1.00, 0.843377, 0.671124, 0.392425, 0.053025, 0.019284, 0.009454, 0.004477, 0.000629, 0.000105, 0.000037, 0.000002, 0.000000, 0.000000, 0.000000, 0.000001]  !  540 eV
prt%p_reflect(:, 53) = [1.00, 0.844758, 0.671612, 0.377846, 0.048218, 0.018055, 0.008719, 0.004043, 0.000605, 0.000094, 0.000032, 0.000002, 0.000000, 0.000000, 0.000001, 0.000001]  !  550 eV
prt%p_reflect(:, 54) = [1.00, 0.846070, 0.671905, 0.360935, 0.044108, 0.016901, 0.008034, 0.003652, 0.000584, 0.000085, 0.000028, 0.000002, 0.000000, 0.000000, 0.000001, 0.000001]  !  560 eV
prt%p_reflect(:, 55) = [1.00, 0.847317, 0.671998, 0.341296, 0.040565, 0.015820, 0.007396, 0.003298, 0.000564, 0.000079, 0.000024, 0.000002, 0.000000, 0.000000, 0.000001, 0.000000]  !  570 eV
prt%p_reflect(:, 56) = [1.00, 0.848501, 0.671887, 0.318607, 0.037478, 0.014804, 0.006804, 0.002980, 0.000546, 0.000075, 0.000021, 0.000002, 0.000000, 0.000000, 0.000000, 0.000000]  !  580 eV
prt%p_reflect(:, 57) = [1.00, 0.849626, 0.671571, 0.292731, 0.034765, 0.013849, 0.006253, 0.002694, 0.000528, 0.000072, 0.000019, 0.000002, 0.000000, 0.000000, 0.000000, 0.000000]  !  590 eV
prt%p_reflect(:, 58) = [1.00, 0.850699, 0.671053, 0.264038, 0.032353, 0.012949, 0.005743, 0.002438, 0.000509, 0.000070, 0.000017, 0.000002, 0.000000, 0.000000, 0.000000, 0.000001]  !  600 eV

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! Table 2: 600 - 1400 eV, 20 eV steps.

prt => surface%table(2)

n_energy = 41;  n_angles = 11
allocate(prt%angle(n_angles), prt%energy(n_energy), prt%p_reflect_scratch(n_angles), prt%int1(n_energy))
allocate(prt%p_reflect(n_angles,n_energy))

prt%angle = [0.0, 0.4, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 90.0]

prt%energy = [(i, i = 600, 1400, 20)]
prt%max_energy = prt%energy(n_energy)

! Angle:                         0        0.4       0.8       1.0       1.5      2.0        2.5      3.0       3.5       4.0       90 
prt%p_reflect(:,  1) = [1.00, 0.940114, 0.880951, 0.850699, 0.769443, 0.671053, 0.530219, 0.264038, 0.068672, 0.032353, 0.000001]  !  600 eV
prt%p_reflect(:,  2) = [1.00, 0.941155, 0.882746, 0.852699, 0.771058, 0.669398, 0.514744, 0.203619, 0.056321, 0.028229, 0.000000]  !  620 eV
prt%p_reflect(:,  3) = [1.00, 0.942127, 0.884409, 0.854534, 0.772371, 0.666912, 0.494590, 0.151924, 0.047344, 0.024806, 0.000000]  !  640 eV
prt%p_reflect(:,  4) = [1.00, 0.943038, 0.885954, 0.856218, 0.773392, 0.663511, 0.467916, 0.115279, 0.040559, 0.021929, 0.000000]  !  660 eV
prt%p_reflect(:,  5) = [1.00, 0.943827, 0.887258, 0.857597, 0.773863, 0.658660, 0.431035, 0.090770, 0.035295, 0.019462, 0.000000]  !  680 eV
prt%p_reflect(:,  6) = [1.00, 0.944466, 0.888270, 0.858610, 0.773693, 0.652129, 0.380397, 0.074279, 0.031097, 0.017323, 0.000000]  !  700 eV
prt%p_reflect(:,  7) = [1.00, 0.945419, 0.889892, 0.860374, 0.774565, 0.646055, 0.316107, 0.062636, 0.027615, 0.015428, 0.000000]  !  720 eV
prt%p_reflect(:,  8) = [1.00, 0.946435, 0.891635, 0.862282, 0.775539, 0.638798, 0.244749, 0.053638, 0.024536, 0.013699, 0.000000]  !  740 eV
prt%p_reflect(:,  9) = [1.00, 0.947375, 0.893230, 0.864001, 0.776100, 0.629252, 0.185451, 0.046495, 0.021842, 0.012156, 0.000000]  !  760 eV
prt%p_reflect(:, 10) = [1.00, 0.948259, 0.894715, 0.865576, 0.776304, 0.616912, 0.144032, 0.040686, 0.019477, 0.010785, 0.000000]  !  780 eV
prt%p_reflect(:, 11) = [1.00, 0.949098, 0.896109, 0.867031, 0.776165, 0.600875, 0.115583, 0.035881, 0.017399, 0.009569, 0.000000]  !  800 eV
prt%p_reflect(:, 12) = [1.00, 0.949893, 0.897416, 0.868369, 0.775659, 0.579615, 0.095364, 0.031839, 0.015563, 0.008492, 0.000000]  !  820 eV
prt%p_reflect(:, 13) = [1.00, 0.950651, 0.898646, 0.869602, 0.774776, 0.550616, 0.080423, 0.028395, 0.013937, 0.007539, 0.000000]  !  840 eV
prt%p_reflect(:, 14) = [1.00, 0.951386, 0.899829, 0.870767, 0.773537, 0.509513, 0.068983, 0.025422, 0.012491, 0.006695, 0.000000]  !  860 eV
prt%p_reflect(:, 15) = [1.00, 0.952098, 0.900965, 0.871863, 0.771898, 0.449163, 0.059953, 0.022829, 0.011201, 0.005948, 0.000000]  !  880 eV
prt%p_reflect(:, 16) = [1.00, 0.952779, 0.902035, 0.872864, 0.769759, 0.366023, 0.052643, 0.020551, 0.010049, 0.005289, 0.000000]  !  900 eV
prt%p_reflect(:, 17) = [1.00, 0.953432, 0.903045, 0.873778, 0.767069, 0.281119, 0.046611, 0.018536, 0.009021, 0.004709, 0.000000]  !  920 eV
prt%p_reflect(:, 18) = [1.00, 0.954069, 0.904020, 0.874634, 0.763786, 0.217636, 0.041556, 0.016748, 0.008102, 0.004198, 0.000000]  !  940 eV
prt%p_reflect(:, 19) = [1.00, 0.954683, 0.904945, 0.875412, 0.759765, 0.174103, 0.037251, 0.015150, 0.007281, 0.003749, 0.000000]  !  960 eV
prt%p_reflect(:, 20) = [1.00, 0.955274, 0.905819, 0.876109, 0.754841, 0.143481, 0.033551, 0.013720, 0.006545, 0.003355, 0.000000]  !  980 eV
prt%p_reflect(:, 21) = [1.00, 0.955843, 0.906643, 0.876723, 0.748825, 0.121036, 0.030326, 0.012434, 0.005889, 0.003010, 0.000000]  ! 1000 eV
prt%p_reflect(:, 22) = [1.00, 0.956387, 0.907410, 0.877245, 0.741378, 0.103975, 0.027505, 0.011277, 0.005302, 0.002709, 0.000000]  ! 1020 eV
prt%p_reflect(:, 23) = [1.00, 0.956917, 0.908142, 0.877701, 0.732152, 0.090617, 0.025009, 0.010232, 0.004778, 0.002446, 0.000000]  ! 1040 eV
prt%p_reflect(:, 24) = [1.00, 0.957454, 0.908885, 0.878147, 0.720636, 0.079869, 0.022788, 0.009288, 0.004311, 0.002218, 0.000000]  ! 1060 eV
prt%p_reflect(:, 25) = [1.00, 0.957980, 0.909596, 0.878529, 0.705721, 0.071028, 0.020797, 0.008434, 0.003896, 0.002019, 0.000000]  ! 1080 eV
prt%p_reflect(:, 26) = [1.00, 0.958486, 0.910260, 0.878819, 0.685598, 0.063646, 0.019007, 0.007660, 0.003526, 0.001847, 0.000000]  ! 1100 eV
prt%p_reflect(:, 27) = [1.00, 0.958975, 0.910879, 0.879019, 0.656947, 0.057404, 0.017391, 0.006958, 0.003198, 0.001697, 0.000000]  ! 1120 eV
prt%p_reflect(:, 28) = [1.00, 0.959444, 0.911448, 0.879117, 0.613065, 0.052052, 0.015927, 0.006323, 0.002907, 0.001567, 0.000000]  ! 1140 eV
prt%p_reflect(:, 29) = [1.00, 0.959896, 0.911972, 0.879114, 0.540959, 0.047417, 0.014594, 0.005747, 0.002650, 0.001454, 0.000000]  ! 1160 eV
prt%p_reflect(:, 30) = [1.00, 0.960332, 0.912450, 0.879003, 0.436458, 0.043377, 0.013381, 0.005225, 0.002422, 0.001357, 0.000000]  ! 1180 eV
prt%p_reflect(:, 31) = [1.00, 0.960752, 0.912883, 0.878780, 0.341217, 0.039814, 0.012270, 0.004752, 0.002223, 0.001272, 0.000000]  ! 1200 eV
prt%p_reflect(:, 32) = [1.00, 0.961158, 0.913271, 0.878434, 0.274941, 0.036660, 0.011252, 0.004324, 0.002047, 0.001198, 0.000000]  ! 1220 eV
prt%p_reflect(:, 33) = [1.00, 0.961554, 0.913622, 0.877968, 0.229167, 0.033839, 0.010316, 0.003937, 0.001894, 0.001134, 0.000000]  ! 1240 eV
prt%p_reflect(:, 34) = [1.00, 0.961945, 0.913943, 0.877377, 0.195997, 0.031312, 0.009454, 0.003587, 0.001760, 0.001077, 0.000000]  ! 1260 eV
prt%p_reflect(:, 35) = [1.00, 0.962331, 0.914232, 0.876641, 0.170961, 0.029030, 0.008658, 0.003271, 0.001644, 0.001027, 0.000000]  ! 1280 eV
prt%p_reflect(:, 36) = [1.00, 0.962709, 0.914481, 0.875733, 0.151438, 0.026949, 0.007921, 0.002988, 0.001545, 0.000982, 0.000000]  ! 1300 eV
prt%p_reflect(:, 37) = [1.00, 0.963070, 0.914667, 0.874592, 0.135821, 0.025055, 0.007239, 0.002734, 0.001459, 0.000941, 0.000000]  ! 1320 eV
prt%p_reflect(:, 38) = [1.00, 0.963416, 0.914786, 0.873179, 0.123102, 0.023322, 0.006606, 0.002506, 0.001385, 0.000904, 0.000000]  ! 1340 eV
prt%p_reflect(:, 39) = [1.00, 0.963747, 0.914835, 0.871442, 0.112591, 0.021728, 0.006018, 0.002303, 0.001323, 0.000868, 0.000000]  ! 1360 eV
prt%p_reflect(:, 40) = [1.00, 0.964060, 0.914798, 0.869281, 0.103837, 0.020264, 0.005470, 0.002123, 0.001271, 0.000835, 0.000000]  ! 1380 eV
prt%p_reflect(:, 41) = [1.00, 0.964356, 0.914668, 0.866596, 0.096476, 0.018906, 0.004959, 0.001964, 0.001228, 0.000803, 0.000000]  ! 1400 eV

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! Table 3: 1400 - 1600 eV, 10 eV steps. There is a resonance here

prt => surface%table(3)

n_energy = 21;  n_angles = 10
allocate(prt%angle(n_angles), prt%energy(n_energy), prt%p_reflect_scratch(n_angles), prt%int1(n_energy))
allocate(prt%p_reflect(n_angles,n_energy))

prt%angle = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 90.0]

prt%energy = [(i, i = 1400, 1600, 10)]
prt%max_energy = prt%energy(n_energy)

! Angle:                       0        0.2       0.4       0.6       0.8       1.0       1.5      2.0        2.5      90 
prt%p_reflect(:,  1) = [1.00, 0.982718, 0.964356, 0.942985, 0.914668, 0.866596, 0.096476, 0.018906, 0.004959, 0.000000]  ! 1400 eV
prt%p_reflect(:,  2) = [1.00, 0.982800, 0.964495, 0.943112, 0.914549, 0.864953, 0.093305, 0.018277, 0.004714, 0.000000]  ! 1410 eV
prt%p_reflect(:,  3) = [1.00, 0.982880, 0.964632, 0.943235, 0.914412, 0.863146, 0.090338, 0.017658, 0.004478, 0.000000]  ! 1420 eV
prt%p_reflect(:,  4) = [1.00, 0.982957, 0.964760, 0.943338, 0.914224, 0.861016, 0.087696, 0.017074, 0.004247, 0.000000]  ! 1430 eV
prt%p_reflect(:,  5) = [1.00, 0.983032, 0.964884, 0.943431, 0.913999, 0.858586, 0.085263, 0.016505, 0.004023, 0.000000]  ! 1440 eV
prt%p_reflect(:,  6) = [1.00, 0.983104, 0.965001, 0.943510, 0.913728, 0.855753, 0.083055, 0.015957, 0.003805, 0.000000]  ! 1450 eV
prt%p_reflect(:,  7) = [1.00, 0.983172, 0.965109, 0.943566, 0.913385, 0.852316, 0.081128, 0.015436, 0.003589, 0.000000]  ! 1460 eV
prt%p_reflect(:,  8) = [1.00, 0.983238, 0.965213, 0.943611, 0.912998, 0.848290, 0.079330, 0.014922, 0.003382, 0.000000]  ! 1470 eV
prt%p_reflect(:,  9) = [1.00, 0.983297, 0.965296, 0.943609, 0.912456, 0.842829, 0.077999, 0.014455, 0.003167, 0.000000]  ! 1480 eV
prt%p_reflect(:, 10) = [1.00, 0.983354, 0.965376, 0.943595, 0.911856, 0.836140, 0.076724, 0.013985, 0.002963, 0.000000]  ! 1490 eV
prt%p_reflect(:, 11) = [1.00, 0.983400, 0.965428, 0.943510, 0.911007, 0.825985, 0.076018, 0.013568, 0.002745, 0.000000]  ! 1500 eV
prt%p_reflect(:, 12) = [1.00, 0.983441, 0.965464, 0.943383, 0.909971, 0.810854, 0.075529, 0.013161, 0.002529, 0.000000]  ! 1510 eV
prt%p_reflect(:, 13) = [1.00, 0.983465, 0.965458, 0.943144, 0.908497, 0.781664, 0.075642, 0.012800, 0.002299, 0.000000]  ! 1520 eV
prt%p_reflect(:, 14) = [1.00, 0.983450, 0.965347, 0.942632, 0.905908, 0.689869, 0.077146, 0.012552, 0.002021, 0.000000]  ! 1530 eV
prt%p_reflect(:, 15) = [1.00, 0.983429, 0.965216, 0.942039, 0.902673, 0.569241, 0.078471, 0.012267, 0.001768, 0.000000]  ! 1540 eV
prt%p_reflect(:, 16) = [1.00, 0.982392, 0.962190, 0.931885, 0.737868, 0.428971, 0.103437, 0.014475, 0.000726, 0.000000]  ! 1550 eV
prt%p_reflect(:, 17) = [1.00, 0.965214, 0.918122, 0.836670, 0.698236, 0.516134, 0.135915, 0.012950, 0.000235, 0.000000]  ! 1560 eV
prt%p_reflect(:, 18) = [1.00, 0.975932, 0.946573, 0.900328, 0.802683, 0.583571, 0.094054, 0.006236, 0.001449, 0.000000]  ! 1570 eV
prt%p_reflect(:, 19) = [1.00, 0.976822, 0.948812, 0.905530, 0.815379, 0.598498, 0.084473, 0.005185, 0.001776, 0.000000]  ! 1580 eV
prt%p_reflect(:, 20) = [1.00, 0.977607, 0.950771, 0.910038, 0.826798, 0.616014, 0.075247, 0.004371, 0.002119, 0.000000]  ! 1590 eV
prt%p_reflect(:, 21) = [1.00, 0.978015, 0.951732, 0.912041, 0.831275, 0.621616, 0.068842, 0.003919, 0.002316, 0.000000]  ! 1600 eV

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! Table 4: 1600 - 5000 eV, 50 eV steps

prt => surface%table(4)

n_energy = 69;  n_angles = 11
allocate(prt%angle(n_angles), prt%energy(n_energy), prt%p_reflect_scratch(n_angles), prt%int1(n_energy))
allocate(prt%p_reflect(n_angles,n_energy))

prt%angle = [0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 90.0]

prt%energy = [(i, i = 1600, 5000, 50)]
prt%max_energy = prt%energy(n_energy)

! Angle:                         0        0.2       0.4       0.5       0.6       0.8       1.0       1.5      2.0        2.5      90 
prt%p_reflect(:,  1) = [1.00, 0.978015, 0.951732, 0.934539, 0.912041, 0.831275, 0.621616, 0.068842, 0.003919, 0.002316, 0.000000]  ! 1600 eV
prt%p_reflect(:,  2) = [1.00, 0.979215, 0.954349, 0.938065, 0.916709, 0.838937, 0.619354, 0.045984, 0.002991, 0.002847, 0.000000]  ! 1650 eV
prt%p_reflect(:,  3) = [1.00, 0.979960, 0.955803, 0.939809, 0.918578, 0.838602, 0.595077, 0.031175, 0.002931, 0.002998, 0.000000]  ! 1700 eV
prt%p_reflect(:,  4) = [1.00, 0.980554, 0.956879, 0.940976, 0.919526, 0.835186, 0.558223, 0.021199, 0.003152, 0.002939, 0.000000]  ! 1750 eV
prt%p_reflect(:,  5) = [1.00, 0.981073, 0.957769, 0.941851, 0.919984, 0.829601, 0.510524, 0.014581, 0.003431, 0.002745, 0.000000]  ! 1800 eV
prt%p_reflect(:,  6) = [1.00, 0.981532, 0.958509, 0.942491, 0.920034, 0.821829, 0.453276, 0.010316, 0.003662, 0.002471, 0.000000]  ! 1850 eV
prt%p_reflect(:,  7) = [1.00, 0.981948, 0.959141, 0.942958, 0.919764, 0.811816, 0.389671, 0.007647, 0.003808, 0.002159, 0.000000]  ! 1900 eV
prt%p_reflect(:,  8) = [1.00, 0.982326, 0.959680, 0.943268, 0.919180, 0.799092, 0.323839, 0.006064, 0.003854, 0.001837, 0.000000]  ! 1950 eV
prt%p_reflect(:,  9) = [1.00, 0.982675, 0.960146, 0.943452, 0.918321, 0.783266, 0.261140, 0.005185, 0.003809, 0.001527, 0.000000]  ! 2000 eV
prt%p_reflect(:, 10) = [1.00, 0.983003, 0.960557, 0.943535, 0.917211, 0.763692, 0.205437, 0.004759, 0.003685, 0.001242, 0.000000]  ! 2050 eV
prt%p_reflect(:, 11) = [1.00, 0.983319, 0.960934, 0.943537, 0.915848, 0.739163, 0.158421, 0.004617, 0.003490, 0.000992, 0.000000]  ! 2100 eV
prt%p_reflect(:, 12) = [1.00, 0.983619, 0.961266, 0.943444, 0.914200, 0.708464, 0.120663, 0.004627, 0.003248, 0.000780, 0.000000]  ! 2150 eV
prt%p_reflect(:, 13) = [1.00, 0.983905, 0.961557, 0.943257, 0.912239, 0.669936, 0.091222, 0.004708, 0.002972, 0.000609, 0.000000]  ! 2200 eV
prt%p_reflect(:, 14) = [1.00, 0.984175, 0.961802, 0.942961, 0.909902, 0.621502, 0.068737, 0.004806, 0.002677, 0.000478, 0.000000]  ! 2250 eV
prt%p_reflect(:, 15) = [1.00, 0.984433, 0.962008, 0.942565, 0.907172, 0.561991, 0.051798, 0.004888, 0.002378, 0.000381, 0.000000]  ! 2300 eV
prt%p_reflect(:, 16) = [1.00, 0.984679, 0.962173, 0.942054, 0.903966, 0.491459, 0.039173, 0.004932, 0.002083, 0.000316, 0.000000]  ! 2350 eV
prt%p_reflect(:, 17) = [1.00, 0.984915, 0.962306, 0.941441, 0.900259, 0.413988, 0.029810, 0.004935, 0.001801, 0.000276, 0.000000]  ! 2400 eV
prt%p_reflect(:, 18) = [1.00, 0.985138, 0.962392, 0.940689, 0.895885, 0.335868, 0.022950, 0.004885, 0.001539, 0.000257, 0.000000]  ! 2450 eV
prt%p_reflect(:, 19) = [1.00, 0.985355, 0.962453, 0.939831, 0.890850, 0.265223, 0.017912, 0.004794, 0.001301, 0.000253, 0.000000]  ! 2500 eV
prt%p_reflect(:, 20) = [1.00, 0.985562, 0.962472, 0.938823, 0.884920, 0.205737, 0.014274, 0.004658, 0.001089, 0.000259, 0.000000]  ! 2550 eV
prt%p_reflect(:, 21) = [1.00, 0.985759, 0.962451, 0.937657, 0.877925, 0.158308, 0.011670, 0.004483, 0.000905, 0.000271, 0.000000]  ! 2600 eV
prt%p_reflect(:, 22) = [1.00, 0.985949, 0.962393, 0.936322, 0.869647, 0.121635, 0.009827, 0.004278, 0.000748, 0.000285, 0.000000]  ! 2650 eV
prt%p_reflect(:, 23) = [1.00, 0.986133, 0.962299, 0.934810, 0.859811, 0.093689, 0.008539, 0.004050, 0.000618, 0.000299, 0.000000]  ! 2700 eV
prt%p_reflect(:, 24) = [1.00, 0.986308, 0.962162, 0.933077, 0.847927, 0.072487, 0.007668, 0.003803, 0.000514, 0.000310, 0.000000]  ! 2750 eV
prt%p_reflect(:, 25) = [1.00, 0.986478, 0.961994, 0.931138, 0.833625, 0.056428, 0.007083, 0.003546, 0.000432, 0.000317, 0.000000]  ! 2800 eV
prt%p_reflect(:, 26) = [1.00, 0.986639, 0.961770, 0.928897, 0.815837, 0.044272, 0.006721, 0.003280, 0.000371, 0.000319, 0.000000]  ! 2850 eV
prt%p_reflect(:, 27) = [1.00, 0.986796, 0.961507, 0.926364, 0.793720, 0.035051, 0.006503, 0.003013, 0.000327, 0.000315, 0.000000]  ! 2900 eV
prt%p_reflect(:, 28) = [1.00, 0.986947, 0.961203, 0.923506, 0.765923, 0.028033, 0.006382, 0.002751, 0.000299, 0.000306, 0.000000]  ! 2950 eV
prt%p_reflect(:, 29) = [1.00, 0.987092, 0.960845, 0.920229, 0.730135, 0.022712, 0.006327, 0.002494, 0.000283, 0.000293, 0.000000]  ! 3000 eV
prt%p_reflect(:, 30) = [1.00, 0.987231, 0.960433, 0.916469, 0.683649, 0.018680, 0.006311, 0.002247, 0.000276, 0.000275, 0.000000]  ! 3050 eV
prt%p_reflect(:, 31) = [1.00, 0.987365, 0.959965, 0.912152, 0.623507, 0.015627, 0.006312, 0.002012, 0.000277, 0.000253, 0.000000]  ! 3100 eV
prt%p_reflect(:, 32) = [1.00, 0.987495, 0.959443, 0.907188, 0.548225, 0.013314, 0.006319, 0.001792, 0.000283, 0.000230, 0.000000]  ! 3150 eV
prt%p_reflect(:, 33) = [1.00, 0.987621, 0.958865, 0.901471, 0.461374, 0.011564, 0.006322, 0.001587, 0.000292, 0.000206, 0.000000]  ! 3200 eV
prt%p_reflect(:, 34) = [1.00, 0.987741, 0.958207, 0.894728, 0.372096, 0.010273, 0.006312, 0.001399, 0.000302, 0.000181, 0.000000]  ! 3250 eV
prt%p_reflect(:, 35) = [1.00, 0.987855, 0.957464, 0.886720, 0.292165, 0.009329, 0.006284, 0.001227, 0.000312, 0.000157, 0.000000]  ! 3300 eV
prt%p_reflect(:, 36) = [1.00, 0.987967, 0.956661, 0.877342, 0.227998, 0.008618, 0.006242, 0.001073, 0.000321, 0.000134, 0.000000]  ! 3350 eV
prt%p_reflect(:, 37) = [1.00, 0.988071, 0.955732, 0.865790, 0.177806, 0.008131, 0.006172, 0.000935, 0.000328, 0.000114, 0.000000]  ! 3400 eV
prt%p_reflect(:, 38) = [1.00, 0.988175, 0.954737, 0.851996, 0.140014, 0.007765, 0.006090, 0.000814, 0.000332, 0.000096, 0.000000]  ! 3450 eV
prt%p_reflect(:, 39) = [1.00, 0.988275, 0.953639, 0.835012, 0.111217, 0.007508, 0.005989, 0.000709, 0.000332, 0.000080, 0.000000]  ! 3500 eV
prt%p_reflect(:, 40) = [1.00, 0.988369, 0.952387, 0.813264, 0.089061, 0.007345, 0.005864, 0.000619, 0.000330, 0.000068, 0.000000]  ! 3550 eV
prt%p_reflect(:, 41) = [1.00, 0.988464, 0.951062, 0.786243, 0.072035, 0.007221, 0.005733, 0.000542, 0.000324, 0.000059, 0.000000]  ! 3600 eV
prt%p_reflect(:, 42) = [1.00, 0.988549, 0.949502, 0.749524, 0.058753, 0.007155, 0.005572, 0.000479, 0.000315, 0.000052, 0.000000]  ! 3650 eV
prt%p_reflect(:, 43) = [1.00, 0.988636, 0.947856, 0.702325, 0.048334, 0.007101, 0.005411, 0.000427, 0.000303, 0.000048, 0.000000]  ! 3700 eV
prt%p_reflect(:, 44) = [1.00, 0.988718, 0.945993, 0.638771, 0.040117, 0.007068, 0.005235, 0.000386, 0.000289, 0.000046, 0.000000]  ! 3750 eV
prt%p_reflect(:, 45) = [1.00, 0.988795, 0.943882, 0.556035, 0.033624, 0.007044, 0.005045, 0.000355, 0.000272, 0.000046, 0.000000]  ! 3800 eV
prt%p_reflect(:, 46) = [1.00, 0.988867, 0.941486, 0.459412, 0.028485, 0.007019, 0.004845, 0.000332, 0.000254, 0.000048, 0.000000]  ! 3850 eV
prt%p_reflect(:, 47) = [1.00, 0.988939, 0.938863, 0.366660, 0.024336, 0.006992, 0.004645, 0.000316, 0.000234, 0.000050, 0.000000]  ! 3900 eV
prt%p_reflect(:, 48) = [1.00, 0.989008, 0.935884, 0.288071, 0.021024, 0.006957, 0.004439, 0.000306, 0.000214, 0.000053, 0.000000]  ! 3950 eV
prt%p_reflect(:, 49) = [1.00, 0.989071, 0.932484, 0.226682, 0.018378, 0.006913, 0.004230, 0.000301, 0.000194, 0.000056, 0.000000]  ! 4000 eV
prt%p_reflect(:, 50) = [1.00, 0.989136, 0.928754, 0.180761, 0.016184, 0.006865, 0.004026, 0.000300, 0.000175, 0.000059, 0.000000]  ! 4050 eV
prt%p_reflect(:, 51) = [1.00, 0.989194, 0.924313, 0.145075, 0.014492, 0.006799, 0.003815, 0.000302, 0.000156, 0.000062, 0.000000]  ! 4100 eV
prt%p_reflect(:, 52) = [1.00, 0.989250, 0.919184, 0.117756, 0.013135, 0.006721, 0.003605, 0.000306, 0.000138, 0.000064, 0.000000]  ! 4150 eV
prt%p_reflect(:, 53) = [1.00, 0.989302, 0.913223, 0.096598, 0.012045, 0.006632, 0.003399, 0.000311, 0.000121, 0.000064, 0.000000]  ! 4200 eV
prt%p_reflect(:, 54) = [1.00, 0.989352, 0.906247, 0.080015, 0.011167, 0.006531, 0.003197, 0.000318, 0.000106, 0.000065, 0.000000]  ! 4250 eV
prt%p_reflect(:, 55) = [1.00, 0.989400, 0.898022, 0.066867, 0.010458, 0.006422, 0.003000, 0.000324, 0.000092, 0.000064, 0.000000]  ! 4300 eV
prt%p_reflect(:, 56) = [1.00, 0.989447, 0.888239, 0.056329, 0.009882, 0.006304, 0.002810, 0.000330, 0.000081, 0.000062, 0.000000]  ! 4350 eV
prt%p_reflect(:, 57) = [1.00, 0.989492, 0.876479, 0.047798, 0.009413, 0.006180, 0.002626, 0.000336, 0.000071, 0.000059, 0.000000]  ! 4400 eV
prt%p_reflect(:, 58) = [1.00, 0.989529, 0.861421, 0.040902, 0.009076, 0.006037, 0.002446, 0.000340, 0.000063, 0.000056, 0.000000]  ! 4450 eV
prt%p_reflect(:, 59) = [1.00, 0.989566, 0.842587, 0.035260, 0.008798, 0.005890, 0.002274, 0.000343, 0.000057, 0.000052, 0.000000]  ! 4500 eV
prt%p_reflect(:, 60) = [1.00, 0.989602, 0.818593, 0.030599, 0.008567, 0.005741, 0.002110, 0.000344, 0.000053, 0.000047, 0.000000]  ! 4550 eV
prt%p_reflect(:, 61) = [1.00, 0.989630, 0.785551, 0.026825, 0.008408, 0.005576, 0.001952, 0.000344, 0.000050, 0.000043, 0.000000]  ! 4600 eV
prt%p_reflect(:, 62) = [1.00, 0.989659, 0.740831, 0.023675, 0.008271, 0.005411, 0.001803, 0.000342, 0.000049, 0.000038, 0.000000]  ! 4650 eV
prt%p_reflect(:, 63) = [1.00, 0.989689, 0.679623, 0.021019, 0.008151, 0.005247, 0.001663, 0.000338, 0.000048, 0.000033, 0.000000]  ! 4700 eV
prt%p_reflect(:, 64) = [1.00, 0.989711, 0.593416, 0.018875, 0.008069, 0.005072, 0.001531, 0.000332, 0.000049, 0.000029, 0.000000]  ! 4750 eV
prt%p_reflect(:, 65) = [1.00, 0.989725, 0.486403, 0.017161, 0.008011, 0.004886, 0.001406, 0.000324, 0.000050, 0.000025, 0.000000]  ! 4800 eV
prt%p_reflect(:, 66) = [1.00, 0.989750, 0.390516, 0.015579, 0.007937, 0.004719, 0.001290, 0.000316, 0.000052, 0.000022, 0.000000]  ! 4850 eV
prt%p_reflect(:, 67) = [1.00, 0.989758, 0.306608, 0.014419, 0.007893, 0.004531, 0.001183, 0.000305, 0.000054, 0.000019, 0.000000]  ! 4900 eV
prt%p_reflect(:, 68) = [1.00, 0.989778, 0.247184, 0.013300, 0.007835, 0.004362, 0.001083, 0.000294, 0.000057, 0.000017, 0.000000]  ! 4950 eV
prt%p_reflect(:, 69) = [1.00, 0.989781, 0.199288, 0.012511, 0.007791, 0.004174, 0.000991, 0.000281, 0.000059, 0.000015, 0.000000]  ! 5000 eV

!---------------------------------------------------------------------------------------------

do i = 1, size(surface%table)
  call finalize_reflectivity_table (surface%table(i), .true.)
enddo

end subroutine photon_reflection_std_surface_init

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine finalize_reflectivity_table (table, in_degrees)
!
! Routine to finalize the construction of the reflectivity tables for a surface.
!
! Input:
!   table       -- photon_reflect_table_struct: Surface tables to be finalized.
!   in_degrees  -- logical: Table angles in degrees?
!
! Output:
!   table       -- photon_reflect_table_struct: Finalized surface tables.
!-

Subroutine finalize_reflectivity_table (table, in_degrees)

implicit none

type (photon_reflect_table_struct), target :: table

real(rp) f, deriv, dprob

integer i, j, k
logical in_degrees

! Take the logiarithm of p_reflect. Where zero, just use 10^-20

if (in_degrees) table%angle = table%angle * pi / 180.0_rp

do j = 1, size(table%energy)
  do k = 1, size(table%angle)
    table%p_reflect(k, j) = log(max(1.0e-20_rp, table%p_reflect(k, j)))
  enddo

  ! First interval interpolation: p = c0 + c1 * angle^n

  deriv = (table%p_reflect(3,j) - table%p_reflect(1,j)) / table%angle(3)
  dprob = table%p_reflect(2,j) - table%p_reflect(1,j)
  table%int1(j)%c0     = table%p_reflect(1,j)
  table%int1(j)%n_exp  = deriv * table%angle(2) / dprob
  table%int1(j)%c1     = dprob / table%angle(2)**table%int1(j)%n_exp

enddo

end subroutine finalize_reflectivity_table

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine read_surface_reflection_file (file_name, surface)
!
! Routine to read the reflection probability data for a given type of surface from a file. 
! 
! Input:
!   file_name -- Character(*): Name of the file.
!
! Output:
!   surface   -- photon_reflect_surface_struct: Surface info.
!-

subroutine read_surface_reflection_file (file_name, surface)

implicit none

type (photon_reflect_surface_struct), target :: surface
type (photon_reflect_table_struct), pointer :: prt

character(*) file_name
character(40) name
character(80) description

integer i, j, iu, n_table, it, n_angles, n_energy, ix_row, ios

real(rp) angles(100), energy_min, energy_max, energy_delta, p_reflect(200), energies(200)
real(rp) surface_roughness_rms, roughness_correlation_len

character(*), parameter :: r_name = 'read_surface_reflection_file'

namelist / general / n_table, surface_roughness_rms, roughness_correlation_len, name, description

namelist / table / angles, energy_min, energy_max, energy_delta, energies
namelist / row / ix_row, p_reflect

! Open file

iu = lunget()
open (iu, file = file_name, status = 'old', action = 'read')

! Allocate the reflection tables global variables.

description = ''
name = '???'
surface_roughness_rms = -1
roughness_correlation_len = -1
read (iu, nml = general, iostat = ios)
if (ios /= 0) then
  call out_io (s_fatal$, r_name, 'CANNOT READ "GENERAL" NAMELIST FROM REFLECTIVITY FILE.')
  read (iu, nml = general)  ! Will generate error message
endif

if (name == '???') then
  call out_io (s_fatal$, r_name, 'SURFACE REFLECTIVITY FILES MUST NOW HAVE A "NAME" IN THE "GENERAL" NAMELIST.')
  stop
endif

if (surface_roughness_rms < 0) then
  call out_io (s_warn$, r_name, 'SURFACE_ROUGHNESS_RMS NOT SET. WILL TREAT AS ZERO.')
  surface_roughness_rms = 0
endif

if (roughness_correlation_len < 0) then
  call out_io (s_warn$, r_name, 'ROUGHNESS_CORRELATION_LEN NOT SET. WILL TREAT AS ZERO.')
  roughness_correlation_len = 0
endif

surface%surface_roughness_rms     = surface_roughness_rms
surface%roughness_correlation_len = roughness_correlation_len
surface%name = name
surface%description = description
surface%reflectivity_file = file_name

if (allocated (surface%table)) deallocate (surface%table)
allocate (surface%table(n_table))

! Fill in each table

do it = 1, n_table
  prt => surface%table(it)

  energy_delta = 0
  energies = -1
  angles = -1

  read (iu, nml = table, iostat = ios)
  if (ios /= 0) then
    call out_io (s_fatal$, r_name, 'ERROR READING TABLE NAMELIST \i0\ FROM SURFACE REFLECTIVITY FILE.', it)
    read (iu, nml = table)
  endif

  do i = 1, size(angles)
    if (angles(i) == -1) then
      n_angles = i - 1
      exit
    endif

    if (i == 1) cycle
    if (angles(i) <= angles(i-1)) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "TABLE" NAMELIST NUMBER \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'ANGLEs MUST BE IN INCREASING ORDER.', i_array = [it])
      if (global_com%exit_on_error) stop
    endif
  enddo

  if (angles(1) /= 0 .or. abs(angles(n_angles) - 90) > 1d-10) then
    call out_io (s_fatal$, r_name, &
            'ERROR IN "TABLE" NAMELIST NUMBER \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
            'FIRST ANGLE MUST BE ZERO AND LAST ANGLE MUST BE 90.', i_array = [it])
    if (global_com%exit_on_error) stop
  endif

  if (energy_delta == 0) then
    do i = 1, size(energies)
      if (energies(i) >= 0) cycle
      n_energy = i - 1
      exit
    enddo
  else
    n_energy = 1 + (energy_max - energy_min) / energy_delta
  endif

  allocate(prt%angle(n_angles), prt%energy(n_energy), prt%p_reflect_scratch(n_angles))
  allocate(prt%p_reflect(n_angles, n_energy))
  allocate(prt%int1(n_energy))

  prt%angle = angles(1:n_angles)

  if (energy_delta == 0) then
    prt%energy = energies(1:n_energy)
  else
    prt%energy = [(energy_min + (i-1) * energy_delta, i = 1, n_energy)]
  endif

  prt%max_energy = prt%energy(n_energy)

  do i = 1, n_energy
    p_reflect = -1
    read (iu, nml = row, iostat = ios)

    if (ios /= 0) then
      call out_io (s_fatal$, r_name, &
                'ERROR READING ROW NAMELIST \i0\ FROM SURFACE REFLECTIVITY FILE: ' // file_name, &
                'FOR TABLE: \i0\ ', i_array = [i, ix_row, it])
      read (iu, nml = row, iostat = ios)
    endif

    if (ix_row /= i) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "ROW" NAMELIST IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'ROW MISMATCH: \2i5\ ', &
              'FOR TABLE: \i0\ ', i_array = [i, ix_row, it])
      if (global_com%exit_on_error) stop
    endif

    if (any(p_reflect(1:n_angles) > 1)) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "ROW" NAMELIST \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'A REFLECTION PROBABILITY IS GREATER THAN 1.', &
              'FOR TABLE: \i0\ ', i_array = [it, ix_row])
      if (global_com%exit_on_error) stop
    endif

    if (any(p_reflect(1:n_angles) == -1)) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "ROW" NAMELIST \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'MISSING REFLECTION PROBABILITY.', &
              'FOR TABLE: \i0\ ', i_array = [it, ix_row])
      if (global_com%exit_on_error) stop
    endif
      
    if (any(p_reflect(1:n_angles) < 0)) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "ROW" NAMELIST \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'NEGATIVE REFLECTION PROBABILITY.', &
              'FOR TABLE: \i0\ ', i_array = [it, ix_row])
      if (global_com%exit_on_error) stop
    endif
      
    if (any(p_reflect(n_angles+1:) /= -1)) then
      call out_io (s_fatal$, r_name, &
              'ERROR IN "ROW" NAMELIST \i0\ IN SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'NUMBER OF PROBABILITIES GREATER THAN NUMBER OF ANGLES.', &
              'FOR TABLE: \i0\ ', i_array = [it, ix_row])
      if (global_com%exit_on_error) stop
    endif

    prt%p_reflect(:, i) = p_reflect(1:n_angles)
  enddo

  if (it > 1) then
    if (energy_min > surface%table(it-1)%max_energy) then
      call out_io (s_fatal$, r_name, &
              'ERROR WITH SURFACE REFLECTION PROBABILITY FILE: ' // file_name, &
              'THE MINIMUM ENERGY \f10.2\ OF TABLE #\i5\ IS LARGER THAN THE PREVIOUS TABLE', &
              r_array = [energy_min], i_array = [it])
      if (global_com%exit_on_error) call err_exit
    endif
  endif

enddo

close (iu)

do i = 1, size(surface%table)
  call finalize_reflectivity_table  (surface%table(i), .true.)
enddo

end subroutine read_surface_reflection_file

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine photon_reflectivity (angle, energy, surface, p_reflect, rel_p_specular)
!
! Routine to evaluate the photon reflectivity.
!   probability of absorption          = 1 - p_reflect
!   probability of reflection          = p_reflect
!   probability of specular reflection = p_reflect * rel_p_specular
!   probability of diffuse reflection  = p_reflect * (1 - rel_p_specular)
!
! Use photon_reflection_std_surface_init or read_surface_reflection_file to get surface info.
!
! Input:
!   angle    -- Real(rp): Incident grazing angle in radians.
!   energy   -- Real(rp): Photon energy in eV.
!   surface  -- photon_reflect_surface_struct: surface info
!
! Output:
!   p_reflect      -- Real(rp): Reflection probability. 
!   rel_p_specular -- Real(rp): Relative specular reflection probability.
!-

subroutine photon_reflectivity (angle, energy, surface, p_reflect, rel_p_specular)

implicit none

type (photon_reflect_surface_struct), target :: surface
type (photon_reflect_table_struct), pointer :: prt

real(rp) angle, energy, e_tot, rel_p_specular, p_reflect
real(rp) max_e, f

integer i, j, ie, ix, n_table, n_energy, ixa, ixa0, ixa1, n_ang, n_a

logical ok

character(*), parameter :: r_name = 'photon_reflectivity'

! Singular cases

if (angle < 1d-10) then
  p_reflect  = 1
  rel_p_specular = 1
  return
endif

if (angle > 1.000001 * pi/2 .or. angle < 0) then
  call out_io (s_fatal$, r_name, 'PHOTON_REFLECTIVITY: ANGLE OUT OF RANGE! \f12.7\ ', angle)
  if (global_com%exit_on_error) call err_exit
endif

! If the energy is less than the minimum of the table then just use the minimum.

e_tot = max(surface%table(1)%energy(1), energy)

! If the energy is greater than what the tables covers assume that things scale as energy*angle

n_table = size(surface%table)
max_e = surface%table(n_table)%max_energy
 
if (e_tot > max_e) then
  angle = angle * e_tot / max_e
  e_tot = max_e
  if (angle > pi/2) then
    p_reflect = 0
    rel_p_specular = 0
    return
  endif
endif

! Find which table to use

do i = 1, n_table 
  prt => surface%table(i)
  if (e_tot <= prt%max_energy) exit
enddo

! Interpolation:
! First do simple linear interpolation in energy.
! e_tot is in the energy interval [prt%energy(ie), prt%energy(ie+1)]

do ie = 1, size(prt%energy)-1
  if (e_tot <= prt%energy(ie+1)) exit
enddo

! Find which angle interval incident angle is in.

n_ang = size(prt%angle)
ixa = bracket_index (angle, prt%angle, 1)

f = (e_tot - prt%energy(ie)) / (prt%energy(ie+1) - prt%energy(ie))

! Interpolate

call this_reflect_prob (prt%int1, prt%p_reflect, angle, p_reflect)
rel_p_specular = exp(-(twopi * surface%surface_roughness_rms * 2 * sin(angle) * energy / (c_light * h_planck))**2)

!----------------------------------------------------------------------------
contains

subroutine this_reflect_prob (int1, p_reflect, angle_incident, reflect_prob)

type (interval1_coef_struct), allocatable :: int1(:)
type (spline_struct) ang_spline(6)

real(rp) :: p_reflect(:,:) 
real(rp) reflect_prob, angle_incident
real(rp) c0, c1, n_exp

! If in the first interval then spline interpolation is not good.
! In this case we use the fact that the probability is monotonic and use the form:
!   prob = c0 + c1 * ang^n

if (ixa == 1) then
  c0    = (1 - f) * int1(ie)%c0 + f * int1(ie)%c0
  c1    = (1 - f) * int1(ie)%c1 + f * int1(ie)%c1
  n_exp = (1 - f) * int1(ie)%n_exp + f * int1(ie)%n_exp
  reflect_prob = c0 + c1 * angle_incident**n_exp
  reflect_prob = exp(reflect_prob)
  return
endif

! Now use Akima spline interpolation in angle.
! Only spline the part of reflect_prob that is needed

prt%p_reflect_scratch = (1 - f) * p_reflect(:, ie) + f * p_reflect(:, ie+1)  ! Linear interpolation

ixa0 = max(1, ixa-2)
ixa1 = min(n_ang, ixa+3)

n_a = ixa1 - ixa0 + 1
ang_spline(1:n_a)%x0 = prt%angle(ixa0:ixa1)
ang_spline(1:n_a)%y0 = prt%p_reflect_scratch(ixa0:ixa1)
call spline_akima(ang_spline(1:n_a), ok)
if (.not. ok) call err_exit

call spline_evaluate (ang_spline(1:n_a), angle_incident, ok, reflect_prob)
if (.not. ok) call err_exit

reflect_prob = exp(reflect_prob)

end subroutine this_reflect_prob

end subroutine photon_reflectivity

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine photon_reflection (graze_angle_in, energy, surface, graze_angle_out, phi_out)
!
! Routine to reflect a photon from a surface including both diffuse and specular reflections.
!
! Input:
!   graze_angle_in  -- Real(rp): Incident grazing (not polar) angle in radians.
!   energy          -- Real(rp): Photon energy in eV.
!   surface         -- photon_reflect_surface_struct: surface info
!
! Output:
!   graze_angle_out -- Real(rp): graze_angle in radians. 
!   phi_out         -- Real(rp): Azimuthal angle in radians.
!-

subroutine photon_reflection (graze_angle_in, energy, surface, graze_angle_out, phi_out)

implicit none

type (photon_reflect_surface_struct), target :: surface

real(rp) graze_angle_in, energy, graze_angle_out,  phi_out
real(rp) p_spec, r, lambda

! If the photon energy is lower than 1eV then use 1eV in the calculation.
! Decide if reflection is specular.

lambda = h_planck * c_light / max(1.0_rp, energy)
P_spec = exp(-(fourpi * surface%surface_roughness_rms * sin(graze_angle_in) / lambda)**2)
call ran_uniform(r)

if (r < P_spec) then   ! Is specular
  graze_angle_out = graze_angle_in
  phi_out = 0

else
  call photon_diffuse_scattering (graze_angle_in, energy, surface, graze_angle_out, phi_out)
endif

end subroutine photon_reflection

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine photon_diffuse_scattering (graze_angle_in, energy, surface, graze_angle_out, phi_out, diffuse_param)
!
! Routine to simulate the diffuse scattering of photons. The outgoing angles are
! choosen using the Dugan distribution.
!
! Also see: photon_reflection.
! Use photon_reflection_std_surface_init or read_surface_reflection_file to get surface info.
!
! Input:
!   graze_angle_in  -- Real(rp): Incident grazing (not polar) angle in radians.
!   energy          -- Real(rp): Photon energy in eV.
!   surface         -- photon_reflect_surface_struct: surface info
!
! Output:
!   graze_angle_out -- Real(rp): graze_angle in radians. 
!   phi_out         -- Real(rp): Azimuthal angle in radians.
!   diffuse_param   -- diffuse_param_struct, optional: Internal parameters used in the calculation.
!                        This is used for diagnostics and is not used in standard simulations.
!-

subroutine photon_diffuse_scattering (graze_angle_in, energy, surface, graze_angle_out, phi_out, diffuse_param)

use random_mod
!! use nr, only: chebft, chint, chebev   !! NR has been removed from Bmad Distributions.
use super_recipes_mod, only: super_rtsafe

implicit none

type (photon_reflect_surface_struct), target :: surface
type (diffuse_param_struct) d_param
type (cheb_diffuse_struct) cheb_param
type (diffuse_param_struct), optional :: diffuse_param

real(rp) graze_angle_in, energy, graze_angle_out,  phi_out
real(rp) sigma, t, ctheta2, sign_phi, tot_integral, rel_integral_err, integral, old_integral
real(rp) fl, fh, df, r, p_spec, ran1, ran2, integral_err(0:50)

integer i, ix, j, n_pt, status

character(*), parameter :: r_name = 'photon_diffuse_scattering'

logical ok

! If the photon energy is lower than 1eV then use 1eV in the calculation.

sigma = surface%surface_roughness_rms
T = surface%roughness_correlation_len
d_param%y = sin(graze_angle_in)
d_param%lambda = h_planck * c_light / max(1.0_rp, energy)

! Pick random numbers

call ran_uniform(ran1)

call ran_uniform(ran2)
if (ran2 > 0.5) then
  sign_phi = 1
  ran2 = 2 * ran2 - 1
else
  sign_phi = -1
  ran2 = 2 * ran2 
endif

! Fit the probability distribution in x = sin(graze_angle_out)
! Also compute the coefficients fo the cumulative distribution in x

if (diffuse_com%use_spline_fit) then

  ! First construct the spline fit for the probability function

  n_pt = -1

  call eval_prob_x (d_param%prob_spline(0), 0.0_rp)
  call eval_prob_x (d_param%prob_spline(1), d_param%y)
  call eval_prob_x (d_param%prob_spline(2), 1.0_rp)

  call insert_spline_point (1)
  call insert_spline_point (0)

  do i = 1, 20
    call spline_akima(d_param%prob_spline(0:n_pt), ok)
    tot_integral = 0
    do j = 0, n_pt-1
      ! Use abs so not to be confused when the integral is negative due to fitting errors.
      ! This should not affect the results significantly.
      tot_integral = tot_integral + abs(spline1(d_param%prob_spline(j), d_param%prob_spline(j+1)%x0, -1))  
    enddo

    do j = 1, n_pt-1
      call integral_err_calc(j)
    enddo

    rel_integral_err = sum(integral_err(1:n_pt-1)) / tot_integral

    ix = maxloc(integral_err(1:n_pt-1), 1)

    if (i < 4) cycle
    if (rel_integral_err < diffuse_com%area_err_tol) exit
    if (i == 20) exit
    call insert_spline_point (ix)
    call insert_spline_point (ix-1)
  enddo

  d_param%chx_norm = tot_integral

  ! Integrate the probability function and find where probability integral from 0 to cos(theta2) equals ran1.
  ! First integrate over the spline intervals to bracket cos(theta2) and then use super_rtsafe to interpolate.

  integral = 0
  do j = 0, n_pt-1
    old_integral = integral
    ! Use abs so not to be confused when the integral is negative due to fitting errors.
    integral = integral + abs(spline1(d_param%prob_spline(j), d_param%prob_spline(j+1)%x0, -1))
    if (integral >= ran1 * tot_integral) exit
  enddo

  ctheta2 = super_rtsafe (d_integral, d_param%prob_spline(j)%x0, d_param%prob_spline(j+1)%x0, &
                                      1.0e-12_rp, 1.0e-5_rp, status)

! Fit the probability distribution to Chebyshev polynomials.
! This is known to produce bad results for smoother surfaces so eventually this
! will be eliminated. Keep for now for cross-check purposes.
! NOTE: Due to the removal of NR, the calls to NR routines have been commented out 2021/12. Look for "!!NR"
else
!!NR  cheb_param%cch = chebft(0.0_rp, 1.0_rp, n_cheb_term$, prob_x_diffuse_vec)
!!NR  cheb_param%cch_int = chint (0.0_rp, 1.0_rp, cheb_param%cch)
!!NR  !  evaluate the normalization constant
!!NR  d_param%chx_norm = chebev(0.0D0, 1.0D0, cheb_param%cch_int, 1.0D0)
!!NR  ! find the value of x for which the cumulative probability equals the random number
!!NR  ctheta2 = super_rtsafe(cumulx, 0.0E0_rp, 1.0E0_rp, 1.0e-12_rp, 1.0D-5, status)
endif

! Evaluate the normalization constant for the cumulative probability in phi, for this x

graze_angle_out = asin(ctheta2)
d_param%x = ctheta2
d_param%c_norm = cos_phi(sigma, T, pi, d_param)

! find the value of phi for which the cumulative probability equals ran2

call cumulr(0.0_rp, fl, df, status)
call cumulr(pi, fh, df, status)
if ((fl > 0 .and. fh > 0) .or. (fl < 0 .and. fh < 0)) then 
  call out_io (s_fatal$, r_name, 'ROOT NOT BRACKETED FOR PHI CALC!', 'fl, fh: \2es14.5\ ', r_array = [fl, fh])
  call output_specular_reflection_input_params(d_param, surface)
endif

phi_out = sign_phi * super_rtsafe(cumulr, 0.0E0_rp, pi, 1.0e-12_rp, 1.0e-5_rp, status)

if (present(diffuse_param)) then
  diffuse_param = d_param
  diffuse_param%n_pt_spline = n_pt
endif

!---------------------
contains

subroutine eval_prob_x(spline, x)
type (spline_struct) spline
real(rp) x
spline%x0  = x
spline%y0 = prob_x_diffuse(spline%x0, d_param, surface)
n_pt = n_pt + 1
end subroutine eval_prob_x

!---------------------
! contains

subroutine insert_spline_point(ix)
integer ix
d_param%prob_spline(ix+2:n_pt+1) = d_param%prob_spline(ix+1:n_pt)
call eval_prob_x(d_param%prob_spline(ix+1), (d_param%prob_spline(ix)%x0 + d_param%prob_spline(ix+2)%x0) / 2)
end subroutine insert_spline_point

!---------------------
! contains

subroutine integral_err_calc(ix)
integer ix
real(rp) x1, y1, x2, y2, a, b, a_spline, a_parab

! The error in the integral when using the spline fit is estimated by the difference in the
! integral using the spline fit and the integral when using a parabolic fit.

if (ix == 0 .or. ix == n_pt) return

x1 = d_param%prob_spline(ix)%x0 - d_param%prob_spline(ix-1)%x0
y1 = d_param%prob_spline(ix)%y0 - d_param%prob_spline(ix-1)%y0

x2 = d_param%prob_spline(ix+1)%x0 - d_param%prob_spline(ix-1)%x0
y2 = d_param%prob_spline(ix+1)%y0 - d_param%prob_spline(ix-1)%y0

a = (y2 * x1 - y1 * x2) / (x1 * x2**2 - x2 * x1**2)
b = (y2 * x1**2 - y1 * x2**2) / (x1**2 * x2 - x2**2 * x1)

a_parab = a * x2**3 / 3 + b * x2**2 / 2 + d_param%prob_spline(ix-1)%y0 * (d_param%prob_spline(ix+1)%x0 - d_param%prob_spline(ix-1)%x0)
a_spline = spline1(d_param%prob_spline(ix-1), d_param%prob_spline(ix)%x0, -1) + spline1(d_param%prob_spline(ix), d_param%prob_spline(ix+1)%x0, -1)

integral_err(ix) = abs(a_parab - a_spline)

end subroutine integral_err_calc

!---------------------
! contains

!+
! Subroutine d_integral (x, fn, df, status)
!
! Wrapper function passed to super_rtsafe.
! Contained routine to calculate integrated probability distribution in x = sin(graze_angle_out).
!-

subroutine d_integral (x, fn, df, status)

implicit none

real(rp), intent(in) :: x
real(rp), intent(out) :: fn, df
integer status

!

fn = (old_integral + abs(spline1(d_param%prob_spline(j), x, -1))) / tot_integral - ran1
df = abs(spline1(d_param%prob_spline(j), x, 0)) / tot_integral

end subroutine d_integral

!---------------------------------------------------------------------------------------------
! contains

!+
! Subroutine cumulr (phi, fn, df, status)
!
! Wrapper function passed to super_rtsafe.
! Contained routine to calculate integrated probability distribution in x = sin(graze_angle_out).
!-

subroutine cumulr (phi, fn, df, status)

implicit none

real(rp), intent(in) :: phi 
real(rp), intent(out) :: fn, df
real(rp) sigma, T
integer status

!

sigma = surface%surface_roughness_rms
T = surface%roughness_correlation_len

fn = cos_phi(sigma, T, phi, d_param) / d_param%c_norm - ran2
df = ptwo(sigma, T, phi, d_param) / d_param%c_norm

end subroutine cumulr

!---------------------------------------------------------------------------------------------
! contains

!+
! Subroutine cumulx (x, fn, df, status)
!
! Wrapper function passed to rtsafe.
! Contained routine to calculate integrated probability distribution in x = sin(graze_angle_out).
!-

subroutine cumulx (x, fn, df, status)

implicit none

real(rp), intent(in) :: x
real(rp), intent(out) :: fn, df
integer status

!

fn = 0; df = 0  ! This to keep the compiler from complaining

!!NR  fn = chebev(0.0E0_rp, 1.0E0_rp, cheb_param%cch_int, x) / d_param%chx_norm - ran1
!!NR  df = chebev(0.0E0_rp, 1.0E0_rp, cheb_param%cch, x) / d_param%chx_norm

end subroutine cumulx

!---------------------------------------------------------------------------------------------
!contains

!+
! Function prob_x_diffuse_vec (x) result (prob_x)
!
! Contained routine to calculate integrated probability distribution in x = sin(graze_angle_out).
! 
! Input:
!   x(:)    -- Real(rp): sin(graze_angle_out) array
!
! Output:
!   prob(:) -- Real(rp): Integrated probability array.
!-

function prob_x_diffuse_vec (x) result (prob_x)

implicit none

real(rp), intent(in) :: x(:)
real(rp) prob_x(size(x))

integer itt

!

do itt = 1, size(x)
  prob_x(itt) = prob_x_diffuse(x(itt), d_param, surface)
enddo

end function prob_x_diffuse_vec

end subroutine photon_diffuse_scattering

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function prob_x_diffuse (x, d_param, surface) result (prob_x)
!
! Contained routine to calculate integrated probability distribution in x = sin(graze_angle_out).
! 
! Input:
!   x    -- Real(rp): sin(graze_angle_out)
!
! Output:
!   prob -- Real(rp): Integrated probability.
!-

function prob_x_diffuse (x, d_param, surface) result (prob_x)

implicit none

type (diffuse_param_struct) d_param
type (photon_reflect_surface_struct), target :: surface

real(rp), intent(in) :: x
real(rp) prob_x
real(rp) pxa, xpy, k, r, sigma, T, xysq, tau, s, g, h, a, y, lambda
real(rp) qexp, bs, fz, b, g0, fexp, factor, pxs, bi, term, ri, xlogi
real(rp) fexpi, fzi, xlogg, a_minus_1

integer i, itt

logical :: xyzero

character(*), parameter :: r_name = 'prob_x_diffuse'
!

y = d_param%y
sigma = surface%surface_roughness_rms
T = surface%roughness_correlation_len
lambda = d_param%lambda

!

pxa = 0.0
xpy = x + y
xysq = x**2 + y**2
k = twopi/lambda
r = sigma*k
tau = T/sigma
s = k*T
g = (r*xpy)**2
h = sqrt((1-x**2)*(1-y**2))
a = h/(1 + x*y)
a_minus_1 = -0.5_rp * xpy*xpy  ! Valid when x and y small
xyzero = .false.

if (x == 1.0 .or. y == 1.0) xyzero = .true.

if (g < gmin) then
  qexp = (2-xysq-2*h)*s**2/4
  bs = s**2*h/2
  fz = zzfi(a, a_minus_1, bs, xyzero)
  prob_x = 0.5*(1 + x*y)**2*r**2*s**2*exp(-qexp-g)*fz
  return
end if

b = h*tau**2/2/xpy**2

if (g > gmax .and. appsw) then
  g0 = tau**2/2/xpy**4*(1 + x*y)**2
  fexp = -(2-xysq-2*h)*tau**2/4/xpy**2
  fz = zzfi(a, a_minus_1, b, xyzero)
  pxa = g0*exp(fexp)*fz
  prob_x = pxa
  return
endif

!

g0 = tau**2*g/2/xpy**4*(1 + x*y)**2
factor = 1.0
pxs = 0.0

do i = 1, maxsum
  bi = b*g/i
  fexpi = -g-(2-xysq-2*h)*tau**2*g/4/i/xpy**2
  fzi = zzfi(a, a_minus_1, bi, xyzero)

  if (i < ismax) then
    factor = factor*i
    term = exp(fexpi)*g**i/factor/i
  else
    ri = i
    xlogi = ri*log(ri)-ri + 0.5*log(twopi*ri)
    xlogg = ri*log(g)-xlogi + fexpi
    term = exp(xlogg)/i
  end if

  term = term*fzi
  pxs = pxs + term

  if (term <= converge * pxs) exit
enddo

prob_x = pxs*g0

if (i > maxsum) then
  call out_io (s_error$, r_name, 'Convergence failed!', &
        'pxs, g0: \2es14.5\ ', r_array = [pxs, g0])
  call output_specular_reflection_input_params(d_param, surface)
endif

end function prob_x_diffuse

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! function ptwo(sigma, T, phi, d_param) result (p_two)
!
! unnormalized two-dimensional probability distribution in x and phi
! polar angles relative to surface normal
! azimuthal angle relative to plane of incidence (plane of incoming ray and surface normal)
! 1/y suppressed
!
! Private routine.
!-

function ptwo(sigma, T, phi, d_param) result (p_two)

implicit none

type (diffuse_param_struct) d_param

real(rp) p_two, sigma, t, lambda, y, x, phi, k, factor, p_twoa, r, tau, s, cp, xpy, a_minus_1
real(rp) xysq, g, h, a, qexp, b, g0, pxa, pxs, xmmax, xtest, fexp, term, ri, xlogi, xlogg
real(rp) :: qlim1 = 2.0D3, qlim2 = -2.0D2

integer i

character(*), parameter :: r_name = 'ptwo'

!

lambda = d_param%lambda
y      = d_param%y
x      = d_param%x

p_twoa = 0.0
p_two = 0.0
k = twopi/lambda
r = sigma*k
tau = T/sigma
s = k*T
cp = cos(phi)
xpy = x+y
xysq = x**2+y**2
g = (r*xpy)**2
h = sqrt((1-x**2)*(1-y**2))
a = h/(1+x*y)
a_minus_1 = -0.5_rp * xpy*xpy  ! Valid when x and y small

if (g < gmin) then
  qexp = (2-xysq-2*h*cp)*s**2/4
  p_two = 0.5/twopi*(1+x*y)**2*(1-a*cp)**2*r**2*s**2*exp(-qexp-g)
  return
end if

b = h*tau**2/2/xpy**2
qexp = (2-xysq)*tau**2/4/xpy**2-b*cp
if(g > gmax.and.appsw) then
  g0 = tau**2/2/twopi/xpy**4*(1+x*y)**2
  pxa = g0*exp(-(2-xysq)*tau**2/4/xpy**2+b*cp)*(1-a*cp)**2
  p_two = pxa
else
  if (g*qexp > qlim1) return
  xmmax = zmmax(qexp,g)
  xtest = xmmax*log(g)-(xmmax*log(xmmax)-xmmax+0.5*log(twopi*xmmax))-g*qexp/xmmax
  if (xtest < qlim2) return
  g0 = tau**2*g/2/twopi/xpy**4*(1+x*y)**2*(1-a*cp)**2
  factor = 1.0
  pxs = 0.0
  i = 0

  do i = 1,maxsum
    fexp = -g*(1.0+qexp/i)
    if (i < ismax) then
      factor = factor*i
      term = exp(fexp)*g**i/factor/i
    else
      ri = i
      xlogi = ri*log(ri)-ri+0.5*log(twopi*ri)
      xlogg = ri*log(g)-xlogi+fexp
      term = exp(xlogg)/i
    end if
    pxs = pxs+term
    if(term < pxs * converge) exit
  enddo

  p_two = pxs*g0

  if (i > maxsum) then
    call out_io (s_error$, r_name, 'Convergence failed!', &
            'phi, pxs, g0: \3es14.5\ ', r_array = [phi, pxs, g0])
    call output_specular_reflection_input_params(d_param)
  endif

end if

end function ptwo

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zmmax(x,g)
!
! Private routine.
!-

function zmmax(x,g) result (zm_max)

implicit none

real(rp) zm_max, x, g
real(rp) :: fk1(3) = [-0.872884, -0.00154422, 3.63838E-7]
real(rp) :: fk2(3) = [6.01394, 0.01683, -3.64794E-6]
real(rp) :: fk3(3) = [0.843621, 0.0006004, -1.33207E-7]
real(rp) :: fk4(3) = [-0.000901106, -1.2543E-6, 0.0]

integer ik

!

zm_max = 0.0
do ik = 1, 3
  zm_max = zm_max + fk1(ik)*x**(ik-1) + fk2(ik)*x**(ik-1)*sqrt(g) + fk3(ik)*x**(ik-1)*g + fk4(ik)*x**(ik-1)*g**2
enddo

if (zm_max < 1) zm_max = 1

end function zmmax

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function cos_phi (sigma, t, phi, d_param) result (cphi)
!
!  computes  unnormalized cumulative distribution function in phi for a given x
!  polar angles relative to surface normal
!  azimuthal angle relative to plane of incidence (plane of incoming ray and surface normal)
!  1/y suppressed
!
! Private routine to calculate integrated probability distribution in x = sin(graze_angle_out).
!-

function cos_phi (sigma, T, phi, d_param) result (cphi)

implicit none

type (diffuse_param_struct) d_param

real(rp) cphi, sigma, t, x, y, phi, a_minus_1
real(rp) k, lambda, cphia, xpy, xysq, r, tau, s, g, h, cphis, bi, fexpi
real(rp) factor, a, fexp, qexp, bs, fz, b, g0, fzi, term, ri, xlogi, xlogg

integer i

logical xyzero

character(*), parameter :: r_name = 'cos_phi'

!

lambda = d_param%lambda
y      = d_param%y
x      = d_param%x

cphia = 0.0
cphi = 0.0
if (phi == 0.0) return
xpy = x+y
xysq = x**2+y**2
k = twopi/lambda
r = sigma*k
tau = T/sigma
s = k*T
g = (r*xpy)**2
h = sqrt((1-x**2)*(1-y**2))
a = h/(1+x*y)
a_minus_1 = -0.5_rp * xpy*xpy  ! Valid when x and y small
xyzero = .false.

if (x == 1.0 .or. y == 1.0) xyzero = .true.

if (g < gmin) then
  qexp = (2-xysq-2*h)*s**2/4
  bs = s**2*h/2
  fz = zzfp(a, a_minus_1, bs, phi, xyzero, d_param)
  cphi = 0.5/twopi*(1+x*y)**2*r**2*s**2*exp(-qexp-g)*fz
  return
end if

b = h*tau**2/2/xpy**2
qexp = (2-xysq)*tau**2/4/xpy**2-b*cos(phi)

if (g > gmax .and. appsw) then
  g0 = tau**2/2/twopi/xpy**4*(1+x*y)**2
  fexp = -(2-xysq-2*h)*tau**2/4/xpy**2
  fz = zzfp(a, a_minus_1, b, phi, xyzero, d_param)
  cphia = g0*exp(fexp)*fz
  cphi = cphia
else
  g0 = tau**2*g/2/twopi/xpy**4*(1+x*y)**2
  factor = 1.0
  cphis = 0.0

  do i = 1, maxsum
    bi = b*g/i
    fexpi = -g-(2-xysq-2*h)*tau**2*g/4/i/xpy**2
    fzi = zzfp(a, a_minus_1, bi, phi, xyzero, d_param)
    if (i < 20) then
      factor = factor*i
      term = exp(fexpi)*g**i/factor/i
    else
      ri = i
      xlogi = ri*log(ri)-ri+0.5*log(twopi*ri)
      xlogg = ri*log(g)-xlogi+fexpi
      term = exp(xlogg)/i
    end if
    term = term*fzi
    cphis = cphis+term
    if (term < converge * cphis) exit
  enddo

  cphi = cphis*g0

  if (i > maxsum) then
    call out_io (s_error$, r_name, 'Convergence failed.', &
              'phi, cphis, g0: \3es14.5\ ', r_array = [phi, cphis, g0])
    call output_specular_reflection_input_params(d_param)
  endif

end if

end function cos_phi

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zzfi (a, a_minus_1, b, xyzero) result (fi)
!
! Private function used in calculating diffuse scattering distribution.
!-

function zzfi (a, a_minus_1, b, xyzero) result (fi)

implicit none

real(rp) a, a_minus_1, b, fi
logical xyzero

!

if (xyzero) then
  fi = (1+a**2)
elseif (b < bmax) then
  fi = ((1+a**2)*zbessi0(b)-a*(a+2*b)*zbessi1(b)/b)/sqrt(twopi*b)
else
  fi = 2*hzz(a, a_minus_1, b, twopi)/twopi
end if

end function zzfi

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zzfp (a, a_minus_1, b, phi, xyzero, d_param) result (fp)
!
! Private function used in calculating diffuse scattering distribution.
!-

function zzfp (a, a_minus_1, b, phi, xyzero, d_param) result (fp)

implicit none

type (diffuse_param_struct) :: d_param

real(rp) a, a_minus_1, b, phi, fp
real(rp) sp, sp2
logical xyzero

!

if(xyzero) then
  sp = sin(phi)
  sp2 = sin(2*phi)
  fp = ((1+a**2/2)*phi-2*a*sp+a**2/4*sp2)
  return
end if
if(b < bmax) then
  fp = zzexp(a,b,phi, d_param)
else
  fp = hzz(a,a_minus_1,b, phi)
end if

end function zzfp
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function hzz(a, a_minus_1, b, phi) result (h)
!
! Private function used in calculating diffuse scattering distribution.
!-

function hzz(a, a_minus_1, b, phi) result (h)

implicit none

real(rp) a, a_minus_1, b, phi, h
real(rp) arg, hzz1, h2zz
real(rp), parameter :: pi_root = sqrt(pi/2)

! This routine computes the phi integral used in the cumulative distribution function. 
! Good for all b, but makes a small phi approximation. Use for b>100.

arg = phi*sqrt(b/2)
hzz1 = -a*phi*exp(-arg**2)*(4*b+a*(3+b*(phi**2-4)))
hzz1 = hzz1/4/b**2
!!h2zz = 4*a*(1-2*b)*b + 4*b**2 + a**2*(3-4*b+4*b**2)
if (abs(a_minus_1) < 1d-6) then
  h2zz = 4 * (a_minus_1 * b)**2 - 4 * a_minus_1 * b + 3*a*a
else
  h2zz =4 * ((a - 1) * b)**2 - 4 * (a - 1) * b + 3*a*a
endif
h2zz = h2zz/4/b**2.5
h2zz = h2zz*pi_root*erf(arg)
h = hzz1+h2zz

end function hzz

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zbessi(n, x) result (zbes)

!
! From numerical recipes: bessel function of index n, multiplied by sqrt(2*Pi*x)*exp(-x).
! Changed IACC from 40 to 400 to improve accuracy for large x.
!
! Private function used in calculating diffuse scattering distribution.
!-

function zbessi(n, x) result (zbes)

implicit none

real(8) zbes, x
real(rp), parameter :: bigno = 1d30, bigni = 1d-30
real(8) bi,bim,bip,tox

integer j, m, n
integer, parameter :: iacc = 400

!  uses zbessi0

if (n < 2) call err_exit ! bad argument n in zbessi

if (x == 0.0) then
  zbes = 0.
else
  tox = 2.0/abs(x)
  bip = 0.0
  bi = 1.0
  zbes = 0.
  m = 2*((n+int(sqrt(float(IACC*n)))))

  do j = m,1,-1
    bim = bip+float(j)*tox*bi
    bip = bi
    bi = bim
    if (abs(bi) > BIGNO) then
      zbes = zbes*BIGNI
      bi = bi*BIGNI
      bip = bip*BIGNI
    endif
    if (j == n) zbes = bip
  enddo
  zbes = zbes*zbessi0(x)/bi
  if (x < 0.0 .and. mod(n,2) == 1) zbes = -zbes
endif

end function zbessi
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zbessi1 (x) result (zbes)
!
! This is the I1 Bessel function from Numerical Recipes multiplied by sqrt(2*Pi*x)*exp(-x)
! Private function used in calculating diffuse scattering distribution.
!-

function zbessi1(x) result (zbes)

implicit none
real(8) zbes, x
real(8) ax, sco, y, sc0
real(8), parameter :: p1 = 0.5d0, p2 = 0.87890594d0, p3 = 0.51498869d0, p4 = 0.15084934d0, &
        p5 = 0.2658733d-1, p6 = 0.301532d-2, p7 = 0.32411d-3
real(8), parameter :: q1 = 0.39894228d0, q2 = -0.3988024d-1, q3 = -0.362018d-2, q4 = 0.163801d-2, &
   q5 = -0.1031555d-1, q6 = 0.2282967d-1, q7 = -0.2895312d-1, q8 = 0.1787654d-1, q9 = -0.420059d-2

!

ax = abs(x)
if (abs(x) < 3.75) then
  y = (x/3.75)**2
  sc0 = sqrt(twopi*ax)*exp(-ax)
  zbes = sc0*(x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))))
else
  y = 3.75/ax
  sc0 = sqrt(twopi)
  zbes = sc0*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
  if (x < 0.0) zbes = -zbes
endif

end function zbessi1

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zbessi0 (x) result (zbes)
!
! This is the I0 Bessel function from Numerical Recipes multiplied by sqrt(2*Pi*x)*exp(-x)
! Private function used in calculating diffuse scattering distribution.
!-

function zbessi0 (x) result (zbes)

implicit none

real(rp) zbes, x
real(8) ax, y, sc0
real(8), parameter :: p1 = 1.0d0, p2 = 3.5156229d0, p3 = 3.0899424d0, p4 = 1.2067492d0, &
                 p5 = 0.2659732d0, p6 = 0.360768d-1, p7 = 0.45813d-2
real(8), parameter :: q1 = 0.39894228d0, q2 = 0.1328592d-1, q3 = 0.225319d-2, q4 = -0.157565d-2, &
        q5 = 0.916281d-2, q6 = -0.2057706d-1, q7 = 0.2635537d-1, q8 = -0.1647633d-1, q9 = 0.392377d-2

!

ax = abs(x)
if (ax < 3.75) then
  y = (x/3.75)**2
  sc0 = sqrt(twopi*ax)*exp(-ax)
  zbes = (p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))*sc0
else
  y = 3.75/ax
  sc0 = sqrt(twopi)
  zbes = sc0*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
endif

end function zbessi0

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Function zzexp (a, b, phi, d_param) result (zexp)
!
! This routine computes the phi integral used in the cumulative distribution function.
! Series does not converge well for b>100.
! Private function used in calculating diffuse scattering distribution.
!-

function zzexp (a, b, phi, d_param) result (zexp)

implicit none 

type (diffuse_param_struct) :: d_param

real(rp) a, b, phi, zexp
real(rp) sp, cp, sp2, cp2, sp3, sp4, r0, r1, r2, spk, cpk, rk, zterm

integer k

character(*), parameter :: r_name = 'zzexp'

!

zexp = 0.0
if (phi == 0.0) return

sp = sin(phi)
cp = cos(phi)
sp2 = sin(2*phi)
cp2 = cos(2*phi)
sp3 = sin(3*phi)
sp4 = sin(4*phi)
r0 = (1+a**2/2)*phi-2*a*sp+a**2/4*sp2
r1 = -a*phi+(1+3*a**2/4)*sp-a*cp*sp+a**2/12*sp3
r2 = a**2/4*phi-a*sp+(12*(a**2+2)*sp2-16*a*sp3+3*a**2*sp4)/48
zexp = r0*zbessi0(b)+2*r1*zbessi1(b)+2*r2*zbessi(2,b)

do k = 3, maxsum
 spk = sin(k*phi)
 cpk = cos(k*phi)
 rk = 2*a*cpk*sp/(k**2-1.0)-2*a*cp*(a*cpk*sp/(k**2-4.0)+k*spk/(k**2-1.0))+ &
         spk*((2+a**2)*(k**2-4.0)+a**2*k**2*cp2)/2/k/(k**2-4.0)
 zterm = 2*rk*zbessi(k,b)
 zexp = zexp+zterm
 if (abs(zterm) < converge * abs(zexp)) exit
enddo

zexp = zexp/sqrt(twopi*b)

if (k > maxsum) then
  call out_io (s_error$, r_name, 'Convergence failed!', &
          'a, b, phi, zzexp: \4es14.4\ ', r_array = [a, b, phi, zexp])
  call output_specular_reflection_input_params(d_param)
endif

end function zzexp

!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!+
! Subroutine output_specular_reflection_input_params(d_param, surface)
!
! Private function used in calculating diffuse scattering distribution.
!-

subroutine output_specular_reflection_input_params (d_param, surface)

implicit none

type (diffuse_param_struct) :: d_param
type (photon_reflect_surface_struct), optional :: surface

real(rp) sigma, t

character(*), parameter :: r_name = 'output_specular_reflection_input_params'

!

sigma = -1
T = -1

if (present(surface)) then
  sigma = surface%surface_roughness_rms
  T     = surface%roughness_correlation_len
endif

call out_io (s_blank$, r_name, 'sigma, T, lambda, x, y: \4es14.5\ ', &
          r_array = [sigma, T, d_param%lambda, d_param%x, d_param%y])

end subroutine output_specular_reflection_input_params

end module