Announcement

Collapse
No announcement yet.

Solar Position Script

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Moskus
    replied
    Thanks for this, sparkman ! I've been meaning to write this script for ages. I found the page you found, and started the cleanup.

    Luckily I found your script when I was looking for som e of those constants, so it saved me the work!

    Leave a comment:


  • sparkman
    replied
    Originally posted by dmurphy View Post
    sparkman brilliant solution, thank you so much for sharing
    Originally posted by rge View Post
    sparkman excellent work, thanks!
    You're welcome, glad you found it useful.

    Cheers
    Al

    Leave a comment:


  • rge
    replied
    sparkman excellent work, thanks!

    Leave a comment:


  • dmurphy
    replied
    sparkman brilliant solution, thank you so much for sharing

    Leave a comment:


  • WayneB
    replied
    Is so appreciate the help of this thread that I was to return the favor of a C# version of the script. Sorry abou the fomatting, don't know how to preserve it yet.

    ///---------------------------------------------------------------------------------------
    /// FileName:
    /// SolarCalculator.cs
    ///
    /// Purpose:
    /// Calculates the current solar position and stores it into two virtual devices.
    ///
    /// Remarks:
    /// Base on work by Sparkman for HomeSeer
    /// From c code posted here:
    /// http://www.psa.es/sdg/sunpos.htm
    /// Original VB.NET Conversion posted here:
    /// http://www.vbforums.com/showthread.p...ion-calculator
    ///---------------------------------------------------------------------------------------

    using System;

    ///---------------------------------------------------------------------------------------
    /// Class (SolarCalculator)
    ///---------------------------------------------------------------------------------------
    public class SolarCalculator
    {
    ///-----------------------------------------------------------------------------------
    ///<summary>
    /// Calculates the suns position. Aguments should be passed as:
    /// HomeLatitude|HomeLongitude|AzmithDeviceName|ZenithDeviceName
    ///
    /// 33.44277|-112.07275|SolarAzimuth|SolarZenith
    ///
    ///</summary>
    ///-----------------------------------------------------------------------------------
    public void CalculatePosition( object parameterObj )
    {
    string parameterString = parameterObj.ToString();

    string[] parameters = parameterString.Split( '|' );

    double latitude = double.Parse( parameters[0] );
    double longitude = double.Parse( parameters[1] );
    string azimuthDevName = parameters[2];
    string zenithDevName = parameters[3];

    DateTime utcDateTime = DateTime.UtcNow;

    double elaspedJulianDays = CalculateElaspedJulianDays( utcDateTime );

    EcllipicCoordinates ecllipicCoordinates = CalculateEcllipicCoordinates( elaspedJulianDays );

    CelestialCoordinates celestialCoordinates = CalculateCelestialCoordinates( ecllipicCoordinates );

    LocalCoordinates solarPosition = CalculateLocalCoordinates( celestialCoordinates,
    latitude,
    longitude,
    elaspedJulianDays,
    utcDateTime );

    hs.SetDeviceValueByName( azimuthDevName, solarPosition.Azimuth );
    hs.SetDeviceValueByName( zenithDevName, solarPosition.Zenith );
    }

    ///-----------------------------------------------------------------------------------
    /// <summary>
    /// Returns the elasped julian days since Julian Day (which is noon January 1, 2000
    /// universal time).
    /// </summary>
    ///-----------------------------------------------------------------------------------
    private double CalculateElaspedJulianDays( DateTime utcDateTime )
    {
    const double Julianday_NOON_JAN_1_2000 = 2451545.0;

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

    // Calculate current Julian Day
    int temp1 = ( utcDateTime.Month - 14 ) / 12;

    int temp2 = ( 1461 * ( utcDateTime.Year + 4800 + temp1 ) ) / 4 +
    ( 367 * ( utcDateTime.Month - 2 - 12 * temp1 ) ) / 12 -
    ( 3 * ( (utcDateTime.Year + 4900 + temp1 ) / 100 ) ) / 4 +
    utcDateTime.Day - 32075;

    // Calculate time of the day into UT decimal hours
    double decimalHours = RetrieveDecimalHours( utcDateTime );

    double julianDate = (double)temp2 - 0.5 + decimalHours / 24.0;

    return julianDate - Julianday_NOON_JAN_1_2000;
    }

    ///-----------------------------------------------------------------------------------
    /// <summary>
    /// Returns the eclipic coordinates for the specified elasped Julian day.
    /// </summary>
    ///-----------------------------------------------------------------------------------
    private EcllipicCoordinates CalculateEcllipicCoordinates( double elapsedJulianDays )
    {
    // 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)
    double omega = 2.1429 - 0.0010394594 * elapsedJulianDays;
    double meanLongitude = 4.895063 + 0.017202791698 * elapsedJulianDays; // Radians
    double meanAnomaly = 6.24006 + 0.0172019699 * elapsedJulianDays;

    EcllipicCoordinates ecllipicCoordinates;

    ecllipicCoordinates.Longitude = meanLongitude + 0.03341607 * Math.Sin( meanAnomaly ) + 0.00034894 * Math.Sin( 2 * meanAnomaly ) - 0.0001134 - 0.0000203 * Math.Sin(omega);
    ecllipicCoordinates.Obliquity = 0.4090928 - 0.000000006214 * elapsedJulianDays + 0.0000396 * Math.Cos(omega);

    return ecllipicCoordinates;
    }

    ///-----------------------------------------------------------------------------------
    /// <summary>
    /// Returns the celestial coordinates for the specified eclipic coordinates.
    /// </summary>
    ///-----------------------------------------------------------------------------------
    private CelestialCoordinates CalculateCelestialCoordinates( EcllipicCoordinates eclipicCoordinates )
    {
    // 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)
    double sin_EclipticLongitude = Math.Sin( eclipicCoordinates.Longitude );
    double x = Math.Cos( eclipicCoordinates.Longitude );
    double y = Math.Cos( eclipicCoordinates.Obliquity ) * sin_EclipticLongitude;

    double rightAscension = Math.Atan2( y, x );

    if ( rightAscension < 0.0 )
    {
    rightAscension = rightAscension + ( 2 * Math.PI );
    }

    CelestialCoordinates celestialCoordinates;

    celestialCoordinates.RightAscension = rightAscension;
    celestialCoordinates.Declination = Math.Asin( Math.Sin( eclipicCoordinates.Obliquity ) * sin_EclipticLongitude );

    return celestialCoordinates;
    }

    ///-----------------------------------------------------------------------------------
    /// <summary>
    /// Returns the local coordinates for the specified eclipic coordinates.
    /// </summary>
    ///-----------------------------------------------------------------------------------
    private LocalCoordinates CalculateLocalCoordinates( CelestialCoordinates celestialCoordinates,
    double latitude,
    double longitude,
    double elapsedJulianDays,
    DateTime utcDateTime )
    {
    double Km_EARTH_MEAN_RADIUS = 6371.01;
    double Km_ASTROMICAL_UNIT = 149597890;
    double Radians_PI = Math.PI / 180;

    double decimalHours = RetrieveDecimalHours( utcDateTime );
    double greenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * elapsedJulianDays + decimalHours;
    double localMeanSiderealTime = ( greenwichMeanSiderealTime * 15 + longitude ) * Radians_PI;
    double hourAngle = localMeanSiderealTime - celestialCoordinates.RightAscension;

    double latitudeInRadians = latitude * Radians_PI;

    double cos_Latitude = Math.Cos( latitudeInRadians );
    double sin_Latitude = Math.Sin( latitudeInRadians );
    double cos_HourAngle = Math.Cos( hourAngle );
    double zenithAngle = Math.Acos( cos_Latitude * cos_HourAngle * Math.Cos( celestialCoordinates.Declination ) + Math.Sin( celestialCoordinates.Declination ) * sin_Latitude );

    double dY = -Math.Sin( hourAngle );
    double dX = Math.Tan( celestialCoordinates.Declination ) * cos_Latitude - sin_Latitude * cos_HourAngle;

    double azimuth = Math.Atan2( dY, dX );
    if ( azimuth < 0.0 )
    {
    azimuth = azimuth + ( 2 * Math.PI );
    }
    azimuth = azimuth / Radians_PI;

    double parallax = ( Km_EARTH_MEAN_RADIUS / Km_ASTROMICAL_UNIT ) * Math.Sin( zenithAngle );
    double zenithAngleParallax = ( zenithAngle + parallax) / Radians_PI;

    double zenith = ( zenithAngleParallax * -1 ) + 90;

    return new LocalCoordinates( azimuth, zenith );
    }

    ///-----------------------------------------------------------------------------------
    /// <summary>
    /// Converts the hours, minutes, and seconds into decimal hour.
    /// </summary>
    ///-----------------------------------------------------------------------------------
    private double RetrieveDecimalHours( DateTime dateTime )
    {
    return dateTime.Hour + ( dateTime.Minute + dateTime.Second / 60.0 ) / 60.0;
    }

    ///-----------------------------------------------------------------------------------
    /// Inner Struct (EcllipicCoordinates)
    ///-----------------------------------------------------------------------------------
    private struct EcllipicCoordinates
    {
    // Fields
    public double Longitude;
    public double Obliquity;

    // Constructor
    public EcllipicCoordinates( double longitude, double obliquity )
    {
    Longitude = longitude;
    Obliquity = obliquity;
    }
    }

    ///-----------------------------------------------------------------------------------
    /// Inner Struct (CelestialCoordinates)
    ///-----------------------------------------------------------------------------------
    public struct CelestialCoordinates
    {
    // Fields
    public double RightAscension;
    public double Declination;

    // Constructor
    public CelestialCoordinates( double rightAscension, double declination )
    {
    RightAscension = rightAscension;
    Declination = declination;
    }
    }

    ///-----------------------------------------------------------------------------------
    /// Inner Struct (LocalCoordinates)
    ///-----------------------------------------------------------------------------------
    private struct LocalCoordinates
    {
    // Fields
    public double Azimuth;
    public double Zenith;

    // Constructor
    public LocalCoordinates( double azimuth, double zenith )
    {
    Azimuth = azimuth;
    Zenith = zenith;
    }
    }
    }

    Leave a comment:


  • langenet
    replied
    @prsmith777 Nice addition. Thanks!

    Leave a comment:


  • prsmith777
    replied
    Just found this. Love it. Thanks much.

    But I couldn't get a grasp on the Azimuth, so I added a bit to make it more understandable with compass direction in the string.
    Code:
    'From c code posted here: http://www.psa.es/sdg/sunpos.htm
    'VB.NET Conversion posted here: http://www.vbforums.com/showthread.php?832645-Solar-position-calculator
    'converted for HomeSeer use by Sparkman v1.0
    
    Imports System.Math
    
    Sub Main(ByVal Parms As String)
    Dim Debug As Boolean = False
    Dim logName As String = "Solar Position"
    Dim dLatitude As Double = xxxxx
    Dim dLongitude As Double = xxxxx
    Dim hsAzimuthDevice As Integer = 5751
    Dim hsAltitudeDevice As Integer = 5750
    
    Dim pi As Double = 3.14159265358979
    Dim rad As Double = pi / 180
    Dim dEarthMeanRadius As Double = 6371.01
    Dim dAstronomicalUnit As Double = 149597890
    
    Dim iYear As Integer = DateTime.UtcNow.Year
    Dim iMonth As Integer = DateTime.UtcNow.Month
    Dim iDay As Integer = DateTime.UtcNow.Day
    Dim dHours As Double = DateTime.UtcNow.Hour
    Dim dMinutes As Double = DateTime.UtcNow.Minute
    Dim dSeconds As Double = DateTime.UtcNow.Second
    
    Dim dZenithAngle As Double
    Dim dZenithAngleParallax As Double
    Dim dAzimuth As Double
    Dim dAltitudeAngle As Double
    Dim dElapsedJulianDays As Double
    Dim dDecimalHours As Double
    Dim dEclipticLongitude As Double
    Dim dEclipticObliquity As Double
    Dim dRightAscension As Double
    Dim dDeclination As Double
    Dim dY As Double
    Dim dX As Double
    Dim dJulianDate As Double
    Dim liAux1 As Integer
    Dim liAux2 As Integer
    Dim dMeanLongitude As Double
    Dim dMeanAnomaly As Double
    Dim dOmega As Double
    Dim dSin_EclipticLongitude As Double
    Dim dGreenwichMeanSiderealTime As Double
    Dim dLocalMeanSiderealTime As Double
    Dim dLatitudeInRadians As Double
    Dim dHourAngle As Double
    Dim dCos_Latitude As Double
    Dim dSin_Latitude As Double
    Dim dCos_HourAngle As Double
    Dim dParallax As Double
    Dim intAzimuth as Integer
    Dim dirCompass as String
    
    Try
    
    ' 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 = dHours + (dMinutes + dSeconds / 60.0) / 60.0
    ' Calculate current Julian Day
    liAux1 = (iMonth - 14) \ 12
    liAux2 = (1461 * (iYear + 4800 + liAux1)) \ 4 + (367 * (iMonth - 2 - 12 * liAux1)) \ 12 - (3 * ((iYear + 4900 + liAux1) \ 100)) \ 4 + iDay - 32075
    dJulianDate = CDbl(liAux2) - 0.5 + dDecimalHours / 24.0
    ' Calculate difference between current Julian Day and JD 2451545.0
    dElapsedJulianDays = dJulianDate - 2451545.0
    If Debug Then hs.writelog(logName,"Elapsed Julian Days Since 2000/01/01: " & CStr(dElapsedJulianDays))
    
    ' 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.895063 + 0.017202791698 * dElapsedJulianDays ' Radians
    dMeanAnomaly = 6.24006 + 0.0172019699 * dElapsedJulianDays
    dEclipticLongitude = dMeanLongitude + 0.03341607 * Math.Sin(dMeanAnomaly) + 0.00034894 * Math.Sin(2 * dMeanAnomaly) - 0.0001134 - 0.0000203 * Math.Sin(dOmega)
    dEclipticObliquity = 0.4090928 - 0.000000006214 * dElapsedJulianDays + 0.0000396 * Math.Cos(dOmega)
    If Debug Then hs.writelog(logName,"Ecliptic Longitude: " & CStr(dEclipticLongitude))
    If Debug Then hs.writelog(logName,"Ecliptic Obliquity: " & CStr(dEclipticObliquity))
    
    ' 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 = Math.Sin(dEclipticLongitude)
    dY = Math.Cos(dEclipticObliquity) * dSin_EclipticLongitude
    dX = Math.Cos(dEclipticLongitude)
    dRightAscension = Math.Atan2(dY, dX)
    If dRightAscension < 0.0 Then
    dRightAscension = dRightAscension + (2 * pi)
    End If
    dDeclination = Math.Asin(Math.Sin(dEclipticObliquity) * dSin_EclipticLongitude)
    If Debug Then hs.writelog(logName,"Declination: " & CStr(dDeclination))
    
    ' Calculate local coordinates ( azimuth and zenith angle ) in degrees
    dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * dElapsedJulianDays + dDecimalHours
    dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + dLongitude) * rad
    dHourAngle = dLocalMeanSiderealTime - dRightAscension
    If Debug Then hs.writelog(logName,"Hour Angle: " & CStr(dHourAngle))
    dLatitudeInRadians = dLatitude * rad
    dCos_Latitude = Math.Cos(dLatitudeInRadians)
    dSin_Latitude = Math.Sin(dLatitudeInRadians)
    dCos_HourAngle = Math.Cos(dHourAngle)
    dZenithAngle = (Math.Acos(dCos_Latitude * dCos_HourAngle * Math.Cos(dDeclination) + Math.Sin(dDeclination) * dSin_Latitude))
    dY = -Math.Sin(dHourAngle)
    dX = Math.Tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle
    dAzimuth = Math.Atan2(dY, dX)
    If dAzimuth < 0.0 Then
    dAzimuth = dAzimuth + (2 * pi)
    End If
    dAzimuth = dAzimuth / rad
    
    If Debug Then hs.writelog(logName,"Azimuth: " & CStr(dAzimuth))
    
    intAzimuth =CInt(dAzimuth)
    Select Case intAzimuth
    Case 0 to 11
    dirCompass = "North"
    Case 12 to 34
    dirCompass = "North North East"
    Case 35 to 56
    dirCompass = "North East"
    Case 57 to 78
    dirCompass = "East North East"
    Case 79 to 101
    dirCompass = "East"
    Case 102 to 124
    dirCompass = "East South East"
    Case 125 to 146
    dirCompass = "South East"
    Case 147 to 167
    dirCompass = "South South East"
    Case 168 to 191
    dirCompass = "South"
    Case 192 to 214
    dirCompass = "South South West"
    Case 215 to 236
    dirCompass = "South West"
    Case 237 to 259
    dirCompass = "West South West"
    Case 260 to 281
    dirCompass = "West"
    Case 282 to 304
    dirCompass = "West North West"
    Case 305 to 325
    dirCompass = "North West"
    Case 326 to 349
    dirCompass = "North North West"
    Case 350 to 360
    dirCompass = "North"
    End Select
    
    hs.setdevicevaluebyref(hsAzimuthDevice,dAzimuth ,True)
    hs.setdevicestring(hsAzimuthDevice, Math.Round(dAzimuth,1).tostring & "° " & dirCompass, True)
    
    ' Parallax Correction
    dParallax = (dEarthMeanRadius / dAstronomicalUnit) * Math.Sin(dZenithAngle)
    dZenithAngleParallax = (dZenithAngle + dParallax) / rad
    dAltitudeAngle = (dZenithAngleParallax * -1) + 90
    
    If Debug Then hs.writelog(logName,"Altitude Angle: " & CStr(dAltitudeAngle))
    hs.setdevicevaluebyref(hsAltitudeDevice,dAltitudeAngle,True)
    
    
    Catch ex As Exception
    hs.WriteLog(logName, "Exception " & ex.ToString)
    End Try
    
    End Sub
    Click image for larger version  Name:	compass.png Views:	0 Size:	81.7 KB ID:	1389592

    Leave a comment:


  • dbrunt
    replied
    That Multisensor 6 looks nice but man, it's pricey! $90 CDN! It's LUX range says 0 to 30,000 though. Is that possible?

    Thought I'd share this link I found that provides a ton of interesting information about Sunrise, Sunset, Daylength and has Astonomical/Nautical/Civil twilight,Solar Noon, Moon Phases, Moon Rise & Moon Set, Azimuths and Altitudes: https://www.timeanddate.com/sun/cana...couver?month=6
    You can enter your longitude/latitude or pick location on a map...

    Leave a comment:


  • dbrunt
    replied
    I too have an ST815. It's mounted on outside patio just a few feet away from my south facing windows so I can differentiate sunny vs cloudy days. I'm currently using it's lux to open/close them rather than sunrise/sunset so they open and close based on morning and evening luminance. Currently I am using <40 and >40 and disabling the triggered event and enabling its twin event so each event only triggers once a day.

    However, I want to reopen the SE blind earlier and as the sun wanes westward reopen the S blind followed later by the SW blind so I can't use the single outside lux, hence my interest in this scripting. I don't want to hard code times and values as I want to automagically tune them from April to October.

    Leave a comment:


  • sparkman
    replied
    Originally posted by dmiller View Post
    I found it easier to put a multisensor 6 between the blinds and the windows. If luminance is over 1500 lux the sun is shining through the widow. I'll probably add some logic to close the blinds at a lower light level if the TV is on.

    I find that trees and partly cloudy days make projecting luminosity problematic.

    The multisensor is the only device I found that could be pointed out the windows and not immediately max out the sensor. Most devices clip before 1000 lux, which doesn't allow differentiation between a sunless sky and direct sunshine.
    I use some Everspring ST815 luminance sensors as they go up to 3000lux.

    Leave a comment:


  • dmiller
    replied
    I found it easier to put a multisensor 6 between the blinds and the windows. If luminance is over 1500 lux the sun is shining through the widow. I'll probably add some logic to close the blinds at a lower light level if the TV is on.

    I find that trees and partly cloudy days make projecting luminosity problematic.

    The multisensor is the only device I found that could be pointed out the windows and not immediately max out the sensor. Most devices clip before 1000 lux, which doesn't allow differentiation between a sunless sky and direct sunshine.

    Leave a comment:


  • dbrunt
    replied
    Yes, that helps greatly. Mine are venetian so 0 is closed down, 50 is flat and 99 is closed up so my adjustment is from 50 to 99 and then back to 50. I may have to tweak the script but I'll get it figured out. Thanks again,
    Daniel

    Leave a comment:


  • sparkman
    replied
    Originally posted by dbrunt View Post
    Hi Al,

    Many thanks for this great solution! I've been struggling for a few days with how to auto-adjust my 3 south facing iBlinds controlled window blinds to shield the sun and just discovered this thread. I've implemented the virtual devices and scripts and will see tomorrow how the blinds react. I have a south facing Everspring outdoor luminance meter to test for LUX as a group condition (>2500) so that the blinds don't close on cloudy days and also only if date is > spring equinox and < fall equinox.

    One question though is how does the third parameter 'value for adjusting the altitude to determine the shade position' affect the blind positioning? I have one that is SE, one is S and one is SW so I want to adjust them differently. SE would tilt up earlier and then back to 50% by solar noon, S would fully close up just after solar noon and then gradually return to 50% somewhat before sunset and SW would fully close later in the day returning to 50% closer to sunset. Maybe if I sleep on it, things will be clearer tomorrow...
    You're welcome. That value is a multiplier that multiplies the solar altitude with whatever you use as a multiplier to create a value to be used as the position value for the shade. For my shades, 0 is closed and 99 open. So if the solar altitude was 33 degrees and the multiplier is 2, then the shade position would be set at 66. For my mid-day events, I use a value of 1.4. The multiplier is determined by trial and error and will not be consistent throughout the day, so you may need different events for different portions of the day. Most of my shade events have conditions for a range of solar azimuths. Also,I don't use the script for all my shade events, most events just set the shade to a specific value when the solar azimuth and altitude are in the ranges set by the event conditions. Hope that helps, if not, let me know.

    Cheers
    Al

    Leave a comment:


  • dbrunt
    replied
    Hi Al,

    Many thanks for this great solution! I've been struggling for a few days with how to auto-adjust my 3 south facing iBlinds controlled window blinds to shield the sun and just discovered this thread. I've implemented the virtual devices and scripts and will see tomorrow how the blinds react. I have a south facing Everspring outdoor luminance meter to test for LUX as a group condition (>2500) so that the blinds don't close on cloudy days and also only if date is > spring equinox and < fall equinox.

    One question though is how does the third parameter 'value for adjusting the altitude to determine the shade position' affect the blind positioning? I have one that is SE, one is S and one is SW so I want to adjust them differently. SE would tilt up earlier and then back to 50% by solar noon, S would fully close up just after solar noon and then gradually return to 50% somewhat before sunset and SW would fully close later in the day returning to 50% closer to sunset. Maybe if I sleep on it, things will be clearer tomorrow...

    Leave a comment:


  • langenet
    replied
    Thanks again Al. I have the living room and office south facing as well. I'm inspired to put this to use soon!

    How often is the event called?

    Leave a comment:

Working...
X