Appendix / Radiative Transfer
|
|
Rama
physics Appendix / Radiative Transfer |
This section explains how we integrate the equations in the Appendix numerically, and how we use this to compute the coupled atmosphere and ground lighting and temperature profiles. The explanations are interspersed with the corresponding source code, using a kind of literate programming style. The source code can be divided in several sections:
We start with the definition of the physical and Rama constants that we will need:
#include <fstream> #include <iostream> #include <cmath> // ----------------------------------------------------------------------------- // Physical and Rama constants // ----------------------------------------------------------------------------- // Stefan-Bolztmann constant, in W.m^-2.K^-4. const double sigma = 5.67e-8; // Specific gas constant or air, in J.kg^-1.K^-1. const double Rs = 8.314 / 0.0289; // Radius of Rama's surface, in m. const double rs = 8000.0; // Angular speed of Rama in rad.s^-1. const double omega = 2.0 * M_PI / 180.0; // Ratio between inner and outer long wave radiations, unitless. // Use 1.3 (resp. 2.9) for the incandescent (resp. energy saving) scenario. const double Delta = 2.9; // Albedo of the ground, unitless. const double alpha = 0.15; // "Albedo" of the ground in the long wave range, unitless. const double alphaBar = Delta / (Delta + 1.0); // Desired average ground and atmosphere temperature, in K. const double Ta = 288.0; // Total radiance of the light sources, per meter, at the ground level, // in W.m^-2, S / 2 pi r_s (computed here to give the average temperature Ta). const double SoverTwoPiRs = sigma * pow(Ta, 4.0) / Delta; // Angular width of each light source, in radians. const double DeltaTheta = M_PI / 18000; // Scale factor used to convert kS and kE values to integers. #define KSCALE 1e6 // Scattering coefficient of the atmosphere at the ground level, in m^-1. const int kS = static_cast<int>(1.50e-5 * KSCALE); // Emissivity coefficient of the atmosphere at the ground level, in m^-1. const int kE = static_cast<int>(2.58e-4 * KSCALE);
// -----------------------------------------------------------------------------
// Scattering / Emissivity coefficients and transmittance functions
// -----------------------------------------------------------------------------
template<int K>
double k(double r) {
const double C = omega * omega / (2.0 * Rs * Ta);
const double k0 = K / (KSCALE * exp(C * rs * rs));
return k0 * exp(C * r * r);
}
We compute the transmittance by using the analytic approximation method
presented in the transmittance section, because
this is much more efficient than a numerical integration, and because this
function is heavily used in the following (note that a precomputation for all
possible parameter values is impractical too, because we have 4 parameters).
However, instead of using the approximation $e^x\approx 1+x+x^2/2$, we use a more precise, third order approximation. Indeed, the second order approximation gives a relative error of 1.5% for $x<1/2$, which is too large compared to the small radiance and temperature variations we expect in the end. The third order Taylor series $e^x\approx 1+x+x^2/2+x^3/6$ gives a 0.16% relative error, which is better. But we can do even better for the same computational cost by fitting a 3rd degree polynomial to $e^x$ for $0 < x < 1/2$, yielding $e^x\approx 1+C_1x+C_2x^2+C_3x^3$ with $C_1=1.00115$, $C_2=0.486507$ and $C_3=0.211682$. The relative error of this approximation is less than +/- 0.03%. Computing the transmittance using this approximation requires the integral of \begin{equation*} 1+C_1(at^2+bt+c)+C_2(at^2+bt+c)^2+C_3(at^2+bt+c)^3 \end{equation*} between 0 and 1, which is given by the expression v below:
template<int K>
double t(double rp, double thetap, double rq, double thetaq) {
const double C = omega * omega / (2.0 * Rs * Ta);
const double k0 = K / (KSCALE * exp(C * rs * rs));
double rp2 = rp * rp;
double rq2 = rq * rq;
double rprq = rp * rq * cos(thetaq - thetap);
double l2 = rp2 + rq2 - 2.0 * rprq;
double a = C * l2;
double b = 2.0 * C * (rprq - rp2);
double c = C * rp2;
double a2 = a * a;
double a3 = a * a2;
double b2 = b * b;
double b3 = b * b2;
double c2 = c * c;
double c3 = c * c2;
const double C1 = 1.00115;
const double C2 = 0.486507;
const double C3 = 0.211682;
double v = a3*C3/7 + a2*b*C3/2 + 3*a2*c*C3/5 + a2*C2/5 + 3*a*b2*C3/5 +
3*a*b*c*C3/2 + a*b*C2/2 + a*c2*C3 + 2*a*c*C2/3 + a*C1/3 + b3*C3/4 +
b2*c*C3 + b2*C2/3 + 3*b*c2*C3/2 + b*c*C2 + b*C1/2 + c3*C3 + c2*C2 + c*C1 +
1;
return exp(-k0 * sqrt(l2) * v);
}
Note that another way to look at this representation is to see it as a linear combination of basis functions, noted here $E^i$, which are equal to 1 at one sample location $\theta_i$ or $r_i,\theta_j$, 0 at all the other sample locations, and are linearly or bilinearly interpolated in between (with the coefficients of the linear combination being the values of the function at each sample location).
Assuming that the values of the function are stored in a valueSamples field of the class, an array of THETA_RES + 1 samples, we can start by defining some methods to initialize a function to a base function, to a constant function or to a copy of an existing function:
// -----------------------------------------------------------------------------
// Ground irradiance functions E
// -----------------------------------------------------------------------------
class Efunction {
public:
void setBaseFunction(int itheta) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] = i == itheta ? 1.0 : 0.0;
}
}
Efunction& set(double k) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] = k;
}
return *this;
}
Efunction& set(const Efunction &E) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] = E.valueSamples[i];
}
return *this;
}
We can continue by adding methods to add, multiply or divide functions by
scalars or other functions, and to compute the average of a function:
Efunction& add(double v) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] += v;
}
return *this;
}
Efunction& add(const Efunction &E) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] += E.valueSamples[i];
}
return *this;
}
Efunction& mult(double k) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] *= k;
}
return *this;
}
Efunction& div(const Efunction &E) {
for (int i = 0; i <= THETA_RES; ++i) {
valueSamples[i] /= E.valueSamples[i];
}
return *this;
}
double average() {
double sum = 0.0;
double area = 0.0;
for (int i = 0; i < THETA_RES; ++i) {
double dtheta = getThetaSample(i + 1) - getThetaSample(i);
double v = valueSamples[i + 1] + valueSamples[i];
sum += v * 0.5 * dtheta;
area += dtheta;
}
return sum / area;
}
We will also need methods to get the location and value of each sample
(the sample locations are stored in a static THETA_SAMPLES array
that will be initialized below):
static inline int getThetaSampleCount() {
return THETA_RES;
}
static inline double getThetaSample(int itheta) {
return THETA_SAMPLES[itheta];
}
inline double getValueSample(int itheta) const {
return valueSamples[itheta];
}
inline void setValueSample(int itheta, double v) {
valueSamples[itheta] = v;
}
Rama, with 3 linear light sources separated from each other by 120°, is
invariant under rotations of 120°. It is also invariant under reflections around
planes containing its axis and a linear light. As a consequence, any physical
function depending on $\theta$ will be periodic of period $2\pi/3$, and each
period will be symmetric around its center. Thus, we don't need to sample this
function over the full $[0,2\pi]$ interval to represent it: it suffice to sample
it over the $[0,\pi/3]$ interval (assuming the light sources are at
$\theta=0,2\pi/3,4\pi/3$). The location and the value of the other samples can
then be deduced by using the symmetries.
This is what the methods below do: they give the total number of samples for the $[0,2\pi]$ interval (i.e. 6 times the number of samples for the $[0,\pi/3]$ interval), the location of these samples, and their values.
static inline int getPeriodicThetaSampleCount() {
return 6 * THETA_RES;
}
static inline double getPeriodicThetaSample(int itheta) {
int part = itheta / THETA_RES;
int i0 = THETA_RES * (part % 2);
int sign = 1 - 2 * (part % 2);
double offset = (part + 1) / 2 * 2.0 * M_PI / 3.0;
return offset + sign * THETA_SAMPLES[i0 + sign * (itheta % THETA_RES)];
}
inline double getPeriodicValueSample(int itheta) const {
int part = itheta / THETA_RES;
int i0 = THETA_RES * (part % 2);
int sign = 1 - 2 * (part % 2);
return valueSamples[i0 + sign * (itheta % THETA_RES)];
}
In order to output the results of our computations, we then define some methods
to save a function to a file:
void dump(const char* file) const {
FILE* f = fopen(file, "w");
for (int i = 0; i <= THETA_RES; ++i) {
fprintf(f, "%d %.9g\n", i, valueSamples[i]);
}
fclose(f);
}
void plot(const char* file) const {
FILE *f = fopen(file, "w");
for (int i = 0; i <= getPeriodicThetaSampleCount(); ++i) {
double theta = getPeriodicThetaSample(i);
fprintf(f, "%.9g %.12g\n", theta,
getPeriodicValueSample(i) / SoverTwoPiRs);
}
fclose(f);
}
void plotTemperature(const char* file) const {
FILE *f = fopen(file, "w");
for (int i = 0; i <= getPeriodicThetaSampleCount(); ++i) {
double theta = getPeriodicThetaSample(i);
double T = pow(getPeriodicValueSample(i) / sigma, 0.25);
fprintf(f, "%.9g %.12g\n", theta, T - 273);
}
fclose(f);
}
The final step is to declare the fields of the class, and to initialize the
location of the samples in the THETA_SAMPLES array. Here we don't use a
regular sampling of the $[0,\pi/3]$ interval, because we would like to have more
samples near the light sources, i.e. near $\theta=0$ (we can suspect that large
gradients could occur here, and this is confirmed by the end results). We
therefore use a sampling scheme of the form $\theta_i=(i/N)^4\pi/3$, with
$0 < i \le N$, which puts many samples near $\theta=0$:
private:
static double* computeThetaSamples() {
double* samples = new double[THETA_RES + 1];
for (int i = 0; i <= THETA_RES; ++i) {
samples[i] = pow(i / double(THETA_RES), THETA_POWER) * M_PI / 3.0;
}
return samples;
}
static const int THETA_RES = 64;
static const double THETA_POWER = 4.0;
static double* THETA_SAMPLES;
double valueSamples[THETA_RES + 1];
};
double* Efunction::THETA_SAMPLES = Efunction::computeThetaSamples();
The definition of the class to represent atmosphere radiance functions
$J(r,\theta)$ is very similar to the above class. Note that we don't use a
regular sampling for the radius samples $r_i$, for the same reason as above. In
order to get more samples near the light sources, and near the surface in
general, we use a sampling scheme of the form $r=(1-(1-i/N)^4)r_s$, with $0 < i
\le N$, which puts many samples near $r=r_s$. Combined with the sampling for
$\theta$, this gives the sample locations shown in Fig. 1
// -----------------------------------------------------------------------------
// Average atmosphere radiance functions J
// -----------------------------------------------------------------------------
class Jfunction {
public:
void setBaseFunction(int ir, int jtheta) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] = i == ir && j == jtheta ? 1.0 : 0.0;
}
}
}
Jfunction& set(double k) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] = k;
}
}
return *this;
}
Jfunction& set(const Jfunction &J) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] = J.valueSamples[i][j];
}
}
return *this;
}
Jfunction& add(double v) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] += v;
}
}
return *this;
}
Jfunction& add(const Jfunction &J) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] += J.valueSamples[i][j];
}
}
return *this;
}
Jfunction& mult(double k) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] *= k;
}
}
return *this;
}
Jfunction& div(const Jfunction &J) {
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
valueSamples[i][j] /= J.valueSamples[i][j];
}
}
return *this;
}
double average() {
double sum = 0.0;
double area = 0.0;
for (int i = 0; i < R_RES; ++i) {
double dr = getRSample(i + 1) - getRSample(i);
for (int j = 0; j < THETA_RES; ++j) {
double dtheta = getThetaSample(j + 1) - getThetaSample(j);
double v = valueSamples[i][j] + valueSamples[i + 1][j] +
valueSamples[i][j + 1] + valueSamples[i + 1][j + 1];
sum += v * 0.25 * dr * dtheta;
area += dr * dtheta;
}
}
return sum / area;
}
static inline int getRSampleCount() {
return R_RES;
}
static inline double getRSample(int i) {
return R_SAMPLES[i];
}
static inline int getThetaSampleCount() {
return THETA_RES;
}
static inline double getThetaSample(int i) {
return THETA_SAMPLES[i];
}
inline double getValueSample(int ir, int itheta) const {
return valueSamples[ir][itheta];
}
inline void setValueSample(int ir, int itheta, double v) {
valueSamples[ir][itheta] = v;
}
static inline int getPeriodicThetaSampleCount() {
return 6 * THETA_RES;
}
static inline double getPeriodicThetaSample(int itheta) {
int part = itheta / THETA_RES;
int i0 = THETA_RES * (part % 2);
int sign = 1 - 2 * (part % 2);
double offset = (part + 1) / 2 * 2.0 * M_PI / 3.0;
return offset + sign * THETA_SAMPLES[i0 + sign * (itheta % THETA_RES)];
}
inline double getPeriodicValueSample(int ir, int itheta) const {
int part = itheta / THETA_RES;
int i0 = THETA_RES * (part % 2);
int sign = 1 - 2 * (part % 2);
return valueSamples[ir][i0 + sign * (itheta % THETA_RES)];
}
void dump(const char* file) const {
FILE* f = fopen(file, "w");
for (int i = 0; i <= R_RES; ++i) {
for (int j = 0; j <= THETA_RES; ++j) {
fprintf(f, "%d %d %.9g\n", i, j, valueSamples[i][j]);
}
fprintf(f, "\n");
}
fclose(f);
}
void plot(const char* file) const {
FILE *f = fopen(file, "w");
for (int i = 0; i <= R_RES; ++i) {
double r = getRSample(i);
for (int j = 0; j <= getPeriodicThetaSampleCount(); ++j) {
double theta = getPeriodicThetaSample(j);
fprintf(f, "%.9g %.9g %.12g\n", r, theta,
getPeriodicValueSample(i, j) / SoverTwoPiRs);
}
fprintf(f, "\n");
}
fclose(f);
}
void plotTemperature(const char* file) const {
FILE *f = fopen(file, "w");
for (int i = 0; i <= R_RES; ++i) {
double r = getRSample(i);
for (int j = 0; j <= getPeriodicThetaSampleCount(); ++j) {
double theta = getPeriodicThetaSample(j);
double T = pow(getPeriodicValueSample(i, j) * M_PI / sigma, 0.25);
fprintf(f, "%.9g %.9g %.12g\n", r, theta, T - 273);
}
fprintf(f, "\n");
}
fclose(f);
}
private:
static double* computeRSamples() {
double* samples = new double[R_RES + 1];
for (int i = 0; i <= R_RES; ++i) {
samples[i] = (1.0 - pow(1.0 - i / double(R_RES), R_POWER)) * rs;
}
return samples;
}
static double* computeThetaSamples() {
double* samples = new double[THETA_RES + 1];
for (int i = 0; i <= THETA_RES; ++i) {
samples[i] = pow(i / double(THETA_RES), THETA_POWER) * M_PI / 3.0;
}
return samples;
}
static const int R_RES = 64;
static const double R_POWER = 4.0;
static double* R_SAMPLES;
static const int THETA_RES = 64;
static const double THETA_POWER = 4.0;
static double* THETA_SAMPLES;
double valueSamples[R_RES + 1][THETA_RES + 1];
};
double* Jfunction::THETA_SAMPLES = Jfunction::computeThetaSamples();
double* Jfunction::R_SAMPLES = Jfunction::computeRSamples();
// -----------------------------------------------------------------------------
// f_n(x) functions
// -----------------------------------------------------------------------------
template<int n>
class f {
public:
double value(double x) {
double ix = x * FN_RES;
int i = static_cast<int>(floor(ix));
double di = ix - i;
double v0 = FN_SAMPLES[i < 0 ? 0 : (i > FN_RES ? FN_RES : i)];
double v1 = FN_SAMPLES[i + 1 < 0 ? 0 : (i + 1 > FN_RES ? FN_RES : i + 1)];
return v0 * (1.0 - di) + v1 * di;
}
private:
static double* computeFnSamples() {
double* samples = new double[FN_RES + 1];
for (int i = 0; i <= FN_RES; ++i) {
double x = i / double(FN_RES);
const int PHI_SAMPLES = 4096;
// Compute the integral in the numerator in value and the integral in the
// denominator in norm, using a single loop, and divide them at the end.
double value = 0.0;
double norm = 0.0;
for (int j = 0; j < PHI_SAMPLES; ++j) {
double phi = (j + 0.5) * M_PI / (2.0 * PHI_SAMPLES);
double cosPhi = cos(phi);
double cosPhiPowN = pow(cosPhi, n);
value += pow(x, 1.0 / cosPhi) * cosPhiPowN;
norm += cosPhiPowN;
}
samples[i] = value / norm;
}
return samples;
}
private:
static const int FN_RES = 256;
static double* FN_SAMPLES;
};
template<> double* f<0>::FN_SAMPLES = f<0>::computeFnSamples();
template<> double* f<1>::FN_SAMPLES = f<1>::computeFnSamples();
template<> double* f<2>::FN_SAMPLES = f<2>::computeFnSamples();
The constructor of the F1 class performs this linear combination, where COEFFICIENTS[i * (n + 1) + j] contains $F_{1,i}^j$:
// -----------------------------------------------------------------------------
// Surface irradiance operator F1
// -----------------------------------------------------------------------------
template<int K>
class F1 : public Efunction {
public:
F1(Efunction& E) {
int n = getThetaSampleCount();
for (int i = 0; i <= n; ++i) {
double v = 0.0;
for (int j = 0; j <= n; ++j) {
v += E.getValueSample(j) * COEFFICIENTS[i * (n + 1) + j];
}
setValueSample(i, v);
}
}
The precomputation of the coefficients is done below. We first try to load them
from disk from a previous execution. If they can't be found we compute them and
then store them on disk. The computation is done inside nested loops over the
$i$ and $j$ indices. For each $i,j$ indices, we must compute the integral
$\mathcal{F}_1[E^j](\theta_i)$. We thus start by initializing a function
base to the base function $E^j$.
private:
static double* computeCoefficients() {
int n = getThetaSampleCount();
int coeffCount = (n + 1) * (n + 1);
double* coefficients = new double[coeffCount];
char name[256];
sprintf(name, "coefficients/F1coefficients_%d.dat", K);
FILE* f = fopen(name, "rb");
if (f) {
fread(coefficients, sizeof(double), coeffCount, f);
fclose(f);
} else {
printf("Computing F1 coefficients...\n");
Efunction base;
for (int i = 0; i <= n; ++i) {
double thetap = getThetaSample(i);
for (int j = 0; j <= n; ++j) {
double v = 0.0;
base.setBaseFunction(j);
By definition, $E^j$ is equal to 0 almost everywhere, except between the samples $j-1$ and $j+1$ and between the intervals obtained by reflection over $\pi/3$ and by translation of $2\pi/3$ and $4\pi/3$ (recall that we store samples only over the $[0,\pi/3]$ interval). To save computations, we therefore define an integral method to compute the integral between two consecutive samples $\theta_i$ and $\theta_{i+1}$, and we call this function for each interval where we know the function is not null:
if (j < n) {
v += integral(thetap, base, j);
}
if (j % n > 0) {
v += integral(thetap, base, j - 1);
v += integral(thetap, base, 2 * n + j - 1);
v += integral(thetap, base, 2 * n + j);
v += integral(thetap, base, 4 * n + j - 1);
v += integral(thetap, base, 4 * n + j);
}
v += integral(thetap, base, 2 * n - j - 1);
v += integral(thetap, base, 2 * n - j);
v += integral(thetap, base, 4 * n - j - 1);
v += integral(thetap, base, 4 * n - j);
v += integral(thetap, base, 6 * n - j - 1);
if (j > 0) {
v += integral(thetap, base, 6 * n - j);
}
coefficients[i * (n + 1) + j] = v;
}
}
f = fopen(name, "wb");
fwrite(coefficients, sizeof(double), coeffCount, f);
fclose(f);
}
return coefficients;
}
We finally need to define the integral method. This method must compute
the integral defining the $\mathcal{F}_1$ operator, but only between $\theta_i$
and $\theta_{i+1}$:
\begin{equation*}
\mathrm{integral}(\theta_p,E,i)\eqdef
\frac{1}{4}\int_{\theta_i}^{\theta_{i+1}}E(\theta)
f_2(\ts(r_s,\theta_p,r_s,\theta))\vert\sin\frac{\theta-\theta_p}{2}\vert
\diff\theta
\end{equation*}
We compute it using the mid-point method:
static double integral(double thetap, const Efunction& E, int itheta) {
double theta0 = getPeriodicThetaSample(itheta);
double theta1 = getPeriodicThetaSample(itheta + 1);
double dtheta = theta1 - theta0;
int stheta = 1;
// Make sure to use small enough dtheta intervals to compute the integral.
while (dtheta > 1e-3) {
dtheta /= 2.0;
stheta *= 2;
}
double value0 = E.getPeriodicValueSample(itheta);
double value1 = E.getPeriodicValueSample(itheta + 1);
double dvalue = (value1 - value0) / stheta;
f<2> f2;
double result = 0.0;
for (int i = 0; i < stheta; ++i) {
double theta = theta0 + (i + 0.5) * dtheta;
double value = value0 + (i + 0.5) * dvalue;
double sinT = fabsf(sin((theta - thetap) * 0.5));
result += value * f2.value(t<K>(rs, thetap, rs, theta)) * sinT;
}
return result * dtheta / 4.0;
}
static double* COEFFICIENTS;
};
template<> double* F1<0>::COEFFICIENTS = F1<0>::computeCoefficients();
template<> double* F1<kS>::COEFFICIENTS = F1<kS>::computeCoefficients();
template<> double* F1<kE>::COEFFICIENTS = F1<kE>::computeCoefficients();
The constructor of the F2 class performs this linear combination, where COEFFICIENTS[[(i * (nr + 1) + j) * (nt + 1) + k] contains $F_{2,i}^{j,k}$:
// -----------------------------------------------------------------------------
// Volume irradiance operator F2
// -----------------------------------------------------------------------------
template<int K>
class F2 : public Efunction {
public:
F2(Jfunction& J) {
int n = getThetaSampleCount();
int nr = Jfunction::getRSampleCount();
int nt = Jfunction::getThetaSampleCount();
for (int i = 0; i <= n; ++i) {
double v = 0.0;
for (int j = 0; j <= nr; ++j) {
for (int k = 0; k <= nt; ++k) {
v += COEFFICIENTS[(i * (nr + 1) + j) * (nt + 1) + k] *
J.getValueSample(j, k);
}
}
setValueSample(i, v);
}
}
The precomputation of the coefficients is done below, using the same method as
for $\mathcal{F}_1$: inside nested loops over the $i,j,k$ indices, we start by
initializing a function base to the base function $J^{j,k}$:
private:
static double* computeCoefficients() {
int n = getThetaSampleCount();
int nr = Jfunction::getRSampleCount();
int nt = Jfunction::getThetaSampleCount();
int coeffCount = (n + 1) * (nr + 1) * (nt + 1);
double* coefficients = new double[coeffCount];
char name[256];
sprintf(name, "coefficients/F2coefficients_%d.dat", K);
FILE* f = fopen(name, "rb");
if (f) {
fread(coefficients, sizeof(double), coeffCount, f);
fclose(f);
} else {
printf("Computing F2 coefficients...\n");
Jfunction base;
for (int i = 0; i <= n; ++i) {
printf("%d / %d\n", i, n);
double thetap = getThetaSample(i);
for (int j = 0; j <= nr; ++j) {
for (int k = 0; k <= nt; ++k) {
double v = 0.0;
base.setBaseFunction(j, k);
we then compute the integral $\mathcal{F}_2[J^{j,k}](\theta_i)$ by decomposing
it into a sum of integrals over the 2D areas where we know the base function
$J^{j,k}$ is not null:
if (j > 0) {
if (k < nt) {
v += integral(thetap, base, j - 1, k);
}
if (k % nt > 0) {
v += integral(thetap, base, j - 1, k - 1);
v += integral(thetap, base, j - 1, 2 * nt + k - 1);
v += integral(thetap, base, j - 1, 2 * nt + k);
v += integral(thetap, base, j - 1, 4 * nt + k - 1);
v += integral(thetap, base, j - 1, 4 * nt + k);
}
v += integral(thetap, base, j - 1, 2 * nt - k - 1);
v += integral(thetap, base, j - 1, 2 * nt - k);
v += integral(thetap, base, j - 1, 4 * nt - k - 1);
v += integral(thetap, base, j - 1, 4 * nt - k);
v += integral(thetap, base, j - 1, 6 * nt - k - 1);
if (k > 0) {
v += integral(thetap, base, j - 1, 6 * nt - k);
}
}
if (j < nr) {
if (k < nt) {
v += integral(thetap, base, j, k);
}
if (k % nt > 0) {
v += integral(thetap, base, j, k - 1);
v += integral(thetap, base, j, 2 * nt + k - 1);
v += integral(thetap, base, j, 2 * nt + k);
v += integral(thetap, base, j, 4 * nt + k - 1);
v += integral(thetap, base, j, 4 * nt + k);
}
v += integral(thetap, base, j, 2 * nt - k - 1);
v += integral(thetap, base, j, 2 * nt - k);
v += integral(thetap, base, j, 4 * nt - k - 1);
v += integral(thetap, base, j, 4 * nt - k);
v += integral(thetap, base, j, 6 * nt - k - 1);
if (k > 0) {
v += integral(thetap, base, j, 6 * nt - k);
}
}
coefficients[(i * (nr + 1) + j) * (nt + 1) + k] = v;
}
}
}
f = fopen(name, "wb");
fwrite(coefficients, sizeof(double), coeffCount, f);
fclose(f);
}
return coefficients;
}
We finally need to define the integral method. This method must compute
the integral defining the $\mathcal{F}_2$ operator, but only between $r_i$ and
$r_{i+1}$ and $\theta_j$ and $\theta_{j+1}$:
\begin{equation*}
\mathrm{integral}(\theta_p,J,i,j)\eqdef
2\int_{\theta_j}^{\theta_{j+1}}\int_{r_i}^{r_{i+1}}k_s(r)J(r,\theta)
f_1(\ts(r_s,\theta_p,r,\theta))\frac{r(r_s-r\cos(\theta-\theta_p))}
{r_s^2+r^2-2rr_s\cos(\theta-\theta_p)}\diff r\diff\theta
\end{equation*}
We compute it using the mid-point method:
static double integral(double thetap, const Jfunction& J, int ir,
int itheta) {
double r0 = J.getRSample(ir);
double r1 = J.getRSample(ir + 1);
double dr = r1 - r0;
int sr = 1;
// Make sure to use small enough dr intervals to compute the integral.
while (sr > 10 && sr < (1 << 20)) {
dr /= 2.0;
sr *= 2;
}
double theta0 = J.getPeriodicThetaSample(itheta);
double theta1 = J.getPeriodicThetaSample(itheta + 1);
double dtheta = theta1 - theta0;
int stheta = 1;
// Make sure to use small enough dtheta intervals to compute the integral.
while (dtheta > 1e-3 && stheta < (1 << 20)) {
dtheta /= 2.0;
stheta *= 2;
}
double value00 = J.getPeriodicValueSample(ir, itheta);
double value01 = J.getPeriodicValueSample(ir, itheta + 1);
double value10 = J.getPeriodicValueSample(ir + 1, itheta);
double value11 = J.getPeriodicValueSample(ir + 1, itheta + 1);
f<1> f1;
double result = 0.0;
for (int i = 0; i < sr; ++i) {
double r = r0 + (i + 0.5) * dr;
double ix = (i + 0.5) / sr;
double value0 = value00 * (1.0 - ix) + value10 * ix;
double value1 = value01 * (1.0 - ix) + value11 * ix;
for (int j = 0; j < stheta; ++j) {
double theta = theta0 + (j + 0.5) * dtheta;
double jx = (j + 0.5) / stheta;
double value = value0 * (1.0 - jx) + value1 * jx;
double rCosT = r * cos(theta - thetap);
result += value * k<K>(r) * f1.value(t<K>(rs, thetap, r, theta)) *
r * (rs - rCosT) / (rs * rs + r * r - 2.0 * rs * rCosT);
}
}
return result * 2.0 * dr * dtheta;
}
static double* COEFFICIENTS;
};
template<> double* F2<kS>::COEFFICIENTS = F2<kS>::computeCoefficients();
template<> double* F2<kE>::COEFFICIENTS = F2<kE>::computeCoefficients();
The constructor of the F3 class performs this linear combination, where COEFFICIENTS[[(i * (nt + 1) + j) * (n + 1) + k] contains $F_{3,i,j}^k$:
// -----------------------------------------------------------------------------
// Average surface radiance operator F3
// -----------------------------------------------------------------------------
template<int K>
class F3 : public Jfunction {
public:
F3(Efunction& E) {
int nr = getRSampleCount();
int nt = getThetaSampleCount();
int n = Efunction::getThetaSampleCount();
for (int i = 0; i <= nr; ++i) {
for (int j = 0; j <= nt; ++j) {
double v = 0.0;
for (int k = 0; k <= n; ++k) {
v += COEFFICIENTS[(i * (nt + 1) + j) * (n + 1) + k] *
E.getValueSample(k);
}
setValueSample(i, j, v);
}
}
}
The precomputation of the coefficients is done below, using the same method as
for $\mathcal{F}_1$ and $\mathcal{F}_2$: inside nested loops over the $i,j,k$
indices, we start by initializing a function base to the base function
$E^k$:
private:
static double* computeCoefficients() {
int nr = getRSampleCount();
int nt = getThetaSampleCount();
int n = Efunction::getThetaSampleCount();
int coeffCount = (nr + 1) * (nt + 1) * (n + 1);
double* coefficients = new double[coeffCount];
char name[256];
sprintf(name, "coefficients/F3coefficients_%d.dat", K);
FILE* f = fopen(name, "rb");
if (f) {
fread(coefficients, sizeof(double), coeffCount, f);
fclose(f);
} else {
printf("Computing F3 coefficients...\n");
Efunction base;
for (int i = 0; i <= nr; ++i) {
printf("%d / %d\n", i, nr);
double rp = std::min(getRSample(i), rs - 0.0001);
for (int j = 0; j <= nt; ++j) {
double thetap = getThetaSample(j);
for (int k = 0; k <= n; ++k) {
double v = 0.0;
base.setBaseFunction(k);
we then compute the integral $\mathcal{F}_3[E^k](r_i,\theta_j)$ by decomposing
it into a sum of integrals over the intervals where we know the base function
$E^k$ is not null:
if (k < n) {
v += integral(rp, thetap, base, k);
}
if (k % n > 0) {
v += integral(rp, thetap, base, k - 1);
v += integral(rp, thetap, base, 2 * n + k - 1);
v += integral(rp, thetap, base, 2 * n + k);
v += integral(rp, thetap, base, 4 * n + k - 1);
v += integral(rp, thetap, base, 4 * n + k);
}
v += integral(rp, thetap, base, 2 * n - k - 1);
v += integral(rp, thetap, base, 2 * n - k);
v += integral(rp, thetap, base, 4 * n - k - 1);
v += integral(rp, thetap, base, 4 * n - k);
v += integral(rp, thetap, base, 6 * n - k - 1);
if (k > 0) {
v += integral(rp, thetap, base, 6 * n - k);
}
coefficients[(i * (nt + 1) + j) * (n + 1) + k] = v;
}
}
}
f = fopen(name, "wb");
fwrite(coefficients, sizeof(double), coeffCount, f);
fclose(f);
}
return coefficients;
}
We finally need to define the integral method. This method must compute
the integral defining the $\mathcal{F}_3$ operator, but only between $\theta_i$
and $\theta_{i+1}$:
\begin{equation*}
\mathrm{integral}(r_p,\theta_p,E,i)\eqdef
\frac{1}{2\pi^2}\int_{\theta_i}^{\theta_{i+1}}E(\theta)
f_1(\ts(r_p,\theta_p,r_s,\theta))
\frac{r_s(r_s-r_p\cos(\theta-\theta_p))}
{r_s^2+r_p^2-2r_pr_s\cos(\theta-\theta_p)}\diff\theta
\end{equation*}
We compute it using the mid-point method, with one subtelty: in order to increase the precision of the result without requiring too much samples, i.e. very long computation times, we use an analytic integral of the fraction term. This means that we approximate the integral over a small sub interval $[\theta_a,\theta_b]$ of $[\theta_i,\theta_{i+1}]$ with: \begin{align*} \int_{\theta_a}^{\theta_b}E(\theta)&f_1(\ts(r_p,\theta_p,r_s,\theta)) \frac{r_s(r_s-r_p\cos(\theta-\theta_p))} {r_s^2+r_p^2-2r_pr_s\cos(\theta-\theta_p)}\diff\theta\\ &\approx E(\frac{\theta_a+\theta_b}{2}) f_1(\ts(r_p,\theta_p,r_s,\frac{\theta_a+\theta_b}{2})) [F(r_p,\theta_b-\theta_p)-F(r_p,\theta_a-\theta_p)] \end{align*} where $F$ is the indefinite integral of the fraction term with respect to $\theta$: \begin{equation*} F(r,\theta)\eqdef \int\frac{r_s(r_s-r\cos\theta)}{r_s^2+r^2-2rr_s\cos\theta}\diff\theta= \tan^{-1}\left(\frac{r-r_s}{r+r_s}\cot\frac{\theta}{2}\right)+ \frac{\theta}{2}+\pi\left\lfloor\frac{\theta}{2\pi}\right\rfloor \end{equation*}
static double integral(double rp, double thetap, const Efunction& E,
int itheta) {
double theta0 = E.getPeriodicThetaSample(itheta);
double theta1 = E.getPeriodicThetaSample(itheta + 1);
double dtheta = theta1 - theta0;
double dthetaMax = rs - rp < 1 ?
(thetap == theta0 || thetap == theta1 ? 1e-5 : 1e-4) : 1e-3;
int stheta = 1;
// Make sure to use small enough dtheta intervals to compute the integral.
while (dtheta > dthetaMax && stheta < (1 << 20)) {
dtheta /= 2.0;
stheta *= 2;
}
double value0 = E.getPeriodicValueSample(itheta);
double value1 = E.getPeriodicValueSample(itheta + 1);
double dvalue = (value1 - value0) / stheta;
f<1> f1;
double rp_ri = (rp - rs) / (rp + rs);
double result = 0.0;
for (int i = 0; i < stheta; ++i) {
double theta = theta0 + (i + 0.5) * dtheta;
double value = value0 + (i + 0.5) * dvalue;
double thetaA = (theta0 - thetap + (i + 1e-9) * dtheta) / 2.0;
double thetaB = (theta0 - thetap + (i + 1 - 1e-9) * dtheta) / 2.0;
double Fa = atan(rp_ri / tan(thetaA)) + thetaA;
double Fb = atan(rp_ri / tan(thetaB)) + thetaB;
double rInt = Fb - Fa + (thetaA * thetaB <= 0.0 ? M_PI : 0.0);
result += value * f1.value(t<K>(rp, thetap, rs, theta)) * rInt;
}
return result / (2.0 * M_PI * M_PI);
}
static double* COEFFICIENTS;
};
template<> double* F3<kS>::COEFFICIENTS = F3<kS>::computeCoefficients();
template<> double* F3<kE>::COEFFICIENTS = F3<kE>::computeCoefficients();
The constructor of the F4 class performs this linear combination, where COEFFICIENTS[((i * (nt + 1) + j) * (nr + 1) + k) * (nt + 1) + l] contains $F_{4,i,j}^{k,l}$:
// -----------------------------------------------------------------------------
// Average volume radiance functional F4.
// -----------------------------------------------------------------------------
template<int K>
class F4 : public Jfunction {
public:
F4(Jfunction& J) {
int nr = getRSampleCount();
int nt = getThetaSampleCount();
for (int i = 0; i <= nr; ++i) {
for (int j = 0; j <= nt; ++j) {
double v = 0.0;
for (int k = 0; k <= nr; ++k) {
for (int l = 0; l <= nt; ++l) {
int index = ((i * (nt + 1) + j) * (nr + 1) + k) * (nt + 1) + l;
v += COEFFICIENTS[index] * J.getValueSample(k, l);
}
}
setValueSample(i, j, v);
}
}
}
The precomputation of the coefficients is done below, using the same method as
above: inside nested loops over the $i,j,k,l$ indices, we start by initializing
a function base to the base function $J^{k,l}$:
private:
static double* computeCoefficients() {
int nr = getRSampleCount();
int nt = getThetaSampleCount();
int coeffCount = (nr + 1) * (nt + 1) * (nr + 1) * (nt + 1);
double* coefficients = new double[coeffCount];
char name[256];
sprintf(name, "coefficients/F4coefficients_%d.dat", K);
FILE* f = fopen(name, "rb");
if (f) {
fread(coefficients, sizeof(double), coeffCount, f);
fclose(f);
} else {
printf("Computing F4 coefficients...\n");
Jfunction base;
for (int i = 0; i <= nr; ++i) {
double rp = std::min(getRSample(i), rs - 0.0001);
for (int j = 0; j <= nt; ++j) {
double thetap = getThetaSample(j);
printf("%d %d / %d %d\n", i, j, nr, nt);
for (int k = 0; k <= nr; ++k) {
for (int l = 0; l <= nt; ++l) {
double v = 0.0;
base.setBaseFunction(k, l);
we then compute the integral $\mathcal{F}_4[J^{k,l}](r_i,\theta_j)$ by
decomposing it into a sum of integrals over the 2D areas where we know the base
function $J^{k,l}$ is not null:
if (k > 0) {
if (l < nt) {
v += integral(rp, thetap, base, k - 1, l);
}
if (l % nt > 0) {
v += integral(rp, thetap, base, k - 1, l - 1);
v += integral(rp, thetap, base, k - 1, 2 * nt + l - 1);
v += integral(rp, thetap, base, k - 1, 2 * nt + l);
v += integral(rp, thetap, base, k - 1, 4 * nt + l - 1);
v += integral(rp, thetap, base, k - 1, 4 * nt + l);
}
v += integral(rp, thetap, base, k - 1, 2 * nt - l - 1);
v += integral(rp, thetap, base, k - 1, 2 * nt - l);
v += integral(rp, thetap, base, k - 1, 4 * nt - l - 1);
v += integral(rp, thetap, base, k - 1, 4 * nt - l);
v += integral(rp, thetap, base, k - 1, 6 * nt - l- 1);
if (l > 0) {
v += integral(rp, thetap, base, k - 1, 6 * nt - l);
}
}
if (k < nr) {
if (l < nt) {
v += integral(rp, thetap, base, k, l);
}
if (l % nt > 0) {
v += integral(rp, thetap, base, k, l - 1);
v += integral(rp, thetap, base, k, 2 * nt + l - 1);
v += integral(rp, thetap, base, k, 2 * nt + l);
v += integral(rp, thetap, base, k, 4 * nt + l - 1);
v += integral(rp, thetap, base, k, 4 * nt + l);
}
v += integral(rp, thetap, base, k, 2 * nt - l - 1);
v += integral(rp, thetap, base, k, 2 * nt - l);
v += integral(rp, thetap, base, k, 4 * nt - l - 1);
v += integral(rp, thetap, base, k, 4 * nt - l);
v += integral(rp, thetap, base, k, 6 * nt - l- 1);
if (l > 0) {
v += integral(rp, thetap, base, k, 6 * nt - l);
}
}
int index = ((i * (nt + 1) + j) * (nr + 1) + k) * (nt + 1) + l;
coefficients[index] = v;
}
}
}
}
f = fopen(name, "wb");
fwrite(coefficients, sizeof(double), coeffCount, f);
fclose(f);
}
return coefficients;
}
We finally need to define the integral method. This method must compute
the integral defining the $\mathcal{F}_4$ operator, but only between $r_i$ and
$r_{i+1}$ and $\theta_j$ and $\theta_{j+1}$:
\begin{align*}
&\mathrm{integral}(r_p,\theta_p,J,i,j)\\
&\quad\eqdef\frac{1}{4}\int_{\theta_j}^{\theta_{j+1}}
\int_{r_i}^{r_{i+1}}k_s(r)J(r,\theta)f_0(\ts(r_p,\theta_p,r,\theta))
\frac{r}{\sqrt{r_p^2+r^2-2rr_p\cos(\theta-\theta_p)}}\diff r\diff\theta
\end{align*}
We compute it using the mid-point method, with one subtelty: in order to increase the precision of the result without requiring too much samples, i.e. very long computation times, we use an analytic integral of the fraction term. This means that we approximate the integral over a small sub area $[r_a,r_b]\times[\theta_a,\theta_b]$ of $[r_i,r_{i+1}]\times[\theta_j,\theta_{j+1}]$ with: \begin{align*} \int_{\theta_a}^{\theta_b}\int_{r_a}^{r_b}&k_s(r)J(r,\theta) f_0(\ts(r_p,\theta_p,r,\theta))\frac{r} {\sqrt{r_p^2+r^2-2rr_p\cos(\theta-\theta_p)}}\diff r\diff\theta\\ &\approx(\theta_b-\theta_a)k_s(r_m)J(r_m,\theta_m) f_0(\ts(r_p,\theta_p,r_m,\theta_m))\\ &\quad\quad\quad\quad\quad[G(r_b,\cos(\theta_m-\theta_p)) -G(r_a,\cos(\theta_m-\theta_p))] \end{align*} where $r_m=(r_a+r_b)/2$, $\theta_m=(\theta_a+\theta_b)/2$ and where $G$ is the indefinite integral of the fraction term with respect to $r$: \begin{align*} G(r,x)&\eqdef \int \frac{r}{\sqrt{r_p^2+r^2-2rr_px)}}\diff r\\ &=\sqrt{r_p^2+r^2-2rr_px}+r_px\ln\left(r+r_px+\sqrt{r_p^2+r^2-2rr_px}\right) \end{align*}
static double integral(double rp, double thetap, const Jfunction& J, int ir,
int itheta) {
double r0 = J.getRSample(ir);
double r1 = J.getRSample(ir + 1);
double dr = r1 - r0;
int sr = 1;
// Make sure to use small enough dr intervals to compute the integral.
while (dr > 100 && sr < (1 << 20)) {
dr /= 2.0;
sr *= 2;
}
double theta0 = J.getPeriodicThetaSample(itheta);
double theta1 = J.getPeriodicThetaSample(itheta + 1);
double dtheta = theta1 - theta0;
int stheta = 1;
// Make sure to use small enough dtheta intervals to compute the integral.
while (dtheta > 1e-3 && stheta < (1 << 20)) {
dtheta /= 2.0;
stheta *= 2;
}
double value00 = J.getPeriodicValueSample(ir, itheta);
double value01 = J.getPeriodicValueSample(ir, itheta + 1);
double value10 = J.getPeriodicValueSample(ir + 1, itheta);
double value11 = J.getPeriodicValueSample(ir + 1, itheta + 1);
f<0> f0;
double result = 0.0;
for (int i = 0; i < sr; ++i) {
double r = r0 + (i + 0.5) * dr;
double rA = r0 + (i + 1e-9) * dr;
double rB = r0 + (i + 1 - 1e-9) * dr;
double ix = (i + 0.5) / sr;
double value0 = value00 * (1.0 - ix) + value10 * ix;
double value1 = value01 * (1.0 - ix) + value11 * ix;
for (int j = 0; j < stheta; ++j) {
double theta = theta0 + (j + 0.5) * dtheta;
double jx = (j + 0.5) / stheta;
double value = value0 * (1.0 - jx) + value1 * jx;
double rpCosT = rp * cos(theta - thetap);
double Ga = sqrt(rp * rp + rA * rA - 2.0 * rA * rpCosT);
double Gb = sqrt(rp * rp + rB * rB - 2.0 * rB * rpCosT);
Ga += rpCosT * log(Ga - rpCosT + rA);
Gb += rpCosT * log(Gb - rpCosT + rB);
result += value * k<K>(r) * f0.value(t<K>(rp, thetap, r, theta)) *
(Gb - Ga);
}
}
return result * dtheta / 4.0;
}
static double* COEFFICIENTS;
};
template<> double* F4<kS>::COEFFICIENTS = F4<kS>::computeCoefficients();
template<> double* F4<kE>::COEFFICIENTS = F4<kE>::computeCoefficients();
// -----------------------------------------------------------------------------
// Initial received ground irradiance function E0
// -----------------------------------------------------------------------------
template<int K>
class E0 : public Efunction {
public:
E0() {
f<2> f2;
for (int i = 0; i <= getThetaSampleCount(); ++i) {
double theta0 = getThetaSample(i);
double theta1 = theta0 + 2.0 * M_PI / 3.0;
double theta2 = theta0 - 2.0 * M_PI / 3.0;
double v0 = SoverTwoPiRs * M_PI / 2.0 *
f2.value(t<K>(rs, theta0, rs, 0.0)) * fabsf(sin(theta0 / 2.0));
double v1 = SoverTwoPiRs * M_PI / 2.0 *
f2.value(t<K>(rs, theta1, rs, 0.0)) * fabsf(sin(theta1 / 2.0));
double v2 = SoverTwoPiRs * M_PI / 2.0 *
f2.value(t<K>(rs, theta2, rs, 0.0)) * fabsf(sin(theta2 / 2.0));
setValueSample(i, (v0 + v1 + v2) / 3.0);
}
}
};
The average radiance $J_0$ received directly from a single linear light source
is equal to
\begin{equation*}
J_0(r,\theta)=\frac{S}{2\pi^2 r_s\Delta\theta}f_1(\ts(r,\theta,r_s,0))
\big(F(r, -\theta+\frac{\Delta\theta}{2})
-F(r, -\theta-\frac{\Delta\theta}{2})\big)
\end{equation*}
where $F$ is the indefinite integral defined above. Its
implementation as a Jfunction subclass, for 3 linear light sources 120°
apart from each other is straightforward:
// -----------------------------------------------------------------------------
// Initial average radiance function J0
// -----------------------------------------------------------------------------
template<int K>
class J0 : public Jfunction {
public:
J0() {
f<1> f1;
for (int i = 0; i <= getRSampleCount(); ++i) {
double r = std::min(getRSample(i), rs - 0.0001);
for (int j = 0; j <= getThetaSampleCount(); ++j) {
double theta0 = getThetaSample(j);
double theta1 = theta0 + 2.0 * M_PI / 3.0;
double theta2 = theta0 - 2.0 * M_PI / 3.0;
double v0 = SoverTwoPiRs / M_PI *
f1.value(t<K>(r, theta0, rs, 0.0)) * helper(r, theta0);
double v1 = SoverTwoPiRs / M_PI *
f1.value(t<K>(r, theta1, rs, 0.0)) * helper(r, theta1);
double v2 = SoverTwoPiRs / M_PI *
f1.value(t<K>(r, theta2, rs, 0.0)) * helper(r, theta2);
setValueSample(i, j, (v0 + v1 + v2) / 3.0);
}
}
}
private:
static double helper(double r, double theta) {
double a = -DeltaTheta / 2.0 - theta;
double b = DeltaTheta / 2.0 - theta;
double Fa = atan((r - rs) / ((r + rs) * tan(a / 2))) + a / 2;
double Fb = atan((r - rs) / ((r + rs) * tan(b / 2))) + b / 2;
return (Fb - Fa + (a * b < 0.0 && r < rs ? M_PI : 0.0)) / DeltaTheta;
}
};
We saw such a case in the Atmosphere Temperature and Pressure section: in the axisymetric case, the ground and atmosphere temperatures are equal to each other and are the same everywhere, and this is true for any scattering and emissivity / absorptivity coefficients. This means $\bar{J}=\bar{E}^{\uparrow}/\pi$. Also, in this case the received and emitted radiances in the long wave range are equal everywhere, as we saw in the Average Lighting and Temperature section. This means $\bar{E}^{\uparrow}=\bar{E}^{\downarrow}$. Reporting this in the equations we found in the Temperature Model section: \begin{align*} \bar{E}^{\uparrow}&=(1-\alpha)\aBar E+\aBar \bar{E}^{\downarrow}\\ \bar{E}^{\downarrow}&= \cF_1[(1-\alpha)\aBar E]+\aBar\cF_1[\bar{E}^{\downarrow}]+\cF_2[\bar{J}]\\ \bar{J}&= \cF_3[(1-\alpha)\aBar E]+\aBar\cF_3[\bar{E}^{\downarrow}]+\cF_4[\bar{J}] \end{align*} we get $(1-\alpha)\aBar E=(1-\aBar)\bar{E}^{\downarrow}$ and from this: \begin{align*} \bar{E}^{\downarrow}&= \cF_1[\bar{E}^{\downarrow}]+\cF_2[\bar{E}^{\downarrow}/\pi]\\ \bar{E}^{\downarrow}/\pi&= \cF_3[\bar{E}^{\downarrow}]+\cF_4[\bar{E}^{\downarrow}/\pi] \end{align*}
In other words, for $E(\theta)=\pi$ and $J(r,\theta)=1$, for instance, we should have $(\cF_1[E]+\cF_2[J])(\theta)/\pi=1$ and $(\cF_3[E]+\cF_4[J])(r,\theta)=1$. This is what we check below, by using the classes defined above.
template<int K>
void validation() {
Efunction e0;
Jfunction j0;
e0.set(M_PI);
j0.set(1.0);
Efunction e1a = F1<K>(e0);
Efunction e1b = F2<K>(j0);
Jfunction j1a = F3<K>(e0);
Jfunction j1b = F4<K>(j0);
char name[256];
sprintf(name, "testF1-%d.txt", K);
e1a.dump(name);
sprintf(name, "testF2-%d.txt", K);
e1b.dump(name);
sprintf(name, "testF3-%d.txt", K);
j1a.dump(name);
sprintf(name, "testF4-%d.txt", K);
j1b.dump(name);
e1a.add(e1b).mult(1.0 / M_PI);
j1a.add(j1b);
sprintf(name, "testF1+F2-%d.txt", K);
e1a.dump(name);
sprintf(name, "testF3+F4-%d.txt", K);
j1a.dump(name);
}
The results are shown in Fig. 2, Fig. 3
and Fig. 4. As can be seen the results are very close to the
constant function 1, with an error which is less than $10^{-4}$ for $k=k_s$ and
about $10^{-3}$ for $k=k_e$. We can also see that the results are always below
1, instead of been centered around 1. Finally, note that the shape of the
functions comes from our adaptive number of samples to compute each integral:
when we use 2 times more samples to get sufficiently small $r$ or $\theta$
intervals, the resulting error suddenly jumps.
int main(int argc, char* argv[]) {
validation<kS>();
validation<kE>();
We then compute $(\cF_1[\pi]+\cF_2[1])/\pi$ and $\cF_3[\pi]+\cF_4[1]$, for
$k=k_s$. We will use them to "normalize" the results below, in order to reduce
approximation errors.
Efunction eCst; Jfunction jCst; eCst.set(M_PI); jCst.set(1.0); Efunction eRef = F1<kS>(eCst).add(F2<kS>(jCst)).mult(1.0 / M_PI); Jfunction jRef = F3<kS>(eCst).add(F4<kS>(jCst));The next step is to initialize the $E$ and $J$ functions we want to compute to the initial functions $E_0$ and $J_0$:
Efunction e; Jfunction j; Efunction eN = E0<kS>(); Jfunction jN = J0<kS>(); e.set(eN); j.set(jN);We can then compute each term of the series $E_{n+1}=\alpha\cF_1[E_n]+cF_2[J_n]$ and $J_{n+1}=\alpha\cF_3[E_n]+\cF_4[J_n]$, and add them to $E$ and $J$ (6 terms are enough). In order to cancel, or at least reduce, the approximation errors (which include a small systematic bias as we saw in the validation above), we "normalize" $E_{n+1}$ and $J_{n+1}$ by dividing them by $(\cF_1[\pi]+\cF_2[1])/\pi$ and $\cF_3[\pi]+\cF_4[1]$, before adding them to $E$ and $J$:
for (int i = 1; i <= 6; ++i) {
Efunction eN1 = F1<kS>(eN).mult(alpha).add(F2<kS>(jN));
Jfunction jN1 = F3<kS>(eN).mult(alpha).add(F4<kS>(jN));
eN.set(eN1.div(eRef));
jN.set(jN1.div(jRef));
e.add(eN);
j.add(jN);
printf("lighting iteration %d\n", i);
printf("E avg (actual / limit): %g\n",
e.average() / (SoverTwoPiRs / (1.0 - alpha)));
printf("J avg (actual / limit): %g\n",
j.average() / (SoverTwoPiRs / (1.0 - alpha) / M_PI));
}
At this point we have the ground irradiance and atmosphere radiance in the short
wave range, and we save them to disk:
e.plot("e.txt");
j.plot("j.txt");
The next step is to compute the same quantities in the long wave range. Again,
we start by computing $(\cF_1[\pi]+\cF_2[1])/\pi$ and $\cF_3[\pi]+\cF_4[1]$, for
$k=k_e$. We will use them to "normalize" the results below, in order to reduce
approximation errors.
eRef = F1<kE>(eCst).add(F2<kE>(jCst)).mult(1.0 / M_PI); jRef = F3<kE>(eCst).add(F4<kE>(jCst));Then we initialize the $\bar{E}^{\downarrow}$ and $\bar{J}$ functions we want to compute to the initial functions $\cF_1[(1-\alpha)\aBar E]$ and $\cF_3[(1-\alpha)\aBar E]$:
Efunction eBarDown; Jfunction jBar; Efunction eBarDownN = F1<kE>(e).mult((1.0 - alpha) * alphaBar); Jfunction jBarN = F3<kE>(e).mult((1.0 - alpha) * alphaBar); eBarDown.set(eBarDownN); jBar.set(jBarN);We can then compute each term of the series $\bar{E}^{\downarrow}_{n+1}=\aBar\cF_1[\bar{E}^{\downarrow}_n]+\cF_2[\bar{J}_n]$ and $\bar{J}_{n+1}=\aBar\cF_3[\bar{E}^{\downarrow}_n]+\cF_4[\bar{J}_n]$, and add them to $\bar{E}^{\downarrow}$ and $\bar{J}$. Note that we need much more terms for the long wave range than for the short wave range: 100 terms are needed instead of 6. One reason is that $\aBar$ is larger than $\alpha$. Thus the terms of the series decrease less rapidly than in the short wave range. The other reason is that $k_e$ is larger than $k_s$. This means that the mean free path is shorter or, in other words, that it takes many bounces for long wave radiations emitted at one side of Rama to reach the other side. And since each term of the series corresponds to an additional bounce, we need many terms.
Due this large number of terms, it is even more important than above to cancel, or at least reduce, the approximation errors (which include a small systematic bias as we saw in the validation above) in each term. For this we "normalize" $\bar{E}^{\downarrow}_{n+1}$ and $\bar{J}_{n+1}$ by dividing them by $(\cF_1[\pi]+\cF_2[1])/\pi$ and $\cF_3[\pi]+\cF_4[1]$, before adding them to $\bar{E}^{\downarrow}$ and $\bar{J}$.
for (int i = 1; i <= 100; ++i) {
Efunction eBarDownN1 = F1<kE>(eBarDownN).mult(alphaBar).add(F2<kE>(jBarN));
Jfunction jBarN1 = F3<kE>(eBarDownN).mult(alphaBar).add(F4<kE>(jBarN));
eBarDownN.set(eBarDownN1.div(eRef));
jBarN.set(jBarN1.div(jRef));
eBarDown.add(eBarDownN);
jBar.add(jBarN);
printf("temperature iteration %d\n", i);
printf("Ebardown avg (actual / limit): %g\n",
eBarDown.average() / (Delta * SoverTwoPiRs));
printf("Jbar avg (actual / limit): %g\n",
jBar.average() / (Delta * SoverTwoPiRs / M_PI));
}
Even with 100 terms, we are still missing the (small) contribution of all the
remaining terms of the series. These terms can be approximated with constant
functions, and we know the theoritical average value of $\bar{E}^{\downarrow}$
and $\bar{J}$, so we account for them by adding to our result the difference
between the theoritical average value and the actual one:
eBarDown.add(Delta * SoverTwoPiRs - eBarDown.average()); jBar.add(Delta * SoverTwoPiRs / M_PI - jBar.average());At this point we have the ground irradiance and atmosphere radiance in the long wave range, and we save them to disk (as well as the corresponding temperatures):
eBarDown.plot("eBarDown.txt");
jBar.plot("jBar.txt");
jBar.plotTemperature("Ta.txt");
// J + Jbar, and the corresponding temperature, gives the temperature a black
// body would have in the atmosphere, depending on its position in Rama.
j.add(jBar).plot("jTotal.txt");
j.plotTemperature("Tt.txt");
Efunction eUp;
eUp.set(eBarDown);
eUp.mult(alphaBar).add(e.mult((1.0 - alpha) * alphaBar));
eUp.plot("eBarUp.txt");
eUp.plotTemperature("Ts.txt");
printf("Eup avg (actual / limit): %g\n",
eUp.average() / (Delta * SoverTwoPiRs));
return 0;
}