replace the SunPosition library for the PSA code

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 <math.h>
+
+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;
+    }
+}

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

src/main.cpp

@@ -9,7 +9,8 @@
 #include <SoftwareSerial.h>
 
 #include <wmm.h>
-#include <SolarPosition.h>
+// #include <SolarPosition.h>
+#include <sunpos.h>
 
 #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();