#ifndef lint
static const char RCSid[] = "$Id: atmos.c,v 2.2 2024/07/19 23:38:28 greg Exp $";
#endif
/*
The Radiance Software License, Version 2.0

Radiance v5.4 Copyright (c) 1990 to 2022, The Regents of the University of
California, through Lawrence Berkeley National Laboratory (subject to receipt
of any required approvals from the U.S. Dept. of Energy).  All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

(1) Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

(2) Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

(3) Neither the name of the University of California, Lawrence Berkeley
National Laboratory, U.S. Dept. of Energy nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

You are under no obligation whatsoever to provide any bug fixes, patches,
or upgrades to the features, functionality or performance of the source
code ("Enhancements") to anyone; however, if you choose to make your
Enhancements available either publicly, or directly to Lawrence Berkeley
National Laboratory, without imposing a separate written license agreement
for such Enhancements, then you hereby grant the following license: a
non-exclusive, royalty-free perpetual license to install, use, modify,
prepare derivative works, incorporate into other computer software,
distribute, and sublicense such enhancements or derivative works thereof,
in binary and source code form.

=================================================================================

Copyright (c) 2017 Eric Bruneton
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors
may be used to endorse or promote products derived from this software without
specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Additional modifications by T. Wang
- N-dimensional data storage and lookup using Radiance data structure
- Mie scattering coefficients precomputed using libRadTrans and interpolated at runtime
- Translation to C from C++ with multi-threading
*/

#include "atmos.h"
#include "data.h"

#if defined(_WIN32) || defined(_WIN64)
#include <windows.h>
#else
#include <pthread.h>
#endif

#if defined(_WIN32) || defined(_WIN64)
#define THREAD HANDLE
#define CREATE_THREAD(thread, func, context)                                   \
  ((*(thread) = CreateThread(NULL, 0, (func),          \
                             (context), 0, NULL)) != NULL)
#define THREAD_RETURN DWORD WINAPI
#define THREAD_JOIN(thread) WaitForSingleObject(thread, INFINITE)
#define THREAD_CLOSE(thread) CloseHandle(thread)
#else
#define THREAD pthread_t
#define CREATE_THREAD(thread, func, context)                                   \
  (pthread_create((thread), NULL, (func), (context)) == 0)
#define THREAD_RETURN void *
#define THREAD_JOIN(thread) pthread_join(thread, NULL)
#endif

typedef struct {
  int start_l;
  int end_l;
  const Atmosphere *atmos;
  DATARRAY *tau_dp;
  DATARRAY *delta_rayleigh_scattering_dp;
  DATARRAY *delta_mie_scattering_dp;
  DATARRAY *scattering_dp;
} Scat1Tdat;

typedef struct {
  int start_l;
  int end_l;
  DATARRAY *tau_dp;
  DATARRAY *delta_scattering_density_dp;
  DATARRAY *delta_multiple_scattering_dp;
  DATARRAY *scattering_dp;
} ScatNTdat;

typedef struct {
  int start_l;
  int end_l;
  const Atmosphere *atmos;
  DATARRAY *tau_dp;
  DATARRAY *delta_rayleigh_scattering_dp;
  DATARRAY *delta_mie_scattering_dp;
  DATARRAY *delta_multiple_scattering_dp;
  DATARRAY *delta_irradiance_dp;
  int scattering_order;
  DATARRAY *delta_scattering_density_dp;
} ScatDenseTdat;

const double ER = 6360000.0; /* Earth radius in meters */
const double AR = 6420000; /* Atmosphere radius in meters */
const double AH = 60000;   /* Atmosphere thickness in meters */

const float START_WVL = 390.0;
const float END_WVL = 770.0;
const double SUNRAD = 0.004625;

/* The cosine of the maximum Sun zenith angle */
const double MU_S_MIN = -0.2;

/* Scale heights (m) */
/* Rayleigh scattering */
const double HR_MS = 8820; /* Midlatitude summer */
const double HR_MW = 8100; /* Midlatitude winter */
const double HR_SS = 8550; /* Subarctic summer */
const double HR_SW = 7600; /* Subarctic winter */
const double HR_T = 9050;  /* Tropical */

const double AOD0_CA = 0.115; /* Broadband AOD for continental average */

const double SOLOMG = 6.7967e-5; /* .533 apex angle */
const int WVLSPAN = 400;         /* 380nm to 780nm */

/* Aerosol optical depth at 550nm for continental clean */
const double AOD_CC_550 = 0.05352;

/** 380nm to 780nm at 20nm intervals **/
/* Extraterrestrial solar W/m^2/nm */
const float EXTSOL[NSSAMP] = {1.11099345, 1.75363199, 1.67126921, 1.97235667,
                              2.01863015, 1.93612482, 1.87310758, 1.88792588,
                              1.86274213, 1.84433248, 1.80370942, 1.73020456,
                              1.67036654, 1.57482149, 1.53646383, 1.45515319,
                              1.38207435, 1.3265141,  1.26648111, 1.20592043};

/* Rayleight scattering coefficients at sea level in m^-1 */
const float BR0_MS[NSSAMP] = {
    4.63889967e-05, 3.76703293e-05, 3.09193474e-05, 2.56223081e-05,
    2.14198405e-05, 1.80483398e-05, 1.53177709e-05, 1.30867298e-05,
    1.12487269e-05, 9.72370830e-06, 8.44932642e-06, 7.37769845e-06,
    6.47116653e-06, 5.70005327e-06, 5.04091603e-06, 4.47430240e-06,
    3.98549839e-06, 3.56178512e-06, 3.19293761e-06, 2.87072599e-06};

const float BR0_MW[NSSAMP] = {
    5.03857564e-05, 4.09159104e-05, 3.35832809e-05, 2.78298619e-05,
    2.32653203e-05, 1.96033395e-05, 1.66375116e-05, 1.42142496e-05,
    1.22178890e-05, 1.05614786e-05, 9.17729921e-06, 8.01334246e-06,
    7.02870602e-06, 6.19115557e-06, 5.47522872e-06, 4.85979708e-06,
    4.32887894e-06, 3.86865960e-06, 3.46803311e-06, 3.11806054e-06};

const float BR0_SS[NSSAMP] = {
    4.73789047e-05, 3.84741872e-05, 3.15791442e-05, 2.61690698e-05,
    2.18769246e-05, 1.84334785e-05, 1.56446412e-05, 1.33659912e-05,
    1.14887667e-05, 9.93120528e-06, 8.62962901e-06, 7.53513326e-06,
    6.60925660e-06, 5.82168833e-06, 5.14848557e-06, 4.56978082e-06,
    4.07054608e-06, 3.63779107e-06, 3.26107261e-06, 2.93198523e-06};

const float BR0_SW[NSSAMP] = {
    5.30623659e-05, 4.30894595e-05, 3.53673035e-05, 2.93082494e-05,
    2.45012287e-05, 2.06447149e-05, 1.75213353e-05, 1.49693438e-05,
    1.28669320e-05, 1.11225291e-05, 9.66481886e-06, 8.43902999e-06,
    7.40208736e-06, 6.52004427e-06, 5.76608570e-06, 5.11796089e-06,
    4.55883914e-06, 4.07417187e-06, 3.65226316e-06, 3.28369923e-06};
const float BR0_T[NSSAMP] = {
    4.55376661e-05, 3.69790036e-05, 3.03519157e-05, 2.51520876e-05,
    2.10267438e-05, 1.77171168e-05, 1.50366593e-05, 1.28465622e-05,
    1.10422904e-05, 9.54525887e-06, 8.29426444e-06, 7.24230298e-06,
    6.35240772e-06, 5.59544593e-06, 4.94840517e-06, 4.39219003e-06,
    3.91235655e-06, 3.49641927e-06, 3.13434084e-06, 2.81804244e-06};

/* Ozone extinction coefficients at sea level */
const float AO0[NSSAMP] = {
    3.5661864e-09, 1.4848362e-08, 4.5415674e-08, 1.2446184e-07, 2.6461576e-07,
    4.8451984e-07, 8.609148e-07,  1.480537e-06,  1.8809e-06,    2.5107328e-06,
    2.5263174e-06, 2.313507e-06,  1.727741e-06,  1.2027012e-06, 7.915902e-07,
    5.0639202e-07, 3.5285684e-07, 2.23021e-07,   1.7395638e-07, 1.5052574e-07};

const double MIE_G = 0.7;  /* Mie eccentricity */

const int TRANSMITTANCE_TEXTURE_WIDTH = 256;
const int TRANSMITTANCE_TEXTURE_HEIGHT = 64;

/* Resolution halved by original impl */
const int SCATTERING_TEXTURE_R_SIZE = 16;
const int SCATTERING_TEXTURE_MU_SIZE = 64;
const int SCATTERING_TEXTURE_MU_S_SIZE = 16;
const int SCATTERING_TEXTURE_NU_SIZE = 4;

const int IRRADIANCE_TEXTURE_WIDTH = 64;
const int IRRADIANCE_TEXTURE_HEIGHT = 16;

static inline double clamp_cosine(double mu) {
  return fmax(-1.0, fmin(1.0, mu));
}

static inline double clamp_distance(double d) { return fmax(d, 0.0); }

static inline double clamp_radius(double r) { return fmax(ER, fmin(r, AR)); }

static inline double safe_sqrt(double a) { return sqrt(fmax(a, 0.0)); }

static inline double distance_to_space(const double radius, const double mu) {
  const double descriminant = radius * radius * (mu * mu - 1.0) + AR * AR;
  return clamp_distance(-radius * mu + safe_sqrt(descriminant));
}

static inline double distance_to_earth(const double radius, const double mu) {
  const double descriminant = radius * radius * (mu * mu - 1.0) + ER * ER;
  return clamp_distance(-radius * mu - safe_sqrt(descriminant));
}

static inline int ray_hits_ground(double radius, double mu) {
  return mu < 0.0 && radius * radius * (mu * mu - 1.0) + ER * ER >= 0.0;
}

static inline double get_layer_density(const DensityProfileLayer *layer,
                                const double altitude) {
  double density = layer->exp_term * exp(layer->exp_scale * altitude) +
                   layer->linear_term * altitude + layer->constant_term;
  return fmax(0.0, fmin(1.0, density));
}

static inline double get_profile_density(const DensityProfile *profile,
                                  const double altitude) {
  if (altitude < profile->layers[0].width) {
    return get_layer_density(&(profile->layers[0]), altitude);
  } else {
    return get_layer_density(&(profile->layers[1]), altitude);
  }
}

static int compute_optical_length_to_space_mie(DATARRAY *mdp, const double radius,
                                               const double mu, double *result) {
  /* Number of intervals for the numerical integration. */
  const int nsamp = 500;
  /* The integration step, i.e. the length of each integration interval. */
  const double dx = distance_to_space(radius, mu) / nsamp;
  /* Integration loop. */
  for (int i = 0; i <= nsamp; ++i) {
    double d_i = i * dx;
    /* Distance between the current sample point and the planet center. */
    double r_i = sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius);
    double pt[1] = {r_i - ER};
    DATARRAY *mie = datavector(mdp, pt);
    double weight_i = i == 0 || i == nsamp ? 0.5 : 1.0;
    for (int j = 0; j < NSSAMP; ++j) {
      result[j] += mie->arr.d[j] * weight_i * dx;
    }
    free(mie);
  }
  return 1;
}

static double compute_optical_length_to_space(const DensityProfile *profile,
                                              const double radius, const double mu) {
  /* Number of intervals for the numerical integration. */
  const int SAMPLE_COUNT = 500;
  /* The integration step, i.e. the length of each integration interval. */
  const double dx = distance_to_space(radius, mu) / SAMPLE_COUNT;
  /* Integration loop. */
  double result = 0.0;
  for (int i = 0; i <= SAMPLE_COUNT; ++i) {
    double d_i = i * dx;
    /* Distance between the current sample point and the planet center. */
    double r_i = sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius);
    /* Number density at the current sample point (divided by the number density */
    /* at the bottom of the atmosphere, yielding a dimensionless number). */
    double y_i = get_profile_density(profile, r_i - ER);
    /* Sample weight (from the trapezoidal rule). */
    double weight_i = i == 0 || i == SAMPLE_COUNT ? 0.5 : 1.0;
    result += y_i * weight_i * dx;
  }
  return result;
}

static void compute_transmittance_to_space(const Atmosphere *atmos,
                                           const double radius, const double mu,
                                           float *result) {
  const double taur =
      compute_optical_length_to_space(&atmos->rayleigh_density, radius, mu);
  double taum[NSSAMP] = {0};
  compute_optical_length_to_space_mie(atmos->beta_m, radius, mu, taum);
  const double tauo =
      compute_optical_length_to_space(&atmos->ozone_density, radius, mu);
  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = exp(-(taur * atmos->beta_r0[i] + taum[i] * atmos->beta_scale +
                      tauo * AO0[i]));
  }
}

static inline double get_texture_coord_from_unit_range(const double x,
                                                       const int texture_size) {
  return 0.5 / (double)texture_size + x * (1.0 - 1.0 / (double)texture_size);
}

static inline double get_unit_range_from_texture_coord(const double u,
                                                       const int texture_size) {
  return (u - 0.5 / (double)texture_size) / (1.0 - 1.0 / (double)texture_size);
}

static void to_transmittance_uv(const double radius, const double mu, double *u,
                                double *v) {
  /* Distance to top atmosphere boundary for a horizontal ray at ground level. */
  double H = sqrt(AR * AR - ER * ER);
  /* Distance to the horizon. */
  double rho = safe_sqrt(radius * radius - ER * ER);
  /* Distance to the top atmosphere boundary for the ray (r,mu), and its minimum */
  /* and maximum values over all mu - obtained for (r,1) and (r,mu_horizon). */
  double d = distance_to_space(radius, mu);
  double d_min = AR - radius;
  double d_max = rho + H;
  double x_mu = (d - d_min) / (d_max - d_min);
  double x_r = rho / H;
  *u = fmax(0.0, fmin(1.0, x_mu));
  *v = fmax(0.0, fmin(1.0, x_r));
}

static void from_transmittance_uv(const double u, const double v, double *radius,
                                  double *mu) {
  double x_mu = u;
  double x_r = v;
  /* Distance to top atmosphere boundary for a horizontal ray at ground level. */
  double H = sqrt(AR * AR - ER * ER);
  /* Distance to the horizon, from which we can compute r: */
  double rho = H * x_r;
  *radius = sqrt(rho * rho + ER * ER);
  /* Distance to the top atmosphere boundary for the ray (r,mu), and its minimum */
  /* and maximum values over all mu - obtained for (r,1) and (r,mu_horizon) - */
  /* from which we can recover mu: */
  double d_min = AR - *radius;
  double d_max = rho + H;
  double d = d_min + x_mu * (d_max - d_min);
  *mu = d == 0.0 ? 1.0 : (H * H - rho * rho - d * d) / (2.0 * *radius * d);
  *mu = clamp_cosine(*mu);
}

static DATARRAY *get_transmittance_to_space(DATARRAY *dp, const double radius,
                                            const double mu) {
  double u, v;
  to_transmittance_uv(radius, mu, &u, &v);

  double pt[2] = {v, u};
  return datavector(dp, pt);
}

static void get_transmittance(DATARRAY *tau_dp, const double radius, const double mu, const double d,
                              const int intersects_ground, double *result) {

  DATARRAY *result1;
  DATARRAY *result2;

  const double r_d = clamp_radius(sqrt(d * d + 2.0 * radius * mu * d + radius * radius));
  const double mu_d = clamp_cosine((radius * mu + d) / r_d);

  if (intersects_ground) {
    result1 = get_transmittance_to_space(tau_dp, r_d, -mu_d);
    result2 = get_transmittance_to_space(tau_dp, radius, -mu);
  } else {
    result1 = get_transmittance_to_space(tau_dp, radius, mu);
    result2 = get_transmittance_to_space(tau_dp, r_d, mu_d);
  }
  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = fmin(result1->arr.d[i] / result2->arr.d[i], 1.0);
  }
  free(result1);
  free(result2);
}

void get_transmittance_to_sun(DATARRAY *tau_dp, const double r,
                                     const double mu_s, double *result) {
  const double sin_theta_h = ER / r;
  const double cos_theta_h = -sqrt(fmax(1.0 - sin_theta_h * sin_theta_h, 0.0));
  DATARRAY *tau_sun = get_transmittance_to_space(tau_dp, r, mu_s);
  /* Smooth step */
  const double edge0 = -sin_theta_h * SUNRAD;
  const double edge1 = edge0 * -1;
  double x = mu_s - cos_theta_h;
  x = (x - edge0) / (edge1 - edge0);
  x = x < 0.0 ? 0.0 : x > 1.0 ? 1.0 : x;
  const double st = x * x * (3.0 - 2.0 * x);
  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = tau_sun->arr.d[i] * st;
  }
  free(tau_sun);
}

static void compute_single_scattering_integrand(
    const Atmosphere *atmos, DATARRAY *tau_dp, const double radius, const double mu,
    const double mu_s, const double nu, const double distance,
    const int r_mu_intersects_ground, double *rayleigh, double *mie) {

  const double radius_d = clamp_radius(sqrt(distance * distance + 2.0 * radius * mu * distance + radius * radius));
  const double mu_s_d = clamp_cosine((radius * mu_s + distance * nu) / radius_d);
  /* double transmittance[NSSAMP]; */
  double tau_r_mu[NSSAMP] = {0};
  double tau_sun[NSSAMP] = {0};
  get_transmittance(tau_dp, radius, mu, distance, r_mu_intersects_ground, tau_r_mu);
  get_transmittance_to_sun(tau_dp, radius_d, mu_s_d, tau_sun);
  const double rayleigh_profile_density =
      get_profile_density(&atmos->rayleigh_density, radius_d - ER);
  DATARRAY *mie_scat;
  double pt[1] = {radius_d - ER};
  mie_scat = datavector(atmos->beta_m, pt);
  for (int i = 0; i < NSSAMP; ++i) {
    double transmittance = tau_r_mu[i] * tau_sun[i];
    rayleigh[i] = transmittance * rayleigh_profile_density;
    mie[i] = transmittance * mie_scat->arr.d[i];
  }
  free(mie_scat);
}

static double
distance_to_nearst_atmosphere_boundary(const double radius, const double mu,
                                       const int hits_ground) {
  if (hits_ground) {
    return distance_to_earth(radius, mu);
  } else {
    return distance_to_space(radius, mu);
  }
}

static void compute_single_scattering(const Atmosphere *atmos, DATARRAY *tau_dp,
                                      const double radius, const double mu,
                                      const double mu_s, const double nu,
                                      const int hits_ground,
                                      float *rayleigh, float *mie) {

  const int nsamp = 50;
  double distance_to_boundary;

  if (hits_ground) {
    distance_to_boundary = distance_to_earth(radius, mu);
  } else {
    distance_to_boundary = distance_to_space(radius, mu);
  }
  const double dx = distance_to_boundary / (double)nsamp;

  double rayleigh_sum[NSSAMP] = {0};
  double mie_sum[NSSAMP] = {0};

  for (int i = 0; i <= nsamp; ++i) {
    const double d_i = i * dx;
    double rayleigh_i[NSSAMP] = {0};
    double mie_i[NSSAMP] = {0};
    compute_single_scattering_integrand(atmos, tau_dp, radius, mu, mu_s, nu, d_i,
                                        hits_ground, rayleigh_i,
                                        mie_i);
    double weight_i = (i == 0 || i == nsamp) ? 0.5 : 1.0;
    for (int j = 0; j < NSSAMP; ++j) {
      rayleigh_sum[j] += rayleigh_i[j] * weight_i;
      mie_sum[j] += mie_i[j] * weight_i;
    }
  }
  for (int i = 0; i < NSSAMP; ++i) {
    rayleigh[i] = rayleigh_sum[i] * dx * EXTSOL[i] * atmos->beta_r0[i];
    mie[i] = dx * EXTSOL[i] * (mie_sum[i] * atmos->beta_scale);
  }
}

inline static double rayleigh_phase_function(const double nu) {
  const double k = 3.0 / (16.0 * PI);
  return k * (1.0 + nu * nu);
}

inline static double mie_phase_function(const double g, const double nu) {
  const double k = 3.0 / (8.0 * PI) * (1.0 - g * g) / (2.0 + g * g);
  return k * (1.0 + nu * nu) / pow(1.0 + g * g - 2.0 * g * nu, 1.5);
}

inline static double cloud_phase_function(const double g, const double nu) {
  const double alpha = 0.5;
  const double g1 = -0.5;
  const double hg0 = mie_phase_function(g, nu);
  const double hg1 = mie_phase_function(g1, nu);
  return alpha * hg0 + (1.0 - alpha) * hg1;
}

static void to_scattering_uvwz(const double radius, const double mu, const double mu_s, const double nu,
                        const int r_mu_hits_ground, double *u, double *v,
                        double *w, double *z) {

  double H = sqrt(AR * AR - ER * ER);
  double rho = safe_sqrt(radius * radius - ER * ER);
  double u_r = rho / H;

  double r_mu = radius * mu;
  double discriminant = r_mu * r_mu - radius * radius + ER * ER;

  double u_mu;
  if (r_mu_hits_ground) {
    double d = -r_mu - safe_sqrt(discriminant);
    double d_min = radius - ER;
    double d_max = rho;
    u_mu = 0.5 - 0.5 * get_texture_coord_from_unit_range(
                           d_max == d_min ? 0.0 : (d - d_min) / (d_max - d_min),
                           SCATTERING_TEXTURE_MU_SIZE / 2);
  } else {
    double d = -r_mu + sqrt(discriminant + H * H);
    double d_min = AR - radius;
    double d_max = rho + H;
    u_mu = 0.5 + 0.5 * get_texture_coord_from_unit_range(
                           (d - d_min) / (d_max - d_min),
                           SCATTERING_TEXTURE_MU_SIZE / 2);
  }

  double d = distance_to_space(ER, mu_s);
  double d_min = AH;
  double d_max = H;
  double a = (d - d_min) / (d_max - d_min);
  double D = distance_to_space(ER, MU_S_MIN);
  double A = (D - d_min) / (d_max - d_min);
  double u_mu_s = get_texture_coord_from_unit_range(
      fmax(1.0 - a / A, 0.0) / (1.0 + a), SCATTERING_TEXTURE_MU_S_SIZE);
  double u_nu = (nu + 1.0) / 2;
  *u = u_nu;
  *v = u_mu_s;
  *w = u_mu;
  *z = u_r;
}

static void from_scattering_uvwz(const double x, const double y, const double z, const double w,
                                 double *radius, double *mu, double *mu_s,
                                 double *nu, int *r_mu_hits_ground) {


  const double H = sqrt(AR * AR - ER * ER);
  const double rho = H * w;
  *radius = sqrt(rho * rho + ER * ER);

  if (z < 0.5) {
    /* Distance to the ground for the ray (r,mu), and its minimum and maximum */
    /* values over all mu - obtained for (r,-1) and (r,mu_horizon) - from which */
    /* we can recover mu: */
    const double d_min = *radius - ER;
    const double d_max = rho;
    const double d = d_min + (d_max - d_min) *
                           get_unit_range_from_texture_coord(
                               1.0 - 2.0 * z, SCATTERING_TEXTURE_MU_SIZE / 2);
    *mu = d == 0.0 ? -1.0 : clamp_cosine(-(rho * rho + d * d) / (2.0 * *radius * d));
    *r_mu_hits_ground = 1;
  } else {
    /* Distance to the top atmosphere boundary for the ray (r,mu), and its */
    /* minimum and maximum values over all mu - obtained for (r,1) and */
    /* (r,mu_horizon) - from which we can recover mu: */
    const double d_min = AR - *radius;
    const double d_max = rho + H;
    const double d = d_min + (d_max - d_min) *
                           get_unit_range_from_texture_coord(
                               2.0 * z - 1.0, SCATTERING_TEXTURE_MU_SIZE / 2);
    *mu = d == 0.0
              ? 1.0
              : clamp_cosine((H * H - rho * rho - d * d) / (2.0 * (*radius) * d));
    *r_mu_hits_ground = 0;
  }

  const double x_mu_s = y;
  const double d_min = AH;
  const double d_max = H;
  const double D = distance_to_space(ER, MU_S_MIN);
  const double A = (D - d_min) / (d_max - d_min);
  const double a = (A - x_mu_s * A) / (1.0 + x_mu_s * A);
  const double d = d_min + fmin(a, A) * (d_max - d_min);
  *mu_s = d == 0.0 ? 1.0 : clamp_cosine((H * H - d * d) / (2.0 * ER * d));
  *nu = clamp_cosine(x * 2.0 - 1.0);
}

static DATARRAY *interpolate_scattering(DATARRAY *dp, const double radius,
                                        const double mu, const double mu_s,
                                        const double nu,
                                        const int r_mu_hits_ground) {
  double u, v, w, z;
  to_scattering_uvwz(radius, mu, mu_s, nu, r_mu_hits_ground, &u, &v, &w,
                     &z);
  double pt[4] = {z, w, v, u};
  return datavector(dp, pt);
}

static void get_scattering(DATARRAY *scat1r, DATARRAY *scat1m, DATARRAY *scat,
                           const double radius, const double mu, const double mu_s,
                           const double nu,
                           const int ray_r_mu_intersects_ground,
                           const int scattering_order, double *result) {
  if (scattering_order == 1) {
    DATARRAY *rayleigh = interpolate_scattering(scat1r, radius, mu, mu_s, nu,
                                                ray_r_mu_intersects_ground);
    DATARRAY *mie = interpolate_scattering(scat1m, radius, mu, mu_s, nu,
                                           ray_r_mu_intersects_ground);
    const double rayleigh_phase = rayleigh_phase_function(nu);
    const double mie_phase = mie_phase_function(MIE_G, nu);
    for (int i = 0; i < NSSAMP; ++i) {
      result[i] =
          rayleigh->arr.d[i] * rayleigh_phase + mie->arr.d[i] * mie_phase;
    }
    free(rayleigh);
    free(mie);
  } else {
    DATARRAY *scattering = interpolate_scattering(scat, radius, mu, mu_s, nu,
                                                  ray_r_mu_intersects_ground);
    for (int i = 0; i < NSSAMP; ++i) {
      result[i] = scattering->arr.d[i];
    }
    free(scattering);
  }
}

static void to_irradiance_uv(const double radius, const double mu_s, double *u,
                             double *v) {
  const double x_r = (radius - ER) / AH;
  const double x_mu_s = mu_s * 0.5 + 0.5;
  *u = x_mu_s;
  *v = x_r;
}

static void from_irradiance_uv(const double u, const double v, double *radius,
                               double *mu_s) {
  *radius = v * AH + ER;
  *mu_s = clamp_cosine(u * 2.0 - 1.0);
}

DATARRAY *get_indirect_irradiance(DATARRAY *dp, const double radius,
                                const double mu_s) {
  double u, v;
  to_irradiance_uv(radius, mu_s, &u, &v);
  double pt[2] = {v, u};
  return datavector(dp, pt);
}

static void
compute_scattering_density(const Atmosphere *atmos, DATARRAY *tau_dp,
                           DATARRAY *scat1r, DATARRAY *scat1m, DATARRAY *scat,
                           DATARRAY *irrad_dp, const double radius, const double mu,
                           const double mu_s, const double nu,
                           const int scattering_order, double *result) {
  /* Compute unit direction vectors for the zenith, the view direction omega and */
  /* and the sun direction omega_s, such that the cosine of the view-zenith */
  /* angle is mu, the cosine of the sun-zenith angle is mu_s, and the cosine of */
  /* the view-sun angle is nu. The goal is to simplify computations below. */
  const double height = radius - ER;
  const double epsilon = 1e-6;
  const FVECT zenith_direction = {0.0, 0.0, 1.0};
  double omega_0 = 1.0 - mu * mu;
  omega_0 = fabs(omega_0) <= epsilon ? 0.0 : sqrt(omega_0);
  const FVECT omega = {omega_0, 0.0, mu};
  const double sun_dir_x = omega[0] == 0.0 ? 0.0 : (nu - mu * mu_s) / omega[0];
  const double sun_dir_y =
      sqrt(fmax(1.0 - sun_dir_x * sun_dir_x - mu_s * mu_s, 0.0));
  const FVECT omega_s = {sun_dir_x, sun_dir_y, mu_s};

  const int SAMPLE_COUNT = 16;
  const double dphi = PI / SAMPLE_COUNT;
  const double dtheta = PI / SAMPLE_COUNT;

  /* Nested loops for the integral over all the incident directions omega_i. */
  for (int l = 0; l < SAMPLE_COUNT; ++l) {
    const double theta = (l + 0.5) * dtheta;
    const double cos_theta = cos(theta);
    const double sin_theta = sin(theta);
    const int r_theta_hits_ground = ray_hits_ground(radius, cos_theta);

    /* The distance and transmittance to the ground only depend on theta, so we */
    /* can compute them in the outer loop for efficiency. */
    double distance_to_ground = 0.0;
    double transmittance_to_ground[NSSAMP] = {0};
    double ground_albedo = 0;
    if (r_theta_hits_ground) {
      distance_to_ground = distance_to_earth(radius, cos_theta);
      /* assert(distance_to_ground >= 0.0 && distance_to_ground <= HMAX); */
      get_transmittance(tau_dp, radius, cos_theta, distance_to_ground, r_theta_hits_ground,
                        transmittance_to_ground);
      ground_albedo = atmos->grefl;
    }

    for (int m = 0; m < 2 * SAMPLE_COUNT; ++m) {
      const double phi = (m + 0.5) * dphi;
      const FVECT omega_i = {cos(phi) * sin_theta, sin(phi) * sin_theta,
                             cos_theta};
      const double domega_i = dtheta * dphi * sin(theta);

      /* The radiance L_i arriving from direction omega_i after n-1 bounces is */
      /* the sum of a term given by the precomputed scattering texture for the */
      /* (n-1)-th order: */
      double nu1 = fdot(omega_s, omega_i);
      if (nu1 <= -1.0) {
        nu1 += 0.1;
      } else if (nu1 >= 1.0) {
        nu1 -= 0.1;
      }
      double incident_radiance[NSSAMP] = {0};
      get_scattering(scat1r, scat1m, scat, radius, omega_i[2], mu_s, nu1,
                     r_theta_hits_ground, scattering_order - 1,
                     incident_radiance);

      /* and of the contribution from the light paths with n-1 bounces and whose */
      /* last bounce is on the ground. This contribution is the product of the */
      /* transmittance to the ground, the ground albedo, the ground BRDF, and */
      /* the irradiance received on the ground after n-2 bounces. */
      FVECT ground_normal = {
          zenith_direction[0] * radius + omega_i[0] * distance_to_ground,
          zenith_direction[1] * radius + omega_i[1] * distance_to_ground,
          zenith_direction[2] * radius + omega_i[2] * distance_to_ground};
      normalize(ground_normal);
      DATARRAY *ground_irradiance =
          get_indirect_irradiance(irrad_dp, ER, fdot(ground_normal, omega_s));
      for (int i = 0; i < NSSAMP; ++i) {
        incident_radiance[i] += transmittance_to_ground[i] * ground_albedo *
                                (1.0 / PI) * ground_irradiance->arr.d[i];
      }
      free(ground_irradiance);

      /* The radiance finally scattered from direction omega_i towards direction */
      /* -omega is the product of the incident radiance, the scattering */
      /* coefficient, and the phase function for directions omega and omega_i */
      /* (all this summed over all particle types, i.e. Rayleigh and Mie). */
      const double nu2 = fdot(omega, omega_i);
      const double rayleigh_density =
          get_profile_density(&atmos->rayleigh_density, height);
      DATARRAY *mie_scat;
      double pt[1] = {height};
      mie_scat = datavector(atmos->beta_m, pt);
      const double rayleigh_phase = rayleigh_phase_function(nu2);
      const double mie_phase = mie_phase_function(MIE_G, nu2);
      for (int j = 0; j < NSSAMP; ++j) {
        result[j] += incident_radiance[j] * domega_i *
                     (atmos->beta_r0[j] * rayleigh_density * rayleigh_phase +
                      mie_scat->arr.d[j] * mie_phase * atmos->beta_scale);
      }
      free(mie_scat);
    }
  }
}

static void compute_multi_scattering(DATARRAY *tau_dp,
                                        DATARRAY *scattering_density_dp,
                                        const double radius, const double mu,
                                        const double mu_s, const double nu,
                                        const int r_mu_hits_ground,
                                        double *result) {


  /* Numerical integration sample count. */
  const int nsamp = 50;
  const double dx = distance_to_nearst_atmosphere_boundary(
                        radius, mu, r_mu_hits_ground) /
                    (double)nsamp;

  for (int i = 0; i <= nsamp; ++i) {
    const double d_i = i * dx;

    /* The r, mu and mu_s parameters at the current integration point (see the */
    /* single scattering section for a detailed explanation). */
    const double r_i =
        clamp_radius(sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius));
    const double mu_i = clamp_cosine((radius * mu + d_i) / r_i);
    const double mu_s_i = clamp_cosine((radius * mu_s + d_i * nu) / r_i);

    /* The Rayleigh and Mie multiple scattering at the current sample point. */
    double transmittance[NSSAMP] = {0};
    DATARRAY *rayleigh_mie_s =
        interpolate_scattering(scattering_density_dp, r_i, mu_i, mu_s_i, nu,
                               r_mu_hits_ground);
    get_transmittance(tau_dp, radius, mu, d_i, r_mu_hits_ground,
                      transmittance);
    /* Sample weight (from the trapezoidal rule). */
    const double weight_i = (i == 0 || i == nsamp) ? 0.5 : 1.0;
    for (int j = 0; j < NSSAMP; ++j) {
      result[j] += rayleigh_mie_s->arr.d[j] * transmittance[j] * dx * weight_i;
    }
    free(rayleigh_mie_s);
  }
}

static void compute_direct_irradiance(DATARRAY *tau_dp, const double r,
                                      const double mu_s, float *result) {

  /* Approximate the average of the cosine factor mu_s over the visible fraction */
  /* of the sun disc */
  const double average_cosine_factor =
      mu_s < -SUNRAD
          ? 0.0
          : (mu_s > SUNRAD
                 ? mu_s
                 : (mu_s + SUNRAD) * (mu_s + SUNRAD) / (4.0 * SUNRAD));

  DATARRAY *transmittance = get_transmittance_to_space(tau_dp, r, mu_s);
  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = EXTSOL[i] * transmittance->arr.d[i] * average_cosine_factor;
  }
  free(transmittance);
}

static void compute_indirect_irradiance(DATARRAY *scat1r, DATARRAY *scat1m,
                                        DATARRAY *scat, const double radius,
                                        const double mu_s,
                                        const int scattering_order,
                                        double *result) {

  const int SAMPLE_COUNT = 32;
  const double dphi = PI / SAMPLE_COUNT;
  const double dtheta = PI / SAMPLE_COUNT;

  const FVECT omega_s = {sqrt(1.0 - mu_s * mu_s), 0.0, mu_s};
  for (int j = 0; j < SAMPLE_COUNT / 2; ++j) {
    const double theta = (j + 0.5) * dtheta;
    for (int i = 0; i < 2 * SAMPLE_COUNT; ++i) {
      const double phi = (i + 0.5) * dphi;
      const FVECT omega = {sin(theta) * cos(phi), sin(theta) * sin(phi),
                           cos(theta)};
      const double domega = dtheta * dphi * sin(theta);
      const double nu = fdot(omega, omega_s);
      double result1[NSSAMP] = {0};
      get_scattering(scat1r, scat1m, scat, radius, omega[2], mu_s, nu, 0,
                     scattering_order, result1);
      for (int k = 0; k < NSSAMP; ++k) {
        result[k] += result1[k] * omega[2] * domega;
      }
    }
  }
}

DATARRAY *allocate_3d_datarray(const char *name, const int ri, const int rj,
                               const int rk) {
  DATARRAY *dp;
  if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
    goto memerr;
  int asize = ri * rj * rk;
  dp->name = savestr(name);
  dp->type = DATATY;
  dp->nd = 3;
  dp->dim[0].org = 0;
  dp->dim[0].siz = 1;
  dp->dim[0].ne = ri;
  dp->dim[0].p = NULL;
  dp->dim[1].org = 0;
  dp->dim[1].siz = 1;
  dp->dim[1].ne = rj;
  dp->dim[1].p = NULL;
  dp->dim[2].org = 0;
  dp->dim[2].siz = rk - 1;
  dp->dim[2].ne = rk;
  dp->dim[2].p = NULL;
  /* dp->arr.p = */
  /* malloc(sizeof(DATATYPE)*dp->dim[0].ne*dp->dim[1].ne*dp->dim[2].ne); */
  if ((dp->arr.d = (DATATYPE *)malloc(asize * sizeof(DATATYPE))) == NULL)
    goto memerr;
  return (dp);

memerr:
  fprintf(stderr, "Memory allocation error in allocate_3d_datarray\n");
  return (NULL);
}

DATARRAY *allocate_5d_datarray(const char *name, const int ri, const int rj,
                               const int rk, const int rl, const int rm) {
  DATARRAY *dp;
  int asize = ri * rj * rk * rl * rm;
  if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
    goto memerr;
  dp->name = savestr(name);
  dp->type = DATATY;
  dp->nd = 5;
  dp->dim[0].org = 0;
  dp->dim[0].siz = 1;
  dp->dim[0].ne = ri;
  dp->dim[0].p = NULL;
  dp->dim[1].org = 0;
  dp->dim[1].siz = 1;
  dp->dim[1].ne = rj;
  dp->dim[1].p = NULL;
  dp->dim[2].org = 0;
  dp->dim[2].siz = 1;
  dp->dim[2].ne = rk;
  dp->dim[2].p = NULL;
  dp->dim[3].org = 0;
  dp->dim[3].siz = 1;
  dp->dim[3].ne = rl;
  dp->dim[3].p = NULL;
  dp->dim[4].org = 0;
  dp->dim[4].siz = rm - 1;
  dp->dim[4].ne = rm;
  dp->dim[4].p = NULL;
  if ((dp->arr.d = (DATATYPE *)malloc(asize * sizeof(DATATYPE))) == NULL)
    goto memerr;
  return (dp);
memerr:
  fprintf(stderr, "Memory allocation error in allocate_5d_datarray\n");
  return (NULL);
}

/* save to a .dat file */
void savedata(DATARRAY *dp) {
  if (dp == NULL)
    return;
  FILE *fp = fopen(dp->name, "w");
  if (fp == NULL) {
    fprintf(stderr, "Error opening file %s\n", dp->name);
    return;
  }
  int i, j;
  int nvals = 1;
  for (i = 0; i < dp->nd; i++) {
    nvals *= dp->dim[i].ne;
  }
  fprintf(fp, "%d\n", dp->nd);
  for (i = 0; i < dp->nd; i++) {
    fprintf(fp, "%f %f %d\n", dp->dim[i].org, dp->dim[i].siz, dp->dim[i].ne);
  }
  for (i = 0; i < nvals; i++) {
    fprintf(fp, "%f\n", dp->arr.d[i]);
  }
  fclose(fp);
}

void increment_dp(DATARRAY *dp1, DATARRAY *dp2) {
  if (dp1 == NULL || dp2 == NULL)
    perror("null pointer in increment_dp\n");

  if (dp1->nd != dp2->nd)
    perror("dimension mismatch in increment_dp\n");

  int i;
  int nvals1 = 1;
  int nvals2 = 1;
  for (i = 0; i < dp1->nd; i++) {
    nvals1 *= dp1->dim[i].ne;
    nvals2 *= dp2->dim[i].ne;
  }
  if (nvals1 != nvals2)
    perror("size mismatch in increment_dp\n");
  for (i = 0; i < nvals1; i++) {
    dp1->arr.d[i] += dp2->arr.d[i];
  }
}

THREAD_RETURN compute_scattering_density_thread(void *arg) {
  ScatDenseTdat *tdata = (ScatDenseTdat *)arg;
  for (unsigned int l = tdata->start_l; l < tdata->end_l; ++l) {
    for (unsigned int k = 0; k < SCATTERING_TEXTURE_MU_SIZE; ++k) {
      for (unsigned int j = 0; j < SCATTERING_TEXTURE_MU_S_SIZE; ++j) {
        for (unsigned int i = 0; i < SCATTERING_TEXTURE_NU_SIZE; ++i) {
          double scattering_density[NSSAMP] = {0};
          double r, mu, mu_s, nu;
          int ray_r_mu_intersects_ground;
          double xr = (double)i / SCATTERING_TEXTURE_NU_SIZE;
          double yr = (double)j / SCATTERING_TEXTURE_MU_S_SIZE;
          double zr = (double)k / SCATTERING_TEXTURE_MU_SIZE;
          double wr = (double)l / SCATTERING_TEXTURE_R_SIZE;
          from_scattering_uvwz(xr, yr, zr, wr, &r, &mu, &mu_s, &nu,
                               &ray_r_mu_intersects_ground);
          nu =
              fmax(mu * mu_s - sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)),
                   fmin(mu * mu_s + sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)),
                        nu));
          compute_scattering_density(
              tdata->atmos, tdata->tau_dp, tdata->delta_rayleigh_scattering_dp,
              tdata->delta_mie_scattering_dp,
              tdata->delta_multiple_scattering_dp, tdata->delta_irradiance_dp,
              r, mu, mu_s, nu, tdata->scattering_order, scattering_density);
          for (int m = 0; m < NSSAMP; ++m) {
            int idx = l * SCATTERING_TEXTURE_MU_SIZE *
                          SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      k * SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      j * SCATTERING_TEXTURE_NU_SIZE * NSSAMP + i * NSSAMP + m;
            tdata->delta_scattering_density_dp->arr.d[idx] =
                scattering_density[m];
          }
        }
      }
    }
  }
  return 0;
}

THREAD_RETURN compute_single_scattering_thread(void *arg) {
  Scat1Tdat *tdata = (Scat1Tdat *)arg;
  for (unsigned int l = tdata->start_l; l < tdata->end_l; ++l) {
    for (int k = 0; k < SCATTERING_TEXTURE_MU_SIZE; ++k) {
      for (int j = 0; j < SCATTERING_TEXTURE_MU_S_SIZE; ++j) {
        for (int i = 0; i < SCATTERING_TEXTURE_NU_SIZE; ++i) {
          float rayleigh[NSSAMP] = {0};
          float mie[NSSAMP] = {0};
          double r, mu, mu_s, nu;
          int r_mu_hits_ground;
          const double xr = (double)i / SCATTERING_TEXTURE_NU_SIZE;
          const double yr = (double)j / SCATTERING_TEXTURE_MU_S_SIZE;
          const double zr = (double)k / SCATTERING_TEXTURE_MU_SIZE;
          const double wr = (double)l / SCATTERING_TEXTURE_R_SIZE;
          from_scattering_uvwz(xr, yr, zr, wr, &r, &mu, &mu_s, &nu,
                               &r_mu_hits_ground);
          const double mu_mu_s = mu * mu_s;
          const double sqrt_term = sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s));
          nu = fmax(mu_mu_s - sqrt_term, fmin(mu_mu_s + sqrt_term, nu));
          compute_single_scattering(tdata->atmos, tdata->tau_dp, r, mu, mu_s,
                                    nu, r_mu_hits_ground, rayleigh,
                                    mie);
          for (int m = 0; m < NSSAMP; ++m) {
            int idx = l * SCATTERING_TEXTURE_MU_SIZE *
                          SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      k * SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      j * SCATTERING_TEXTURE_NU_SIZE * NSSAMP + i * NSSAMP + m;
            tdata->delta_rayleigh_scattering_dp->arr.d[idx] = rayleigh[m];
            tdata->delta_mie_scattering_dp->arr.d[idx] = mie[m];
            tdata->scattering_dp->arr.d[idx] = rayleigh[m];
          }
        }
      }
    }
  }
  return 0;
}

THREAD_RETURN compute_multi_scattering_thread(void *arg) {
  ScatNTdat *tdata = (ScatNTdat *)arg;
  for (unsigned int l = tdata->start_l; l < tdata->end_l; ++l) {
    for (unsigned int k = 0; k < SCATTERING_TEXTURE_MU_SIZE; ++k) {
      for (unsigned int j = 0; j < SCATTERING_TEXTURE_MU_S_SIZE; ++j) {
        for (unsigned int i = 0; i < SCATTERING_TEXTURE_NU_SIZE; ++i) {
          double delta_multiple_scattering[NSSAMP] = {0};
          double r, mu, mu_s, nu;
          int r_mu_hits_ground;
          const double xr = (double)i / SCATTERING_TEXTURE_NU_SIZE;
          const double yr = (double)j / SCATTERING_TEXTURE_MU_S_SIZE;
          const double zr = (double)k / SCATTERING_TEXTURE_MU_SIZE;
          const double wr = (double)l / SCATTERING_TEXTURE_R_SIZE;
          from_scattering_uvwz(xr, yr, zr, wr, &r, &mu, &mu_s, &nu,
                               &r_mu_hits_ground);
          const double mu_mu_s = mu * mu_s;
          const double sqrt_term = sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s));
          nu = fmax(mu_mu_s - sqrt_term, fmin(mu_mu_s + sqrt_term, nu));
          compute_multi_scattering(
              tdata->tau_dp, tdata->delta_scattering_density_dp, r, mu, mu_s,
              nu, r_mu_hits_ground, delta_multiple_scattering);
          const double rayleigh_phase = rayleigh_phase_function(nu);
          for (int m = 0; m < NSSAMP; ++m) {
            int idx = l * SCATTERING_TEXTURE_MU_SIZE *
                          SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      k * SCATTERING_TEXTURE_MU_S_SIZE *
                          SCATTERING_TEXTURE_NU_SIZE * NSSAMP +
                      j * SCATTERING_TEXTURE_NU_SIZE * NSSAMP + i * NSSAMP + m;
            tdata->delta_multiple_scattering_dp->arr.d[idx] =
                delta_multiple_scattering[m];
            tdata->scattering_dp->arr.d[idx] +=
                delta_multiple_scattering[m] * (1.0 / rayleigh_phase);
          }
        }
      }
    }
  }
  return 0;
}

int precompute(const int sorder, const DpPaths dppaths, const Atmosphere *atmos,
               int num_threads) {
  num_threads = (num_threads < SCATTERING_TEXTURE_R_SIZE)
                    ? num_threads
                    : SCATTERING_TEXTURE_R_SIZE;
  const int tchunk = SCATTERING_TEXTURE_R_SIZE / num_threads;
  const int tremainder = SCATTERING_TEXTURE_R_SIZE % num_threads;

  if (sorder < 2) {
    fprintf(stderr, "scattering order must be at least 2\n");
    return 0;
  }

  DATARRAY *tau_dp =
      allocate_3d_datarray(dppaths.tau, TRANSMITTANCE_TEXTURE_HEIGHT,
                           TRANSMITTANCE_TEXTURE_WIDTH, NSSAMP);
  DATARRAY *delta_irradiance_dp = allocate_3d_datarray(
      "ssky_delta_irradiance.dat", IRRADIANCE_TEXTURE_HEIGHT,
      IRRADIANCE_TEXTURE_WIDTH, NSSAMP);
  DATARRAY *irradiance_dp =
      allocate_3d_datarray(dppaths.irrad, IRRADIANCE_TEXTURE_HEIGHT,
                           IRRADIANCE_TEXTURE_WIDTH, NSSAMP);

  DATARRAY *delta_rayleigh_scattering_dp = allocate_5d_datarray(
      "ssky_delta_rayleigh_scattering_dp.dat", SCATTERING_TEXTURE_R_SIZE,
      SCATTERING_TEXTURE_MU_SIZE, SCATTERING_TEXTURE_MU_S_SIZE,
      SCATTERING_TEXTURE_NU_SIZE, NSSAMP);
  DATARRAY *delta_mie_scattering_dp = allocate_5d_datarray(
      dppaths.scat1m, SCATTERING_TEXTURE_R_SIZE, SCATTERING_TEXTURE_MU_SIZE,
      SCATTERING_TEXTURE_MU_S_SIZE, SCATTERING_TEXTURE_NU_SIZE, NSSAMP);
  DATARRAY *scattering_dp = allocate_5d_datarray(
      dppaths.scat, SCATTERING_TEXTURE_R_SIZE, SCATTERING_TEXTURE_MU_SIZE,
      SCATTERING_TEXTURE_MU_S_SIZE, SCATTERING_TEXTURE_NU_SIZE, NSSAMP);
  DATARRAY *delta_multiple_scattering_dp = allocate_5d_datarray(
      "ssky_delta_multiple_scattering.dat", SCATTERING_TEXTURE_R_SIZE,
      SCATTERING_TEXTURE_MU_SIZE, SCATTERING_TEXTURE_MU_S_SIZE,
      SCATTERING_TEXTURE_NU_SIZE, NSSAMP);
  DATARRAY *delta_scattering_density_dp = allocate_5d_datarray(
      "ssky_delta_scattering_density.dat", SCATTERING_TEXTURE_R_SIZE,
      SCATTERING_TEXTURE_MU_SIZE, SCATTERING_TEXTURE_MU_S_SIZE,
      SCATTERING_TEXTURE_NU_SIZE, NSSAMP);

  int idx = 0;
  /* Computing transmittance... */
  for (int j = 0; j < TRANSMITTANCE_TEXTURE_HEIGHT; ++j) {
    for (int i = 0; i < TRANSMITTANCE_TEXTURE_WIDTH; ++i) {
      double r, mu;
      float tau[NSSAMP] = {0};
      const float xr = (float)i / TRANSMITTANCE_TEXTURE_WIDTH;
      const float yr = (float)j / TRANSMITTANCE_TEXTURE_HEIGHT;
      from_transmittance_uv(xr, yr, &r, &mu);
      compute_transmittance_to_space(atmos, r, mu, tau);
      for (int k = 0; k < NSSAMP; ++k) {
        tau_dp->arr.d[idx] = tau[k];
        idx += 1;
      }
    }
  }
  savedata(tau_dp);

  /* Computing direct irradiance... */
  idx = 0;
  for (int j = 0; j < IRRADIANCE_TEXTURE_HEIGHT; ++j) {
    for (int i = 0; i < IRRADIANCE_TEXTURE_WIDTH; ++i) {
      double r, mu_s;
      float result[NSSAMP] = {0};
      const double xr = (double)i / IRRADIANCE_TEXTURE_WIDTH;
      const double yr = (double)j / IRRADIANCE_TEXTURE_HEIGHT;
      from_irradiance_uv(xr, yr, &r, &mu_s);
      compute_direct_irradiance(tau_dp, r, mu_s, result);
      for (int k = 0; k < NSSAMP; ++k) {
        delta_irradiance_dp->arr.d[idx] = result[k];
        irradiance_dp->arr.d[idx] = 0.0;
        idx += 1;
      }
    }
  }

  /* Computing single scattering... */
  THREAD *threads = (THREAD *)malloc(num_threads * sizeof(THREAD));
  Scat1Tdat *tdata = (Scat1Tdat *)malloc(num_threads * sizeof(Scat1Tdat));
  for (int i = 0; i < num_threads; i++) {
    tdata[i].atmos = atmos;
    tdata[i].tau_dp = tau_dp;
    tdata[i].delta_rayleigh_scattering_dp = delta_rayleigh_scattering_dp;
    tdata[i].delta_mie_scattering_dp = delta_mie_scattering_dp;
    tdata[i].scattering_dp = scattering_dp;
    tdata[i].start_l = i * tchunk;
    tdata[i].end_l = (i + 1) * tchunk;
    if (i == num_threads - 1)
      tdata[i].end_l += tremainder;
    CREATE_THREAD(&threads[i], compute_single_scattering_thread,
                  (void *)&tdata[i]);
  }

  for (int i = 0; i < num_threads; i++) {
    THREAD_JOIN(threads[i]);
#if defined(_WIN32) || defined(_WIN64)
    THREAD_CLOSE(threads[i]);
#endif
  }
  free(tdata);
  free(threads);
  savedata(delta_mie_scattering_dp);


  for (int scattering_order = 2; scattering_order <= sorder; ++scattering_order) {
    THREAD *threads = (THREAD *)malloc(num_threads * sizeof(THREAD));
    ScatDenseTdat *tdata = (ScatDenseTdat *)malloc(num_threads * sizeof(ScatDenseTdat));
    for (int i = 0; i < num_threads; i++) {
      tdata[i].atmos = atmos;
      tdata[i].tau_dp = tau_dp;
      tdata[i].delta_rayleigh_scattering_dp = delta_rayleigh_scattering_dp;
      tdata[i].delta_mie_scattering_dp = delta_mie_scattering_dp;
      tdata[i].delta_multiple_scattering_dp = delta_multiple_scattering_dp;
      tdata[i].delta_irradiance_dp = delta_irradiance_dp;
      tdata[i].scattering_order = scattering_order;
      tdata[i].delta_scattering_density_dp = delta_scattering_density_dp;
      tdata[i].start_l = i * tchunk;
      tdata[i].end_l = (i + 1) * tchunk;
      if (i == num_threads - 1)
        tdata[i].end_l += tremainder;
      CREATE_THREAD(&threads[i], compute_scattering_density_thread,
                    (void *)&tdata[i]);
    }

    for (int i = 0; i < num_threads; i++) {
      THREAD_JOIN(threads[i]);
#if defined(_WIN32) || defined(_WIN64)
      THREAD_CLOSE(threads[i]);
#endif
    }
    free(tdata);
    free(threads);

    idx = 0;
    for (unsigned int j = 0; j < IRRADIANCE_TEXTURE_HEIGHT; ++j) {
      for (unsigned int i = 0; i < IRRADIANCE_TEXTURE_WIDTH; ++i) {
        double delta_irradiance[NSSAMP] = {0};
        double r, mu_s;
        const double xr = (double)i / IRRADIANCE_TEXTURE_WIDTH;
        const double yr = (double)j / IRRADIANCE_TEXTURE_HEIGHT;
        from_irradiance_uv(xr, yr, &r, &mu_s);
        compute_indirect_irradiance(delta_rayleigh_scattering_dp,
                                    delta_mie_scattering_dp,
                                    delta_multiple_scattering_dp, r, mu_s,
                                    scattering_order - 1, delta_irradiance);
        for (int k = 0; k < NSSAMP; ++k) {
          delta_irradiance_dp->arr.d[idx] = delta_irradiance[k];
          idx += 1;
        }
      }
    }

    increment_dp(irradiance_dp, delta_irradiance_dp);

    THREAD *threads2 = (THREAD *)malloc(num_threads * sizeof(THREAD));
    ScatNTdat *tdata2 = (ScatNTdat *)malloc(num_threads * sizeof(ScatNTdat));
    for (int i = 0; i < num_threads; i++) {
      tdata2[i].tau_dp = tau_dp;
      tdata2[i].delta_multiple_scattering_dp = delta_multiple_scattering_dp;
      tdata2[i].delta_scattering_density_dp = delta_scattering_density_dp;
      tdata2[i].scattering_dp = scattering_dp;
      tdata2[i].start_l = i * tchunk;
      tdata2[i].end_l = (i + 1) * tchunk;
      if (i == num_threads - 1)
        tdata2[i].end_l += tremainder;
      CREATE_THREAD(&threads2[i], compute_multi_scattering_thread,
                    (void *)&tdata2[i]);
    }

    for (int i = 0; i < num_threads; i++) {
      THREAD_JOIN(threads2[i]);
#if defined(_WIN32) || defined(_WIN64)
      THREAD_CLOSE(threads2[i]);
#endif
    }
    free(tdata2);
    free(threads2);
  }
  savedata(scattering_dp);
  savedata(irradiance_dp);

  freedata(tau_dp);
  freedata(delta_rayleigh_scattering_dp);
  freedata(delta_mie_scattering_dp);
  freedata(scattering_dp);
  freedata(delta_multiple_scattering_dp);
  freedata(delta_scattering_density_dp);
  freedata(delta_irradiance_dp);
  freedata(irradiance_dp);
  return 1;
}

int compute_sundir(const int year, const int month, const int day,
                   const double hour, const int tsolar, double sundir[3]) {
  double sd, st = hour;

  if (year) {
    /* Michalsky algorithm? */
    double mjd = mjdate(year, month, day, hour);
    if (tsolar)
      sd = msdec(mjd, NULL);
    else
      sd = msdec(mjd, &st);
  } else {
    int jd = jdate(month, day); /* Julian date */
    sd = sdec(jd);              /* solar declination */
    if (!tsolar)                /* get solar time? */
      st = hour + stadj(jd);
  }
  const double altitude = salt(sd, st);
  const double azimuth = sazi(sd, st);
  sundir[0] = -sin(azimuth) * cos(altitude);
  sundir[1] = -cos(azimuth) * cos(altitude);
  sundir[2] = sin(altitude);
  return 1;
}

void get_sky_transmittance(DATARRAY *tau, const double radius, const double mu, float *result) {
  DATARRAY *trans = get_transmittance_to_space(tau, radius, mu);
  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = trans->arr.d[i];
  }
  free(trans);
}

void get_sky_radiance(DATARRAY *scat, DATARRAY *scat1m, const double radius, 
                      const double mu, const double mu_s, const double nu, float *result) {
  double u, v, w, z;
  to_scattering_uvwz(radius, mu, mu_s, nu, 0, &u, &v, &w, &z);
  double pt[4] = {z, w, v, u};

  DATARRAY *scattering = datavector(scat, pt);
  DATARRAY *single_mie_scattering = datavector(scat1m, pt);

  const double rayleigh_phase = rayleigh_phase_function(nu);
  const double mie_phase = mie_phase_function(MIE_G, nu);

  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = scattering->arr.d[i] * rayleigh_phase +
                single_mie_scattering->arr.d[i] * mie_phase;
  }
  free(single_mie_scattering);
  free(scattering);
}

static void get_irradiance(DATARRAY *tau, DATARRAY *irrad, const double radius, 
                             const FVECT point, const FVECT normal, 
                             const FVECT sundir, double *result) {
  double mu_s = fdot(point, sundir) / radius;

  double point_trans_sun[NSSAMP] = {0};
  get_transmittance_to_sun(tau, radius, mu_s, point_trans_sun);

  DATARRAY *indirect_irradiance = get_indirect_irradiance(irrad, radius, mu_s);
  const double sun_ct = fmax(fdot(normal, sundir), 0.0);

  for (int i = 0; i < NSSAMP; ++i) {
    result[i] = indirect_irradiance->arr.d[i] + point_trans_sun[i] * EXTSOL[i] * sun_ct;
  }

  free(indirect_irradiance);
}


void get_ground_radiance(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, DATARRAY *irrad, 
                         const FVECT view_point, const FVECT view_direction, 
                         const double radius, const double mu, const double sun_ct, 
                         const double nu, const double grefl, const FVECT sundir, 
                         float *ground_radiance) {

  const int intersect = ray_hits_ground(radius, mu);

  if (intersect) {
    FVECT point, normal;
    const double distance = distance_to_earth(radius, mu);

    /* direct + indirect irradiance */
    VSUM(point, view_point, view_direction, distance);
    VCOPY(normal, point);
    normalize(normal);
    double irradiance[NSSAMP] = {0};
    get_irradiance(tau, irrad, ER, point, normal, sundir, irradiance);

    /* transmittance between view point and ground point */
    double trans[NSSAMP] = {0};
    get_transmittance(tau, radius, mu, distance, intersect, trans);
    if (trans[0] == 0) {printf("trans 0\n");}

    /* inscattering */
    float inscatter[NSSAMP] = {0}; 
    get_sky_radiance(scat, scat1m, radius, mu, sun_ct, nu, inscatter);

    for (int i = 0; i < NSSAMP; ++i) {
      ground_radiance[i] = inscatter[i] + irradiance[i] * trans[i] * grefl / PI;
    }
  }
}

void get_solar_radiance(DATARRAY *tau, DATARRAY *scat, DATARRAY *scat1m, const FVECT sundir, const double radius, const double sun_ct, double *sun_radiance) {
  double trans_sun[NSSAMP] = {0};
  float sky_radiance_sun[NSSAMP] = {0};
  double u, v, w, z;
  const double nu = fdot(sundir, sundir);
  get_transmittance_to_sun(tau, radius, sun_ct, trans_sun);
  get_sky_radiance(scat, scat1m, radius, sun_ct, sun_ct, nu, sky_radiance_sun);
  for (int i = 0; i < NSSAMP; ++i) {
    sun_radiance[i] = sky_radiance_sun[i] + trans_sun[i] * EXTSOL[i] / SOLOMG;
  }
}