1 // oursun.hxx -- model earth's sun
2 //
3 // Written by Durk Talsma. Originally started October 1997, for distribution
4 // with the FlightGear project. Version 2 was written in August and
5 // September 1998. This code is based upon algorithms and data kindly
6 // provided by Mr. Paul Schlyter. (pausch@saaf.se).
7 //
8 // Separated out rendering pieces and converted to ssg by Curt Olson,
9 // March 2000
10 // This library is free software; you can redistribute it and/or
11 // modify it under the terms of the GNU Library General Public
12 // License as published by the Free Software Foundation; either
13 // version 2 of the License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
18 // Library General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with this program; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
23 //
24 // $Id$
25 
26 
27 #ifdef HAVE_CONFIG_H
28 #  include <simgear_config.h>
29 #endif
30 
31 #include <simgear/compiler.h>
32 
33 #include <osg/AlphaFunc>
34 #include <osg/BlendFunc>
35 #include <osg/Fog>
36 #include <osg/Geode>
37 #include <osg/Geometry>
38 #include <osg/Material>
39 #include <osg/ShadeModel>
40 #include <osg/TexEnv>
41 #include <osg/Texture2D>
42 #include <osgDB/ReadFile>
43 
44 #include <simgear/math/SGMath.hxx>
45 #include <simgear/scene/util/SGReaderWriterOptions.hxx>
46 #include <simgear/screen/colors.hxx>
47 #include <simgear/scene/model/model.hxx>
48 #include "oursun.hxx"
49 
50 using namespace simgear;
51 
52 // Constructor
SGSun(void)53 SGSun::SGSun( void ) :
54     visibility(-9999.0), prev_sun_angle(-9999.0), path_distance(60000.0),
55     sun_exp2_punch_through(7.0e-06), horizon_angle(0.0)
56 {
57 
58 }
59 
60 
61 // Destructor
~SGSun(void)62 SGSun::~SGSun( void ) {
63 }
64 
65 
66 // initialize the sun object and connect it into our scene graph root
67 osg::Node*
build(SGPath path,double sun_size,SGPropertyNode * property_tree_Node)68 SGSun::build( SGPath path, double sun_size, SGPropertyNode *property_tree_Node ) {
69 
70     env_node = property_tree_Node;
71 
72     osg::ref_ptr<SGReaderWriterOptions> options;
73     options = SGReaderWriterOptions::fromPath(path);
74     // build the ssg scene graph sub tree for the sky and connected
75     // into the provide scene graph branch
76     sun_transform = new osg::MatrixTransform;
77     osg::StateSet* stateSet = sun_transform->getOrCreateStateSet();
78 
79     osg::TexEnv* texEnv = new osg::TexEnv;
80     texEnv->setMode(osg::TexEnv::MODULATE);
81     stateSet->setTextureAttribute(0, texEnv, osg::StateAttribute::ON);
82 
83     osg::Material* material = new osg::Material;
84     material->setColorMode(osg::Material::AMBIENT_AND_DIFFUSE);
85     material->setEmission(osg::Material::FRONT_AND_BACK, osg::Vec4(0,0,0,1));
86     material->setSpecular(osg::Material::FRONT_AND_BACK, osg::Vec4(0,0,0,1));
87     stateSet->setAttribute(material);
88 
89     osg::ShadeModel* shadeModel = new osg::ShadeModel;
90     shadeModel->setMode(osg::ShadeModel::SMOOTH);
91     stateSet->setAttributeAndModes(shadeModel);
92 
93     osg::AlphaFunc* alphaFunc = new osg::AlphaFunc;
94     alphaFunc->setFunction(osg::AlphaFunc::ALWAYS);
95     stateSet->setAttributeAndModes(alphaFunc);
96 
97     osg::BlendFunc* blendFunc = new osg::BlendFunc;
98     blendFunc->setSource(osg::BlendFunc::SRC_ALPHA);
99     blendFunc->setDestination(osg::BlendFunc::ONE_MINUS_SRC_ALPHA);
100     stateSet->setAttributeAndModes(blendFunc);
101 
102     stateSet->setMode(GL_FOG, osg::StateAttribute::OFF);
103     stateSet->setMode(GL_LIGHTING, osg::StateAttribute::OFF);
104     stateSet->setMode(GL_CULL_FACE, osg::StateAttribute::OFF);
105     stateSet->setMode(GL_DEPTH_TEST, osg::StateAttribute::OFF);
106 
107     osg::Geode* geode = new osg::Geode;
108     stateSet = geode->getOrCreateStateSet();
109 
110     stateSet->setRenderBinDetails(-6, "RenderBin");
111 
112     // set up the sun-state
113     osg::Texture2D* texture = SGLoadTexture2D("sun.png", options.get());
114     stateSet->setTextureAttributeAndModes(0, texture);
115 
116     // Build scenegraph
117     sun_cl = new osg::Vec4Array;
118     sun_cl->push_back(osg::Vec4(1, 1, 1, 1));
119 
120     scene_cl = new osg::Vec4Array;
121     scene_cl->push_back(osg::Vec4(1, 1, 1, 1));
122 
123     osg::Vec3Array* sun_vl = new osg::Vec3Array;
124     sun_vl->push_back(osg::Vec3(-sun_size, 0, -sun_size));
125     sun_vl->push_back(osg::Vec3(sun_size, 0, -sun_size));
126     sun_vl->push_back(osg::Vec3(-sun_size, 0, sun_size));
127     sun_vl->push_back(osg::Vec3(sun_size, 0, sun_size));
128 
129     osg::Vec2Array* sun_tl = new osg::Vec2Array;
130     sun_tl->push_back(osg::Vec2(0, 0));
131     sun_tl->push_back(osg::Vec2(1, 0));
132     sun_tl->push_back(osg::Vec2(0, 1));
133     sun_tl->push_back(osg::Vec2(1, 1));
134 
135     osg::Geometry* geometry = new osg::Geometry;
136     geometry->setUseDisplayList(false);
137     geometry->setVertexArray(sun_vl);
138     geometry->setColorArray(sun_cl.get(), osg::Array::BIND_OVERALL);
139     geometry->setNormalBinding(osg::Geometry::BIND_OFF);
140     geometry->setTexCoordArray(0, sun_tl, osg::Array::BIND_PER_VERTEX);
141     geometry->addPrimitiveSet(new osg::DrawArrays(GL_TRIANGLE_STRIP, 0, 4));
142     geode->addDrawable(geometry);
143 
144     sun_transform->addChild( geode );
145 
146     // set up the inner-halo state
147     geode = new osg::Geode;
148     stateSet = geode->getOrCreateStateSet();
149     stateSet->setRenderBinDetails(-7, "RenderBin");
150 
151     texture = SGLoadTexture2D("inner_halo.png", options.get());
152     stateSet->setTextureAttributeAndModes(0, texture);
153 
154     // Build ssg structure
155     ihalo_cl = new osg::Vec4Array;
156     ihalo_cl->push_back(osg::Vec4(1, 1, 1, 1));
157 
158     float ihalo_size = sun_size * 2.0;
159     osg::Vec3Array* ihalo_vl = new osg::Vec3Array;
160     ihalo_vl->push_back(osg::Vec3(-ihalo_size, 0, -ihalo_size));
161     ihalo_vl->push_back(osg::Vec3(ihalo_size, 0, -ihalo_size));
162     ihalo_vl->push_back(osg::Vec3(-ihalo_size, 0, ihalo_size));
163     ihalo_vl->push_back(osg::Vec3(ihalo_size, 0, ihalo_size));
164 
165     osg::Vec2Array* ihalo_tl = new osg::Vec2Array;
166     ihalo_tl->push_back(osg::Vec2(0, 0));
167     ihalo_tl->push_back(osg::Vec2(1, 0));
168     ihalo_tl->push_back(osg::Vec2(0, 1));
169     ihalo_tl->push_back(osg::Vec2(1, 1));
170 
171     geometry = new osg::Geometry;
172     geometry->setUseDisplayList(false);
173     geometry->setVertexArray(ihalo_vl);
174     geometry->setColorArray(ihalo_cl.get(), osg::Array::BIND_OVERALL);
175     geometry->setNormalBinding(osg::Geometry::BIND_OFF);
176     geometry->setTexCoordArray(0, ihalo_tl, osg::Array::BIND_PER_VERTEX);
177     geometry->addPrimitiveSet(new osg::DrawArrays(GL_TRIANGLE_STRIP, 0, 4));
178     geode->addDrawable(geometry);
179 
180     sun_transform->addChild( geode );
181 
182     // set up the outer halo state
183 
184     geode = new osg::Geode;
185     stateSet = geode->getOrCreateStateSet();
186     stateSet->setRenderBinDetails(-8, "RenderBin");
187 
188     texture = SGLoadTexture2D("outer_halo.png", options.get());
189     stateSet->setTextureAttributeAndModes(0, texture);
190 
191     // Build ssg structure
192     ohalo_cl = new osg::Vec4Array;
193     ohalo_cl->push_back(osg::Vec4(1, 1, 1, 1));
194 
195     double ohalo_size = sun_size * 8.0;
196     osg::Vec3Array* ohalo_vl = new osg::Vec3Array;
197     ohalo_vl->push_back(osg::Vec3(-ohalo_size, 0, -ohalo_size));
198     ohalo_vl->push_back(osg::Vec3(ohalo_size, 0, -ohalo_size));
199     ohalo_vl->push_back(osg::Vec3(-ohalo_size, 0, ohalo_size));
200     ohalo_vl->push_back(osg::Vec3(ohalo_size, 0, ohalo_size));
201 
202     osg::Vec2Array* ohalo_tl = new osg::Vec2Array;
203     ohalo_tl->push_back(osg::Vec2(0, 0));
204     ohalo_tl->push_back(osg::Vec2(1, 0));
205     ohalo_tl->push_back(osg::Vec2(0, 1));
206     ohalo_tl->push_back(osg::Vec2(1, 1));
207 
208     geometry = new osg::Geometry;
209     geometry->setUseDisplayList(false);
210     geometry->setVertexArray(ohalo_vl);
211     geometry->setColorArray(ohalo_cl.get(), osg::Array::BIND_OVERALL);
212     geometry->setNormalBinding(osg::Geometry::BIND_OFF);
213     geometry->setTexCoordArray(0, ohalo_tl, osg::Array::BIND_PER_VERTEX);
214     geometry->addPrimitiveSet(new osg::DrawArrays(GL_TRIANGLE_STRIP, 0, 4));
215     geode->addDrawable(geometry);
216 
217     sun_transform->addChild( geode );
218 
219 
220     // set up the brilliance state
221 
222     geode = new osg::Geode;
223     stateSet = geode->getOrCreateStateSet();
224     stateSet->setRenderBinDetails(-9, "RenderBin");
225 
226 
227     texture = SGLoadTexture2D("brilliance.png", options.get());
228     stateSet->setTextureAttributeAndModes(0, texture);
229 
230 
231     // Build ssg structure
232     brilliance_cl = new osg::Vec4Array;
233     brilliance_cl->push_back(osg::Vec4(1, 1, 1, 1));
234 
235     double brilliance_size = sun_size * 12.0;
236     osg::Vec3Array* brilliance_vl = new osg::Vec3Array;
237     brilliance_vl->push_back(osg::Vec3(-brilliance_size, 0, -brilliance_size));
238     brilliance_vl->push_back(osg::Vec3(brilliance_size, 0, -brilliance_size));
239     brilliance_vl->push_back(osg::Vec3(-brilliance_size, 0, brilliance_size));
240     brilliance_vl->push_back(osg::Vec3(brilliance_size, 0, brilliance_size));
241 
242     osg::Vec2Array* brilliance_tl = new osg::Vec2Array;
243     brilliance_tl->push_back(osg::Vec2(0, 0));
244     brilliance_tl->push_back(osg::Vec2(1, 0));
245     brilliance_tl->push_back(osg::Vec2(0, 1));
246     brilliance_tl->push_back(osg::Vec2(1, 1));
247 
248     geometry = new osg::Geometry;
249     geometry->setUseDisplayList(false);
250     geometry->setVertexArray(brilliance_vl);
251     geometry->setColorArray(brilliance_cl.get(), osg::Array::BIND_OVERALL);
252     geometry->setNormalBinding(osg::Geometry::BIND_OFF);
253     geometry->setTexCoordArray(0, brilliance_tl, osg::Array::BIND_PER_VERTEX);
254     geometry->addPrimitiveSet(new osg::DrawArrays(GL_TRIANGLE_STRIP, 0, 4));
255     geode->addDrawable(geometry);
256 
257     sun_transform->addChild( geode );
258 
259     // force a repaint of the sun colors with arbitrary defaults
260     repaint( 0.0, 1.0 );
261 
262     return sun_transform.get();
263 }
264 
265 
266 // repaint the sun colors based on current value of sun_angle in
267 // degrees relative to verticle
268 // 0 degrees = high noon
269 // 90 degrees = sun rise/set
270 // 180 degrees = darkest midnight
repaint(double sun_angle,double new_visibility)271 bool SGSun::repaint( double sun_angle, double new_visibility ) {
272 
273     if ( visibility != new_visibility ) {
274         if (new_visibility < 100.0) new_visibility = 100.0;
275         else if (new_visibility > 45000.0) new_visibility = 45000.0;
276         visibility = new_visibility;
277         sun_exp2_punch_through = 2.0/log(visibility);
278     }
279 
280     if ( prev_sun_angle != sun_angle ) {
281         prev_sun_angle = sun_angle;
282 
283         // determine how much aerosols are in the air (rough guess)
284         double aerosol_factor;
285         if ( visibility < 100 ) {
286             aerosol_factor = 8000;
287         }
288         else {
289             aerosol_factor = 80.5 / log( visibility / 99.9 );
290         }
291 
292         // get environmental data from property tree or use defaults
293         double rel_humidity, density_avg;
294 
295         if ( !env_node ) {
296             rel_humidity = 0.5;
297             density_avg = 0.7;
298         }
299         else {
300             rel_humidity = env_node->getFloatValue( "relative-humidity" );
301             density_avg =  env_node->getFloatValue( "atmosphere/density-tropo-avg" );
302         }
303 
304         // ok, now let's go and generate the sun and scene color
305         osg::Vec4 i_halo_color, o_halo_color, scene_color, sun_color, brilliance_color;
306 
307         // Some comments:
308         // * When the sunangle changes, light has to travel a longer
309         //   distance through the atmosphere. So it's scattered more due
310         //   to raleigh scattering, which affects blue more than green
311         //   light.
312         // * Red is almost not scattered and effectively only get's
313         //   touched when the sun is near the horizon.
314         // * Visability also affects suncolor inasmuch as more particles
315         //   are in the air that cause more scattering.
316         // * We base our calculation on the halo's color, which is most
317         //   scattered.
318         double red_scat_f, red_scat_corr_f, green_scat_f, blue_scat_f;
319 
320         // Red - is almost not scattered
321         // Lambda is 700 nm
322 
323         red_scat_f = (aerosol_factor * path_distance * density_avg)/5E+07;
324         red_scat_corr_f = sun_exp2_punch_through / (1 - red_scat_f);
325         sun_color[0] = 1;
326         scene_color[0] = 1 - red_scat_f;
327 
328         // Green - 546.1 nm
329         green_scat_f = (aerosol_factor * path_distance * density_avg)/8.8938E+06;
330         sun_color[1] = 1 - green_scat_f * red_scat_corr_f;
331         scene_color[1] = 1 - green_scat_f;
332 
333         // Blue - 435.8 nm
334         blue_scat_f = (aerosol_factor * path_distance * density_avg)/3.607E+06;
335         sun_color[2] = 1 - blue_scat_f * red_scat_corr_f;
336         scene_color[2] = 1 - blue_scat_f;
337 
338         // Alpha
339         sun_color[3] = 1;
340         scene_color[3] = 1;
341 
342         // Now that we have the color calculated
343         // let's consider the saturation which is produced by mie scattering
344         double saturation = 1 - ( rel_humidity / 200 );
345         scene_color[1] += (( 1 - saturation ) * ( 1 - scene_color[1] ));
346         scene_color[2] += (( 1 - saturation ) * ( 1 - scene_color[2] ));
347 
348         if (sun_color[0] > 1.0) sun_color[0] = 1.0;
349         if (sun_color[0] < 0.0) sun_color[0] = 0.0;
350         if (sun_color[1] > 1.0) sun_color[1] = 1.0;
351         if (sun_color[1] < 0.0) sun_color[1] = 0.0;
352         if (sun_color[2] > 1.0) sun_color[2] = 1.0;
353         if (sun_color[2] < 0.0) sun_color[2] = 0.0;
354 
355         if (scene_color[0] > 1.0) scene_color[0] = 1.0;
356         if (scene_color[0] < 0.0) scene_color[0] = 0.0;
357         if (scene_color[1] > 1.0) scene_color[1] = 1.0;
358         if (scene_color[1] < 0.0) scene_color[1] = 0.0;
359         if (scene_color[2] > 1.0) scene_color[2] = 1.0;
360         if (scene_color[2] < 0.0) scene_color[2] = 0.0;
361 
362         double scene_f = 0.5 * (1 / (1 - red_scat_f));
363         double sun_f = 1.0 - scene_f;
364         i_halo_color[0] = sun_f * sun_color[0] + scene_f * scene_color[0];
365         i_halo_color[1] = sun_f * sun_color[1] + scene_f * scene_color[1];
366         i_halo_color[2] = sun_f * sun_color[2] + scene_f * scene_color[2];
367         i_halo_color[3] = 1;
368 
369         o_halo_color[0] = 0.2 * sun_color[0] + 0.8 * scene_color[0];
370         o_halo_color[1] = 0.2 * sun_color[1] + 0.8 * scene_color[1];
371         o_halo_color[2] = 0.2 * sun_color[2] + 0.8 * scene_color[2];
372         o_halo_color[3] = blue_scat_f;
373         if ((visibility < 10000) && (blue_scat_f > 1)) {
374             o_halo_color[3] = 2 - blue_scat_f;
375         }
376         if (o_halo_color[3] > 1) o_halo_color[3] = 1;
377         if (o_halo_color[3] < 0) o_halo_color[3] = 0;
378 
379 	brilliance_color[0] = i_halo_color[0];
380 	brilliance_color[1] = i_halo_color[1];
381 	brilliance_color[2] = i_halo_color[2];
382 
383 	double norm = (i_halo_color[0] * i_halo_color[0] + i_halo_color[1] * i_halo_color[1] + i_halo_color[2] * i_halo_color[2])/1.732;
384 
385 	brilliance_color[3] = pow(norm, 6.0);
386 	if (brilliance_color[3] < 0.0) {brilliance_color[3] = 0.0;}
387 
388 
389 
390         gamma_correct_rgb( i_halo_color._v );
391         gamma_correct_rgb( o_halo_color._v );
392         gamma_correct_rgb( scene_color._v );
393         gamma_correct_rgb( sun_color._v );
394         gamma_correct_rgb( brilliance_color._v );
395 
396 	if (sun_angle >91.0 * 3.1415/180.0 + horizon_angle)
397 		{
398 		sun_color[3] = 0;
399 		o_halo_color[3]=0;
400 		i_halo_color[3]=0;
401 		brilliance_color[3]=0;
402 		}
403 
404         (*sun_cl)[0] = sun_color;
405         sun_cl->dirty();
406         (*scene_cl)[0] = scene_color;
407         scene_cl->dirty();
408         (*ihalo_cl)[0] = i_halo_color;
409         ihalo_cl->dirty();
410         (*ohalo_cl)[0] = o_halo_color;
411         ohalo_cl->dirty();
412 	(*brilliance_cl)[0] = brilliance_color;
413         brilliance_cl->dirty();
414     }
415 
416     return true;
417 }
418 
419 
420 // reposition the sun at the specified right ascension and
421 // declination, offset by our current position (p) so that it appears
422 // fixed at a great distance from the viewer.  Also add in an optional
423 // rotation (i.e. for the current time of day.)
424 // Then calculate stuff needed for the sun-coloring
reposition(double rightAscension,double declination,double sun_dist,double lat,double alt_asl,double sun_angle)425 bool SGSun::reposition( double rightAscension, double declination,
426                         double sun_dist, double lat, double alt_asl, double sun_angle)
427 {
428     // GST - GMT sidereal time
429     osg::Matrix T2, RA, DEC;
430     // xglRotatef( ((SGD_RADIANS_TO_DEGREES * rightAscension)- 90.0),
431     //             0.0, 0.0, 1.0);
432     RA.makeRotate(rightAscension - 90*SGD_DEGREES_TO_RADIANS, osg::Vec3(0, 0, 1));
433     // xglRotatef((SGD_RADIANS_TO_DEGREES * declination), 1.0, 0.0, 0.0);
434     DEC.makeRotate(declination, osg::Vec3(1, 0, 0));
435     // xglTranslatef(0,sun_dist);
436     T2.makeTranslate(osg::Vec3(0, sun_dist, 0));
437     sun_transform->setMatrix(T2*DEC*RA);
438     // Suncolor related things:
439     if ( prev_sun_angle != sun_angle ) {
440       if ( sun_angle == 0 ) sun_angle = 0.1;
441          const double r_earth_pole = 6356752.314;
442          const double r_tropo_pole = 6356752.314 + 8000;
443          const double epsilon_earth2 = 6.694380066E-3;
444          const double epsilon_tropo2 = 9.170014946E-3;
445 
446          double r_tropo = r_tropo_pole / sqrt ( 1 - ( epsilon_tropo2 * pow ( cos( lat ), 2 )));
447          double r_earth = r_earth_pole / sqrt ( 1 - ( epsilon_earth2 * pow ( cos( lat ), 2 )));
448 
449          double position_radius = r_earth + alt_asl;
450 
451          double gamma =  SG_PI - sun_angle;
452          double sin_beta =  ( position_radius * sin ( gamma )  ) / r_tropo;
453          if (sin_beta > 1.0) sin_beta = 1.0;
454          double alpha =  SG_PI - gamma - asin( sin_beta );
455 
456          // OK, now let's calculate the distance the light travels
457          path_distance = sqrt( pow( position_radius, 2 ) + pow( r_tropo, 2 )
458                         - ( 2 * position_radius * r_tropo * cos( alpha ) ));
459 
460          double alt_half = sqrt( pow ( r_tropo, 2 ) + pow( path_distance / 2, 2 ) - r_tropo * path_distance * cos( asin( sin_beta )) ) - r_earth;
461 
462          if ( alt_half < 0.0 ) alt_half = 0.0;
463 
464 	 //angle at which the sun is visible below the horizon
465 	 horizon_angle = acos(min(r_earth/position_radius, 1.0));
466 
467          // Push the data to the property tree, so it can be used in the enviromental code
468          if ( env_node ){
469             env_node->setDoubleValue( "atmosphere/altitude-troposphere-top", r_tropo - r_earth );
470             env_node->setDoubleValue( "atmosphere/altitude-half-to-sun", alt_half );
471       }
472     }
473 
474     return true;
475 }
476 
477 SGVec4f
get_color()478 SGSun::get_color()
479 {
480     return SGVec4f((*sun_cl)[0][0], (*sun_cl)[0][1], (*sun_cl)[0][2], (*sun_cl)[0][3]);
481 }
482 
483 SGVec4f
get_scene_color()484 SGSun::get_scene_color()
485 {
486     return SGVec4f((*scene_cl)[0][0], (*scene_cl)[0][1], (*scene_cl)[0][2], (*scene_cl)[0][3]);
487 }
488