From 7f2c635454cfb6ed563697944f627cbe7a77b2ba Mon Sep 17 00:00:00 2001 From: Luis Rodil-Fernandez Date: Mon, 1 Aug 2022 12:44:59 +0200 Subject: [PATCH] replace the SunPosition library for the PSA code --- lib/sunpos/sunpos.cpp | 103 ++++++++++++++++++++++++++++++++++++++++++ lib/sunpos/sunpos.h | 38 ++++++++++++++++ src/main.cpp | 38 ++++++++++++---- 3 files changed, 171 insertions(+), 8 deletions(-) create mode 100644 lib/sunpos/sunpos.cpp create mode 100644 lib/sunpos/sunpos.h diff --git a/lib/sunpos/sunpos.cpp b/lib/sunpos/sunpos.cpp new file mode 100644 index 0000000..4f73a0d --- /dev/null +++ b/lib/sunpos/sunpos.cpp @@ -0,0 +1,103 @@ +// This file is available in electronic form at http://www.psa.es/sdg/sunpos.htm + +#include "sunpos.h" +#include + +void sunpos(cTime udtTime,cLocation udtLocation, cSunCoordinates *udtSunCoordinates) +{ + // Main variables + double dElapsedJulianDays; + double dDecimalHours; + double dEclipticLongitude; + double dEclipticObliquity; + double dRightAscension; + double dDeclination; + + // Auxiliary variables + double dY; + double dX; + + // Calculate difference in days between the current Julian Day + // and JD 2451545.0, which is noon 1 January 2000 Universal Time + { + double dJulianDate; + long int liAux1; + long int liAux2; + // Calculate time of the day in UT decimal hours + dDecimalHours = udtTime.dHours + (udtTime.dMinutes + + udtTime.dSeconds / 60.0 ) / 60.0; + // Calculate current Julian Day + liAux1 =(udtTime.iMonth-14)/12; + liAux2=(1461*(udtTime.iYear + 4800 + liAux1))/4 + (367*(udtTime.iMonth + - 2-12*liAux1))/12- (3*((udtTime.iYear + 4900 + + liAux1)/100))/4+udtTime.iDay-32075; + dJulianDate=(double)(liAux2)-0.5+dDecimalHours/24.0; + // Calculate difference between current Julian Day and JD 2451545.0 + dElapsedJulianDays = dJulianDate-2451545.0; + } + + // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the + // ecliptic in radians but without limiting the angle to be less than 2*Pi + // (i.e., the result may be greater than 2*Pi) + { + double dMeanLongitude; + double dMeanAnomaly; + double dOmega; + dOmega=2.1429-0.0010394594*dElapsedJulianDays; + dMeanLongitude = 4.8950630+ 0.017202791698*dElapsedJulianDays; // Radians + dMeanAnomaly = 6.2400600+ 0.0172019699*dElapsedJulianDays; + dEclipticLongitude = dMeanLongitude + 0.03341607*sin( dMeanAnomaly ) + + 0.00034894*sin( 2*dMeanAnomaly )-0.0001134 + -0.0000203*sin(dOmega); + dEclipticObliquity = 0.4090928 - 6.2140e-9*dElapsedJulianDays + +0.0000396*cos(dOmega); + } + + // Calculate celestial coordinates ( right ascension and declination ) in radians + // but without limiting the angle to be less than 2*Pi (i.e., the result may be + // greater than 2*Pi) + { + double dSin_EclipticLongitude; + dSin_EclipticLongitude= sin( dEclipticLongitude ); + dY = cos( dEclipticObliquity ) * dSin_EclipticLongitude; + dX = cos( dEclipticLongitude ); + dRightAscension = atan2( dY,dX ); + if( dRightAscension < 0.0 ) dRightAscension = dRightAscension + twopi; + dDeclination = asin( sin( dEclipticObliquity )*dSin_EclipticLongitude ); + } + + // Calculate local coordinates ( azimuth and zenith angle ) in degrees + { + double dGreenwichMeanSiderealTime; + double dLocalMeanSiderealTime; + double dLatitudeInRadians; + double dHourAngle; + double dCos_Latitude; + double dSin_Latitude; + double dCos_HourAngle; + double dParallax; + dGreenwichMeanSiderealTime = 6.6974243242 + + 0.0657098283*dElapsedJulianDays + + dDecimalHours; + dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime*15 + + udtLocation.dLongitude)*rad; + dHourAngle = dLocalMeanSiderealTime - dRightAscension; + dLatitudeInRadians = udtLocation.dLatitude*rad; + dCos_Latitude = cos( dLatitudeInRadians ); + dSin_Latitude = sin( dLatitudeInRadians ); + dCos_HourAngle= cos( dHourAngle ); + udtSunCoordinates->dZenithAngle = (acos( dCos_Latitude*dCos_HourAngle + *cos(dDeclination) + sin( dDeclination )*dSin_Latitude)); + dY = -sin( dHourAngle ); + dX = tan( dDeclination )*dCos_Latitude - dSin_Latitude*dCos_HourAngle; + udtSunCoordinates->dAzimuth = atan2( dY, dX ); + if ( udtSunCoordinates->dAzimuth < 0.0 ) + udtSunCoordinates->dAzimuth = udtSunCoordinates->dAzimuth + twopi; + udtSunCoordinates->dAzimuth = udtSunCoordinates->dAzimuth/rad; + // Parallax Correction + dParallax=(dEarthMeanRadius/dAstronomicalUnit) + *sin(udtSunCoordinates->dZenithAngle); + udtSunCoordinates->dZenithAngle=(udtSunCoordinates->dZenithAngle + + dParallax)/rad; + } +} \ No newline at end of file diff --git a/lib/sunpos/sunpos.h b/lib/sunpos/sunpos.h new file mode 100644 index 0000000..f099e60 --- /dev/null +++ b/lib/sunpos/sunpos.h @@ -0,0 +1,38 @@ +// This file is available in electronic form at http://www.psa.es/sdg/sunpos.htm + +#ifndef __SUNPOS_H +#define __SUNPOS_H + + +// Declaration of some constants +#define pi 3.14159265358979323846 +#define twopi (2*pi) +#define rad (pi/180) +#define dEarthMeanRadius 6371.01 // In km +#define dAstronomicalUnit 149597890 // In km + +struct cTime +{ + int iYear; + int iMonth; + int iDay; + double dHours; + double dMinutes; + double dSeconds; +}; + +struct cLocation +{ + double dLongitude; + double dLatitude; +}; + +struct cSunCoordinates +{ + double dZenithAngle; + double dAzimuth; +}; + +void sunpos(cTime udtTime, cLocation udtLocation, cSunCoordinates *udtSunCoordinates); + +#endif \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index da64e8f..500c9fb 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,7 +9,8 @@ #include #include -#include +// #include +#include #define SCREEN_WIDTH 128 #define SCREEN_HEIGHT 64 @@ -76,10 +77,10 @@ void draw_sunpos() { display.setCursor(0, 0); //display.print(F("lat ------> ")); display.print(F("lat ....... ")); - display.print((double)sun.latitude, 2); + display.print((double)sun.latitude, 4); display.println(); display.print(F("lon ....... ")); - display.print((double)sun.longitude, 2); + display.print((double)sun.longitude, 4); display.println(); display.print(F("WMM ....... ")); display.print((double)sun.declination, 2); @@ -140,18 +141,39 @@ void loop() { sun.latitude = GPS.latitudeDegrees; sun.longitude = GPS.longitudeDegrees; + // method 2 - using SolarPosition library // create a fixed UNIX time to test fixed time method - tmElements_t someTime = {GPS.seconds, GPS.minute, GPS.hour, 0, GPS.day, GPS.month, CalendarYrToTm(GPS.year) }; - sun.epoch = makeTime(someTime); + // tmElements_t someTime = {GPS.seconds, GPS.minute, GPS.hour, 0, GPS.day, GPS.month, CalendarYrToTm(GPS.year) }; + // sun.epoch = makeTime(someTime); if (GPS.fix) { - SolarPosition location(GPS.latitudeDegrees, GPS.longitudeDegrees); + // SolarPosition location(GPS.latitudeDegrees, GPS.longitudeDegrees); // calculate magnetic declination E0000(GPS.latitudeDegrees, GPS.longitudeDegrees, wmmtime, &sun.declination); - sun.elevation = location.getSolarPosition(sun.epoch).elevation; - sun.azimuth = location.getSolarPosition(sun.epoch).azimuth; + cLocation loc; + loc.dLatitude = GPS.latitudeDegrees; + loc.dLongitude = GPS.longitudeDegrees; + + cTime t; + t.dHours = GPS.hour; + t.dMinutes = GPS.minute; + t.dSeconds = GPS.seconds; + t.iDay = GPS.day; + t.iMonth = GPS.month; + t.iYear = GPS.year; + + cSunCoordinates sunloc; + + sunpos(t, loc, &sunloc); + + sun.azimuth = sunloc.dAzimuth; + sun.elevation = 90 - sunloc.dZenithAngle; + + // method 2 - using SolarPosition library + // sun.elevation = location.getSolarPosition(sun.epoch).elevation; + // sun.azimuth = location.getSolarPosition(sun.epoch).azimuth; debug_solar(); draw_sunpos();