ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/atmos.c
(Generate patch)

Comparing ray/src/gen/atmos.c (file contents):
Revision 2.1 by greg, Fri Jul 5 18:04:36 2024 UTC vs.
Revision 2.2 by greg, Fri Jul 19 23:38:28 2024 UTC

# Line 1 | Line 1
1 + #ifndef lint
2 + static const char RCSid[] = "$Id$";
3 + #endif
4 + /*
5 + The Radiance Software License, Version 2.0
6 +
7 + Radiance v5.4 Copyright (c) 1990 to 2022, The Regents of the University of
8 + California, through Lawrence Berkeley National Laboratory (subject to receipt
9 + of any required approvals from the U.S. Dept. of Energy).  All rights reserved.
10 +
11 + Redistribution and use in source and binary forms, with or without
12 + modification, are permitted provided that the following conditions are met:
13 +
14 + (1) Redistributions of source code must retain the above copyright notice,
15 + this list of conditions and the following disclaimer.
16 +
17 + (2) Redistributions in binary form must reproduce the above copyright
18 + notice, this list of conditions and the following disclaimer in the
19 + documentation and/or other materials provided with the distribution.
20 +
21 + (3) Neither the name of the University of California, Lawrence Berkeley
22 + National Laboratory, U.S. Dept. of Energy nor the names of its contributors
23 + may be used to endorse or promote products derived from this software
24 + without specific prior written permission.
25 +
26 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
27 + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
29 + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
30 + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
31 + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
32 + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
33 + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
34 + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
35 + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
36 + POSSIBILITY OF SUCH DAMAGE.
37 +
38 + You are under no obligation whatsoever to provide any bug fixes, patches,
39 + or upgrades to the features, functionality or performance of the source
40 + code ("Enhancements") to anyone; however, if you choose to make your
41 + Enhancements available either publicly, or directly to Lawrence Berkeley
42 + National Laboratory, without imposing a separate written license agreement
43 + for such Enhancements, then you hereby grant the following license: a
44 + non-exclusive, royalty-free perpetual license to install, use, modify,
45 + prepare derivative works, incorporate into other computer software,
46 + distribute, and sublicense such enhancements or derivative works thereof,
47 + in binary and source code form.
48 +
49 + =================================================================================
50 +
51 + Copyright (c) 2017 Eric Bruneton
52 + All rights reserved.
53 +
54 + Redistribution and use in source and binary forms, with or without
55 + modification, are permitted provided that the following conditions are met:
56 +
57 + 1. Redistributions of source code must retain the above copyright notice, this
58 + list of conditions and the following disclaimer.
59 +
60 + 2. Redistributions in binary form must reproduce the above copyright notice,
61 + this list of conditions and the following disclaimer in the documentation
62 + and/or other materials provided with the distribution.
63 +
64 + 3. Neither the name of the copyright holder nor the names of its contributors
65 + may be used to endorse or promote products derived from this software without
66 + specific prior written permission.
67 +
68 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
69 + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
70 + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
71 + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
72 + FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
73 + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
74 + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
75 + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
76 + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
77 + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
78 +
79 + Additional modifications by T. Wang
80 + - N-dimensional data storage and lookup using Radiance data structure
81 + - Mie scattering coefficients precomputed using libRadTrans and interpolated at runtime
82 + - Translation to C from C++ with multi-threading
83 + */
84 +
85   #include "atmos.h"
86   #include "data.h"
87  
# Line 55 | Line 139 | typedef struct {
139    DATARRAY *delta_scattering_density_dp;
140   } ScatDenseTdat;
141  
142 < const double ER = 6360000.0; // Earth radius in meters
143 < const double AR = 6420000; // Atmosphere radius in meters
144 < const double AH = 60000;   // Atmosphere thickness in meters
142 > const double ER = 6360000.0; /* Earth radius in meters */
143 > const double AR = 6420000; /* Atmosphere radius in meters */
144 > const double AH = 60000;   /* Atmosphere thickness in meters */
145  
146   const float START_WVL = 390.0;
147   const float END_WVL = 770.0;
148   const double SUNRAD = 0.004625;
149  
150 < // The cosine of the maximum Sun zenith angle
150 > /* The cosine of the maximum Sun zenith angle */
151   const double MU_S_MIN = -0.2;
152  
153 < // Scale heights (m)
154 < // Rayleigh scattering
155 < const double HR_MS = 8820; // Midlatitude summer
156 < const double HR_MW = 8100; // Midlatitude winter
157 < const double HR_SS = 8550; // Subarctic summer
158 < const double HR_SW = 7600; // Subarctic winter
159 < const double HR_T = 9050;  // Tropical
153 > /* Scale heights (m) */
154 > /* Rayleigh scattering */
155 > const double HR_MS = 8820; /* Midlatitude summer */
156 > const double HR_MW = 8100; /* Midlatitude winter */
157 > const double HR_SS = 8550; /* Subarctic summer */
158 > const double HR_SW = 7600; /* Subarctic winter */
159 > const double HR_T = 9050;  /* Tropical */
160  
161 < const double AOD0_CA = 0.115; // Broadband AOD for continental average
161 > const double AOD0_CA = 0.115; /* Broadband AOD for continental average */
162  
163 < const double SOLOMG = 6.7967e-5; // .533 apex angle
164 < const int WVLSPAN = 400;         // 380nm to 780nm
163 > const double SOLOMG = 6.7967e-5; /* .533 apex angle */
164 > const int WVLSPAN = 400;         /* 380nm to 780nm */
165  
166 < // Aerosol optical depth at 550nm for continental clean
166 > /* Aerosol optical depth at 550nm for continental clean */
167   const double AOD_CC_550 = 0.05352;
168  
169   /** 380nm to 780nm at 20nm intervals **/
170 < // Extraterrestrial solar W/m^2/nm
170 > /* Extraterrestrial solar W/m^2/nm */
171   const float EXTSOL[NSSAMP] = {1.11099345, 1.75363199, 1.67126921, 1.97235667,
172                                2.01863015, 1.93612482, 1.87310758, 1.88792588,
173                                1.86274213, 1.84433248, 1.80370942, 1.73020456,
174                                1.67036654, 1.57482149, 1.53646383, 1.45515319,
175                                1.38207435, 1.3265141,  1.26648111, 1.20592043};
176  
177 < // Rayleight scattering coefficients at sea level in m^-1
177 > /* Rayleight scattering coefficients at sea level in m^-1 */
178   const float BR0_MS[NSSAMP] = {
179      4.63889967e-05, 3.76703293e-05, 3.09193474e-05, 2.56223081e-05,
180      2.14198405e-05, 1.80483398e-05, 1.53177709e-05, 1.30867298e-05,
# Line 125 | Line 209 | const float BR0_T[NSSAMP] = {
209      6.35240772e-06, 5.59544593e-06, 4.94840517e-06, 4.39219003e-06,
210      3.91235655e-06, 3.49641927e-06, 3.13434084e-06, 2.81804244e-06};
211  
212 < // Ozone extinction coefficients at sea level
212 > /* Ozone extinction coefficients at sea level */
213   const float AO0[NSSAMP] = {
214      3.5661864e-09, 1.4848362e-08, 4.5415674e-08, 1.2446184e-07, 2.6461576e-07,
215      4.8451984e-07, 8.609148e-07,  1.480537e-06,  1.8809e-06,    2.5107328e-06,
216      2.5263174e-06, 2.313507e-06,  1.727741e-06,  1.2027012e-06, 7.915902e-07,
217      5.0639202e-07, 3.5285684e-07, 2.23021e-07,   1.7395638e-07, 1.5052574e-07};
218  
219 < const double MIE_G = 0.7;  // Mie eccentricity
219 > const double MIE_G = 0.7;  /* Mie eccentricity */
220  
221   const int TRANSMITTANCE_TEXTURE_WIDTH = 256;
222   const int TRANSMITTANCE_TEXTURE_HEIGHT = 64;
223  
224 < // Resolution halved by original impl
224 > /* Resolution halved by original impl */
225   const int SCATTERING_TEXTURE_R_SIZE = 16;
226   const int SCATTERING_TEXTURE_MU_SIZE = 64;
227   const int SCATTERING_TEXTURE_MU_S_SIZE = 16;
# Line 188 | Line 272 | static inline double get_profile_density(const Density
272  
273   static int compute_optical_length_to_space_mie(DATARRAY *mdp, const double radius,
274                                                 const double mu, double *result) {
275 <  // Number of intervals for the numerical integration.
275 >  /* Number of intervals for the numerical integration. */
276    const int nsamp = 500;
277 <  // The integration step, i.e. the length of each integration interval.
277 >  /* The integration step, i.e. the length of each integration interval. */
278    const double dx = distance_to_space(radius, mu) / nsamp;
279 <  // Integration loop.
279 >  /* Integration loop. */
280    for (int i = 0; i <= nsamp; ++i) {
281      double d_i = i * dx;
282 <    // Distance between the current sample point and the planet center.
282 >    /* Distance between the current sample point and the planet center. */
283      double r_i = sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius);
284      double pt[1] = {r_i - ER};
285      DATARRAY *mie = datavector(mdp, pt);
# Line 210 | Line 294 | static int compute_optical_length_to_space_mie(DATARRA
294  
295   static double compute_optical_length_to_space(const DensityProfile *profile,
296                                                const double radius, const double mu) {
297 <  // Number of intervals for the numerical integration.
297 >  /* Number of intervals for the numerical integration. */
298    const int SAMPLE_COUNT = 500;
299 <  // The integration step, i.e. the length of each integration interval.
299 >  /* The integration step, i.e. the length of each integration interval. */
300    const double dx = distance_to_space(radius, mu) / SAMPLE_COUNT;
301 <  // Integration loop.
301 >  /* Integration loop. */
302    double result = 0.0;
303    for (int i = 0; i <= SAMPLE_COUNT; ++i) {
304      double d_i = i * dx;
305 <    // Distance between the current sample point and the planet center.
305 >    /* Distance between the current sample point and the planet center. */
306      double r_i = sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius);
307 <    // Number density at the current sample point (divided by the number density
308 <    // at the bottom of the atmosphere, yielding a dimensionless number).
307 >    /* Number density at the current sample point (divided by the number density */
308 >    /* at the bottom of the atmosphere, yielding a dimensionless number). */
309      double y_i = get_profile_density(profile, r_i - ER);
310 <    // Sample weight (from the trapezoidal rule).
310 >    /* Sample weight (from the trapezoidal rule). */
311      double weight_i = i == 0 || i == SAMPLE_COUNT ? 0.5 : 1.0;
312      result += y_i * weight_i * dx;
313    }
# Line 257 | Line 341 | static inline double get_unit_range_from_texture_coord
341  
342   static void to_transmittance_uv(const double radius, const double mu, double *u,
343                                  double *v) {
344 <  // Distance to top atmosphere boundary for a horizontal ray at ground level.
344 >  /* Distance to top atmosphere boundary for a horizontal ray at ground level. */
345    double H = sqrt(AR * AR - ER * ER);
346 <  // Distance to the horizon.
346 >  /* Distance to the horizon. */
347    double rho = safe_sqrt(radius * radius - ER * ER);
348 <  // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
349 <  // and maximum values over all mu - obtained for (r,1) and (r,mu_horizon).
348 >  /* Distance to the top atmosphere boundary for the ray (r,mu), and its minimum */
349 >  /* and maximum values over all mu - obtained for (r,1) and (r,mu_horizon). */
350    double d = distance_to_space(radius, mu);
351    double d_min = AR - radius;
352    double d_max = rho + H;
# Line 276 | Line 360 | static void from_transmittance_uv(const double u, cons
360                                    double *mu) {
361    double x_mu = u;
362    double x_r = v;
363 <  // Distance to top atmosphere boundary for a horizontal ray at ground level.
363 >  /* Distance to top atmosphere boundary for a horizontal ray at ground level. */
364    double H = sqrt(AR * AR - ER * ER);
365 <  // Distance to the horizon, from which we can compute r:
365 >  /* Distance to the horizon, from which we can compute r: */
366    double rho = H * x_r;
367    *radius = sqrt(rho * rho + ER * ER);
368 <  // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
369 <  // and maximum values over all mu - obtained for (r,1) and (r,mu_horizon) -
370 <  // from which we can recover mu:
368 >  /* Distance to the top atmosphere boundary for the ray (r,mu), and its minimum */
369 >  /* and maximum values over all mu - obtained for (r,1) and (r,mu_horizon) - */
370 >  /* from which we can recover mu: */
371    double d_min = AR - *radius;
372    double d_max = rho + H;
373    double d = d_min + x_mu * (d_max - d_min);
# Line 328 | Line 412 | void get_transmittance_to_sun(DATARRAY *tau_dp, const
412    const double sin_theta_h = ER / r;
413    const double cos_theta_h = -sqrt(fmax(1.0 - sin_theta_h * sin_theta_h, 0.0));
414    DATARRAY *tau_sun = get_transmittance_to_space(tau_dp, r, mu_s);
415 <  // Smooth step
415 >  /* Smooth step */
416    const double edge0 = -sin_theta_h * SUNRAD;
417    const double edge1 = edge0 * -1;
418    double x = mu_s - cos_theta_h;
# Line 348 | Line 432 | static void compute_single_scattering_integrand(
432  
433    const double radius_d = clamp_radius(sqrt(distance * distance + 2.0 * radius * mu * distance + radius * radius));
434    const double mu_s_d = clamp_cosine((radius * mu_s + distance * nu) / radius_d);
435 <  // double transmittance[NSSAMP];
435 >  /* double transmittance[NSSAMP]; */
436    double tau_r_mu[NSSAMP] = {0};
437    double tau_sun[NSSAMP] = {0};
438    get_transmittance(tau_dp, radius, mu, distance, r_mu_intersects_ground, tau_r_mu);
# Line 485 | Line 569 | static void from_scattering_uvwz(const double x, const
569    *radius = sqrt(rho * rho + ER * ER);
570  
571    if (z < 0.5) {
572 <    // Distance to the ground for the ray (r,mu), and its minimum and maximum
573 <    // values over all mu - obtained for (r,-1) and (r,mu_horizon) - from which
574 <    // we can recover mu:
572 >    /* Distance to the ground for the ray (r,mu), and its minimum and maximum */
573 >    /* values over all mu - obtained for (r,-1) and (r,mu_horizon) - from which */
574 >    /* we can recover mu: */
575      const double d_min = *radius - ER;
576      const double d_max = rho;
577      const double d = d_min + (d_max - d_min) *
# Line 496 | Line 580 | static void from_scattering_uvwz(const double x, const
580      *mu = d == 0.0 ? -1.0 : clamp_cosine(-(rho * rho + d * d) / (2.0 * *radius * d));
581      *r_mu_hits_ground = 1;
582    } else {
583 <    // Distance to the top atmosphere boundary for the ray (r,mu), and its
584 <    // minimum and maximum values over all mu - obtained for (r,1) and
585 <    // (r,mu_horizon) - from which we can recover mu:
583 >    /* Distance to the top atmosphere boundary for the ray (r,mu), and its */
584 >    /* minimum and maximum values over all mu - obtained for (r,1) and */
585 >    /* (r,mu_horizon) - from which we can recover mu: */
586      const double d_min = AR - *radius;
587      const double d_max = rho + H;
588      const double d = d_min + (d_max - d_min) *
# Line 588 | Line 672 | compute_scattering_density(const Atmosphere *atmos, DA
672                             DATARRAY *irrad_dp, const double radius, const double mu,
673                             const double mu_s, const double nu,
674                             const int scattering_order, double *result) {
675 <  // Compute unit direction vectors for the zenith, the view direction omega and
676 <  // and the sun direction omega_s, such that the cosine of the view-zenith
677 <  // angle is mu, the cosine of the sun-zenith angle is mu_s, and the cosine of
678 <  // the view-sun angle is nu. The goal is to simplify computations below.
675 >  /* Compute unit direction vectors for the zenith, the view direction omega and */
676 >  /* and the sun direction omega_s, such that the cosine of the view-zenith */
677 >  /* angle is mu, the cosine of the sun-zenith angle is mu_s, and the cosine of */
678 >  /* the view-sun angle is nu. The goal is to simplify computations below. */
679    const double height = radius - ER;
680    const double epsilon = 1e-6;
681    const FVECT zenith_direction = {0.0, 0.0, 1.0};
# Line 607 | Line 691 | compute_scattering_density(const Atmosphere *atmos, DA
691    const double dphi = PI / SAMPLE_COUNT;
692    const double dtheta = PI / SAMPLE_COUNT;
693  
694 <  // Nested loops for the integral over all the incident directions omega_i.
694 >  /* Nested loops for the integral over all the incident directions omega_i. */
695    for (int l = 0; l < SAMPLE_COUNT; ++l) {
696      const double theta = (l + 0.5) * dtheta;
697      const double cos_theta = cos(theta);
698      const double sin_theta = sin(theta);
699      const int r_theta_hits_ground = ray_hits_ground(radius, cos_theta);
700  
701 <    // The distance and transmittance to the ground only depend on theta, so we
702 <    // can compute them in the outer loop for efficiency.
701 >    /* The distance and transmittance to the ground only depend on theta, so we */
702 >    /* can compute them in the outer loop for efficiency. */
703      double distance_to_ground = 0.0;
704      double transmittance_to_ground[NSSAMP] = {0};
705      double ground_albedo = 0;
706      if (r_theta_hits_ground) {
707        distance_to_ground = distance_to_earth(radius, cos_theta);
708 <      // assert(distance_to_ground >= 0.0 && distance_to_ground <= HMAX);
708 >      /* assert(distance_to_ground >= 0.0 && distance_to_ground <= HMAX); */
709        get_transmittance(tau_dp, radius, cos_theta, distance_to_ground, r_theta_hits_ground,
710                          transmittance_to_ground);
711        ground_albedo = atmos->grefl;
# Line 633 | Line 717 | compute_scattering_density(const Atmosphere *atmos, DA
717                               cos_theta};
718        const double domega_i = dtheta * dphi * sin(theta);
719  
720 <      // The radiance L_i arriving from direction omega_i after n-1 bounces is
721 <      // the sum of a term given by the precomputed scattering texture for the
722 <      // (n-1)-th order:
720 >      /* The radiance L_i arriving from direction omega_i after n-1 bounces is */
721 >      /* the sum of a term given by the precomputed scattering texture for the */
722 >      /* (n-1)-th order: */
723        double nu1 = fdot(omega_s, omega_i);
724        if (nu1 <= -1.0) {
725          nu1 += 0.1;
# Line 647 | Line 731 | compute_scattering_density(const Atmosphere *atmos, DA
731                       r_theta_hits_ground, scattering_order - 1,
732                       incident_radiance);
733  
734 <      // and of the contribution from the light paths with n-1 bounces and whose
735 <      // last bounce is on the ground. This contribution is the product of the
736 <      // transmittance to the ground, the ground albedo, the ground BRDF, and
737 <      // the irradiance received on the ground after n-2 bounces.
734 >      /* and of the contribution from the light paths with n-1 bounces and whose */
735 >      /* last bounce is on the ground. This contribution is the product of the */
736 >      /* transmittance to the ground, the ground albedo, the ground BRDF, and */
737 >      /* the irradiance received on the ground after n-2 bounces. */
738        FVECT ground_normal = {
739            zenith_direction[0] * radius + omega_i[0] * distance_to_ground,
740            zenith_direction[1] * radius + omega_i[1] * distance_to_ground,
# Line 664 | Line 748 | compute_scattering_density(const Atmosphere *atmos, DA
748        }
749        free(ground_irradiance);
750  
751 <      // The radiance finally scattered from direction omega_i towards direction
752 <      // -omega is the product of the incident radiance, the scattering
753 <      // coefficient, and the phase function for directions omega and omega_i
754 <      // (all this summed over all particle types, i.e. Rayleigh and Mie).
751 >      /* The radiance finally scattered from direction omega_i towards direction */
752 >      /* -omega is the product of the incident radiance, the scattering */
753 >      /* coefficient, and the phase function for directions omega and omega_i */
754 >      /* (all this summed over all particle types, i.e. Rayleigh and Mie). */
755        const double nu2 = fdot(omega, omega_i);
756        const double rayleigh_density =
757            get_profile_density(&atmos->rayleigh_density, height);
# Line 694 | Line 778 | static void compute_multi_scattering(DATARRAY *tau_dp,
778                                          double *result) {
779  
780  
781 <  // Numerical integration sample count.
781 >  /* Numerical integration sample count. */
782    const int nsamp = 50;
783    const double dx = distance_to_nearst_atmosphere_boundary(
784                          radius, mu, r_mu_hits_ground) /
# Line 703 | Line 787 | static void compute_multi_scattering(DATARRAY *tau_dp,
787    for (int i = 0; i <= nsamp; ++i) {
788      const double d_i = i * dx;
789  
790 <    // The r, mu and mu_s parameters at the current integration point (see the
791 <    // single scattering section for a detailed explanation).
790 >    /* The r, mu and mu_s parameters at the current integration point (see the */
791 >    /* single scattering section for a detailed explanation). */
792      const double r_i =
793          clamp_radius(sqrt(d_i * d_i + 2.0 * radius * mu * d_i + radius * radius));
794      const double mu_i = clamp_cosine((radius * mu + d_i) / r_i);
795      const double mu_s_i = clamp_cosine((radius * mu_s + d_i * nu) / r_i);
796  
797 <    // The Rayleigh and Mie multiple scattering at the current sample point.
797 >    /* The Rayleigh and Mie multiple scattering at the current sample point. */
798      double transmittance[NSSAMP] = {0};
799      DATARRAY *rayleigh_mie_s =
800          interpolate_scattering(scattering_density_dp, r_i, mu_i, mu_s_i, nu,
801                                 r_mu_hits_ground);
802      get_transmittance(tau_dp, radius, mu, d_i, r_mu_hits_ground,
803                        transmittance);
804 <    // Sample weight (from the trapezoidal rule).
804 >    /* Sample weight (from the trapezoidal rule). */
805      const double weight_i = (i == 0 || i == nsamp) ? 0.5 : 1.0;
806      for (int j = 0; j < NSSAMP; ++j) {
807        result[j] += rayleigh_mie_s->arr.d[j] * transmittance[j] * dx * weight_i;
# Line 729 | Line 813 | static void compute_multi_scattering(DATARRAY *tau_dp,
813   static void compute_direct_irradiance(DATARRAY *tau_dp, const double r,
814                                        const double mu_s, float *result) {
815  
816 <  // Approximate the average of the cosine factor mu_s over the visible fraction
817 <  // of the sun disc
816 >  /* Approximate the average of the cosine factor mu_s over the visible fraction */
817 >  /* of the sun disc */
818    const double average_cosine_factor =
819        mu_s < -SUNRAD
820            ? 0.0
# Line 795 | Line 879 | DATARRAY *allocate_3d_datarray(const char *name, const
879    dp->dim[2].siz = rk - 1;
880    dp->dim[2].ne = rk;
881    dp->dim[2].p = NULL;
882 <  // dp->arr.p =
883 <  // malloc(sizeof(DATATYPE)*dp->dim[0].ne*dp->dim[1].ne*dp->dim[2].ne);
882 >  /* dp->arr.p = */
883 >  /* malloc(sizeof(DATATYPE)*dp->dim[0].ne*dp->dim[1].ne*dp->dim[2].ne); */
884    if ((dp->arr.d = (DATATYPE *)malloc(asize * sizeof(DATATYPE))) == NULL)
885      goto memerr;
886    return (dp);
# Line 1053 | Line 1137 | int precompute(const int sorder, const DpPaths dppaths
1137        SCATTERING_TEXTURE_NU_SIZE, NSSAMP);
1138  
1139    int idx = 0;
1140 <  // Computing transmittance...
1140 >  /* Computing transmittance... */
1141    for (int j = 0; j < TRANSMITTANCE_TEXTURE_HEIGHT; ++j) {
1142      for (int i = 0; i < TRANSMITTANCE_TEXTURE_WIDTH; ++i) {
1143        double r, mu;
# Line 1070 | Line 1154 | int precompute(const int sorder, const DpPaths dppaths
1154    }
1155    savedata(tau_dp);
1156  
1157 <  // Computing direct irradiance...
1157 >  /* Computing direct irradiance... */
1158    idx = 0;
1159    for (int j = 0; j < IRRADIANCE_TEXTURE_HEIGHT; ++j) {
1160      for (int i = 0; i < IRRADIANCE_TEXTURE_WIDTH; ++i) {
# Line 1088 | Line 1172 | int precompute(const int sorder, const DpPaths dppaths
1172      }
1173    }
1174  
1175 <  // Computing single scattering...
1175 >  /* Computing single scattering... */
1176    THREAD *threads = (THREAD *)malloc(num_threads * sizeof(THREAD));
1177    Scat1Tdat *tdata = (Scat1Tdat *)malloc(num_threads * sizeof(Scat1Tdat));
1178    for (int i = 0; i < num_threads; i++) {
# Line 1116 | Line 1200 | int precompute(const int sorder, const DpPaths dppaths
1200    savedata(delta_mie_scattering_dp);
1201  
1202  
1119 –  // Compute the 2nd, 3rd and 4th order of scattering, in sequence.
1203    for (int scattering_order = 2; scattering_order <= sorder; ++scattering_order) {
1121 –    // Compute the scattering density, and store it in
1122 –    // delta_scattering_density_texture.
1204      THREAD *threads = (THREAD *)malloc(num_threads * sizeof(THREAD));
1205      ScatDenseTdat *tdata = (ScatDenseTdat *)malloc(num_threads * sizeof(ScatDenseTdat));
1206      for (int i = 0; i < num_threads; i++) {
# Line 1148 | Line 1229 | int precompute(const int sorder, const DpPaths dppaths
1229      free(tdata);
1230      free(threads);
1231  
1151 –    // Compute the indirect irradiance, store it in delta_irradiance_texture
1152 –    // and accumulate it in irradiance_texture_.
1232      idx = 0;
1233      for (unsigned int j = 0; j < IRRADIANCE_TEXTURE_HEIGHT; ++j) {
1234        for (unsigned int i = 0; i < IRRADIANCE_TEXTURE_WIDTH; ++i) {
# Line 1171 | Line 1250 | int precompute(const int sorder, const DpPaths dppaths
1250  
1251      increment_dp(irradiance_dp, delta_irradiance_dp);
1252  
1174 –    // Computing multiple scattering...
1253      THREAD *threads2 = (THREAD *)malloc(num_threads * sizeof(THREAD));
1254      ScatNTdat *tdata2 = (ScatNTdat *)malloc(num_threads * sizeof(ScatNTdat));
1255      for (int i = 0; i < num_threads; i++) {
# Line 1294 | Line 1372 | void get_ground_radiance(DATARRAY *tau, DATARRAY *scat
1372      FVECT point, normal;
1373      const double distance = distance_to_earth(radius, mu);
1374  
1375 <    // direct + indirect irradiance
1375 >    /* direct + indirect irradiance */
1376      VSUM(point, view_point, view_direction, distance);
1377      VCOPY(normal, point);
1378      normalize(normal);
1379      double irradiance[NSSAMP] = {0};
1380      get_irradiance(tau, irrad, ER, point, normal, sundir, irradiance);
1381  
1382 <    // transmittance between view point and ground point
1382 >    /* transmittance between view point and ground point */
1383      double trans[NSSAMP] = {0};
1384      get_transmittance(tau, radius, mu, distance, intersect, trans);
1385      if (trans[0] == 0) {printf("trans 0\n");}
1386  
1387 <    // inscattering
1387 >    /* inscattering */
1388      float inscatter[NSSAMP] = {0};
1389      get_sky_radiance(scat, scat1m, radius, mu, sun_ct, nu, inscatter);
1390  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines