next topic »

[on:]

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:

#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.

[Reply #1 on:]

//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.

[Reply #2 on:]

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)

[Reply #3 on:]

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-)

[Reply #4 on:]

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.

[Reply #5 on:]

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:

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

[Reply #6 on:]

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!

[Reply #7 on:]

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.

[Reply #8 on:]

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

[Reply #9 on:]

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?

[Reply #10 on:]

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

[Reply #11 on:]

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
 

[Reply #12 on:]

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?  :

[Reply #13 on:]

@leo72:

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

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

[Reply #14 on:]

Poffarbacco!! Non mi ero accorto di quel passaggio  

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.