Estimation of clear-sky solar radiation

Contributions to this software library are always welcome. Please ensure that you post program listings rather than .raw files. They give a reasonable idea of what your program does without having to load them into a DM42 and you can also include comments in your code. Check out the following link for a decoder/encoder: https://technical.swissmicros.com/decoders/dm42/

You can then copy/paste the listing and post it in "code" tags.
Post Reply
User avatar
Patrizia_Favaron
Posts: 12
Joined: Thu Mar 11, 2021 9:44 am
Location: Lomagna

Estimation of clear-sky solar radiation

Post by Patrizia_Favaron »

Moderator note: 42 code placed in "code" tags for compactness

1.Background

Incoming solar radiation \(R_{g}\) at a point on the Earth surface is a very important environmental parameter: it is indeed the far dominant source of energy driving Earth natural processes, life among them.

Meanwhile, unfortunately, \(R_{g}\) is quite difficult to measure precisely. The sensor used to measure inbound solar radiation, the pyranometer, is prone to misalignment (should be placed horizontally) and other maintenance mishaps, and its components descriptive parameters may drift with time. To make things worse, it is really difficult to calibrate: many operators prefer to perform field-comparison campaigns rather than testing outcome against some "reference" artificial signal: as far as I know, we humans have not yet reproduced the Sun behaviour exactly.

It is then important to be able to estimate the clear-sky inbound solar radiation, \(R_{so}\): would a pyranometer reading \(R_{g}\) exceed \(R_{so}\) significantly, we would be in serious troubles...

If you, like me, are curious of what happens to inbound solar radiation, a reading of the current edition of IPCC technical report might be an interesting starter. Incidentally, this is my current research interest, and if you like I may give you more references to get deeper.

But for the purposes of estimation I've chosen to stay simple, and use a method amenable to a simple implementation on the DM42. "Simple" here means using a minimum of parameters to enter, so that would you be on field you could check sensor readings with relative ease. More specifically, I selected the approximate method described in

R.G. Allen, I.A. Walter, R.L. Elliott, T.A. Howell, D. Itenfisu, M.E. Jensen, R.L. Snyder (eds.), The ASCE Standardized Reference Evapotranspiration Equation, ASCE, 2005, ISBN 0-7844-0805-X

The "ASCE model" is in particular non-refractive (it does not take into account the effect of refraction due to decreasing air density with altitude). Although approximate, it is accurate enough for engineering applications and for gross pyranometer malfunction/mis-calibration check.

2.The ASCE clear-sky solar radiation model

In this section I'll list the model's formulae, using as much as possible the original ASCE notation. Of the many relationships present in the original reference you will find here only the ones relevant to estimating the clear-sky solar radiation.

We start with the formula for the clear-sky radiation itself:

\(R_{so} = \left( 0.75 + 2\cdot 10^{-5} z \right) \cdot R_{a}\)

where \(z\) is the height above mean sea level (m) and \(R_{a}\) the so-called "extraterrestrial radiation", that is the solar radiation available at the top of the atmosphere (W/m^2). The extraterrestrial radiation is in turn estimated as

\( R_{a} = \frac{\pi}{12} G_{s} d_{r} \left[ (\omega_{2}-\omega_{1}) \cdot \sin \phi \cdot \sin \delta + \left( \sin \omega_{2} - \sin \omega_{1} \right) \cdot \cos \phi \cdot \cos \delta\right] \)

with \(G_{s}\) the solar "constant" (1367 W/m^2). \(d_r\) is the squared inverse relative distance factor for the Earth-Sun system (unitless), \(\delta\) is the solar declination (radians), \(\phi\) is latitude (radians), \(\omega_{1}\) is the solar azimuth at beginning of averaging period (radians), and \(\omega_{2}\) the solar azimuth at end of the averaging period (radians).

The two solar azimuths can be expressed preliminarily as

\(\omega_{1} = \omega - \frac{\pi \Delta t}{24}\)

and

\(\omega_{2} = \omega + \frac{\pi \Delta t}{24}\)

respectively, with \(\Delta t\) the averaging time. The mid-period azimuth \(\omega\) is computed through the formula

\(\omega = \frac{\pi}{12} \left[ t + 0.06667 \cdot \left( L_{z} - L_{m} \right) + S_{c} - 12 \right]\)

with \(t\) the mid-time-interval instant in decimal hours, \(L_{z}\) the longitude of the central meridian of the time zone (degrees, positive towards West), \(L_{m}\) is the longitude of the location where the estimation is made (degrees, positive towards West), and \(S_{c}\) the seasonal correction for solar time (decimal hours).

Of course the extraterrestrial (and clear-sky) solar radiation is non-zero only on daytime, so the two preliminary evaluations of \(\omega{1}\) and \(\omega_{2}\) should be corrected by limiting the actual range to sunrise and sunset, that is, when sunlight is actually present. The sunset azimuth is given by

\(\omega_{s} = \cos^{-1} \left( -\tan \phi \tan \delta \right)\)

Of course the sunrise is \(-\omega_{s}\). Given the sunrise and sunset angles, we may then calculate the final \(\omega_{1}\) and \(\omega_{2}\):

\(\omega_{1} = \max \left( \omega_{1}, -\omega_{s} \right)\)
\(\omega_{2} = \max \left( \omega_{2}, -\omega_{s} \right)\)
\(\omega_{1} = \min \left( \omega_{1}, \omega_{s} \right)\)
\(\omega_{2} = \min \left( \omega_{2}, \omega_{s} \right)\)
\(\omega_{1} = \min \left( \omega_{1}, \omega_{2} \right)\)

In order to completely specify the model, we have to assign values to \(S_{c}\), \(d_{r}\) and \(\delta\). The seasonal correction for the solar time is

\(S_{c} = 0.1645 \sin(2b) - 0.1255 \cos b - 0.025 \sin b\)

with

\(b = 2 \pi \frac{J - 81}{364}\)

where \(J\) is the number of day in year.

The inverse squared distance factor is

\(d_{r} = 1 + 0.033 \cos \left( \frac{2 \pi}{365} J \right)\)

The solar declination is

\(\delta = 0.409 \sin \left( \frac{2 \pi}{365} J - 1.39 \right)\)

All the three last formulae depend on the number of day within year:

\(J = D - 32 \left\lfloor 275 \frac{M}{9} \right\rfloor + 2 \left\lfloor \frac{3}{M+1} \right\rfloor \left\lfloor \frac{M}{100} - \frac{Y \mod 4}{4} + 0.975 \right\rfloor \)

This concludes the ASCE model presentation. As a comment, we may notice all relations are algebraic: overall the formulae may seem "complicate", but the model itself is "mathematically simple". This is aligned with the desire I mentioned.

3.Implementation

3.1.Foreword

Any implementation is to some extent subjective, and mine is not an exception. Even a simple machine like the DM42 permits many coding paths. The one I decided to follow is through use of registers, namely R00 to R19, to hold input data and intermediate results.

Additionally, I've chosen to divide the code in two functional parts. The first, using the same heterogeneous units of the original ASCE model, contains all the formulae (one per program), written in a form which is basically compatible with the DM41X. The second acts as an interface between input data, in forms which are uniform and customary, and the final result.

3.2.Register use

Here is the (inelegant, I admit) way I used the registers:

R00 - J (day in the year - computed)
R01 - Y (year, input)
R02 - M (month, input)
R03 - D (day, input)
R04 - \(d_{r}\) (inverse squared distance factor, computed)
R05 - \(\delta\) (solar declination, computed)
R06 - \(S_{c}\) (Hour correction term, computed)
R07 - \(t\) (Decimal hour, input)
R08 - \(L_{z}\) (Time zone central meridian longitude, degrees, positive towards West, input)
R09 - \(L_{m}\) (Measurement site longitude, degrees, positive westwards, input)
R10 - \(\omega\) (Solar azimuth, computed)
R11 - \(\phi\) (Measurement site latitude, radians, positive northwards, input)
R12 - \(\omega_{s}\) (Azimuth at sunset, computed)
R13 - \(\Delta_{t}\) (Time step, decimal hours, input)
R14 - \(\omega_{1}\) (Initial step azimuth, computed)
R15 - \(\omega_{2}\) (Final step azimuth, computed)
R16 - \(R_{a}\) (Extraterrestrial radiation, computed)
R17 - \(z\) (Measurement site altitude above mean sea level, meters, input)
R18 - \(R_{so}\) (Clear-sky solar radiation estimate, W/m^2, computed)
R19 - Temporary storage.

3.3.Computational programs

Here are the programs constituting the computational part, in HP-42S form (revised a little bit for Windows keyboard...). They should be compatible with DM41X.

In the descriptions, I'll roughly present the listings in reverse order respect to the formulae.

The first program computes the day in year:

Code: Select all

01 LBL "JVAL"
02 RCL 03
03 32
04 -
05 RCL 02
06 9
07 /
08 275
09 x
10 IP
11 +
12 3
13 RCL 02
14 1
15 +
16 /
17 IP
18 2
19 x
20 +
21 RCL 02
22 100
23 /
24 RCL 01
25 4
26 MOD
27 4
28 /
29 -
30 0.975
31 +
32 IP
33 +
34 STO 00
35 RTN
The next program computes the solar declination, seasonal time correction factor and the distance factor. All depend on the day in year, J:

Code: Select all

01 LBL "DECL"
02 RAD
03 RCL 00
04 2
05 x
06 PI
07 x
08 365
09 /
10 COS
11 0.033
12 x
13 1
14 +
15 STO 04
16 RCL 00
17 2
18 x
19 PI
20 x
21 365
22 /
23 1.39
24 -
25 SIN
26 0.409
27 x
28 STO 05
29 RCL 00
30 81
31 -
32 2
33 x
34 PI
35 x
36 364
37 /
38 STO 06
39 2
40 x
41 SIN
42 0.1645
43 x
44 RCL 06
45 COS
46 0.1255
47 x
48 -
49 RCL 06
50 SIN
51 0.025
52 x
53 -
54 STO 06
55 RTN
Then we have the program for computing the current solar azimut:

Code: Select all

01 LBL "OMEGA"
02 RCL 08
03 RCL 09
04 -
05 0.06667
06 x
07 RCL 07
08 +
09 RCL 06
10 +
11 12
12 -
13 PI
14 x
15 12
16 /
17 STO 10
18 RTN
Here we compute the sunset azimuth:

Code: Select all

01 LBL "OMS"
02 RCL 11
03 TAN
04 RCL 05
05 TAN
06 x
07 +/-
08 ACOS
09 STO 12
10 RTN
The sunset azimuth allows computing and delimiting the initial and final averaging time azimuths:

Code: Select all

01 LBL "OM12"
02 RCL 13
03 PI
04 x
05 24
06 /
07 ENTER
08 +/-
09 RCL 10
10 +
11 STO 14
12 X<>Y
13 RCL 10
14 +
15 STO 15
16 RCL 14
17 RCL 12
18 +/-
19 XEQ 00
20 STO 14
21 RCL 15
22 RCL 12
23 +/-
24 XEQ 00
25 STO 15
26 RCL 14
27 RCL 12
28 XEQ 01
29 STO 14
30 RCL 15
31 RCL 12
32 XEQ 01
33 STO 15
34 RCL 14
35 RCL 15
36 XEQ 01
37 STO 14
38 RTN
39 LBL 00      ; Min(X,Y)
40 X<Y?
41 X<>Y
42 RTN
43 LBL 01      ; Max(X,Y)
44 X>Y?
45 X<>Y
46 RTN
We have all the elements to compute the extraterrestrial radiation:

Code: Select all

01 LBL "RA"
02 RCL 15
03 SIN
04 RCL 14
05 SIN
06 -
07 RCL 11
08 COS
09 x
10 RCL 05
11 COS
12 x
13 RCL 15
14 RCL 14
15 -
16 RCL 11
17 SIN
18 x
19 RCL 05
20 SIN
21 x
22 +
23 RCL 04
24 x
25 12
26 x
27 PI
28 /
29 1367
30 x
31 RCL 13
32 /
33 STO 16
34 RTN
Notice lines 31 and 32: the result is divided by \(\Delta t\), in order to rescale it to a "hourly" form. The original formula does not contain this term, and calculates an energy flux whose magnitude does, for example, halves if the averaging interval is also halved. In case the energy flux is what you need, then you may multiply the result back by \(\Delta t\).

The final program in the computational section calculates the actual estimate of inbound solar radiation.

Code: Select all

01 LBL "RSOL"
02 RCL 17
03 2E-5
04 x
05 0.75
06 +
07 RCL 16
08 x
09 STO 18
10 RTN
3.4.User interface

By paying much attention, one might use the computational programs in section 3.3 directly, as they are. But fortunately there is another, more humane approach: building a user interface which converts between the user world and actual values. The user interface may well be not only performing conversions: it might attempt making the user's lives simpler.

For instance, parameter input may be made uniform, by adopting a unique measurement unit and form, of common use; as an European, I would also like longitudes to be positive towards East; and as a reasonably lazy person I would also be grateful if date is inserted in YYYYMMDD form, instead of inputting its three parts separately.

Access and storage of all input parameters would benefit greatly, if it is performed using a programmable menu. This is of noteworthy practical importance as it allows to store the parameters in the order you desire, and if needed to change one of them only. This would be impossible using a conventional interactive program, but unfortunately restricts the system type to HP-42S and DM42. This may be quite unhappy, and yet from time to time we have to make some choices in our life. :)

The program, listed here, is commented to allow following the unfolding of conversions.

Code: Select all

001 LBL "SOLAR"
002 LBL A
003 1                    ; Default value for [latex]\Delta t[/latex] = 1 hour
004 STO 13
005 CLST
006 "DATE"
007 KEY 1 XEQ 01  ; Enter date in YYYYMMDD form
008 "TIME"
009 KEY 2 XEQ 02  ; Enter time as decimal hour
010 "DT"
011 KEY 3 XEQ 03  ; Enter averaging time as decimal hour (default = 1 hour)
012 "ZONE"
013 KEY 4 XEQ 04  ; Enter time zone (1 for Paris and Rome) and deduce central longitude from it
014 "LON"
015 KEY 5 XEQ 05  ; Enter longitude in ° ' " sexagesimal form
016 "LAT"
017 KEY 6 XEQ 06  ; Enter latitude in ° ' " sexagesimal form
018 KEY 7 GTO B   ; Key "up"
019 KEY 8 GTO B   ; Key "down"
020 KEY 9 GTO 99  ; Key "exit"
021 MENU              ; Show menu
022 LBL 20            ; Manage successive keypress of R/S
023 STOP              ; by "doing nothing"
024 GTO 20

025 LBL B
026 "Z"
027 KEY 1 XEQ 07  ; Enter altitude in m above mean sea level
028 "RUN"
029 KEY 2 XEQ 08  ; Perform calculations in the proper sequence
030 KEY 7 GTO A   ; Key "up"
031 KEY 8 GTO A   ; Key "down"
032 LBL 21            ; Manage successive keypresses of R/S
033 STOP              ; by doing "nothing"
034 GTO 21

035 LBL 99            ; Perform a graceful exit by removing the menu
036 CLMENU
037 EXITALL
038 RTN

039 LBL 01            ; Split date from YYYYMMDD into components
040 STO 19
041 100
042 /
043 FP
044 100
045 x
046 STO 03             ; Day
047 RCL 19
048 100
049 /
050 IP
051 100
052 /
053 FP
054 100
055 x
056 STO 02             ; Month
057 RCL 19
058 1E4
059 /
060 IP
061 STO 01             ; Year
062 CLST
063 RTN

064 LBL 02             ; Conversion of time at start of interval to time at center interval
065 >HR
066 RCL 13
067 2
068 /
069 +
070 STO 07
071 CLST
072 RTN

073 LBL 03             ; Change the value of [latex]\Delta t[/latex], and update time at center interval accordingly
074 RCL 07
075 RCL 13
076 2
077 /
078 -
079 X<>Y
080 >HR
081 STO 13
082 2
083 /
084 +
085 STO 07
086 CLST
087 RTN

088 LBL 04             ; Convert to time zone to west-bound longitude
089 -15
090 x
091 STO 08
092 CLST
093 RTN

094 LBL 05             ; Convert east-bound sexagesimal longitude to west-bound decimal
095 +/-
096 >HR
097 STO 09
098 CLST
099 RTN

100 LBL 06             ; Convert latitude from sexagesimal degrees to radians
101 >HR
102 >RAD
103 STO 11
104 CLST
105 RTN

106 LBL 07             ; Store altitude as it is
107 STO 17
108 CLST
109 RTN

110 LBL 08             ; Execute calculations in proper order, assuming all input data have been assigned
111 XEQ "JVAL"
112 XEQ "DECL"
113 XEQ "OMEGA"
114 XEQ "OMS"
115 XEQ "OM12"
116 XEQ "RA"
117 XEQ "RSOL"
118 CLST
119 RCL 18             ; The actual clear-sky solar radiation value
120 RTN

4.What am I using this program for

Sure all the formulae I wanted to code would have found a much more compact form in a high level language like Fortran or Julia (and they will soon). But the very high precision of the DM42 allows comparing their results to conventional programs, with the additional benefit of a fixed and well defined instruction execution order. Would any significant difference be found, then a delicate point in implementation would have been found, possibly due to a numerical instability: something to work on in the "final" high-level implementations.

Another of my uses is on-field, to check pyranometers during met station maintenance. Being it possible to perform calculations on a calculator is for me an advantage: the calculator fits in my purse, but a computer would decidedly not.

Of course many other uses are possible. Including e.g. to understand the behaviour of some astronomical quantities like the solar declination.

5.What next?

I've done my best to provide a correct and accurate program, but can not guarantee some error is still present in my code. I'll fix them and post the results as comments to this short report. If you find one, of course, please tell me about it: I'll be glad to correct.

For now, sorry so sloppy.

Patrizia
User avatar
akaTB
Posts: 794
Joined: Tue May 02, 2017 1:56 pm
Location: Milan, Italy

Re: Estimation of clear-sky solar radiation

Post by akaTB »

so sloppy?
I'd call this high quality documentation.

Hats off, Patrizia. Thank you.
Greetings,
    Massimo
ajcaton
-+×÷ left is right and right is wrong :twisted: Casted in gold
User avatar
Patrizia_Favaron
Posts: 12
Joined: Thu Mar 11, 2021 9:44 am
Location: Lomagna

Re: Estimation of clear-sky solar radiation

Post by Patrizia_Favaron »

Thanks, too generous!
Post Reply