Solar Position Script

[sparkman]

Al this is very interesting. So you call this via an event and pass in the three parameters? How often is this run? How is the multiplier parameter determined? Hit/miss.
I will be purchasing some Bali blinds very soon so I'm keenly interested in implementing this.. Would you mind showing your event?

As well, do you use Alexa and do you know if you can tell Alexa to open/close the blinds natively. I suppose using Jon's helper that's possible but wasn't sure if Alexa supports this.

Thanks,

Robert

Hi Robert, yes, called from events with 3 parameters. See the example events for what triggers them. The multiplier can be a bit of trial and error, but the numbers I picked initially are pretty close. There's also a minimum change that's hard coded in the script (4%) so that the blinds don't adjust too often. I don't use Alexa but assume that it will work as the child devices look like switches and dimmers inside HS.

This event deals with the morning and early afternoon sun in a south facing window. I'll modify it when it gets a bit warmer out to close as well if no one is home, whereas in the winter, I typically want it open to get more of the solar heat. I may just base it on forecasted high temps.



This event deal with sun later in the afternoon for the same window. There's a TV in the corner of the room and this prevents the sun shining on it, so one of the conditions is that the TV is on (right now it's assumed that when the receiver is on the TV is on).

[sparkman]

This works great! Thanks for posting it. Now, I have to figure out how to use it to close my motorized blinds. Elliott

Glad that it's working for you. I just posted some examples on how I use it, along with the script in this post: https://forums.homeseer.com/forum/de...42#post1286742

[langenet]

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?

[dbrunt]

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

[sparkman]

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

[dbrunt]

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

[dmiller]

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.

[sparkman]

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.

[dbrunt]

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.

[dbrunt]

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

[prsmith777]

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.

'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

[langenet]

@prsmith777 Nice addition. Thanks!

[WayneB]

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;
}
}
}

[dmurphy]

sparkman brilliant solution, thank you so much for sharing

[rge]

sparkman excellent work, thanks!