On Monday 27 April 2009 02:31:48 Davide Bettio wrote: > moved moonphase code to time dataengine from luna plasmoid/moonphase kicker > applet (as requested by petri). phases.cpp is missing but it comes from the > plasmoid.
I have made moonrise / moonset calculation to the code. There is also moon phase there (phase calculation seems simple after sun and moon values for rise / set are calculated). Currently phase is 0.0 - 1.0 (0.5 meaning full moon). How should we continue from here? Petri
Index: timeengine.cpp
===================================================================
--- timeengine.cpp (revision 959786)
+++ timeengine.cpp (working copy)
@@ -106,7 +106,7 @@
// this is relatively cheap - KSTZ::local() is cached
timezone = KSystemTimeZones::local().name();
} else {
- KTimeZone newTz = KSystemTimeZones::zone(tz);
+ KTimeZone newTz = KSystemTimeZones::zone(timezone);
if (!newTz.isValid()) {
return false;
}
@@ -118,7 +118,7 @@
QString trTimezone = i18n(timezone.toUtf8());
data[I18N_NOOP("Timezone")] = trTimezone;
-
+
QStringList tzParts = trTimezone.split("/");
data[I18N_NOOP("Timezone Continent")] = tzParts.value(0);
data[I18N_NOOP("Timezone City")] = tzParts.value(1);
@@ -126,7 +126,8 @@
if (args.count() > 0) {
static const QString latitude = I18N_NOOP("Latitude");
static const QString longitude = I18N_NOOP("Longitude");
-
+ bool sun = false, moon = false;
+
data[latitude] = args[latitude].toDouble();
data[longitude] = args[longitude].toDouble();
QDateTime dt = QDateTime::fromString(args[I18N_NOOP("DateTime")], Qt::ISODate);
@@ -137,10 +138,16 @@
data[I18N_NOOP("Time")] = dt.time();
}
if (args.contains(I18N_NOOP("Solar"))) {
+ sun = true;
+ }
+ if (args.contains(I18N_NOOP("Lunar"))) {
+ moon = true;
+ }
+ if (sun || moon) {
if (!solarPosition) {
solarPosition = new SolarPosition;
}
- solarPosition->appendData(data);
+ solarPosition->appendData(data, sun, moon);
}
if (args.contains(I18N_NOOP("Moon"))) {
appendMoonphase(data);
Index: solarposition.cpp
===================================================================
--- solarposition.cpp (revision 959786)
+++ solarposition.cpp (working copy)
@@ -20,48 +20,106 @@
#include <QDateTime>
/*
- * This class is ported from public domain javascript code:
+ * Mathematics, ideas, public domain code used for these classes from:
+ * http://www.stjarnhimlen.se/comp/tutorial.html
+ * http://www.stjarnhimlen.se/comp/riset.html
* http://www.srrb.noaa.gov/highlights/solarrise/azel.html
* http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html
- * Calculation details:
- * http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
+ * http://bodmas.org/astronomy/riset.html
+ * moontool.c by John Walker
+ * Wikipedia
*/
-class NOAASolarCalc
+class SolarSystemObject
{
public:
- static void calc(const QDateTime &dt, double longitude, double latitude,
- double zone, double *jd, double *century, double *eqTime,
- double *solarDec, double *azimuth, double *zenith);
- static double radToDeg(double angleRad);
- static double degToRad(double angleDeg);
- static double calcJD(double year, double month, double day);
- static QDateTime calcDateFromJD(double jd, double minutes, double zone);
- static double calcTimeJulianCent(double jd);
- static double calcJDFromJulianCent(double t);
- static double calcGeomMeanLongSun(double t);
- static double calcGeomMeanAnomalySun(double t);
- static double calcEccentricityEarthOrbit(double t);
- static double calcSunEqOfCenter(double t);
- static double calcSunTrueLong(double t);
- static double calcSunTrueAnomaly(double t);
- static double calcSunRadVector(double t);
- static double calcSunApparentLong(double t);
- static double calcMeanObliquityOfEcliptic(double t);
- static double calcObliquityCorrection(double t);
- static double calcSunRtAscension(double t);
- static double calcSunDeclination(double t);
- static double calcEquationOfTime(double t);
- static void calcAzimuthAndZenith(QDateTime now, double eqTime, double zone,
- double solarDec, double latitude, double longitude,
- double *zenith, double *azimuth);
- static double calcElevation(double zenith);
- static double calcHourAngle(double zenith, double solarDec, double latitude);
- static void calcTimeUTC(double zenith, bool rise, double *jd, double *minutes,
- double latitude, double longitude);
- static double calcSolNoonUTC(double t, double longitude);
+ SolarSystemObject(double latitude, double longitude);
+ virtual ~SolarSystemObject();
+
+ double meanLongitude() { return L; };
+ double meanAnomaly() { return M; };
+ double siderealTime();
+ double altitude() { return m_altitude; };
+ double azimuth() { return m_azimuth; };
+ double calcElevation();
+ QDateTime dateTime() { return m_local; };
+ double lambda() { return m_lambda; };
+ double eclipticLongitude() { return m_eclipticLongitude; };
+
+ virtual void calcForDateTime(const QDateTime& local, int offset);
+ QList< QPair<QDateTime, QDateTime> > timesForAngles(const QList<double>& angles,
+ int offset);
+
+ protected:
+ void calc();
+ virtual bool calcPerturbations(double*, double*, double*) { return false; };
+ virtual void rotate(double*, double*) { };
+ virtual void topocentricCorrection(double*, double*) { };
+
+ inline double rev(double x);
+ inline double asind(double x);
+ inline double sind(double x);
+ inline double cosd(double x);
+ inline double atand(double x);
+ inline double tand(double x);
+ inline double atan2d(double y, double x);
+ void toRectangular(double lo, double la, double r, double *x, double *y, double *z);
+ void toSpherical(double x, double y, double z, double *lo, double *la, double *r);
+ QPair<double, double> zeroPoints(QPointF p1, QPointF p2, QPointF p3);
+
+ double N;
+ double i;
+ double w;
+ double a;
+ double e;
+ double M;
+ double m_obliquity;
+
+ QDateTime m_utc;
+ QDateTime m_local;
+ double m_day;
+ double m_latitude;
+ double m_longitude;
+
+ double L;
+ double rad;
+ double RA;
+ double dec;
+ double HA;
+ double m_altitude;
+ double m_azimuth;
+ double m_eclipticLongitude;
+ double m_lambda;
};
+class Sun : public SolarSystemObject
+{
+ public:
+ Sun(double latitude, double longitude);
+
+ virtual void calcForDateTime(const QDateTime& local, int offset);
+
+ protected:
+ virtual void rotate(double*, double*);
+};
+
+class Moon : public SolarSystemObject
+{
+ public:
+ Moon(double latitude, double longitude, Sun *sun);
+
+ virtual void calcForDateTime(const QDateTime& local, int offset);
+ double phase();
+
+ protected:
+ virtual bool calcPerturbations(double *RA, double *dec, double *r);
+ virtual void rotate(double*, double*);
+ virtual void topocentricCorrection(double*, double*);
+
+ private:
+ Sun *m_sun;
+};
+
SolarPosition::SolarPosition()
{
}
@@ -70,649 +128,352 @@
{
}
-void SolarPosition::appendData(Plasma::DataEngine::Data &data)
+void SolarPosition::appendData(Plasma::DataEngine::Data &data, bool fillSun, bool fillMoon)
{
//QTime time;
//time.start();
double longitude = data["Longitude"].toDouble();
double latitude = data["Latitude"].toDouble();
- double zone = data["Offset"].toDouble() / -3600.0;
- QDateTime dt(data["Date"].toDate(), data["Time"].toTime());
-
- double jd;
- double century;
- double eqTime;
- double solarDec;
- double azimuth;
- double zenith;
-
- NOAASolarCalc::calc(dt, longitude, latitude, zone, &jd, ¢ury, &eqTime,
- &solarDec, &azimuth, &zenith);
- data["Equation of Time"] = eqTime;
- data["Solar Declination"] = solarDec;
- data["Azimuth"] = azimuth;
- data["Zenith"] = zenith;
- data["Corrected Elevation"] = NOAASolarCalc::calcElevation(zenith);
-
+ int offset = data["Offset"].toInt();
+ QDateTime local(data["Date"].toDate(), data["Time"].toTime());
const QString cacheKey = QString("%1|%2|%3")
- .arg(latitude).arg(longitude).arg(dt.date().toString(Qt::ISODate));
- if (!m_cache.contains(cacheKey)) {
- double minutes;
- double jd2;
-
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(90.833, true, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Apparent Sunrise"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(90.833, false, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Apparent Sunset"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
+ .arg(latitude).arg(longitude).arg(local.date().toString(Qt::ISODate));
+ Sun sun(latitude, longitude);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(90.0, true, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Sunrise"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(90.0, false, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Sunset"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
+ if (fillSun) {
+ sun.calcForDateTime(local, offset);
+ data["Sun Azimuth"] = sun.azimuth();
+ data["Sun Zenith"] = 90.0 - sun.altitude();
+ data["Sun Corrected Elevation"] = sun.calcElevation();
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(96.0, true, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Civil Dawn"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(96.0, false, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Civil Dusk"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
+ if (!m_cache.contains(cacheKey) || !m_cache[cacheKey].contains("Sunrise")) {
+ QList< QPair<QDateTime, QDateTime> > times = sun.timesForAngles(
+ QList<double>() << -0.833 << -6.0 << -12.0 << -18.0, offset);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(102.0, true, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Nautical Dawn"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(102.0, false, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Nautical Dusk"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
+ m_cache[cacheKey]["Sunrise"] = times[0].first;
+ m_cache[cacheKey]["Sunset"] = times[0].second;
+ m_cache[cacheKey]["Civil Dawn"] = times[1].first;
+ m_cache[cacheKey]["Civil Dusk"] = times[1].second;
+ m_cache[cacheKey]["Nautical Dawn"] = times[2].first;
+ m_cache[cacheKey]["Nautical Dusk"] = times[2].second;
+ m_cache[cacheKey]["Astronomical Dawn"] = times[3].first;
+ m_cache[cacheKey]["Astronomical Dusk"] = times[3].second;
+ }
+ }
+ if (fillMoon) {
+ Moon moon(latitude, longitude, &sun);
+ moon.calcForDateTime(local, offset);
+ data["Moon Azimuth"] = moon.azimuth();
+ data["Moon Zenith"] = 90 - moon.altitude();
+ data["Moon Corrected Elevation"] = moon.calcElevation();
+ data["Moon Phase"] = moon.phase();
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(108.0, true, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Astronomical Dawn"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
- jd2 = jd;
- NOAASolarCalc::calcTimeUTC(108.0, false, &jd2, &minutes, latitude, longitude);
- m_cache[cacheKey]["Astronomical Dusk"] = NOAASolarCalc::calcDateFromJD(jd2, minutes, zone);
-
- century = NOAASolarCalc::calcTimeJulianCent(jd);
- minutes = NOAASolarCalc::calcSolNoonUTC(century, longitude);
- dt = NOAASolarCalc::calcDateFromJD(jd, minutes, zone);
- NOAASolarCalc::calc(dt, longitude, latitude, zone, &jd, ¢ury, &eqTime,
- &solarDec, &azimuth, &zenith);
- m_cache[cacheKey]["Solar Noon"] = dt;
- m_cache[cacheKey]["Min Zenith"] = zenith;
- m_cache[cacheKey]["Max Corrected Elevation"] = NOAASolarCalc::calcElevation(zenith);
+ if (!m_cache.contains(cacheKey) || !m_cache[cacheKey].contains("Moonrise")) {
+ QList< QPair<QDateTime, QDateTime> > times = moon.timesForAngles(
+ QList<double>() << -0.833, offset);
+ m_cache[cacheKey]["Moonrise"] = times[0].first;
+ m_cache[cacheKey]["Moonset"] = times[0].second;
+ }
}
data.unite(m_cache[cacheKey]);
//kDebug() << "Calculation took:" << time.elapsed() << "ms";
//kDebug() << data;
}
-void NOAASolarCalc::calc(const QDateTime &dt, double longitude, double latitude,
- double zone, double *jd, double *century, double *eqTime,
- double *solarDec, double *azimuth, double *zenith)
+Sun::Sun(double latitude, double longitude)
+ : SolarSystemObject(latitude, longitude)
{
- double timenow = dt.time().hour() + dt.time().minute() / 60.0 +
- dt.time().second() / 3600.0 + zone;
- *jd = NOAASolarCalc::calcJD(dt.date().year(), dt.date().month(), dt.date().day());
- *century = NOAASolarCalc::calcTimeJulianCent(*jd + timenow / 24.0);
- //double earthRadVec = NOAASolarCalc::calcSunRadVector(*century);
- //double alpha = NOAASolarCalc::calcSunRtAscension(*century);
- *eqTime = NOAASolarCalc::calcEquationOfTime(*century);
- *solarDec = NOAASolarCalc::calcSunDeclination(*century);
- NOAASolarCalc::calcAzimuthAndZenith(dt, *eqTime, zone, *solarDec, latitude, longitude,
- zenith, azimuth);
}
-// Convert radian angle to degrees
-double NOAASolarCalc::radToDeg(double angleRad)
+void Sun::calcForDateTime(const QDateTime& local, int offset)
{
- return (180.0 * angleRad / M_PI);
+ SolarSystemObject::calcForDateTime(local, offset);
+
+ N = 0.0;
+ i = 0.0;
+ w = rev(282.9404 + 4.70935E-5 * m_day);
+ a = 1.0;
+ e = rev(0.016709 - 1.151E-9 * m_day);
+ M = rev(356.0470 + 0.9856002585 * m_day);
+
+ calc();
}
-// Convert degree angle to radians
-double NOAASolarCalc::degToRad(double angleDeg)
+void Sun::rotate(double* y, double* z)
{
- return (M_PI * angleDeg / 180.0);
+ *y *= cosd(m_obliquity);
+ *z *= sind(m_obliquity);
}
-//***********************************************************************/
-//* Name: calcJD */
-//* Type: Function */
-//* Purpose: Julian day from calendar day */
-//* Arguments: */
-//* year : 4 digit year */
-//* month: January = 1 */
-//* day : 1 - 31 */
-//* Return value: */
-//* The Julian day corresponding to the date */
-//* Note: */
-//* Number is returned for start of day. Fractional days should be */
-//* added later. */
-//***********************************************************************/
-double NOAASolarCalc::calcJD(double year, double month, double day)
+Moon::Moon(double latitude, double longitude, Sun *sun)
+ : SolarSystemObject(latitude, longitude)
+ , m_sun(sun)
{
- if (month <= 2) {
- year -= 1;
- month += 12;
- }
- double A = floor(year / 100.0);
- double B = 2 - A + floor(A / 4.0);
-
- return floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1))
- + day + B - 1524.5;
}
-//***********************************************************************/
-//* Name: calcDateFromJD */
-//* Type: Function */
-//* Purpose: Calendar date from Julian Day */
-//* Arguments: */
-//* jd : Julian Day */
-//* Return value: */
-//* QDateTime */
-//* Note: */
-//***********************************************************************/
-QDateTime NOAASolarCalc::calcDateFromJD(double jd, double minutes, double zone)
+void Moon::calcForDateTime(const QDateTime& local, int offset)
{
- QDateTime result;
-
- minutes -= 60 * zone;
- if (minutes > 1440.0) {
- minutes -= 1440.0;
- jd += 1.0;
+ if (m_sun->dateTime() != local) {
+ m_sun->calcForDateTime(local, offset);
}
- if (minutes < 0.0) {
- minutes += 1440.0;
- jd -= 1.0;
- }
- double floatHour = minutes / 60.0;
- double hour = floor(floatHour);
- double floatMinute = 60.0 * (floatHour - floor(floatHour));
- double minute = floor(floatMinute);
- double floatSec = 60.0 * (floatMinute - floor(floatMinute));
- double second = floor(floatSec + 0.5);
- result.setTime(QTime(hour, minute, second));
- double z = floor(jd + 0.5);
- double f = (jd + 0.5) - z;
- double A;
- if (z < 2299161.0) {
- A = z;
- } else {
- double alpha = floor((z - 1867216.25) / 36524.25);
- A = z + 1.0 + alpha - floor(alpha / 4.0);
- }
- double B = A + 1524.0;
- double C = floor((B - 122.1) / 365.25);
- double D = floor(365.25 * C);
- double E = floor((B - D) / 30.6001);
- double day = B - D - floor(30.6001 * E) + f;
- double month = (E < 14.0) ? E - 1 : E - 13.0;
- double year = (month > 2.0) ? C - 4716.0 : C - 4715.0;
- result.setDate(QDate(year, month, day));
+ SolarSystemObject::calcForDateTime(local, offset);
- return result;
-}
+ N = rev(125.1228 - 0.0529538083 * m_day);
+ i = 5.1454;
+ w = rev(318.0634 + 0.1643573223 * m_day);
+ a = 60.2666;
+ e = 0.054900;
+ M = rev(115.3654 + 13.0649929509 * m_day);
-//***********************************************************************/
-//* Name: calcTimeJulianCent */
-//* Type: Function */
-//* Purpose: convert Julian Day to centuries since J2000.0. */
-//* Arguments: */
-//* jd : the Julian Day to convert */
-//* Return value: */
-//* the T value corresponding to the Julian Day */
-//***********************************************************************/
-double NOAASolarCalc::calcTimeJulianCent(double jd)
-{
- return (jd - 2451545.0) / 36525.0;
+ calc();
}
-//***********************************************************************/
-//* Name: calcJDFromJulianCent */
-//* Type: Function */
-//* Purpose: convert centuries since J2000.0 to Julian Day. */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* the Julian Day corresponding to the t value */
-//***********************************************************************/
-double NOAASolarCalc::calcJDFromJulianCent(double t)
+bool Moon::calcPerturbations(double *lo, double *la, double *r)
{
- return t * 36525.0 + 2451545.0;
+ double Ms = m_sun->meanAnomaly();
+ double D = L - m_sun->meanLongitude();
+ double F = L - N;
+
+ *lo += -1.274 * sind(M - 2 * D)
+ +0.658 * sind(2 * D)
+ -0.186 * sind(Ms)
+ -0.059 * sind(2 * M - 2 * D)
+ -0.057 * sind(M - 2 * D + Ms)
+ +0.053 * sind(M + 2 * D)
+ +0.046 * sind(2 * D - Ms)
+ +0.041 * sind(M - Ms)
+ -0.035 * sind(D)
+ -0.031 * sind(M + Ms)
+ -0.015 * sind(2 * F - 2 * D)
+ +0.011 * sind(M - 4 * D);
+ *la += -0.173 * sind(F - 2 * D)
+ -0.055 * sind(M - F - 2 * D)
+ -0.046 * sind(M + F - 2 * D)
+ +0.033 * sind(F + 2 * D)
+ +0.017 * sind(2 * M + F);
+ *r += -0.58 * cosd(M - 2 * D)
+ -0.46 * cosd(2 * D);
+ return true;
}
-//***********************************************************************/
-//* Name: calcGeomMeanLongSun */
-//* Type: Function */
-//* Purpose: calculate the Geometric Mean Longitude of the Sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* the Geometric Mean Longitude of the Solar in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcGeomMeanLongSun(double t)
+void Moon::topocentricCorrection(double* RA, double* dec)
{
- double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
- while (L0 > 360.0) {
- L0 -= 360.0;
- }
- while(L0 < 0.0) {
- L0 += 360.0;
- }
- return L0; // in degrees
+ double HA = rev(siderealTime() - *RA);
+ double gclat = m_latitude - 0.1924 * sind(2 * m_latitude);
+ double rho = 0.99833 + 0.00167 * cosd(2 * m_latitude);
+ double mpar = asind(1 / rad);
+ double g = atand(tand(gclat) / cosd(HA));
+
+ *RA -= mpar * rho * cosd(gclat) * sind(HA) / cosd(*dec);
+ *dec -= mpar * rho * sind(gclat) * sind(g - *dec) / sind(g);
}
-//***********************************************************************/
-//* Name: calcGeomAnomalySun */
-//* Type: Function */
-//* Purpose: calculate the Geometric Mean Anomaly of the Sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* the Geometric Mean Anomaly of the Solar in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcGeomMeanAnomalySun(double t)
+double Moon::phase()
{
- return 357.52911 + t * (35999.05029 - 0.0001537 * t); // in degrees
+ return rev(m_eclipticLongitude - m_sun->lambda()) / 360.0;
}
-//***********************************************************************/
-//* Name: calcEccentricityEarthOrbit */
-//* Type: Function */
-//* Purpose: calculate the eccentricity of earth's orbit */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* the unitless eccentricity */
-//***********************************************************************/
-double NOAASolarCalc::calcEccentricityEarthOrbit(double t)
+void Moon::rotate(double* y, double* z)
{
- return 0.016708634 - t * (0.000042037 + 0.0000001267 * t); // unitless
+ double t = *y;
+ *y = t * cosd(m_obliquity) - *z * sind(m_obliquity);
+ *z = t * sind(m_obliquity) + *z * cosd(m_obliquity);
}
-//***********************************************************************/
-//* Name: calcSunEqOfCenter */
-//* Type: Function */
-//* Purpose: calculate the equation of center for the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunEqOfCenter(double t)
+void SolarSystemObject::calc()
{
- double m = calcGeomMeanAnomalySun(t);
+ double x, y, z;
+ double la, r;
- double mrad = degToRad(m);
- double sinm = sin(mrad);
- double sin2m = sin(mrad+mrad);
- double sin3m = sin(mrad+mrad+mrad);
+ L = rev(N + w + M);
+ double E0 = 720.0;
+ double E = M + (180.0 / M_PI) * e * sind(M) * (1.0 + e * cosd(M));
+ for (int j = 0; fabs(E0 - E) > 0.005 && j < 10; ++j) {
+ E0 = E;
+ E = E0 - (E0 - (180.0 / M_PI) * e * sind(E0) - M) / (1 - e * cosd(E0));
+ }
+ x = a * (cosd(E) - e);
+ y = a * sind(E) * sqrt(1.0 - e * e);
+ r = sqrt(x * x + y * y);
+ double v = rev(atan2d(y, x));
+ m_lambda = rev(v + w);
+ x = r * (cosd(N) * cosd(m_lambda) - sind(N) * sind(m_lambda) * cosd(i));
+ y = r * (sind(N) * cosd(m_lambda) + cosd(N) * sind(m_lambda) * cosd(i));
+ z = r * sind(m_lambda);
+ if (!qFuzzyCompare(i, 0.0)) {
+ z *= sind(i);
+ }
+ toSpherical(x, y, z, &m_eclipticLongitude, &la, &r);
+ if (calcPerturbations(&m_eclipticLongitude, &la, &r)) {
+ toRectangular(m_eclipticLongitude, la, r, &x, &y, &z);
+ }
+ rotate(&y, &z);
+ toSpherical(x, y, z, &RA, &dec, &rad);
+ topocentricCorrection(&RA, &dec);
- // in degrees
- return sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) +
- sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
+ HA = rev(siderealTime() - RA);
+ x = cosd(HA) * cosd(dec) * sind(m_latitude) - sind(dec) * cosd(m_latitude);
+ y = sind(HA) * cosd(dec);
+ z = cosd(HA) * cosd(dec) * cosd(m_latitude) + sind(dec) * sind(m_latitude);
+ m_azimuth = atan2d(y, x) + 180.0;
+ m_altitude = asind(z);
}
-//***********************************************************************/
-//* Name: calcSunTrueLong */
-//* Type: Function */
-//* Purpose: calculate the true longitude of the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar's true longitude in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunTrueLong(double t)
+double SolarSystemObject::siderealTime()
{
- double l0 = calcGeomMeanLongSun(t);
- double c = calcSunEqOfCenter(t);
-
- return l0 + c; // in degrees
+ double UT = m_utc.time().hour() + m_utc.time().minute() / 60.0 +
+ m_utc.time().second() / 3600.0;
+ double GMST0 = rev(282.9404 + 4.70935E-5 * m_day + 356.0470 + 0.9856002585 * m_day + 180.0);
+ return GMST0 + UT * 15.0 + m_longitude;
}
-//***********************************************************************/
-//* Name: calcSunTrueAnomaly */
-//* Type: Function */
-//* Purpose: calculate the true anamoly of the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar's true anamoly in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunTrueAnomaly(double t)
+void SolarSystemObject::calcForDateTime(const QDateTime& local, int offset)
{
- double m = calcGeomMeanAnomalySun(t);
- double c = calcSunEqOfCenter(t);
-
- return m + c; // in degrees
+ m_local = local;
+ m_utc = local.addSecs(-offset);
+ m_day = 367 * m_utc.date().year() - (7 * (m_utc.date().year() +
+ ((m_utc.date().month() + 9) / 12))) / 4 +
+ (275 * m_utc.date().month()) / 9 + m_utc.date().day() - 730530;
+ m_day += m_utc.time().hour() / 24.0 + m_utc.time().minute() / (24.0 * 60.0) +
+ m_utc.time().second() / (24.0 * 60.0 * 60.0);
+ m_obliquity = 23.4393 - 3.563E-7 * m_day;
}
-//***********************************************************************/
-//* Name: calcSunRadVector */
-//* Type: Function */
-//* Purpose: calculate the distance to the sun in AU */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar radius vector in AUs */
-//***********************************************************************/
-double NOAASolarCalc::calcSunRadVector(double t)
+SolarSystemObject::SolarSystemObject(double latitude, double longitude)
+ : m_latitude(latitude)
+ , m_longitude(-longitude)
{
- double v = calcSunTrueAnomaly(t);
- double e = calcEccentricityEarthOrbit(t);
-
- // in AUs
- return (1.000001018 * (1 - e * e)) / (1 + e * cos(degToRad(v)));
}
-//***********************************************************************/
-//* Name: calcSunApparentLong */
-//* Type: Function */
-//* Purpose: calculate the apparent longitude of the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar's apparent longitude in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunApparentLong(double t)
+SolarSystemObject::~SolarSystemObject()
{
- double o = calcSunTrueLong(t);
-
- double omega = 125.04 - 1934.136 * t;
- // in degrees
- return o - 0.00569 - 0.00478 * sin(degToRad(omega));
}
-//***********************************************************************/
-//* Name: calcMeanObliquityOfEcliptic */
-//* Type: Function */
-//* Purpose: calculate the mean obliquity of the ecliptic */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* mean obliquity in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcMeanObliquityOfEcliptic(double t)
+double SolarSystemObject::rev(double x)
{
- double seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * (0.001813)));
- return 23.0 + (26.0 + (seconds/60.0)) / 60.0; // in degrees
+ return x - floor(x / 360.0) * 360.0;
}
-//***********************************************************************/
-//* Name: calcObliquityCorrection */
-//* Type: Function */
-//* Purpose: calculate the corrected obliquity of the ecliptic */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* corrected obliquity in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcObliquityCorrection(double t)
+double SolarSystemObject::asind(double x)
{
- double e0 = calcMeanObliquityOfEcliptic(t);
-
- double omega = 125.04 - 1934.136 * t;
- return e0 + 0.00256 * cos(degToRad(omega)); // in degrees
+ return asin(x) * 180.0 / M_PI;
}
-//***********************************************************************/
-//* Name: calcSunRtAscension */
-//* Type: Function */
-//* Purpose: calculate the right ascension of the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar's right ascension in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunRtAscension(double t)
+double SolarSystemObject::sind(double x)
{
- double e = calcObliquityCorrection(t);
- double lambda = calcSunApparentLong(t);
-
- double tananum = (cos(degToRad(e)) * sin(degToRad(lambda)));
- double tanadenom = (cos(degToRad(lambda)));
- return radToDeg(atan2(tananum, tanadenom)); // in degrees
+ return sin(x * M_PI / 180.0);
}
-//***********************************************************************/
-//* Name: calcSunDeclination */
-//* Type: Function */
-//* Purpose: calculate the declination of the sun */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* solar's declination in degrees */
-//***********************************************************************/
-double NOAASolarCalc::calcSunDeclination(double t)
+double SolarSystemObject::cosd(double x)
{
- double e = calcObliquityCorrection(t);
- double lambda = calcSunApparentLong(t);
-
- double sint = sin(degToRad(e)) * sin(degToRad(lambda));
- return radToDeg(asin(sint)); // in degrees
+ return cos(x * M_PI / 180.0);
}
-//***********************************************************************/
-//* Name: calcHourAngle */
-//* Type: Function */
-//* Purpose: calculate the hour angle of the sun at sunset for the */
-//* latitude */
-//* Arguments: */
-//* zenith : zenith angle in degrees */
-//* lat : latitude of observer in degrees */
-//* solarDec : declination angle of sun in degrees */
-//* Return value: */
-//* hour angle of sunset in radians */
-//***********************************************************************/
-double NOAASolarCalc::calcHourAngle(double zenith, double solarDec, double latitude)
+double SolarSystemObject::tand(double x)
{
- double latRad = degToRad(latitude);
- double sdRad = degToRad(solarDec);
-
- //double HAarg = (cos(degToRad(90.833)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad));
-
- double HA = (acos(cos(degToRad(zenith)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad)));
-
- return HA; // in radians
+ return tan(x * M_PI / 180.0);
}
-//***********************************************************************/
-//* Name: calcEquationOfTime */
-//* Type: Function */
-//* Purpose: calculate the difference between true solar time and mean */
-//* solar time */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* Return value: */
-//* equation of time in minutes of time */
-//***********************************************************************/
-double NOAASolarCalc::calcEquationOfTime(double t)
+double SolarSystemObject::atan2d(double y, double x)
{
- double epsilon = calcObliquityCorrection(t);
- double l0 = calcGeomMeanLongSun(t);
- double e = calcEccentricityEarthOrbit(t);
- double m = calcGeomMeanAnomalySun(t);
+ return atan2(y, x) * 180.0 / M_PI;
+}
- double y = tan(degToRad(epsilon) / 2.0);
- y *= y;
-
- double sin2l0 = sin(2.0 * degToRad(l0));
- double sinm = sin(degToRad(m));
- double cos2l0 = cos(2.0 * degToRad(l0));
- double sin4l0 = sin(4.0 * degToRad(l0));
- double sin2m = sin(2.0 * degToRad(m));
-
- double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
- - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
-
- return radToDeg(Etime) * 4.0; // in minutes of time
+double SolarSystemObject::atand(double x)
+{
+ return atan(x) * 180.0 / M_PI;
}
-//***********************************************************************/
-//* Name: calcSolNoonUTC */
-//* Type: Function */
-//* Purpose: calculate the Universal Coordinated Time (UTC) of solar */
-//* noon for the given day at the given location on earth */
-//* Arguments: */
-//* t : number of Julian centuries since J2000.0 */
-//* longitude : longitude of observer in degrees */
-//* Return value: */
-//* time in minutes from zero Z */
-//***********************************************************************/
-double NOAASolarCalc::calcSolNoonUTC(double t, double longitude)
+void SolarSystemObject::toRectangular(double lo, double la, double r, double *x, double *y, double *z)
{
- // First pass uses approximate solar noon to calculate eqtime
- double tnoon = calcTimeJulianCent(calcJDFromJulianCent(t) + longitude / 360.0);
- double eqTime = calcEquationOfTime(tnoon);
- double solNoonUTC = 720 + (longitude * 4) - eqTime; // min
+ *x = r * cosd(lo) * cosd(la);
+ *y = r * sind(lo) * cosd(la);
+ *z = r * sind(la);
+}
- double newt = calcTimeJulianCent(calcJDFromJulianCent(t) -0.5 + solNoonUTC / 1440.0);
-
- eqTime = calcEquationOfTime(newt);
- // double solarNoonDec = calcSunDeclination(newt);
- solNoonUTC = 720 + (longitude * 4) - eqTime; // min
- return solNoonUTC;
+void SolarSystemObject::toSpherical(double x, double y, double z, double *lo, double *la, double *r)
+{
+ *r = sqrt(x * x + y * y + z * z);
+ *la = asind(z / *r);
+ *lo = rev(atan2d(y, x));
}
-//***********************************************************************/
-//* Name: calcTimeUTC */
-//* Type: Function */
-//* Purpose: calculate the Universal Coordinated Time (UTC) of sunset */
-//* for the given day at the given location on earth */
-//* Arguments: */
-//* zenith : zenith in degrees */
-//* rise : rise or set */
-//* JD : julian day */
-//* latitude : latitude of observer in degrees */
-//* longitude : longitude of observer in degrees */
-//* Return value: */
-//* time in minutes from zero Z */
-//***********************************************************************/
-void NOAASolarCalc::calcTimeUTC(double zenith, bool rise, double *jd, double *minutes,
- double latitude, double longitude)
+QPair<double, double> SolarSystemObject::zeroPoints(QPointF p1, QPointF p2, QPointF p3)
{
- forever {
- double t = calcTimeJulianCent(*jd);
- double m = (rise) ? 1.0 : -1.0;
+ double a = ((p2.y() - p1.y()) * (p1.x() - p3.x()) + (p3.y() - p1.y()) * (p2.x() - p1.x())) /
+ ((p1.x() - p3.x()) * (p2.x() * p2.x() - p1.x() * p1.x()) + (p2.x() - p1.x()) *
+ (p3.x() * p3.x() - p1.x() * p1.x()));
+ double b = ((p2.y() - p1.y()) - a * (p2.x() * p2.x() - p1.x() * p1.x())) / (p2.x() - p1.x());
+ double c = p1.y() - a * p1.x() * p1.x() - b * p1.x();
+ double discriminant = b * b - 4.0 * a * c;
+ double z1 = -1.0, z2 = -1.0;
- // *** Find the time of solar noon at the location, and use
- // that declination. This is better than start of the
- // Julian day
-
- double noonmin = calcSolNoonUTC(t, longitude);
- double tnoon = calcTimeJulianCent(*jd + noonmin / 1440.0);
-
- // First calculates sunrise and approx length of day
-
- double eqTime = calcEquationOfTime(tnoon);
- double solarDec = calcSunDeclination(tnoon);
- double hourAngle = m * calcHourAngle(zenith, solarDec, latitude);
-
- double delta = longitude - radToDeg(hourAngle);
- double timeDiff = 4 * delta;
- *minutes = 720 + timeDiff - eqTime;
-
- // first pass used to include fractional day in gamma calc
-
- double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + *minutes / 1440.0);
- eqTime = calcEquationOfTime(newt);
- solarDec = calcSunDeclination(newt);
- hourAngle = m * calcHourAngle(zenith, solarDec, latitude);
-
- delta = longitude - radToDeg(hourAngle);
- timeDiff = 4 * delta;
- *minutes = 720 + timeDiff - eqTime; // in minutes
- if (isnan(*minutes)) {
- *jd -= m;
- } else {
- return;
- }
+ if (discriminant >= 0.0) {
+ z1 = (-b + sqrt(discriminant)) / (2 * a);
+ z2 = (-b - sqrt(discriminant)) / (2 * a);
}
+ return QPair<double, double>(z1, z2);
}
-void NOAASolarCalc::calcAzimuthAndZenith(QDateTime dt, double eqTime, double zone,
- double solarDec, double latitude, double longitude,
- double *zenith, double *azimuth)
+QList< QPair<QDateTime, QDateTime> > SolarSystemObject::timesForAngles(const QList<double>& angles,
+ int offset)
{
- double solarTimeFix = eqTime - 4.0 * longitude + 60.0 * zone;
- double trueSolarTime = dt.time().hour() * 60.0 + dt.time().minute() +
- dt.time().second() / 60.0 + solarTimeFix;
- // in minutes
- while (trueSolarTime > 1440) {
- trueSolarTime -= 1440;
+ QList<double> altitudes;
+ QDate d = m_local.date();
+ QDateTime local(d, QTime(0, 0));
+ for (int j = 0; j <= 24; ++j) {
+ calcForDateTime(local, offset);
+ altitudes.append(altitude());
+ local = local.addSecs(60 * 60);
}
- //double hourAngle = calcHourAngle(timenow, m_longitude, eqTime);
- double hourAngle = trueSolarTime / 4.0 - 180.0;
- // Thanks to Louis Schwarzmayr for finding our error,
- // and providing the following 4 lines to fix it:
- if (hourAngle < -180.0) {
- hourAngle += 360.0;
- }
-
- double haRad = degToRad(hourAngle);
-
- double csz = sin(degToRad(latitude)) * sin(degToRad(solarDec)) +
- cos(degToRad(latitude)) * cos(degToRad(solarDec)) * cos(haRad);
- if (csz > 1.0) {
- csz = 1.0;
- } else if (csz < -1.0) {
- csz = -1.0;
- }
- *zenith = radToDeg(acos(csz));
- double azDenom = (cos(degToRad(latitude)) * sin(degToRad(*zenith)));
-
- if (fabs(azDenom) > 0.001) {
- double azRad = ((sin(degToRad(latitude)) * cos(degToRad(*zenith))) -
- sin(degToRad(solarDec))) / azDenom;
- if (fabs(azRad) > 1.0) {
- if (azRad < 0) {
- azRad = -1.0;
- } else {
- azRad = 1.0;
+ QList< QPair<QDateTime, QDateTime> > result;
+ QTime rise, set;
+ foreach (double angle, angles) {
+ for (int j = 3; j <= 24; j += 2) {
+ QPointF p1((j - 2) * 60 * 60, altitudes[j - 2] - angle);
+ QPointF p2((j - 1) * 60 * 60, altitudes[j - 1] - angle);
+ QPointF p3(j * 60 * 60, altitudes[j] - angle);
+ QPair<double, double> z = zeroPoints(p1, p2, p3);
+ if (z.first > p1.x() && z.first < p3.x()) {
+ if (p1.y() < 0.0) {
+ rise = QTime(0, 0).addSecs(z.first);
+ } else {
+ set = QTime(0, 0).addSecs(z.first);
+ }
}
+ if (z.second > p1.x() && z.second < p3.x()) {
+ if (p3.y() < 0.0) {
+ set = QTime(0, 0).addSecs(z.second);
+ } else {
+ rise = QTime(0, 0).addSecs(z.second);
+ }
+ }
}
-
- *azimuth = 180.0 - radToDeg(acos(azRad));
-
- if (hourAngle > 0.0) {
- *azimuth = -*azimuth;
- }
- } else {
- if (latitude > 0.0) {
- *azimuth = 180.0;
- } else {
- *azimuth = 0.0;
- }
+ result.append(QPair<QDateTime, QDateTime>(QDateTime(d, rise), QDateTime(d, set)));
}
- if (*azimuth < 0.0) {
- *azimuth += 360.0;
- }
+ return result;
}
-double NOAASolarCalc::calcElevation(double zenith)
+double SolarSystemObject::calcElevation()
{
- double exoatmElevation = 90.0 - zenith;
double refractionCorrection;
- if (exoatmElevation > 85.0) {
+
+ if (m_altitude > 85.0) {
refractionCorrection = 0.0;
} else {
- double te = tan(degToRad(exoatmElevation));
- if (exoatmElevation > 5.0) {
+ double te = tand(m_altitude);
+ if (m_altitude > 5.0) {
refractionCorrection = 58.1 / te - 0.07 / (te * te * te) +
0.000086 / (te * te * te * te * te);
- } else if (exoatmElevation > -0.575) {
- refractionCorrection = 1735.0 + exoatmElevation *
- (-518.2 + exoatmElevation * (103.4 + exoatmElevation *
- (-12.79 + exoatmElevation * 0.711) ) );
+ } else if (m_altitude > -0.575) {
+ refractionCorrection = 1735.0 + m_altitude *
+ (-518.2 + m_altitude * (103.4 + m_altitude *
+ (-12.79 + m_altitude * 0.711) ) );
} else {
refractionCorrection = -20.774 / te;
}
refractionCorrection = refractionCorrection / 3600.0;
}
- double solarZen = zenith - refractionCorrection;
- return 90.0 - solarZen;
+ return m_altitude + refractionCorrection;
}
-
Index: solarposition.h
===================================================================
--- solarposition.h (revision 959786)
+++ solarposition.h (working copy)
@@ -26,7 +26,7 @@
SolarPosition();
virtual ~SolarPosition();
- void appendData(Plasma::DataEngine::Data &data);
+ void appendData(Plasma::DataEngine::Data &data, bool sun, bool moon);
private:
QHash<QString, Plasma::DataEngine::Data> m_cache;
signature.asc
Description: This is a digitally signed message part.
_______________________________________________ Plasma-devel mailing list Plasma-devel@kde.org https://mail.kde.org/mailman/listinfo/plasma-devel
