diff --git a/simgear/ephemeris/ephemeris.cxx b/simgear/ephemeris/ephemeris.cxx index f918c3f7..0cc9b996 100644 --- a/simgear/ephemeris/ephemeris.cxx +++ b/simgear/ephemeris/ephemeris.cxx @@ -68,7 +68,8 @@ SGEphemeris::~SGEphemeris( void ) { void SGEphemeris::update( double mjd, double lst, double lat ) { // update object positions our_sun->updatePosition( mjd ); - moon->updatePosition( mjd, lst, lat, our_sun ); + // moon->updatePositionTopo( mjd, lst, lat, our_sun ); + moon->updatePosition( mjd, our_sun ); mercury->updatePosition( mjd, our_sun ); venus->updatePosition( mjd, our_sun ); mars->updatePosition( mjd, our_sun ); diff --git a/simgear/ephemeris/ephemeris.hxx b/simgear/ephemeris/ephemeris.hxx index abe5a8a0..0da2d1cb 100644 --- a/simgear/ephemeris/ephemeris.hxx +++ b/simgear/ephemeris/ephemeris.hxx @@ -51,7 +51,7 @@ * Written by Durk Talsma and Curtis Olson * * - * Introduction + * Introduction * * The SGEphemeris class computes and stores the positions of the Sun, * the Moon, the planets, and the brightest stars. These positions @@ -146,6 +146,11 @@ public: return moon->getDeclination(); } + /** @return the geocentric distance to the Moon in unit of its semi-mayor axis. */ + inline double getMoonDistanceInMayorAxis() const { + return moon->getDistanceInMayorAxis(); + } + /** @return the numbers of defined planets. */ inline int getNumPlanets() const { return nplanets; } diff --git a/simgear/ephemeris/moonpos.cxx b/simgear/ephemeris/moonpos.cxx index 94f1b36c..3948fb5d 100644 --- a/simgear/ephemeris/moonpos.cxx +++ b/simgear/ephemeris/moonpos.cxx @@ -1,9 +1,9 @@ /************************************************************************** * moonpos.cxx - * Written by Durk Talsma. Originally started October 1997, for distribution - * with the FlightGear project. Version 2 was written in August and - * September 1998. This code is based upon algorithms and data kindly - * provided by Mr. Paul Schlyter. (pausch@saaf.se). + * Written by Durk Talsma. Originally started October 1997, for distribution + * with the FlightGear project. Version 2 was written in August and + * September 1998. This code is based upon algorithms and data kindly + * provided by Mr. Paul Schlyter. (pausch@saaf.se). * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public @@ -37,10 +37,10 @@ /************************************************************************* * MoonPos::MoonPos(double mjd) - * Public constructor for class MoonPos. Initializes the orbital elements and + * Public constructor for class MoonPos. Initializes the orbital elements and * sets up the moon texture. * Argument: The current time. - * the hard coded orbital elements for MoonPos are passed to + * the hard coded orbital elements for MoonPos are passed to * CelestialBody::CelestialBody(); ************************************************************************/ MoonPos::MoonPos(double mjd) : @@ -70,21 +70,23 @@ MoonPos::~MoonPos() /***************************************************************************** - * void MoonPos::updatePosition(double mjd, Star *ourSun) - * this member function calculates the actual topocentric position (i.e.) - * the position of the moon as seen from the current position on the surface - * of the moon. + * void MoonPos::updatePositionTopo(double mjd, double lst, double + * lat, Star *ourSun) this member function calculates the actual + * topocentric position (i.e.) the position of the moon as seen from + * the current position on the surface of the earth. This include + * parallax effects, the moon appears on different stars background + * for different observers on the surface of the earth. ****************************************************************************/ -void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) +void MoonPos::updatePositionTopo(double mjd, double lst, double lat, Star *ourSun) { - double + double eccAnom, ecl, actTime, xv, yv, v, r, xh, yh, zh, zg, xe, Ls, Lm, D, F, mpar, gclat, rho, HA, g, geoRa, geoDec, cosN, sinN, cosvw, sinvw, sinvw_cosi, cosecl, sinecl, rcoslatEcl, FlesstwoD, MlesstwoD, twoD, twoM, twolat, alpha; - + double max_loglux = -0.504030345621; double min_loglux = -4.39964634562; double conv = 1.0319696543787917; // The log foot-candle to log lux conversion factor. @@ -99,7 +101,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) yv = a * (sqrt(1.0 - e*e) * sin(eccAnom)); v = atan2(yv, xv); // the moon's true anomaly r = sqrt (xv*xv + yv*yv); // and its distance - + // repetitive calculations, minimised for speed cosN = cos(N); sinN = sin(N); @@ -118,7 +120,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) lonEcl = atan2 (yh, xh); latEcl = atan2(zh, sqrt(xh*xh + yh*yh)); - /* Calculate a number of perturbation, i.e. disturbances caused by the + /* Calculate a number of perturbation, i.e. disturbances caused by the * gravitational influence of the sun and the other major planets. * The largest of these even have a name */ Ls = ourSun->getM() + ourSun->getw(); @@ -130,7 +132,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) twoM = 2 * M; FlesstwoD = F - twoD; MlesstwoD = M - twoD; - + lonEcl += SGD_DEGREES_TO_RADIANS * (-1.274 * sin(MlesstwoD) +0.658 * sin(twoD) -0.186 * sin(ourSun->getM()) @@ -154,12 +156,13 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) -0.46 * cos(twoD) ); distance = r; + distance_in_a = r/a; // SG_LOG(SG_GENERAL, SG_INFO, "Running moon update"); rcoslatEcl = r * cos(latEcl); xg = cos(lonEcl) * rcoslatEcl; yg = sin(lonEcl) * rcoslatEcl; zg = r * sin(latEcl); - + xe = xg; ye = yg * cosecl -zg * sinecl; ze = yg * sinecl +zg * cosecl; @@ -167,17 +170,17 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) geoRa = atan2(ye, xe); geoDec = atan2(ze, sqrt(xe*xe + ye*ye)); - /* SG_LOG( SG_GENERAL, SG_INFO, - "(geocentric) geoRa = (" << (SGD_RADIANS_TO_DEGREES * geoRa) + /* SG_LOG( SG_GENERAL, SG_INFO, + "(geocentric) geoRa = (" << (SGD_RADIANS_TO_DEGREES * geoRa) << "), geoDec= (" << (SGD_RADIANS_TO_DEGREES * geoDec) << ")" ); */ - // Given the moon's geocentric ra and dec, calculate its + // Given the moon's geocentric ra and dec, calculate its // topocentric ra and dec. i.e. the position as seen from the // surface of the earth, instead of the center of the earth - // First calculate the moon's parallax, that is, the apparent size of the - // (equatorial) radius of the earth, as seen from the moon + // First calculate the moon's parallax, that is, the apparent size of the + // (equatorial) radius of the earth, as seen from the moon mpar = asin ( 1 / r); // SG_LOG( SG_GENERAL, SG_INFO, "r = " << r << " mpar = " << mpar ); // SG_LOG( SG_GENERAL, SG_INFO, "lat = " << f->get_Latitude() ); @@ -188,12 +191,12 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) rho = 0.99883 + 0.00167 * cos(twolat); // SG_LOG( SG_GENERAL, SG_INFO, "rho = " << rho ); - + if (geoRa < 0) geoRa += SGD_2PI; - + HA = lst - (3.8197186 * geoRa); - /* SG_LOG( SG_GENERAL, SG_INFO, "t->getLst() = " << t->getLst() + /* SG_LOG( SG_GENERAL, SG_INFO, "t->getLst() = " << t->getLst() << " HA = " << HA ); */ g = atan (tan(gclat) / cos ((HA / 3.8197186))); @@ -212,8 +215,140 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun) // << SGD_RADIANS_TO_DEGREES * (geoDec - declination) << endl; } - /* SG_LOG( SG_GENERAL, SG_INFO, - "Ra = (" << (SGD_RADIANS_TO_DEGREES *rightAscension) + /* SG_LOG( SG_GENERAL, SG_INFO, + "Ra = (" << (SGD_RADIANS_TO_DEGREES *rightAscension) + << "), Dec= (" << (SGD_RADIANS_TO_DEGREES *declination) << ")" ); */ + + // Moon age and phase calculation + age = lonEcl - ourSun->getlonEcl(); + phase = (1 - cos(age)) / 2; + + // The log of the illuminance of the moon outside the atmosphere. + // This is the base 10 log of equation 20 from Krisciunas K. and Schaefer B.E. + // (1991). A model of the brightness of moonlight, Publ. Astron. Soc. Pacif. + // 103(667), 1033-1039 (DOI: http://dx.doi.org/10.1086/132921). + alpha = SGD_RADIANS_TO_DEGREES * SGMiscd::normalizeAngle(age + SGMiscd::pi()); + log_I = -0.4 * (3.84 + 0.026*fabs(alpha) + 4e-9*pow(alpha, 4.0)); + + // Convert from foot-candles to lux. + log_I += conv; + + // The moon's illuminance factor, bracketed between 0 and 1. + I_factor = (log_I - max_loglux) / (max_loglux - min_loglux) + 1.0; + I_factor = SGMiscd::clip(I_factor, 0, 1); +} + + + +/***************************************************************************** + * void MoonPos::updatePosition(double mjd, Star *ourSun) this member + * function calculates the geocentric position (i.e.) the position of + * the moon as seen from the center of earth. As such, it does not + * include any parallax effects. These are taken into account during + * the rendering. + ****************************************************************************/ +void MoonPos::updatePosition(double mjd, Star *ourSun) +{ + double + eccAnom, ecl, actTime, + xv, yv, v, r, xh, yh, zh, zg, xe, + Ls, Lm, D, F, geoRa, geoDec, + cosN, sinN, cosvw, sinvw, sinvw_cosi, cosecl, sinecl, rcoslatEcl, + FlesstwoD, MlesstwoD, twoD, twoM, twolat, alpha; + + double max_loglux = -0.504030345621; + double min_loglux = -4.39964634562; + double conv = 1.0319696543787917; // The log foot-candle to log lux conversion factor. + updateOrbElements(mjd); + actTime = sgCalcActTime(mjd); + + // calculate the angle between ecliptic and equatorial coordinate system + // in Radians + ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 * actTime); + eccAnom = sgCalcEccAnom(M, e); // Calculate the eccentric anomaly + xv = a * (cos(eccAnom) - e); + yv = a * (sqrt(1.0 - e*e) * sin(eccAnom)); + v = atan2(yv, xv); // the moon's true anomaly + r = sqrt (xv*xv + yv*yv); // and its distance + + // repetitive calculations, minimised for speed + cosN = cos(N); + sinN = sin(N); + cosvw = cos(v+w); + sinvw = sin(v+w); + sinvw_cosi = sinvw * cos(i); + cosecl = cos(ecl); + sinecl = sin(ecl); + + // estimate the geocentric rectangular coordinates here + xh = r * (cosN * cosvw - sinN * sinvw_cosi); + yh = r * (sinN * cosvw + cosN * sinvw_cosi); + zh = r * (sinvw * sin(i)); + + // calculate the ecliptic latitude and longitude here + lonEcl = atan2 (yh, xh); + latEcl = atan2(zh, sqrt(xh*xh + yh*yh)); + + /* Calculate a number of perturbation, i.e. disturbances caused by the + * gravitational influence of the sun and the other major planets. + * The largest of these even have a name */ + Ls = ourSun->getM() + ourSun->getw(); + Lm = M + w + N; + D = Lm - Ls; + F = Lm - N; + + twoD = 2 * D; + twoM = 2 * M; + FlesstwoD = F - twoD; + MlesstwoD = M - twoD; + + lonEcl += SGD_DEGREES_TO_RADIANS * (-1.274 * sin(MlesstwoD) + +0.658 * sin(twoD) + -0.186 * sin(ourSun->getM()) + -0.059 * sin(twoM - twoD) + -0.057 * sin(MlesstwoD + ourSun->getM()) + +0.053 * sin(M + twoD) + +0.046 * sin(twoD - ourSun->getM()) + +0.041 * sin(M - ourSun->getM()) + -0.035 * sin(D) + -0.031 * sin(M + ourSun->getM()) + -0.015 * sin(2*F - twoD) + +0.011 * sin(M - 4*D) + ); + latEcl += SGD_DEGREES_TO_RADIANS * (-0.173 * sin(FlesstwoD) + -0.055 * sin(M - FlesstwoD) + -0.046 * sin(M + FlesstwoD) + +0.033 * sin(F + twoD) + +0.017 * sin(twoM + F) + ); + r += (-0.58 * cos(MlesstwoD) + -0.46 * cos(twoD) + ); + // from the orbital elements, the unit of the distance is in Earth radius, around 60. + distance = r; + distance_in_a = r/a; + // SG_LOG(SG_GENERAL, SG_INFO, "Running moon update"); + rcoslatEcl = r * cos(latEcl); + xg = cos(lonEcl) * rcoslatEcl; + yg = sin(lonEcl) * rcoslatEcl; + zg = r * sin(latEcl); + + xe = xg; + ye = yg * cosecl -zg * sinecl; + ze = yg * sinecl +zg * cosecl; + + geoRa = atan2(ye, xe); + geoDec = atan2(ze, sqrt(xe*xe + ye*ye)); + + if (geoRa < 0) + geoRa += SGD_2PI; + + rightAscension = geoRa; + declination = geoDec; + + + /* SG_LOG( SG_GENERAL, SG_INFO, + "Ra = (" << (SGD_RADIANS_TO_DEGREES *rightAscension) << "), Dec= (" << (SGD_RADIANS_TO_DEGREES *declination) << ")" ); */ // Moon age and phase calculation diff --git a/simgear/ephemeris/moonpos.hxx b/simgear/ephemeris/moonpos.hxx index b601d425..078a55f6 100644 --- a/simgear/ephemeris/moonpos.hxx +++ b/simgear/ephemeris/moonpos.hxx @@ -39,6 +39,7 @@ private: double xg, yg; // the moon's rectangular geocentric coordinates double ye, ze; // the moon's rectangular equatorial coordinates double distance; // the moon's distance to the earth + double distance_in_a; // the moon's distance to the earth in unit of its semi-mayor axis a double age; // the moon's age from 0 to 2pi double phase; // the moon's phase double log_I; // the moon's illuminance outside the atmosphere (logged) @@ -59,8 +60,9 @@ public: MoonPos(double mjd); MoonPos(); ~MoonPos(); - void updatePosition(double mjd, double lst, double lat, Star *ourSun); - // void newImage(); + void updatePositionTopo(double mjd, double lst, double lat, Star *ourSun); + void updatePosition(double mjd, Star *ourSun); + // void newImage(); double getM() const; double getw() const; double getxg() const; @@ -68,6 +70,7 @@ public: double getye() const; double getze() const; double getDistance() const; + double getDistanceInMayorAxis() const; double getAge() const; double getPhase() const; double getLogIlluminance() const; @@ -109,6 +112,11 @@ inline double MoonPos::getDistance() const return distance; } +inline double MoonPos::getDistanceInMayorAxis() const +{ + return distance_in_a; +} + inline double MoonPos::getAge() const { return age; diff --git a/simgear/math/simd.hxx b/simgear/math/simd.hxx index 162820ac..b5cab2e5 100644 --- a/simgear/math/simd.hxx +++ b/simgear/math/simd.hxx @@ -82,7 +82,7 @@ inline T magnitude2(const simd4_t& vi) { template inline simd4_t interpolate(T tau, const simd4_t& v1, const simd4_t& v2) { - return (T(1)-tau)*v1 + tau*v2; + return v1 + tau*(v2-v1); } template diff --git a/simgear/scene/sky/moon.cxx b/simgear/scene/sky/moon.cxx index 363c1335..665115e4 100644 --- a/simgear/scene/sky/moon.cxx +++ b/simgear/scene/sky/moon.cxx @@ -1,9 +1,9 @@ // moon.hxx -- model earth's moon // -// Written by Durk Talsma. Originally started October 1997, for distribution -// with the FlightGear project. Version 2 was written in August and -// September 1998. This code is based upon algorithms and data kindly -// provided by Mr. Paul Schlyter. (pausch@saaf.se). +// Written by Durk Talsma. Originally started October 1997, for distribution +// with the FlightGear project. Version 2 was written in August and +// September 1998. This code is based upon algorithms and data kindly +// provided by Mr. Paul Schlyter. (pausch@saaf.se). // // Separated out rendering pieces and converted to ssg by Curt Olson, // March 2000 @@ -69,7 +69,7 @@ SGMoon::~SGMoon( void ) { osg::Node* SGMoon::build( SGPath path, double moon_size ) { - osg::Node* orb = SGMakeSphere(moon_size, 15, 15); + osg::Node* orb = SGMakeSphere(moon_size, 40, 20); osg::StateSet* stateSet = orb->getOrCreateStateSet(); stateSet->setRenderBinDetails(-5, "RenderBin"); @@ -141,18 +141,18 @@ bool SGMoon::repaint( double moon_angle ) { prev_moon_angle = moon_angle; float moon_factor = 4*cos(moon_angle); - + if (moon_factor > 1) moon_factor = 1.0; if (moon_factor < -1) moon_factor = -1.0; moon_factor = moon_factor/2 + 0.5; - + osg::Vec4 color; color[1] = sqrt(moon_factor); color[0] = sqrt(color[1]); color[2] = moon_factor * moon_factor; color[2] *= color[2]; color[3] = 1.0; - + gamma_correct_rgb( color._v ); orb_material->setDiffuse(osg::Material::FRONT_AND_BACK, color); @@ -162,22 +162,74 @@ bool SGMoon::repaint( double moon_angle ) { // reposition the moon at the specified right ascension and -// declination, offset by our current position (p) so that it appears -// fixed at a great distance from the viewer. Also add in an optional -// rotation (i.e. for the current time of day.) -bool SGMoon::reposition( double rightAscension, double declination, - double moon_dist ) -{ - osg::Matrix T2, RA, DEC; +// declination from the center of Earth. Because the view is actually +// offset by our current position (p), we first evaluate our current +// position w.r.t the Moon and then shift to the articial center of +// earth before shifting to the rendered moon distance. This allows to +// implement any parallax effects. moon_dist_bare is expected to not +// change during the rendering, it gives us the normalisation factors +// between real distances and units used in the +// rendering. moon_dist_fact is any extra factors to put the moon +// further or closer. + +bool SGMoon::reposition( double rightAscension, double declination, + double moon_dist_bare, double moon_dist_factor, + double lst, double lat, double alt ) +{ + osg::Matrix TE, T2, RA, DEC; + + // semi mayor axis + const double moon_a_in_rearth = 60.266600; + // average (we could also account for equatorial streching like in + // oursun.cxx if required) + const double earth_radius_in_meters = 6371000.0; + //local hour angle in radians + double lha; + // in unit of the rendering + double moon_dist; + double earth_radius, viewer_radius; + double xp,yp,zp; + + //rendered earth radius according to what has been specified by + //moon_dist_bare + earth_radius = moon_dist_bare/moon_a_in_rearth; + + //how far are we from the center of Earth + viewer_radius = (1.0 + alt/earth_radius_in_meters)*earth_radius; + + // the local hour angle of the moon, .i.e. its angle with respect + // to the meridian of the viewer + lha = lst - rightAscension; + + // the shift vector of the observer w.r.t. earth center (funny + // convention on x?) + zp = viewer_radius * sin(lat); + yp = viewer_radius * cos(lat)*cos(lha); + xp = viewer_radius * cos(lat)*sin(-lha); + + //rotate along the z axis RA.makeRotate(rightAscension - 90.0 * SGD_DEGREES_TO_RADIANS, osg::Vec3(0, 0, 1)); - + //rotate along the rotated x axis DEC.makeRotate(declination, osg::Vec3(1, 0, 0)); + //move to the center of Earth + TE.makeTranslate(osg::Vec3(-xp,-yp,-zp)); + + //move the moon from the center of Earth to moon_dist + moon_dist = moon_dist_bare * moon_dist_factor; T2.makeTranslate(osg::Vec3(0, moon_dist, 0)); - moon_transform->setMatrix(T2*DEC*RA); + // cout << " viewer radius= " << viewer_radius << endl; + // cout << " xp yp zp= " << xp <<" " << yp << " " << zp << endl; + // cout << " lha= " << lha << endl; + // cout << " moon_dist_bare= " << moon_dist_bare << endl; + // cout << " moon_dist_factor= " << moon_dist_factor << endl; + // cout << " moon_dist= " << moon_dist << endl; + + moon_transform->setMatrix(T2*TE*DEC*RA); + return true; } diff --git a/simgear/scene/sky/moon.hxx b/simgear/scene/sky/moon.hxx index 5efd5575..95474ecc 100644 --- a/simgear/scene/sky/moon.hxx +++ b/simgear/scene/sky/moon.hxx @@ -65,11 +65,12 @@ public: bool repaint( double moon_angle ); // reposition the moon at the specified right ascension and - // declination, offset by our current position (p) so that it - // appears fixed at a great distance from the viewer. Also add in - // an optional rotation (i.e. for the current time of day.) + // declination, at moon_dist_bare*moon_dist_factor from the center + // of Earth. lst, lat and alt are need to estimate where is the + // center of Earth from the current view). bool reposition( double rightAscension, double declination, - double moon_dist ); + double moon_dist_bare, double moon_dist_factor, + double lst, double lat, double alt ); }; diff --git a/simgear/scene/sky/sky.cxx b/simgear/scene/sky/sky.cxx index 1c7f13e5..e2cb1efb 100644 --- a/simgear/scene/sky/sky.cxx +++ b/simgear/scene/sky/sky.cxx @@ -152,7 +152,7 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt double angleRad = SGMiscd::deg2rad(angle); SGVec3f zero_elev, view_up; - double lon, lat, alt; + double lon, lat, alt, lst; SGGeod geodZeroViewPos = SGGeod::fromGeodM(st.pos_geod, 0); zero_elev = toVec3f( SGVec3d::fromGeod(geodZeroViewPos) ); @@ -165,7 +165,9 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt lon = st.pos_geod.getLongitudeRad(); lat = st.pos_geod.getLatitudeRad(); alt = st.pos_geod.getElevationM(); - + // Local sidereal time + lst = angleRad + lon; + dome->reposition( zero_elev, alt, lon, lat, st.spin ); osg::Matrix m = osg::Matrix::rotate(angleRad, osg::Vec3(0, 0, -1)); @@ -178,7 +180,12 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt double moon_ra = eph.getMoonRightAscension(); double moon_dec = eph.getMoonDeclination(); - moon->reposition( moon_ra, moon_dec, st.moon_dist ); + double moon_r = eph.getMoonDistanceInMayorAxis(); + + //this allows to render the moon closer to the viewer when the + //moon is closer to the center of Earth, times any articial extra factors + double moon_dist_factor = moon_r * st.moon_dist_factor; + moon->reposition( moon_ra, moon_dec, st.moon_dist_bare, moon_dist_factor, lst, lat, alt ); for ( unsigned i = 0; i < cloud_layers.size(); ++i ) { if ( cloud_layers[i]->getCoverage() != SGCloudLayer::SG_CLOUD_CLEAR || diff --git a/simgear/scene/sky/sky.hxx b/simgear/scene/sky/sky.hxx index a07d9d0c..9799280c 100644 --- a/simgear/scene/sky/sky.hxx +++ b/simgear/scene/sky/sky.hxx @@ -63,7 +63,8 @@ struct SGSkyState double gst; //!< GMT side real time. double sun_dist; //!< the sun's distance from the current view point // (to keep it inside your view volume). - double moon_dist;//!< The moon's distance from the current view point. + double moon_dist_bare ;//!< The moon's semi-mayor axis in the rendering (constant) + double moon_dist_factor ;//!< Any factor that are needed to artificially change the moon distance double sun_angle; };