3
$\begingroup$

Good Afternoon, I have calculated the Incident beam radiation for a given locations local sunlight hours.

Rise = Sunrise[GeoPosition[{57.7053, -3.33917}], 
   DateRange[DateObject[{2019, 1, 1}], DateObject[{2019, 12, 31}]]];
Sunrises = DateString[#, {"Hour", "Minute"}] & /@ Rise["Values"]
Sets = Sunset[GeoPosition[{57.7053, -3.33917}], 
   DateRange[DateObject[{2019, 1, 1}], DateObject[{2019, 12, 31}]]];
Sunsets = DateString[#, {"Hour", "Minute"}] & /@ Sets["Values"]; 

Sunrises[[1]]
Sunsets[[1]]

n = Table[0 + i, {i, 1, 1}];(*Defines what day of the year is being used*)

LSH = Table[i + 0.01, {i, 8.93, 15.616667, 0.01}];
(*Local Solar Hour,this shows the daylight hours of the location*)
H = 360/24*(LSH - 12);(*Hour Angle*)
L = 57.7053;(*Latitude*)

\[Delta] = 
 23.45*Sin[360/365*(n - 81) Degree](*Solar Declination of the sun*);

\[Beta] = (Cos[L Degree]*Cos[First[\[Delta]] Degree ]*
 Cos[H Degree ]) + ((Sin[L Degree]*Sin[First[\[Delta]] Degree]));

Arc\[Beta] = ArcSin[\[Beta] ]*180/Pi;(*Shows the altitude angle of the sun*)
A = 1160 + 75 Sin[360/365*(n - 275) Degree](*Extraterrestial Flux*);

m = 1/(Sin[
 Arc\[Beta] Degree])(*Air Mass Ratio for every hour of the day*);
k = 0.174 + 0.035*Sin[360/365*(n - 100) Degree];(*Optical Depth*)

IB = First[A]*E^(-First[k]*m); (*Direct Beam radiation Wm^-2*)
test = Transpose[{LSH, IB}];

ListLinePlot[test, PlotRange -> {{8.5, 16}, {0, 600}}]

This code shows the sunrise and sunset times of a given location for a specified day of the year It then calculates the angles of the sun that are required to calculate the Incident Beam radiation; or how much sunlight is incident on a location as trh sun traces a path through the sky. It then transposes the LocalSolarHour and Incident Radiation datasets so that it can be plotted to show how much radiation is incident over the the daylight hours of the location for a given day. From the ListLinePlot, it should show solar radiation for the daylight hours of the location. Where sunrise is 0856 and sunset is 1537 (Considering time of year) However the graph shows a sudden increase before the given sunset time where the graph should have a bell shape resemblence between 0856 and 1537(8.933333 and 15.616667 in numerical format) I've been going over this code for days trying to solve this so any help would be appreciated.

Luke

Solar Radiation vs Daylight Hours

$\endgroup$
  • $\begingroup$ Would you be able to add in some more detail what your code is meant to do? It's a little unclear to me what the variables are and so on. I'm curious what the algorithm is doing. Thanks! $\endgroup$ – Carl Lange Sep 21 '19 at 15:11
  • 1
    $\begingroup$ Of course, sorry for being sparse with details! hope this is clearer! $\endgroup$ – Luke4737 Sep 21 '19 at 15:19
  • $\begingroup$ I would guess this is because of β becoming negative when the sun goes down. If you clip β to 0 I think this will work fine (eg by using Ramp[β]) $\endgroup$ – Carl Lange Sep 21 '19 at 15:21
2
$\begingroup$

There it will be better to put a limit on m, then there will be no division by 0 and you can calculate the solar flux at any time of the day

Rise = Sunrise[GeoPosition[{57.7053, -3.33917}], 
   DateRange[DateObject[{2019, 1, 1}], DateObject[{2019, 12, 31}]]];
Sunrises = DateString[#, {"Hour", "Minute"}] & /@ Rise["Values"];
Sets = Sunset[GeoPosition[{57.7053, -3.33917}], 
   DateRange[DateObject[{2019, 1, 1}], DateObject[{2019, 12, 31}]]];
Sunsets = DateString[#, {"Hour", "Minute"}] & /@ Sets["Values"];


n = Table[
  0 + i, {i, 1, 
   1}];(*Defines what day of the year is being used*)LSH = 
 Table[i, {i, 8.93, 15.616667, 0.01}];
(*Local Solar Hour,this shows the daylight hours of the location*)
H = 360/24*(LSH - 12);(*Hour Angle*)L = 57.7053;(*Latitude*)\[Delta] =
  23.45*Sin[360/365*(n - 81) Degree](*Solar Declination of the sun*);

\[Beta] = (Cos[L Degree]*Cos[First[\[Delta]] Degree]*
     Cos[H Degree]) + ((Sin[L Degree]*Sin[First[\[Delta]] Degree]));

(*Arc\[Beta]=ArcSin[\[Beta]]*180/Pi;*)Arc\[Beta] = 
 ArcSin[Abs[\[Beta]]]*180/
   Pi;(*Shows the altitude angle of the sun*)A = 
 1160 + 75 Sin[360/365*(n - 275) Degree](*Extraterrestial Flux*);

m = Table[
   If[0 < 1/(Sin[ArcSin[\[Beta][[i]]]*180/Pi Degree]) < 100, 
    1/(Sin[ArcSin[\[Beta][[i]]]*180/Pi Degree]), 100], {i, 1, 
    Length[\[Beta]]}](*Air Mass Ratio for every hour of the day*);
k = 0.174 + 0.035*Sin[360/365*(n - 100) Degree];(*Optical Depth*)IB = 
 First[A]*E^(-First[k]*m);(*Direct Beam radiation Wm^-2*)test = 
 Transpose[{LSH, IB}];

The solar flux from sunrise to sunset and during the day

{ListLinePlot[test, PlotRange -> {{8.5, 16}, {0, 600}}], 
 ListLinePlot[test, PlotRange -> {{0, 24}, {0, 600}}]}

Figure 1

$\endgroup$
7
$\begingroup$

This is because of the value of β becoming negative. If we replace the Arcβ definition with the following:

Arc\[Beta] = ArcSin[Ramp[\[Beta]] ]*180/Pi;

your code appears to work ok (although it complains a little because now there are divisions by zero. It's unclear to me how best to handle these given your algorithm, so I would simply Quiet them - however, you may have a better alternative.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.