26 #include "NumericalApproximation.h"
75 double trap[maxits+1];
86 for(
int i = 0; i < maxits; i++) {
89 trap[i] = RefineExtendedTrap(am, sub, a, b, trap[i], i + 1);
91 interp.
AddData(5, &h[i-4], &trap[i-4]);
92 ss = interp.
Evaluate(0.0, NumericalApproximation::Extrapolate);
96 if(fabs(dss) <= epsilon * fabs(ss))
return ss;
97 if(fabs(dss) <= epsilon2)
return ss;
100 h[i+1] = 0.25 * h[i];
111 "NumericalAtmosApprox::RombergsMethod() - Caught the following error: ",
115 "NumericalAtmosApprox::RombergsMethod() - Failed to converge in "
116 +
IString(maxits) +
" iterations.",
152 double NumericalAtmosApprox::RefineExtendedTrap(
AtmosModel *am,
IntegFunc sub,
double a,
double b,
double s,
unsigned int n) {
160 if(sub == InnerFunction) {
161 begin = InrFunc2Bint(am, a);
162 end = InrFunc2Bint(am, b);
165 begin = OutrFunc2Bint(am, a);
166 end = OutrFunc2Bint(am, b);
168 return (0.5 * (b - a) * (begin + end));
172 double delta, tnm, x, sum;
173 it = (int)(pow(2.0, (
double)(n - 2)));
175 delta = (b - a) / tnm;
178 for(
int i = 0; i < it; i++) {
179 if(sub == InnerFunction) {
180 sum = sum + InrFunc2Bint(am, x);
183 sum = sum + OutrFunc2Bint(am, x);
187 return (0.5 * (s + (b - a) * sum / tnm));
193 "NumericalAtmosApprox::RefineExtendedTrap() - Caught the following error: ",
217 double NumericalAtmosApprox::OutrFunc2Bint(
AtmosModel *am,
double phi) {
220 sub = NumericalAtmosApprox::InnerFunction;
222 am->p_atmosPhi = phi;
223 am->p_atmosCosphi = cos((
PI / 180.0) * phi);
233 "NumericalAtmosApprox::OutrFunc2Bint() - Caught the following error: ",
266 double NumericalAtmosApprox::InrFunc2Bint(
AtmosModel *am,
double mu) {
280 ema = acos(mu) * (180.0 /
PI);
281 sine = sin(ema * (
PI / 180.0));
282 if((am->p_atmosAtmSwitch == 0) || (am->p_atmosAtmSwitch == 2)) {
283 alpha = am->p_atmosSini * sine * am->p_atmosCosphi +
284 am->p_atmosMunot * mu;
287 alpha = am->p_atmosSini * sine * am->p_atmosCosphi -
288 am->p_atmosMunot * mu;
290 phase = acos(alpha) * (180.0 /
PI);
293 if(am->p_atmosAtmSwitch == 0) {
295 result = mu * am->p_atmosPM->
CalcSurfAlbedo(phase, am->p_atmosInc, ema);
299 xx = -am->
AtmosTau() / max(am->p_atmosMunot, 1.0e-30);
309 xx = -am->
AtmosTau() / max(mu, 1.0e-30);
320 if(mu == am->p_atmosMunot) {
321 tfac = am->
AtmosTau() * emunot / (am->p_atmosMunot * am->p_atmosMunot);
324 tfac = (emunot - emu) / (am->p_atmosMunot - mu);
326 if(am->p_atmosAtmSwitch == 1) {
327 result = mu * (phasefn - 1.0) * tfac;
329 else if(am->p_atmosAtmSwitch == 2) {
330 result = am->p_atmosMunot * mu * (phasefn - 1.0) * (1.0 - emunot * emu) / (am->p_atmosMunot + mu);
332 else if(am->p_atmosAtmSwitch == 3) {
333 result = -sine * am->p_atmosCosphi * (phasefn - 1.0) * tfac;
336 string msg =
"NumericalAtmosApprox::InrFunc2Bint() - Invalid value of atmospheric ";
337 msg +=
"switch used as argument to this function";