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!
Announcement
Collapse
No announcement yet.
Solar Position Script
Collapse
X
-
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;
}
}
}
- Likes 1
Leave a comment:
-
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
- Likes 1
Leave a comment:
-
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:
-
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:
-
Originally posted by dmiller View PostI 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:
-
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:
-
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:
-
Originally posted by dbrunt View PostHi 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...
Cheers
Al
Leave a comment:
-
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:
-
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:
Leave a comment: