WSL/SLF GitLab Repository

Commit 3086970d authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Since the _USE_MATH_DEFINES for Windows was making trouble (it should have...

Since the _USE_MATH_DEFINES for Windows was making trouble (it should have enabled using M_PI), the basic math constants have been added to meteolaws/meteoconst.h. Now, please use things like Cst::PI and include meteoconst.h.
parent 51703c49
......@@ -27,9 +27,6 @@ IF(CMAKE_COMPILER_IS_GNUCXX)
SET(WARNINGS "-Wall")
SET(PROFILING "-pg -fprofile-arcs")
SET(EXTRA_WARNINGS "-Wextra -ansi -pedantic")
IF(WIN32)
LIST(APPEND CFLAGS " -D_USE_MATH_DEFINES") #USE_MATH_DEFINES needed for Win32
ENDIF(WIN32)
SET(OPTIM "-g -O3 -DNDEBUG -DNOSAFECHECKS")
EXECUTE_PROCESS(COMMAND ${CMAKE_CXX_COMPILER} -dumpversion OUTPUT_VARIABLE GCC_VERSION)
IF (GCC_VERSION VERSION_GREATER 4.2 OR GCC_VERSION VERSION_EQUAL 4.2)
......@@ -48,7 +45,7 @@ IF(MSVC)
SET(ARCH_OPTIM "/arch:SSE2")
SET(ARCH_SAFE "")
SET(DEBUG "/Z7 /Od /D__DEBUG /MDd")
LIST(APPEND CFLAGS " /D_CRT_SECURE_NO_WARNINGS /D_USE_MATH_DEFINES") #Za: strict ansi
LIST(APPEND CFLAGS " /D_CRT_SECURE_NO_WARNINGS") #Za: strict ansi
SET(_VERSION "/D_VERSION=\\\"${_versionString}\\\"")
ENDIF(MSVC)
......@@ -271,7 +268,7 @@ IF(UNIX)
SET(CPACK_DEBIAN_PACKAGE_ARCHITECTURE "i386") #dpkg --print-architecture
SET(CPACK_RPM_PACKAGE_NAME "meteoio")
SET(CPACK_RPM_PACKAGE_LICENSE "LPGLv3")
SET(CPACK_RPM_PACKAGE_REQUIRES "libstdc++6, libproj4")
SET(CPACK_RPM_PACKAGE_REQUIRES "libstdc++6, libproj4") #simply libproj for OpenSuse
SET(CPACK_RPM_PACKAGE_ARCHITECTURE "i386")
SET(CPACK_RPM_PACKAGE_GROUP "Development/Libraries")
ENDIF(APPLE)
......
......@@ -18,6 +18,7 @@
#include <cmath>
#include <meteoio/Coords.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
#ifdef PROJ4
#include <proj_api.h>
......@@ -982,7 +983,7 @@ void Coords::WGS84_to_UTM(double lat_in, double long_in, double& east_out, doubl
//also http://www.uwgb.edu/dutchs/usefuldata/UTMFormulas.HTM
//also http://www.oc.nps.edu/oc2902w/maps/utmups.pdf or Chuck Gantz (http://www.gpsy.com/gpsinfo/geotoutm/)
//Geometric constants
const double to_rad = M_PI / 180.0;
const double to_rad = Cst::PI / 180.0;
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis
const double e2 = (a*a-b*b) / (a*a); //ellispoid eccentricity, squared
......@@ -1047,7 +1048,7 @@ void Coords::UTM_to_WGS84(double east_in, double north_in, double& lat_out, doub
//also http://www.uwgb.edu/dutchs/usefuldata/UTMFormulas.HTM
//also http://www.oc.nps.edu/oc2902w/maps/utmups.pdf or Chuck Gantz (http://www.gpsy.com/gpsinfo/geotoutm/)
//Geometric constants
const double to_deg = 180.0 / M_PI;
const double to_deg = 180.0 / Cst::PI;
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis
const double e2 = (a*a-b*b) / (a*a); //ellispoid eccentricity, squared
......@@ -1205,7 +1206,7 @@ void Coords::distance(const Coords& destination, double& distance, double& beari
//HACK: this is the 2D distance, it does not work in 3D!!
if(isSameProj(destination)) {
//we can use simple cartesian grid arithmetic
const double to_deg = 180.0 / M_PI;
const double to_deg = 180.0 / Cst::PI;
distance = sqrt( IOUtils::pow2(easting - destination.getEasting()) + IOUtils::pow2(northing - destination.getNorthing()) );
bearing = atan2( northing - destination.getNorthing() , easting - destination.getEasting() );
bearing = fmod( bearing*to_deg+360. , 360. );
......@@ -1233,7 +1234,7 @@ void Coords::distance(const Coords& destination, double& distance, double& beari
void Coords::WGS84_to_local(double lat_in, double long_in, double& east_out, double& north_out) const
{
double alpha;
const double to_rad = M_PI / 180.0;
const double to_rad = Cst::PI / 180.0;
double distance;
if((ref_latitude==IOUtils::nodata) || (ref_longitude==IOUtils::nodata)) {
......@@ -1267,7 +1268,7 @@ void Coords::WGS84_to_local(double lat_in, double long_in, double& east_out, dou
*/
void Coords::local_to_WGS84(double east_in, double north_in, double& lat_out, double& long_out) const
{
const double to_deg = 180.0 / M_PI;
const double to_deg = 180.0 / Cst::PI;
const double distance = sqrt( IOUtils::pow2(east_in) + IOUtils::pow2(north_in) );
const double bearing = fmod( atan2(east_in, north_in)*to_deg+360. , 360.);
......@@ -1312,8 +1313,8 @@ void Coords::WGS84_to_NULL(double /*lat_in*/, double /*long_in*/, double& /*east
void Coords::cosineInverse(const double& lat_ref, const double& lon_ref, const double& distance, const double& bearing, double& lat, double& lon) const
{
const double Rearth = 6371.e3;
const double to_rad = M_PI / 180.0;
const double to_deg = 180.0 / M_PI;
const double to_rad = Cst::PI / 180.0;
const double to_deg = 180.0 / Cst::PI;
const double lat_ref_rad = lat_ref*to_rad;
const double bearing_rad = bearing*to_rad;
......@@ -1328,7 +1329,7 @@ void Coords::cosineInverse(const double& lat_ref, const double& lon_ref, const d
cos(lat_ref_rad)*sin(distance/Rearth)*cos(bearing_rad) );
lon = lon_ref*to_rad + atan2( sin(bearing_rad)*sin(distance/Rearth)*cos(lat_ref_rad) ,
cos(distance/Rearth) - sin(lat_ref_rad)*sin(lat) );
lon = fmod(lon+M_PI, 2.*M_PI) - M_PI;
lon = fmod(lon+Cst::PI, 2.*Cst::PI) - Cst::PI;
lat *= to_deg;
lon *= to_deg;
......@@ -1353,7 +1354,7 @@ double Coords::cosineDistance(const double& lat1, const double& lon1, const doub
}
const double Rearth = 6371.e3;
const double to_rad = M_PI / 180.0;
const double to_rad = Cst::PI / 180.0;
const double d = acos(
sin(lat1*to_rad) * sin(lat2*to_rad)
+ cos(lat1*to_rad) * cos(lat2*to_rad) * cos((lon2-lon1)*to_rad)
......@@ -1386,7 +1387,7 @@ double Coords::VincentyDistance(const double& lat1, const double& lon1, const do
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis
const double f = (a - b) / a; //ellispoid flattening
const double to_rad = M_PI / 180.0;
const double to_rad = Cst::PI / 180.0;
const double L = (lon1 - lon2)*to_rad;
const double U1 = atan( (1.-f)*tan(lat1*to_rad) );
......@@ -1462,8 +1463,8 @@ void Coords::VincentyInverse(const double& lat_ref, const double& lon_ref, const
const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis, value for wgs84
const double b = ellipsoids[E_WGS84].b; //minor ellipsoid semi-axis, value for wgs84
const double f = (a - b) / a; //ellispoid flattening
const double to_rad = M_PI / 180.0;
const double to_deg = 180.0 / M_PI;
const double to_rad = Cst::PI / 180.0;
const double to_deg = 180.0 / Cst::PI;
const double alpha1 = bearing*to_rad;
const double tanU1 = (1.-f)*tan(lat_ref*to_rad);
......@@ -1477,7 +1478,7 @@ void Coords::VincentyInverse(const double& lat_ref, const double& lon_ref, const
const double B = u2/1024. * (256. + u2*(-128.+u2*(74.-47.*u2)));
double sigma = distance / (b*A);
double sigma_p = 2.*M_PI;
double sigma_p = 2.*Cst::PI;
double cos2sigma_m = cos( 2.*sigma1 + sigma ); //required to avoid uninitialized value
while (fabs(sigma - sigma_p) > thresh) {
......
......@@ -19,6 +19,7 @@
#include <limits.h>
#include <meteoio/DEMObject.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
/**
* @file DEMObject.cc
......@@ -593,7 +594,7 @@ double DEMObject::getHorizon(const Coords& point, const double& bearing) {
}
//returning the angle matching the highest tangent
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
return ( atan(max_tangent)*to_deg );
}
......@@ -607,7 +608,7 @@ double DEMObject::getHorizon(const Coords& point, const double& bearing) {
void DEMObject::getHorizon(const Coords& point, const double& increment, std::vector<double>& horizon) {
for(double bearing=0.0; bearing <360.; bearing += increment) {
const double alpha = getHorizon(point, bearing * M_PI/180.);
const double alpha = getHorizon(point, bearing * Cst::PI/180.);
horizon.push_back(alpha);
}
}
......@@ -712,7 +713,7 @@ double DEMObject::CalculateAspect(const double& Nx, const double& Ny, const doub
return (0.); // north facing
}
} else { //there is a E-W slope
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
if ( Nx > 0. ) {
return (90. - atan(Ny/Nx)*to_deg);
} else {
......@@ -738,7 +739,7 @@ void DEMObject::CalculateHick(double A[4][4], double& slope, double& Nx, double&
Nz = IOUtils::nodata;
slope_failures++;
} else {
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
slope = atan(smax)*to_deg;
//Nx and Ny: x and y components of the normal pointing OUT of the surface
......@@ -770,7 +771,7 @@ void DEMObject::CalculateFleming(double A[4][4], double& slope, double& Nx, doub
Nx = 0.5 * (A[2][1] - A[2][3]) / cellsize;
Ny = 0.5 * (A[3][2] - A[1][2]) / cellsize;
Nz = 1.;
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * to_deg;
} else {
CalculateHick(A, slope, Nx, Ny, Nz);
......@@ -789,7 +790,7 @@ void DEMObject::CalculateHorn(double A[4][4], double& slope, double& Nx, double&
//There is no difference between slope = acos(n_z/|n|) and slope = atan(sqrt(sx*sx+sy*sy))
//slope = acos( (Nz / sqrt( Nx*Nx + Ny*Ny + Nz*Nz )) );
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * to_deg;
} else {
//steepest slope method (Dunn and Hickey, 1998)
......@@ -806,7 +807,7 @@ void DEMObject::CalculateCorripio(double A[4][4], double& slope, double& Nx, dou
Nz = 1.;
//There is no difference between slope = acos(n_z/|n|) and slope = atan(sqrt(sx*sx+sy*sy))
//slope = acos( (Nz / sqrt( Nx*Nx + Ny*Ny + Nz*Nz )) );
const double to_deg = 180./M_PI;
const double to_deg = 180./Cst::PI;
slope = atan( sqrt(Nx*Nx+Ny*Ny) ) * to_deg;
} else {
//steepest slope method (Dunn and Hickey, 1998)
......
......@@ -22,6 +22,7 @@
#include <meteoio/Array2D.h>
#include <meteoio/Grid2DObject.h>
#include <meteoio/IOUtils.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
#include <cmath>
#include <limits>
......@@ -68,16 +69,16 @@ class DEMObject : public Grid2DObject {
NORMAL=2, ///< update the normals
CURVATURE=4 ///< update the curvatures
} update_type;
DEMObject(const slope_type& i_algorithm=DFLT);
DEMObject(const unsigned int& ncols_in, const unsigned int& nrows_in,
const double& cellsize_in, const Coords& llcorner_in, const slope_type& i_algorithm=DFLT);
DEMObject(const unsigned int& ncols_in, const unsigned int& nrows_in,
const double& cellsize_in, const Coords& llcorner_in, const Array2D<double>& altitude_in,
const bool& i_update=true, const slope_type& i_algorithm=DFLT);
DEMObject(const Grid2DObject& dem_in, const bool& i_update=true, const slope_type& i_algorithm=DFLT);
DEMObject (const DEMObject& i_dem,
......@@ -104,7 +105,7 @@ class DEMObject : public Grid2DObject {
private:
void CalculateAziSlopeCurve(slope_type algorithm);
double CalculateAspect(const double& Nx, const double& Ny, const double& Nz, const double& slope, const double no_slope=M_PI);
double CalculateAspect(const double& Nx, const double& Ny, const double& Nz, const double& slope, const double no_slope=Cst::PI);
void CalculateHick(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void CalculateFleming(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
void CalculateHorn(double A[4][4], double& slope, double& Nx, double& Ny, double& Nz);
......
......@@ -17,6 +17,7 @@
*/
#include <cmath>
#include <meteoio/FilterAlgorithms.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
using namespace std;
......@@ -1042,14 +1043,14 @@ void FilterAlgorithms::WindAvgProcess(const std::vector<MeteoData>& vecM, const
//calculate ve and vn
double ve=0.0, vn=0.0;
for (unsigned int jj=0; jj<vecSize; jj++){
ve += vecWindowVW[jj] * sin(vecWindowDW[jj] * M_PI / 180.); //turn into radians
vn += vecWindowVW[jj] * cos(vecWindowDW[jj] * M_PI / 180.); //turn into radians
ve += vecWindowVW[jj] * sin(vecWindowDW[jj] * Cst::PI / 180.); //turn into radians
vn += vecWindowVW[jj] * cos(vecWindowDW[jj] * Cst::PI / 180.); //turn into radians
}
ve /= vecSize;
vn /= vecSize;
meanspeed = sqrt(ve*ve + vn*vn);
meandirection = fmod( atan2(ve,vn) * 180. / M_PI + 360. , 360.); // turn into degrees [0;360)
meandirection = fmod( atan2(ve,vn) * 180. / Cst::PI + 360. , 360.); // turn into degrees [0;360)
}
vecWindowM[ii](MeteoData::VW) = meanspeed;
......
/***********************************************************************************/
/* Copyright 2009 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/* Copyright 2011 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
......
......@@ -18,6 +18,7 @@
#include <cmath>
#include <meteoio/IOUtils.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
#include <meteoio/Config.h> // to avoid forward declaration hell
#include <meteoio/MeteoData.h> // to avoid forward declaration hell
......@@ -81,12 +82,12 @@ double IOUtils::pow4(const double& val)
}
double IOUtils::bearing_to_angle(const double& bearing) {
const double to_rad = M_PI/180.;
const double to_rad = Cst::PI/180.;
return (fmod(360.-bearing+90., 360.)*to_rad);
}
double IOUtils::angle_to_bearing(const double& angle) {
const double to_deg = 180.0 / M_PI;
const double to_deg = 180.0 / Cst::PI;
return (fmod( 90.-angle*to_deg+360. , 360. ));
}
......
......@@ -84,14 +84,14 @@ double FilterWindAvg::calc_avg(const unsigned int& index, const std::vector<cons
//calculate ve and vn
double ve=0.0, vn=0.0;
for (size_t jj=0; jj<vecSize; jj++){
ve += vec_window[jj]->operator()(MeteoData::VW) * sin(vec_window[jj]->operator()(MeteoData::DW) * M_PI / 180.); //turn into radians
vn += vec_window[jj]->operator()(MeteoData::VW) * cos(vec_window[jj]->operator()(MeteoData::DW) * M_PI / 180.); //turn into radians
ve += vec_window[jj]->operator()(MeteoData::VW) * sin(vec_window[jj]->operator()(MeteoData::DW) * Cst::PI / 180.); //turn into radians
vn += vec_window[jj]->operator()(MeteoData::VW) * cos(vec_window[jj]->operator()(MeteoData::DW) * Cst::PI / 180.); //turn into radians
}
ve /= vecSize;
vn /= vecSize;
meanspeed = sqrt(ve*ve + vn*vn);
meandirection = fmod( atan2(ve,vn) * 180. / M_PI + 360. , 360.); // turn into degrees [0;360)
meandirection = fmod( atan2(ve,vn) * 180. / Cst::PI + 360. , 360.); // turn into degrees [0;360)
}
if(index==MeteoData::VW)
......
......@@ -44,6 +44,21 @@ namespace Cst {
const double earth_R0 = 6356766.0; // (m)
const double solcon = 1366.1; // (W/m^2)
//Math constants
const double e = 2.71828182845904523536; // e
const double Log2e = 1.44269504088896340736; // log2(e)
const double Log10e = 0.434294481903251827651; // log10(e)
const double Ln2 = 0.693147180559945309417; // ln(2)
const double Ln10 = 2.30258509299404568402; // ln(10)
const double PI = 3.14159265358979323846; // pi
const double PI2 = 1.57079632679489661923; // pi/2
const double PI4 = 0.785398163397448309616; // pi/4
const double InvPI = 0.318309886183790671538; // 1/pi
const double TwoOverPI = 0.636619772367581343076; // 2/pi
const double TwoOverSqrtPI = 1.12837916709551257390; // 2/sqrt(pi)
const double Sqrt2 = 1.41421356237309504880; // sqrt(2)
const double InvSqrt2 = 0.707106781186547524401; // 1/sqrt(2)
} //end CST namespace
} //end namespace
......
......@@ -116,7 +116,7 @@ void SunObject::getBeamPotential(const double& sun_elevation, const double& Ecce
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)
const double beta = 0.03;
const double to_rad = M_PI/180.;
const double to_rad = Cst::PI/180.;
if(ta==IOUtils::nodata || rh==IOUtils::nodata || pressure==IOUtils::nodata || mean_albedo==IOUtils::nodata) {
R_toa = IOUtils::nodata;
......@@ -269,7 +269,7 @@ void SunObject::getSlopeRadiation(const double& slope_azi, const double& slope_e
*/
double SunObject::getSplitting(const double& iswr_modeled, const double& iswr_measured) const
{
const double to_rad = M_PI/180.;
const double to_rad = Cst::PI/180.;
double splitting_coef;
double azimuth, elevation;
position.getHorizontalCoordinates(azimuth, elevation);
......
......@@ -18,11 +18,12 @@
#include <cmath>
#include <string>
#include <meteoio/meteolaws/Suntrajectory.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
namespace mio {
//class SunTrajectory
const double SunTrajectory::to_deg = 180./M_PI;
const double SunTrajectory::to_rad = M_PI/180.;
const double SunTrajectory::to_deg = 180./Cst::PI;
const double SunTrajectory::to_rad = Cst::PI/180.;
/**
* @brief Compute the solar incidence (rad), i.e. the angle between the incident sun beam
......@@ -246,14 +247,14 @@ void SunMeeus::getEquatorialCoordinates(double& right_ascension, double& declina
}
void SunMeeus::getEquatorialSunVector(double& sunx, double& suny, double& sunz) {
const double to_rad = M_PI/180.;
const double to_rad = Cst::PI/180.;
double azi_Sacw;
// Convert to angle measured from South, counterclockwise (rad)
if ( SolarAzimuthAngle <= 90. ) {
azi_Sacw = M_PI - SolarAzimuthAngle*to_rad;
azi_Sacw = Cst::PI - SolarAzimuthAngle*to_rad;
} else {
azi_Sacw = 3.*M_PI - SolarAzimuthAngle*to_rad;
azi_Sacw = 3.*Cst::PI - SolarAzimuthAngle*to_rad;
}
// derived as shown in Funk (1984) p. 107, but for a y-coordinate increasing northwards
......
......@@ -19,6 +19,7 @@
#include <meteoio/meteostats/libinterpol2D.h>
#include <meteoio/meteolaws/Atmosphere.h>
#include <meteoio/meteolaws/Meteoconst.h> //for math constants
using namespace std;
......@@ -609,7 +610,7 @@ void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW
double Ww; // Wind weighting
double Od; // Diverting factor
const double to_rad = M_PI/180.;
const double to_rad = Cst::PI/180.;
const double dem_min_slope=dem.min_slope*to_rad, dem_range_slope=(dem.max_slope-dem_min_slope)*to_rad;
const double dem_min_curvature=dem.min_curvature, dem_range_curvature=(dem.max_curvature-dem_min_curvature);
......@@ -627,7 +628,7 @@ void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW
DW.grid2D(i, j) = IOUtils::nodata;
} else {
//convert direction to rad
dir *= ((M_PI) / 180.);
dir *= ((Cst::PI) / 180.);
//Speed and direction converted to zonal et meridional
//components
u = (-1.) * (speed * sin(dir));
......@@ -635,7 +636,7 @@ void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW
// Converted back to speed and direction
speed = sqrt(u*u + v*v);
dir = (1.5 * M_PI) - atan(v/u);
dir = (1.5 * Cst::PI) - atan(v/u);
//normalize curvature and beta.
//Note: it should be slopeDir instead of beta, but beta is more efficient
......@@ -656,7 +657,7 @@ void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW
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));
DW.grid2D(i, j) = (dir + Od) * (180. / (Cst::PI));
if( DW.grid2D(i, j)>360. ) {
DW.grid2D(i, j) -= 360.;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment