← Back to branch summary

~ubuntu-branches/ubuntu/trusty/xplanet/trusty

« back to all changes in this revision

Viewing changes to src/libmultiple/RayleighScattering.cpp

  • Committer: Package Import Robot
  • Author(s): Steve McIntyre
  • Date: 2014-01-30 18:54:10 UTC
  • mfrom: (4.1.8 sid)
  • Revision ID: package-import@ubuntu.com-20140130185410-o2u51emh707gbevv
Tags: 1.3.0-1
http://bugs.debian.org/678464
http://bugs.debian.org/736052
http://bugs.debian.org/662568
* New maintainer
* Move to new upstream version 1.3.0. (Closes: #678464)
  + Upstream highlights:
    - add "outlined" keyword to marker files
    - update JPL ephemeris code for 64 bit machines
    - add bump_shade config file parameter
    - add opacity keyword for markers
    - implement Rayleigh scattering
  + Remove old Debian patches now obsoleted:
    - possible-buffer-overflow (included upstream)
    - orbit-step-bug (included upstream)
    - gcc44-includes (included upstream)
    - hyphen-used-as-minus-sign (included upstream)
    - tsc-aspect-ratio (included upstream)
  + Update (and include!) FreeMonoBold
* Update libtiff build-deps (Closes: #736052)
* Update libpng build-deps (Closes: #662568)
* Update Standards-Version to 3.9.5
* Update to debhelper 9 to enable hardening build.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libmultiple/RayleighScattering.cpp
 
1
#include <algorithm>
 
2
#include <cmath>
 
3
#include <cstring>
 
4
#include <fstream>
 
5
#include <sstream>
 
6
using namespace std;
 
7
 
 
8
#include "findFile.h"
 
9
#include "parse.h"
 
10
#include "xpUtil.h"
 
11
 
 
12
#include "libimage/Image.h"
 
13
#include "libmultiple/RayleighScattering.h"
 
14
 
 
15
// Calculate the airmass using the expressions in Smith and Smith, JGR
 
16
// 77(19), 3592--3597, 1972
 
17
// x = distance from planet center, units of atmospheric scale height
 
18
// chi = zenith angle
 
19
static double
 
20
chapman(double x, double chi)
 
21
{
 
22
    double chapman;
 
23
    double y = sqrt(x/2) * fabs(cos(chi));
 
24
 
 
25
    if (chi < M_PI_2)
 
26
        chapman = sqrt(x * M_PI_2) * exp(y*y) * erfc(y);
 
27
    else
 
28
        chapman = sqrt(2 * M_PI * x) 
 
29
            * (sqrt(sin(chi))*exp(x*(1-sin(chi)))
 
30
               -0.5*exp(y*y)*erfc(y));
 
31
 
 
32
    return chapman;
 
33
}
 
34
 
 
35
// Find the altitude (units of planetary radius) where the sun sets
 
36
// for a given zenith angle
 
37
static double
 
38
shadowHeight(double sza)
 
39
{
 
40
    double shadowHeight;
 
41
    if (cos(sza) > 0)
 
42
        shadowHeight = 0;
 
43
    else
 
44
        shadowHeight = 1/sin(sza) - 1;
 
45
 
 
46
    return shadowHeight;
 
47
}
 
48
 
 
49
// find the height above the ground given the zenith angle at the ground
 
50
// and the slant path length to the ground
 
51
static double 
 
52
slantHeight(double chi, double slant)
 
53
{
 
54
    double a = 1;
 
55
    double b = 2;
 
56
    double c = -slant * (slant + 2 * cos(chi));
 
57
 
 
58
    return (-b + sqrt(b*b-4*a*c))/2*a; 
 
59
}
 
60
 
 
61
// find the slant path length to the ground given the zenith angle at the ground
 
62
// and the height above the ground
 
63
static double
 
64
slantLength(double chi, double height)
 
65
{
 
66
    double a = 1;
 
67
    double b = 2 * cos(chi);
 
68
    double c = -height * (height + 2);
 
69
 
 
70
    return (-b + sqrt(b*b-4*a*c))/2*a; 
 
71
}
 
72
 
 
73
// find the zenith angle along the slant path given the
 
74
// zenith angle at the ground and the height above the ground
 
75
static double
 
76
zenithAlongSlant(double chi, double height)
 
77
{
 
78
    return (asin(sin(chi) / (1+height)));
 
79
}
 
80
 
 
81
// find the zenith angle at a point along the tangent given the zenith
 
82
// angle at the tangent point.  Set isNear to true if the point is
 
83
// closer than the tangent, false if it is on the far side.  Tangent
 
84
// point is assumed to be at the ground.
 
85
static double
 
86
zenithAlongTangent(double chi, double height, bool isNear)
 
87
{
 
88
    double angle1 = asin(1/(1+height));
 
89
    double angle2 = M_PI_2 - chi;
 
90
 
 
91
    if (isNear)
 
92
    {
 
93
        return angle1 - angle2;
 
94
    }
 
95
    else
 
96
    {
 
97
        return M_PI - angle1 - angle2;
 
98
    }
 
99
}
 
100
 
 
101
void 
 
102
RayleighScattering::createTables()
 
103
{
 
104
    double deg_to_rad = M_PI/180;
 
105
 
 
106
    double planck1 = 1.19104e-6;
 
107
    double planck2 = 0.0143883;
 
108
    double temp = 5700;
 
109
 
 
110
    double peakLambda = 510e-9;
 
111
    double peakRadiance = planck1/(pow(peakLambda,5) 
 
112
                                  * (exp(planck2/(peakLambda*temp))-1));
 
113
 
 
114
    double htop = tanHeight_[tanHeight_.size()-1];
 
115
    double min_dh = scaleHeight_ / 5;
 
116
    
 
117
    // assume tables for each degree between 0 and 180
 
118
    vector<double> phase_d;
 
119
    vector<double> phase;
 
120
    for (int i = 0; i < 181; i++)
 
121
    {
 
122
        phase_d.push_back(i);
 
123
        phase.push_back(i * deg_to_rad);
 
124
    }
 
125
 
 
126
    int area_th = incidence_.size() * tanHeight_.size();
 
127
    int area_em = incidence_.size() * emission_.size();
 
128
        
 
129
    double n2m1 = indexOfRefraction_*indexOfRefraction_-1;
 
130
 
 
131
    // build the limb tables
 
132
    for (unsigned int ip = 0; ip < phase.size(); ip++)
 
133
    {
 
134
        double th_array[lambda_.size() * area_th];
 
135
        for (unsigned int i = 0; i < lambda_.size() * area_th; i++) 
 
136
            th_array[i] = 0;
 
137
 
 
138
        for (unsigned int il = 0; il < lambda_.size(); il++)
 
139
        {
 
140
            double lambda4 = pow(lambda_[il], 4);
 
141
            double rayleighScale = 1/peakRadiance;
 
142
            rayleighScale *= planck1/(pow(lambda_[il],5) 
 
143
                                      * (exp(planck2/(lambda_[il]*temp))-1));
 
144
            rayleighScale *= (2 * pow(M_PI*n2m1, 2)) / (3*numberDensity_);
 
145
            double crossRayleigh = 4*M_PI*rayleighScale / (numberDensity_*lambda4);
 
146
 
 
147
            // fudge factor - scale the red down at low phase angles
 
148
            if (il == 0) rayleighScale *= 0.5 * (1 + (1-cos(phase[ip]/2)));
 
149
 
 
150
            for (unsigned int ii = 0; ii < incidence_.size(); ii++)
 
151
            {
 
152
                if (phase[ip] > (M_PI_2+incidence_[ii])+2*deg_to_rad) continue;
 
153
                if (phase[ip] < fabs(M_PI_2-incidence_[ii])-2*deg_to_rad) continue;
 
154
 
 
155
                double hmin = shadowHeight(incidence_[ii]) * radius_;
 
156
                if (htop < hmin) continue;
 
157
 
 
158
                double cos2phase = cos(phase[ip]) * cos(phase[ip]);
 
159
 
 
160
                for (unsigned int it = 0; it < tanHeight_.size(); it++)
 
161
                {
 
162
                    unsigned int ipos = il*area_th + it*incidence_.size() + ii;
 
163
 
 
164
                    double rayleigh = 0;
 
165
                    // slant distance from the tangent point to htop
 
166
                    double xfar = slantLength(M_PI_2,
 
167
                                              (htop-tanHeight_[it])
 
168
                                              /(radius_+tanHeight_[it]));
 
169
                    xfar *= (radius_+tanHeight_[it]);
 
170
                    int n_h = (int) ceil(xfar/min_dh);
 
171
                    double dx = xfar / (n_h-1);
 
172
 
 
173
                    // the far side of the tangent point
 
174
                    for (int ix = n_h-2; ix >= 0; ix--)
 
175
                    {
 
176
                        double x = (ix + 0.5)* dx;
 
177
                        double h = slantHeight(M_PI_2, 
 
178
                                               x/(radius_+tanHeight_[it]));
 
179
                        double thisIncidence = zenithAlongTangent(incidence_[ii], h, 
 
180
                                                                  (phase[ip]>M_PI_2));
 
181
                        double thisEmission = zenithAlongTangent(M_PI_2, h, false);
 
182
                        h *= (radius_ + tanHeight_[it]);
 
183
                        h += tanHeight_[it];
 
184
 
 
185
                        hmin = shadowHeight(thisIncidence) * radius_;
 
186
                        if (h < hmin) continue;
 
187
 
 
188
                        double atten = exp(-h/scaleHeight_);
 
189
                        double tauScale = numberDensity_*atten*scaleHeight_*crossRayleigh;
 
190
                        double tauSun = chapman((h+radius_)/scaleHeight_, thisIncidence);
 
191
                        double tauView = chapman((h+radius_)/scaleHeight_, thisEmission);
 
192
 
 
193
                        rayleigh += dx * atten * exp(-(tauSun+tauView)*tauScale);
 
194
                    }
 
195
                    // the near side of the tangent point
 
196
                    for (int ix = 0; ix < n_h-1; ix++)
 
197
                    {
 
198
                        double x = (ix+0.5) * dx;
 
199
                        double h = slantHeight(M_PI_2, x/(radius_+tanHeight_[it]));
 
200
                        double thisIncidence = zenithAlongTangent(incidence_[ii], h, 
 
201
                                                                  (phase[ip]<M_PI_2));
 
202
                        double thisEmission = zenithAlongTangent(M_PI_2, h, true);
 
203
                        h *= (radius_ + tanHeight_[it]);
 
204
                        h += tanHeight_[it];
 
205
                        
 
206
                        hmin = shadowHeight(thisIncidence) * radius_;
 
207
                        if (h < hmin) continue;
 
208
 
 
209
                        double atten = exp(-h/scaleHeight_);
 
210
                        double tauScale = numberDensity_*atten*scaleHeight_*crossRayleigh;
 
211
                        double tauSun = chapman((h+radius_)/scaleHeight_, thisIncidence);
 
212
                        double tauView = chapman((h+radius_)/scaleHeight_, thisEmission);
 
213
 
 
214
                        rayleigh += dx * atten * exp(-(tauSun+tauView)*tauScale);
 
215
                    }
 
216
                    rayleigh *= (0.75*(1+cos2phase)*rayleighScale/lambda4);
 
217
 
 
218
                    th_array[ipos] = rayleigh;
 
219
                } // tangent height
 
220
            } // incidence
 
221
        } // lambda
 
222
        char buffer[64];
 
223
        snprintf(buffer, 64, limbTemplate_.c_str(), (int) (phase_d[ip]));
 
224
 
 
225
        writeTable(buffer, th_array, lambda_.size(), incidence_.size(), 
 
226
                   tanHeight_.size(), limbPNG_);
 
227
    }
 
228
 
 
229
    // Now build the tables for the disk
 
230
    for (unsigned int ip = 0; ip < phase.size(); ip++)
 
231
    {
 
232
        double em_array[lambda_.size() * area_em];
 
233
        for (unsigned int i = 0; i < lambda_.size() * area_em; i++) 
 
234
            em_array[i] = 0;
 
235
 
 
236
        double cos2phase = cos(phase[ip]) * cos(phase[ip]);
 
237
        for (unsigned int il = 0; il < lambda_.size(); il++)
 
238
        {
 
239
            double lambda4 = pow(lambda_[il], 4);
 
240
            double rayleighScale = 1/peakRadiance;
 
241
            rayleighScale *= planck1/(pow(lambda_[il],5) 
 
242
                                      * (exp(planck2/(lambda_[il]*temp))-1));
 
243
            rayleighScale *= (2 * pow(M_PI*n2m1, 2)) / (3*numberDensity_);
 
244
            double crossRayleigh = 4*M_PI*rayleighScale 
 
245
                / (numberDensity_*lambda4);
 
246
 
 
247
            for (unsigned int ii = 0; ii < incidence_.size(); ii++)
 
248
            {
 
249
                double hmin = shadowHeight(incidence_[ii]) * radius_;
 
250
                if (hmin > htop) continue;
 
251
 
 
252
                for (unsigned int ie = 0; ie < emission_.size(); ie++)
 
253
                {
 
254
                    if (phase[ip] > (emission_[ie]+incidence_[ii])+2*deg_to_rad) continue;
 
255
                    if (phase[ip] < fabs(emission_[ie]-incidence_[ii])-2*deg_to_rad) continue;
 
256
 
 
257
 
 
258
                    unsigned int ipos = il * area_em + ie * incidence_.size() + ii;
 
259
 
 
260
                    double rayleigh = 0;
 
261
 
 
262
                    // slant distance from the tangent point to htop
 
263
                    double xfar = slantLength(emission_[ie], htop/radius_);
 
264
                    xfar *= radius_;
 
265
                    int n_h = (int) ceil(xfar/min_dh);
 
266
                    double dx = xfar / (n_h-1);
 
267
 
 
268
                    // integrate along the slant path
 
269
                    for (int ix = 0; ix < n_h-1; ix++)
 
270
                    {
 
271
                        double x = (ix+0.5) * dx;
 
272
                        double h = slantHeight(emission_[ie], x/radius_);
 
273
                        double thisIncidence = zenithAlongSlant(incidence_[ii], h);
 
274
                        double thisEmission = zenithAlongSlant(emission_[ie], h);
 
275
                        h *= radius_;
 
276
 
 
277
                        if (h < hmin) continue;
 
278
 
 
279
                        double atten = exp(-h/scaleHeight_);
 
280
                        double tauSun = chapman((h+radius_)/scaleHeight_, thisIncidence);
 
281
                        double tauView = chapman((h+radius_)/scaleHeight_,thisEmission);
 
282
                        double tauScale = numberDensity_ * atten * 
 
283
                            scaleHeight_ * crossRayleigh;
 
284
 
 
285
                        rayleigh += dx * atten * exp(-(tauSun+tauView)*tauScale);
 
286
                    }
 
287
                    rayleigh *= (0.75*(1+cos2phase)*rayleighScale/lambda4);
 
288
 
 
289
                    em_array[ipos] = rayleigh;
 
290
                } // emission
 
291
            } // incidence
 
292
        } // lambda
 
293
 
 
294
        char buffer[64];
 
295
        snprintf(buffer, 64, diskTemplate_.c_str(), (int) (phase_d[ip]));
 
296
 
 
297
        writeTable(buffer, em_array, lambda_.size(), incidence_.size(), 
 
298
                   emission_.size(), diskPNG_);
 
299
    }
 
300
}
 
301
 
 
302
RayleighScattering::RayleighScattering(string configFile) :
 
303
    diskPNG_(false),
 
304
    diskTemplate_(""),
 
305
    indexOfRefraction_(0.),
 
306
    limbPNG_(false),
 
307
    limbTemplate_(""),
 
308
    numberDensity_(0.),
 
309
    radius_(0.),
 
310
    scaleHeight_(0.)
 
311
{
 
312
    bool foundFile = findFile(configFile, "scattering");
 
313
    if (!foundFile)
 
314
    {
 
315
        ostringstream errStr;
 
316
        errStr << "Can't load scattering file " << configFile << "\n";
 
317
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
318
    }
 
319
 
 
320
    phaseDeg_.clear();
 
321
    for (int i = 0; i < 181; i++)
 
322
    {
 
323
        phaseDeg_.push_back((double) i);
 
324
        PNGTable_.insert(make_pair(i, (Image *) NULL));
 
325
        BINTable_.insert(make_pair(i, (double *) NULL));
 
326
    }
 
327
 
 
328
    for (int i = 0; i < 3; i++) scattering_[i] = 0;
 
329
 
 
330
    readConfigFile(configFile);
 
331
}
 
332
 
 
333
void 
 
334
RayleighScattering::readConfigFile(string configFile)
 
335
{
 
336
    ifstream inFile(configFile.c_str());
 
337
    char line[MAX_LINE_LENGTH];
 
338
 
 
339
    vector<double> values;
 
340
    if (!readBlock(inFile, "INCIDENCE %d", values))
 
341
    {
 
342
        ostringstream errStr;
 
343
        errStr << "INCIDENCE block not found in " << configFile << "\n";
 
344
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
345
    }
 
346
    incidence_.clear();
 
347
    for (unsigned int ii = 0; ii < values.size(); ii++)
 
348
        incidence_.push_back(values[ii] * deg_to_rad);
 
349
 
 
350
    if (!readBlock(inFile, "EMISSION %d", values))
 
351
    {
 
352
        ostringstream errStr;
 
353
        errStr << "EMISSION block not found in " << configFile << "\n";
 
354
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
355
    }
 
356
    emission_.clear();
 
357
    for (unsigned int ii = 0; ii < values.size(); ii++)
 
358
        emission_.push_back(values[ii] * deg_to_rad);
 
359
 
 
360
    if (!readBlock(inFile, "TANGENT_HEIGHT %d", values))
 
361
    {
 
362
        ostringstream errStr;
 
363
        errStr << "TANGENT_HEIGHT block not found in " << configFile << "\n";
 
364
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
365
    }
 
366
    tanHeight_.clear();
 
367
    for (unsigned int ii = 0; ii < values.size(); ii++)
 
368
        tanHeight_.push_back(values[ii] * 1e3);
 
369
 
 
370
    diskTemplate_.clear();
 
371
    limbTemplate_.clear();
 
372
    while (inFile.getline(line, MAX_LINE_LENGTH, '\n') != NULL)
 
373
    {
 
374
        int i = 0;
 
375
        while (isDelimiter(line[i]))
 
376
        {
 
377
            i++;
 
378
            if (static_cast<unsigned int> (i) > strlen(line)) break;
 
379
        }
 
380
        if (static_cast<unsigned int> (i) > strlen(line)) continue;
 
381
        if (isEndOfLine(line[i])) continue;
 
382
 
 
383
        char templ1[MAX_LINE_LENGTH];
 
384
        char templ2[MAX_LINE_LENGTH];
 
385
        sscanf(line, "TEMPLATES %s %s", templ1, templ2);
 
386
 
 
387
        diskTemplate_.assign(templ1);
 
388
        limbTemplate_.assign(templ2);
 
389
 
 
390
        size_t found = diskTemplate_.find(".png");
 
391
        if (found == string::npos) found = diskTemplate_.find(".PNG");
 
392
        diskPNG_ = (found != string::npos);
 
393
            
 
394
        found = limbTemplate_.find(".png");
 
395
        if (found == string::npos) found = limbTemplate_.find(".PNG");
 
396
        limbPNG_ = (found != string::npos);
 
397
           
 
398
        break;
 
399
    }
 
400
    if (diskTemplate_.length() == 0) 
 
401
    {
 
402
        ostringstream errStr;
 
403
        errStr << "TEMPLATE disk file not found in " << configFile << "\n";
 
404
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
405
    }
 
406
    if (limbTemplate_.length() == 0) 
 
407
    {
 
408
        ostringstream errStr;
 
409
        errStr << "TEMPLATE limb file not found in " << configFile << "\n";
 
410
        xpExit(errStr.str(), __FILE__, __LINE__);       
 
411
    }
 
412
 
 
413
    // remaining parameters are only required for table generation
 
414
    readBlock(inFile, "WAVELENGTHS %d", values);
 
415
    lambda_.clear();
 
416
    for (unsigned int ii = 0; ii < values.size(); ii++)
 
417
        lambda_.push_back(values[ii] * 1e-9);
 
418
 
 
419
    readValue(inFile, "RADIUS %lf", radius_);
 
420
    radius_ *= 1e3;
 
421
 
 
422
    readValue(inFile, "SCALE_HEIGHT %lf", scaleHeight_);
 
423
 
 
424
    readValue(inFile, "INDEX_OF_REFRACTION %lf", indexOfRefraction_);
 
425
 
 
426
    readValue(inFile, "DENSITY %lf", numberDensity_);
 
427
}
 
428
 
 
429
RayleighScattering::~RayleighScattering()
 
430
{
 
431
    clearTables();
 
432
}
 
433
 
 
434
bool
 
435
RayleighScattering::readBlock(ifstream &inFile, 
 
436
                              const char *format, 
 
437
                              vector<double> &values)
 
438
{
 
439
    values.clear();
 
440
 
 
441
    char line[MAX_LINE_LENGTH];
 
442
    while (inFile.getline(line, MAX_LINE_LENGTH, '\n') != NULL)
 
443
    {
 
444
        int i = 0;
 
445
        while (isDelimiter(line[i]))
 
446
        {
 
447
            i++;
 
448
            if (static_cast<unsigned int> (i) > strlen(line)) break;
 
449
        }
 
450
        if (static_cast<unsigned int> (i) > strlen(line)) continue;
 
451
        if (isEndOfLine(line[i])) continue;
 
452
 
 
453
        int size;
 
454
        if (sscanf(line, format, &size) == 0) break;
 
455
 
 
456
        for (int ii = 0; ii < size; ii++) 
 
457
        {
 
458
            double thisValue;
 
459
            inFile >> thisValue;
 
460
            values.push_back(thisValue);
 
461
        }
 
462
        return true;
 
463
    }
 
464
    return false;
 
465
}
 
466
 
 
467
bool
 
468
RayleighScattering::readValue(ifstream &inFile, 
 
469
                              const char *format, 
 
470
                              double &value)
 
471
{
 
472
    char line[MAX_LINE_LENGTH];
 
473
    while (inFile.getline(line, MAX_LINE_LENGTH, '\n') != NULL)
 
474
    {
 
475
        int i = 0;
 
476
        while (isDelimiter(line[i]))
 
477
        {
 
478
            i++;
 
479
            if (static_cast<unsigned int> (i) > strlen(line)) break;
 
480
        }
 
481
        if (static_cast<unsigned int> (i) > strlen(line)) continue;
 
482
        if (isEndOfLine(line[i])) continue;
 
483
 
 
484
        if (sscanf(line, format, &value) == 0) break;
 
485
 
 
486
        return true;
 
487
    }
 
488
 
 
489
    return false;
 
490
}
 
491
 
 
492
void
 
493
RayleighScattering::clear()
 
494
{ 
 
495
    for (int i = 0; i < 3; i ++) scattering_[i] = 0; 
 
496
}
 
497
 
 
498
void
 
499
RayleighScattering::clearTables()
 
500
{ 
 
501
    for (map<int, Image *>::iterator it = PNGTable_.begin(); 
 
502
         it != PNGTable_.end(); it++)
 
503
    {
 
504
        delete it->second;
 
505
        it->second = NULL;
 
506
    }
 
507
    for (map<int, double *>::iterator it = BINTable_.begin(); 
 
508
         it != BINTable_.end(); it++)
 
509
    {
 
510
        delete [] it->second;
 
511
        it->second = NULL;
 
512
    }
 
513
}
 
514
 
 
515
// store a double (0 < x < 1) as a four byte array of unsigned char
 
516
// big endian order: 
 
517
// four bytes ABCD -> alpha = A, red = B, green = C, blue = D
 
518
void 
 
519
RayleighScattering::doubleToARGB(double x, unsigned char argb[4])
 
520
{
 
521
    static const unsigned long maxValue = 0xffffffff;
 
522
    unsigned long ul = (unsigned long) (x * maxValue);
 
523
    
 
524
    argb[0] = (ul >> 24) & 0xff;
 
525
    argb[1] = (ul >> 16) & 0xff;
 
526
    argb[2] = (ul >> 8) & 0xff;
 
527
    argb[3] =  ul & 0xff;
 
528
}
 
529
 
 
530
// turn a four byte unsigned char array to a double (0 < x < 1)
 
531
// Note: four bytes is only good enough for a float
 
532
double
 
533
RayleighScattering::ARGBToDouble(unsigned char argb[4])
 
534
{
 
535
    static const unsigned long maxValue = 0xffffffff;
 
536
    unsigned long ul = argb[0] << 24 | argb[1] << 16 | argb[2] << 8 | argb[3];
 
537
 
 
538
    double x = ul;
 
539
    x /= maxValue;
 
540
 
 
541
    return x;
 
542
}
 
543
void
 
544
RayleighScattering::writeTable(const char *buffer, double *array, 
 
545
                               size_t dim0, size_t dim1, size_t dim2, 
 
546
                               bool usePNG)
 
547
{
 
548
    if (usePNG)
 
549
    {
 
550
        unsigned int area = dim1 * dim2;
 
551
            
 
552
        unsigned char rgb[dim0 * 3 * area];
 
553
        unsigned char alpha[dim0 * area];
 
554
        memset(rgb, 0, dim0 * 3 * area);
 
555
        memset(alpha, 0, dim0 * area);
 
556
            
 
557
        unsigned char argb[4];
 
558
 
 
559
        for (unsigned int il = 0; il < dim0; il++)
 
560
        {
 
561
            for (unsigned int it = 0; it < dim2; it++)
 
562
            {
 
563
                for (unsigned int ii = 0; ii < dim1; ii++)
 
564
                {
 
565
                    unsigned int ipos = il * area + it * dim1 + ii;
 
566
                    doubleToARGB(array[ipos], argb);
 
567
                    alpha[il * area + it * dim1 + ii] = argb[0];
 
568
                    memcpy(rgb+3*(il * area + it*dim1 + ii), 
 
569
                           argb + 1, 3);
 
570
                }
 
571
            }
 
572
        }
 
573
        Image pngImage(dim1, dim0*dim2, rgb, alpha);
 
574
        if (!pngImage.Write(buffer))
 
575
        {
 
576
            ostringstream errStr;
 
577
            errStr << "Can't write scattering file " << buffer << "\n";
 
578
            xpWarn(errStr.str(), __FILE__, __LINE__);       
 
579
        }
 
580
        else
 
581
        {
 
582
            string message("Wrote ");
 
583
            message.append(buffer);
 
584
            message.append("\n");
 
585
            xpMsg(message, __FILE__, __LINE__);
 
586
        }  
 
587
    }
 
588
    else
 
589
    {
 
590
        FILE *outfile = fopen(buffer, "wb");
 
591
        if (outfile == NULL)
 
592
        {
 
593
            ostringstream errStr;
 
594
            errStr << "Can't write scattering file " << buffer << "\n";
 
595
            xpWarn(errStr.str(), __FILE__, __LINE__);       
 
596
        }
 
597
        else
 
598
        {
 
599
            fwrite(&dim0, sizeof(size_t), 1, outfile);
 
600
            fwrite(&dim1, sizeof(size_t), 1, outfile);
 
601
            fwrite(&dim2, sizeof(size_t), 1, outfile);
 
602
            fwrite(array, sizeof(double), dim0*dim1*dim2, outfile);
 
603
            fclose(outfile);
 
604
            string message("Wrote ");
 
605
            message.append(buffer);
 
606
            message.append("\n");
 
607
            xpMsg(message, __FILE__, __LINE__);
 
608
        }
 
609
    }
 
610
}
 
611
 
 
612
double *
 
613
RayleighScattering::readBinaryTable(const char *filename)
 
614
{
 
615
    double *dblArray = NULL;
 
616
 
 
617
    FILE *inFile = fopen(filename, "rb");
 
618
 
 
619
    if (inFile != NULL)
 
620
    {
 
621
        size_t dim0, dim1, dim2;
 
622
        fread(&dim0, sizeof(size_t), 1, inFile);
 
623
        fread(&dim1, sizeof(size_t), 1, inFile);
 
624
        fread(&dim2, sizeof(size_t), 1, inFile);
 
625
 
 
626
        size_t size = dim0*dim1*dim2;
 
627
 
 
628
        dblArray = new double[size];
 
629
        fread(dblArray, sizeof(double), size, inFile);
 
630
    }
 
631
 
 
632
    fclose(inFile);
 
633
 
 
634
    return dblArray;
 
635
}
 
636
 
 
637
double
 
638
RayleighScattering::getColor(int index)
 
639
{
 
640
    return scattering_[index];
 
641
}
 
642
 
 
643
void
 
644
RayleighScattering::calcScatteringLimb(double inc, double tanht, double phase)
 
645
{
 
646
    for (int ic = 0; ic < 3; ic++)
 
647
        calcScattering(inc, tanht, phase, ic, tanHeight_, limbTemplate_, 
 
648
                       limbPNG_);
 
649
}
 
650
 
 
651
void
 
652
RayleighScattering::calcScatteringDisk(double inc, double ems, double phase)
 
653
{
 
654
    if (phase > M_PI_2) phase = M_PI - phase;
 
655
    for (int ic = 0; ic < 3; ic++)
 
656
        calcScattering(inc, ems, phase, ic, emission_, diskTemplate_, diskPNG_);
 
657
}
 
658
 
 
659
void
 
660
RayleighScattering::readTableValue(string thisTemplate, bool usePNG, 
 
661
                                   int color, int x, int y, 
 
662
                                   vector<double> yaxis, 
 
663
                                   int phase, double values[4])
 
664
{
 
665
    if (usePNG)
 
666
    {
 
667
        map<int, Image *>::iterator it = PNGTable_.find(phase);
 
668
        Image *image = NULL;
 
669
        if (it->second == NULL)
 
670
        {
 
671
            char buffer[64];
 
672
            snprintf(buffer, 64, thisTemplate.c_str(), phase);
 
673
            string thisFile(buffer);
 
674
            if (!findFile(thisFile, "scattering")) 
 
675
            {
 
676
                ostringstream errStr;
 
677
                errStr << "Can't load scattering file " << thisFile << "\n";
 
678
                xpExit(errStr.str(), __FILE__, __LINE__);
 
679
            }
 
680
            image = new Image();
 
681
            if (!image->Read(thisFile.c_str()))
 
682
            {
 
683
                scattering_[color] = 0;
 
684
                return;
 
685
            }
 
686
            it->second = image;
 
687
        }
 
688
 
 
689
        unsigned char argb[4];
 
690
        image = it->second;
 
691
        image->getPixel(x, y, argb+1, argb);
 
692
        values[0] = ARGBToDouble(argb);
 
693
        image->getPixel(x+1, y, argb+1, argb);
 
694
        values[1] = ARGBToDouble(argb);
 
695
        image->getPixel(x, y+1, argb+1, argb);
 
696
        values[2] = ARGBToDouble(argb);
 
697
        image->getPixel(x+1, y+1, argb+1, argb);
 
698
        values[3] = ARGBToDouble(argb);
 
699
    }
 
700
    else
 
701
    {
 
702
        map<int, double *>::iterator it = BINTable_.find(phase);
 
703
        double *array = NULL;
 
704
        if (it->second == NULL)
 
705
        {
 
706
            char buffer[64];
 
707
            snprintf(buffer, 64, thisTemplate.c_str(), phase);
 
708
            string thisFile(buffer);
 
709
            if (!findFile(thisFile, "scattering"))
 
710
            {
 
711
                ostringstream errStr;
 
712
                errStr << "Can't load scattering file " << thisFile << "\n";
 
713
                xpExit(errStr.str(), __FILE__, __LINE__);
 
714
            }
 
715
            array = readBinaryTable(thisFile.c_str());
 
716
            if (array == NULL)
 
717
            {
 
718
                scattering_[color] = 0;
 
719
                return;
 
720
            }
 
721
            it->second = array;
 
722
        }
 
723
 
 
724
        array = it->second;
 
725
        unsigned int ipos[4];
 
726
        ipos[0] = y * incidence_.size() + x;
 
727
        ipos[1] = y * incidence_.size() + x+1;
 
728
        ipos[2] = (y+1) * incidence_.size() + x;
 
729
        ipos[3] = (y+1) * incidence_.size() + x+1;
 
730
        unsigned int isize = 3 * yaxis.size() * incidence_.size();
 
731
        for (int i = 0; i < 4; i++) 
 
732
        {
 
733
            if (isize > ipos[i]) values[i] = array[ipos[i]];
 
734
            else values[i] = 0;
 
735
        }
 
736
    }
 
737
}
 
738
 
 
739
void
 
740
RayleighScattering::calcScattering(double inc, double yValue, double phase, 
 
741
                                   int color, vector<double> yaxis, 
 
742
                                   string thisTemplate, bool usePNG)
 
743
{
 
744
    double deg_to_rad = M_PI / 180;
 
745
 
 
746
    if (inc < incidence_.front() || inc > incidence_.back() \
 
747
        || yValue < yaxis.front() || yValue > yaxis.back())
 
748
    {
 
749
        scattering_[color] = 0;
 
750
        return;
 
751
    }
 
752
 
 
753
    double phase_deg = phase / deg_to_rad;
 
754
 
 
755
    int phase_lo = floor(phase_deg);
 
756
    if (phase_lo < 0) phase_lo = 0;
 
757
 
 
758
    int phase_hi = ceil(phase_deg);
 
759
    if (phase_hi - phase_deg < 1e-3) phase_hi = phase_lo+1;
 
760
 
 
761
    vector<double>::iterator inc_hi, inc_lo;
 
762
    // lower_bound returns first element greater than or equal to value
 
763
    inc_lo = lower_bound(incidence_.begin(), incidence_.end(), inc);
 
764
    if (inc_lo > incidence_.begin()) inc_lo--;
 
765
    // upper_bound returns first element greater than value
 
766
    inc_hi = upper_bound(incidence_.begin(), incidence_.end(), inc);
 
767
 
 
768
    vector<double>::iterator y_hi, y_lo;
 
769
    y_lo = lower_bound(yaxis.begin(), yaxis.end(), yValue);
 
770
    if (y_lo > yaxis.begin()) y_lo--;
 
771
    y_hi = upper_bound(yaxis.begin(), yaxis.end(), yValue);
 
772
 
 
773
    int x = inc_lo - incidence_.begin();
 
774
    double xFrac = (*inc_hi == *inc_lo ? 0 : (inc - *inc_lo) / (*inc_hi - *inc_lo));
 
775
    int y = y_lo - yaxis.begin();
 
776
    double yFrac = (*y_hi == *y_lo ? 0 : (yValue - *y_lo) / (*y_hi - *y_lo));
 
777
 
 
778
    double weight[4];
 
779
    getWeights(xFrac, 1-yFrac, weight);
 
780
 
 
781
    y += color * yaxis.size();
 
782
 
 
783
    double values[4];
 
784
 
 
785
    readTableValue(thisTemplate, usePNG, color, x, y, yaxis, phase_lo, values);
 
786
    double loPhase = 0;
 
787
    for (int i = 0; i < 4; i++) 
 
788
        loPhase += weight[i] * values[i];
 
789
 
 
790
    readTableValue(thisTemplate, usePNG, color, x, y, yaxis, phase_hi, values);
 
791
    double hiPhase = 0;
 
792
    for (int i = 0; i < 4; i++) 
 
793
        hiPhase += weight[i] * values[i];
 
794
 
 
795
    scattering_[color] = loPhase;
 
796
    if (phase_hi - phase_lo != 0)
 
797
        scattering_[color] += (phase_deg-phase_lo)/(phase_hi-phase_lo) 
 
798
            * (hiPhase - loPhase);
 
799
 
 
800
    if (scattering_[color] < 0) scattering_[color] = 0;
 
801
}