Here’s a link to NREL’s PvWatts:
http://rredc.nrel.gov/solar/calculators/PVWATTS/version2/
It mainly a program for PV panel sizing but also gives heat energy in
kWh/m^2 for many locations in the US, which is handy for sizing solar
concentrators. For us old Americans 1kWh = 3412.142 BTUs/h and 1
meter^2 = 10.76391 ft^2, so 3412 / 10.76391 = 317 BTU/ft^2/h.
And here is the source code for two solar tracking programs, the first
was written in C by NREL and includes a VB port I did, this NREL
program has a high degree of accuracy but is more CPU intensive, while
the second program (not written by NREL) is both less accurate and CPU
intensive.
Have fun,
Curbie
Option Compare Database
Option Explicit
'=============================================================================
' Contains:
' S_solpos (computes solar position and intensity from time and
place)
' INPUTS: via posdata struct) year, daynum, hour, minute,
second, latitude, longitude, timezone, intervl
' OPTIONAL: (via posdata struct) month, day, press, temp, tilt,
aspect, function
' OUTPUTS: EVERY variable in the struct posdata (defined in
solpos.h)
'
' NOTE: Certain conditions exist during which some of
the output variables are undefined or cannot be
' calculated. In these cases, the variables are
returned with flag values indicating such. In other
' cases, the variables may return a realistic, though
invalid, value. These variables and the flag values
' or invalid conditions are listed below:
'
' amass -1.0 at zenetr angles greater than 93.0
degrees
' ampress -1.0 at zenetr angles greater than 93.0
degrees
' azim invalid at zenetr angle 0.0 or latitude
+/-90.0 or at night elevetr limited to -9 degrees at
' night
' etr 0.0 at night
' etrn 0.0 at night
' etrtilt 0.0 when cosinc is less than 0
' prime invalid at zenetr angles greater than 93.0
degrees
' sretr +/- 2999.0 during periods of 24 hour sunup
or sundown
' ssetr +/- 2999.0 during periods of 24 hour sunup
or sundown
' ssha invalid at the North and South Poles unprime
invalid at zenetr angles greater than 93.0 degrees
' zenetr limited to 99.0 degrees at night
'
' S_init (optional initialization for all input parameters in
the posdata struct)
' INPUTS: struct posdata*
' OUTPUTS: struct posdata*
'
' (Note: initializes the required S_solpos INPUTS above
to out-of-bounds conditions, forcing the user to
' supply the parameters; initializes the OPTIONAL
S_solpos inputs above to nominal values.)
'
' S_decode (optional utility for decoding the S_solpos return
code)
' INPUTS: long integer S_solpos return value, struct posdata*
' OUTPUTS: text to stderr
'
' Usage: In calling program, just after other 'includes',
insert:
'
' #include "solpos00.h"
'
' Function calls:
' S_init(struct posdata*) [optional]
' .
' .
' [set time and location parameters before S_solpos
call]
' .
' .
' int retval = S_solpos(struct posdata*)
' S_decode(int retval, struct posdata*) [optional]
' (Note: you should always look at the S_solpos return
value, which contains error codes. S_decode is one
' option for examining these codes. It can also serve
as a template for building your own application-
' specific decoder.)
'
' Martin Rymes
' National Renewable Energy Laboratory
' 25 March 1998
'
' 27 April 1999 REVISION: Corrected leap year in S_date.
' 13 January 2000 REVISION: SMW converted to structure posdata
parameter
' and subdivided into functions.
' 01 February 2001 REVISION: SMW corrected ecobli calculation
' (changed sign). Error is small (max
0.015 deg
' in calculation of declination angle)
'-----------------------------------------------------------------------------
'#include <math.h>
'#include <string.h>
'#include <stdio.h>
'#include "solpos00.h"
'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
'
' Structures defined for this module
'
'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Private Type trigdata ' used to pass
calculated values locally
cd As Single ' cosine of the
declination
ch As Single ' cosine of the hour
angle
cl As Single ' cosine of the latitude
sd As Single ' sine of the
declination
Sl As Single ' sine of the latitude
End Type
'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
'
' Temporary global variables used only in this file:
'
'+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
'cumulative number of days prior to beginning of month 0=non-leap
year, 1=leap year
Private month_days(1, 12) As Integer ' day of month array
Private array_init As Boolean ' day of month array
initalization flag
Private Const degrad As Single = 57.295779513 ' converts from
radians to degrees
Private Const raddeg As Single = 0.0174532925 ' converts from
degrees to radians
'=============================================================================
' Long integer function S_solpos, adapted from the VAX solar libraries
'
' This function calculates the apparent solar position and the
intensity of the sun (theoretical maximum solar energy)
' from time and place on Earth.
'
' Requires (from the struct posdata parameter):
' Date and time:
' year
' daynum (requirement depends on the S_DOY switch)
' month (requirement depends on the S_DOY switch)
' day (requirement depends on the S_DOY switch)
' hour
' minute
' second
' interval DEFAULT 0
' Location:
' latitude
' longitude
' Location/time adjuster:
' timezone
' Atmospheric pressure and temperature:
' press DEFAULT 1013.0 mb
' temp DEFAULT 10.0 degrees C
' Tilt of flat surface that receives solar energy:
' aspect DEFAULT 180 (South)
' tilt DEFAULT 0 (Horizontal)
' Function Switch (codes defined in solpos.h)
' function DEFAULT S_ALL
'
' Returns (via the struct posdata parameter):
' everything defined in the struct posdata in solpos.h.
'-----------------------------------------------------------------------------
Public Function S_solpos(pdat As posdata) As Long
Dim retval As Long ' return value
Dim tdat As trigdata ' trig data
' initialize the trig structure
tdat.sd = -999# ' flag to force
calculation of trig data
tdat.cd = 1# '
tdat.ch = 1# ' set the rest of these
to something safe
tdat.cl = 1# '
tdat.Sl = 1# '
If Not array_init Then InitArray ' if day of month array
not initalized, do it
If ((retval = validate(pdat)) <> 0) Then ' if position data
validates
S_solpos = retval ' return validate flags
End If ' if position data
validates
If (pdat.function And L_DOY) Then ' if day of year
conversion required
doy2dom pdat ' convert input doy to
month-day
Else ' else convert day of
month
dom2doy pdat ' convert input
month-day to doy
End If ' if day of year
conversion required
If (pdat.function And L_GEOM) Then ' if basic geometry
calculations required
geometry pdat ' do basic geometry
calculations
End If ' if basic geometry
calculations required
If (pdat.function And L_ZENETR) Then ' if etr at
non-refracted zenith angle required
zen_no_ref pdat, tdat ' do etr at
non-refracted zenith angle
End If ' if etr at
non-refracted zenith angle required
If (pdat.function And L_SSHA) Then ' if Sunset hour
calculation required
ssha pdat, tdat ' do Sunset hour
calculation required
End If ' if Sunset hour
calculation required
If (pdat.function And L_SBCF) Then ' if ShadoWBRSd
correction factor required
sbcf pdat, tdat ' do ShadoWBRSd
correction factor
End If ' if ShadoWBRSd
correction factor required
If (pdat.function And L_TST) Then ' if true solar time
required
Tst pdat ' do true solar time
End If ' if true solar time
required
If (pdat.function And L_SRSS) Then ' if sunrise/sunset
calculations required
srss pdat ' if sunrise/sunset
calculations
End If ' if sunrise/sunset
calculations required
If (pdat.function And L_SOLAZM) Then ' if solar azimuth
calculations required
SAzm pdat, tdat ' do solar azimuth
calculations
End If ' if solar azimuth
calculations required
If (pdat.function And L_REFRAC) Then ' if atmospheric
refraction calculations required
Refrac pdat ' do atmospheric
refraction calculations
End If ' if atmospheric
refraction calculations required
If (pdat.function And L_AMASS) Then ' if airmass
calculations required
amass pdat ' do airmass
calculations required
End If ' if airmass
calculations required
If (pdat.function And L_PRIME) Then ' if kt-prime/unprime
calculations required
prime pdat ' do kt-prime/unprime
calculations
End If ' if kt-prime/unprime
calculations required
If (pdat.function And L_ETR) Then ' if ETR and ETRN
(refracted) required
etr pdat ' do ETR and ETRN
(refracted)
End If ' if ETR and ETRN
(refracted) required
If (pdat.function And L_TILT) Then ' if tilt calculations
required
Tilt pdat ' do ETR and ETRN
(refracted)
End If ' if tilt calculations
required
S_solpos = 0 ' done
End Function
'============================================================================
' Void function S_init
'
' This function initiates all of the input parameters in the struct
posdata passed to S_solpos(). Initialization is
' either to nominal values or to out of range values, which forces
the calling program to specify parameters.
'
' NOTE: This function is optional if you initialize ALL input
parameters in your calling code. Note that the required
' parameters of date and location are deliberately initialized out of
bounds to force the user to enter real-world
' values.
'
' Requires: Pointer to a posdata structure, members of which are
initialized.
'
' Returns: Void
'-----------------------------------------------------------------------------
Public Sub S_init(pdat As posdata)
pdat.Day = -99 ' Day of month (May 27 =
27, etc.)
pdat.daynum = -999 ' Day number (day of
year; Feb 1 = 32 )
pdat.Hour = -99 ' Hour of day, 0 - 23
pdat.Minute = -99 ' Minute of hour, 0 - 59
pdat.Month = -99 ' Month number (Jan = 1,
Feb = 2, etc.)
pdat.second = -99 ' Second of minute, 0 -
59
pdat.Year = -99 ' 4-digit year
pdat.interval = 0 ' instantaneous
measurement interval
pdat.aspect = 180# ' Azimuth of panel
surface (direction it faces) N=0, E=90, S=180, W=270
pdat.Latitude = -99# ' Latitude, degrees
north (south negative)
pdat.Longitude = -999# ' Longitude, degrees
east (west negative)
pdat.press = 1013# ' Surface pressure,
millibars
pdat.solcon = 1367# ' Solar constant, 1367
W/sq m
pdat.Temp = 15# ' Ambient dry-bulb
temperature, degrees C
pdat.Tilt = 0# ' Degrees tilt from
horizontal of panel
pdat.TimeZone = -99# ' Time zone, east (west
negative).
pdat.sbwid = 7.6 ' Eppley shadow band
width
pdat.sbrad = 31.7 ' Eppley shadow band
radius
pdat.sbsky = 0.04 ' Drummond factor for
partly cloudy skies
pdat.function = S_ALL ' compute all parameters
End Sub
'============================================================================
' Local long int function validate
'
' Validates the input parameters
'----------------------------------------------------------------------------
Private Function validate(pdat As posdata) As Long
Dim retval As Long ' return value
retval = 0 ' start return value
with no errors
If (pdat.function And L_GEOM) Then ' if basic geometry
calculations required
If ((pdat.Year < 1950) Or (pdat.Year > 2050)) Then
retval = retval Or ShL32F(1, S_YEAR_ERROR)
End If ' if invalid algoritm
year limits
If (Not (pdat.function And S_DOY) And ((pdat.Month < 1) Or
(pdat.Month > 12))) Then
retval = retval Or ShL32F(1, S_MONTH_ERROR)
End If ' if invalid month
If (Not (pdat.function And S_DOY) And ((pdat.Day < 1) Or (pdat.Day
retval = retval Or ShL32F(1, S_DAY_ERROR)
End If ' if invalid day month
If ((pdat.function And S_DOY) And ((pdat.daynum < 1) Or
(pdat.daynum > 366))) Then
retval = retval Or ShL32F(1, S_DOY_ERROR)
End If ' if invalid day of year
If ((pdat.Hour < 0) Or (pdat.Hour > 24)) Then
retval = retval Or ShL32F(1, S_HOUR_ERROR)
End If ' if vainvalidlid hour
of day
If ((pdat.Minute < 0) Or (pdat.Minute > 59)) Then
retval = retval Or ShL32F(1, S_MINUTE_ERROR)
End If ' if invalid mintue of
hour
If ((pdat.second < 0) Or (pdat.second > 59)) Then
retval = retval Or ShL32F(1, S_SECOND_ERROR)
End If ' if invalid second of
minute
If ((pdat.Hour = 24) And (pdat.Minute > 0)) Then
retval = retval Or (ShL32F(1, S_HOUR_ERROR) Or ShL32F(1,
S_MINUTE_ERROR))
End If ' if more than 24 hrs
(addditional minutes)
If ((pdat.Hour = 24) And (pdat.second > 0)) Then
retval = retval Or (ShL32F(1, S_HOUR_ERROR) Or ShL32F(1,
S_SECOND_ERROR))
End If ' if more than 24 hrs
(addditional seconds)
If (CSng(Abs(pdat.TimeZone)) > 12#) Then
retval = retval Or ShL32F(1, S_TZONE_ERROR)
End If ' if invalid time zone
If ((pdat.interval < 0) Or (pdat.interval > 28800)) Then
retval = retval Or ShL32F(1, S_INTRVL_ERROR)
End If ' if invalid interval
If (CSng(Abs(pdat.Longitude)) > 180#) Then
retval = retval Or ShL32F(1, S_LON_ERROR)
End If ' if invalid longitude
If (CSng(Abs(pdat.Latitude)) > 90#) Then
retval = retval Or ShL32F(1, S_LAT_ERROR)
End If ' if invalid latitude
End If ' if basic geometry
calculations required
If ((pdat.function And L_REFRAC) And (CSng(Abs(pdat.Temp)) > 100#))
Then
retval = retval Or ShL32F(1, S_TEMP_ERROR)
End If ' if invalid
temperatures
If ((pdat.function And L_REFRAC) And (pdat.press < 0#) Or
(pdat.press > 2000#)) Then
retval = retval Or ShL32F(1, S_PRESS_ERROR)
End If ' if invalid pressures
If ((pdat.function And L_TILT) And (CSng(Abs(pdat.Tilt)) > 180#))
Then
retval = retval Or ShL32F(1, S_TILT_ERROR)
End If ' if invalid tilt
If ((pdat.function And L_TILT) And (CSng(Abs(pdat.aspect)) > 360#))
Then
retval = retval Or ShL32F(1, S_ASPECT_ERROR)
End If ' if invalid aspect
If ((pdat.function And L_SBCF) And (pdat.sbwid < 1#) Or (pdat.sbwid
retval = retval Or ShL32F(1, S_SBWID_ERROR)
End If ' if invalid shadoWBRSds
width
If ((pdat.function And L_SBCF) And (pdat.sbrad < 1#) Or (pdat.sbrad
retval = retval Or ShL32F(1, S_SBRAD_ERROR)
End If ' if invalid shadoWBRSds
radians
If ((pdat.function And L_SBCF) And (CSng(Abs(pdat.sbsky)) > 1#))
Then
retval = retval Or ShL32F(1, S_SBSKY_ERROR)
End If ' if invalid shadoWBRSds
sky
validate = retval ' return error status
End Function
'*****************************************************************************
' Private Sub InitArray TWW porting construct
'
' loads day-of-month array with non-leap year and leap year days until
start
' of month
'
' Requires: nothing
'*****************************************************************************
Private Sub InitArray()
month_days(0, 2) = 31 ' Febuary non-leap
year
month_days(1, 2) = 31 ' Febuary leap year
month_days(0, 3) = 59 ' March non-leap
year
month_days(1, 3) = 60 ' March leap year
month_days(0, 4) = 90 ' April non-leap
year
month_days(1, 4) = 91 ' April leap year
month_days(0, 5) = 120 ' May non-leap
year
month_days(1, 5) = 121 ' May leap year
month_days(0, 6) = 151 ' June non-leap
year
month_days(1, 6) = 152 ' June leap year
month_days(0, 7) = 181 ' July non-leap
year
month_days(1, 7) = 182 ' July leap year
month_days(0, 8) = 212 ' August non-leap
year
month_days(1, 8) = 213 ' August leap year
month_days(0, 9) = 243 ' September non-leap
year
month_days(1, 9) = 244 ' September leap year
month_days(0, 10) = 273 ' October non-leap
year
month_days(1, 10) = 273 ' October leap year
month_days(0, 11) = 304 ' November non-leap
year
month_days(1, 11) = 305 ' November leap year
month_days(0, 12) = 334 ' December non-leap
year
month_days(1, 12) = 335 ' December leap year
array_init = True ' month array
initilized, we don't need to do it again
End Sub
'============================================================================
' Local Void function dom2doy
'
' Converts day-of-month to day-of-year
'
' Requires (from struct posdata parameter):
' year
' month
' day
'
' Returns (via the struct posdata parameter):
' year
' daynum
'-----------------------------------------------------------------------------
Private Sub dom2doy(pdat As posdata)
' cumulative number of days prior to beginning of month for leap year
pdat.daynum = pdat.Day + month_days(0, pdat.Month)
' Corrected by NREL (TW)
' If (((pdat.year Mod 4) = 0) And (((pdat.year Mod 100) <> 0) Or
((pdat.year Mod 400) = 0)) And (pdat.month > 2)) Then
If (((pdat.Year Mod 4) = 0) And (((pdat.Year Mod 100) <> 0) Or
((pdat.Year Mod 400) = 0))) Then
pdat.daynum = pdat.daynum + 1 ' adjust for leap year
End If ' if leap year
End Sub
'============================================================================
' Local void function doy2dom
'
' This function computes the month/day from the day number.
'
' Requires (from struct posdata parameter):
' Year and day number:
' year
' daynum
'
' Returns (via the struct posdata parameter):
' year
' month
' day
'-----------------------------------------------------------------------------
Private Sub doy2dom(pdat As posdata)
Dim Leap As Integer ' leap year switch
Dim imon As Integer ' Month (month_days)
array counte
If (((pdat.Year Mod 4) = 0) And (((pdat.Year Mod 100) <> 0) Or
((pdat.Year Mod 400) = 0))) Then
Leap = 1 ' set the leap year
switch
Else ' if not leap year
Leap = 0 ' set the non-leap year
switch
End If ' if leap year
imon = 12 ' default to Decamber
While pdat.daynum <= month_days(Leap, imon) ' loop to find number of
days prior to beginning of month
imon = imon - 1 ' decrement month
pointer
Wend ' loop to find number of
days prior to beginning of month
pdat.Month = imon ' set the month and day
of month
pdat.Day = pdat.daynum - month_days(Leap, imon)
End Sub
'============================================================================
' Local Void function geometry
'
' Does the underlying geometry for a given time and location
'-----------------------------------------------------------------------------
Private Sub geometry(pdat As posdata)
Dim bottom As Single ' denominator (bottom)
of the fraction
Dim c2 As Single ' cosine of d2
Dim cd As Single ' cosine of the day
angle or delination
Dim d2 As Single ' pdat.dayang times two
Dim Delta As Single ' difference between
current year and 1949
Dim s2 As Single ' sine of d2
Dim sd As Single ' sine of the day angle
Dim top As Single ' numerator (top) of the
fraction
Dim Leap As Integer ' (month_days) array
switch
' Day angle
' Iqbal, M. 1983. An Introduction to Solar Radiation. Academic
Press, NY., page 3
pdat.dayang = 360# * (pdat.daynum - 1) / 365#
' Earth radius vector * solar constant = solar energy Spencer, J. W.
1971. Fourier series representation of the
' position of the sun. Search 2 (5), page 172
sd = Sin(raddeg * pdat.dayang) ' get the sine of the
day angle
cd = Cos(raddeg * pdat.dayang) ' get the cosine of the
day angle or delination
d2 = 2# * pdat.dayang ' get the day angle
times two
c2 = Cos(raddeg * d2) ' get the cosine of d2
s2 = Sin(raddeg * d2) ' get the sine of d2
pdat.erv = 1.00011 + 0.034221 * cd + 0.00128 * sd
pdat.erv = pdat.erv + 0.000719 * c2 + 0.000077 * s2
' Universal Coordinated (Greenwich standard) time Michalsky, J.
1988. The Astronomical Almanac's algorithm for
' approximate solar position (1950-2050). Solar Energy 40 (3),
pp. 227-235.
pdat.utime = pdat.Hour * 3600# + pdat.Minute * 60# + pdat.second -
CSng(pdat.interval) / 2#
pdat.utime = pdat.utime / 3600# - pdat.TimeZone
' Julian Day minus 2,400,000 days (to eliminate roundoff errors)
Michalsky, J. 1988. The Astronomical Almanac's
' algorithm for approximate solar position (1950-2050). Solar
Energy 40 (3), pp. 227-235.
' No adjustment for century non-leap years since this function is
bounded by 1950 - 2050
Delta = pdat.Year - 1949 ' get the difference
between current year and 1949
Leap = CInt((Delta / 4#)) ' get (month_days) array
switch
pdat.julday = 32916.5 + Delta * 365# + Leap + pdat.daynum +
pdat.utime / 24#
' Time used in the calculation of ecliptic coordinates Noon 1 JAN
2000 = 2,400,000 + 51,545 days Julian Date
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar
' Energy 40 (3), pp. 227-235.
pdat.ectime = pdat.julday - 51545# ' get time of ecliptic
calculations
' Mean longitude
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar
' Energy 40 (3), pp. 227-235.
pdat.MnLong = 280.46 + 0.9856474 * pdat.ectime
' (dump the multiples of 360, so the answer is between 0 and 360) */
pdat.MnLong = pdat.MnLong - 360# * CInt((pdat.MnLong / 360#))
If (pdat.MnLong < 0#) Then '
pdat.MnLong = pdat.MnLong = 360# '
End If '
' Mean anomaly
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar
' Energy 40 (3), pp. 227-235.
pdat.MnAnom = 357.528 + 0.9856003 * pdat.ectime
' (dump the multiples of 360, so the answer is between 0 and 360)
pdat.MnAnom = pdat.MnAnom - 360# * CInt((pdat.MnAnom / 360#))
If (pdat.MnAnom < 0#) Then '
pdat.MnAnom = pdat.MnAnom + 360# '
End If '
' Ecliptic longitude */
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar
' Energy 40 (3), pp. 227-235.
pdat.EcLong = pdat.MnLong + 1.915 * Sin(pdat.MnAnom * raddeg) + 0.02
* Sin(2# * pdat.MnAnom * raddeg)
' (dump the multiples of 360, so the answer is between 0 and 360)
pdat.EcLong = pdat.EcLong - 360# * CInt(pdat.EcLong / 360#)
If (pdat.EcLong < 0#) Then
pdat.EcLong = pdat.EcLong + 360#
End If '
' Obliquity of the ecliptic
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3) pp. 227-235.
pdat.ecobli = 23.439 - 0.0000004 * pdat.ectime
' Declination
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3), pp. 227-235.
pdat.declin = degrad * Asin(Sin(pdat.ecobli * raddeg) *
Sin(pdat.EcLong * raddeg))
' Right ascension
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3), pp. 227-235.
top = Cos(raddeg * pdat.ecobli) * Sin(raddeg * pdat.EcLong)
bottom = Cos(raddeg * pdat.EcLong) '
pdat.rascen = degrad * Atan(top, bottom) '
' (make it a positive angle)
If (pdat.rascen < 0#) Then '
pdat.rascen = pdat.rascen + 360# '
End If '
' Greenwich mean sidereal time
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3), pp. 227-235.
pdat.Gmst = 6.697375 + 0.0657098242 * pdat.ectime + pdat.utime
' (dump the multiples of 24, so the answer is between 0 and 24)
pdat.Gmst = pdat.Gmst - 24# * CInt(pdat.Gmst / 24#)
If (pdat.Gmst < 0#) Then '
pdat.Gmst = pdat.Gmst + 24# '
End If '
' Local mean sidereal time
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3), pp. 227-235.
pdat.Lmst = pdat.Gmst * 15# + pdat.Longitude
' (dump the multiples of 360, so the answer is between 0 and 360)
pdat.Lmst = pdat.Lmst - 360# * CInt(pdat.Lmst / 360#)
If (pdat.Lmst < 0#) Then '
pdat.Lmst = pdat.Lmst + 360# '
End If '
' Hour angle
' Michalsky, J. 1988. The Astronomical Almanac's algorithm for
approximate solar position (1950-2050). Solar Energy
' 40 (3), pp. 227-235.
pdat.HrAng = pdat.Lmst - pdat.rascen '
' (force it between -180 and 180 degrees)
If (pdat.HrAng < -180#) Then '
pdat.HrAng = pdat.HrAng + 360# '
Else
If (pdat.HrAng > 180#) Then '
pdat.HrAng = pdat.HrAng - 360# '
End If '
End If '
End Sub
'=============================================================================
' Local Void function zen_no_ref
'
' ETR solar zenith angle
' Iqbal, M. 1983. An Introduction to Solar Radiation.
' Academic Press, NY., page 15
'-----------------------------------------------------------------------------
Private Sub zen_no_ref(pdat As posdata, tdat As trigdata)
Dim Cz As Single ' cosine of the solar
zenith angle
localtrig pdat, tdat '
Cz = tdat.sd * tdat.Sl + tdat.cd * tdat.cl * tdat.ch
' (watch out for the roundoff errors)
If CSng(Abs(Cz)) > 1# Then '
' If Abs(CZ) > 1# Then '
If (Cz >= 0#) Then '
Cz = 1# '
Else '
Cz = -1# '
End If '
End If '
pdat.zenetr = Acos(Cz) * degrad ''
' (limit the degrees below the horizon to 9 [+90 -> 99])
If (pdat.zenetr > 99#) Then '
pdat.zenetr = 99# '
End If '
pdat.elevetr = 90# - pdat.zenetr '
End Sub
'============================================================================
' Local Void function ssha
'
' Sunset hour angle, degrees
' Iqbal, M. 1983. An Introduction to Solar Radiation.
' Academic Press, NY., page 16
'-----------------------------------------------------------------------------
Private Sub ssha(pdat As posdata, tdat As trigdata)
Dim cssha As Single ' cosine of the sunset
hour angle
Dim cdcl As Single ' ( cd * cl )
localtrig pdat, tdat '
cdcl = tdat.cd * tdat.cl '
If CSng(Abs(cdcl)) >= 0.001 Then '
' If Abs(cdcl) >= 0.001 Then '
cssha = -tdat.Sl * tdat.sd / cdcl '
' This keeps the cosine from blowing on roundoff
If (cssha < -1#) Then '
pdat.ssha = 180# '
ElseIf (cssha > 1#) Then '
pdat.ssha = 0# '
Else '
pdat.ssha = degrad * Acos(cssha) '
End If '
ElseIf (((pdat.declin >= 0#) And (pdat.Latitude > 0#)) Or
((pdat.declin < 0#) And (pdat.Latitude < 0#))) Then
pdat.ssha = 180# '
Else '
pdat.ssha = 0# '
End If '
End Sub '
'=============================================================================
' Local Void function sbcf
'
' ShadoWBRSd correction factor
' Drummond, A. J. 1956. A contribution to absolute pyrheliometry.
' Q. J. R. Meteorol. Soc. 82, pp. 481-493
'-----------------------------------------------------------------------------
Private Sub sbcf(pdat As posdata, tdat As trigdata)
Dim P As Single ' used to compute sbcf
Dim t1 As Single ' used to compute sbcf
Dim t2 As Single ' used to compute sbcf
localtrig pdat, tdat '
P = 0.6366198 * pdat.sbwid / pdat.sbrad * tdat.cd ^ 3
t1 = tdat.Sl * tdat.sd * pdat.ssha * raddeg '
t2 = tdat.cl * tdat.cd * Sin(pdat.ssha * raddeg)
pdat.sbcf = pdat.sbsky + 1# / (1# - P * (t1 + t2))
End Sub
'=============================================================================
' Local Void function tst
'
' TST -> True Solar Time = local standard time + TSTfix, time
' in minutes from midnight.
' Iqbal, M. 1983. An Introduction to Solar Radiation.
' Academic Press, NY., page 13
'-----------------------------------------------------------------------------
Private Sub Tst(pdat As posdata)
pdat.Tst = (180# + pdat.HrAng) * 4# '
' add back half of the interval
pdat.tstfix = pdat.Tst - CSng(pdat.Hour) * 60# - pdat.Minute -
CSng(pdat.second) / 60# + CSng(pdat.interval) / 120#
While pdat.tstfix > 720# ' bound tstfix to this
day
pdat.tstfix = pdat.tstfix - 1440# '
Wend '
While pdat.tstfix < -720# '
pdat.tstfix = pdat.tstfix + 1440# '
Wend '
pdat.eqntim = pdat.tstfix + 60# * pdat.TimeZone - 4# *
pdat.Longitude
End Sub
'=============================================================================
' Local Void function srss
'
' Sunrise and sunset times (minutes from midnight)
'-----------------------------------------------------------------------------
Private Sub srss(pdat As posdata)
If (pdat.ssha <= 1#) Then '
pdat.sretr = 2999# '
pdat.ssetr = -2999# '
ElseIf (pdat.ssha >= 179#) Then '
pdat.sretr = -2999# '
pdat.ssetr = 2999# '
Else '
pdat.sretr = 720# - 4# * pdat.ssha - pdat.tstfix '
pdat.ssetr = 720# + 4# * pdat.ssha - pdat.tstfix '
End If '
End Sub
'=============================================================================
' Local Void function sazm
'
' Solar azimuth angle
' Iqbal, M. 1983. An Introduction to Solar Radiation.
' Academic Press, NY., page 15
'-----------------------------------------------------------------------------
Private Sub SAzm(pdat As posdata, tdat As trigdata)
Dim ca As Single ' cosine of the solar
azimuth angle
Dim ce As Single ' cosine of the solar
elevation
Dim cecl As Single ' ( ce * cl )
Dim se As Single ' sine of the solar
elevation
localtrig pdat, tdat '
ce = Cos(raddeg * pdat.elevetr) '
se = Sin(raddeg * pdat.elevetr) '
pdat.azim = 180# '
cecl = ce * tdat.cl '
If Abs(cecl) >= 0.001 Then '
ca = (se * tdat.Sl - tdat.sd) / cecl '
If (ca > 1#) Then '
ca = 1# '
ElseIf (ca < -1#) Then '
ca = -1# '
End If '
pdat.azim = 180# - Acos(ca) * degrad '
If (pdat.HrAng > 0) Then '
pdat.azim = 360# - pdat.azim '
End If '
End If '
End Sub
'=============================================================================
' Local Int function refrac
'
' Refraction correction, degrees
' Zimmerman, John C. 1981. Sun-pointing programs and their
accuracy.
' SAND81-0761, Experimental Systems Operation Division 4721,
' Sandia National Laboratories, Albuquerque, NM.
'-----------------------------------------------------------------------------
Private Sub Refrac(pdat As posdata)
Dim prestemp As Single ' temporary
pressure/temperature correction
Dim refcor As Single ' temporary refraction
correction
Dim tanelev As Single ' tangent of the solar
elevation angle
Dim Part1 As Single ' expression part 1
Dim Part2 As Single ' expression part 2
Dim Part3 As Single ' expression part 3
If (pdat.elevetr > 85#) Then ' If the sun is near
zenith, the algorithm bombs; refraction near 0
refcor = 0# '
Else ' Otherwise, we have
refraction
tanelev = Tan(raddeg * pdat.elevetr) '
If (pdat.elevetr >= 5#) Then '
refcor = 58.1 / tanelev - 0.07 / (tanelev ^ 3) + 0.000086 /
(tanelev ^ 5)
ElseIf (pdat.elevetr >= -0.575) Then '
' refcor = 1735# + pdat.elevetr * (-518.2 + pdat.elevetr * (103.4
+ pdat.elevetr * (-12.79 + pdat.elevetr * 0.711)))
Part1 = -12.79 + pdat.elevetr * 0.711 '
Part2 = 103.4 + pdat.elevetr * Part1 '
Part3 = -518.2 + pdat.elevetr * Part2 '
refcor = 1735# + pdat.elevetr * Part3 '
Else '
refcor = -20.774 / tanelev '
End If '
prestemp = (pdat.press * 283#) / (1013# * (273# + pdat.Temp))
refcor = refcor * prestemp / 3600# '
End If ' If the sun is near
zenith, the algorithm bombs; refraction near 0
pdat.elevref = pdat.elevetr + refcor ' Refracted solar
elevation angle
If (pdat.elevref < -9#) Then ' (limit the degrees
below the horizon to 9)
pdat.elevref = -9# '
End If ' (limit the degrees
below the horizon to 9)
pdat.zenref = 90# - pdat.elevref ' Refracted solar zenith
angle
pdat.coszen = Cos(raddeg * pdat.zenref) '
End Sub
'=============================================================================
' Local Void function amass
'
' Airmass
' Kasten, F. and Young, A. 1989. Revised optical air mass tables
and approximation formula. Applied Optics 28 (22),
' pp. 4735-4738
'-----------------------------------------------------------------------------
Private Sub amass(pdat As posdata)
If (pdat.zenref > 93#) Then '
pdat.amass = -1# '
pdat.ampress = -1# '
Else '
pdat.amass = 1# / (Cos(raddeg * pdat.zenref) + 0.50572 * (96.07995
- pdat.zenref) ^ -1.6364)
pdat.ampress = pdat.amass * pdat.press / 1013#
End If '
End Sub
'=============================================================================
' Local Void function prime
'
' Prime and Unprime
' Prime converts Kt to normalized Kt', etc.
' Unprime deconverts Kt' to Kt, etc. Perez, R., P. Ineichen, Seals,
R., & Zelenka, A. 1990. Making
' full use of the clearness index for parameterizing hourly
insolation conditions. Solar Energy 45 (2), pp. 111-114
'-----------------------------------------------------------------------------
Private Sub prime(pdat As posdata)
pdat.unprime = 1.031 * Exp(-1.4 / (0.9 + 9.4 / pdat.amass)) + 0.1
pdat.prime = 1# / pdat.unprime '
End Sub
'=============================================================================
' Local Void function etr
'
' Extraterrestrial (top-of-atmosphere) solar irradiance
'-----------------------------------------------------------------------------
Private Sub etr(pdat As posdata)
If (pdat.coszen > 0#) Then '
pdat.etrn = pdat.solcon * pdat.erv '
pdat.etr = pdat.etrn * pdat.coszen '
Else '
pdat.etrn = 0# '
pdat.etr = 0# '
End If '
End Sub
'=============================================================================
' Local Void function localtrig
'
' Does trig on internal variable used by several functions
'-----------------------------------------------------------------------------
Private Sub localtrig(pdat As posdata, tdat As trigdata)
' define masks to prevent calculation of uninitialized variables
Const SD_MASK As Long = (L_ZENETR Or L_SSHA Or S_SBCF Or
S_SOLAZM)
Const SL_MASK As Long = (L_ZENETR Or L_SSHA Or S_SBCF Or
S_SOLAZM)
Const CL_MASK As Long = (L_ZENETR Or L_SSHA Or S_SBCF Or
S_SOLAZM)
Const CD_MASK As Long = (L_ZENETR Or L_SSHA Or S_SBCF)
Const CH_MASK As Long = (L_ZENETR)
If (tdat.sd < -900#) Then ' sd was initialized
-999 as flag
tdat.sd = 1# ' reflag as having
completed calculations
If (pdat.function Or CD_MASK) Then tdat.cd = Cos(raddeg *
pdat.declin)
If (pdat.function Or CH_MASK) Then tdat.ch = Cos(raddeg *
pdat.HrAng)
If (pdat.function Or CL_MASK) Then tdat.cl = Cos(raddeg *
pdat.Latitude)
If (pdat.function Or SD_MASK) Then tdat.sd = Sin(raddeg *
pdat.declin)
If (pdat.function Or SL_MASK) Then tdat.Sl = Sin(raddeg *
pdat.Latitude)
End If '
End Sub
'=============================================================================
' Local Void function tilt
'
' ETR on a tilted surface
'-----------------------------------------------------------------------------
Private Sub Tilt(pdat As posdata)
Dim ca As Single ' cosine of the solar
azimuth angle
Dim cp As Single ' cosine of the panel
aspect
Dim ct As Single ' cosine of the panel
tilt
Dim sa As Single ' sine of the solar
azimuth angle
Dim sp As Single ' sine of the panel
aspect
Dim st As Single ' sine of the panel tilt
Dim sz As Single ' sine of the refraction
corrected solar zenith angle
' Cosine of the angle between the sun and a tipped flat surface,
' useful for calculating solar energy on tilted surfaces
ca = Cos(raddeg * pdat.azim) '
cp = Cos(raddeg * pdat.aspect) '
ct = Cos(raddeg * pdat.Tilt) '
sa = Sin(raddeg * pdat.azim) '
sp = Sin(raddeg * pdat.aspect) '
st = Sin(raddeg * pdat.Tilt) '
sz = Sin(raddeg * pdat.zenref) '
pdat.CosInc = pdat.coszen * ct + sz * st * (ca * cp + sa * sp)
If (pdat.CosInc > 0#) Then '
pdat.etrtilt = pdat.etrn * pdat.CosInc '
Else '
pdat.etrtilt = 0# '
End If
End Sub
'=============================================================================
' Void function S_decode
'
' This function decodes the error codes from S_solpos return value
'
' Requires the long integer return value from S_solpos
'
' Returns descriptive text to stderr
'-----------------------------------------------------------------------------
Public Sub S_decode(Code As Long, pdat As posdata)
If Code = 0 Then Exit Sub ' if no errors, return
(added by Tommy)
If (Code And (ShL32F(1, S_YEAR_ERROR))) Then '
MsgBox "S_decode ==> Please fix the year: " & pdat.Year & "
[1950-2050]"
End If '
If (Code And (ShL32F(1, S_MONTH_ERROR))) Then '
MsgBox "S_decode ==> Please fix the month: " & pdat.Month
End If '
If (Code And (ShL32F(1, S_DAY_ERROR))) Then '
MsgBox "S_decode ==> Please fix the month: " & pdat.Day
End If '
If (Code And (ShL32F(1, S_DOY_ERROR))) Then '
MsgBox "S_decode ==> Please fix the day-of-year: " & pdat.daynum
End If '
If (Code And (ShL32F(1, S_HOUR_ERROR))) Then '
MsgBox "S_decode ==> Please fix the hour: " & pdat.Hour
End If '
If (Code And (ShL32F(1, S_MINUTE_ERROR))) Then '
MsgBox "S_decode ==> Please fix the minute: " & pdat.Minute
End If '
If (Code And (ShL32F(1, S_SECOND_ERROR))) Then '
MsgBox "S_decode ==> Please fix the second: " & pdat.second
End If '
If (Code And (ShL32F(1, S_TZONE_ERROR))) Then '
MsgBox "S_decode ==> Please fix the time zone: " & pdat.TimeZone
End If '
If (Code And (ShL32F(1, S_INTRVL_ERROR))) Then '
MsgBox "S_decode ==> Please fix the interval: " & pdat.interval
End If '
If (Code And (ShL32F(1, S_LAT_ERROR))) Then '
MsgBox "S_decode ==> Please fix the latitude: " & pdat.Latitude
End If '
If (Code And (ShL32F(1, S_LON_ERROR))) Then '
MsgBox "S_decode ==> Please fix the longitude: " & pdat.Longitude
End If '
If (Code And (ShL32F(1, S_TEMP_ERROR))) Then '
MsgBox "S_decode ==> Please fix the temperature: " & pdat.Temp
End If '
If (Code And (ShL32F(1, S_PRESS_ERROR))) Then '
MsgBox "S_decode ==> Please fix the pressure: " & pdat.press
End If '
If (Code And (ShL32F(1, S_TILT_ERROR))) Then '
MsgBox "S_decode ==> Please fix the tilt: " & pdat.Tilt
End If '
If (Code And (ShL32F(1, S_ASPECT_ERROR))) Then '
MsgBox "S_decode ==> Please fix the aspect: " & pdat.aspect
End If '
If (Code And (ShL32F(1, S_SBWID_ERROR))) Then '
MsgBox "S_decode ==> Please fix the shadoWBRSd width: " &
pdat.sbwid
End If '
If (Code And (ShL32F(1, S_SBRAD_ERROR))) Then '
MsgBox "S_decode ==> Please fix the shadoWBRSd radius: " &
pdat.sbrad
End If '
If (Code And (ShL32F(1, S_SBSKY_ERROR))) Then '
MsgBox "S_decode ==> Please fix the shadoWBRSd sky factor: " &
pdat.sbsky
End If '
End Sub
*****************************************************************************
Option Compare Database
Option Explicit
Public Const Pi = 3.14159265 ' define pi
Public Function Radians(Degree As Double) As Double
' convert degrees to radians
Radians = Degree * Pi / 180 ' calculate and
return degrees converted to radians
End Function
Public Function Degrees(Radian As Double) As Double
' convert radians to degrees
Degrees = Radian * 180 / Pi ' calculate and
return radians converted to degrees
End Function
Public Function EquationOfTime(DayNumber As Double) As Double
' Equation of Time using day number
Dim R As Double ' radians
R = Radians((360 * (DayNumber - 1)) / 365.242) ' load radians
EquationOfTime = 0.258 * Cos(R) - 7.416 * Sin(R) - 3.648 * Cos(2 *
R) - 9.228 * Sin(2 * R)
End Function
Public Function LongitudeCorrection(TimeZone As Double, Longitude As
Double) As Double
' Longitude Correction using time zone and longitude
LongitudeCorrection = (15 * TimeZone - Longitude) / 15 ' calculate and
return solar time
End Function
Public Function LocalTime(ST As Double, DST As Double, EOT As Double,
LC As Double) As Double
' Local Time using solar time, daylight saving time, equation of time,
and longitude correction
LocalTime = ST - (EOT / 60) + LC + DST ' calculate and
return local time
End Function
Public Function SolarTime(LTD As Double, DST As Double, EOT As Double,
LC As Double) As Double
' Solar Time using local time decmal, daylight saving time, equation
of time, and longitude correction
SolarTime = LTD + (EOT / 60) - LC - DST ' calculate and
return solar time
End Function
Public Function DecimalToHMS(X As Double) As String
' Return time in 24Hr. string formatted 19:29:12
Dim Hrs As Integer ' hours
Dim Mins As Integer ' minutes
Dim Secs As Integer ' seconds
Hrs = Int(X) ' calculate hours
Mins = Int(60 * (X - Hrs)) ' calculate
minutes
Secs = Int(3600 * (X - Hrs - Mins / 60)) ' calculate
seconds
DecimalToHMS = Format(Hrs, "00") & ":" & Format(Mins, "00") & ":" &
Format(Secs, "00")
End Function
' arc sine
' error if value is outside the range [-1,1]
Function ASin(value As Double) As Double
If Abs(value) <> 1 Then
ASin = Atn(value / Sqr(1 - value * value))
Else
ASin = 1.5707963267949 * Sgn(value)
End If
End Function
' arc cosine
' error if NUMBER is outside the range [-1,1]
Function ACos(ByVal number As Double) As Double
If Abs(number) <> 1 Then
ACos = 1.5707963267949 - Atn(number / Sqr(1 - number * number))
ElseIf number = -1 Then
ACos = 3.14159265358979
End If
'elseif number=1 --> Acos=0 (implicit)
End Function
' arc cotangent
' error if NUMBER is zero
Function ACot(value As Double) As Double
ACot = Atn(1 / value)
End Function
' arc secant
' error if value is inside the range [-1,1]
Function ASec(value As Double) As Double
' NOTE: the following lines can be replaced by a single call
' ASec = ACos(1 / value)
If Abs(value) <> 1 Then
ASec = 1.5707963267949 - Atn((1 / value) / Sqr(1 - 1 / (value *
value)))
Else
ASec = 3.14159265358979 * Sgn(value)
End If
End Function
' arc cosecant
' error if value is inside the range [-1,1]
Function ACsc(value As Double) As Double
' NOTE: the following lines can be replaced by a single call
' ACsc = ASin(1 / value)
If Abs(value) <> 1 Then
ACsc = Atn((1 / value) / Sqr(1 - 1 / (value * value)))
Else
ACsc = 1.5707963267949 * Sgn(value)
End If
End Function
]
