Sun Angle Script

The following rhinoscript code calculates and draws a series of sun angles based on the time, day and latitude.
NOTE: in my test there’s something happening to the angles around noon causing them to be inaccurate because the calculations are based on cos values. All other times besides noon should be fine.

Updated version with nicer UI from Ezio. Thank you!
Download:sunbatchrender_ezb.rvb

Option Explicit
'Script written by cwwang.com augmented by algorithmicdesign.net
'Script version Sunday, March 01, 2009 9:39:09 PM

Call Main()
Sub Main()

	Dim i,j,k
	Dim hour
	Dim arrSettings : arrSettings   = getSunSettings()
	Dim longitude	: longitude		= CDbl(arrSettings(0))
	Dim latitude	: latitude		= CDbl(arrSettings(1))
	Dim offset		: offset		= CDbl(arrSettings(2))
	'day number --> January 1 = day 1, January 2 = day 2, etc.
	'March equinox 78
	'September equinox 264
	'June solstice 171
	'Dec solstice 354
	Dim startDay	: startDay		= CDbl(arrSettings(3))
	Dim endDay		: endDay		= CDbl(arrSettings(4))
	Dim dayIncrement: dayIncrement 	= CDbl(arrSettings(5))
	Dim startHour   : startHour     = CDbl(arrSettings(6))
	Dim endHour		: endHour		= CDbl(arrSettings(7))
	Dim increment	: increment		= CDbl(arrSettings(8))
	'Dim increment	: increment		=.1666666666666666666666666666666' 10/60=.16666 every ten  minutes
	'Dim increment	: increment		=.01666666666666666666666666666666' 1/60=.016666 every minute
	Dim pathFolder 	: pathFolder 	= Rhino.BrowseForFolder(,"Pick a destination folder", "Output Folder")
	Dim nTimes		: nTimes     	=(endHour-startHour)/increment
	Dim nDays		: nDays      	=(endDay-startDay)/dayIncrement

	ReDim analemmaPoints(CInt(nTimes),CInt(nDays))
	'ReDim analemmaPoints(150,50)

	Call Rhino.EnableRedraw(False)

	Dim count:count=0
	For i=startDay To endDay

		Dim count2:count2=0
		Dim paz:paz=0
		Dim switchaz:switchaz=False

		hour=startHour
		For k=0 To (endHour-hour)/increment

			Dim N:N=i 

			Dim g:g = (360/365.242199)*(N + hour/24.0)
			'Rhino.Print(g)

			'calculation of the declination of the sun
			Dim	D:D = 0.396372-22.91327*cos(Rhino.ToRadians(g))+4.02543*sin(Rhino.ToRadians(g))-0.387205*cos(Rhino.ToRadians(2*g))+0.051967*sin(Rhino.ToRadians(2*g))-0.154527*cos(Rhino.ToRadians(3*g)) + 0.084798*sin(Rhino.ToRadians(3*g))
			'Rhino.Print("Declination "&d)

			'calculate the time correction For solar angle:
			Dim TC:TC = 0.004297+0.107029*cos(Rhino.ToRadians(g))-1.837877*sin(Rhino.ToRadians(g))-0.837378*cos(Rhino.ToRadians(2*g))-2.340475*sin(Rhino.ToRadians(2*g))
			'Rhino.Print("Time Correction " &tc)

			'calculate the Solar Hour Angle (SHA)
			Dim SHA:SHA = (hour+offset-12)*15 + Longitude + TC
			'rhino.Print("Solar Hour Angle "&sha)

			'Note that If SHA Is greater than 180, Then you must add (-360) To the
			'result And If SHA Is lower than -180, Then you must add 360 To the
			If (SHA>180) Then
				SHA=SHA-360
			ElseIf(SHA<-180)Then
				SHA=SHA+360
			End If
			'calculate the Sun Zenith Angle (SZA):
			Dim SZA:SZA = (sin(Rhino.ToRadians(Latitude))*sin(Rhino.ToRadians(D)) )+ (cos(Rhino.ToRadians(Latitude))*cos(Rhino.ToRadians(D))*cos(Rhino.ToRadians(SHA)))
			'rhino.Print("cos(sza) "&sza)

			'NOTE: If cos(SZA) formula gives a figure greater than 1 use 1 And If
			'it gives you a figure lower than -1 use -1.
			'Dim SSZA:SSZA=SZA
			If(SZA>=1) Then
				SZA=1
			ElseIf(SZA<=-1) Then
				SZA=-1
			End If

			sza=ArcCos(sza)
			'SZA Is the complementary angle of the Sun Elevation Angle Or
			Dim sea:sea=90-sza

			'calculate the Azimuth Angle (AZ):
			Dim AZ
			AZ = (sin(Rhino.ToRadians(D))-sin(Rhino.ToRadians(Latitude))*cos(Rhino.ToRadians(SZA) ) )/(cos(Rhino.ToRadians(Latitude))*sin(Rhino.ToRadians(SZA))) 

			az=arccos(az)

			If(paz>az And switchaz=False)Then
				switchAz=True
			End If
			paz=az

			If(switchaz=True)Then
				az=0-az
			End If

			'rhino.Print("Azimuth Angle "& az)
			'rhino.Print("Zenith Angle "& sza)
			'rhino.Print("Altitude Angle "& sea)

			Dim ray:ray=Rhino.AddLine( Array(0,0,0),array(0,40,0) )
			ray= Rhino.RotateObject(ray,array(0,0,0),sea,array(1,0,0),False)
			ray= Rhino.RotateObject(ray,array(0,0,0),az,array(0,0,-1),False)

			analemmaPoints(count2,count)=rhino.CurveEndPoint(ray)
			Dim sun : sun = Rhino.AddDirectionalLight (rhino.CurveEndPoint(ray), rhino.CurveStartPoint(ray))
			Call BatchRender(i, k, pathFolder)
			Call Rhino.DeleteObject(sun)
			'Call rhino.AddText(cstr(hour),analemmaPoints(count),.04)
			count2=count2+1

			hour=hour+increment
			'Rhino.Print(hour)

		Next

		i=i+dayIncrement

		count=count+1

	Next
	'Call rhino.AddLine(analemmaPoints(0,0),analemmaPoints(0,1) )

	For i=0 To count2-1
		For j=0 To count-2			

			'Dim analemmaCurve:analemmaCurve=rhino.AddCurve(analemmaPoints(3),4 )
			Call rhino.AddLine(analemmaPoints(i,j),analemmaPoints(i,j+1) )
		Next	

	Next

	Call Rhino.EnableRedraw(True)
	Rhino.Print("---")
End Sub

Function ArcCos(ByVal X)
	If X <> 1 Then
		ArcCos =rhino.ToDegrees( Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1) )
	Else
		ArcCos = 0
	End If
End Function

Sub BatchRender(strDay, strHour, pathFolder)

	Dim FileName : FileName = Chr(34) & pathFolder & strDay & "_" & strHour & ".jpg" & Chr(34)
	Rhino.Command "-_Render "
	Rhino.Command "-_SaveRenderWindowAs " & FileName & " _Enter"
	Rhino.Command "-_CloseRenderWindow "
End Sub

Function getSunSettings()

	getSunSettings	= Null

	Dim arrSettings(8), arrPrompt(8)
	arrPrompt(0)   = "Longitude"
	arrPrompt(1)   = "Latitude"
	arrPrompt(2)   = "GMT"
	arrPrompt(3)   = "StartDay"
	arrPrompt(4)   = "EndDay"
	arrPrompt(5)   = "DayIncrement"
	arrPrompt(6)   = "StartHour"
	arrPrompt(7)   = "EndHour"
	arrPrompt(8)   = "HourIncrement"

	arrSettings(0) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(0))
	If IsNull(arrSettings(0)) Then arrSettings(0) = -73.997
	arrSettings(1) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(1))
	If IsNull(arrSettings(1)) Then arrSettings(1) = 40.73
	arrSettings(2) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(2))
	If IsNull(arrSettings(2)) Then arrSettings(2) = 5
	arrSettings(3) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(3))
	If IsNull(arrSettings(3)) Then arrSettings(3) = 0
	arrSettings(4) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(4))
	If IsNull(arrSettings(4)) Then arrSettings(4) = 1
	arrSettings(5) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(5))
	If IsNull(arrSettings(5)) Then arrSettings(5) = 5
	arrSettings(6) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(6))
	If IsNull(arrSettings(6)) Then arrSettings(6) = 6.0
	arrSettings(7) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(7))
	If IsNull(arrSettings(7)) Then arrSettings(7) = 18.0
	arrSettings(8) = Rhino.GetSettings(Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(8))
	If IsNull(arrSettings(8)) Then arrSettings(8) = 1
	'Dim increment	: increment		=.1666666666666666666666666666666' 10/60=.16666 every ten  minutes
	'Dim increment	: increment		=.01666666666666666666666666666666' 1/60=.016666 every minute
	'day number --> January 1 = day 1, January 2 = day 2, etc.
	'March equinox 78
	'September equinox 264
	'June solstice 171
	'Dec solstice 354

	Dim arrReturnSettings : arrReturnSettings =	Rhino.PropertyListBox (arrPrompt,  arrSettings, "Please enter the appropriate values" ,"Sun Batch Render Settings")
	Dim i
	For i = 0 To Ubound(arrSettings)
		'write these settings to a file so that you don't have to re-enter them next time
		Rhino.SaveSettings Rhino.InstallFolder & "SunRender.ini", "SunRender", arrPrompt(i), arrReturnSettings(i)
	Next

	getSunSettings = arrReturnSettings

End Function