Calculating Sunrise and Sunset 2

[flotulopex]

Hi all,

I'm looking for an algorithm calculating the sunrise/sunset time. The few threads/posts found here didn't help, unfortunately.

My wall clock using a DS3231 RTC is updated once daily by the means of a GPS module. So I do have all coordinates et precise time information.

Does anyone have done this yet in PICBasic?

[Jerson]

I have done this in C many years back. However, you need to use double precision floating point numbers to get accurate results. PicBasic code? I doubt it as Picbasic does not support this data type.

[flotulopex]

Yes, this is what I found out too.

Is there any chip that would do this math for me I could embed in my circuit, at least as a plan B?

[HenrikOlsson]

I don't know how to calculate it but do you really have to? I mean, do you have to calculate it on the PIC, at runtime? It's not like your wall clock is changing location on a day to day basis, is it?

You could calculate it (or simply look it up) on your PC, store the results in your PBP program and look it up at runtime. If calculating it is "heavy" then a pre-calculated table is probably smaller in size and definetly faster.

[flotulopex]

Henrik,

Right, no need for runtime data.

Any idea where a kind of data table would be to be found?

[Ioannis]

I think it is best to be calculated in a Excel file and then the results stored in an External EEPROM?

Ioannis

[pedja089]

You can store lot of data in program memory. And get it with read code.

[HenrikOlsson]

Henrik,

Right, no need for runtime data.

Any idea where a kind of data table would be to be found?

Let me Google that for you....
https://www.timeanddate.com/sun/swed...h=10&year=2024

It defaulted to Gothenburg for me so obviously you'll need to change that.

[amgen]

Isn't the changing sunrise and sunset times mostly a sine wave function, with the 2 times opposite of each other. Basic has a sin function. and the sin values would synchronize with the date and the sin multiplier for minutes change/day value would depend on the latitude of the location of the clock....

[Acetronics2]

I do know it's " another language " ( and written in French ! ) ... but might be useful

volatile int f_wdt=1;
volatile int Nb_Inter;

const int Fuseau_horaire=0; // Pour la france +0
//const float Latitude = 46.309; // latiture du poulailler pour le calcul du lever et du coucher du soleil
//const float Longitude = 5.165; // longitude idem ""
const float Latitude = 49.916549; // latiture du poulailler pour le calcul du lever et du coucher du soleil
const float Longitude = 1.661536; // longitude idem ""
const float TensionMaxA0= 22.90; // Tension Max mesurable pour 1023 points dépend du pont de résistance
const float SeuilAlerteA0 = 10; // Alerte batterie faible sur mesure A0
const float TensionMaxA1= 5; // Tension Max mesurable pour 1023 points dépend du pont de résistance
const float SeuilAlerteA1 = -1; // Alerte batterie faible sur mesue A1
const int Nb_Boucles_Mini = 0; // Nb d'interruption Mini 1 = 8s
const int Nb_Boucles_Maxi = 0; // Nb d'interruption Maxi 1 = 8s valeur par defaut 8
const unsigned long TempsCycleStd = 1000; //Temps de cyle standard en ms valeur 60000
const unsigned long TempsCycleAnomalie = 1000; //Temps de cyle anomalie en ms
const int on_off_moteur=2; // Sortie N°2 pour on_off moteur
const int sens_moteur=3; // Sortie N°3 pour le sens
const int Porte_ON=4; // Entrée N° 4 pour le contact de porte fermée
const int Porte_OFF=5; // Entrée N° 5 pour le contat de porte ouverte
const int M_A=6; // Entrée N° 6 pour la fonction marche /arrêt
const int Led_Nuit=7; // Sortie N°7 pour Flash Nuit
const int DC=8; // Sortie N°8 pour activer convertisseur DC-DC
const int Beeper_Out=12; // Sortie N°12 pour Beeper
const int Seuil_Nuit = 500; // Seuil de détection de la nuit à calibrer en fonction de la cellule photorésitance
const int Seuil_Jour = 700; // Seuil de détection du jour à calibrer en fonction de la cellule photorésitance
const int Nb_Anomalie_Max=1; // Nb d'essai en cas d'anomalie
const unsigned long TempsMaxConst= 35; // Limite de temps ouveture ou fermeture avant anomalie en secondes
const unsigned long TempsSecurite = 10; //  % temps complémentaire se sécurité après temps étalonnage
const int Marge=0.50; // latitude en heures pour tenir compte de la photorésistance par rapport aux dates prévues de fermeture et d'ouverture
const float Aube=0.50; // temps en heures entre l'aube et le levé du soleil
const float Crepuscule=1.00; // temps en heures entre le couché du soleil et le crépuscule
const int Nb_Cycles_Detection=5; // Nb de cycles pour détection jour-nuit cellule 
const boolean Debug_2=false; // flag pour le débug avancé
const boolean Calcul_Soleil_OK = true; // flag pour utiliser le calcul du soleil sinon utilisation du tabelau des heures
const int Heure_chgt_horaire = 3; // à 3h le dernier dimache de mars, il sera 4h (+1) idem avec -1 le dernier dimanche d'octobre
const boolean Etalonnage_OK = true; // détermine si le système engage un cycle d'étalonnage pour mesurer le temps d'ouverture
const int Bauds =9600;// détermine la vitesse de communication de la RS
const int Tps_secu =10;// Temps de sécurité de maintien de la descente lors de la fermeture de la porte en ms
const boolean Debugage = false; // à mettre à true pour faire des modification et activer le mode débug

//Variables

boolean Etalonnage_Fait = false;
unsigned long TempsMax;
int Heure=21;
int Minute=30;
int Jour=24;
int Mois=01;
int Annee=2015;
boolean Date; // Mode réception Date via RS
int IndDate=0; // Indice de réception format Date via RS sous Forme HHMMJJMMAAAA
int Nb_Anomalie=0; // Nb d'anomalie constatées 
boolean Debug=false; // active l'envoi par ligne série des infos de débug
unsigned long prevtime;  // Variable qui stock le temps
unsigned long MemoTemps; // Variable qui stock en temps  en secondes
unsigned long TempsDerMvt; // Variable qui stock le temps en secondes du dernier mvt
float Temps; // Variable qui stock le temps en hh.mm
float nuit_mois[13]; // Tableau Heure de fermeture hh.mm
float jour_mois[13]; // Tableau Heure d'ouverture hh.mm
unsigned long TempsCycle; // Temps de cycle
float Heure_Ouverture;
float Heure_Fermeture;
float Heure_Couche_Soleil;
float Heure_Leve_Soleil;
boolean RTC_OK=true; // utilise ou non le module RTC Tiny à mettre en constante si taille mémoire trop grande


int Led_Arduino = 13; // N° Led interne Arduino
boolean Memo_Led_Nuit = false; // Pilote declignotement de la led Nuit
boolean Memo_Led_Arduino = false; // Pilote declignotement de la led Nuit
int Etat_Porte=-1; // Etat de la porte 0=Ouvertr 1=fermée -1 en mouvement (indéterminé)
boolean Anomalie = false; // Gère les anomalies reset par M_A
boolean Beeper = true; // active la gestion des alertes beep des niveau bas batterie
char Carlu; // Caratère lu
boolean Interruption = true; // 
boolean Autoris_Inter = true; // autorise le mode interrution en marche
int Ana0; // Pour lecture Analogique 0
int Ana1; // Pour lecture Analogique 1
int Ana2; // Pour lecture Analogique 2
int Nb_Boucles; // Compte le nb d'interruption
int Nb_Cycles_Jour; // Compte le nb de cycles avec détection Jour
int Nb_Cycles_Nuit; // Compte le nb de cycles avec détection Nuit
int Rang_Jour; // Rang d'une journée dans l'année pour calculer lever et coucher du soleil
int Memo_Rang_Jour; // Mémorise le rang jour
int Rang_Jour_date_chgt_hr_ete; // Mémorise le rang jour pour la date de chgt d'heure d'été
int Rang_Jour_date_chgt_hr_hiver; // Mémorise le rang jour pour la date de chgt d'heure d'été
int Heure_ete_OK = true; 
int Memo_Annee =0;
boolean Chgt_hr_ete_fait=false;
boolean Chgt_hr_hiver_fait=false;


// convertit un angle degré en radian
float conv_rad(float angle_degre)
{
  return( (angle_degre * 71.0) / 4068.0);
}

// convertit un angle radian en degré
float conv_deg(float angle_radian)
{
  return ( (angle_radian * 4068.0) / 71.0);
}

// calcul de N° du jour dans la semaine 0=Dimanche, 1= Lundi etc..
// Jour de semaine D = { [(23m)/9] + d + 4 + y + [z/4] - [z/100] + [z/400] - 2 (si m >= 3) } mod 7 
// algo de Mike Keith
int calcul_jour_semaine(int jour,int mois,int annee)
{
  int num_jour_semaine;
  int z;
  if (mois >=3 )
  {
    z= annee;
    num_jour_semaine = ((23*mois)/9+ jour+ 4+ annee+(z/4)-(z/100)+(z/400)-2) %7;
  }
  else
  {
    z= annee-1;
    num_jour_semaine = ((23*mois)/9+ jour+ 4+ annee+(z/4)-(z/100)+(z/400)) %7;
  }
  // num_jour_semaine = (jour + annee + annee/4 - annee/100 + annee/400 + ((31*mois)/12) ) % 7; ancienne formule

  if (Debug_2)
  {
    Serial.print("num jour semaine:");
    Serial.print(num_jour_semaine);
    Serial.println(); 
  }
  return(num_jour_semaine); 
}

// calcul de dernier dimanche du mois
void calcul_dates_chgt_hr (int annee)
{
  boolean trouver_dimanche = false;
  // on recherche le dernier dimanche de mars
  int num_jour=31; // pour le mois de mars
  while (!trouver_dimanche)
  {
    if (calcul_jour_semaine(num_jour,3,annee)==0)
      trouver_dimanche = true;
    else
      num_jour--;
  }
  Rang_Jour_date_chgt_hr_ete =calcul_rang_date(num_jour,3,annee);

  // on recherche le dernier dimanche d'octobre
  trouver_dimanche = false;
  num_jour=31; // pour le mois d'octobre
  while (!trouver_dimanche)
  {
    if (calcul_jour_semaine(num_jour,10,annee)==0)
      trouver_dimanche = true;
    else
      num_jour--;
  }
  Rang_Jour_date_chgt_hr_hiver =calcul_rang_date(num_jour,10,annee);
  if (Debug)
  {
    Serial.print("Rang Jour :");
    Serial.print(Rang_Jour);
    Serial.print(" Rang chtg ete:");
    Serial.print(Rang_Jour_date_chgt_hr_ete);
    Serial.print(" Rang chtg hiver:");
    Serial.print(Rang_Jour_date_chgt_hr_hiver);
    Serial.println(); 
  }
}



// calcul le rang d'un journée pour le calcul du coucher du soleil http://jean-paul.cornec.pagesperso-orange.fr/formule_rang.htm
// vérifier les calculs sur tables http://jean-paul.cornec.pagesperso-orange.fr/rang_des_jours.htm
int calcul_rang_date(int jour,int mois,int annee) 
{
  int N1;
  int N2;
  int K;
  int Rang;

  N1= int(( mois* 275 ) / 9);
  N2 = int( (mois + 9) / 12 );
  K = 1 + int( ( annee - 4 * int (annee / 4) + 2 ) / 3 );
  Rang = N1 - N2*K + jour -30;
  if (Debug_2)
  {
    Serial.print("Rang :");
    Serial.print(Rang);
    Serial.println(); 
  }
  return(Rang);
}

// calcul l'heure de lever et du coucher du soleil
// basé sur le calcul simplifié de jean-paul Cornec http://jean-paul.cornec.pagesperso-orange.fr/heures_lc.htm

void calcul_soleil()
{

  float M;
  float C;
  float L;
  float R;
  float Et;
  float Dec;
  float Ho;
  float VL;
  float TL;
  float HL;
  float VC;
  float TC;
  float HC;

  M = 357.0 + (0.9856 * Rang_Jour);
  C = (1.914 * sin(conv_rad(M))) + (0.02 * sin(conv_rad(2.0*M)));
  L = 280.0 + C + (0.9856 * Rang_Jour);

  R = (-2.465 * sin(conv_rad(2.0*L))) + (0.053 * sin(conv_rad(4.0*L)));
  Et = 4.0*(C + R);
  Dec =conv_deg(asin(0.3978 * sin(conv_rad(L))));
  Ho = conv_deg(acos(( -0.01454 - sin ( conv_rad(Dec) ) * sin ( conv_rad(Latitude) ) ) / ( cos ( conv_rad(Dec) ) * cos ( conv_rad(Latitude) ) )));
  VL = 12 - Ho/15;
  TL = 12 - Ho/15 + Et/60 + Longitude/15;
  if (Heure_ete_OK)
    HL = 12 - Ho/15 + Et/60 + Longitude/15 + Fuseau_horaire +1; //(+ 1 pour heure été)
  else 
    HL = 12 - Ho/15 + Et/60 + Longitude/15 + Fuseau_horaire;

  Heure_Leve_Soleil=HL;

  VC = 12 + Ho/15;
  TC = 12 + Ho/15 + Et/60 + Longitude/15;
  if (Heure_ete_OK)
    HC = 12 + Ho/15 + Et/60 + Longitude/15 + Fuseau_horaire +1; //(+ 1 pour heure été)
  else
    HC = 12 + Ho/15 + Et/60 + Longitude/15 + Fuseau_horaire;

  Heure_Couche_Soleil=HC;
  if (Debug_2)
  {
    Serial.print("M:");
    Serial.println(M);
    Serial.print("C:");
    Serial.println(C);
    Serial.print("L:");
    Serial.println(L);
    Serial.print("R:");
    Serial.println(R);
    Serial.print("Et:");
    Serial.println(Et);
    Serial.print("Dec:");
    Serial.println(Dec);
    Serial.print("Ho:");
    Serial.println(Ho);
    Serial.print("VL:");
    Serial.println(VL);
    Serial.print("HL:");
    Serial.println(HL);
    Serial.print("VC:");
    Serial.println(VC);
    Serial.print("TC:");
    Serial.println(TC);
    Serial.print("HC:");
    Serial.println(HC);
  }
}

Alain

[HenrikOlsson]

Like I said, for a fixed location I'd use a lookup table instead of expensive and complicated computations at runtime.
For example, in the southern part of Sweden (Malmö to be precise) the sunrise occurs between 04:23 and 08:37 (a window of 4h14min or 254 minutes). Sunset occurs between 15:35 and 21:55 (a window of 6h20min or 380min).

If we can accept 2 minute resolution then we can store sunrise and sunset as "2 minute ticks from 04:20 for sunrise and 2 minute ticks from 15:30 for sunset using 8bits for each day using only 16bits so 730bytes worth of flash memory. I'd be surprised if you can get better resolution and/or smaller code doing it at runtime.

The drawback, obviously, is that you need to adapt and encode the data for your specific location and changing it is not as simple as punching in GPS coordinates. But for a one-off wall clock it's certainly means to an end IMHO.

[amgen]

but also....... with a few known data points like Hennrik states, a simple calculation is possible. The change in rise and set times per day is a known from charts for a given latitude. So at some given latitude, the overall change in rise and set in hours and minutes can be estimated. And that happens for 180 days at a time, longer for 180 days and then going shorter for 180 days. So with Hennric's example, sunrise changes for 254 minutes ....or 254/180.... on average adds 1.4 minutes a day. Similarly, sunset (with it's minutes calculation) goes X minutes lower each day. Could keep a byte counter to track days of 1 to 182 and when 182 is reached, the additions and subtractions are switched. Those days are Summer Solstice, about June 21 ..... and Winter Solstice ABOUT December 21. Those are the dates to switch the additions and subtractions for rise and set... So only check for dates of change (182 days apart), keep track of mode (season... getting longer/shorter) and add or subtract minutes daily from some basic set point.

[flotulopex]

Nice idea Henrik 👍

[flotulopex]

Just a thought: from what I understand from Amgen's and Henrik's suggestions, is that the sunrise and sunset would be roughtly linearly calculable.

Is it so?

[Jerson]

Depending on your application, the approach you take to solve your problem could differ.

For a changeable lattitude, longitude, Height above sea level, DST the calculation is not trivial. The astronomical calculations used require double precision (15digits floating point) math to get accurate results.

If your location is fixed, then perhaps a table driven approach with corrections applied to the table for different dates of the year may work - I am not sure of this.

You may find C code for the calculations if you search for such calculations.

IMO : this is most certainly not easy to undetake in PicBasic.

[Acetronics2]

table example here ... https://www.suntoday.org/sunrise-sunset/2024.html

Alain

[amgen]

Just a thought: from what I understand from Amgen's and Henrik's suggestions, is that the sunrise and sunset would be roughtly linearly calculable

black is my attempt at a straight line.... follows actual times fairly close enough to actual exact times. It is just adding and subtracting X minutes each day for 1/2 a year then switching the add-subtract.

[rocket_troy]

And there's always NBit math libraries if you really wanted to do with with a PIC. I've used them a lot with PICs and PBP. The workarounds can be painful if you're dealing with >32bit numbers that could either be positive or negative as the numbers aren't signed.

Troy