TerrainRadiationComplex.h

Go to the documentation of this file.

/***********************************************************************************/
/*  Copyright 2009-2015 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of Alpine3D.
    Alpine3D 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.
 
    Alpine3D 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 Alpine3D.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef TERRAINRadiationComplex_H
#define TERRAINRadiationComplex_H
 
#include <meteoio/MeteoIO.h>
 
#include <alpine3d/ebalance/TerrainRadiationAlgorithm.h>
#include <alpine3d/ebalance/RadiationField.h>
#include <alpine3d/ebalance/SolarPanel.h>
#include <alpine3d/ebalance/SnowBRDF.h>
#include <array>
 
class TerrainRadiationComplex : public TerrainRadiationAlgorithm
{
 
public:
        TerrainRadiationComplex(const mio::Config &cfg, const mio::DEMObject &dem_in, const std::string &method);
 
        virtual void getRadiation(mio::Array2D<double>& direct, mio::Array2D<double>& diffuse,
                            mio::Array2D<double>& terrain, const mio::Array2D<double>& direct_unshaded_horizontal,
                            const mio::Array2D<double>& total_ilwr, mio::Array2D<double>& sky_ilwr,
                            mio::Array2D<double>& terrain_ilwr, double solarAzimuth, double solarElevation) override;
        virtual void setMeteo(const mio::Array2D<double>& albedo, const mio::Array2D<double>& ta) override;
 
        void getSkyViewFactor(mio::Array2D<double> &o_sky_vf) override;
        void setSP(const mio::Date timestamp, const double solarAzimuth, const double solarElevation) override;
        void writeSP(const unsigned int max_steps) override;
 
private:
        typedef std::array<double, 3> Vec3D;
 
        // Initialisation Functions
        void initBasicSetHorizontal();
        void initBasicSetRotated();
        void initViewList();
        void initRList();
        void initSortList();
        void WriteViewList();
        bool ReadViewList();
 
        // auxiliary functions
        void TriangleNormal(size_t ii_dem, size_t jj_dem, int which_triangle, Vec3D &v_out);
        double IntersectionRayTriangle(const Vec3D &v_view, size_t ii_0, size_t jj_0, size_t ii_dem, size_t jj_dem, size_t which_triangle);
        size_t vectorToSPixel(const Vec3D &vec_in, size_t ii_dem, size_t jj_dem, size_t which_triangle);
        void initSkyViewFactor();
        double computeSkyViewFactor(size_t ii_dem, size_t jj_dem, size_t which_triangle);
        double getLandViewFactor(size_t ii_dem, size_t jj_dem, size_t which_triangle);
        double getSkyViewFactor(size_t ii_dem, size_t jj_dem, size_t which_triangle);
        void getVectorSun(double solarAzimuth, double solarElevation, Vec3D &v_out);
        double TerrainBiggestDifference(const mio::Array3D<double> &terrain_old, const mio::Array3D<double> &terrain_new);
 
        // Standard Vector operations
        double NormOfVector(const Vec3D &vec1);
        void normalizeVector(const Vec3D &vec1, Vec3D &v_out);
        double VectorScalarProduct(const Vec3D &vec1, const Vec3D &vec2);
        void VectorCrossProduct(const Vec3D &vec1, const Vec3D &vec2, Vec3D &v_out);
        void VectorSum(const Vec3D &vec1, const Vec3D &vec2, Vec3D &v_out);
        void VectorDifference(const Vec3D &vec1, const Vec3D &vec2, Vec3D &v_out);
        void VectorStretch(const Vec3D &vec1, double factor, Vec3D &v_out);
        void RotN(const Vec3D &axis, const Vec3D &vec_in, double rad, Vec3D &v_out);
        void ProjectVectorToPlane(const Vec3D &vec1, const Vec3D &plane_normal, Vec3D &v_out);
        double AngleBetween2Vectors(const Vec3D &vec1, const Vec3D &vec2);
 
        // PVP functions
        void readSP();
 
        //
        void resizeArrays();
 
 
        // Output functions
        void PrintProgress(double percentage);
 
        // Variables
        const size_t dimx, dimy;
        size_t dimx_process;
        size_t startx, endx;
        mio::DEMObject dem;
        const mio::Config &cfg;
 
        SnowBRDF BRDFobject;
        SolarPanel SP;
        std::vector<std::vector<double> > pv_points;
 
        mio::Array3D<std::vector<double>> SortList;                       // Used for speedup in Terrain Iterations
        std::vector<Vec3D> BasicSet_Horizontal; // Horizontal Basic Set [MT 2.1.1 Basic Set]
        mio::Array4D<Vec3D> BasicSet_rotated;     // Basic Set rotated in Triangular Pixel Plane [MT 2.1.3 View List, eq. 2.38]
        mio::Array4D<std::vector<double>> ViewList;                       // Stores all information of network between pixels [MT 2.1.3 View List, eq. 2.47]
        mio::Array2D<double> RList;                                                       // List pre-storage of BRDF values
        mio::Array2D<double> albedo_grid;                                         // Albedo value for each square Pixel
 
        mio::Array2D<double> sky_vf_mean;                 // Grid to store view factor
        std::vector<mio::Array2D<double>> sky_vf; // Array to stor grid of view factor for both values of which_triangle
        unsigned int M_epsilon;                                   // Number of small circles in Basic Set [MT fig. 2.1]
        unsigned int M_phi;                                               // Number of vectors per small circle of Basic Set [MT fig. 2.1]
        unsigned int S;                                                   // Number of vectors per Basic Set [MT fig. 2.1]
        double delta_F_max;                                               // Stopping Treshold for Iteration in W/m2 [MT eq. 2.100]
 
        // Keys from io-file
        bool if_anisotropy = false;             // Anisotropic or Isotropic Snow Model ?
        bool if_multiple = false;               // Do Multiple Scattering in Terrain or not ?
        bool if_write_view_list = true; // Write ViewList to file ?
        bool if_read_view_list = false; // Read existing View-list file? -> Speeds up Initialisation by factor ~200
};
 
#endif