From: curt Date: Mon, 22 Mar 1999 12:08:57 +0000 (+0000) Subject: Initial revision. X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=50aabc3bc9953cc633532a42df1100242b829c51;p=flightgear.git Initial revision. --- diff --git a/Time/moonpos.cxx b/Time/moonpos.cxx new file mode 100644 index 000000000..7247cbebf --- /dev/null +++ b/Time/moonpos.cxx @@ -0,0 +1,608 @@ +// moonpos.cxx (basically, this is a slightly modified version of the 'sunpos.cxx' file, adapted from XEarth) + +// kirk johnson +// july 1993 +// +// code for calculating the position on the earth's surface for which +// the moon is directly overhead (adapted from _practical astronomy +// with your calculator, third edition_, peter duffett-smith, +// cambridge university press, 1988.) +// +// Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson +// +// Parts of the source code (as marked) are: +// Copyright (C) 1989, 1990, 1991 by Jim Frost +// Copyright (C) 1992 by Jamie Zawinski +// +// Permission to use, copy, modify and freely distribute xearth for +// non-commercial and not-for-profit purposes is hereby granted +// without fee, provided that both the above copyright notice and this +// permission notice appear in all copies and in supporting +// documentation. +// +// The author makes no representations about the suitability of this +// software for any purpose. It is provided "as is" without express or +// implied warranty. +// +// THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +// INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, +// IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT +// OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM +// LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, +// NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN +// CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. +// +// $Id$ +// (Log is kept at end of this file) + + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include "Include/compiler.h" +#ifdef FG_HAVE_STD_INCLUDES +# include +# include +# include +#else +# include +# include +# include +#endif + + +//#include +#include +#include +#include +#include
+#include +#include +#include +#include +#include +#include + +#include "fg_time.hxx" +#include "moonpos.hxx" + +extern SolarSystem *solarSystem; + +#undef E + + +/* + * the epoch upon which these astronomical calculations are based is + * 1990 january 0.0, 631065600 seconds since the beginning of the + * "unix epoch" (00:00:00 GMT, Jan. 1, 1970) + * + * given a number of seconds since the start of the unix epoch, + * DaysSinceEpoch() computes the number of days since the start of the + * astronomical epoch (1990 january 0.0) + */ + +#define EpochStart (631065600) +#define DaysSinceEpoch(secs) (((secs)-EpochStart)*(1.0/(24*3600))) + +/* + * assuming the apparent orbit of the moon about the earth is circular, + * the rate at which the orbit progresses is given by RadsPerDay -- + * FG_2PI radians per orbit divided by 365.242191 days per year: + */ + +#define RadsPerDay (FG_2PI/365.242191) + +/* + * details of moon's apparent orbit at epoch 1990.0 (after + * duffett-smith, table 6, section 46) + * + * Epsilon_g (ecliptic longitude at epoch 1990.0) 279.403303 degrees + * OmegaBar_g (ecliptic longitude of perigee) 282.768422 degrees + * Eccentricity (eccentricity of orbit) 0.016713 + */ + +#define Epsilon_g (279.403303*(FG_2PI/360)) +#define OmegaBar_g (282.768422*(FG_2PI/360)) +#define Eccentricity (0.016713) + +/* + * MeanObliquity gives the mean obliquity of the earth's axis at epoch + * 1990.0 (computed as 23.440592 degrees according to the method given + * in duffett-smith, section 27) + */ +#define MeanObliquity (23.440592*(FG_2PI/360)) + +/* static double solve_keplers_equation(double); */ +/* static double moon_ecliptic_longitude(time_t); */ +static void ecliptic_to_equatorial(double, double, double *, double *); +static double julian_date(int, int, int); +static double GST(time_t); + +/* + * solve Kepler's equation via Newton's method + * (after duffett-smith, section 47) + */ +/* +static double solve_keplers_equation(double M) { + double E; + double delta; + + E = M; + while (1) { + delta = E - Eccentricity*sin(E) - M; + if (fabs(delta) <= 1e-10) break; + E -= delta / (1 - Eccentricity*cos(E)); + } + + return E; +} +*/ + + +/* compute ecliptic longitude of moon (in radians) (after + * duffett-smith, section 47) */ +/* +static double moon_ecliptic_longitude(time_t ssue) { + // time_t ssue; // seconds since unix epoch + double D, N; + double M_moon, E; + double v; + + D = DaysSinceEpoch(ssue); + + N = RadsPerDay * D; + N = fmod(N, FG_2PI); + if (N < 0) N += FG_2PI; + + M_moon = N + Epsilon_g - OmegaBar_g; + if (M_moon < 0) M_moon += FG_2PI; + + E = solve_keplers_equation(M_moon); + v = 2 * atan(sqrt((1+Eccentricity)/(1-Eccentricity)) * tan(E/2)); + + return (v + OmegaBar_g); +} +*/ + + +/* convert from ecliptic to equatorial coordinates (after + * duffett-smith, section 27) */ + +static void ecliptic_to_equatorial(double lambda, double beta, + double *alpha, double *delta) { + /* double lambda; ecliptic longitude */ + /* double beta; ecliptic latitude */ + /* double *alpha; (return) right ascension */ + /* 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_l*cos_e - tan(beta)*sin_e, cos_l); + *delta = asin(sin(beta)*cos_e + cos(beta)*sin_e*sin_l); +} + + +/* computing julian dates (assuming gregorian calendar, thus this is + * only valid for dates of 1582 oct 15 or later) (after duffett-smith, + * section 4) */ + +static double julian_date(int y, int m, int d) { + /* int y; year (e.g. 19xx) */ + /* int m; month (jan=1, feb=2, ...) */ + /* int d; day of month */ + + int A, B, C, D; + double JD; + + /* lazy test to ensure gregorian calendar */ + if (y < 1583) { + FG_LOG( FG_EVENT, FG_ALERT, + "WHOOPS! Julian dates only valid for 1582 oct 15 or later" ); + } + + if ((m == 1) || (m == 2)) { + y -= 1; + m += 12; + } + + A = y / 100; + B = 2 - A + (A / 4); + C = (int)(365.25 * y); + D = (int)(30.6001 * (m + 1)); + + JD = B + C + D + d + 1720994.5; + + return JD; +} + + +/* compute greenwich mean sidereal time (GST) corresponding to a given + * number of seconds since the unix epoch (after duffett-smith, + * section 12) */ +static double GST(time_t ssue) { + /* time_t ssue; seconds since unix epoch */ + + double JD; + double T, T0; + double UT; + struct tm *tm; + + tm = gmtime(&ssue); + + JD = julian_date(tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday); + T = (JD - 2451545) / 36525; + + T0 = ((T + 2.5862e-5) * T + 2400.051336) * T + 6.697374558; + + T0 = fmod(T0, 24.0); + if (T0 < 0) T0 += 24; + + UT = tm->tm_hour + (tm->tm_min + tm->tm_sec / 60.0) / 60.0; + + T0 += UT * 1.002737909; + T0 = fmod(T0, 24.0); + if (T0 < 0) T0 += 24; + + return T0; +} + + +/* given a particular time (expressed in seconds since the unix + * epoch), compute position on the earth (lat, lon) such that moon is + * directly overhead. (lat, lon are reported in radians */ + +void fgMoonPosition(time_t ssue, 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 = moon_ecliptic_longitude(ssue); */ + /* ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta); */ + //ecliptic_to_equatorial (solarPosition.lonMoon, 0.0, &alpha, &delta); + + /* ********************************************************************** + * NOTE: in the next function, each time the moon's position is updated, the + * the moon's longitude is returned from solarSystem->moon. Note that the + * moon'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 fgMoonPosition we're talking about here concerns the changing + * position of the moon due to the daily rotation of the earth. + * The ecliptic longitude, however, represents the position of the moon with + * respect to the stars, and completes just one cycle over the course of a + * year. Its therefore pretty safe to update the moon's longitude only once + * every ten minutes. (Comment added by Durk Talsma). + ************************************************************************/ + + ecliptic_to_equatorial( SolarSystem::theSolarSystem->getMoon()->getLon(), + 0.0, &alpha, &delta ); + tmp = alpha - (FG_2PI/24)*GST(ssue); + 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; +} + + +/* given a particular time expressed in side real time at prime + * meridian (GST), compute position on the earth (lat, lon) such that + * moon is directly overhead. (lat, lon are reported in radians */ + +static void fgMoonPositionGST(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 = moon_ecliptic_longitude(ssue); */ + /* ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta); */ + //ecliptic_to_equatorial (solarPosition.lonMoon, 0.0, &alpha, &delta); + ecliptic_to_equatorial( SolarSystem::theSolarSystem->getMoon()->getLon(), + SolarSystem::theSolarSystem->getMoon()->getLat(), + &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 moon position +void fgUpdateMoonPos( void ) { + fgLIGHT *l; + fgTIME *t; + FGView *v; + MAT3vec nup, nmoon, v0, surface_to_moon; + Point3D p, rel_moonpos; + double dot, east_dot; + double moon_gd_lat, sl_radius; + double ntmp; + + l = &cur_light_params; + t = &cur_time_params; + v = ¤t_view; + + FG_LOG( FG_EVENT, FG_INFO, " Updating Moon position" ); + + // (not sure why there was two) + // fgMoonPosition(t->cur_time, &l->moon_lon, &moon_gd_lat); + fgMoonPositionGST(t->gst, &l->moon_lon, &moon_gd_lat); + + fgGeodToGeoc(moon_gd_lat, 0.0, &sl_radius, &l->moon_gc_lat); + + p = Point3D( l->moon_lon, l->moon_gc_lat, sl_radius ); + l->fg_moonpos = fgPolarToCart3d(p); + + FG_LOG( FG_EVENT, FG_INFO, " t->cur_time = " << t->cur_time ); + FG_LOG( FG_EVENT, FG_INFO, + " Moon Geodetic lat = " << moon_gd_lat + << " Geocentric lat = " << l->moon_gc_lat ); + + // I think this will work better for generating the moon light vector + l->moon_vec[0] = l->fg_moonpos.x(); + l->moon_vec[1] = l->fg_moonpos.y(); + l->moon_vec[2] = l->fg_moonpos.z(); + MAT3_NORMALIZE_VEC(l->moon_vec, ntmp); + MAT3_SCALE_VEC(l->moon_vec_inv, l->moon_vec, -1.0); + + // make sure these are directional light sources only + l->moon_vec[3] = 0.0; + l->moon_vec_inv[3] = 0.0; + + // printf(" l->moon_vec = %.2f %.2f %.2f\n", l->moon_vec[0], l->moon_vec[1], + // l->moon_vec[2]); + + // calculate the moon's relative angle to local up + MAT3_COPY_VEC(nup, v->get_local_up()); + nmoon[0] = l->fg_moonpos.x(); + nmoon[1] = l->fg_moonpos.y(); + nmoon[2] = l->fg_moonpos.z(); + MAT3_NORMALIZE_VEC(nup, ntmp); + MAT3_NORMALIZE_VEC(nmoon, ntmp); + + l->moon_angle = acos(MAT3_DOT_PRODUCT(nup, nmoon)); + // printf(" MOON ANGLE relative to current location = %.3f rads.\n", + // l->moon_angle); + + // calculate vector to moon's position on the earth's surface + rel_moonpos = l->fg_moonpos - (v->get_view_pos() + scenery.center); + v->set_to_moon( rel_moonpos.x(), rel_moonpos.y(), rel_moonpos.z() ); + // printf( "Vector to moon = %.2f %.2f %.2f\n", + // v->to_moon[0], v->to_moon[1], v->to_moon[2]); + + // make a vector to the current view position + Point3D view_pos = v->get_view_pos(); + MAT3_SET_VEC(v0, view_pos.x(), view_pos.y(), view_pos.z()); + + // Given a vector from the view position to the point on the + // earth's surface the moon is directly over, map into onto the + // local plane representing "horizontal". + map_vec_onto_cur_surface_plane( v->get_local_up(), v0, v->get_to_moon(), + surface_to_moon ); + MAT3_NORMALIZE_VEC(surface_to_moon, ntmp); + v->set_surface_to_moon( surface_to_moon[0], surface_to_moon[1], + surface_to_moon[2] ); + // printf("Surface direction to moon is %.2f %.2f %.2f\n", + // v->surface_to_moon[0], v->surface_to_moon[1], v->surface_to_moon[2]); + // printf("Should be close to zero = %.2f\n", + // MAT3_DOT_PRODUCT(v->local_up, v->surface_to_moon)); + + // calculate the angle between v->surface_to_moon 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( surface_to_moon, v->get_surface_east() ); + // printf(" East dot product = %.2f\n", east_dot); + + // calculate the angle between v->surface_to_moon and + // v->surface_south. this is how much we have to rotate the sky + // for it to align with the moon + dot = MAT3_DOT_PRODUCT( surface_to_moon, v->get_surface_south() ); + // printf(" Dot product = %.2f\n", dot); + if ( east_dot >= 0 ) { + l->moon_rotation = acos(dot); + } else { + l->moon_rotation = -acos(dot); + } + // printf(" Sky needs to rotate = %.3f rads = %.1f degrees.\n", + // angle, angle * RAD_TO_DEG); */ +} + + +// $Log$ +// Revision 1.1 1999/03/22 12:08:57 curt +// Initial revision. +// +// Revision 1.19 1999/01/07 20:25:37 curt +// Portability changes and updates from Bernie Bright. +// +// Revision 1.18 1998/12/09 18:50:36 curt +// Converted "class fgVIEW" to "class FGView" and updated to make data +// members private and make required accessor functions. +// +// Revision 1.17 1998/11/09 23:41:53 curt +// Log message clean ups. +// +// Revision 1.16 1998/11/07 19:07:14 curt +// Enable release builds using the --without-logging option to the configure +// script. Also a couple log message cleanups, plus some C to C++ comment +// conversion. +// +// Revision 1.15 1998/10/18 01:17:24 curt +// Point3D tweaks. +// +// 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 moonrise/moonset 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? +// moonpos.cxx: Caculate local direction to moon here. (what compass direction +// do we need to face to point directly at moon) +// +// 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. +// +// Revision 1.7 1998/04/30 12:36:05 curt +// C++-ifying a couple source files. +// +// Revision 1.6 1998/04/28 01:22:18 curt +// Type-ified fgTIME and fgVIEW. +// +// Revision 1.5 1998/04/26 05:10:05 curt +// "struct fgLIGHT" -> "fgLIGHT" because fgLIGHT is typedef'd. +// +// Revision 1.4 1998/04/25 22:06:34 curt +// Edited cvs log messages in source files ... bad bad bad! +// +// Revision 1.3 1998/04/25 20:24:03 curt +// Cleaned up initialization sequence to eliminate interdependencies +// between moon position, lighting, and view position. This creates a +// valid single pass initialization path. +// +// Revision 1.2 1998/04/24 00:52:31 curt +// Wrapped "#include " in "#ifdef HAVE_CONFIG_H" +// Fog color fixes. +// Separated out lighting calcs into their own file. +// +// Revision 1.1 1998/04/22 13:24:07 curt +// C++ - ifiing the code a bit. +// Starting to reorginize some of the lighting calcs to use a table lookup. +// +// Revision 1.27 1998/04/03 22:12:57 curt +// Converting to Gnu autoconf system. +// Centralized time handling differences. +// +// Revision 1.26 1998/02/23 19:08:00 curt +// Incorporated Durk's Astro/ tweaks. Includes unifying the moon position +// calculation code between moon display, and other FG sections that use this +// for things like lighting. +// +// Revision 1.25 1998/02/09 15:07:53 curt +// Minor tweaks. +// +// Revision 1.24 1998/01/27 00:48:07 curt +// Incorporated Paul Bleisch's new debug message +// system and commandline/config file processing code. +// +// Revision 1.23 1998/01/19 19:27:21 curt +// Merged in make system changes from Bob Kuehne +// This should simplify things tremendously. +// +// Revision 1.22 1998/01/19 18:40:40 curt +// Tons of little changes to clean up the code and to remove fatal errors +// when building with the c++ compiler. +// +// Revision 1.21 1997/12/30 23:10:19 curt +// Calculate lighting parameters here. +// +// Revision 1.20 1997/12/30 22:22:43 curt +// Further integration of event manager. +// +// Revision 1.19 1997/12/30 20:47:59 curt +// Integrated new event manager with subsystem initializations. +// +// Revision 1.18 1997/12/23 04:58:40 curt +// Tweaked the sky coloring a bit to build in structures to allow finer rgb +// control. +// +// Revision 1.17 1997/12/15 23:55:08 curt +// Add xgl wrappers for debugging. +// Generate terrain normals on the fly. +// +// Revision 1.16 1997/12/11 04:43:57 curt +// Fixed moon vector and lighting problems. I thing the moon is now lit +// correctly. +// +// Revision 1.15 1997/12/10 22:37:55 curt +// Prepended "fg" on the name of all global structures that didn't have it yet. +// i.e. "struct WEATHER {}" became "struct fgWEATHER {}" +// +// Revision 1.14 1997/12/09 04:25:39 curt +// Working on adding a global lighting params structure. +// +// Revision 1.13 1997/11/25 19:25:42 curt +// Changes to integrate Durk's moon/moon code updates + clean up. +// +// Revision 1.12 1997/11/15 18:15:39 curt +// Reverse direction of moon vector, so object normals can be more normal. +// +// Revision 1.11 1997/10/28 21:07:21 curt +// Changed GLUT/ -> Main/ +// +// Revision 1.10 1997/09/13 02:00:09 curt +// Mostly working on stars and generating sidereal time for accurate star +// placement. +// +// Revision 1.9 1997/09/05 14:17:31 curt +// More tweaking with stars. +// +// Revision 1.8 1997/09/05 01:36:04 curt +// Working on getting stars right. +// +// Revision 1.7 1997/09/04 02:17:40 curt +// Shufflin' stuff. +// +// Revision 1.6 1997/08/27 03:30:37 curt +// Changed naming scheme of basic shared structures. +// +// Revision 1.5 1997/08/22 21:34:41 curt +// Doing a bit of reorganizing and house cleaning. +// +// Revision 1.4 1997/08/19 23:55:09 curt +// Worked on better simulating real lighting. +// +// Revision 1.3 1997/08/13 20:23:49 curt +// The interface to moonpos now updates a global structure rather than returning +// current moon position. +// +// Revision 1.2 1997/08/06 00:24:32 curt +// Working on correct real time moon lighting. +// +// Revision 1.1 1997/08/01 15:27:56 curt +// Initial revision. +// diff --git a/Time/moonpos.hxx b/Time/moonpos.hxx new file mode 100644 index 000000000..dad36e880 --- /dev/null +++ b/Time/moonpos.hxx @@ -0,0 +1,60 @@ +/* + * moonpos.hxx (taken from XEarth) + * kirk johnson + * july 1993 + * + * code for calculating the position on the earth's surface for which + * the sun is directly overhead (adapted from _practical astronomy + * 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: + * Copyright (C) 1989, 1990, 1991 by Jim Frost + * Copyright (C) 1992 by Jamie Zawinski + * + * Permission to use, copy, modify and freely distribute xearth for + * non-commercial and not-for-profit purposes is hereby granted + * without fee, provided that both the above copyright notice and this + * permission notice appear in all copies and in supporting + * documentation. + * + * The author makes no representations about the suitability of this + * software for any purpose. It is provided "as is" without express or + * implied warranty. + * + * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, + * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT + * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM + * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, + * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN + * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + + +#ifndef _MOONPOS_HXX +#define _MOONPOS_HXX + + +#ifndef __cplusplus +# error This library requires C++ +#endif + +#include "Include/compiler.h" +#ifdef FG_HAVE_STD_INCLUDES +# include +#else +# include +#endif + +/* update the cur_time_params structure with the current moon position */ +void fgUpdateMoonPos( void ); + +void fgMoonPosition(time_t ssue, double *lon, double *lat); + + +#endif /* _MOONPOS_H */