Skip to content

Instantly share code, notes, and snippets.

@ksasao

ksasao/sunset.h

Last active March 10, 2025 00:40
Show Gist options
  • Save ksasao/af9e6a274ed44da678a11efd13e70b95 to your computer and use it in GitHub Desktop.
Save ksasao/af9e6a274ed44da678a11efd13e70b95 to your computer and use it in GitHub Desktop.
Download ZIP
現在地の日没時刻/日の出時刻を表示します https://x.com/ksasao/status/1898305387969012007
Raw
sunset.h
#ifndef SUNSET_H_
#define SUNSET_H_
#include<math.h>
// 計算式は下記のWebサイトを参考に実装しています
// 日出日没計算、やってみよう
// https://hhsprings.pinoko.jp/site-hhs/2015/02/%e6%97%a5%e5%87%ba%e6%97%a5%e6%b2%a1%e8%a8%88%e7%ae%97%e3%80%81%e3%82%84%e3%81%a3%e3%81%a6%e3%81%bf%e3%82%88%e3%81%86/
typedef struct _DateTime {
int year;
int month;
int day;
int hour;
int minute;
int second;
int msec;
} DateTime;
double Sind(double d)
{
return sin(d * M_PI / 180);
}
// cos function using degree
double Cosd(double d)
{
return cos(d * M_PI / 180);
}
// tan function using degree
double tand(double d)
{
return tan(d * M_PI / 180);
}
// calculate Julius year (year from 2000/1/1, for variable "t")
double jy(int yy, int mm, int dd, int h, double m, int s, int i)
{ // yy/mm/dd h:m:s, i: time difference
yy -= 2000;
if (mm <= 2)
{
mm += 12;
yy--;
}
double k = 365.0 * yy + 30.0 * mm + dd - 33.5 - i / 24.0 + floor(3.0 * (mm + 1) / 5.0)
+ floor(yy / 4.0) - floor(yy / 100.0) + floor(yy / 400.0);
k += ((s / 60.0 + m) / 60.0 + h) / 24.0; // plus time
k += (65 + yy) / 86400; // plus delta T
return k / 365.25;
}
// solar position1 (celestial longitude, degree)
double spls(double t)
{ // t: Julius year
double l = 280.4603 + 360.00769 * t
+ (1.9146 - 0.00005 * t) * Sind(357.538 + 359.991 * t)
+ 0.0200 * Sind(355.05 + 719.981 * t)
+ 0.0048 * Sind(234.95 + 19.341 * t)
+ 0.0020 * Sind(247.1 + 329.640 * t)
+ 0.0018 * Sind(297.8 + 4452.67 * t)
+ 0.0018 * Sind(251.3 + 0.20 * t)
+ 0.0015 * Sind(343.2 + 450.37 * t)
+ 0.0013 * Sind(81.4 + 225.18 * t)
+ 0.0008 * Sind(132.5 + 659.29 * t)
+ 0.0007 * Sind(153.3 + 90.38 * t)
+ 0.0007 * Sind(206.8 + 30.35 * t)
+ 0.0006 * Sind(29.8 + 337.18 * t)
+ 0.0005 * Sind(207.4 + 1.50 * t)
+ 0.0005 * Sind(291.2 + 22.81 * t)
+ 0.0004 * Sind(234.9 + 315.56 * t)
+ 0.0004 * Sind(157.3 + 299.30 * t)
+ 0.0004 * Sind(21.1 + 720.02 * t)
+ 0.0003 * Sind(352.5 + 1079.97 * t)
+ 0.0003 * Sind(329.7 + 44.43 * t);
while (l >= 360) { l -= 360; }
while (l < 0) { l += 360; }
return l;
}
// solar position2 (distance, AU)
double spds(double t)
{ // t: Julius year
double r = (0.007256 - 0.0000002 * t) * Sind(267.54 + 359.991 * t)
+ 0.000091 * Sind(265.1 + 719.98 * t)
+ 0.000030 * Sind(90.0)
+ 0.000013 * Sind(27.8 + 4452.67 * t)
+ 0.000007 * Sind(254 + 450.4 * t)
+ 0.000007 * Sind(156 + 329.6 * t);
r = pow(10, r);
return r;
}
// solar position3 (declination, degree)
double spal(double t)
{ // t: Julius year
double ls = spls(t);
double ep = 23.439291 - 0.000130042 * t;
double al = atan(tand(ls) * Cosd(ep)) * 180 / M_PI;
if ((ls >= 0) && (ls < 180))
{
while (al < 0) { al += 180; }
while (al >= 180) { al -= 180; }
}
else
{
while (al < 180) { al += 180; }
while (al >= 360) { al -= 180; }
}
return al;
}
// solar position4 (the right ascension, degree)
double spdl(double t)
{ // t: Julius year
double ls = spls(t);
double ep = 23.439291 - 0.000130042 * t;
double dl = asin(Sind(ls) * Sind(ep)) * 180 / M_PI;
return dl;
}
// Calculate sidereal hour (degree)
double sh(double t, int h, double m, int s, double l, int i)
{ // t: julius year, h: hour, m: minute, s: second,
// l: longitude, i: time difference
double d = ((s / 60.0 + m) / 60.0 + h) / 24.0; // elapsed hour (from 0:00 a.m.)
double th = 100.4606 + 360.007700536 * t + 0.00000003879 * t * t - 15 * i;
th += l + 360.0 * d;
while (th >= 360) { th -= 360; }
while (th < 0) { th += 360; }
return th;
}
// Calculating the seeming horizon altitude "sa"(degree)
double eandp(double alt, double ds)
{ // subfunction for altitude and parallax
double e = 0.035333333 * sqrt(alt);
double p = 0.002442818 / ds;
return p - e;
}
double sa(double alt, double ds)
{ // alt: altitude (m), ds: solar distance (AU)
double s = 0.266994444 / ds;
double r = 0.585555555;
double k = eandp(alt, ds) - s - r;
return k;
}
// Calculating solar alititude (degree) {
double soal(double la, double th, double al, double dl)
{ // la: latitude, th: sidereal hour,
// al: solar declination, dl: right ascension
double h = Sind(dl) * Sind(la) + Cosd(dl) * Cosd(la) * Cosd(th - al);
h = asin(h) * 180 / M_PI;
return h;
}
// Calculating solar direction (degree) {
double sodr(double la, double th, double al, double dl)
{ // la: latitude, th: sidereal hour,
// al: solar declination, dl: right ascension
double t = th - al;
double dc = -Cosd(dl) * Sind(t);
double dm = Sind(dl) * Sind(la) - Cosd(dl) * Cosd(la) * Cosd(t);
double dr = 0;
if (dm == 0)
{
double st = Sind(t);
if (st > 0) dr = -90;
if (st == 0) dr = 9999;
if (st < 0) dr = 90;
}
else
{
dr = atan(dc / dm) * 180 / M_PI;
if (dm < 0) dr += 180;
}
if (dr < 0) dr += 360;
return dr;
}
double FineTune(int yy, int mm, int dd, int hh, double m, int s, int i, double la, double lo, double alt)
{
m = m - 1.0;
double t = jy(yy, mm, dd, hh, m, 0, i);
double th = sh(t, hh, m, 0, lo, i);
double ds = spds(t);
double ls = spls(t);
double alp = spal(t);
double dlt = spdl(t);
double pht = soal(la, th, alp, dlt);
for (double mmm = m; mmm < m+1; mmm += 0.00197)
{
t = jy(yy, mm, dd, hh, mmm, 0, i);
th = sh(t, hh, mmm, 0, lo, i);
ds = spds(t);
ls = spls(t);
alp = spal(t);
dlt = spdl(t);
double ht = soal(la, th, alp, dlt);
double dr = sodr(la, th, alp, dlt);
double t4 = sa(alt, ds);
if( ((pht - t4) * (ht - t4)) < 0)
{
return mmm;
}
pht = ht;
}
return m;
}
DateTime MinToSec(double m,DateTime d)
{
int min = (int)floor(m);
int sec = (int)floor((m - min) * 60.0);
int milli = (int)floor((((m - min) * 60.0)-sec)*1000.0);
d.minute = min;
d.second = sec;
d.msec = milli;
return d;
}
DateTime Sunset(int yy, int mm, int dd, int i, double la, double lo, double alt, bool isSunrise = false)
{ // main routine
DateTime result;
double t = jy(yy, mm, dd - 1, 23, 59, 0, i);
double th = sh(t, 23, 59, 0, lo, i);
double ds = spds(t);
double ls = spls(t);
double alp = spal(t);
double dlt = spdl(t);
double pht = soal(la, th, alp, dlt);
if(alt<0){
alt = 0;
}
double t4 = sa(alt, ds);
for (int j = 0; j < 24 * 60; j++) {
// 最初は粗い粒度（1分単位）で探索
int hh = j / 60;
double m = (double)(j % 60);
t = jy(yy, mm, dd, hh, m, 0, i);
th = sh(t, hh, m, 0, lo, i);
ds = spds(t);
ls = spls(t);
alp = spal(t);
dlt = spdl(t);
double ht = soal(la, th, alp, dlt);
if ( ((pht > t4) && (ht < t4) && !isSunrise) || ((pht < t4) && (ht > t4) && isSunrise)) {
// 条件が満たされる粗い領域を見つけたら、より細かい粒度で再探索
j--;
double start = j;
double end = j + 60.0;
for (double fi = start; fi <= end; fi += 0.1) {
hh = (int)(fi / 60.0);
double fm = (double)(fi - 60.0 * hh);
t = jy(yy, mm, dd, hh, fm, 0, i);
th = sh(t, hh, fm, 0, lo, i);
ds = spds(t);
ls = spls(t);
alp = spal(t);
dlt = spdl(t);
double fht = soal(la, th, alp, dlt);
if ( ((pht > t4) && (fht < t4) && !isSunrise) || ((pht < t4) && (fht > t4) && isSunrise)) {
double m2 = FineTune(yy, mm, dd, hh, fm, 0, i, la, lo, alt);
result.year = yy;
result.month = mm;
result.day = dd;
result.hour = hh;
result = MinToSec(m2, result);
return result;
}
pht = fht;
}
}
pht = ht;
}
return result;
}
#endif
Raw
sunsettime.ino
// 現在地の日没時刻/日の出時刻を表示します
// M5StickC Plus2 / M5Stack GPS Unit v1.1 (AT6668)
// Board: M5Stack 2.1.3 / Lib: M5StickCPlus2 1.0.1
#include <M5GFX.h>
#include <M5StickCPlus2.h>
#include <TinyGPS++.h>
#include <Preferences.h>
#include "sunset.h"
#define TZ_OFFSET (9)
bool isSunrise = false;
//static const uint32_t GPSBaud = 9600; // for M5Stack GPS Unit
static const uint32_t GPSBaud = 115200; // for M5Stack GPS Unit v1.1 (AT6668)
TinyGPSPlus gps;
M5GFX display;
char txt[64];
long oldTime = -1;
void drawSunset(int fix){
StickCP2.Display.setTextColor(TFT_WHITE, TFT_BLACK);
StickCP2.Display.setTextWrap(false);
StickCP2.Display.setTextDatum(textdatum_t::top_right);
StickCP2.Display.setTextSize(0.75,1);
int year = gps.date.year();
int month = gps.date.month();
int day = gps.date.day();
int hour = gps.time.hour();
int minute = gps.time.minute();
int second = gps.time.second();
double lat = gps.location.lat();
double lng = gps.location.lng();
double alt = gps.altitude.meters();
if(year < 2022 || lng < -180 || lng > 180 || hour < 0 || minute < 0 || second < 0 || hour > 23 || minute > 59 || second > 59){
return;
}
DateTime sunset = Sunset(year,month,day, TZ_OFFSET, lat, lng, alt, isSunrise);
// 時刻描画
StickCP2.Display.setTextSize(0.5,0.75);
StickCP2.Display.setFont(&fonts::Font8);
StickCP2.Display.setTextDatum(textdatum_t::top_left);
sprintf(txt,"%02d:%02d:%02d", sunset.hour,sunset.minute,sunset.second);
StickCP2.Display.drawString(txt, 0,40);
M5.Log.printf(txt);
// 時刻の小数点以下を小さく出力
StickCP2.Display.setTextSize(0.5/2.0,0.75/2.0);
StickCP2.Display.setFont(&fonts::Font8);
sprintf(txt,".%03d\n", sunset.msec);
StickCP2.Display.drawString(txt, 190,68);
M5.Log.printf(txt);
// 単位描画
StickCP2.Display.setTextSize(1,1);
StickCP2.Display.setFont(&fonts::lgfxJapanGothic_32);
StickCP2.Display.setTextDatum(textdatum_t::top_left);
if(isSunrise){
StickCP2.Display.drawString("日の出時刻　", 0,0);
}else{
StickCP2.Display.drawString("日没時刻　", 0,0);
}
// 衛星数描画
sprintf(txt,"Sat.: %d ", fix);
StickCP2.Display.drawString(txt, 5,100);
}
void setup(){
auto cfg = M5.config();
StickCP2.begin(cfg);
Serial.begin(115200);
M5.delay(500);
StickCP2.Display.init();
StickCP2.Display.setRotation(1);
Serial2.begin(GPSBaud, SERIAL_8N1, 33, 32);
Serial2.flush();
updateStatus();
M5.Log.println("Started.");
}
void updateStatus(){
int fix = gps.satellites.value();
M5.Log.print("Location: ");
M5.Log.printf("%f,%f,%f\n",gps.location.lat(),gps.location.lng(),gps.altitude.meters());
if (!StickCP2.Display.displayBusy())
{
StickCP2.Display.startWrite();
drawSunset(fix);
StickCP2.Display.endWrite();
}
}
void loop(){
StickCP2.update(); // for Btn input
long now = gps.time.value();
if(now != oldTime){
updateStatus();
oldTime = now;
}
if (millis() > 5000 && gps.charsProcessed() < 10){
M5.Log.println("No GPS data received: check wiring");
}
if(StickCP2.BtnA.wasPressed()){
isSunrise = !isSunrise;
StickCP2.Display.clear(0);
}
smartDelay(50);
}
char* buf[256];
// This custom version of delay() ensures that the gps object
// is being "fed".
static void smartDelay(unsigned long ms)
{
unsigned long start = millis();
do
{
while (Serial2.available()){
int c = Serial2.read();
gps.encode(c);
}
while (Serial.available()){
int c = Serial.read();
Serial2.write((char)c);
}
} while (millis() - start < ms);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment