WSL/SLF GitLab Repository

libinterpol2D.cc 23.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/***********************************************************************************/
/*  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/>.
*/
//This is the two 2D meteo interpolation library.
#include <meteoio/libinterpol2D.h>

using namespace std;

namespace mio {

25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
//Quake3 fast 1/x² approximation
// For Magic Derivation see: Chris Lomont http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
// Credited to Greg Walsh.
// 32  Bit float magic number 
//#define SQRT_MAGIC_F 0x5f3759df
#define SQRT_MAGIC_F 0x5f375a86

//maximum relative error is <1.7% while computation time for sqrt is <1/4. At 0, returns a large number
//on a large scale interpolation test on TA, max relative error is 1e-6
inline float invSqrt(const float x) {
	const float xhalf = 0.5f*x;
	
	union {
		// get bits for floating value
		float x;
		int i;
	} u;
	u.x = x;
	u.i = SQRT_MAGIC_F - (u.i >> 1);  // gives initial guess y0
	return u.x*(1.5f - xhalf*u.x*u.x);// Newton step, repeating increases accuracy
}

inline float fastSqrt_Q3(const float x) {
	return x * invSqrt(x);
}

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
const double Interpol2D::dflt_temperature_lapse_rate = -0.0065;
const double Interpol2D::wind_ys = 0.58;
const double Interpol2D::wind_yc = 0.42;

double Interpol2D::getReferenceAltitude(const DEMObject& dem)
{
	double ref_altitude = 1500.;

	if(dem.min_altitude!=IOUtils::nodata && dem.max_altitude!=IOUtils::nodata) {
		//we use the median elevation as the reference elevation for reprojections
		ref_altitude = 0.5 * (dem.min_altitude+dem.max_altitude); 
	} else {
		//since there is nothing else that we can do, we use an arbitrary elevation
		ref_altitude = 1500.;
	}
	return ref_altitude;
}

//Usefull functions
/**
* @brief Computes the horizontal distance between points, given by coordinates in a geographic grid
* @param X1 (const double) first point's X coordinate
* @param Y1 (const double) first point's Y coordinate
* @param X2 (const double) second point's X coordinate
* @param Y2 (const double) second point's Y coordinate
* @return (double) distance in m
*/
double Interpol2D::HorizontalDistance(const double& X1, const double& Y1, const double& X2, const double& Y2)
{
	//This function computes the horizontaldistance between two points
	//coordinates are given in a square, metric grid system
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
	const double DX=(X1-X2), DY=(Y1-Y2);
	return sqrt( DX*DX + DY*DY );
}

/**
* @brief Computes the 1/horizontal distance between points, given by coordinates in a geographic grid
* @param X1 (const double) first point's X coordinate
* @param Y1 (const double) first point's Y coordinate
* @param X2 (const double) second point's X coordinate
* @param Y2 (const double) second point's Y coordinate
* @return (double) 1/distance in m
*/
double Interpol2D::InvHorizontalDistance(const double& X1, const double& Y1, const double& X2, const double& Y2)
{
	//This function computes 1/horizontaldistance between two points
	//coordinates are given in a square, metric grid system
	const double DX=(X1-X2), DY=(Y1-Y2);
	return invSqrt( DX*DX + DY*DY ); //we use the optimized approximation for 1/sqrt
100
101
102
103
104
105
106
107
108
109
110
}

/**
* @brief Computes the horizontal distance between points, given by their cells indexes
* @param X1 (const double) first point's i index
* @param Y1 (const double) first point's j index
* @param X2 (const double) second point's X coordinate
* @param Y2 (const double) second point's Y coordinate
* @return (double) distance in m
*/
double Interpol2D::HorizontalDistance(const DEMObject& dem, const int& i, const int& j, const double& X2, const double& Y2)
111
{//TODO: store DEM easting/northing, etc as private members
112
113
114
115
116
	//This function computes the horizontal distance between two points
	//coordinates are given in a square, metric grid system
	//for grid points toward real coordinates
	const double X1 = (dem.llcorner.getEasting()+i*dem.cellsize);
	const double Y1 = (dem.llcorner.getNorthing()+j*dem.cellsize);
117
118
	const double DX=(X1-X2), DY=(Y1-Y2);
	return sqrt( DX*DX + DY*DY );
119
120
}

121
122
123
124
125
126
127
//these weighting functions take the square of a distance as an argument and return a weight
double Interpol2D::weightInvDist(const double& d2)
{
	return invSqrt( d2 ); //we use the optimized approximation for 1/sqrt
}
double Interpol2D::weightInvDistSqrt(const double& d2)
{
128
	return fastSqrt_Q3( invSqrt(d2) ); //we use the optimized approximation for 1/sqrt
129
130
131
132
133
134
135
136
137
}
double Interpol2D::weightInvDist2(const double& d2)
{
	return 1./d2; //we use the optimized approximation for 1/sqrt
}
double Interpol2D::weightInvDistN(const double& d2)
{
	return (double)pow( invSqrt(d2) , (float)dist_pow); //we use the optimized approximation for 1/sqrt
}
138
139
140
141
142
143
144
145
146
147

//Data regression models
/**
* @brief Computes the linear regression coefficients fitting the points given as X and Y in two vectors
* the linear regression has the form Y = aX + b with a regression coefficient r. If the regression coefficient is not good enough, a bad point is looked removed.
* @param X (vector\<double\>) vector of X coordinates
* @param Y (vector\<double\>) vector of Y coordinates (same order as X)
* @param coeffs (vector\<double\>) a,b,r coefficients in a vector
* @return (int) EXIT_SUCCESS or EXIT_FAILURE
*/
148
int Interpol2D::LinRegression(const std::vector<double>& in_X, const std::vector<double>& in_Y, std::vector<double>& coeffs)
149
150
151
{
	//finds the linear regression for points (x,y,z,Value)
	const double r_thres=0.7;
152
	//we want at least 4 points AND 85% of the initial data set kept in the regression
153
	const unsigned int min_dataset=(unsigned int)floor(0.85*(double)in_X.size());
154
	const unsigned int min_pts=(min_dataset>4)?min_dataset:4;
155
	const unsigned int nb_pts = in_X.size();
156
	double a,b,r;
157

158
159
	if (nb_pts==2) {
		std::cout << "[W] only two points for linear regression!\n";
160
	}
161
162
	if(nb_pts<2) { //this should not be needed, we should have refrained from calling LinRegression in such a case
		std::cerr << "[E] Not enough data point for linear regression!\n";
163
		coeffs[1]=0.;
164
		coeffs[2]=in_X[1];
165
166
167
168
		coeffs[3]=1.;
		return EXIT_FAILURE;
	}

169
	Interpol1D::LinRegression(in_X, in_Y, coeffs[1], coeffs[2], coeffs[3]);
170
171
172
173
	if(fabs(coeffs[3])>=r_thres)
		return EXIT_SUCCESS;

	std::vector<double> X(in_X), Y(in_Y);
174
	unsigned int nb_valid_pts=nb_pts;
175

176
	while(fabs(coeffs[3])<r_thres && nb_valid_pts>min_pts) {
177
178
179
180
		//we try to remove the one point in the data set that is the worst
		coeffs[3]=0.;
		unsigned int index_bad=0;
		for (unsigned int i=0; i<nb_pts; i++) {
181
182
183
			//invalidating alternatively each point
			const double Y_tmp=Y[i]; Y[i]=IOUtils::nodata;
			Interpol1D::LinRegression(X, Y, a, b, r);
184
185
			Y[i]=Y_tmp;

186
187
188
189
			if (fabs(r)>fabs(coeffs[3])) {
				coeffs[1]=a;
				coeffs[2]=b;
				coeffs[3]=r;
190
				index_bad=i;
191
192
			}
		}
193
		//the worst point has been found, we overwrite it
194
195
		Y[index_bad]=IOUtils::nodata;
		nb_valid_pts--;
196
	}
197

198
199
	//check if r is reasonnable
	if (fabs(coeffs[3])<r_thres) {
200
		std::cout << "[W] Poor regression coefficient: " << std::setprecision(4) << coeffs[3] << "\n";
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	}

	return EXIT_SUCCESS;
}


//Now, the core interpolation functions: they project a given parameter to a reference altitude, given a constant lapse rate
//example: Ta projected to 1500m with a rate of -0.0065K/m
/**
* @brief Projects a given parameter to another elevation: 
* This implementation keeps the value constant as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::ConstProject(const double& value, const double&, const double&, const std::vector<double>&)
{
	return value;
}

/**
* @brief Projects a given parameter to another elevation: 
* This implementation assumes a linear dependency of the value as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
* @param altitude altitude of the original value
* @param new_altitude altitude of the reprojected value
* @param coeffs coefficients to use for the projection
* @return reprojected value
*/
double Interpol2D::LinProject(const double& value, const double& altitude, const double& new_altitude, const std::vector<double>& coeffs)
{
	//linear lapse: coeffs must have been already computed
	if (coeffs.size()<1) {
		throw IOException("Linear regression coefficients not initialized", AT);
	}
	return (value + coeffs[1] * (new_altitude - altitude));
}

//Filling Functions
/**
* @brief Grid filling function: 
* This implementation builds a standard air pressure as a function of the elevation
* @param dem array of elevations (dem)
* @param grid 2D array to fill
*/
void Interpol2D::stdPressureGrid2DFill(const DEMObject& dem, Grid2DObject& grid) {
251
                                       grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
252
253

	//provide each point with an altitude dependant pressure... it is worth what it is...
254
255
	for (unsigned int j=0; j<grid.nrows; j++) {
		for (unsigned int i=0; i<grid.ncols; i++) {
256
			const double& cell_altitude=dem.grid2D(i,j);
257
258
			if (cell_altitude!=IOUtils::nodata) {
				grid.grid2D(i,j) = lw_AirPressure(cell_altitude);
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

/**
* @brief Grid filling function:
* This implementation fills the grid with a constant value
* @param value value to put in the grid
* @param dem array of elevations (dem). This is needed in order to know if a point is "nodata"
* @param grid 2D array to fill
*/
void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, Grid2DObject& grid)
{
	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);

	//fills a data table with constant values
278
279
	for (unsigned int j=0; j<grid.nrows; j++) {
		for (unsigned int i=0; i<grid.ncols; i++) {
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
			if (dem.grid2D(i,j)!=IOUtils::nodata) {
				grid.grid2D(i,j) = value;
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

/**
* @brief Grid filling function: 
* This implementation fills a flat grid with a constant value and then reprojects it to the terrain's elevation.
* for example, the air temperature measured at one point at 1500m would be given as value, the 1500m as altitude and the dem would allow to reproject this temperature on the full DEM using the detrending function provided as pointer (with its previously calculated coefficients).
* @param value value to put in the grid
* @param altitude altitude of the "value"
* @param dem array of elevations (dem)
* @param vecCoefficients vector of detrending coefficients
* @param funcptr detrending function pointer (that uses the detrending coefficients)
* @param grid 2D array to fill
*/
void Interpol2D::constantLapseGrid2DFill(const double& value, const double& altitude, 
                                         const DEMObject& dem, const std::vector<double>& vecCoefficients,
                                         const LapseRateProjectPtr& funcptr, Grid2DObject& grid)
{
	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);

	//fills a data table with constant values and then reprojects it to the DEM's elevation from a given altitude
	//the laspe rate parameters must have been set before
308
309
	for (unsigned int j=0; j<grid.nrows; j++) {
		for (unsigned int i=0; i<grid.ncols; i++) {
310
311
312
			const double cell_altitude=dem.grid2D(i,j);
			if (cell_altitude!=IOUtils::nodata) {
				grid.grid2D(i,j) = funcptr(value, altitude, cell_altitude, vecCoefficients);
313
314
315
316
317
318
319
320
321
322
323
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<double>& vecData_in,
                           const std::vector<StationData>& vecStations_in)
{
	//The value at any given cell is the sum of the weighted contribution from each source
324
	const unsigned int n_stations=vecStations_in.size();
325
	double parameter=0., norm=0.;
326
	const double scale = 1.e6;
327
	
328
	for (unsigned int i=0; i<n_stations; i++) {
329
330
		/*const double weight=1./(HorizontalDistance(x, y, vecStations_in[i].position.getEasting(),
		       vecStations_in[i].position.getNorthing()) + 1e-6);*/
331
		const Coords& position = vecStations_in[i].position;
332
		const double DX=x-position.getEasting(); //optimization: precompute and store in a vector?
333
334
335
		const double DY=y-position.getNorthing();
		const double weight = invSqrt( DX*DX + DY*DY + scale ); //use the optimized 1/sqrt approximation
		//const double weight = weightInvDistSqrt( (DX*DX + DY*DY) );
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
		parameter += weight*vecData_in[i];
		norm += weight;
	}
	return (parameter/norm); //normalization
}

/**
* @brief Grid filling function: 
* This implementation fills a flat grid using Inverse Distance Weighting and then reproject it to the terrain's elevation.
* for example, the air temperatures measured at several stations would be given as values, the stations altitude and positions
* as positions and projected to a flat grid. Afterward, the grid would be reprojected to the correct elevation as given
* by the dem would using the detrending function provided as pointer (with its previously calculated coefficients).
* @param vecData_in input values to use for the IDW
* @param vecStations_in position of the "values" (altitude and coordinates)
* @param dem array of elevations (dem)
* @param vecCoefficients vector of detrending coefficients
* @param funcptr detrending function pointer (that uses the detrending coefficients)
* @param grid 2D array to fill
*/
void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vector<StationData>& vecStations_in,
                          const DEMObject& dem, const std::vector<double>& vecCoefficients,
                          const LapseRateProjectPtr& funcptr,
                          Grid2DObject& grid)
{	//multiple source stations: lapse rate projection, IDW Krieging, re-projection
	const double ref_altitude = getReferenceAltitude(dem);
361
	const unsigned int n_stations=vecStations_in.size();
362
363
364
365

	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
	std::vector<double> vecTref(vecStations_in.size(), 0.0); // init to 0.0
	
366
	for (unsigned int i=0; i<n_stations; i++) {
367
368
369
370
371
372
373
		vecTref[i] = funcptr(vecData_in[i], vecStations_in[i].position.getAltitude(), 
		                     ref_altitude, vecCoefficients);
	}

	const double xllcorner = dem.llcorner.getEasting();
	const double yllcorner = dem.llcorner.getNorthing();
	const double cellsize = dem.cellsize;
374
375
	for (unsigned int j=0; j<grid.nrows; j++) {
		for (unsigned int i=0; i<grid.ncols; i++) {
376
377
			const double cell_altitude=dem.grid2D(i,j);
			if (cell_altitude!=IOUtils::nodata) {
378
379
380
				grid.grid2D(i,j) = IDWCore((xllcorner+i*cellsize),
				                           (yllcorner+j*cellsize), vecTref, vecStations_in);
				grid.grid2D(i,j) = funcptr(grid.grid2D(i,j), ref_altitude,
381
				                           cell_altitude, vecCoefficients);
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

/**
* @brief Grid filling function: 
* This implementation fills a grid using Inverse Distance Weighting.
* for example, the air temperatures measured at several stations would be given as values, the stations positions
* as positions and projected to a grid. No elevation detrending is performed, the DEM is only used for checking if a grid point is "nodata".
* @param vecData_in input values to use for the IDW
* @param vecStations_in position of the "values" (altitude and coordinates)
* @param dem array of elevations (dem). This is needed in order to know if a point is "nodata"
* @param grid 2D array to fill
*/
void Interpol2D::IDW(const std::vector<double>& vecData_in, const std::vector<StationData>& vecStations_in,
                     const DEMObject& dem, Grid2DObject& grid)
{
	grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);

	//multiple source stations: simple IDW Krieging
	const double xllcorner = dem.llcorner.getEasting();
	const double yllcorner = dem.llcorner.getNorthing();
	const double cellsize = dem.cellsize;
408
409
	for (unsigned int j=0; j<grid.nrows; j++) {
		for (unsigned int i=0; i<grid.ncols; i++) {
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
			if (dem.grid2D(i,j)!=IOUtils::nodata) {
				grid.grid2D(i,j) = IDWCore((xllcorner+i*cellsize), (yllcorner+j*cellsize),
				                           vecData_in, vecStations_in);
			} else {
				grid.grid2D(i,j) = IOUtils::nodata;
			}
		}
	}
}

/**
* @brief Grid filling function: 
* This implementation fills a grid using a curvature and slope algorithm, as described in "A Meteorological
* Distribution System for High-Resolution Terrestrial Modeling (MicroMet)", Liston and Elder, 2006.
* @param dem array of elevations (dem). The slope must have been updated as it is required for the DEM analysis.
* @param VW 2D array of Wind Velocity to fill
* @param DW 2D array of Wind Direction to fill
*/
void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW, Grid2DObject& DW)
{
	if ((!VW.isSameGeolocalization(DW)) || (!VW.isSameGeolocalization(dem))){
		throw IOException("Requested grid VW and grid DW don't match the geolocalization of the DEM", AT);
	}

	//This method computes the speed of the wind and returns a table in 2D with this values
	double speed;		// Wind speed (m s-1)
	double dir;		// Wind direction
	double u;		// Zonal component u (m s-1)
	double v;		// Meridional component v (m s-1)
	double beta;		// Terrain slope
	double azi;		// Topographic slope azimuth
	double curvature;	// Topographic curvature
	double slopeDir;	// Slope in the direction of the wind
	double Ww;		// Wind weighting
	double Od;		// Diverting factor
445
446
447

	const double dem_min_slope=dem.min_slope, dem_range_slope=(dem.max_slope-dem_min_slope);
	const double dem_min_curvature=dem.min_curvature, dem_range_curvature=(dem.max_curvature-dem_min_curvature);
448
	
449
450
	for (unsigned int j=0;j<VW.nrows-1;j++) {
		for (unsigned int i=0;i<VW.ncols-1;i++){
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
			// Get input data
			speed = VW.grid2D(i,j);
			dir = DW.grid2D(i,j);
			beta = dem.slope(i, j);
			azi = dem.azi(i, j);
			curvature = dem.curvature(i, j);

			if(speed==IOUtils::nodata || dir==IOUtils::nodata || beta==IOUtils::nodata || azi==IOUtils::nodata || curvature==IOUtils::nodata) {
				VW.grid2D(i, j) = IOUtils::nodata;
				DW.grid2D(i, j) = IOUtils::nodata;
			} else {
				//convert direction to rad
				dir *= ((M_PI) / 180.);
				//Speed and direction converted to zonal et meridional
				//components 
				u = (-1.) * (speed * sin(dir));
				v = (-1.) * (speed * cos(dir));

				// Converted back to speed and direction
				speed = sqrt(u*u + v*v);
				dir = (1.5 * M_PI) - atan(v/u);

				//normalize curvature and beta. 
				//Note: it should be slopeDir instead of beta, but beta is more efficient
				//to compute (only once for each dem) and it should not be that different...
476
477
				beta = (beta - dem_min_slope)/dem_range_slope - 0.5;
				curvature = (curvature - dem_min_curvature)/dem_range_curvature - 0.5;
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610

				// Calculate the slope in the direction of the wind
				slopeDir = beta * cos(dir - azi);
		
				// Calculate the wind weighting factor
				Ww = 1. + wind_ys * slopeDir + wind_yc * curvature;

				// Modify the wind direction by a diverting factor
				Od = -0.5 * slopeDir * sin(2.*(azi - dir));

				// Calculate the terrain-modified wind speed
				VW.grid2D(i, j) = Ww * speed;

				// Add the diverting factor to the wind direction and convert to degrees
				DW.grid2D(i, j) = (dir + Od) * (180. / (M_PI));
				if( DW.grid2D(i, j)>360. ) {
					DW.grid2D(i, j) -= 360.;
				}
			}
		}
	}
}

double Interpol2D::RhtoDewPoint(double RH, double TA, const short int& force_water)
{
	//Convert a Relative Humidity into a dew point temperature
	//TA is in Kelvins, RH between 0 and 1, returns Td in Kelvins
	TA = K_TO_C(TA);
	double Es, E, Tdw, Tdi; //saturation and current water vapor pressure
	const double Aw = 611.21, Bw = 17.502, Cw = 240.97;	//parameters for water
	const double Ai = 611.15, Bi = 22.452, Ci = 272.55;	//parameters for ice
	const double Tfreeze = 0.;			//freezing temperature
	const double Tnucl = -16.0;			//nucleation temperature
	const double di = 1. / ((TA - Tnucl) * (TA - Tnucl) + 1e-6);		//distance to pure ice
	const double dw = 1. / ((Tfreeze - TA) * (Tfreeze - TA) + 1e-6);	//distance to pure water

	//in order to avoid getting NaN if RH=0
	RH += 0.0001;
	assert(RH>0.);
	if (TA >= Tfreeze || force_water==1) {//above freezing point, water
		Es = Aw * exp( (Bw * TA) / (Cw + TA) );
		E = RH * Es;
		Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) );
		return C_TO_K(Tdw);
	}
	if (TA < Tnucl) { //below nucleation, ice
		Es = Ai * exp( (Bi * TA) / (Ci + TA) );
		E = RH * Es;
		Tdi = ( Ci * log(E / Ai) ) / ( Bi - log(E / Ai) );
		return C_TO_K(Tdi);
	}

	//no clear state, we do a smooth interpolation between water and ice
	Es = Ai * exp( (Bi*TA) / (Ci + TA) );
	E = RH * Es;
	Tdi = ( Ci * log(E / Ai) ) / ( Bi - log(E / Ai) );

	Es = Aw * exp( (Bw * TA) / (Cw + TA) );
	E = RH * Es;
	Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) );

	return C_TO_K( (di / (di + dw) * Tdi + dw / (di + dw) * Tdw) );
}

double Interpol2D::DewPointtoRh(double TD, double TA, const short int& force_water)
{
	//Convert a dew point temperature into a Relative Humidity
	//TA, TD are in Kelvins, RH is returned between 0 and 1
	TA = K_TO_C(TA);
	TD = K_TO_C(TD);
	double Es, E, Rhi, Rhw, Rh;			//saturation and current water vapro pressure
	const double Aw = 611.21, Bw = 17.502, Cw = 240.97;	//parameters for water
	const double Ai = 611.15, Bi = 22.452, Ci = 272.55;	//parameters for ice
	const double Tfreeze = 0.;			//freezing temperature
	const double Tnucl = -16.0;			//nucleation temperature
	const double di = 1. / ((TA - Tnucl) * (TA - Tnucl) + 1e-6);		//distance to pure ice
	const double dw = 1. / ((Tfreeze - TA) * (Tfreeze - TA) + 1e-6);	//distance to pure water

	if (TA >= Tfreeze || force_water==1) {
		//above freezing point, water
		Es = Aw * exp( (Bw * TA) / (Cw + TA) );
		E  = Aw * exp( (Bw * TD) / (Cw + TD) );
		Rhw = (E / Es);
		if (Rhw > 1.) {
			return 1.;
		} else {
			return Rhw;
		}
	}
	if (TA < Tnucl) {
		//below nucleation, ice
		Es = Ai * exp( (Bi * TA) / (Ci + TA) );
		E  = Ai * exp( (Bi * TD) / (Ci + TD) );
		Rhi = (E / Es);
		if (Rhi > 1.) {
			return 1.;
		} else {
			return Rhi;
		}
	}

	//no clear state, we do a smooth interpolation between water and ice
	Es = Ai * exp( (Bi * TA) / (Ci + TA) );
	E  = Ai * exp( (Bi * TD) / (Ci + TD) );
	Rhi = E / Es;

	Es = Aw * exp( (Bw * TA) / (Cw + TA) );
	E  = Aw * exp( (Bw * TD) / (Cw + TD) );
	Rhw = E / Es;

	Rh = (di / (di + dw) * Rhi + dw / (di + dw) * Rhw);
	if(Rh > 1.) {
		return 1.;
	} else {
		return Rh;
	}
}

double Interpol2D::lw_AirPressure(const double& altitude)
{
	double p;
	const double p0 = 101325.; 		// Air and standard pressure in Pa
	const double lapse_rate = 0.0065;	// K m-1
	const double sea_level_temp = 288.15;	// K
	const double expo = GRAVITY / (lapse_rate * GAS_CONSTANT_AIR);
	const double R0 = 6356766.0;		// Earth's radius in m
	
	p = p0 * pow( 1. - ( (lapse_rate * R0 * altitude) / (sea_level_temp * (R0 + altitude)) ), expo );
	
	return(p);
}

} //namespace