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