Search
The Renewable Energy site for Do-It-Yourselfers
Sun Class Source
Code
This is the source code for Sun
Class.
This is just to take a look at before
you download -- if you want to download it, then download all the solar classes
here: Solar Tools Source
Sun Class:
'Solar Analysis Tools -
'Copyright (C) 2005 Gary Reysa (gary@BuildItSolar.com)
' This program is free software; you can redistribute it and/or modify it
' under the terms of the GNU General Public License as published by the Free
' Software Foundation; either version 2 of the License, or (at your option)
' any later version.
' This program is distributed in the hope that it will be useful and encourage
' further development of solar tools, but WITHOUT ANY WARRANTY; without even
' the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
' See the GNU General Public License for more details.
' Sun Class -- Instances of Sun represent the sun
' Given the inputs: Date, Time, Longitude, Latitude and Altitude (abv sea level)
' The sun object will provide solar position and clear day solar radiation data.
' It is assumed that the time inputs to sun class are LOCAL SOLAR TIME.
'
' The calcs are based on the Lunde book, Solar Thermal Engineering.
'
' Azimuth is measured from true South.
' + azimuths are to East, - to West
' (see note below)
'
' NOTE: Azimuth is measured from Due South.
' Internally, Sun Class uses the convention of + azimuth to West, and - to East
' (this is the convention Lunde used in this calc in his book)
' Since this is counter to + rotations being CCW when viewed from the + Z direction,
' I always change the Azimuth angles sign, and reverse the Y component on the
' ToSunVector calcs before making them available as properties.
'
' Checked a number of dates and times -- seems to work well -- returns data
' that agrees well with tables in books.
' Should take a look at whether the "year" variable is needed at all -- probably not
' The calcs are based on the Lunde book, Solar Thermal Engineering.
Imports System.Math
Public Class SunCls
Public Const RadPerDeg As Single = 1 / 57.296
Public Const DegPerRad As Single = 57.296
Public Const Isc As Single = 429.0 ' Solar constant - solar inten at top of atmosphere BTU/hr-ft^2
' Constants used in solar calcs:
' Parameters used to estimate solar radiation intensity
' From Lunde table 3-2, pg 65
' Data is by month, jan index is 0, dec index is 11
Public ReadOnly DayNo() As Double = _
{19.85, 54.06, 80.0, 110.47, 140.15, 172.5, 201.84, 232.49, 265.0, 292.34, 324.2, 357.5}
Public ReadOnly Decl() As Double = _
{-20.0, -10.0, 0.0, 11.6, 20.0, 23.45, 20.6, 12.3, 0.0, -10.5, -19.8, -23.45}
Public ReadOnly A() As Double = _
{390, 385, 376, 360, 350, 345, 344, 351, 365, 378, 387, 391}
Public ReadOnly B() As Double = _
{0.142, 0.144, 0.156, 0.18, 0.196, 0.205, 0.207, 0.201, 0.177, 0.16, 0.149, 0.142}
Public ReadOnly C() As Double = _
{0.058, 0.06, 0.071, 0.097, 0.121, 0.134, 0.136, 0.122, 0.092, 0.073, 0.063, 0.057}
Public ReadOnly EqnT() As Double = _
{-11.2, -13.9, -7.5, 1.1, 3.3, -1.4, -6.2, -2.4, 7.5, 15, 4, 13.8, 1.6}
'Input vars
' Changing any of these vars requrires "recalc"
Private mDayOfMonth As Integer ' 1 to 31
Private mMonth As Integer ' 1 is jan, 2 feb ...
Private mYear As Integer ' the year - pick a non-leap year for default
Private mASTHr As Double ' Solar time, 12 is solar noon
Private mLongitude As Double ' Longitude -- West is + -- radians
Private mLatitude As Double ' Latitude -- + in NA -- radians
Private mAltitude As Double ' Altitude in ft
'Output vars
' All of the these vars are changed by a "recalc"
Private mDayOfYear As Integer ' day of the year
Private mHrAng As Double ' hour angle in radians, is measured from south, + to West (pm)
Private mDeclination As Double ' Declination - apparent angle of earths axis of rot to ecliptic
Private mZenithAng As Double ' The complement of the solar altitude -- i.e. 90 - AltAng radians
Private mAzimuthAng As Double ' Azimuth angle of sun -- measured from South, + to West
' NOTE this is + to West internally (to match Lunde book), but the
' property AzimuthAngle is reported as + to east (CCW about Z azis)
Private mAltAng As Double ' Solar altitude -- radians
Private mSunVec As New VectorCls ' A vector pointing toward sun (determined from AltAng and Azimuth
Private mIoext As Double ' solar radiation at top of atmosphere
Private mIdn As Double ' Direct normal radiation BTU/hr-ft^2
Private mIdiffuse As Double ' Diffuse radiation BTU/hr=ft^2
Private mItotal As Double ' Total solar radiation (Idn + Idiffuse)
Private mIHorz As Double ' Total solar radiation on horz surface
Private mSunRise As Double ' Hour of sunrise in solar time
Private mSunSet As Double ' Hour of sunset in solar time
Private UpToDate As Boolean
' Constructors
Sub New(ByVal Day As Integer, ByVal Month As Integer, ByVal Year As Integer, _
ByVal ASThr As Double, ByVal Longitude As Double, ByVal Latitude As Double, _
ByVal Altitude As Double)
mDayOfMonth = Day
mMonth = Month
mYear = Year
mASTHr = ASThr
mLongitude = Longitude
mLatitude = Latitude
mAltitude = Altitude
UpToDate = False
End Sub
Sub New()
mDayOfMonth = 21
mMonth = 12
mYear = 2003
mASTHr = 12
mLongitude = 0
mLatitude = 45 * RadPerDeg
mAltitude = 0
UpToDate = False
End Sub
' ---- Input Variable Gets and Sets for inputs
' ReCalc is called if any of these inputs are set
Public Property Year() As Integer
Get
Year = mYear
End Get
Set(ByVal Value As Integer)
mYear = Value
UpToDate = False
End Set
End Property
' Month number -- Jan = 1, ... Dec = 12
Public Property Month() As Integer
Get
Month = mMonth
End Get
Set(ByVal Value As Integer)
If (Value < 13) And (Value > 0) Then
mMonth = Value
UpToDate = False
End If
End Set
End Property
Public Property DayOfMonth() As Integer
Get
DayOfMonth = mDayOfMonth
End Get
Set(ByVal Value As Integer)
If (Value < 32) And (Value > 0) Then
mDayOfMonth = Value
UpToDate = False
End If
End Set
End Property
' ASThr -- solar time in hours -- 12hr = solar noon (sun is due south)
Public Property ASThr() As Double
Get
ASThr = mASTHr
End Get
Set(ByVal Value As Double)
mASTHr = Value
UpToDate = False
End Set
End Property
Public Property Longitude() As Double
Get
Longitude = mLongitude
End Get
Set(ByVal Value As Double)
mLongitude = Value
UpToDate = False
End Set
End Property
Public Property Latitude() As Double
Get
Latitude = mLatitude
End Get
Set(ByVal Value As Double)
mLatitude = Value
UpToDate = False
End Set
End Property
' Altitude above sea level in ft
Public Property Altitude() As Double
Get
Altitude = mAltitude
End Get
Set(ByVal Value As Double)
mAltitude = Value
UpToDate = False
End Set
End Property
' ---- Output Variable Gets ------------
' Hour Angle -- measured from South, + to East in radians (2 pi radians per day)
' NOTE: that internally HrAng is + to West, but this property is reported
' out as + to east (CCW about Z axis)
Public ReadOnly Property HrAng() As Double
Get
If Not UpToDate Then ReCalc()
HrAng = -mHrAng
End Get
End Property
' DayOfYear-- )
Public ReadOnly Property DayOfYear() As Double
Get
If Not UpToDate Then ReCalc()
DayOfYear = mDayOfYear
End Get
End Property
'Solar Declination angle (rad) -- apparent angle of earchs axis of rotation to Earth Sun line
' 23.45 deg at winter solstice, -23.45 at summer solstice, 0 at equinoxs
Public ReadOnly Property Declination() As Double
Get
If Not UpToDate Then ReCalc()
Declination = mDeclination
End Get
End Property
' Zenith Angle (rad) -- Is the complement of Altitude angle
Public ReadOnly Property ZenithAng() As Double
Get
If Not UpToDate Then ReCalc()
ZenithAng = mZenithAng
End Get
End Property
' Azimuth Angle (rad) -- Measure from South, + to East
' NOTE: that internally mAzimuthAngle is + to West, but this property is reported
' out as + to east (CCW about Z axis)
Public ReadOnly Property AzimuthAng() As Double
Get
If Not UpToDate Then ReCalc()
AzimuthAng = -mAzimuthAng
End Get
End Property
' Altitude Angle -- Elevation angle of sun (rad)
Public ReadOnly Property AltAng() As Double
Get
If Not UpToDate Then ReCalc()
AltAng = mAltAng
End Get
End Property
' ToSunVector -- a unit vector in the direction of the sun
' The Y comonpent of this vector is flipped to keep the convention
' + azimuth to East externally, and + az to West internally
Public ReadOnly Property ToSunVector() As VectorCls
Get
If Not UpToDate Then ReCalc()
Dim V As New VectorCls(mSunVec.X, -mSunVec.Y, mSunVec.Z)
ToSunVector = V
End Get
End Property
' Ioext -- solar radiation at top of atmosphere
Public ReadOnly Property Ioext() As Double
Get
If Not UpToDate Then ReCalc()
Ioext = mIoext
End Get
End Property
' Direct solar radiation -- (adjusted for date, altitude, and airmass) BTU/hr ft^2
Public ReadOnly Property Idn() As Double
Get
If Not UpToDate Then ReCalc()
Idn = mIdn
End Get
End Property
' Diffuse radiation -- BTU/hr ft^2
Public ReadOnly Property Idiffuse() As Double
Get
If Not UpToDate Then ReCalc()
Idiffuse = mIdiffuse
End Get
End Property
' Total solar radiation -- sum of Idn + Idiffuse BTU/hr ft^2
Public ReadOnly Property Itotal() As Double
Get
If Not UpToDate Then ReCalc()
Itotal = mItotal
End Get
End Property
' Total solar radiation on Horz Surface -- sum of Idn + Idiffuse BTU/hr ft^2
Public ReadOnly Property Ihorz() As Double
Get
If Not UpToDate Then ReCalc()
Ihorz = mIHorz
End Get
End Property
'Sunrise hour -- time of sunrise in solar time
Public ReadOnly Property Sunrise() As Double
Get
If Not UpToDate Then ReCalc()
Sunrise = mSunRise
End Get
End Property
' Sunset hour -- time of sunset in solar time
Public ReadOnly Property Sunset() As Double
Get
If Not UpToDate Then ReCalc()
Sunset = mSunSet
End Get
End Property
' ReCalc -- this routine must be called after any of the input vars in changed by user.
' It recalculates all of the output solar parameters to agree with current inputs
' The calculations are based on the book "Solar Thermal Engineering", Lunde
Public Sub ReCalc()
UpToDate = True ' moved this here to prevent recursive calls to ReCalc
' Day of Year
Dim mDate As Date
mDate = DateSerial(mYear, mMonth, mDayOfMonth)
mDayOfYear = mDate.DayOfYear
' Ioext is solar radiation at top of atmosphere --
' it is based on the solar constant corrected for earth to sun distance
mIoext = Isc * (1 + 0.033 * Cos((360 * mDayOfYear / 370.0) * RadPerDeg))
' Declination -- The apparent tilt of the earths axis of rotation to the ecliptic
' +23.45 at summer solstice, -23.45 at winter solstice, and 0 at the equinoxes
' returns declination for the day of year in radians
mDeclination = RadPerDeg * 23.45 * Sin(((mDayOfYear - 80) / 370) * 360 * RadPerDeg)
'Hour Angle --
' found by taking the time between solar noon and current AST, and
' multiplying by 15 deg/hr, then converting to radians.
mHrAng = (mASTHr - 12.0) * (360 / 24) * RadPerDeg
' Altitude Angle -- the solar altitude in radians
mAltAng = Asin(Cos(mLatitude) * Cos(mDeclination) * Cos(mHrAng) _
+ Sin(mLatitude) * Sin(mDeclination))
' Zenith Angle
mZenithAng = 90 * RadPerDeg - mAltAng
' Solar Azimuth -- angle of sun from due south in radians -- + to west
' this is the one that is dual valued
' Solar Azimuth -- angle of sun from due south in radians -- + to west
' this is the one that is dual valued
Dim sinAzAng As Single = Cos(mDeclination) * Sin(mHrAng) / Cos(mAltAng)
mAzimuthAng = Asin(sinAzAng)
' Alt at azimuth of 90
Dim alt90 As Double
alt90 = Asin(Sin(mDeclination) / Sin(mLatitude))
' if the sun higher than alt90, then mAzimuth is ok, else is 180 - mAz
If (mAltAng < alt90) Then
If mAzimuthAng < 0 Then
mAzimuthAng = -180 * RadPerDeg - mAzimuthAng
Else
mAzimuthAng = 180 * RadPerDeg - mAzimuthAng
End If
End If
' Calc the ToSun Vector
mSunVec = VectorCls.AzElevVector(mAzimuthAng, mAltAng)
' Direct Normal Irradiation -- in BTU/hr - ft^2
' This is direct normal irradation on sunny day corrected for air mass and
' for altitude
Dim A, B As Double
Dim Palt_to_Psl As Double
Palt_to_Psl = Exp(-0.0000361 * mAltitude)
A = Aparm()
B = Bparm()
' Check for nightime
If mAltAng > 0 Then
mIdn = A * Exp(-Palt_to_Psl * B / Cos(mZenithAng))
Else
mIdn = 0
End If
' Diffuse solar radiantion
' Calced as table factor "C" times the direct normal
Dim C As Double
C = Cparm()
mIdiffuse = Idn() * Cparm()
' Toal solar radiation -- is sum of direct normal (Idn) and diffuse (Idiffuse)
' BTU/hr - ft^2
mItotal = mIdn + mIdiffuse
' Total radiation radiation on Horz surface
mIHorz = mIdn * Sin(mAltAng) + mIdiffuse
' Sunrise hour
Dim SRHourAng As Double
SRHourAng = Acos(-Tan(mLatitude) * Tan(mDeclination))
mSunRise = 12.0 - SRHourAng / (15 * RadPerDeg)
mSunSet = 12.0 + SRHourAng / (15 * RadPerDeg)
End Sub ' ReCalc
' interpolations for A as function of day of year
' these routines are called by ReCalc to get interpolated values for
' the parms in the table
Friend Function Aparm() As Double
Dim H, L As Double
Dim frac As Double
If (mDayOfYear < 19.85) Then
L = 391
H = 390
frac = mDayOfYear / 19.85
ElseIf (mDayOfYear >= 357) Then
L = 391
H = 391
frac = 0
Else
Dim i As Integer = 0
While (mDayOfYear > DayNo(i))
i = i + 1
End While
L = A(i - 1)
H = A(i)
frac = (mDayOfYear - DayNo(i - 1)) / (DayNo(i) - DayNo(i - 1))
End If
Aparm = L + (H - L) * frac
End Function
' interpolations for B as function of day of year
' B is Air Mass factor from the B array above
Friend Function Bparm() As Double
Dim H, L As Double
Dim frac As Double
If (mDayOfYear < 19.85) Then
L = 0.142
H = 0.142
frac = 0
ElseIf (mDayOfYear > 357) Then
L = 0.142
H = 0.142
frac = 0
Else
Dim i As Integer = 0
While (mDayOfYear > DayNo(i))
i = i + 1
End While
L = B(i - 1)
H = B(i)
frac = (mDayOfYear - DayNo(i - 1)) / (DayNo(i) - DayNo(i - 1))
End If
Bparm = L + (H - L) * frac
End Function
' interpolations for C as function of day of year
' C is a factor from the C array above
Friend Function Cparm() As Double
Dim H, L As Double
Dim frac As Double
If (mDayOfYear < 19.85) Then
L = 0.058
H = 0.057
frac = mDayOfYear / 19.85
ElseIf (mDayOfYear > 357) Then
L = 0.057
H = 0.057
frac = 0
Else
Dim i As Integer = 0
While (mDayOfYear > DayNo(i))
i = i + 1
End While
L = C(i - 1)
H = C(i)
frac = (mDayOfYear - DayNo(i - 1)) / (DayNo(i) - DayNo(i - 1))
End If
Cparm = L + (H - L) * frac
End Function
' interpolations for Equation of Time as function of day of year
Friend Function EqnTparm() As Double
Dim H, L As Double
Dim frac As Double
If (mDayOfYear < 19.85) Then
L = 1.6
H = -11.2
frac = mDayOfYear / 19.85
ElseIf (mDayOfYear > 357) Then
L = 1.6
H = 1.6
frac = 0
Else
Dim i As Integer = 0
While (mDayOfYear > DayNo(i))
i = i + 1
End While
L = EqnT(i - 1)
H = EqnT(i)
frac = (mDayOfYear - DayNo(i - 1)) / (DayNo(i) - DayNo(i - 1))
End If
EqnTparm = L + (H - L) * frac
End Function
Public Sub Dump()
If Not UpToDate Then ReCalc()
Console.WriteLine("Sun Object:")
Console.WriteLine("Date, DayOfYr {0,2} {1,2} {2,4} DOY {3,4}", Month, DayOfMonth, Year, DayOfYear) ' As Integer = 21 ' 1 to 31
Console.WriteLine("ASThr {0,4:F1}", ASThr) ' As Double = 12.0 ' Solar time, 12 is solar noon
Console.WriteLine("Long, Lat, Altitude {0,4:F1} {1,4:F1} {2,4:F1}", _
Longitude * DegPerRad, Latitude * DegPerRad, Altitude) ' As Double = 111.0 * RadPerDeg ' Longitude -- West is + -- radians
'Output vars
' All of the these vars are changed by a "recalc"
Console.WriteLine("Declination, HrAngle {0,4:F1} {1,4:F1}", Declination * DegPerRad, HrAng * DegPerRad) ' hour angle in radians, is measured from south, + to West (pm)
Console.WriteLine("Azimuth, Altitude {0,5:F1} {1,5:F1}", AzimuthAng * DegPerRad, AltAng * DegPerRad) ' Azimuth angle of sun -- measured from South, + to West internally,
' but reported out as + to East
Console.WriteLine("ToSunVec X,Y,Z {0,6:F3} {1,6:F3} {2,6:F3}", ToSunVector.X, ToSunVector.Y, ToSunVector.Z) ' to sun vector
Console.WriteLine("Ioext {0,6:F1}", Ioext) ' solar radiation at top of atmosphere
Console.WriteLine("Idn, Idif, Itot {0,6:F1} {1,6:F1} {2,6:F1}", Idn, Idiffuse, Itotal) ' Direct normal radiation BTU/hr-ft^2
Console.WriteLine("Ihorz {0,6:F1}", Ihorz) ' Total solar radiation on horz surface
Console.WriteLine("SunRise, Sunset {0,4:F1} {1,4:F1}", Sunrise, Sunset) ' Hour of sunrise in solar time
End Sub
End Class ' sun