-// sunpos.c (taken from XEarth)
+// sunpos.cxx (adapted from XEarth)
// kirk johnson
// july 1993
//
// with your calculator, third edition_, peter duffett-smith,
// cambridge university press, 1988.)
//
-// RCS $Id$
-//
// Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson
//
// Parts of the source code (as marked) are:
#include <stdio.h>
#include <time.h>
-#include <Astro/orbits.hxx>
+//#include <Astro/orbits.hxx>
+#include <Astro/solarsystem.hxx>
#include <Include/fg_constants.h>
#include <Main/views.hxx>
-#include <Math/fg_geodesy.h>
+#include <Math/fg_geodesy.hxx>
#include <Math/mat3.h>
-#include <Math/polar3d.h>
-#include <Math/vector.h>
+#include <Math/point3d.hxx>
+#include <Math/polar3d.hxx>
+#include <Math/vector.hxx>
#include <Scenery/scenery.hxx>
#include "fg_time.hxx"
#include "sunpos.hxx"
+extern SolarSystem *solarSystem;
#undef E
/* double *delta; (return) declination */
double sin_e, cos_e;
+ double sin_l, cos_l;
sin_e = sin(MeanObliquity);
cos_e = cos(MeanObliquity);
+ sin_l = sin(lambda);
+ cos_l = cos(lambda);
- *alpha = atan2(sin(lambda)*cos_e - tan(beta)*sin_e, cos(lambda));
- *delta = asin(sin(beta)*cos_e + cos(beta)*sin_e*sin(lambda));
+ *alpha = atan2(sin_l*cos_e - tan(beta)*sin_e, cos_l);
+ *delta = asin(sin(beta)*cos_e + cos(beta)*sin_e*sin_l);
}
/* lambda = sun_ecliptic_longitude(ssue); */
/* ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta); */
- ecliptic_to_equatorial (solarPosition.lonSun, 0.0, &alpha, &delta);
-
+ //ecliptic_to_equatorial (solarPosition.lonSun, 0.0, &alpha, &delta);
+
+ /* **********************************************************************
+ * NOTE: in the next function, each time the sun's position is updated, the
+ * the sun's longitude is returned from solarSystem->sun. Note that the
+ * sun's position is updated at a much higher frequency than the rate at
+ * which the solar system's rebuilds occur. This is not a problem, however,
+ * because the fgSunPosition we're talking about here concerns the changing
+ * position of the sun due to the daily rotation of the earth.
+ * The ecliptic longitude, however, represents the position of the sun with
+ * respect to the stars, and completes just one cycle over the course of a
+ * year. Its therefore pretty safe to update the sun's longitude only once
+ * every ten minutes. (Comment added by Durk Talsma).
+ ************************************************************************/
+
+ ecliptic_to_equatorial( SolarSystem::theSolarSystem->getSun()->getLon(),
+ 0.0, &alpha, &delta );
tmp = alpha - (FG_2PI/24)*GST(ssue);
if (tmp < -FG_PI) {
do tmp += FG_2PI;
}
+/* given a particular time expressed in side real time at prime
+ * meridian (GST), compute position on the earth (lat, lon) such that
+ * sun is directly overhead. (lat, lon are reported in radians */
+
+static void fgSunPositionGST(double gst, double *lon, double *lat) {
+ /* time_t ssue; seconds since unix epoch */
+ /* double *lat; (return) latitude */
+ /* double *lon; (return) longitude */
+
+ /* double lambda; */
+ double alpha, delta;
+ double tmp;
+
+ /* lambda = sun_ecliptic_longitude(ssue); */
+ /* ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta); */
+ //ecliptic_to_equatorial (solarPosition.lonSun, 0.0, &alpha, &delta);
+ ecliptic_to_equatorial( SolarSystem::theSolarSystem->getSun()->getLon(),
+ 0.0, &alpha, &delta );
+
+// tmp = alpha - (FG_2PI/24)*GST(ssue);
+ tmp = alpha - (FG_2PI/24)*gst;
+ if (tmp < -FG_PI) {
+ do tmp += FG_2PI;
+ while (tmp < -FG_PI);
+ } else if (tmp > FG_PI) {
+ do tmp -= FG_2PI;
+ while (tmp < -FG_PI);
+ }
+
+ *lon = tmp;
+ *lat = delta;
+}
+
+
// update the cur_time_params structure with the current sun position
void fgUpdateSunPos( void ) {
fgLIGHT *l;
fgTIME *t;
fgVIEW *v;
MAT3vec nup, nsun, v0;
- fgPolarPoint3d p;
+ Point3D p, rel_sunpos;
+ double dot, east_dot;
double sun_gd_lat, sl_radius;
double ntmp;
printf(" Updating Sun position\n");
- fgSunPosition(t->cur_time, &l->sun_lon, &sun_gd_lat);
- fgSunPosition(t->cur_time, &l->sun_lon, &sun_gd_lat);
+ // (not sure why there was two)
+ // fgSunPosition(t->cur_time, &l->sun_lon, &sun_gd_lat);
+ fgSunPositionGST(t->gst, &l->sun_lon, &sun_gd_lat);
fgGeodToGeoc(sun_gd_lat, 0.0, &sl_radius, &l->sun_gc_lat);
- p.lon = l->sun_lon;
- p.lat = l->sun_gc_lat;
- p.radius = sl_radius;
+ p.setvals( l->sun_lon, l->sun_gc_lat, sl_radius );
l->fg_sunpos = fgPolarToCart3d(p);
printf( " t->cur_time = %ld\n", t->cur_time);
sun_gd_lat, l->sun_gc_lat);
// I think this will work better for generating the sun light vector
- l->sun_vec[0] = l->fg_sunpos.x;
- l->sun_vec[1] = l->fg_sunpos.y;
- l->sun_vec[2] = l->fg_sunpos.z;
+ l->sun_vec[0] = l->fg_sunpos.x();
+ l->sun_vec[1] = l->fg_sunpos.y();
+ l->sun_vec[2] = l->fg_sunpos.z();
MAT3_NORMALIZE_VEC(l->sun_vec, ntmp);
MAT3_SCALE_VEC(l->sun_vec_inv, l->sun_vec, -1.0);
l->sun_vec[3] = 0.0;
l->sun_vec_inv[3] = 0.0;
- printf(" l->sun_vec = %.2f %.2f %.2f\n", l->sun_vec[0], l->sun_vec[1],
- l->sun_vec[2]);
+ // printf(" l->sun_vec = %.2f %.2f %.2f\n", l->sun_vec[0], l->sun_vec[1],
+ // l->sun_vec[2]);
// calculate the sun's relative angle to local up
MAT3_COPY_VEC(nup, v->local_up);
- nsun[0] = l->fg_sunpos.x;
- nsun[1] = l->fg_sunpos.y;
- nsun[2] = l->fg_sunpos.z;
+ nsun[0] = l->fg_sunpos.x();
+ nsun[1] = l->fg_sunpos.y();
+ nsun[2] = l->fg_sunpos.z();
MAT3_NORMALIZE_VEC(nup, ntmp);
MAT3_NORMALIZE_VEC(nsun, ntmp);
l->sun_angle = acos(MAT3_DOT_PRODUCT(nup, nsun));
- printf(" SUN ANGLE relative to current location = %.3f rads.\n",
- l->sun_angle);
+ // printf(" SUN ANGLE relative to current location = %.3f rads.\n",
+ // l->sun_angle);
// calculate vector to sun's position on the earth's surface
- v->to_sun[0] = l->fg_sunpos.x - (v->view_pos.x + scenery.center.x);
- v->to_sun[1] = l->fg_sunpos.y - (v->view_pos.y + scenery.center.y);
- v->to_sun[2] = l->fg_sunpos.z - (v->view_pos.z + scenery.center.z);
+ rel_sunpos = l->fg_sunpos - (v->view_pos + scenery.center);
+ v->to_sun[0] = rel_sunpos.x();
+ v->to_sun[1] = rel_sunpos.y();
+ v->to_sun[2] = rel_sunpos.z();
// printf( "Vector to sun = %.2f %.2f %.2f\n",
// v->to_sun[0], v->to_sun[1], v->to_sun[2]);
// make a vector to the current view position
- MAT3_SET_VEC(v0, v->view_pos.x, v->view_pos.y, v->view_pos.z);
+ MAT3_SET_VEC(v0, v->view_pos.x(), v->view_pos.y(), v->view_pos.z());
// Given a vector from the view position to the point on the
// earth's surface the sun is directly over, map into onto the
// v->surface_to_sun[0], v->surface_to_sun[1], v->surface_to_sun[2]);
// printf("Should be close to zero = %.2f\n",
// MAT3_DOT_PRODUCT(v->local_up, v->surface_to_sun));
+
+ // calculate the angle between v->surface_to_sun and
+ // v->surface_east. We do this so we can sort out the acos()
+ // ambiguity. I wish I could think of a more efficient way ... :-(
+ east_dot = MAT3_DOT_PRODUCT(v->surface_to_sun, v->surface_east);
+ // printf(" East dot product = %.2f\n", east_dot);
+
+ // calculate the angle between v->surface_to_sun and
+ // v->surface_south. this is how much we have to rotate the sky
+ // for it to align with the sun
+ dot = MAT3_DOT_PRODUCT(v->surface_to_sun, v->surface_south);
+ // printf(" Dot product = %.2f\n", dot);
+ if ( east_dot >= 0 ) {
+ l->sun_rotation = acos(dot);
+ } else {
+ l->sun_rotation = -acos(dot);
+ }
+ // printf(" Sky needs to rotate = %.3f rads = %.1f degrees.\n",
+ // angle, angle * RAD_TO_DEG); */
}
// $Log$
+// Revision 1.14 1998/10/17 01:34:32 curt
+// C++ ifying ...
+//
+// Revision 1.13 1998/10/16 00:56:12 curt
+// Converted to Point3D class.
+//
+// Revision 1.12 1998/09/15 04:27:50 curt
+// Changes for new astro code.
+//
+// Revision 1.11 1998/08/12 21:13:22 curt
+// Optimizations by Norman Vine.
+//
+// Revision 1.10 1998/07/22 21:45:39 curt
+// fg_time.cxx: Removed call to ctime() in a printf() which should be harmless
+// but seems to be triggering a bug.
+// light.cxx: Added code to adjust fog color based on sunrise/sunset effects
+// and view orientation. This is an attempt to match the fog color to the
+// sky color in the center of the screen. You see discrepancies at the
+// edges, but what else can be done?
+// sunpos.cxx: Caculate local direction to sun here. (what compass direction
+// do we need to face to point directly at sun)
+//
+// Revision 1.9 1998/07/08 14:48:39 curt
+// polar3d.h renamed to polar3d.hxx
+//
// Revision 1.8 1998/05/02 01:53:18 curt
// Fine tuning mktime() support because of varying behavior on different
// platforms.
projects / flightgear.git / blobdiff