Go Down

Topic: Orologio con Alba/Tramonto da calcoli astronomici (Read 7226 times) previous topic - next topic

Federico Vanzati

Dec 21, 2010, 02:41 pm Last Edit: Dec 21, 2010, 02:57 pm by Federico_Vanzati Reason: 1
in questi giorni mi sto cimentando a creare (con arduino :)) un orologio che sia in grado di calcolare gli orari di alba e tramonto, così da far accendere e spegnere delle luci in automatico.

Il codice per adesso è in grado di far dialogare arduino con un RTC esterno (il classico DS1307) e tramite una funzione che contiene qualche formula di conversione di calcolare gli orari di alba e tramonto, non essendo un esperto di astronomia, non sono certo dell'esattezza del procedimento ma, sta di fatto che con delle simulazioni i risultati ottenuti sono piuttosto soddisfacenti (manca però il calcolo dei crepuscoli).

Prima di passare alle fasi di introduzione delle condizioni che gestiranno un relè, sto cercando di verificare l'affidabilità del chip RTC e del codice che ho scritto fino a questo punto.

Ed ecco il codice:

prima parte:
Quote

#include <Math.h>
#include <WProgram.h>
#include <Wire.h>
#include <DS1307.h> // written by  mattt on the Arduino forum and modified by D. Sjunnesson


int command = 0;
int prevDay = 0;
int actDay = 0;

float deg_rad = 57.2958; // conversione gradi radianti
float omega = 0.9863; //(= 360/365) [gradi/giorno] -> gradi per ogni giro / giorni_anno
double mezzo_di_medio = 12.1667; //corrisponde alle 12h10m -> 50m indietro rispetto a Greenwich
float lat = 45.47;
double alba, tramonto;

int alba_tramonto_HM[4];

void setup()
{
  
  Serial.begin(57600);
  omega = omega/deg_rad; //few conversion

  for(int i=0; i<5; i++)
    alba_tramonto_HM = 0;
}

void loop()
{
  if(Serial.available()){
    command = Serial.read();
    if (command == 84) {      //If command = "T" Set Date
      setRTC();
    }
  }

  actDay = RTC.get(DS1307_DATE,false);
  if( actDay != prevDay ){
    calcAlba_Tramonto( mezzo_di_medio );
    prevDay = actDay;
  }

  Serial.print(RTC.get(DS1307_HR,true)); //read the hour and also update all the values by pushing in true
  Serial.print(":");
  Serial.print(RTC.get(DS1307_MIN,false));//read minutes without update (false)
  Serial.print(":");
  Serial.print(RTC.get(DS1307_SEC,false));//read seconds
  Serial.print("    ");                 // some space for a more happy life
  Serial.print(RTC.get(DS1307_DATE,false));//read date
  Serial.print("/");
  Serial.print(RTC.get(DS1307_MTH,false));//read month
  Serial.print("/");
  Serial.print(RTC.get(DS1307_YR,false)); //read year
  Serial.println();

  Serial.println(calcDOY(RTC.get(DS1307_DATE,false), RTC.get(DS1307_MTH,false), RTC.get(DS1307_YR,false)));

  Serial.print("alba: ");
  Serial.print(alba_tramonto_HM[0], DEC);
  Serial.print("h ");
  Serial.print(alba_tramonto_HM[1], DEC);
  Serial.println("m ");

  Serial.print("tramonto: ");
  Serial.print( alba_tramonto_HM[2], DEC);
  Serial.print("h ");
  Serial.print( alba_tramonto_HM[3], DEC);
  Serial.println("m ");

  double ora_dec = conv_sessagimali( RTC.get(DS1307_HR,true), RTC.get(DS1307_MIN,true) );


  delay(1000);

}

void setRTC()
{
  byte second = (byte) ((Serial.read() - 48) * 10 + (Serial.read() - 48)); // Use of (byte) type casting and ascii math to achieve result.  
  byte minute = (byte) ((Serial.read() - 48) *10 +  (Serial.read() - 48));
  byte hour  = (byte) ((Serial.read() - 48) *10 +  (Serial.read() - 48));
  byte dayOfWeek = (byte) (Serial.read() - 48);
  byte dayOfMonth = (byte) ((Serial.read() - 48) *10 +  (Serial.read() - 48));
  byte month = (byte) ((Serial.read() - 48) *10 +  (Serial.read() - 48));
  byte year= (byte) ((Serial.read() - 48) *10 +  (Serial.read() - 48));

  RTC.stop();
  RTC.set(DS1307_SEC,second);          //set the seconds
  RTC.set(DS1307_MIN,minute);          //set the minutes
  RTC.set(DS1307_HR,hour);             //set the hours
  RTC.set(DS1307_DOW,dayOfWeek);       //set the day of the week
  RTC.set(DS1307_DATE,dayOfMonth);     //set the date
  RTC.set(DS1307_MTH,month);           //set the month
  RTC.set(DS1307_YR,year);             //set the year
  RTC.start();
}



La data nell'RTC viene impostata tramite seriale inviando il comando "T(sec)(min)(ore)(giorno della settimana)(giorno)(mese)(anno)"

Per prima cosa non so se l'inizializzazione dell'RTC è fatta nel modo ottimale, perchè succede che a volte non riesco più a settare l'ora, riesco solo a farlo se faccio l'upload di uno sketch con l'impostazione contenuta nel Setup.
F

Federico Vanzati

Quote

//Questa funzione una volta ricevuti come parametri in ingresso la data
//restiturisce il numero del giorno dell'anno...da 0 a 365!
//Day Of Year
int calcDOY(int day, int month, int year)

  int n1 = floor(275 * month /9);
  int n2 = floor((month + 9) / 12);
  int n3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
  int n = n1 - (n2 * n3) + day - 30;

  return n;
}

void calcAlba_Tramonto( double mezzo_di_medio )
{
  float doy = (float)calcDOY(RTC.get(DS1307_DATE,false), RTC.get(DS1307_MTH,false), RTC.get(DS1307_YR,false));
  float t = doy + 284.0;

  double dec = 23.45 * sin(omega * t); //declinazione del sole
  double E1 = sin(2.0 * omega * (doy - 81.0));
  double E2 = sin(omega * (doy - 1.0));
  double E = -9.87 * E1 + 7.67 * E2; //equazione del tempo

  Serial.println((float)E,4);

  double mezzo_di = mezzo_di_medio + E/10.0;
  Serial.println((float)mezzo_di,4);

  double T = 12 * (1 + (dec * tan(lat / deg_rad) / 90.0));
  Serial.println((float)T,4);

  double mezzoT = T / 2.0;
  alba = mezzo_di - (mezzoT);
  tramonto = mezzo_di + (mezzoT);

  alba_tramonto_HM[0] = (int)floor(alba);
  alba_tramonto_HM[1] = (int)floor((alba - alba_tramonto_HM[0]) * 60.0);

  alba_tramonto_HM[2] = (int)floor(tramonto);
  alba_tramonto_HM[3] = (int)floor((tramonto - alba_tramonto_HM[2]) * 60.0);

  alba_tramonto_HM[4]++;
}

double conv_sessagimali( int ore, int minuti ){
  double min_dec = minuti / 60;

  return ore + min_dec;
}



e queste sono le due funzioni che contengono le formule astronomiche per calcolare alba e tramonto.

Se qualcuno trova degli errori, delle sviste o vuole semplicemente aiutarmi ad andare avanti io ne sarei ben felice!
A progetto concluso penso di farne una libreria (anche se non ne sono ancora in grado), e magari anche un tutorial compelto.

P.S.: nel caso qualcuno la trovasse, so che esiste una libreria già fatta che fa calcoli di questo tipo, si chiama TimeLord, ho anche provato ad aprire il file .CPP per vedere che procedimento usava e mi è sembrato più macchinoso, in più non mi sembra neanche ben documentata.
F

Federico

La libreria che hai usato per il ds1307, se e' quella che nasce da una lunga discussione sul forum di arduino, l'ho usata io parecchio e ho un oggetto che quotidianamente da parecchi mesi sembra funzionare, quindi potremmo dirla affidabile anche se ogni tanto ha delle magagne (ma che secondo me alcune volte sono dovute al cablaggio del ds1307 stesso, la mia configurazione ottimale vede un condensatore qualsiasi sulla corrente in ingresso)

Ho letto il codice ma visto che non sono esperto di calcoli solari non ho capito come si fa a capire la durata dell'alba/tramonto. Ad esempio, e' possibile capire che l'alba in quel giorno durera' mettiamo, mezz'ora? (sparo a caso, non ho idea di quanto duri una fase di alba)
Federico - Sideralis
Arduino &C: http://www.sideralis.org
Foto: http://blackman.amicofigo.com

Federico

Se posso farei anche un'altro appunto, utilizzerei variabili piu' chiare (o commentate) ed eviterei di salvare valori double per quelle di chiara comprensione, per risparmiare la memoria dell'atmega, come ad esempio mezzoT o min_dec. E' solo un appunto  8-)
Federico - Sideralis
Arduino &C: http://www.sideralis.org
Foto: http://blackman.amicofigo.com

leo72

Non ho provato l'accuratezza del tuo software.
Per aiutarti però posso intanto suggerirti questo programma in BASIC che calcola gli istanti del sorgere e del tramonto del Sole, così magari hai un modo per verificare l'accuratezza dei calcoli.

L'unica modifica da fare, secondo me, sarebbe quella di utilizzare i dati astronomici dell'orbita solare relativi all'anno 2000 dato che (sicuramente) il programma che ti ho linkato userà le effemeridi del 1960.

@federico:
non esiste una durata dell'alba o del tramonto. Essi sono "istanti": gli istanti in cui il sole compare all'orizzonte e scompare dietro di esso, rispettivamente.
Casomai esiste il "crepuscolo", quel periodo in cui la luce aumenta (o diminuisce) prima dell'alba (o dopo il tramonto).
Diciamo che una buona approssimazione è il crepuscolo civile, che si ha quando il sole è a meno di 6 gradi dall'orizzonte. Civilmente viene riconosciuto come iniziante circa 30 minuti prima dell'alba (o terminante 30 minuti dopo il tramonto): nel crepuscolo civile si ha già una quantità di luce sufficiente a compiere attività senza l'ausilio di illuminazioni artificiali, per cui in genere si potrebbe anche spegnere un eventuale impianto.

Cmq un'idea per migliorare il progetto sarebbe quella di usare in abbinamento un fotodiodo così che si regoli l'illuminazione anche in relazione alla luce effettivamente presente. Potrebbe esserci un grosso temporale che ritarda la visibilità mattutina, ad esempio.

dadebo1

ciao nn ho tempo di analizzare i calcoli .. ma ho un po' di tempo fa avevo fatto quessto programmino in vb .. che calcola la posizione del sole secondo per secondo a  seconda di coordinate date e anche l'ora di alba e tramonto (in realtà c'è uno sfasamento di pochi minuti con quelli ufficiali) ....
lo posto ..
Code: [Select]

Option Explicit
' cordinate long, lat 9.261474459325999,45.47665555555555
Const pi = 3.14159265358979
Const pi2 = 2 * pi

Dim lat As Double
Dim lon As Double
Dim anno1, annog, annoh, annom, ha, alba, albam, tram, ht, minuto, h1, d, latR, lonR, W, eTempo, difH, dd As Double
Dim aa As Double
Dim hh As Double
Dim azimut As Double
Dim azimutR As Double
Dim elevation As Double
Dim elevationR As Double
Dim declinazioneR As Double
Dim declinazione As Double
Dim elevationG, azimutG
Dim NN
Dim AB, AC, AH, AN, ASS
Dim LL As Data
Dim AGG As Double
Dim NNN As Date

Function Gdeg2dec(d, m, s As String)

 Gdeg2dec = (d + (m / 60) + (s / 3600))

End Function
Function Convert_Degree(Decimal_Deg) As Variant
  Dim degrees As Variant
  Dim minutes As Variant
  Dim seconds As Variant
 
       'Set degree to Integer of Argument Passed
       degrees = Left(Decimal_Deg, Len(Int(Decimal_Deg)))
       'Set minutes to 60 times the number to the right
       'of the decimal for the variable Decimal_Deg
       minutes = (Abs(Decimal_Deg) - Abs(degrees)) * 60
       'Set seconds to 60 times the number to the right of the
       'decimal for the variable Minute
       seconds = Format(((minutes - Int(minutes)) * 60), "0")
       'Returns the Result of degree conversion
      '(for example, 10.46 = 10~ 27  ' 36")
       Convert_Degree = " " & degrees & "° " & Int(minutes) & "' " _
           & seconds + Chr(34)

End Function




Function Arccos(x As Double) As Double

   Arccos = Atn(-x / Sqr((-x * x) + 1)) + 2 * Atn(1)

End Function

Function Arcsin(x As Double)
 
  Arcsin = Atn(x / Sqr(-x * x + 1))

End Function

'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::  this function converts radians  to decimal degrees          :::
'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Function rad2deg(rad As Double)
   rad2deg = CDbl(rad * 180 / pi)
End Function

'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
':::  this function converts decimal degrees to radians             :::
'::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Function deg2rad(deg As Double) As Double
   deg2rad = CDbl(deg * pi / 180)
End Function




Private Sub Command1_Click()
Check2 = 0
Call orologio
Call coordinate

sole.Visible = False
AH = Second(Label1) + 2
       
       Do Until Second(Label1) = AH
           Call orologio
           DoEvents
       Loop
       
Call orologio
Call coordinate
Shape1(1).Visible = False
sole.Visible = True


Do Until Check2 = 1
DoEvents
Call orologio
   If Right(Minute(Label1), 1) = 5 Or Right(Minute(Label1), 1) = 0 Then
       If Second(Label1) < 1 Then

           Call coordinate
       
       AH = Second(Label1) + 1
       
       Do Until Second(Label1) = AH
           Call orologio
           DoEvents
       Loop
       
       End If
   End If


Loop
End Sub

Private Sub Command2_Click()
Check2 = 1

End Sub


Private Sub Command3_Click()
Check2 = 1
DoEvents
AH = Second(Label1) + 2
DoEvents
       Do Until Second(Label1) = AH
           DoEvents
           Call orologio
       Loop

Unload Me

End Sub


Public Sub coordinate()
   DoEvents
       
       NN = 0
       If Check1 = 1 Then
          NN = 1
       End If
       
       lat = Val(Text8)
       lon = Val(Text9)
       
       d = 0
       anno1 = ("31/12/" & Year(Label1) - 1) 'indica l'ultimo giorno dell'anno precedente (last day of last year)
       annog = (DateDiff("D", anno1, Format(Label1, "DD/MM/YYYY")) + 1) 'giorno progressivo dell'anno (progressive day)
       annoh = (((annog - 1)) + (Hour(Label1) - NN) / 24) + (Minute(Label1) / 24 / 60) + (Second(Label1) / 24 / 60 / 60) ' ora esatta progressiva dell'anno inore copresi i sec sec (progressive time i h with minute and sec)
       
       latR = deg2rad(lat) ' in rad
       lonR = deg2rad(lon) ' in rad
       
       declinazione = 23.45 * Sin(((284) + annoh) * pi2 / (365))
       
       Text1 = declinazione
       
       declinazioneR = deg2rad(declinazione)
       DoEvents
       
       
       
       W = annoh * pi / 180
       
       eTempo = 0.42 * Cos(W) - 3.23 * Cos(2 * W) - 0.09 * Cos(3 * W) - 7.35 * Sin(W) - 9.39 * Sin(2 * W) - 0.34 * Sin(3 * W) ' equazione del tempo , corregge la nn regolarità dell'orbita terrestre (equal time)
       
       difH = 12 - ((Hour(Label1) - NN) + (Minute(Label1) / 60) + (Second(Label1) / 60 / 60)) 'ORA CORRENTE
       
       hh = 15 * (difH) - 0.25 * (eTempo - 4 * (15 - lon)) 'angolo orario (hour 's angle)
       
       
       hh = deg2rad(hh)
       
       elevationR = Arcsin((Sin(latR) * Sin(declinazioneR)) + (Cos(latR) * Cos(declinazioneR) * Cos(hh)))
       
       If elevationR < 0 Then elevationR = 0
       
       elevation = rad2deg(elevationR)
       
       
       
       elevationG = Convert_Degree(elevation)
         
       
       
       
       aa = (((Sin(elevationR) * Sin(latR)) - Sin(declinazioneR)) / (Cos(elevationR) * Cos(latR)))
       If (-aa * aa + 1) <> 0 Then
           azimutR = Arccos(aa)
       Else
           azimutR = 0
       End If
       
       'elevationG = Left(elevation, 2) & " . " & Int(Mid(elevation, 4, 2) * 0.6) & "'' . " & Mid(elevation, 6, Len(elevation))
         
          'elevationG = Int(elevation) & "° : " & (Int((elevation - Int(elevation)) * 6))
       
       
           Text3 = hh
           azimutR = azimutR
           azimut = rad2deg(azimutR)
          'If azimut <= Abs(Val(Replace(Text5, ",", "."))) Then
           If hh > 0 Then
                   
                   AN = -1
           
               Else
                    AN = 1
       
           End If
         
       azimutR = azimutR * AN
       azimut = azimut * AN
           
           
     
       azimutG = Convert_Degree(azimut)
     
     
     
     
       Text4 = elevation
       Text5 = azimut
       Text6 = elevationG
       Text7 = azimutG
       
       
       DoEvents
       
       sole.Top = (1800 - (elevation * 20)) - 240
       
       sole.Left = ((2400 + (azimut * 20))) - 240 + 150
       
       
       AB = Shape1.Count
       AC = AB + 1
     
       Do Until AB = AC
           Load Shape1(AB)
           AB = AB + 1
       Loop
     
   
       Text2 = Shape1.Count
       
       Shape1(AB - 1).Visible = True
       
       Shape1(AB - 1).Top = 1800 - (elevation * 20) - 37
       
       Shape1(AB - 1).Left = (2400 + (azimut * 20)) - 37 + 150
       
     DoEvents

End Sub

spero possa aiutare

Federico Vanzati

si, avete ragione, forse ho avuto un po' di fretta di creare discussione, i commenti ci vogliono e sono sacrosanti.

Ho creato variabili double perchè ad esempio le funzioni trigonometriche accettano float e restituiscono double, mentre la funzione floor, ho letto su nonguru che accetta in ingresso double e restituisce double.
Ho iniziato a convertirle tutte e fare alcuni cast perchè avevo notato che i valori calcolati si discostavano di parecchio rispetto a quelli simulati in matlab, così i risultati sono simili.

il procedimento che ho usato è:
- calcolare la declinazione del sole nel dato giorno
- calcolo dell'equazione del tempo, che come risultato da i minuti da sommare o sottrarre al mezzo dì medio, che per semplicità ho approssimato alle 12:10 (lo si potrebbe calcolare con precisione partendo alla longitudine)
- T invece è una formula che consente di calcolare la lunghezza del giorno, espressa in ore di luce.
- una volta conosciuto il mezzogiorno del giorno di cui si sta facendo il calcolo e la lunghezza del giorno, basta dividere per due la lunghezza del giorno e sottrarre e sommare all'istante di mezzogiorno per ottenere gli orari di alba e tramonto.

il mio calcolo (non ne sono sicurissimo) considera gli istanti di alba e tramonto quando il sole passa l'orizzonte...quindi si trova a 0°.
In realtà l'atmosfera terreste funziona tipo una lente e quindi crea un effetto ottico per cui la posizione del sole che vediamo è spostata di qualche grado...per questo motivo si usano delle correzioni da -6°  fino  a -12° per considerare i crepuscoli, a seconda che sia civile, nautico, militare ecc.

@leo72: ti ringrazio per il link, qualche paragone l'ho già fatto e so che si può migliorare parecchio, tipo applicare qualche formula di correzione tipo Gauss, solo che non so se appesantirebbe troppo arduino, pensavo che per accendere e spegnere le lucine di natale fosse accettabile anche un errore di un quarto d'ora.

@federico: grazie, i tuoi consigli sono sempre ben accetti e vedrò di metterli in pratica...come si dice...stay tuned!
F

dadebo1

si ... io avevo fatto le stesse considerazioni, ma nn trovando ( o nn volendo spendere) un dispositivo che mi dessè l'orario avevo deciso di far fare  i conti al pc poi ....
ma comunque si .. il tuo approccio è corretto.

dadebo1

ultima cosa occhio ai gradi e ai radianti .... ;)

Federico Vanzati

#9
Dec 22, 2010, 12:04 pm Last Edit: Dec 22, 2010, 12:04 pm by Federico_Vanzati Reason: 1
grazie anche a te dadebo1 ho notato che il procedimento che hai messo nel tuo programma è molto simile al mio...ma lo guarderò molto meglio stasera a casa. ;-)

la conversione in radianti è stata fatta dividendo per 57.2958 (180/pi)...m l'hai detto perchè hai trovato qualche svista?
F

dadebo1

no .. perchè mi ero scontrato anche io .
se riesco riguardo meglio il tuo listato ... e vedo se posso aiutarti ..  :)

dadebo1

mi sono perso nelle fomule.. ..ma  OMEGA? che dovrebbe essere 2pgreco/365 quando la definisci?
poi cosa intendi per equazione del tempo ?
sarebbe lo scostamento di orari che abbiamo rispetto al meridiano di riferimento ecc....
boh .... hai provato a confrontarlo con le ore di alba e tramonto.. ?
ciao


leo72

Per l'equazione del tempo vedete qui.

Anch'io ai tempi del GW-BASIC feci un programma di calcolo dei tempi di sorgere e del tramonto del sole che teneva conto di una serie di fattori quali anche (cosa che nessuno ha citato qui) la rifrazione atmosferica, vale a dire lo scostamento dell'effettivo momento di sorgere e tramonto dato dallo spostamento dei raggi luminosi per via dell'atmosfera.

Secondo me però qui stiamo entrando nel regno delle paranoie dato che per un semplice progetto per accendere/spengere delle luci tutta questa accuratezza mi pare eccessiva.  ;)
Quando hai i tempi "indicativi" del sorgere e del tramonto, approssimati anche a qualche minuto, mi pare che sia più che sufficiente, non trovate?  ::)

Federico Vanzati

@leo72:

l'equazione del tempo l'ho presa proprio da wikipedia  :)

Quote
In realtà l'atmosfera terreste funziona tipo una lente e quindi crea un effetto ottico per cui la posizione del sole che vediamo è spostata di qualche grado


non ho usato la parola rifrazione, ma comunque se ne era discusso
F

leo72

Poffarbacco!! Non mi ero accorto di quel passaggio  :D

Cmq ciò non toglie che tutta questa accuratezza mi pare eccessiva per il tipo di funzione del progetto  ::)

Io credo che se già l'approssimazione sia vicina al minuto, il risultato è più che buono.

Go Up