1 #include "Map.h"
2 #include "Options.h"
3 #include "PlanetProperties.h"
4 #include "View.h"
5 #include "xpUtil.h"
6
7 #include "libdisplay/libdisplay.h"
8 #include "libmultiple/libmultiple.h"
9 #include "libmultiple/RayleighScattering.h"
10 #include "libplanet/Planet.h"
11
12 void
drawSphere(const double pX,const double pY,const double pR,const double oX,const double oY,const double oZ,const double X,const double Y,const double Z,DisplayBase * display,const View * view,const Map * map,Planet * planet,PlanetProperties * planetProperties)13 drawSphere(const double pX, const double pY, const double pR,
14 const double oX, const double oY, const double oZ,
15 const double X, const double Y, const double Z,
16 DisplayBase *display,
17 const View *view, const Map *map, Planet *planet,
18 PlanetProperties *planetProperties)
19 {
20 double lat, lon;
21 unsigned char color[3];
22
23 const int j0 = 0;
24 const int j1 = display->Height();
25 const int i0 = 0;
26 const int i1 = display->Width();
27
28 // P1 (Observer) is at (oX, oY, oZ), or (0, 0, 0) in view
29 // coordinates
30 // P2 (current pixel) is on the line from P1 to (vX, vY, vZ)
31 // P3 (Planet center) is at (X, Y, Z) in heliocentric rectangular
32 // Now find the intersection of the line with the planet's sphere.
33 // This algorithm is from
34 // http://astronomy.swin.edu.au/~pbourke/geometry/sphereline
35
36 const double p1X = 0, p1Y = 0, p1Z = 0;
37 double p2X, p2Y, p2Z;
38 double p3X, p3Y, p3Z;
39
40 view->RotateToViewCoordinates(X, Y, Z, p3X, p3Y, p3Z);
41
42 const double planetRadius = planet->Radius() * planetProperties->Magnify();
43 const double c = 2 * (dot(p1X - p3X, p1Y - p3Y, p1Z - p3Z,
44 p1X - p3X, p1Y - p3Y, p1Z - p3Z)
45 - planetRadius * planetRadius);
46
47 Options *options = Options::getInstance();
48
49 // compute the value of the determinant at the center of the body
50 view->PixelToViewCoordinates(options->CenterX() - pX,
51 options->CenterY() - pY,
52 p2X, p2Y, p2Z);
53
54 const double centerA = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z,
55 p2X - p1X, p2Y - p1Y, p2Z - p1Z));
56 const double centerB = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z,
57 p1X - p3X, p1Y - p3Y, p1Z - p3Z));
58 const double centerDet = centerB * centerB - centerA * c;
59
60 double plX, plY, plZ;
61 planet->getPosition(plX, plY, plZ);
62
63 RayleighScattering *rayleighDisk = NULL;
64 RayleighScattering *rayleighLimb = NULL;
65 double rayleighScale = planetProperties->RayleighScale();
66 if (rayleighScale > 0)
67 {
68 rayleighDisk = new RayleighScattering(planetProperties->RayleighFile());
69 rayleighLimb = new RayleighScattering(planetProperties->RayleighFile());
70
71 double radiansPerPixel = options->FieldOfView() / display->Width();
72 double dX = plX - oX;
73 double dY = plY - oY;
74 double dZ = plZ - oZ;
75 double targetDist = sqrt(dX*dX + dY*dY + dZ*dZ);
76
77 double kmPerPixel = radiansPerPixel * targetDist * AU_to_km;
78 double minRes = 2 * rayleighLimb->getScaleHeightKm()
79 * planetProperties->RayleighLimbScale();
80 if (kmPerPixel > minRes)
81 {
82 delete rayleighLimb;
83 rayleighLimb = NULL;
84 }
85 }
86
87 for (int j = j0; j < j1; j++)
88 {
89 for (int i = i0; i < i1; i++)
90 {
91 const double dX = options->CenterX() - i;
92 const double dY = options->CenterY() - j;
93
94 view->PixelToViewCoordinates(dX, dY, p2X, p2Y, p2Z);
95
96 const double a = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z,
97 p2X - p1X, p2Y - p1Y, p2Z - p1Z));
98 const double b = 2 * (dot(p2X - p1X, p2Y - p1Y, p2Z - p1Z,
99 p1X - p3X, p1Y - p3Y, p1Z - p3Z));
100 const double determinant = b*b - a * c;
101
102 double u;
103 double iX, iY, iZ;
104
105 if (rayleighLimb != NULL)
106 {
107 u = -b/a;
108 iX = p1X + u * (p2X - p1X);
109 iY = p1Y + u * (p2Y - p1Y);
110 iZ = p1Z + u * (p2Z - p1Z);
111 view->RotateToXYZ(iX, iY, iZ, iX, iY, iZ);
112
113 double lon, lat, rad;
114 planet->XYZToPlanetographic(iX, iY, iZ, lat, lon, rad);
115 if (rad >= 1 || pR * determinant/centerDet < 10)
116 {
117 double incidence = acos(ndot(iX-plX, iY-plY, iZ-plZ,
118 -iX, -iY, -iZ));
119 double phase = acos(ndot(oX-iX, oY-iY, oZ-iZ,
120 -iX, -iY, -iZ));
121
122 double tanht = planet->Radius() * AU_to_km * 1e3 * (rad-1);
123 if (tanht < 0) tanht = 0;
124 if (planetProperties->RayleighLimbScale() > 0)
125 tanht /= planetProperties->RayleighLimbScale();
126 rayleighLimb->calcScatteringLimb(incidence, tanht, phase);
127
128 display->getPixel(i, j, color);
129 double opacity[3] = { 0, 0, 0 };
130 for (int ic = 0; ic < 3; ic++)
131 {
132 double thisColor = rayleighLimb->getColor(ic);
133 thisColor = (rayleighScale * 255 * thisColor);
134 if (thisColor > 255) thisColor = 255;
135 color[ic] = thisColor;
136 opacity[ic] = thisColor / 255;
137 }
138 display->setPixel(i, j, color, opacity);
139 }
140 rayleighLimb->clear();
141 }
142
143 if (determinant < 0) continue;
144
145 u = -(b + sqrt(determinant));
146 u /= a;
147
148 // if the intersection point is behind the observer, don't
149 // plot it
150 if (u < 0) continue;
151
152 // coordinates of the intersection point
153 iX = p1X + u * (p2X - p1X);
154 iY = p1Y + u * (p2Y - p1Y);
155 iZ = p1Z + u * (p2Z - p1Z);
156
157 view->RotateToXYZ(iX, iY, iZ, iX, iY, iZ);
158 planet->XYZToPlanetographic(iX, iY, iZ, lat, lon);
159
160 map->GetPixel(lat, lon, color);
161
162 if (rayleighDisk != NULL)
163 {
164 double incidence = acos(ndot(iX-plX, iY-plY, iZ-plZ,
165 -iX, -iY, -iZ));
166 double emission = acos(ndot(iX-plX, iY-plY, iZ-plZ,
167 oX-iX, oY-iY, oZ-iZ));
168 double phase = acos(ndot(oX-iX, oY-iY, oZ-iZ,
169 -iX, -iY, -iZ));
170
171 double emsScale = 1;
172 if (planetProperties->RayleighEmissionWeight() > 0)
173 emsScale = pow(sin(emission),
174 planetProperties->RayleighEmissionWeight());
175
176 rayleighDisk->calcScatteringDisk(incidence, emission, phase);
177 for (int ic = 0; ic < 3; ic++)
178 {
179 double thisColor = rayleighDisk->getColor(ic) * emsScale;
180 thisColor = (rayleighScale * 255 * thisColor + color[ic]);
181 if (thisColor > 255) thisColor = 255;
182 color[ic] = thisColor;
183 }
184 rayleighDisk->clear();
185 }
186
187 double darkening = 1;
188 if (planet->Index() != SUN && rayleighScale <= 0)
189 {
190 darkening = ndot(X - iX, Y - iY, Z - iZ,
191 X - oX, Y - oY, Z - oZ);
192 if (darkening < 0)
193 darkening = 0;
194 else
195 darkening = photoFunction(darkening);
196 }
197
198 for (int k = 0; k < 3; k++)
199 color[k] = static_cast<unsigned char> (color[k]
200 * darkening);
201 double opacity = 1;
202 if (pR * determinant/centerDet < 10)
203 opacity = 1 - pow(1-determinant/centerDet, pR);
204 display->setPixel(i, j, color, opacity);
205 }
206 }
207
208 delete rayleighDisk;
209 delete rayleighLimb;
210 }
211