Browse Source

first commit

master
Luis Rodil-Fernandez 2 years ago
commit
6fd89c7ee1
  1. 5
      README.md
  2. 46
      lib/README
  3. 3
      lib/wmm/README.md
  4. 302
      lib/wmm/wmm.cpp
  5. 44
      lib/wmm/wmm.h
  6. 26
      platformio.ini
  7. 162
      src/main.cpp

5
README.md

@ -0,0 +1,5 @@
## Heliograph prototype
ESP32 with a GPS module calculates sun position and magnetic declination.
- Is this enough for an earth bound object to know where to find the sun?

46
lib/README

@ -0,0 +1,46 @@
This directory is intended for project specific (private) libraries.
PlatformIO will compile them to static libraries and link into executable file.
The source code of each library should be placed in a an own separate directory
("lib/your_library_name/[here are source files]").
For example, see a structure of the following two libraries `Foo` and `Bar`:
|--lib
| |
| |--Bar
| | |--docs
| | |--examples
| | |--src
| | |- Bar.c
| | |- Bar.h
| | |- library.json (optional, custom build options, etc) https://docs.platformio.org/page/librarymanager/config.html
| |
| |--Foo
| | |- Foo.c
| | |- Foo.h
| |
| |- README --> THIS FILE
|
|- platformio.ini
|--src
|- main.c
and a contents of `src/main.c`:
```
#include <Foo.h>
#include <Bar.h>
int main (void)
{
...
}
```
PlatformIO Library Dependency Finder will find automatically dependent
libraries scanning project source files.
More information about PlatformIO Library Dependency Finder
- https://docs.platformio.org/page/librarymanager/ldf.html

3
lib/wmm/README.md

@ -0,0 +1,3 @@
World Magnetic Model library slightly modified
from original written by Rene Herrero Gomez
https://github.com/sevenseas-io/wmm

302
lib/wmm/wmm.cpp

@ -0,0 +1,302 @@
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include "wmm.h"
#define PI_CONST 3.14159265359f
#define RADIANS_TO_DEGREES 0.017453292f
#define DEGREES_TO_RADIANS (PI_CONST / 180.0f)
#define A_CONST 6378.137f
#define A2_CONST (A_CONST * A_CONST)
#define B_CONST 6356.7523142f
#define B2_CONST (B_CONST * B_CONST)
#define RE_CONST 6371.2f
#define A4_CONST (A2_CONST * A2_CONST)
#define B4_CONST (B2_CONST * B2_CONST)
#define C2_CONST (A2_CONST - B2_CONST)
#define C4_CONST (A4_CONST - B4_CONST)
#define COEFFICIENTS_COUNT 90U
static float c[13][13];
static float cd[13][13];
static float k[13][13];
static float snorm[169];
static float fn[13];
static float fm[13];
const uint8_t wmm_cof_entries_encoded[] =
{0xDD, 0xF2, 0x23, 0x00, 0x83, 0x01, 0x00, 0xEB, 0xE2, 0x01, 0x81, 0xD7, 0x05, 0x8D, 0x01, 0xFB, 0x03, 0xE8, 0x86, 0x03, 0x00, 0xF3, 0x01, 0x00, 0xBC, 0xD1, 0x03, 0xDC, 0xD3, 0x03, 0xC7, 0x01, 0xEE, 0x04, 0x80, 0x86, 0x02, 0xF4, 0x72, 0x56, 0xEF, 0x03, 0x87, 0xD5, 0x01, 0x00, 0x1C, 0x00, 0xC2, 0xF4, 0x02, 0xF6, 0x0C, 0x7E, 0x39, 0x8A, 0xC1, 0x01, 0xB2, 0x25, 0x22, 0x4A, 0x89, 0x52, 0xF5, 0x54, 0xFA, 0x01, 0x0B, 0x87, 0x8D, 0x01, 0x00, 0x4B, 0x00, 0x9E, 0x7E, 0x84, 0x2C, 0x50, 0x02, 0x9E, 0x0D, 0xF0, 0x18, 0x7C, 0x85, 0x01, 0xD6, 0x30, 0x8E, 0x1F, 0x36, 0x25, 0x9F, 0x07, 0xED, 0x36, 0x77, 0x78, 0xE8, 0x24, 0x00, 0x43, 0x00, 0xAF, 0x38, 0x9D, 0x07, 0x06, 0x01, 0x96, 0x1D, 0xA4, 0x20, 0x47, 0x19, 0xFF, 0x15, 0xFD, 0x12, 0x01, 0x49, 0xE8, 0x17, 0x82, 0x05, 0x0C, 0x1E, 0x89, 0x02, 0x9F, 0x0F, 0x0A, 0x05, 0x93, 0x0A, 0x00, 0x46, 0x00, 0x90, 0x0A, 0xFF, 0x02, 0x44, 0x01, 0x9A, 0x0B, 0xBA, 0x03, 0x05, 0x52, 0xFF, 0x12, 0x8F, 0x08, 0x0E, 0x4E, 0xEA, 0x05, 0xC4, 0x0A, 0x4E, 0x09, 0x87, 0x02, 0x9A, 0x01, 0x00, 0x01, 0xC7, 0x0A, 0xA9, 0x0A, 0x08, 0x0A, 0xA6, 0x0C, 0x00, 0x41, 0x00, 0xC0, 0x0C, 0xC2, 0x08, 0x43, 0x05, 0xD3, 0x01, 0xE8, 0x02, 0x41, 0x06, 0xB5, 0x08, 0x17, 0x07, 0x47, 0x9E, 0x02, 0xAB, 0x03, 0x02, 0x42, 0x80, 0x01, 0x56, 0x45, 0x4C, 0xC8, 0x01, 0xD0, 0x04, 0x48, 0x02, 0xA2, 0x01, 0x53, 0x0A, 0x03, 0xAC, 0x03, 0x00, 0x41, 0x00, 0xA2, 0x01, 0x94, 0x01, 0x01, 0x43, 0xEF, 0x02, 0xD9, 0x02, 0x41, 0x07, 0x44, 0x80, 0x02, 0x05, 0x42, 0xD3, 0x03, 0xF6, 0x01, 0x41, 0x05, 0x99, 0x02, 0x95, 0x02, 0x04, 0x43, 0x89, 0x02, 0x24, 0x05, 0x45, 0xE5, 0x02, 0xC5, 0x01, 0x00, 0x04, 0x43, 0x1C, 0x04, 0x01, 0x32, 0x00, 0x41, 0x00, 0x92, 0x01, 0xE9, 0x03, 0x42, 0x43, 0x1D, 0xAF, 0x01, 0x00, 0x02, 0x4E, 0xA2, 0x01, 0x04, 0x44, 0x4B, 0x73, 0x43, 0x04, 0xC5, 0x02, 0x7E, 0x00, 0x01, 0x0B, 0x8E, 0x01, 0x03, 0x00, 0x99, 0x01, 0x04, 0x00, 0x42, 0xDD, 0x01, 0x4F, 0x00, 0x05, 0xF7, 0x01, 0xA1, 0x01, 0x44, 0x02, 0x53, 0x00, 0x00, 0x00, 0x7E, 0x22, 0x00, 0x00, 0x41, 0x42, 0x00, 0x01, 0x11, 0x23, 0x02, 0x43, 0x49, 0x30, 0x41, 0x01, 0x06, 0xD6, 0x01, 0x42, 0x42, 0x49, 0x41, 0x00, 0x01, 0x13, 0x6A, 0x41, 0x00, 0x0E, 0x62, 0x42, 0x41, 0x58, 0x41, 0x41, 0x02, 0x67, 0xD8, 0x01, 0x00, 0x00, 0x1E, 0x00, 0x00, 0x00, 0x4E, 0x00, 0x41, 0x00, 0x59, 0x1A, 0x00, 0x01, 0x18, 0x45, 0x00, 0x00, 0x49, 0x44, 0x00, 0x02, 0x03, 0x06, 0x41, 0x00, 0x47, 0x42, 0x00, 0x00, 0x41, 0x51, 0x00, 0x01, 0x0E, 0x50, 0x41, 0x00, 0x46, 0x5E, 0x41, 0x41, 0x02, 0x54, 0x41, 0x00, 0x1F, 0x5A, 0x41, 0x00, 0x54, 0x00, 0x00, 0x00, 0x41, 0x4C, 0x00, 0x00, 0x05, 0x05, 0x00, 0x00, 0x0D, 0x0D, 0x00, 0x41, 0x4C, 0x52, 0x00, 0x01, 0x07, 0x01, 0x00, 0x00, 0x03, 0x07, 0x00, 0x00, 0x05, 0x41, 0x00, 0x00, 0x42, 0x06, 0x00, 0x01, 0x45, 0x02, 0x00, 0x00, 0x01, 0x49, 0x00, 0x00, 0x4B, 0x00, 0x00, 0x00, 0x43, 0x05, 0x41, 0x41};
static wmm_cof_record_t wmm_cof_entries[COEFFICIENTS_COUNT];
static float convert_varint_to_float(char **bytes);
float wmm_get_date(uint8_t year, uint8_t month, uint8_t date)
{
return (float)year + 2000.0f + (float)(month - 1U) / 12.0f + (float)(date - 1U) / (365.0f);
}
static float convert_varint_to_float(char **bytes)
{
float result;
int32_t result_int;
bool negative = false;
bool first_byte = true;
uint8_t shift;
do
{
if (first_byte)
{
if (**bytes & 0x40)
{
negative = true;
}
result_int = **bytes & 0x3f;
shift = 6U;
first_byte = false;
}
else
{
result_int += (uint32_t)(**bytes & 0x7f) << shift;
shift += 7U;
}
if ((**bytes & 0x80) == 0U)
{
(*bytes)++;
break;
}
(*bytes)++;
} while (true);
result = ((float)result_int) / 10.0f;
if (negative)
{
result = -result;
}
return result;
}
void wmm_init(void)
{
uint8_t j;
uint8_t m;
uint8_t n;
uint8_t D2;
float gnm;
float hnm;
float dgnm;
float dhnm;
float flnmj;
uint8_t i;
char *bytes = (char *)&wmm_cof_entries_encoded[0];
// unpack coefficients
for (i = 0U; i < COEFFICIENTS_COUNT; i++)
{
wmm_cof_entries[i].gnm = convert_varint_to_float(&bytes);
wmm_cof_entries[i].hnm = convert_varint_to_float(&bytes);
wmm_cof_entries[i].dgnm = convert_varint_to_float(&bytes);
wmm_cof_entries[i].dhnm = convert_varint_to_float(&bytes);
}
c[0][0] = 0.0f;
cd[0][0] = 0.0f;
j = 0U;
for (n = 1U; n <= 12U; n++)
{
for (m = 0U; m <= n; m++)
{
gnm = wmm_cof_entries[j].gnm;
hnm = wmm_cof_entries[j].hnm;
dgnm = wmm_cof_entries[j].dgnm;
dhnm = wmm_cof_entries[j].dhnm;
j++;
if (m <= n)
{
c[m][n] = gnm;
cd[m][n] = dgnm;
if (m != 0U)
{
c[n][m - 1U] = hnm;
cd[n][m - 1U] = dhnm;
}
}
}
}
// CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED
*snorm = 1.0f;
for (n = 1U; n <= 12U; n++)
{
*(snorm + n) = *(snorm + n - 1U) * (float)(2U * n - 1U) / (float)n;
j = 2U;
m = 0U;
for (D2 = n - m + 1U; D2 > 0U; D2--)
{
k[m][n] = (float)(((n - 1U) * (n - 1U)) - (m * m)) / (float)((2U * n - 1U) * (2U * n - 3U));
if (m > 0U)
{
flnmj = (float)((n - m + 1U) * j) / (float)(n + m);
*(snorm + n + m * 13U) = *(snorm + n + (m - 1U) * 13U) * sqrt(flnmj);
j = 1U;
c[n][m - 1U] = *(snorm + n + m * 13U) * c[n][m - 1U];
cd[n][m - 1U] = *(snorm + n + m * 13U) * cd[n][m - 1U];
}
c[m][n] = *(snorm + n + m * 13U) * c[m][n];
cd[m][n] = *(snorm + n + m *13U) * cd[m][n];
m += 1U;
}
fn[n] = (float)(n + 1U);
fm[n] = (float)n;
}
k[1][1] = 0.0f;
}
void E0000(float glat, float glon, float time_years, float *dec)
{
static float tc[13][13];
static float sp[13];
static float cp[13];
static float dp[13][13];
static float pp[13];
float dt = time_years - WMM_EPOCH;
float rlon = glon * DEGREES_TO_RADIANS;
float rlat = glat * DEGREES_TO_RADIANS;
float srlon = sinf(rlon);
float srlat = sinf(rlat);
float crlon = cosf(rlon);
float crlat = cosf(rlat);
float srlat2 = srlat * srlat;
float crlat2 = crlat * crlat;
sp[0] = 0.0f;
sp[1] = srlon;
cp[0] = 1.0f;
cp[1] = crlon;
dp[0][0] = 0.0f;
pp[0] = 1.0f;
// CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS
float q = sqrtf(A2_CONST - C2_CONST * srlat2);
float q2 = (A2_CONST / (B2_CONST)) * (A2_CONST / B2_CONST);
float ct = srlat / sqrtf(q2 * crlat2 + srlat2);
float st = sqrtf(1.0f - (ct * ct));
float r2 = (A4_CONST - C4_CONST * srlat2) / (q * q);
float r = sqrtf(r2);
float d = sqrtf(A2_CONST * crlat2 + B2_CONST * srlat2);
float ca = d / r;
float sa = C2_CONST * crlat * srlat / (r * d);
for (uint8_t m = 2U; m <= 12U; m++)
{
sp[m] = sp[1] * cp[m - 1U] + cp[1] * sp[m - 1U];
cp[m] = cp[1] * cp[m - 1U] - sp[1] * sp[m - 1U];
}
float aor = RE_CONST / r;
float ar = aor * aor;
float br = 0.0f;
float bt = 0.0f;
float bp = 0.0f;
float bpp = 0.0f;
for (uint16_t n = 1U; n <= 12U; n++)
{
ar = ar * aor;
uint8_t m = 0U;
for (uint8_t D4 = n + 1U; D4 > 0U; D4--)
{
// COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS AND DERIVATIVES VIA RECURSION RELATIONS
if (n == m)
{
*(snorm + n + m * 13U) = st * *(snorm + n - 1U + (m - 1U) * 13U);
dp[m][n] = st * dp[m - 1U][n - 1U] + ct * *(snorm + n - 1U + (m - 1U) * 13U);
goto S50;
}
if (n == 1U && m == 0U)
{
*(snorm + n + m * 13U) = ct * *(snorm + n - 1U + m * 13U);
dp[m][n] = ct * dp[m][n - 1U] - st * *(snorm + n - 1U + m * 13U);
goto S50;
}
if (n > 1U && n != m)
{
if (m > n - 2U)
{
*(snorm + n - 2U + m * 13U) = 0.0f;
}
if (m > n - 2U)
{
dp[m][n - 2U] = 0.0f;
}
*(snorm + n + m * 13U) = ct * *(snorm + n - 1U + m * 13U) - k[m][n] * *(snorm + n - 2U + m * 13U);
dp[m][n] = ct * dp[m][n - 1U] - st * *(snorm + n - 1U + m * 13U) - k[m][n] * dp[m][n - 2U];
}
S50:
// TIME ADJUST THE GAUSS COEFFICIENTS
tc[m][n] = c[m][n] + dt * cd[m][n];
if (m != 0U)
{
tc[n][m - 1U] = c[n][m - 1U] + dt * cd[n][m - 1U];
}
// ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
float par = ar * *(snorm + n + m * 13U);
float temp1;
float temp2;
if (m == 0)
{
temp1 = tc[m][n] * cp[m];
temp2 = tc[m][n] * sp[m];
}
else
{
temp1 = tc[m][n] * cp[m] + tc[n][m - 1U] * sp[m];
temp2 = tc[m][n] * sp[m] - tc[n][m - 1U] * cp[m];
}
bt = bt - ar * temp1 * dp[m][n];
bp += (fm[m] * temp2 * par);
br += (fn[n] * temp1 * par);
// SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES
if (st == 0.0f && m == 1U)
{
if (n == 1U)
{
pp[n] = pp[n - 1U];
}
else
{
pp[n] = ct * pp[n - 1U] - k[m][n] * pp[n - 2U];
}
bpp += (fm[m] * temp2 * ar * pp[n]);
}
m += 1U;
}
}
if (st == 0.0f)
{
bp = bpp;
}
else
{
bp /= st;
}
// ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO GEODETIC COORDINATES
float bx = -bt * ca - br * sa;
float by = bp;
// COMPUTE DECLINATION
*dec = atan2f(by, bx) / DEGREES_TO_RADIANS;
}

44
lib/wmm/wmm.h

@ -0,0 +1,44 @@
#ifndef WMM_H
#define WMM_H
#include <stdint.h>
#define WMM_EPOCH 2020.0f
typedef struct
{
float gnm;
float hnm;
float dgnm;
float dhnm;
} wmm_cof_record_t;
/**
* Initialize the WMM. Needs calling only once.
*/
void wmm_init(void);
/**
* Get the date in WMM format
*
* @param year Year in 2 digit format of 21st centuary, i.e. 20 represents 2020
* @param month Month, 1 to 12
* @param date Date of month, 1 to 31
* @return Date in WMM format
* @note No checking of illegal dates is done
*/
float wmm_get_date(uint8_t year, uint8_t month, uint8_t date);
/**
* Get the magnetic variation at a point on the earth's surface
*
* @param glat Latitude in degrees and fractional degrees, negative west
* @param glon Longitude in degrees and fractional degrees, negative west
* @param time_years The date as returned from wmm_get_date
* @param dec Pointer to float holding calculated magnetic variation (also known as declination). Negative is west.
* @note The altitude used is the ellipsoid at the supplied latitude/longitude, not the earth's surface. This will
* give very small errors in some parts of the world comapred to sea level.
*/
void E0000(float glat, float glon, float time_years, float *dec);
#endif

26
platformio.ini

@ -0,0 +1,26 @@
; PlatformIO Project Configuration File
;
; Build options: build flags, source filter
; Upload options: custom upload port, speed and extra flags
; Library options: dependencies, extra library storages
; Advanced options: extra scripting
;
; Please visit documentation for the other options and examples
; https://docs.platformio.org/page/projectconf.html
[env:lolin32]
platform = espressif32
board = lolin32
framework = arduino
monitor_speed = 115200
lib_deps =
SSD1306=adafruit/Adafruit SSD1306@^2.5.6
Adafruit_GPS=https://github.com/adafruit/Adafruit_GPS/archive/refs/tags/1.7.0.zip
SoftwareSerial=https://github.com/plerup/espsoftwareserial/archive/refs/tags/6.16.1.zip
Time=paulstoffregen/Time@^1.6.1
SolarPosition=https://github.com/KenWillmott/SolarPosition.git
; NeoGPS=slashdevin/NeoGPS@^4.2.9

162
src/main.cpp

@ -0,0 +1,162 @@
#include <Arduino.h>
#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>
#include <Adafruit_GPS.h>
#include <SoftwareSerial.h>
#include <wmm.h>
#include <SolarPosition.h>
#define SCREEN_WIDTH 128
#define SCREEN_HEIGHT 64
#define OLED_RESET -1
#define SCREEN_ADDRESS 0x3C
Adafruit_SSD1306 display(SCREEN_WIDTH, SCREEN_HEIGHT, &Wire, OLED_RESET);
SoftwareSerial altserial(14, 12);
Adafruit_GPS GPS(&altserial);
uint32_t tstamp = millis();
typedef struct {
time_t epoch;
float latitude;
float longitude;
float declination;
float elevation;
float azimuth;
} tSolar;
tSolar sun;
void init_display() {
Wire.begin(5, 4);
if(!display.begin(SSD1306_SWITCHCAPVCC, SCREEN_ADDRESS, false, false)) {
Serial.println(F("SSD1306 allocation failed"));
while(true) ; // do not proceed
}
display.clearDisplay();
}
void init_gps() {
altserial.begin(9600);
delay(2000);
// Serial.println("Software Serial GPS Test Echo Test");
// you can send various commands to get it started
//mySerial.println(PMTK_SET_NMEA_OUTPUT_RMCGGA);
altserial.println(PMTK_SET_NMEA_OUTPUT_ALLDATA);
altserial.println(PMTK_SET_NMEA_UPDATE_1HZ);
// Serial.println("Get version!");
altserial.println(PMTK_Q_RELEASE);
}
// ///////////////////////////////////////////////////////////////////////////////////////
// ///////////////////////////////////////////////////////////////////////////////////////
// ///////////////////////////////////////////////////////////////////////////////////////
void setup() {
Serial.begin(115200);
init_display();
init_gps();
wmm_init();
}
void draw_sunpos() {
display.clearDisplay();
display.setTextSize(1);
display.setTextColor(SSD1306_WHITE);
display.setCursor(0, 0);
//display.print(F("lat ------> "));
display.print(F("lat ....... "));
display.print((double)sun.latitude, 2);
display.println();
display.print(F("lon ....... "));
display.print((double)sun.longitude, 2);
display.println();
display.print(F("WMM ....... "));
display.print((double)sun.declination, 2);
display.println();
display.print(F("elevation .. "));
display.print((double)sun.elevation, 4);
display.println();
display.print(F("azimuth .... "));
display.print((double)sun.azimuth, 4);
display.println();
// refresh
display.display();
}
void debug_solar() {
Serial.print("Location: ");
Serial.print(sun.latitude, 4); //Serial.print(GPS.lat);
Serial.print(", ");
Serial.print(sun.longitude, 4); //Serial.println(GPS.lon);
Serial.println();
Serial.print("WMM declination "); Serial.print(sun.declination, 4);
Serial.println();
Serial.print(F("elevation: "));
Serial.print(sun.elevation, 4);
Serial.print(F(" azimuth: "));
Serial.println(sun.azimuth, 4);
Serial.println();
// Serial.print("Speed (knots): "); Serial.println(GPS.speed);
// Serial.print("Angle: "); Serial.println(GPS.angle);
// Serial.print("Altitude: "); Serial.println(GPS.altitude);
// Serial.print("Satellites: "); Serial.println((int)GPS.satellites);
// Serial.print("Antenna status: "); Serial.println((int)GPS.antenna);
}
void loop() {
char c = GPS.read();
// // if you want to debug, this is a good time to do it!
//if (c) Serial.write(c);
// if a sentence is received, we can check the checksum, parse it...
if (GPS.newNMEAreceived()) {
if (!GPS.parse(GPS.lastNMEA())) // this also sets the newNMEAreceived() flag to false
return; // we can fail to parse a sentence in which case we should just wait for another
}
// approximately every 2 seconds or so, print out the current stats
if (millis() - tstamp > 2000) {
tstamp = millis(); // reset the timer
float wmmtime = wmm_get_date(GPS.year, GPS.month, GPS.day);
sun.latitude = GPS.latitudeDegrees;
sun.longitude = GPS.longitudeDegrees;
// 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);
if (GPS.fix) {
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;
debug_solar();
draw_sunpos();
}
}
// draw_sunpos();
} // loop()
Loading…
Cancel
Save