Слежение за солнцем SUNPOS

cheese
Offline
Зарегистрирован: 08.05.2012
Помогите разобраться. Часы идут, а расчет стоит на месте.

//Sun Position Calculation by Mowcius (mowcius.co.uk)
//Provides sun position (relative) from static variables

#include <math.h>
#include <Wire.h>
#include "RTClib.h"
#define pi    3.14159265358979323846
#define twopi (2*pi)
#define rad   (pi/180)
#define EarthMeanRadius     6371.01      // In km
#define AstronomicalUnit    149597890      // In km

RTC_DS1307 RTC;

//Input Variables --------------------- TIME HAS TO BE IN UT (UNIVERSAL TIME)! NO TIME ZONES OR SUMMER TIMES --------
//My last modifications were probably at this time on this date!
     

     float Longitude = 5.311; //enter longitude here
     float Latitude = 5.007; //enter latitude here
//--------

//Program Variables
     float ZenithAngle;
     float Azimuth;
       float RightAscension;
     float Declination;
       float Parallax;
       float ElevationAngle;

     float ElapsedJulianDays;
     float DecimalHours;
     float EclipticLongitude;
     float EclipticObliquity;
//--------

void setup() {
Serial.begin(57600);
Wire.begin();
RTC.begin();
RTC.adjust(DateTime(__DATE__, __TIME__));

}

void sunPos() {
       DateTime now = RTC.now();

       int Year = (now.year(), DEC); //year
     int Month = (now.month(), DEC); //month
     int Day = (now.day(), DEC); //day
     float Hours = (now.hour(), DEC); //hour
     float Minutes = (now.minute(), DEC); //minutes      
        
     // Auxiliary variables
     float dY;
     float dX;

     // Calculate difference in days between the current Julian Day 
     // and JD 2451545.0, which is noon 1 January 2000 Universal Time
     
           float JulianDate;
           long int liAux1;
           long int liAux2;
           // Calculate time of the day in UT decimal hours
           DecimalHours = Hours + (Minutes / 60.0);
           // Calculate current Julian Day
           liAux1 =(Month-14)/12;
           liAux2=(1461*(Year + 4800 + liAux1))/4 + (367*(Month 
                 - 2-12*liAux1))/12- (3*((Year + 4900 
           + liAux1)/100))/4+Day-32075;
           JulianDate=(float)(liAux2)-0.5+DecimalHours/24.0;
           // Calculate difference between current Julian Day and JD 2451545.0 
           ElapsedJulianDays = JulianDate-2451545.0;
     
     // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the 
     // ecliptic in radians but without limiting the angle to be less than 2*Pi 
     // (i.e., the result may be greater than 2*Pi)
     
           float MeanLongitude;
           float MeanAnomaly;
           float Omega;
           Omega=2.1429-0.0010394594*ElapsedJulianDays;
           MeanLongitude = 4.8950630+ 0.017202791698*ElapsedJulianDays; // Radians
           MeanAnomaly = 6.2400600+ 0.0172019699*ElapsedJulianDays;
           EclipticLongitude = MeanLongitude + 0.03341607*sin( MeanAnomaly ) 
                 + 0.00034894*sin( 2*MeanAnomaly )-0.0001134
                 -0.0000203*sin(Omega);
           EclipticObliquity = 0.4090928 - 6.2140e-9*ElapsedJulianDays
                 +0.0000396*cos(Omega);
     
     // Calculate celestial coordinates ( right ascension and declination ) in radians 
     // but without limiting the angle to be less than 2*Pi (i.e., the result may be 
     // greater than 2*Pi)
     
           float Sin_EclipticLongitude;
           Sin_EclipticLongitude= sin( EclipticLongitude );
           dY = cos( EclipticObliquity ) * Sin_EclipticLongitude;
           dX = cos( EclipticLongitude );
           RightAscension = atan2( dY,dX );
           if( RightAscension < 0.0 ) RightAscension = RightAscension + twopi;
           Declination = asin( sin( EclipticObliquity )*Sin_EclipticLongitude );
     
     // Calculate local coordinates ( azimuth and zenith angle ) in degrees
     
           float GreenwichMeanSiderealTime;
           float LocalMeanSiderealTime;
           float LatitudeInRadians;
           float HourAngle;
           float Cos_Latitude;
           float Sin_Latitude;
           float Cos_HourAngle;
           GreenwichMeanSiderealTime = 6.6974243242 + 
                 0.0657098283*ElapsedJulianDays 
                 + DecimalHours;
           LocalMeanSiderealTime = (GreenwichMeanSiderealTime*15 
                 + Longitude)*rad;
           HourAngle = LocalMeanSiderealTime - RightAscension;
           LatitudeInRadians = Latitude*rad;
           Cos_Latitude = cos( LatitudeInRadians );
           Sin_Latitude = sin( LatitudeInRadians );
           Cos_HourAngle= cos( HourAngle );
           ZenithAngle = (acos( Cos_Latitude*Cos_HourAngle
                 *cos(Declination) + sin( Declination )*Sin_Latitude));
           dY = -sin( HourAngle );
           dX = tan( Declination )*Cos_Latitude - Sin_Latitude*Cos_HourAngle;
           Azimuth = atan2( dY, dX );
           if ( Azimuth < 0.0 ) 
             Azimuth = Azimuth + twopi;
           Azimuth = Azimuth/rad;
           // Parallax Correction
           Parallax=(EarthMeanRadius/AstronomicalUnit)
                 *sin(ZenithAngle);
           ZenithAngle=(ZenithAngle //Zenith angle is from the top of the visible sky (thanks breaksbassbleeps)
                 + Parallax)/rad;
               ElevationAngle = (90-ZenithAngle); //Retrieve useful elevation angle from Zenith angle
}

void loop(){
 DateTime now = RTC.now();
 sunPos(); //Run sun position calculations
 Serial.print("Elevation Angle:  ");
 Serial.println(ElevationAngle); //Print Elevation (Vertical) with no decimal places as accuracy is not really great enough
 Serial.print("Azimuth:  ");
 Serial.println(Azimuth); //Print Azimuth (Horizontal) with no decimal places
// Serial.print("Seconds:  ");
// Serial.println(now.second(), DEC);
// печать времени в сериал
   Serial.print(now.year(), DEC);
    Serial.print('/');
    Serial.print(now.month(), DEC);
    Serial.print('/');
    Serial.print(now.day(), DEC);
    Serial.print(' ');
    Serial.print(now.hour(), DEC);
    Serial.print(':');
    Serial.print(now.minute(), DEC);
    Serial.print(':');
    Serial.print(now.second(), DEC);
    Serial.println();
 
 if(ElevationAngle < 0)
 Serial.println("The sun has set. Get some sleep!");
 delay(1000); //Delay 10 seconds - Values aren't going to have changed anyway as they are currently static variables!
}

 

leshak
Offline
Зарегистрирован: 29.09.2011

 "Расчет, стой! ать-два".

Какой расчет, расчет чего именно?

Вообще если есть проблемма к с какой-то длинной логикой, то общее решение: натыкать Serial.print где только можно, выводить не только конечный результат расчета, но и промежуточные решения. Продублировать расчет руками на бумажки, ну и смотреть с какого момента "ожидаемое" не совпадает с фактическим.

P.S. А зачем так сложно?  Не проще несколько датчиков освещения и следить за фактическим положением солнца. Где светит - там и солнце ;)

 

AK54
Offline
Зарегистрирован: 15.03.2013

Пример программы для расчета высоты солнца над горизонтом по указанному времени.
Данные вводятся при помощи монитора порта. Формат данных описан в тексте программы.
Расчёт производят функции Nday и heightSun.
Результат выдаётся функцией showAngle отрицательное значение – темно, положительное светло.
Алгоритм отлаживался на Borland C++ 5.5.
Данные совпали с калькулятором http://www.hmn.ru/index.php?index=41&value=1
При переводе на Arduino 1.0.1 разница составила 4 или 5-ть минут, что может быть связано с использованием двойной точности на Borland C++.
Недостатки :
Большое количество констант занимает много оперативной памяти.

 

// отладочная программа расчёта высоты солнца над горизонтом
// ПО АЛГОРИТМУ ИЗЛОЖЕННОМУ В http://www.astronomy.ru/forum/index.php/topic,21684.0.html
/* команды
 D(01-31)(01-12)(10-99) устанавливает дату на которую считаем
 D121112 12 декабря 2012 года 

 T(00-24).(00-59) местное время
 Н4 4-ый часовой пояс
 X00.00 широта (latitude) только для северной широты
 У00.00 долгота (longitude) только для восточной долготы
 X55.75  Москва
 У37.62
 R сбросить все буфера и установить умолчание
 P печать исходных данных
 С вычислить угол над горизонтом
*/

byte day, month, year ; // ден месяц и две последнии цифры года расчёта. Проверки на коректность ввода нет !
float LocalTime ;   // местное время в часах минуты десятичные
float TimeLocal ;   // местное время в часах минуты обычные
byte timeZone  ;    // временная зона
float gorizont ;    // угол на котором считается что расвет и закат наступил
                     // с учётом рефракции и углового размера солнца, а также
                     // с учётом угла погружения солнца при гражданских, административных и астрономических сумерках              
float latitude ;    // введёная широта
float longitude ;   // введёная  долгота в градусах
float angle ;       // угол солнца над горизонтом
//
//
void readBuf() ;       // ввод дданных через терминал и разбор ввода
void setDefault() ;    // установим умолчание
void listPar() ;       // печатаем в монитор исходные данные 
long Nday(const byte &D, byte M, const byte &Y)  ;// Вычисление модифицированной  юлианской даты
float heightSun(const  long &MD , const int  &TZ , const float &LT , const float &LON , const float &LAT ) ;  // расчёт высоты солнца на горизонтом
void showAngle(float angle , float gorizont) ; // показывает угол солнца над горизонтом

//
//    разбор ввода
// 

void readBuf() // ввод дданных через терминал и разбор ввода
  {
    byte buf=0 ;
    byte command=0 ;
    if (Serial.available())
    {
    command = Serial.read();
    if  (command == 68)     // D- дата  
        {
         day=(Serial.read() - 48) * 10 + (Serial.read() - 48);
         month=(Serial.read() - 48) * 10 + (Serial.read() - 48);
         year=(Serial.read() - 48) * 10 + (Serial.read() - 48);
         Serial.print("DAY= "); Serial.print(day, DEC); Serial.print("  MONTH= "); Serial.print(month, DEC) ; Serial.print("  YEA= ");  Serial.println(year, DEC); 
        }
    if  (command == 84)  // T ввели время
        {     
         byte h= (Serial.read() - 48) * 10 + (Serial.read() - 48);
         buf=Serial.read(); // убрали разделитель
         byte m= (Serial.read() - 48) * 10 + (Serial.read() - 48);
         LocalTime=h+(m/60.0);
         TimeLocal=h+(m/100.0) ;
         Serial.print("LocalTime= "); Serial.println(LocalTime, DEC);
         Serial.print("TimeLocal= "); Serial.println(TimeLocal, DEC);
        } 
    if  (command == 72)  // H- ввели часовой пояс  ( не учитываем западную долготу, которая должна быть с минусом и только один знак от 0 до 9 )     
        {
         timeZone=(Serial.read() - 48);
         Serial.print("TIMEZONE= "); Serial.println(timeZone, DEC);
        }
    if  (command == 71)  // G- угол над горизонтом
        {
         gorizont = (Serial.read() - 48) * 10 + (Serial.read() - 48); 
         buf =  Serial.read(); // убрали разделитель
         buf = (Serial.read() - 48) * 10 + (Serial.read() - 48);
         gorizont=gorizont+(buf/100.0);
         Serial.print("GORIZONT= "); Serial.println(gorizont, DEC); 
        }
    if  (command == 88)  // X введена широта
        {     
         byte x= (Serial.read() - 48) * 10 + (Serial.read() - 48);
         buf=Serial.read(); // убрали разделитель
         byte xm= (Serial.read() - 48) * 10 + (Serial.read() - 48);
         latitude=x+(xm/100.0);
         Serial.print(" LATITUDE= "); Serial.println(latitude, DEC);
        } 
      if (command ==89 ) // Y введена долгота
         {byte y=(Serial.read() - 48) * 10 + (Serial.read() - 48);
         buf=Serial.read() ; // убрали разделитель
         byte ym=(Serial.read() - 48) * 10 +(Serial.read() - 48);
         longitude = y +(ym/100.0);
         Serial.print(" LONGITUDE= "); Serial.println(longitude,DEC);
         }
     if  (command == 82)     // R установим умолчание
          { 
          setDefault ();
          listPar();
          }
     if  (command == 80)  listPar() ; // P печать исходных данных 
         //
     if  (command == 67)   // 'C' вычислить УГОЛ НАД ГОРИЗОНТОМ если угол отрицательный, то темно
            {  
             long MD=Nday( day,  month,  year) ; // вычислим порядковый день года
             // Serial.println(""),Serial.print("otladka MD = ");  Serial.println(MD , DEC ) ; Serial.println("") ; // отладка
             angle = heightSun ( MD , timeZone , LocalTime , longitude , latitude ) ; // вычислим угол солнца над горизонтом
             showAngle(angle , gorizont ) ;
     } 
         }
     }
void setDefault() // установим умолчание
   {
    latitude = 55.83  ; // широта Москвы 
    longitude=  36.61  ; // долгота Москвы 
    day = 22 ;
    month = 6 ;
    year = 12 ;
    timeZone=4 ;  // часовой пояс Москва
    gorizont=90.86 ; // угол от зенита с учётом угловоно размера солнца и рефракции. Расвет или закат наступает когла солнце неже горизонта на 0.86 градуса
    LocalTime=12.50 ; 
    TimeLocal=12.30 ;
  }
void listPar()     // печатаем в монитор исходные данные 
   {
    Serial.print("VER = 1.2 ");  Serial.println(" ");
    Serial.print("day : month : year = ");  Serial.print(day, DEC) ; Serial.print(" : "); 
                                            Serial.print(month, DEC) ; Serial.print(" : ");
    Serial.print("20");                     Serial.print(year, DEC) ; Serial.println(" ");  
    Serial.print("         LocalTime = ");  Serial.println(LocalTime, DEC) ;
    Serial.print("         TimeLocal = ");  Serial.println(TimeLocal, DEC) ;
    Serial.print("          timeZone = ");  Serial.println(timeZone, DEC) ;  
    Serial.print("          latitude = ");  Serial.println(latitude, DEC) ; 
    Serial.print("         longitude = ");  Serial.println(longitude, DEC) ;
    Serial.print("          gorizont = ");  Serial.println(gorizont, DEC) ;
   }
//
// *************************************************************************************** 
//  
void showAngle(float angle, float gorizont)
   {     
    float a ;
    a=gorizont-angle ;
    Serial.print("height of the sun above the horizon = "); Serial.println(a,DEC);   
   }
//
// *************************************************************************************** 
//  
long Nday(const byte &D, byte M, const byte &Y) //    Вычисление модифицированной юлианской даты 
    {
      long Ye =2000+Y ;            //        храним две последние цифры года
      if (M<=2) { M=M+12 ; Ye=Ye-1; }
      //
      long VAR0 = 365 * Ye ;
      long VAR1 = (VAR0 - 679004L );
      long VAR2 = (Ye / 400 - Ye / 100 + Ye / 4) ;
      long VAR3 = 306001L * (M + 1) / 10000L ;
      long MD= VAR1+VAR2+VAR3+D ;
      //
      /*   отладка 
      Serial.print("M    = ");  Serial.println(M, DEC) ;
      Serial.print("Ye   = ");  Serial.println(Ye, DEC) ;
      Serial.print("VAR0 = ");  Serial.println(VAR0, DEC) ;
      Serial.print("VAR1 = ");  Serial.println(VAR1, DEC) ;
      Serial.print("VAR2 = ");  Serial.println(VAR2, DEC) ;
      Serial.print("VAR3 = ");  Serial.println(VAR3, DEC) ;
      Serial.print("MD   = ");  Serial.println(MD, DEC) ;  
      return (365 * Ye - 679004L ) + (Ye / 400 - Ye / 100 + Ye / 4) + 306001L * (M + 1) / 10000L + D ; \\ так не работает из-за преобразования в типах данных 
      */
      return MD ;
    }
//
float heightSun(const long &MD , const int  &TZ , const float &LT , const float &LON , const float &LAT )
      {
 // в процедуре много констант которые будут занимать оперативную память желательно перенести их а память программ 
        //
        float Rg=PI/180.0 ; // можно вычислять в радианах, но пример был в градусах 
        float To ,SO , ST ;
          {        
             float a1 = 24110.54841;  float a2 = 8640184.812 ; float a3 = 0.093104 ; float a4 = 0.0000062 ;
             To = (MD - 51544.5) / 36525.0 ; //  мод.юл.дата на начало суток в юлианских столетиях
             SO = a1 + a2 * To + a3 * pow(To , 2) - a4 * pow(To , 3) ;
      //
             //  UT = (LT - TZ) ;                 всемирное время в часах от полуночи/N
             //  Nsec= (LT - TZ)*3600 ;      Nsec =UT * 3600 ‘ количество секунд, прошедших
             //                              от начала суток до момента наблюдения всемирноe времени
             //                              NsecS = Nsec * 366.2422 / 365.2422 ;     количество звездных секунд
      float NsecS = (LT - TZ) * 3609.8564 ; 
      //                       NsecS = (LT - TZ) * 3609,85647
             //                                                 3609.8564 = 3600 * 366.2422 / 365.2422  
      //         SG = (SO + NsecS) / 3600.0 * 15.0 ;    гринвическое среднее звездное время в градусах
             //  
      float SG= (SO + NsecS)*0.004166666 ;  //  15/3600.0 = 0,004166666
      ST=fmod((SG+LON),360.0) ;
   }
 //
        //                         3. Вычисление эклиптических координат Солнца  
        //
 float M = 357.528 + 35999.05 * To + 0.04107 * (LT - TZ) ; // средняя аномалия
 M=fmod(M,360.0) ;
 float LO = 280.46 + 36000.772 * To + 0.04107 * (LT - TZ) ;
        float L = LO + (1.915 - 0.0048 * To) * sin(M*Rg) + 0.02 * sin(2 *M*Rg) ; //   долгота Солнца      
 L=fmod(L,360.0) ; // приведение угла в диапазон 0-360 градусов
        //  
 //                                       раcчёт высоты солнца на горизонтом 
 //
        //                                    X = cos(L)     Y = sin(L)   Z = 0.0       
        //
        //                          4.Координаты Cолнца в прямоугольной экваториальной системе координат
        //
        float Ra , Dec ;
          {
          float Eps = 23.439281*Rg ; //    наклон эклиптики к экватору в радианах
          float Xprim = cos(L*Rg) ; //    вектор в экваториальной системе координат
          float Yprim = sin(L*Rg) * cos(Eps) ; // - Z * sin(Eps) ' Z=0   
          float Zprim = sin(L*Rg) * sin(Eps) ; // + Z * cos(Eps) ' Z=0   
   Ra=atan2( Yprim ,Xprim ) *180.0/PI ;  // 
          Dec=atan2(Zprim,sqrt(Xprim*Xprim + Yprim*Yprim))*180.0/PI ;  
   }  
        //                          5.Экваториальные геоцентрические координаты Солнца
 //    float Ra=atan(sin(L*Rg) / cos(L*Rg) ) *180/PI ;
 //    float Dec=atan(sin(L*Rg)/cos(L*Rg))*180/PI ;
 //
 // Можно и топоцентрические вычислить, но разница несущественна - несколько угловых секунд.
        //
        //                          6. Азимутальные координаты Солнца
 //
        float Th = ST - Ra ; //  часовой угол
 float Cosz=sin(LAT*Rg)*sin(Dec*Rg) + cos(LAT*Rg)*cos(Dec*Rg)*cos(Th*Rg) ;
        return acos(Cosz)*180.0/PI ;
      }
//
//---------------------------------------------------------------------------------------------------------------------------//
//
void setup()
  {
   Serial.begin(9600); 
   setDefault();
   listPar(); 
  }
  //
void loop()
  {
  readBuf();
  delay(50) ;   
  } 

 

storm134
Offline
Зарегистрирован: 14.04.2017

А как тут еще азимутальный угол относительно севера высчитать?