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 changes 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 :-)

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