WSL/SLF GitLab Repository

Sun.cc 16.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <cmath>
#include <string>
#include <iostream> //for "fixed"
#include <iomanip> //for "setprecision"
22
23
24
#include <meteoio/meteoLaws/Sun.h>
#include <meteoio/meteoLaws/Atmosphere.h>
#include <meteoio/meteoLaws/Meteoconst.h>
25
26
27
28
29

namespace mio {

const double SunObject::elevation_dftlThreshold = 5.; //in degrees

30
31
32
33
34
SunObject::SunObject(SunObject::position_algo /*alg*/)
           : position(), julian_gmt(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata), altitude(IOUtils::nodata),
             elevation_threshold(elevation_dftlThreshold),
             beam_toa(IOUtils::nodata), beam_direct(IOUtils::nodata), beam_diffuse(IOUtils::nodata)
{
35
36
}

37
38
39
40
41
SunObject::SunObject(const double& i_latitude, const double& i_longitude, const double& i_altitude)
           : position(), julian_gmt(IOUtils::nodata), latitude(i_latitude), longitude(i_longitude), altitude(i_altitude),
             elevation_threshold(elevation_dftlThreshold),
             beam_toa(IOUtils::nodata), beam_direct(IOUtils::nodata), beam_diffuse(IOUtils::nodata)
{
42
43
}

44
45
46
47
48
SunObject::SunObject(const double& i_latitude, const double& i_longitude, const double& i_altitude, const double& i_julian, const double& TZ)
           : position(), julian_gmt(i_julian - TZ/24.), latitude(i_latitude), longitude(i_longitude), altitude(i_altitude),
             elevation_threshold(elevation_dftlThreshold),
             beam_toa(IOUtils::nodata), beam_direct(IOUtils::nodata), beam_diffuse(IOUtils::nodata)
{
49
50
51
	update();
}

52
53
void SunObject::setDate(const double& i_julian, const double& TZ)
{
54
55
56
57
58
59
60
61
62
63
	const double i_julian_gmt = i_julian - TZ/24.;

	if(i_julian_gmt!=julian_gmt) { //invalidate all fields if receiving a new date
		position.reset();
		julian_gmt = i_julian_gmt;
		beam_toa = beam_direct = beam_diffuse = IOUtils::nodata;
	}

	//if the date was new or if the previous date had not lead to an update -> update now
	if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata && altitude!=IOUtils::nodata &&
64
65
	  (beam_toa==IOUtils::nodata || beam_direct==IOUtils::nodata || beam_diffuse==IOUtils::nodata)) {
		update();
66
67
68
	}
}

69
70
void SunObject::setLatLon(const double& i_latitude, const double& i_longitude, const double& i_altitude)
{
71
	position.reset();
72
73
74
	latitude = i_latitude;
	longitude = i_longitude;
	altitude = i_altitude;
75
76
77
78
79
80
	beam_toa = beam_direct = beam_diffuse = IOUtils::nodata;
	if(julian_gmt!=IOUtils::nodata) {
		update();
	}
}

81
82
void SunObject::setElevationThresh(const double& i_elevation_threshold)
{
83
	elevation_threshold = i_elevation_threshold;
84
85
86
	beam_toa = beam_direct = beam_diffuse = IOUtils::nodata;
}

87
88
89
90
double SunObject::getElevationThresh() const {
	return elevation_threshold;
}

91
double SunObject::getJulian(const double& TZ) const {
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
	return (julian_gmt+TZ/24.);
}

void SunObject::calculateRadiation(const double& ta, const double& rh, const double& pressure, const double& mean_albedo) {
//ta in K, rh in [0,1], pressure in Pa, altitude in m
	double azimuth, elevation, eccentricity;
	position.getHorizontalCoordinates(azimuth, elevation, eccentricity);
	getBeamPotential(elevation, eccentricity, ta, rh, pressure, mean_albedo, beam_toa, beam_direct, beam_diffuse);
}

void SunObject::calculateRadiation(const double& ta, const double& rh, const double& mean_albedo) {
//ta in K, rh in [0,1], pressure in Pa, altitude in m
	double azimuth, elevation, eccentricity;
	position.getHorizontalCoordinates(azimuth, elevation, eccentricity);

	const double p = Atmosphere::stdAirPressure(altitude);
	getBeamPotential(elevation, eccentricity, ta, rh, p, mean_albedo, beam_toa, beam_direct, beam_diffuse);
}

void SunObject::update() {
	position.setAll(latitude, longitude, julian_gmt);
}

115
//see http://www.meteoexploration.com/products/solarcalc.php for a validation calculator
116
void SunObject::getBeamPotential(const double& sun_elevation, const double& Eccentricity_corr,
117
                                 const double& ta, const double& rh, const double& pressure, const double& ground_albedo,
118
                                 double& R_toa, double& R_direct, double& R_diffuse) const
119
{
120
	if(ta==IOUtils::nodata || rh==IOUtils::nodata || pressure==IOUtils::nodata || ground_albedo==IOUtils::nodata) {
121
122
123
124
125
		R_toa = IOUtils::nodata;
		R_direct = IOUtils::nodata;
		R_diffuse = IOUtils::nodata;
		return;
	}
126
127
128
129
	if(sun_elevation<0.) { //the Sun is below the horizon, our formulas don't apply
		R_toa = 0.;
		R_direct = 0.;
		R_diffuse = 0.;
130
	} else {
131
		R_toa = Cst::solcon * (1.+Eccentricity_corr);
132
133
134
		getClearSky(sun_elevation, R_toa, ta, rh, pressure, ground_albedo, R_direct, R_diffuse);
	}
}
135

136
137
138
139
void SunObject::getClearSky(const double& sun_elevation, const double& R_toa,
                            const double& ta, const double& rh, const double& pressure, const double& ground_albedo,
                            double& R_direct, double& R_diffuse) const
{
140
	//these pow cost us a lot here, but replacing them by fastPow() has a large impact on accuracy (because of the exp())
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
	const double olt = 0.32;   //ozone layer thickness (cm) U.S.standard = 0.34 cm
	const double w0 = 0.9;     //fraction of energy scattered to total attenuation by aerosols (Bird and Hulstrom(1981))
	const double fc = 0.84;    //fraction of forward scattering to total scattering (Bird and Hulstrom(1981))
	const double alpha = 1.3;  //wavelength exponent (Iqbal(1983) p.118). Good average value: 1.3+/-0.5. Related to the size distribution of the particules
	const double beta = 0.03;  //amount of particules index (Iqbal(1983) p.118). Between 0 & .5 and above.
	const double zenith = 90. - sun_elevation; //this is the TRUE zenith because the elevation is the TRUE elevation
	const double cos_zenith = cos(zenith*Cst::to_rad); //this uses true zenith angle

	// relative optical air mass Young (1994), see http://en.wikipedia.org/wiki/Airmass
	//const double mr = 1. / (cos_zenith + 0.50572 * pow( 96.07995-zenith , -1.6364 )); //pbl: this should use apparent zenith angle, and we only get true zenith angle here...
	// relative optical air mass, Young, A. T. 1994. Air mass and refraction. Applied Optics. 33:1108–1110.
	const double mr = ( 1.002432*cos_zenith*cos_zenith + 0.148386*cos_zenith + 0.0096467) /
	                  ( cos_zenith*cos_zenith*cos_zenith + 0.149864*cos_zenith*cos_zenith
	                  + 0.0102963*cos_zenith +0.000303978);

	// actual air mass: because mr is applicable for standard pressure
	// it is modified for other pressures (in Iqbal (1983), p.100)
	// pressure in Pa
	const double ma = mr * (pressure/101325.);

	// the equations for all the transmittances of the individual atmospheric constituents
	// are from Bird and Hulstrom (1980, 1981) and can be found summarized in Iqbal (1983)
	// on the quoted pages

	// broadband transmittance by Rayleigh scattering (Iqbal (1983), p.189)
	const double taur = exp( -0.0903 * pow(ma,0.84) * (1. + ma - pow(ma,1.01)) );

	// broadband transmittance by ozone (Iqbal (1983), p.189)
	const double u3 = olt * mr; // ozone relative optical path length
	const double alpha_oz = 0.1611 * u3 * pow(1. + 139.48 * u3, -0.3035) -
	                        0.002715 * u3 / ( 1. + 0.044  * u3 + 0.0003 * u3 * u3); //ozone absorbance
	const double tauoz = 1. - alpha_oz;

	// broadband transmittance by uniformly mixed gases (Iqbal (1983), p.189)
	const double taug = exp( -0.0127 * pow(ma, 0.26) );

	// saturation vapor pressure in Pa
	//const double Ps = exp( 26.23 - 5416./ta ); //as used for the parametrization
	const double Ps = Atmosphere::waterSaturationPressure(ta);

	// Leckner (1978) (in Iqbal (1983), p.94), reduced precipitable water
	const double w = 0.493 * rh * Ps / ta;

	// pressure corrected relative optical path length of precipitable water (Iqbal (1983), p.176)
	// pressure and temperature correction not necessary since it is included in its numerical constant
	const double u1 = w * mr;

	// broadband transmittance by water vapor (in Iqbal (1983), p.189)
	const double tauw = 1. - 2.4959 * u1  / (pow(1.0 + 79.034 * u1, 0.6828) + 6.385 * u1);

	// broadband total transmittance by aerosols (in Iqbal (1983), pp.189-190)
	// using Angstroem's turbidity formula Angstroem (1929, 1930) for the aerosol thickness
	// in Iqbal (1983), pp.117-119
	// aerosol optical depth at wavelengths 0.38 and 0.5 micrometer
	const double ka1 = beta * pow(0.38, -alpha);
	const double ka2 = beta * pow(0.5, -alpha);

	// broadband aerosol optical depth:
	const double ka  = 0.2758 * ka1 + 0.35 * ka2;

	// total aerosol transmittance function for the two wavelengths 0.38 and 0.5 micrometer:
	const double taua = exp( -pow(ka, 0.873) * (1. + ka - pow(ka, 0.7088)) * pow(ma, 0.9108) );

	// broadband transmittance by aerosols due to absorption only (Iqbal (1983) p. 190)
	const double tauaa = 1. - (1. - w0) * (1. - ma + pow(ma, 1.06)) * (1. - taua);

	// broadband transmittance function due to aerosols scattering only
	// Iqbal (1983) p. 146 (Bird and Hulstrom (1981))
	const double tauas = taua / tauaa;

	// direct normal solar irradiance in range 0.3 to 3.0 micrometer (Iqbal (1983) ,p.189)
	// 0.9751 is for this wavelength range.
	// Bintanja (1996) (see Corripio (2002)) introduced a correction beta_z for increased
	// transmittance with altitude that is linear up to 3000 m and than fairly constant up to 5000 - 6000 m
	const double beta_z = (altitude<3000.)? 2.2*1.e-5*altitude : 2.2*1.e-5*3000.;

	//Now calculating the radiation
	//Top of atmosphere radiation (it will always be positive, because we check for sun elevation before)
	const double tau_commons = tauoz * taug * tauw * taua;

	// Diffuse radiation from the sky
222
	const double factor = 0.79 * R_toa * tau_commons / (1. - ma + pow( ma,1.02 ));  //avoid recomputing pow() twice
223
	// Rayleigh-scattered diffuse radiation after the first pass through atmosphere (Iqbal (1983), p.190)
224
	const double Idr = factor * 0.5 * (1. - taur );
225
226

	// aerosol scattered diffuse radiation after the first pass through atmosphere (Iqbal (1983), p.190)
227
	const double Ida = factor * fc  * (1. - tauas);
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247

	// cloudless sky albedo Bird and Hulstrom (1980, 1981) (in Iqbal (1983) p. 190)
	//in Iqbal, it is recomputed with ma=1.66*pressure/101325.; and alb_sky=0.0685+0.17*(1.-taua_p)*w0;
	const double alb_sky = 0.0685 + (1. - fc) * (1. - tauas);


	//Now, we compute the direct and diffuse radiation components
	//Direct radiation. All transmitances, including Rayleigh scattering (Iqbal (1983), p.189)
	R_direct = 0.9751*( taur * tau_commons + beta_z ) * R_toa ;

	// multiple reflected diffuse radiation between surface and sky (Iqbal (1983), p.154)
	const double Idm = (Idr + Ida + R_direct) * ground_albedo * alb_sky / (1. - ground_albedo * alb_sky);
	R_diffuse = (Idr + Ida + Idm)*cos_zenith; //Iqbal always "project" diffuse radiation on the horizontal

	if( sun_elevation < elevation_threshold ) {
		//if the Sun is too low on the horizon, we put all the radiation as diffuse
		//the splitting calculation that might take place later on will reflect this
		//instead point radiation, it becomes the radiation of a horizontal sky above the domain
		R_diffuse += R_direct*cos_zenith; //HACK
		R_direct = 0.;
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
	}
}

void SunObject::getBeamRadiation(double& R_toa, double& R_direct, double& R_diffuse) const
{
	R_toa = beam_toa;
	R_direct = beam_direct;
	R_diffuse = beam_diffuse;
}

//diffuse remains beam_diffuse
void SunObject::getHorizontalRadiation(double& R_toa, double& R_direct, double& R_diffuse) const
{
	R_toa = position.getRadiationOnHorizontal(beam_toa);
	R_direct = position.getRadiationOnHorizontal(beam_direct);
	R_diffuse = beam_diffuse;
}

void SunObject::getSlopeRadiation(const double& slope_azi, const double& slope_elev, double& R_toa, double& R_direct, double& R_diffuse) const
{
	R_toa = position.getRadiationOnSlope(slope_azi, slope_elev, beam_toa);
	R_direct = position.getRadiationOnSlope(slope_azi, slope_elev, beam_direct);
	R_diffuse = beam_diffuse;
}

/**
 * @brief Evaluate the splitting coefficient between direct and diffuse components of the
275
276
277
278
 * incoming short wave radiation. Splitting is based on "clearness of the sky", ie. the ratio of
 * measured incoming global radiation to top of the atmosphere radiation toa_h,
 * see
 * D. G. Erbs, S.A. Klein, J.A. Duffie, <i>"Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation"</i>, Solar Energy, <b>28</b>, 4, 1982, Pages 293-302 and
279
280
281
 * M. Iqbal, <i>"An introduction to solar radiation"</i>, 1983, Academic Press,  ISBN: 0-12-373750-8 and
 * D. T. Reindl, W. A. Beckman, J. A. Duffle, <i>"Diffuse fraction correlations</i>, Solar Energy, <b>45</b>, 1990, pp 1-7.
 * @param iswr_modeled modelled radiation, it should be horizontal Top Of Atmosphere Radiation (W/m²)
282
283
284
285
286
287
288
289
290
291
 * @param iswr_measured measured Incoming Short Wave Radiation on the ground (W/m²)
 * @return splitting coefficient (between 0 and 1, 1 being 100% diffuse radiation)
 */
double SunObject::getSplitting(const double& iswr_modeled, const double& iswr_measured) const
{
	double splitting_coef;
	double azimuth, elevation;
	position.getHorizontalCoordinates(azimuth, elevation);

	if( elevation < elevation_threshold ) {
292
		//when the Sun is low above the horizon, Mt is getting abnormaly too large pretending
293
294
295
296
297
		// this is a clear sky day when almost all the radiation should be diffuse
		// no matter how the sky is
		splitting_coef = 1.0;
	} else {
		// clear sky index (ratio global measured to top of atmosphere radiation)
298
		const double Mt = iswr_measured / iswr_modeled; // should be <=1.2, aka clearness index Kt
299
		const double clear_sky = 0.147;
300

301
302
303
304
305
306
307
308
309
		// diffuse fraction: hourly ratio of diffuse to global radiation incident on a horizontal surface
		// splitting according to a combination of Reindl et al.(1990)'s models (Mt-model and Mt&Psolar->elev-model):
		if( Mt >= 0.78 ) { // Mt in [0.78;1] -> clear day
			 splitting_coef = clear_sky;
		} else {
			if( Mt <= 0.3 ) { // Mt in [0;0.3] -> overcast
				splitting_coef = 1.02 - 0.248*Mt;
				if(splitting_coef>1.) splitting_coef=1.;
			} else {           // Mt in ]0.3;0.78[ -> cloudy
310
				splitting_coef = 1.4 - 1.749*Mt + 0.177*sin(elevation*Cst::to_rad);
311
312
313
314
315
316
				if(splitting_coef>0.97) splitting_coef = 0.97;
				if(splitting_coef<clear_sky) splitting_coef = clear_sky;
			}
		}
	}

317
	return (splitting_coef); // should be <=1.1; diff/toa should be <=0.8
318
319
}

320
321
double SunObject::getSplitting(const double& iswr_measured) const
{
322
323
324
	double toa_h, direct_h, diffuse;
	getHorizontalRadiation(toa_h, direct_h, diffuse);
	return getSplitting(toa_h, iswr_measured);
325
326
}

327
const std::string SunObject::toString() const
328
{
329
	std::ostringstream os;
330
	const std::streamsize old_prec = os.precision();
331
	os << "<SunObject>\n";
332
	os << position.toString();
333
	os << std::fixed << std::setprecision(4);
334
335
336
	os << "Julian (gmt)\t" << julian_gmt << "\n";
	os << "Lat/Long/Alt\t" << std::setw(7) << latitude << "° " << std::setw(7) << longitude << "° " << std::setprecision(0) << std::setw(4) << altitude << "\n";
	os << "Elev. thresh.\t" << std::setprecision(1) << elevation_threshold << \n";
337

338
	const int colw=10;
339
340
341
342
343
	os << std::setw(colw) << "" << std::fixed << std::setw(colw) << std::setprecision(1) << "toa";
	os << std::fixed << std::setw(colw) << std::setprecision(1) << "direct";
	os << std::fixed << std::setw(colw) << std::setprecision(1) << "diffuse";
	os << std::fixed << std::setw(colw) << std::setprecision(1) << "sum\n";

344
345
346
347
	os << std::setw(colw) << "Beam" << std::fixed << std::setw(colw) << std::setprecision(1) << beam_toa;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << beam_direct;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << beam_diffuse;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << beam_direct+beam_diffuse << "\n";
348
349

	double R_toa, R_direct, R_diffuse;
350
	getHorizontalRadiation(R_toa, R_direct, R_diffuse);
351
352
353
354
355
	os << std::setw(colw) << "Horizontal" << std::fixed << std::setw(colw) << std::setprecision(1) << R_toa;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << R_direct;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << R_diffuse;
	os << std::fixed << std::setw(colw) << std::setprecision(1) << R_direct+R_diffuse << "\n";

356
357
	os << "</SunObject>\n";
	os.precision(old_prec);
358
	return os.str();
359
360
361
}

} //end namespace