1 // tmp.cxx -- stuff I don't know what to do with at the moment
3 // Written by Curtis Olson, started July 2000.
5 // Copyright (C) 2000 Curtis L. Olson - http://www.flightgear.org/~curt
7 // This program is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU General Public License as
9 // published by the Free Software Foundation; either version 2 of the
10 // License, or (at your option) any later version.
12 // This program is distributed in the hope that it will be useful, but
13 // WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // General Public License for more details.
17 // You should have received a copy of the GNU General Public License
18 // along with this program; if not, write to the Free Software
19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 #include <simgear/math/SGMath.hxx>
29 #include <simgear/math/vector.hxx>
30 #include <simgear/misc/sg_path.hxx>
31 #include <simgear/timing/sg_time.hxx>
33 #include <Main/fg_props.hxx>
34 #include <Main/globals.hxx>
35 #include <Main/viewer.hxx>
36 #include <Scenery/scenery.hxx>
39 #include "sunsolver.hxx"
43 // periodic time updater wrapper
44 void fgUpdateLocalTime() {
45 static const SGPropertyNode *longitude
46 = fgGetNode("/position/longitude-deg");
47 static const SGPropertyNode *latitude
48 = fgGetNode("/position/latitude-deg");
50 SGPath zone( globals->get_fg_root() );
51 zone.append( "Timezone" );
53 SG_LOG(SG_GENERAL, SG_INFO, "updateLocal("
54 << longitude->getDoubleValue() * SGD_DEGREES_TO_RADIANS
56 << latitude->getDoubleValue() * SGD_DEGREES_TO_RADIANS
57 << ", " << zone.str() << ")");
58 globals->get_time_params()->updateLocal( longitude->getDoubleValue()
59 * SGD_DEGREES_TO_RADIANS,
60 latitude->getDoubleValue()
61 * SGD_DEGREES_TO_RADIANS,
66 // update the cur_time_params structure with the current sun position
67 void fgUpdateSunPos( void ) {
69 // This only works at lat,lon = 0,0
70 // need to find a way to get it working at other locations
72 FGLight *light = (FGLight *)(globals->get_subsystem("lighting"));
73 FGViewer *viewer = globals->get_current_view();
74 SGTime *time_now = globals->get_time_params();
76 SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" );
77 SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << time_now->getGst() );
79 double sun_lon, sun_lat;
80 fgSunPositionGST(time_now->getGst(), &sun_lon, &sun_lat);
81 light->set_sun_lon(sun_lon);
82 light->set_sun_lat(sun_lat);
84 // update the sun light vector
85 // calculations are in the horizontal normal plane: x-north, y-east, z-down
86 static SGQuatd q = SGQuatd::fromLonLat(SGGeod::fromRad(0,0));
89 SGGeod geodSunPos = SGGeod::fromRad(sun_lon, sun_lat);
90 SGQuatd sunOr = SGQuatd::fromLonLat(geodSunPos);
92 // scenery orientation
93 SGGeod geodViewPos = SGGeod::fromCart(viewer->getViewPosition());
94 SGQuatd hlOr = SGQuatd::fromLonLat(geodViewPos);
95 SGVec3d localAt = hlOr.backTransform(SGVec3d::e3());
97 // transpose the sun direction from (lat,lon) to (0,0)
98 SGVec3d transSunDir = (q*sunOr).transform(-localAt);
99 SGQuatd sunDirOr = SGQuatd::fromRealImag(0, transSunDir);
101 // transpose the calculated sun vector back to (lat,lon)
102 SGVec3d sunDirection = sunDirOr.transform(localAt);
103 light->set_sun_rotation( acos(sunDirection[1])-SGD_PI_2 );
104 light->set_sun_angle( acos(-sunDirection[2]) );
106 SGVec3d sunPos = SGVec3d::fromGeod(geodSunPos);
107 light->sun_vec() = SGVec4f(toVec3f(normalize(sunPos)), 0);
108 light->sun_vec_inv() = -light->sun_vec();
111 FGLight *l = (FGLight *)(globals->get_subsystem("lighting"));
112 SGTime *t = globals->get_time_params();
113 FGViewer *v = globals->get_current_view();
115 SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" );
116 SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t->getGst() );
120 fgSunPositionGST(t->getGst(), &sun_l, &sun_gd_lat);
121 l->set_sun_lon(sun_l);
122 l->set_sun_lat(sun_gd_lat);
123 SGVec3d sunpos(SGVec3d::fromGeod(SGGeod::fromRad(sun_l, sun_gd_lat)));
125 SG_LOG( SG_EVENT, SG_DEBUG, " t->cur_time = " << t->get_cur_time() );
126 SG_LOG( SG_EVENT, SG_DEBUG,
127 " Sun Geodetic lat = " << sun_gd_lat
128 << " Geodetic lat = " << sun_gd_lat );
130 // update the sun light vector
131 l->sun_vec() = SGVec4f(toVec3f(normalize(sunpos)), 0);
132 l->sun_vec_inv() = - l->sun_vec();
134 // calculate the sun's relative angle to local up
135 SGVec3d viewPos = v->get_view_pos();
136 SGQuatd hlOr = SGQuatd::fromLonLat(SGGeod::fromCart(viewPos));
137 SGVec3f world_up = toVec3f(hlOr.backTransform(-SGVec3d::e3()));
138 SGVec3f nsun = toVec3f(normalize(sunpos));
139 // cout << "nup = " << nup[0] << "," << nup[1] << ","
140 // << nup[2] << endl;
141 // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
142 // << nsun[2] << endl;
144 l->set_sun_angle( acos( dot ( world_up, nsun ) ) );
145 SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
146 << l->get_sun_angle() );
148 // calculate vector to sun's position on the earth's surface
149 SGVec3d rel_sunpos = sunpos - v->get_view_pos();
150 // vector in cartesian coordinates from current position to the
151 // postion on the earth's surface the sun is directly over
152 SGVec3f to_sun = toVec3f(rel_sunpos);
153 // printf( "Vector to sun = %.2f %.2f %.2f\n",
154 // v->to_sun[0], v->to_sun[1], v->to_sun[2]);
156 // Given a vector from the view position to the point on the
157 // earth's surface the sun is directly over, map into onto the
158 // local plane representing "horizontal".
160 // surface direction to go to head towards sun
161 SGVec3f surface_to_sun;
162 SGVec3f view_pos = toVec3f(v->get_view_pos());
163 sgmap_vec_onto_cur_surface_plane( world_up.data(), view_pos.data(),
164 to_sun.data(), surface_to_sun.data() );
165 surface_to_sun = normalize(surface_to_sun);
166 // cout << "(sg) Surface direction to sun is "
167 // << surface_to_sun[0] << ","
168 // << surface_to_sun[1] << ","
169 // << surface_to_sun[2] << endl;
170 // cout << "Should be close to zero = "
171 // << sgScalarProductVec3(nup, surface_to_sun) << endl;
173 // calculate the angle between surface_to_sun and
174 // v->get_surface_east(). We do this so we can sort out the
175 // acos() ambiguity. I wish I could think of a more efficient
177 SGVec3f surface_east(toVec3f(hlOr.backTransform(SGVec3d::e2())));
178 float east_dot = dot( surface_to_sun, surface_east );
179 // cout << " East dot product = " << east_dot << endl;
181 // calculate the angle between v->surface_to_sun and
182 // v->surface_south. this is how much we have to rotate the sky
183 // for it to align with the sun
184 SGVec3f surface_south(toVec3f(hlOr.backTransform(-SGVec3d::e1())));
185 float dot_ = dot( surface_to_sun, surface_south );
186 // cout << " Dot product = " << dot << endl;
189 SG_LOG( SG_ASTRO, SG_INFO,
190 "Dot product = " << dot_ << " is greater than 1.0" );
193 else if (dot_ < -1.0) {
194 SG_LOG( SG_ASTRO, SG_INFO,
195 "Dot product = " << dot_ << " is less than -1.0" );
199 if ( east_dot >= 0 ) {
200 l->set_sun_rotation( acos(dot_) );
202 l->set_sun_rotation( -acos(dot_) );
204 // cout << " Sky needs to rotate = " << angle << " rads = "
205 // << angle * SGD_RADIANS_TO_DEGREES << " degrees." << endl;