← Back to branch summary

~ubuntu-branches/ubuntu/warty/xplanet/warty

« back to all changes in this revision

Viewing changes to src/libmultiple/drawEllipsoid.cpp

  • Committer: Bazaar Package Importer
  • Author(s): LaMont Jones
  • Date: 2004-08-24 07:14:00 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040824071400-2dr4qnjbjmm8z3ia
Tags: 1.0.6-1ubuntu1
Build-depend: libtiff4-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libmultiple/drawEllipsoid.cpp
 
1
#include "Map.h"
 
2
#include "Options.h"
 
3
#include "View.h"
 
4
#include "xpUtil.h"
 
5
 
 
6
#include "libdisplay/libdisplay.h"
 
7
#include "libplanet/Planet.h"
 
8
 
 
9
void
 
10
drawEllipsoid(const double pX, const double pY, const double pR, 
 
11
              const double oX, const double oY, const double oZ, 
 
12
              const double X, const double Y, const double Z,
 
13
              DisplayBase *display, 
 
14
              const View *view, const Map *map, Planet *planet,
 
15
              const double planetRadius)
 
16
{
 
17
    double lat, lon;
 
18
    unsigned char color[3];
 
19
 
 
20
    const int j0 = 0;
 
21
    const int j1 = display->Height();
 
22
    const int i0 = 0;
 
23
    const int i1 = display->Width();
 
24
 
 
25
    // P1 (Observer) is at (oX, oY, oZ), or (0, 0, 0) in view
 
26
    // coordinates
 
27
    // P2 (current pixel) is on the line from P1 to (vX, vY, vZ)
 
28
    // P3 (Planet center) is at (X, Y, Z) in heliocentric rectangular
 
29
    // Now find the intersection of the line with the planet's ellipsoid. 
 
30
 
 
31
    // Define a new coordinate system: planetary XYZ
 
32
    // X = 0 lon, Z = north pole, Y = Z x X
 
33
    // Convert P1, P2, P3 to planetary XYZ:
 
34
    // 1) Convert to planetocentric
 
35
    // 2) Convert to planetary XYZ (units of planetary radius)
 
36
 
 
37
    double p1X = 0, p1Y = 0, p1Z = 0;
 
38
    double p2X = 0, p2Y = 0, p2Z = 0;
 
39
    double p3X = 0, p3Y = 0, p3Z = 0;
 
40
            
 
41
    const double ratio = 1/(1 - planet->Flattening());
 
42
 
 
43
    planet->XYZToPlanetaryXYZ(oX, oY, oZ, p1X, p1Y, p1Z);
 
44
    p1Z *= ratio;
 
45
 
 
46
    const double c = 2 * (dot(p1X - p3X, p1Y - p3Y, p1Z - p3Z, 
 
47
                              p1X - p3X, p1Y - p3Y, p1Z - p3Z) 
 
48
                          - planetRadius * planetRadius);
 
49
 
 
50
    Options *options = Options::getInstance();
 
51
 
 
52
    // compute the value of the determinant at the center of the body
 
53
    view->PixelToViewCoordinates(options->getCenterX() - pX, 
 
54
                                 options->getCenterY() - pY, 
 
55
                                 p2X, p2Y, p2Z);
 
56
    
 
57
    view->RotateToXYZ(p2X, p2Y, p2Z, p2X, p2Y, p2Z);
 
58
    planet->XYZToPlanetaryXYZ(p2X, p2Y, p2Z, p2X, p2Y, p2Z);
 
59
    p2Z *= ratio;
 
60
    const double centerA = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z, 
 
61
                                    p2X - p1X, p2Y - p1Y, p2Z - p1Z));
 
62
    const double centerB = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z, 
 
63
                                    p1X - p3X, p1Y - p3Y, p1Z - p3Z));
 
64
    const double centerDet = centerB * centerB - centerA * c;
 
65
 
 
66
    for (int j = j0; j < j1; j++)
 
67
    {
 
68
        for (int i = i0; i < i1; i++)
 
69
        {
 
70
            const double dX = options->getCenterX() - i;
 
71
            const double dY = options->getCenterY() - j;
 
72
 
 
73
            view->PixelToViewCoordinates(dX, dY, p2X, p2Y, p2Z);
 
74
 
 
75
            view->RotateToXYZ(p2X, p2Y, p2Z, p2X, p2Y, p2Z);
 
76
            planet->XYZToPlanetaryXYZ(p2X, p2Y, p2Z, p2X, p2Y, p2Z);
 
77
            p2Z *= ratio;
 
78
 
 
79
            const double a = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z, 
 
80
                                      p2X - p1X, p2Y - p1Y, p2Z - p1Z));
 
81
            const double b = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z, 
 
82
                                      p1X - p3X, p1Y - p3Y, p1Z - p3Z));
 
83
            const double determinant = b*b - a * c;
 
84
 
 
85
            if (determinant < 0) continue;
 
86
 
 
87
            double u = -(b + sqrt(determinant));
 
88
            u /= a;
 
89
                    
 
90
            // coordinates of the intersection point
 
91
            double iX, iY, iZ;
 
92
            iX = p1X + u * (p2X - p1X);
 
93
            iY = p1Y + u * (p2Y - p1Y);
 
94
            iZ = p1Z + u * (p2Z - p1Z);
 
95
 
 
96
            iZ /= ratio;
 
97
                    
 
98
            planet->PlanetaryXYZToXYZ(iX, iY, iZ, iX, iY, iZ);
 
99
            planet->XYZToPlanetographic(iX, iY, iZ, lat, lon);
 
100
 
 
101
            map->GetPixel(lat, lon, color);
 
102
            double darkening = ndot(X - iX, Y - iY, Z - iZ, 
 
103
                                    X - oX, Y - oY, Z - oZ);
 
104
            if (darkening < 0) 
 
105
                darkening = 0;
 
106
            else
 
107
                darkening = sqrt(darkening);
 
108
 
 
109
            for (int k = 0; k < 3; k++) 
 
110
                color[k] = static_cast<unsigned char> (color[k] 
 
111
                                                       * darkening);
 
112
            double opacity = 1;
 
113
            if (pR * determinant/centerDet < 10)
 
114
            {
 
115
                opacity = 1 - pow(1-determinant/centerDet, pR);
 
116
            }
 
117
            display->setPixel(i, j, color, opacity);
 
118
        }
 
119
    }
 
120
}