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