File content as of revision 0:ad31da30ae64:

/*
  Helios.cpp-
  Copyright (c) 2011 Hannes Hassler.  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

  This library can be used for calculating the solar position on Arduino.

  The algorithm is an adaption from
  the "PSA" solar positioning algorithm, as documented in:

  Blanco-Muriel et al.: Computing the Solar Vector. Solar Energy Vol 70 No 5 pp 431-441.
  http://dx.doi.org/10.1016/S0038-092X(00)00156-0

  According to the paper, "The algorithm allows .. the true solar vector
  to be determined with an accuracy of 0.5
  minutes of arc for the period 1999–2015.
  The original code has been downloaded from
  http://www.psa.es/sdg/sunpos.htm

  Adaptions:
  Modified calculation of number of Days since 1.Jan 2000 (dJulianDate-2451545.0)
  Neccessary because of the limited double precision on Arduino
  (double has the same precision as float on the current Arduino (2011))
  It should be used only for dates between 1.1.2000 and 31.12.2100
  (PSA itself has garantueed accuracy only until 2015)

*/

#include "Helios.h"
#include <math.h>
#include <mbed.h>


Helios::Helios(double latitude, double longitude, int tzOffset)
{
    setLocalLatitude(latitude);
    setLocalLongitude(longitude);
    setLocalTimeZoneOffset(tzOffset);
}


void Helios::updatePosition()
{
    // rtc time
    time_t rtcTime;
    struct tm *rtcTimeInfo;
    
    double dElapsedJulianDays;
    double dDecimalHours;
    double dEclipticLongitude;
    double dEclipticObliquity;
    double dRightAscension;
    double dDeclination;

    double dSin_EclipticLongitude;
    double dMeanLongitude;
    double dMeanAnomaly;
    double dOmega;

    double dGreenwichMeanSiderealTime;
    double dLocalMeanSiderealTime;
    
    double dHourAngle;
    double dCos_HourAngle;
    double dParallax;
    double dZenithAngle;
    
    // get current time from RTC
    time(&rtcTime);
    rtcTimeInfo = localtime(&rtcTime);

    // Calculate difference in days between the current Julian Day
    // and JD 2451545.0, which is noon 1 January 2000 Universal Time

    // Calculate time of the day in UT decimal hours
    dDecimalHours = (rtcTimeInfo->tm_hour - iTzOffset) + (rtcTimeInfo->tm_min
                    + rtcTimeInfo->tm_sec / 60.0 ) / 60.0;

    // Calculate current Julian Day
    long int iYfrom2000=rtcTimeInfo->tm_year-100;
    long int iA=(14-(rtcTimeInfo->tm_mon+1))/12;
    long int iM=(rtcTimeInfo->tm_mon+1)+12*iA-3;
    
    long int liAux3=(153*iM+2)/5;
    long int liAux4=365*(iYfrom2000-iA);
    long int liAux5=(iYfrom2000-iA)/4;
    
    
    dElapsedJulianDays=(double)(rtcTimeInfo->tm_mday+liAux3+liAux4+liAux5+59)+
                        -0.5+dDecimalHours/24.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)
    dOmega=2.1429-0.0010394594*dElapsedJulianDays;
    dMeanLongitude = 4.8950630+ 0.017202791698*dElapsedJulianDays; // Radians
    dMeanAnomaly = 6.2400600+ 0.0172019699*dElapsedJulianDays;
    dEclipticLongitude = dMeanLongitude + 0.03341607*sin( dMeanAnomaly )
                         + 0.00034894*sin( 2*dMeanAnomaly )-0.0001134
                         -0.0000203*sin(dOmega);
    dEclipticObliquity = 0.4090928 - 6.2140e-9*dElapsedJulianDays
                         +0.0000396*cos(dOmega);

    // 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)
    dSin_EclipticLongitude= sin( dEclipticLongitude );
    double dY1 = cos( dEclipticObliquity ) * dSin_EclipticLongitude;
    double dX1 = cos( dEclipticLongitude );
    dRightAscension = atan2( dY1,dX1 );
    if( dRightAscension < 0.0 ) dRightAscension = dRightAscension + twopi;
    dDeclination = asin( sin( dEclipticObliquity )*dSin_EclipticLongitude );

    // Calculate local coordinates ( azimuth and zenith angle ) in degrees
    dGreenwichMeanSiderealTime = 6.6974243242 +
                                 0.0657098283*dElapsedJulianDays
                                 + dDecimalHours;

    dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime*15
                              + dLongitude)*rad;
    dHourAngle = dLocalMeanSiderealTime - dRightAscension;

    dCos_HourAngle= cos( dHourAngle );
    dZenithAngle = (acos( dCos_Latitude*dCos_HourAngle
                          *cos(dDeclination) + sin( dDeclination )*dSin_Latitude));
    double dY = -sin( dHourAngle );
    double dX = tan( dDeclination )*dCos_Latitude - dSin_Latitude*dCos_HourAngle;
    dAzimuth=atan2( dY, dX );
    if ( dAzimuth < 0.0 )
        dAzimuth = dAzimuth + twopi;
    dAzimuth = dAzimuth/rad;
    // Parallax Correction
    dParallax=(dEarthMeanRadius/dAstronomicalUnit)
              *sin(dZenithAngle);
    dZenithAngle=(dZenithAngle
                  + dParallax)/rad;
    dElevation=90-dZenithAngle;
}

void Helios::setLocalLatitude(double latitude)
{
    double dLatitudeInRadians;

    dLatitude = latitude;
    dLatitudeInRadians = dLatitude*rad;
    dCos_Latitude = cos( dLatitudeInRadians );
    dSin_Latitude = sin( dLatitudeInRadians );
}

void Helios::setLocalLongitude(double longitude)
{
    dLongitude = longitude;
}

void Helios::setLocalTimeZoneOffset(int tzOffset)
{
    iTzOffset = tzOffset;
}

double Helios::azimuth()
{
    return dAzimuth;
}

double Helios::elevation()
{
    return dElevation;
}