models/johnpye/sunpos.a4c

REQUIRE "atoms.a4l";

(*
	Using the equations given in Duffie & Beckman (1980)
	the following model allows calculation of the angle of incidence of
	sunlight onto an inclined surface

	The equation of time has been changed to the Duffie & Beckman 2nd edition
	form, cited as coming from Iqbal 1983. The equation for 'B' has also
	been changed to the 2nd Ed D&B form.

	For the longitude correction we use the spherical geometry convention
	of making east positive. D&B are americocentric, obviously :-)

	See further notes and references at the end of the file.

	by John Pye, 2006
*)
MODEL sunpos;
	t IS_A time; (* starting at zero at midnight on 1 Jan *)
		
	L_st IS_A angle; (* standard meridian for the current time zone *)
	L_loc IS_A angle; (* local longitude *)
	phi IS_A angle; (* latitude, north positive *)

	E IS_A time; (* 'Equation of time' correction *)
	E = 229.2{min}*(0.000075 + 0.001868*cos(B) - 0.032077*sin(B)
			- 0.014615*cos(2*B) - 0.04089*sin(2*B) );
	(* where *) B IS_A angle;
	B = t*360{deg}/365{day};

	t_solar IS_A time; (* solar time *)
	t_solar = t - 4{min/deg}*(L_st - L_loc) + E;
	t_long IS_A time;
	t_long = 4{min/deg}*(L_st - L_loc);
	t_corr IS_A time;
	t_corr = E - t_long;

	delta IS_A angle; (* solar declination *)
	delta = 23.45{deg}*sin(360{deg}*(285{day}+t)/(365{day}));

	theta IS_A angle; (* angle of incidence of sun onto inclined surface *)
	omega IS_A angle; (* hour angle, noon=0, afternoon positive *)
	gamma IS_A angle; (* surface azimuth angle, south=0, west positive *)
	beta IS_A angle; (* surface inclination, horiz=0, +90deg=vertical *)	

	omega = (t_solar - 12{h})*(360{deg}/1{d});
	omega_RHP IS_A angle;
	omega_RHP = arccos(cos(omega));

	omega_s IS_A angle;
	omega_s = arccos( - tan(phi)*tan(delta) );

	theta = arccos( 
		sin(delta)*sin(phi)*cos(beta) - sin(delta)*cos(phi)*sin(beta)*cos(gamma)
		+ cos(delta)*cos(phi)*cos(beta)*cos(omega)
		+ cos(delta)*sin(phi)*sin(beta)*cos(gamma)*cos(omega)
		+ cos(delta)*sin(beta)*sin(gamma)*sin(omega)
	);

	theta_z IS_A angle; (* zenith angle of the sun *)
	theta_z = arccos(
		cos(delta)*cos(phi)*cos(omega) + sin(delta)*sin(phi)
	);

	costheta IS_A factor;
	costheta = cos(theta);
	costhetaz IS_A factor;
	costhetaz = cos(theta_z);

	Gsc IS_A real_constant; (* solar constant *)
	Gsc :== 1353 {W/m^2};

	(* extraterrestrial irradiance on a normal surface*)
	Gon IS_A irradiance;
	Gon = Gsc*(1 + 0.033*cos(
		300{deg}/365{day}*(1{day}+t))
	);

	R_b IS_A factor; (* ratio of beam radiation, inclined : horizontal *)
	R_b = cos(theta) / cos(theta_z);

METHODS
METHOD default_self;
	t := 0{day};
	L_st := -120{deg}; (* USA Pacific time *)
	L_loc := -(116{deg}+47{arcmin}); (* Daggett, California, home of SEGS *)
END default_self;
METHOD on_load;
	RUN default_self;
	RUN reset;
	RUN values;
END on_load;
METHOD specify;
	FIX t, L_st, L_loc, phi; (* time and location *)
	FIX beta, gamma; (* surface orientation *)
END specify;
METHOD values;
	(* values for Duffie & Beckman examples 1.5.1 ff *)
	L_st := -90{deg}; (* USA Central time*)
	L_loc := -89.4{deg};
	phi := +43{deg};
	(* t := 32.4375 {d}; *)
	t := 32{d} + 10{h}+30{min};

	(* surface orientation *)
	beta := 45{deg};
	gamma := 15{deg};

END values;
END sunpos;

(*
	Total daily extraterrestrial radiation, from sunrise to sunset.
	
	Set 't' to the start of the day in question.
*)
MODEL dailyextraterrestrial REFINES sunpos;
	(* extraterrestrial irradiance on a horizontal surface -- problems coz it goes negative! *)
	(* Go IS_A irradiance;
	Go = Gsc*(1 + 0.033*cos( 300{deg}/365{day}*(1{day}+t)) )*cos(theta_z); *)

	(* day's-total extraterrestrial radiation on a horizontal surface *)
	Ho IS_A irradiation;
	Ho = 1{d}*Gon*(
		1/1{PI} * (cos(phi)*cos(delta)*sin(omega_s) + omega_s*sin(phi)*sin(delta))
	);
END dailyextraterrestrial;

(*------------------------------------------------------------------------------
	EXAMPLES

	The following examples come from chapter one of Duffie and Beckman (1980)
*)

(*
	For Madison (Wisconsin), calculate the angle of incidence of beam radiation
	on a surface at 10:30 AM solar time on February 13, if the surface is
	tilted 45 from the horizontal and pointed 15 degrees west of south.

	checked, this looks OK -- JP 
*)
MODEL example_1_6_1 REFINES sunpos;
METHODS
METHOD specify;
	RUN sunpos::specify;
	FREE t;
	FIX t_solar;
END specify;
METHOD values;
	RUN sunpos::values;
	t_solar := 43{d} + 10{h} + 30{min};
	beta := 45 {deg};
	gamma := 15 {deg};
	L_st := -90{deg}; (* USA Central time*)
	L_loc := -89.4{deg};
	phi := +43{deg};
END values;
METHOD self_test;
	ASSERT abs(theta-35.0{deg}) < 0.15{deg};
	ASSERT abs(delta-(-13.80{deg})) < 0.02{deg};
END self_test;
END example_1_6_1;

(*
	Calculate the zenith angle of the sun at 9:30 AM solar time in Madison on
	Feb 13.

	checked, this looks OK -- JP
*)
MODEL example_1_6_2 REFINES sunpos;
METHODS
METHOD specify;
	RUN sunpos::specify;
	FREE t;
	FIX t_solar;
END specify;
METHOD values;
	RUN sunpos::values;
	t_solar := 43{d} + 9{h} + 30{min};
	L_st := -90{deg}; (* USA Central time*)
	L_loc := -89.4{deg};
	phi := +43{deg};
END values;
METHOD self_test;
	ASSERT abs(theta_z-66.0{deg}) < 0.4{deg};
END self_test;
END example_1_6_2;

(*
	What is the ratio of beam radiiation for the surface and time specified in
	Example 1.6.1 to that on a horizontal surface?

	checked, this looks OK -- JP
*)
MODEL example_1_7_1 REFINES example_1_6_1;
METHODS
METHOD self_test;
	ASSERT abs(R_b-1.6577) < 0.013;
END self_test;
END example_1_7_1;

(* 
	Calculate Rb for a surface at latitude 40 N, at a tilt 30 degrees toward 
	the south, for the hour 9 to 10 (solar time) on Feb 16.

	checked, this looks OK -- JP
*)
MODEL example_1_7_2 REFINES example_1_6_1;
METHODS
METHOD values;
	phi := 40{deg};
	beta := 30{deg};
	gamma := 0{deg};
	t_solar := 46{d} + 9.5{h};
END values;
METHOD self_test;
	ASSERT abs(delta - (-13{deg})) < 0.5{deg};
	ASSERT abs(R_b-1.61) < 0.005;
END self_test;
END example_1_7_2;

(*
	As for Example 1.7.2, but with a surface inclined at 50 degrees

	checked, this looks OK
*)
MODEL example_1_7_3 REFINES example_1_7_2;
METHODS
METHOD values;
	RUN example_1_7_2::values;
	beta := 50{deg};
END values;
METHOD self_test;
	ASSERT abs(delta - (-13{deg})) < 0.5{deg};
	ASSERT abs(R_b-1.80) < 0.05;
END self_test;
END example_1_7_3;

(* checked, looks OK (although the error in Ho is a little high) *)
MODEL example_1_8_1 REFINES dailyextraterrestrial;
METHODS
METHOD values;
	phi := 43{deg};
	t := 104{d};
END values;
METHOD self_test;
	ASSERT abs(omega_s - 98.9{deg}) < 0.05 {deg};
	ASSERT abs(Ho - 33.4{MJ/m^2}) < 0.4{MJ/m^2};
END self_test;
END dailyextraterrestrial;

(* 
	NOTE
	although these calcs are very useful, we are now starting to see why it 
	might be a good idea to do them outside ASCEND in an external library.

	The reason is that when integrating around the clock, some of these
	functions go unphysically negative, eg radiation on a horizontal surface.
	Using external functions we can easily apply such constraints without
	even necessarily needing to cause nonsmoothness in the 'G' values.

	Another reason is that using the sun position calculations enable smarter
	time-wise interpolation of TMY data, preventing non-physical situations
	from arising, eg sun after sunset.

	--

	REFERENCES

	Duffie & Beckman (1980) Solar Engineering of Thermal Processes, 
	1st ed., Wiley

	Duffie & Beckman (1991) Solar Engineering of Thermal Processes, 
	2nd ed., Wiley

	Iqbal (1983) An Introduction to Solar Radiation, Academic Press, Toronto.
*)
1	johnpye	822	REQUIRE "atoms.a4l";
2
3			(*
4			Using the equations given in Duffie & Beckman (1980)
5			the following model allows calculation of the angle of incidence of
6			sunlight onto an inclined surface
7
8	jpye	1974	The equation of time has been changed to the Duffie & Beckman 2nd edition
9	johnpye	822	form, cited as coming from Iqbal 1983. The equation for 'B' has also
10			been changed to the 2nd Ed D&B form.
11
12			For the longitude correction we use the spherical geometry convention
13			of making east positive. D&B are americocentric, obviously :-)
14
15	jpye	1974	See further notes and references at the end of the file.
16
17	johnpye	822	by John Pye, 2006
18			*)
19			MODEL sunpos;
20			t IS_A time; (* starting at zero at midnight on 1 Jan *)
21
22			L_st IS_A angle; (* standard meridian for the current time zone *)
23			L_loc IS_A angle; (* local longitude *)
24			phi IS_A angle; (* latitude, north positive *)
25
26			E IS_A time; (* 'Equation of time' correction *)
27			E = 229.2{min}(0.000075 + 0.001868cos(B) - 0.032077*sin(B)
28			- 0.014615cos(2B) - 0.04089sin(2B) );
29			(* where *) B IS_A angle;
30			B = t*360{deg}/365{day};
31
32			t_solar IS_A time; (* solar time *)
33			t_solar = t - 4{min/deg}*(L_st - L_loc) + E;
34			t_long IS_A time;
35			t_long = 4{min/deg}*(L_st - L_loc);
36			t_corr IS_A time;
37			t_corr = E - t_long;
38
39			delta IS_A angle; (* solar declination *)
40			delta = 23.45{deg}sin(360{deg}(285{day}+t)/(365{day}));
41
42			theta IS_A angle; (* angle of incidence of sun onto inclined surface *)
43			omega IS_A angle; (* hour angle, noon=0, afternoon positive *)
44			gamma IS_A angle; (* surface azimuth angle, south=0, west positive *)
45			beta IS_A angle; (* surface inclination, horiz=0, +90deg=vertical *)
46
47			omega = (t_solar - 12{h})*(360{deg}/1{d});
48			omega_RHP IS_A angle;
49			omega_RHP = arccos(cos(omega));
50
51			omega_s IS_A angle;
52			omega_s = arccos( - tan(phi)*tan(delta) );
53
54			theta = arccos(
55			sin(delta)sin(phi)cos(beta) - sin(delta)cos(phi)sin(beta)*cos(gamma)
56			+ cos(delta)cos(phi)cos(beta)*cos(omega)
57			+ cos(delta)sin(phi)sin(beta)cos(gamma)cos(omega)
58			+ cos(delta)sin(beta)sin(gamma)*sin(omega)
59			);
60
61			theta_z IS_A angle; (* zenith angle of the sun *)
62			theta_z = arccos(
63			cos(delta)cos(phi)cos(omega) + sin(delta)*sin(phi)
64			);
65
66			costheta IS_A factor;
67			costheta = cos(theta);
68			costhetaz IS_A factor;
69			costhetaz = cos(theta_z);
70
71			Gsc IS_A real_constant; (* solar constant *)
72			Gsc :== 1353 {W/m^2};
73
74			(* extraterrestrial irradiance on a normal surface*)
75			Gon IS_A irradiance;
76			Gon = Gsc(1 + 0.033cos(
77			300{deg}/365{day}*(1{day}+t))
78			);
79
80			R_b IS_A factor; (* ratio of beam radiation, inclined : horizontal *)
81			R_b = cos(theta) / cos(theta_z);
82
83			METHODS
84			METHOD default_self;
85			t := 0{day};
86			L_st := -120{deg}; (* USA Pacific time *)
87			L_loc := -(116{deg}+47{arcmin}); (* Daggett, California, home of SEGS *)
88			END default_self;
89			METHOD on_load;
90			RUN default_self;
91			RUN reset;
92			RUN values;
93			END on_load;
94			METHOD specify;
95			FIX t, L_st, L_loc, phi; (* time and location *)
96			FIX beta, gamma; (* surface orientation *)
97			END specify;
98			METHOD values;
99			(* values for Duffie & Beckman examples 1.5.1 ff *)
100			L_st := -90{deg}; (* USA Central time*)
101			L_loc := -89.4{deg};
102			phi := +43{deg};
103			(* t := 32.4375 {d}; *)
104			t := 32{d} + 10{h}+30{min};
105
106			(* surface orientation *)
107			beta := 45{deg};
108			gamma := 15{deg};
109
110			END values;
111			END sunpos;
112
113			(*
114			Total daily extraterrestrial radiation, from sunrise to sunset.
115
116			Set 't' to the start of the day in question.
117			*)
118			MODEL dailyextraterrestrial REFINES sunpos;
119			(* extraterrestrial irradiance on a horizontal surface -- problems coz it goes negative! *)
120			(* Go IS_A irradiance;
121			Go = Gsc(1 + 0.033cos( 300{deg}/365{day}(1{day}+t)) )cos(theta_z); *)
122
123			(* day's-total extraterrestrial radiation on a horizontal surface *)
124			Ho IS_A irradiation;
125			Ho = 1{d}Gon(
126			1/1{PI} * (cos(phi)cos(delta)sin(omega_s) + omega_ssin(phi)sin(delta))
127			);
128			END dailyextraterrestrial;
129
130			(*------------------------------------------------------------------------------
131			EXAMPLES
132
133			The following examples come from chapter one of Duffie and Beckman (1980)
134			*)
135
136			(*
137			For Madison (Wisconsin), calculate the angle of incidence of beam radiation
138			on a surface at 10:30 AM solar time on February 13, if the surface is
139			tilted 45 from the horizontal and pointed 15 degrees west of south.
140
141			checked, this looks OK -- JP
142			*)
143			MODEL example_1_6_1 REFINES sunpos;
144			METHODS
145			METHOD specify;
146			RUN sunpos::specify;
147			FREE t;
148			FIX t_solar;
149			END specify;
150			METHOD values;
151			RUN sunpos::values;
152			t_solar := 43{d} + 10{h} + 30{min};
153			beta := 45 {deg};
154			gamma := 15 {deg};
155			L_st := -90{deg}; (* USA Central time*)
156			L_loc := -89.4{deg};
157			phi := +43{deg};
158			END values;
159			METHOD self_test;
160			ASSERT abs(theta-35.0{deg}) < 0.15{deg};
161			ASSERT abs(delta-(-13.80{deg})) < 0.02{deg};
162			END self_test;
163			END example_1_6_1;
164
165			(*
166			Calculate the zenith angle of the sun at 9:30 AM solar time in Madison on
167			Feb 13.
168
169			checked, this looks OK -- JP
170			*)
171			MODEL example_1_6_2 REFINES sunpos;
172			METHODS
173			METHOD specify;
174			RUN sunpos::specify;
175			FREE t;
176			FIX t_solar;
177			END specify;
178			METHOD values;
179			RUN sunpos::values;
180			t_solar := 43{d} + 9{h} + 30{min};
181			L_st := -90{deg}; (* USA Central time*)
182			L_loc := -89.4{deg};
183			phi := +43{deg};
184			END values;
185			METHOD self_test;
186			ASSERT abs(theta_z-66.0{deg}) < 0.4{deg};
187			END self_test;
188			END example_1_6_2;
189
190			(*
191			What is the ratio of beam radiiation for the surface and time specified in
192			Example 1.6.1 to that on a horizontal surface?
193
194			checked, this looks OK -- JP
195			*)
196			MODEL example_1_7_1 REFINES example_1_6_1;
197			METHODS
198			METHOD self_test;
199			ASSERT abs(R_b-1.6577) < 0.013;
200			END self_test;
201			END example_1_7_1;
202
203			(*
204			Calculate Rb for a surface at latitude 40 N, at a tilt 30 degrees toward
205			the south, for the hour 9 to 10 (solar time) on Feb 16.
206
207			checked, this looks OK -- JP
208			*)
209			MODEL example_1_7_2 REFINES example_1_6_1;
210			METHODS
211			METHOD values;
212			phi := 40{deg};
213			beta := 30{deg};
214			gamma := 0{deg};
215			t_solar := 46{d} + 9.5{h};
216			END values;
217			METHOD self_test;
218			ASSERT abs(delta - (-13{deg})) < 0.5{deg};
219			ASSERT abs(R_b-1.61) < 0.005;
220			END self_test;
221			END example_1_7_2;
222
223			(*
224			As for Example 1.7.2, but with a surface inclined at 50 degrees
225
226			checked, this looks OK
227			*)
228			MODEL example_1_7_3 REFINES example_1_7_2;
229			METHODS
230			METHOD values;
231			RUN example_1_7_2::values;
232			beta := 50{deg};
233			END values;
234			METHOD self_test;
235			ASSERT abs(delta - (-13{deg})) < 0.5{deg};
236	johnpye	974	ASSERT abs(R_b-1.80) < 0.05;
237	johnpye	822	END self_test;
238			END example_1_7_3;
239
240			(* checked, looks OK (although the error in Ho is a little high) *)
241			MODEL example_1_8_1 REFINES dailyextraterrestrial;
242			METHODS
243			METHOD values;
244			phi := 43{deg};
245			t := 104{d};
246			END values;
247			METHOD self_test;
248			ASSERT abs(omega_s - 98.9{deg}) < 0.05 {deg};
249			ASSERT abs(Ho - 33.4{MJ/m^2}) < 0.4{MJ/m^2};
250			END self_test;
251			END dailyextraterrestrial;
252
253			(*
254			NOTE
255			although these calcs are very useful, we are now starting to see why it
256			might be a good idea to do them outside ASCEND in an external library.
257
258			The reason is that when integrating around the clock, some of these
259			functions go unphysically negative, eg radiation on a horizontal surface.
260			Using external functions we can easily apply such constraints without
261			even necessarily needing to cause nonsmoothness in the 'G' values.
262
263			Another reason is that using the sun position calculations enable smarter
264			time-wise interpolation of TMY data, preventing non-physical situations
265			from arising, eg sun after sunset.
266
267			--
268
269			REFERENCES
270
271			Duffie & Beckman (1980) Solar Engineering of Thermal Processes,
272			1st ed., Wiley
273
274			Duffie & Beckman (1991) Solar Engineering of Thermal Processes,
275			2nd ed., Wiley
276
277			Iqbal (1983) An Introduction to Solar Radiation, Academic Press, Toronto.
278			*)