.. _program_listing_file_src_modules_solar.cpp: Program Listing for File solar.cpp ================================== |exhale_lsh| :ref:`Return to documentation for file ` (``src/modules/solar.cpp``) .. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS .. code-block:: cpp // // Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel // modular unstructured mesh based approach for hydrological modelling // Copyright (C) 2018 Christopher Marsh // // This file is part of Canadian Hydrological Model. // // Canadian Hydrological Model is free software: you can redistribute it and/or // modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Canadian Hydrological Model 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 General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Canadian Hydrological Model. If not, see // . // #include "solar.hpp" REGISTER_MODULE_CPP(solar); solar::solar(config_file cfg) : module_base("solar", parallel::data, cfg) { provides("solar_el"); provides("solar_az"); provides_parameter("svf"); } solar::~solar() { } void solar::run(mesh_elem &face) { double Lon = 0; double Lat = 0; if(global_param->is_geographic()) { Lon = face->center().x(); Lat = face->center().y(); } else{ auto& data = face->get_module_data(ID); Lon = data.lng; Lat = data.lat; } //Following the RA DEC to Az Alt conversion sequence explained here: //http://www.stargazing.net/kepler/altaz.html //UTC offset. Don't know how to use datetime's UTC converter yet.... boost::posix_time::time_duration UTC_offset = boost::posix_time::hours(global_param->_utc_offset); std::tm tm = boost::posix_time::to_tm(global_param->posix_time()+UTC_offset); double year = tm.tm_year + 1900.; //convert from epoch double month = tm.tm_mon + 1.;//conert jan == 0 double day = tm.tm_mday; //starts at 1, ok double hour = tm.tm_hour; // 0 = midnight, ok double min = tm.tm_min; // 0, ok double sec = tm.tm_sec; // [0,60] in c++11, ok http://en.cppreference.com/w/cpp/chrono/c/tm double Alt = face->center().z();//0.; //TODO: fix this? if (month <= 2.0) { year = year -1.0; month = month +12.0; } double jd = floor( 365.25*(year + 4716.0)) + floor( 30.6001*( month + 1.0)) + 2.0 - \ floor( year/100.0 ) + floor( floor( year/100.0 )/4.0 ) + day - 1524.5 + \ (hour + min/60. + sec/3600.)/24.; double d = jd-2451543.5; // Keplerian Elements for the Sun (geocentric) double w = 282.9404+4.70935*pow(10,-5)*d; // (longitude of perihelion degrees) double a = 1.000000; // (mean distance, a.u.) double e = 0.016709- 1.151*pow(10.,-9.)*d; // (eccentricity) double M = fmod(356.0470+0.9856002585*d,360.0); // (mean anomaly degrees) double L = w + M; //(Sun's mean longitude degrees) double oblecl = 23.4393-3.563e-7*d; //(Sun's obliquity of the ecliptic) //auxiliary angle double E = M+(180./M_PI)*e*sin(M*(M_PI/180.))*(1+e*cos(M*(M_PI/180.))); //rectangular coordinates in the plane of the ecliptic (x axis toward //perhilion) double x = cos(E*(M_PI/180.))-e; double y = sin(E*(M_PI/180.))*sqrt(1.-e*e); //find the distance and true anomaly double r = sqrt(x*x + y*y); double v = atan2(y,x)*(180./M_PI); //find the longitude of the sun double lon = v + w; //compute the ecliptic rectangular coordinates double xeclip = r*cos(lon*(M_PI/180.)); double yeclip = r*sin(lon*(M_PI/180.)); double zeclip = 0.0; //rotate these coordinates to equitorial rectangular coordinates double xequat = xeclip; double yequat = yeclip*cos(oblecl*(M_PI/180.))+zeclip*sin(oblecl*(M_PI/180.)); double zequat = yeclip*sin(23.4406*(M_PI/180.))+zeclip*cos(oblecl*(M_PI/180.)); //convert equatorial rectangular coordinates to RA and Decl: r = sqrt(xequat*xequat + yequat*yequat + zequat*zequat)-(Alt/149598000.0); //roll up the altitude correction double RA = atan2(yequat,xequat)*(180./M_PI); double delta = asin(zequat/r)*(180./M_PI); //Find the J2000 value double J2000 = jd - 2451545.0; double UTH = hour+min/60.0+sec/3600.0; //Calculate local siderial time double GMST0=fmod(L+180.,360.)/15.; double SIDTIME = GMST0 + UTH + Lon/15.; //Replace RA with hour angle HA double HA = (SIDTIME*15. - RA); //convert to rectangular coordinate system x = cos(HA*(M_PI/180.))*cos(delta*(M_PI/180.)); y = sin(HA*(M_PI/180.))*cos(delta*(M_PI/180.)); double z = sin(delta*(M_PI/180.)); //rotate this along an axis going east-west. double xhor = x*cos((90.-Lat)*(M_PI/180.))-z*sin((90.-Lat)*(M_PI/180.)); double yhor = y; double zhor = x*sin((90.-Lat)*(M_PI/180.))+z*cos((90.-Lat)*(M_PI/180.)); //Find the h and AZ double Az = atan2(yhor,xhor)*(180./M_PI) + 180.; double El = asin(zhor)*(180./M_PI); (*face)["solar_az"_s]=Az; (*face)["solar_el"_s]=El; } void solar::init(mesh& domain) { //number of steps along the search vector to check for a higher point int steps = cfg.get("svf.steps",10); //max distance to search double max_distance = cfg.get("svf.max_distance",1000.0); //size of the step to take double size_of_step = max_distance / steps; //number of azimuthal sections int N = cfg.get("svf.nsectors", 12); double azimuthal_width = 360./(double)N; // in degrees bool svf_compute = cfg.get("svf.compute",true); #pragma omp parallel { OGRSpatialReference monUtm, monGeo; OGRCoordinateTransformation* coordTrans = nullptr; // we are UTM and need to convert internally to lat long to calc the solar position if(!domain->is_geographic()) { // the proj is likely x/y ordering but force it to be x/y so we know what goes where in the assignment below monUtm.importFromProj4(domain->proj4().c_str()); monUtm.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); // "EPSG:4326" is lat/lon (y/x) so force it to be x,y so we know how to do the assignment below monGeo.SetWellKnownGeogCS("EPSG:4326"); monGeo.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); coordTrans = OGRCreateCoordinateTransformation(&monUtm, &monGeo); } #pragma omp for for (size_t i = 0; i < domain->size_local_faces(); i++) { auto face = domain->face(i); // we are UTM and need to convert internally to lat long to calc the solar position if (!domain->is_geographic()) { double x = face->center().x(); double y = face->center().y(); // do the transform with the enforce x/y ordering if(!coordTrans->Transform(1, &x, &y)) { CHM_THROW_EXCEPTION(module_error, "Unable to covert face coordinates to lat long for solar rad calcuation."); } auto& d = face->make_module_data(ID); d.lat = y; d.lng = x; } double svf = 0.0; if (svf_compute) { Point_3 me = face->center(); auto cosSlope = cos(face->slope()); auto sinSlope = sin(face->slope()); // for each search azimuthal sector for (int k = 0; k < N; k++) { double phi = 0.; // search along each azimuth in j step increments to find horizon angle for (int j = 1; j <= steps; ++j) { double distance = j * size_of_step; auto f = domain->find_closest_face(math::gis::point_from_bearing(me, k * azimuthal_width, distance)); double z_diff = (f->center().z() - me.z()); if (z_diff > 0) { double dist = math::gis::distance(f->center(), me); phi = std::max(atan(z_diff / dist), phi); } } auto cosPhi = cos(phi); auto sinPhi = sin(phi); auto azi_in_rad = (k * azimuthal_width * M_PI / 180.); svf += cosSlope * cosPhi * cosPhi + sinSlope * cos(azi_in_rad - face->aspect()) * (M_PI_2 - phi - sinPhi * cosPhi); } svf /= (double)N; } else { svf = 1.; } face->parameter("svf"_s) = std::max(0.0, svf); } OGRCoordinateTransformation::DestroyCT(coordTrans); } }