Source code: C++ Solar Irradiation and Daily Light Integral Classes đź› 

C++ source code for calculating solar irradiation and daily light integrals for your location. This class builds upon and utilizes the previously posted classes for the calculation of solar position, found at

and at

This is the first pass version that utilizes basic models for calculating the estimated irradiation for a particular location, time, and date. Right now, it’s just for solar application. As you can imagine, there are quite a few possibilities here including the interrogation of on-line weather data for the calculation of diffuse and direct (beam) radiation, the utilization of real-time spectrometer and pyranometer data, inclusion of supplemental or sole sourced lighting fixtures (as opposed to the sun), etc, etc… This initial version provides basic but ballpark accurate models for “clear sky” days as a starting point. A follow-on class will utilize this class along with the VPD classes to model plant energy balance, transpiration, etc. But, this class, as with the other classes, can be used for whatever ideas that you might conceive.

This is for the DIY types that have some software background. Questions/comment are great but I’d like to ask forum members to avoid the “what’s this good for” type questions within this thread to avoid clutter.

//
// Daily_Light.cpp
// Version timestamp: 9-30-2018, 7:42 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# include "Daily_Light.h"

using namespace std::chrono;
using namespace date;
Daily_Light::Daily_Light()
{

}

Daily_Light::~Daily_Light()
{
    
}

double Daily_Light::calculate_dli(zoned_sun_time time, sun *sun_element)
{
	calculate_dli(time, sun_element, ATOM_SIMPLE);
}


double Daily_Light::calculate_dli(zoned_sun_time time, sun *sun_element, unsigned int atmosphere)
{
    
	sun sun_time;
    std::pair<double, double> element;
    boost::container::flat_set<std::pair<double, double>> the_set;
    unsigned int day_start = (unsigned int)(floor(sun_element->sunrise * 24.0));
    unsigned int day_end = (unsigned int)(floor(sun_element->sunset * 24.0));


    /* Skip over hours where there is no daylight locally, we are not calculating the moonlight contribution. */
	auto local_time = time.get_local_time() + hours(day_start);
	time = date::make_zoned(time.get_time_zone(), local_time);
    
    for (unsigned int hour = day_start;hour <= day_end;hour++)
    {
        
	    sun_time(sun_element->name, sun_element->observer_latitude, sun_element->observer_longitude, time, false);

	    double irrad = calculate_irrad(time, &sun_time, atmosphere);
	    if (((hour == day_start) && (irrad == 0.0)) || (hour == day_end))
	    {
		    /* Find the first and last daylight in the hour. */
		    zoned_sun_time temp_time = time;
		    unsigned int minute_stride = 5;
		    for (unsigned int minute = 1; minute <= 59; minute += minute_stride)
		    {
			    temp_time = time.get_local_time() + minutes(minute);
			    
			    sun_time(sun_element->name, sun_element->observer_latitude, sun_element->observer_longitude, temp_time, false);
			    irrad = calculate_irrad(temp_time, &sun_time, atmosphere);
		    
			    if (irrad > 0.0)
			    {
				    element = { (1.0*hour + minute / MINUTES_HOUR), irrad };
				    the_set.insert(element);
			    }
			    else if (hour == day_end)
			    {
				    element = { (1.0*hour + minute / MINUTES_HOUR), irrad };
				    the_set.insert(element);
				    break;
			    }
		    }
	    }
	    else
	    {
		    element = { hour, irrad };
		    the_set.insert(element);
	    }

	    time = time.get_local_time() + hours(1);
    }
    
    integrate integral(&the_set);
	
	/* SECONDS_HOUR converts the solar constant from unit J / m^2 s (or W m^2) to J m^2 hour */
	this->irrad = SOLAR_CONSTANT * ((integral.integral) * (SECONDS_HOUR)) / 1000000.0; /* MJ/m^2 */
	this->ppfd  = 4.6242 * SOLAR_CONSTANT * ((integral.integral) * (SECONDS_HOUR)) / 1000000.0;
		
    return (irrad);
}

/* Extraterrestrial radiation at 0 degree angle to incident sun radiation (normal plane) */
double Daily_Light::calculate_extraterestrial_radiation(double rad_vector)
{
	double ic_prime = 1 / (rad_vector * rad_vector);
    return (ic_prime);
}


/* Extraterrestrial radiation at surface angle (horizontal plane) to inbound radiation */
double Daily_Light::calculate_incident_surface_irradiation(double normal_radiation, double solar_elevation)
{
	double irrad;
	if (solar_elevation > 0.0)
	{
		double solar_zenith = 90.0 - solar_elevation;
		irrad = /* normal_radiation * */fabs(TO_RADIAN(calculate_plane_surface_angle(solar_zenith)));
	}
	else
	{
		irrad = 0;	
	}
	
    return (irrad);
}


/* Projection */
/* surface_tilt = surface tilt angle from horizontal. */
/* surface_azimuth = surface azimuth angle from the north. */
/* solar_zenith = apparent solar zenith of the sun/light source. */
/* solar_azimuth = apparent solar azimuth of the sun/light source. */
double Daily_Light::calculate_projection(double surface_tilt, double solar_zenith, double surface_azimuth, double solar_azimuth)
{
    double angle = cos(TO_RADIAN(surface_tilt)) * cos(TO_RADIAN(solar_zenith));
    angle += sin(TO_RADIAN(surface_tilt)) * sin(TO_RADIAN(solar_zenith));
    angle *= cos(TO_RADIAN(solar_azimuth - surface_azimuth));
    angle = TO_DEGREE(angle);
    
    return (angle);
}


/* Incident plane angle to inbound radiation */
/* Incident radiation is on a plane that is tangent to the outer surface of the atmosphere. E.g. laying flat on the surface. */
double Daily_Light::calculate_plane_surface_angle(double solar_zenith)
{
    double angle = TO_DEGREE(cos(TO_RADIAN(solar_zenith)) );
    return (angle);
}

/* Incident angle to inbound radiation */
/* surface_tilt = surface tilt angle from horizontal. */
/* surface_azimuth = surface azimuth angle from the north. */
/* solar_zenith = apparent solar zenith of the sun/light source. */
/* solar_azimuth = apparent solar azimuth of the sun/light source. */
double Daily_Light::calculate_incident_angle(double surface_tilt, double solar_zenith, double surface_azimuth, double solar_azimuth)
{
    double angle = calculate_projection(surface_tilt, solar_zenith, surface_azimuth, solar_azimuth);
    angle = TO_DEGREE(acos(TO_RADIAN(angle)));
    
    return (angle);
}


/* Clear day beam transmittance. */
/* The solar radiation received from the sun without having been scattered by the atmosphere. */
double Daily_Light::calculate_beam_transmittance_clear(double altitude_km, std::tuple<double, double, double> climate, double solar_elevation)
{
    double r0 = std::get<0>(climate);
    double r1 = std::get<1>(climate);
    double rk = std::get<2>(climate);
    
    /*  Hottel, ASHRAE 1985  */
	double a0 = r0 * (0.4237 - 0.00821 * (6.0 - altitude_km) * (6.0 - altitude_km));
	double a1 = r1 * (0.5055 - 0.00595 * (6.5 - altitude_km) * (6.5 - altitude_km));
	double k = (-1.0) * rk * (0.2711 - 0.01858 * (2.5 - altitude_km) * (2.5 - altitude_km));
   
    
    double tb = 0.0;
	double zenith = 90.0 - solar_elevation;
    
	if (solar_elevation > 0.0)
	{
		tb = a0 + a1 * exp(k / (TO_DEGREE(sin(TO_RADIAN(fabs(solar_elevation)))) + 0.0000001));
	}
	else
	{
		tb = 0.0;
	}

    return (tb);

}

/* The solar radiation received from the sun having been scattered by the atmosphere. */
double Daily_Light::calculate_diffuse_radiation(double beam_transmittance)
{
	double td = 0.0;
	
	/*  Liu and Jordan (1960)  */
	if (beam_transmittance > 0)
	{
		td = 0.2710 - (0.2939 * beam_transmittance);
	}
    
    return (td);
}

double Daily_Light::calculate_irrad(zoned_sun_time time, sun *sun_element, unsigned int atmosphere)
{
	double solar_elevation = sun_element->solar_elevation_corrected;
    std::tuple<double, double, double> climate = std::tuple<double, double, double>(0.97, 0.99, 1.02);
    
    this->ic_prime = calculate_extraterestrial_radiation(sun_element->rad_vector);
	this->irrad_tilted = calculate_incident_surface_irradiation(this->ic_prime, solar_elevation);
	switch (atmosphere)
	{
	case (ATOM_NONE):
		{
			this->direct_normal = 1.0;
			this->beam_transmittance = 1.0;
			this->diffuse_radiation = 0.0;
		
			/* Total solar radiation. */
			this->irrad = this->irrad_tilted * (this->beam_transmittance + this->diffuse_radiation);
			
			break;
		}
	case (ATOM_SIMPLE):
		{
			this->relative_optical_path = optical_air_mass(solar_elevation, sun_element->sec_z_inverse, 0.0, AM_PICKERING);
			this->beam_transmittance = calculate_beam_transmittance_clear(0.0, climate, solar_elevation);
			this->diffuse_radiation = calculate_diffuse_radiation(this->beam_transmittance);
			
			this->irrad = this->irrad_tilted * (this->beam_transmittance + this->diffuse_radiation);
			
			break;
		}
	case (ATOM_BIRD):
		{
			this->relative_optical_path = optical_air_mass(solar_elevation, sun_element->sec_z_inverse, 0.0, AM_PICKERING);
			bird bird_calculation(solar_elevation, ic_prime, this->relative_optical_path, TO_MILLIBAR(ATM1_AVE), 0.3, 1.5, 0.15, 0.1, 0.85, 0.076370, 0.2);
			this->global_horizontal = bird_calculation.global_horizontal;
			this->irrad = this->global_horizontal;
			
			break;
		}
	case (ATOM_PEREZ):
		{
			this->irrad = 0.0;
			break;
		}
	case (ATOM_INEICHEN):
		{
			this->irrad = 0.0;
			break;
		}
	default:
		{
			this->direct_normal = 1.0;
			this->beam_transmittance = 1.0;
			this->diffuse_radiation = 0.0;
		
			/* Total solar radiation. */
			this->irrad = this->irrad_tilted * (this->beam_transmittance + this->diffuse_radiation);
			break;
		}
	}
	
	
    return (this->irrad);
}

/* Air mass, relative optical path. Zenith is apparent zenith. */
double Daily_Light::optical_air_mass(double solar_elevation, double inverse_secant_zenith, double altitude, unsigned int model)
{
    double relative_optical_path;
    double altitude_h = 1.0;
	double zenith = 90.0 - solar_elevation;
    
    if (altitude > 0.0)
    {
        altitude_h = exp(-0.0001184 * altitude); 
    }
    
    switch (model)
    {
    case AM_KASTENYOUNG:
        {
            double power = 0.50572*pow((6.07995 + (90 - zenith)), -1.6364);
            relative_optical_path = altitude_h / (cos(TO_RADIAN(zenith)) + power);
            break;
        }	
    case AM_KASTEN:
        {
            double power = 0.15*pow((93.885 - zenith), -1.253);
            relative_optical_path = altitude_h / (cos(TO_RADIAN(zenith)) + power);
            break;
        }	
    case AM_PICKERING:
        {
            double power = pow(165.0 + 47.0 * (90.0 - zenith), 1.1);
            relative_optical_path = (altitude_h / (sin(TO_RADIAN(90.0 - zenith + 244.0 / power))));
            break;
        }
    case AM_GUEYMARD:
        {
            double power = zenith * 0.00176759 * pow((94.37515 - zenith), -1.21563);
            relative_optical_path = altitude_h / (cos(TO_RADIAN(zenith)) + power);
            break;
            break;
        }
    default:
        break;
    }
	
	if (isnanf(relative_optical_path))
	{
		relative_optical_path = 0.0;
	}

    return (relative_optical_path);
}

/* Air mass, absolute air mass at input atmospheric pressure. Return is in pascals */
double Daily_Light::absolute_air_mass(double relative_optical_path, double air_pressure)
{
    double absolute_airmass = (relative_optical_path*air_pressure) / (ATM1_AVE);
    
    return (absolute_airmass);
}

/* Air mass, absolute air mass at sea level. Return is in pascals */
double Daily_Light::absolute_air_mass(double relative_optical_path)
{
    double absolute_airmass = (relative_optical_path);
    
    return (absolute_airmass);
}
//
// Daily_Light.h
// Version timestamp: 9-30-2018, 7:42 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# pragma once
# include "integrate.h"
# include "sun.h"
# include "bird.h"

enum
{
	AM_HARDIE,
	AM_KASTENYOUNG,
	AM_KASTEN,
	AM_PICKERING,
	AM_GUEYMARD
};

enum 
{
	ATOM_NONE,
        ATOM_SIMPLE,
	ATOM_BIRD,
	ATOM_PEREZ,
	ATOM_INEICHEN
};

class Daily_Light : public integrate, sun
{
public:
	Daily_Light();
	~Daily_Light();
	double calculate_dli(zoned_sun_time time, sun *sun_element);
	double calculate_dli(zoned_sun_time time, sun *sun_element, unsigned int atmosphere);
	double calculate_irrad(zoned_sun_time time, sun *sun_element, unsigned int atmosphere);
	double optical_air_mass(double zenith, double atmospheric_pressue);
	double calculate_extraterestrial_radiation(double rad_vector);
	double calculate_incident_surface_irradiation(double normal_radiation, double solar_elevation);
	double calculate_plane_surface_angle(double solar_zenith);
	double calculate_projection(double surface_tilt, double solar_zenith, double surface_azimuth, double solar_azimuth);
	double calculate_incident_angle(double surface_tilt, double solar_zenith, double surface_azimuth, double solar_azimuth);

	double calculate_beam_transmittance_clear(double optical_air_mass, std::tuple<double, double, double> climate, double solar_elevation);
	double calculate_diffuse_radiation(double beam_transmittance);
	double optical_air_mass(double solar_elevation, double inverse_secant_zenith, double altitude, unsigned int model);
	double absolute_air_mass(double relative_optical_path);
	double absolute_air_mass(double relative_optical_path, double air_pressure);

	double ic_prime;
	double relative_optical_path;
	double beam_transmittance;
	double diffuse_radiation;
	double irrad_tilted;
	double irrad;
	double ppfd;
	
	double direct_normal;
	double direct_horizontal;
	double global_horizontal;
};
//
// integrate.cpp
// Version timestamp: 9-26-2018, 10:56 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# include "integrate.h"

using namespace boost::container;

integrate::integrate()
{
}

/* In this case, a user supplies a set of unique (x,y) values to integrate. */
integrate::integrate(flat_set<std::pair<double, double>> *by_point)
{
	integrate_by_point(by_point);
}

integrate::~integrate()
{
}

void integrate::integrate_by_point(flat_set<std::pair<double, double>> *by_point)
{
	/* Integrate the current set of values */
	if (by_point->size() > 1)
	{
		bool first = true;
		
		double x1, x2, y1, y2;
		for (std::pair<double, double> element : *by_point)
		{
			if (first)
			{
				first = false;
				x2 = element.first;
				y2 = element.second;
			}
			else
			{
				x1 = x2;
				y1 = y2;
				x2 = element.first;
				y2 = element.second;
				running_trapezodal_integral(x1, y1, x2, y2);
			}

		}
	}
	else
	{
		this->integral = 0.0;
	}
}

void integrate::operator()(flat_set<std::pair<double, double>> *by_point)
{
	this->integral = 0.0;
	integrate_by_point(by_point);
}

void integrate::clear()
{
	this->integral = 0.0;
}

void integrate::running_trapezodal_integral(double x1, double y1, double x2, double y2)
{
	/* Trapazodal integration. */
	double dt = (x2 - x1);
	this->integral += ((y2 + y1)*dt) * 0.5;
}
//
// integrate.h
// Version timestamp: 9-26-2018, 10:56 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# pragma once

# include <cstdlib>
# include <iostream>
# include <sstream>
# include <boost/container/flat_set.hpp>

class integrate
{
public:
	integrate();
	integrate(boost::container::flat_set<std::pair<double, double>> *by_point);
	~integrate();
	void clear();
	void integrate_by_point(boost::container::flat_set<std::pair<double, double>> *by_point);
	void running_trapezodal_integral(double x1, double y1, double x2, double y2);
	void operator()(boost::container::flat_set<std::pair<double, double>> *by_point);
	
	double integral;
};

CC BY-SA 4.0

3 Likes

Example


*********************** Extraterrestrial radiation normal to solar radiation ***************************************


The following example routine utilizes the after-mentioned classes to determine the solar radiation on a plane normal to the sun at the surface of the earths atmosphere (or on earth with no atmosphere). This is the starting point to calculating the solar radiation after it passes through the air, clouds, reflections, etc…

/* Calculate the extraterrestrial radiation per each day of a year normal to the plane of the sun. */
	std::cout << "********************************************************************************************************************" << std::endl;
	std::cout << "*********************** Extraterrestrial radiation normal to solar radiation ***************************************" << std::endl;
	std::cout << "********************************************************************************************************************" << std::endl;
	do
	{
		Daily_Light dli;
		zoned_sun_time newyork_area;
		sun sun_element_newyork;
		std::list<std::tuple<double, double, zoned_sun_time>> rad_map;
		
		double solar_noon;
		
		double distmax, distmax_copy, distmin_copy, distmin = 5000.0;
		for (int day = 0; day < 365; day++)
		{
			/* Find the day with the maximum and minimum amounts of daylight. */
			newyork_area = date::make_zoned("America/Phoenix", date::local_days{ date::January / 01 / 2018 } + std::chrono::hours(12) + date::days(day));
			sun_element_newyork.set("Pheonix", 33.448376, -112.074036, newyork_area, false);
			solar_noon = sun_element_newyork.solar_noon;
			
			/* Calculate sun position at solar noon. */
			auto time_hms = sun_element_newyork.fp_days_to_chrono(solar_noon);
			auto time_hmsd = date::local_days{ date::January / 01 / 2018 } + time_hms + date::days(day);
			newyork_area = date::make_zoned("America/Phoenix", time_hmsd);
			sun_element_newyork.set("Pheonix", 40.730610, -73.935242, newyork_area, false);
			
			double distance = dli.calculate_extraterestrial_radiation(sun_element_newyork.rad_vector);
			
			/* Save this information. */
			rad_map.push_back(std::tuple<double, double, zoned_sun_time>((double)day, distance, newyork_area));
        
			if (distance > distmax) distmax = distance;
			if (distance < distmin) distmin = distance;
		}
    
		distmax_copy = distmax;
		int distmax_ceil = ceil(distmax);
		distmin_copy = distmin;
		int distmin_floor = floor(distmin);	
		double chart_x_quant = 96.0 / (distmax_ceil - distmin_floor);

		/* Plot the chart */
		int day_stride = 4;
		
		for (std::tuple<double, double, zoned_sun_time> element : rad_map)
		{
			int day = std::get<0>(element);
			double distance = std::get<1>(element);
			newyork_area = std::get<2>(element);		
        
			if ((distance) >= (distmax_copy) || (distance) <= (distmin_copy) || (day % day_stride == 0))
			{

				std::cout << date::format("%b %d %Y", newyork_area) << ": " << std::string(round(((distance) - distmin)*chart_x_quant), '*');
        
				if ((distance) >= (distmax_copy))
				{
					std::cout << " <-max- " << distance  << " W/m^2";
				}
            
				if ((distance) <= (distmin_copy))
				{
					std::cout << " <-min- " << distance << " W/m^2";
				}
        
				std::cout << std::endl;
			}	    
		}
		
		std::cout << std::endl;
	
	} while (0);
	

Produces the following:

CC BY-SA 4.0

2 Likes

Example


*********************** Extraterrestrial radiation tangent to surface **********************************************


	/*************/
	/* Calculate the extraterrestrial radiation per each day of a year tangent to the surface. */
	std::cout << "********************************************************************************************************************" << std::endl;
	std::cout << "*********************** Extraterrestrial radiation tangent to surface **********************************************" << std::endl;
	std::cout << "********************************************************************************************************************" << std::endl;
	do
	{
		Daily_Light dli;
		zoned_sun_time newyork_area;
		sun sun_element_newyork;
		std::list<std::tuple<double, double, zoned_sun_time>> rad_map;
		
		double distmax, distmax_copy, distmin_copy, distmin = 5000.0;
		for (int day = 0; day < 365; day++)
		{
			/* Find the solar noons and the incident radiation max/mins throughout the year. */
			/* Find solar noon for today. */
			newyork_area = date::make_zoned("America/Phoenix", date::local_days{ date::January / 01 / 2018 } + std::chrono::hours(12) + date::days(day));
			sun_element_newyork.set("Pheonix", 33.448376, -112.074036, newyork_area, false);
			
			double solar_noon = sun_element_newyork.solar_noon;
			
			/* Calculate sun position at solar noon. */
			auto time_hms = sun_element_newyork.fp_days_to_chrono(solar_noon);
			auto time_hmsd = date::local_days{ date::January / 01 / 2018 } + time_hms + date::days(day);
			newyork_area = date::make_zoned("America/Phoenix", time_hmsd);
			sun_element_newyork.set("Pheonix", 33.448376, -112.074036, newyork_area, false);
			
			/* Determine radiation. */
			double radiation = dli.calculate_extraterestrial_radiation(sun_element_newyork.rad_vector);
			double distance = dli.calculate_incident_surface_irradiation(radiation, sun_element_newyork.solar_elevation_corrected);
			/* Save this information. */
			rad_map.push_back(std::tuple<double, double, zoned_sun_time>((double)day, distance, newyork_area));
        
			if (distance > distmax) distmax = distance;
			if (distance < distmin) distmin = distance;
		}
    
		distmax_copy = distmax;
		int distmax_ceil = ceil(distmax);
		distmin_copy = distmin;
		int distmin_floor = floor(distmin);	
		double chart_x_quant = 96.0 / (distmax_ceil - distmin_floor);

		/* Plot the chart */
		int day_stride = 4;

		for (std::tuple<double, double, zoned_sun_time> element : rad_map)
		{
			int day = std::get<0>(element);
			double distance = std::get<1>(element);
			newyork_area = std::get<2>(element);

			if ((distance) >= (distmax_copy) || (distance) <= (distmin_copy) || (day % day_stride == 0))
			{

				std::cout << date::format("%b %d %Y %R", newyork_area) << " : " << FG_GREEN << std::string(round(((distance) - distmin)*chart_x_quant), '*') << FG_DEFAULT << " : " << FG_YELLOW << distance << FG_DEFAULT;
				
				if ((distance) >= (distmax_copy))
				{
					std::cout << FG_GREEN << " <-max- W/m^2" << FG_DEFAULT;
				}
            
				if ((distance) <= (distmin_copy))
				{
					std::cout << FG_GREEN << " <-min- W/m^2" << FG_DEFAULT;
				}
				std::cout << std::endl;
			}	    
        
		}
		std::cout << std::endl;
	
	} while (0) ;

Produces the following:

CC BY-SA 4.0

1 Like

Example.

The following example generates graphs for several locations for the solar radiation received at the ground as if there were no atmosphere.


*********************** Min / Max DLI solar radiation days *********************************************************


	/*************/
	/* Calculate the dli on the min/max solar radiation days. */
	do
	{
		zoned_sun_time shanghai1 = date::make_zoned("Asia/Shanghai", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time losangeles1 = date::make_zoned("America/Los_Angeles", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time sydney1 = date::make_zoned("Australia/Sydney", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time newyork1 = date::make_zoned("America/New_York", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time moscow1 = date::make_zoned("Europe/Moscow", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time chicago1 = date::make_zoned("America/Chicago", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time newzealand1 = date::make_zoned("Pacific/Auckland", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time quito1 = date::make_zoned("America/Guayaquil", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time Reykjavik = date::make_zoned("Atlantic/Reykjavik", date::local_days{ date::January / 01 / 2018 });
	
		std::vector<sun> sun_info_vec1;
		sun sun_element1("Reykjavik, Iceland", 64.135666, -21.862675, Reykjavik);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Shanghai", 31.22222, 121.45806, shanghai1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Los Angeles", 34.052235, -118.243683, losangeles1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Sydney", -33.865143, 151.209900, sydney1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("New York", 40.730610, -73.935242, newyork1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Moscow", 55.751244, 37.618423, moscow1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Chicago", 41.881832, -87.623177, chicago1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Auckland,NZ", -36.848461, 174.763336, newzealand1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Quito, Ecuador", -0.22985, -78.52495, quito1);
		sun_info_vec1.push_back(sun_element1);

	
		for (sun sun_element_location : sun_info_vec1)
		{
		
			std::cout << "********************************************************************************************************************" << std::endl;
			std::cout << "*********************** Min / Max DLI solar Radiation days without Atmosphere **************************************" << std::endl;
			std::cout << "********************************************************************************************************************" << std::endl;

			Daily_Light dli;
			std::list<std::tuple<double, double, double, zoned_sun_time>> rad_map;
			zoned_sun_time sun_time = *sun_element_location.zone_time;
		
			std::cout << "Calculating : ";
			double distmax = 0;
			double distmax_copy = 0;
			double distmin_copy = 0;
			double distmin = 5000.0; 
			double julian_last = 0;
			
			for (int day = 0; day < 365; day++)
			{
				if (day % 5 == 0) std::cout << "." << std::flush;
			
				/* Determine radiation. */
				double distance = dli.calculate_dli(sun_time, &sun_element_location, false);
			
				/* Save this information. */
				rad_map.push_back(std::tuple<double, double, double, zoned_sun_time>((double)day, distance, sun_element_location.solar_elevation_corrected, sun_time));
				if (distance > distmax)
				{
					distmax = distance;
				}
				if (distance < distmin) 
				{
					distmin = distance;
				}
				
				/* Update the time. */
				sun_time = sun_time.get_local_time() + date::days{ 1 };
				sun_element_location.set(sun_element_location.name, sun_element_location.observer_latitude, sun_element_location.observer_longitude, sun_time);
			}
			std::cout << std::endl;
			std::cout << "Location : " << sun_element_location.name << std::endl;
    
			distmax_copy = distmax;
			int distmax_ceil = ceil(distmax);
			distmin_copy = distmin;
			int distmin_floor = floor(distmin);	
			double chart_x_quant = 96.0 / (distmax_ceil - distmin_floor);

		
			/* Plot the chart */
			int day_stride = 3;
			double yearly_DLI = 0.0; julian_last = 0;
			for (std::tuple<double, double, double, zoned_sun_time> element : rad_map)
			{
				int day = std::get<0>(element);
				double distance = std::get<1>(element);
				double temp = std::get<2>(element);
				zoned_sun_time sun_time_location = std::get<3>(element);
				yearly_DLI += distance;
			
			
				if ((distance) >= (distmax_copy) || (distance) <= (distmin_copy) || (day % day_stride == 0))
				{

					std::cout << date::format("%b %d %Y", sun_time_location) << " : " << FG_GREEN << std::string(round(((distance) - distmin)*chart_x_quant), '*') << FG_DEFAULT << " : " << FG_YELLOW << distance << FG_DEFAULT;
				
					if ((distance) >= (distmax_copy))
					{
						std::cout << FG_GREEN << " <-max- DLI (MJ / m^2 day)" << FG_DEFAULT;
					}
            
					if ((distance) <= (distmin_copy))
					{
						std::cout << FG_GREEN << " <-min- DLI (MJ / m^2 day)" << FG_DEFAULT;
					}
					std::cout << std::endl;
				}	    
        
			}
			std::cout << std::endl;	
			std::cout << "Yearly DLI : " << yearly_DLI << " (MJ / m^2 year)" << std::endl;
			rad_map.clear();
		}
	
	} while (0);

Produces the following:

and, so on…

In a similar fashion, estimating the DLI taking into account direct beam radiation, diffuse radiation, and the effect of the atmosphere:

	/*************/
	/* Calculate the dli on the min/max solar radiation days. */
	do
	{
		zoned_sun_time shanghai1 = date::make_zoned("Asia/Shanghai", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time losangeles1 = date::make_zoned("America/Los_Angeles", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time sydney1 = date::make_zoned("Australia/Sydney", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time newyork1 = date::make_zoned("America/New_York", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time moscow1 = date::make_zoned("Europe/Moscow", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time chicago1 = date::make_zoned("America/Chicago", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time newzealand1 = date::make_zoned("Pacific/Auckland", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time quito1 = date::make_zoned("America/Guayaquil", date::local_days{ date::January / 01 / 2018 });
		zoned_sun_time Reykjavik = date::make_zoned("GMT", date::local_days{ date::January / 01 / 2018 });
	
		std::vector<sun> sun_info_vec1;
		sun sun_element1("Shanghai", 31.22222, 121.45806, shanghai1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Los Angeles", 34.052235, -118.243683, losangeles1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Sydney", -33.865143, 151.209900, sydney1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("New York", 40.730610, -73.935242, newyork1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Moscow", 55.751244, 37.618423, moscow1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Chicago", 41.881832, -87.623177, chicago1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Auckland,NZ", -36.848461, 174.763336, newzealand1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Quito, Ecuador", -0.22985, -78.52495, quito1);
		sun_info_vec1.push_back(sun_element1);
		sun_element1.set("Reykjavik, Iceland", 64.135666, -21.862675, Reykjavik);
		sun_info_vec1.push_back(sun_element1);
	
		for (sun sun_element_location : sun_info_vec1)
		{
		
			std::cout << "********************************************************************************************************************" << std::endl;
			std::cout << "*********************** Min / Max DLI solar Radiation days with Atmosphere *****************************************" << std::endl;
			std::cout << "********************************************************************************************************************" << std::endl;

			Daily_Light dli;
			std::list<std::tuple<double, double, double, zoned_sun_time>> rad_map;
			zoned_sun_time sun_time = *sun_element_location.zone_time;
		
			std::cout << "Calculating : ";
			double distmax = 0;
			double distmax_copy = 0;
			double distmin_copy = 0;
			double distmin = 5000.0; 
			double julian_last = 0;
			
			for (int day = 0; day < 365; day++)
			{
				if (day % 5 == 0) std::cout << "." << std::flush;
			
				/* Determine radiation. */
				double distance = dli.calculate_dli(sun_time, &sun_element_location, true);
			
				/* Save this information. */
				rad_map.push_back(std::tuple<double, double, double, zoned_sun_time>((double)day, distance, sun_element_location.solar_elevation_corrected, sun_time));
				if (distance > distmax)
				{
					distmax = distance;
				}
				if (distance < distmin) 
				{
					distmin = distance;
				}
				
				/* Update the time. */
				sun_time = sun_time.get_local_time() + date::days{ 1 };
				sun_element_location.set(sun_element_location.name, sun_element_location.observer_latitude, sun_element_location.observer_longitude, sun_time);
			}
			std::cout << std::endl;
			std::cout << "Location : " << sun_element_location.name << std::endl;
    
			distmax_copy = distmax;
			int distmax_ceil = ceil(distmax);
			distmin_copy = distmin;
			int distmin_floor = floor(distmin);	
			double chart_x_quant = 96.0 / (distmax_ceil - distmin_floor);

		
			/* Plot the chart */
			int day_stride = 5;
			double yearly_DLI = 0.0; julian_last = 0;
			for (std::tuple<double, double, double, zoned_sun_time> element : rad_map)
			{
				int day = std::get<0>(element);
				double distance = std::get<1>(element);
				double temp = std::get<2>(element);
				zoned_sun_time sun_time_location = std::get<3>(element);
				yearly_DLI += distance;
			
			
				if ((distance) >= (distmax_copy) || (distance) <= (distmin_copy) || (day % day_stride == 0))
				{

					std::cout << date::format("%b %d %Y", sun_time_location) << " : " << FG_GREEN << std::string(round(((distance) - distmin)*chart_x_quant), '*') << FG_DEFAULT << " : " << FG_YELLOW << distance << FG_DEFAULT;
				
					if ((distance) >= (distmax_copy))
					{
						std::cout << FG_GREEN << " <-max- DLI (MJ / m^2 day)" << FG_DEFAULT;
					}
            
					if ((distance) <= (distmin_copy))
					{
						std::cout << FG_GREEN << " <-min- DLI (MJ / m^2 day)" << FG_DEFAULT;
					}
					std::cout << std::endl;
				}	    
        
			}
			std::cout << std::endl;	
			std::cout << "Yearly DLI : " << yearly_DLI << " (MJ / m^2 year)" << std::endl;
			rad_map.clear();
		}
	
	} while (0);

CC BY-SA 4.0

2 Likes

Example

The following example generates a map of DLI by latitude for the solar radiation received at the ground as if there were no atmosphere. Atmosphere and weather affect the actual DLI received, the above class, currently, can also calculate the DLI with a first order estimation of attenuation.


******************************************* DLI by latitude and year day ***********************************************


		
	/*************/
	/* Calculate the dli on the min/max solar radiation days. */
	std::cout << "****************************************************************************************************************************************" << std::endl;
	std::cout << "******************************************* DLI by latitude and year day **************************************************************" << std::endl;
	std::cout << "****************************************************************************************************************************************" << std::endl;
	do
	{
		Daily_Light dli;
		zoned_sun_time location;
		sun sun_element_newyork;
		std::list<std::tuple<double, double, double>> rad_map;
		
		int day_stride = 3;
		std::cout << "LAT     1" <<  std::string((round(365 / day_stride) - 10) / 2, ' ') << "day" <<  std::string((round(365 / day_stride) - 1) / 2, ' ') << "365   Total(yr)" << std::endl;
		
		for (int calc_latitude = 84; calc_latitude >= -84; calc_latitude -= 2)
		{
			double total_irrad = 0.0;
			
			for (int day = 0; day < 365; day++)
			{
				/* Find the solar noons and the incident radiation max/mins throughout the year. */
				/* Find solar noon for today. */
				location = date::make_zoned("UTC", date::local_days{ date::January / 01 / 2018 } + std::chrono::hours(0) + date::days(day));
				sun_element_newyork.set("UTC", calc_latitude, -73.935242, location);
			
				/* Determine radiation. */
				double distance = dli.calculate_dli(location, &sun_element_newyork, false);
			
				/* Save this information. */
				rad_map.push_back(std::tuple<double, double, double>((double)day, distance, calc_latitude));
			}
		
			/* Plot the chart */

			for (std::tuple<double, double, double> element : rad_map)
			{
				int day = std::get<0>(element);
				double distance = std::get<1>(element);
				double calc_latitude = std::get<2>(element);
				
				if (day == 0) std::cout << boost::format("%-5.1f")  % calc_latitude << " : ";
				
				total_irrad += distance;
				
				if ((day % day_stride == 0))
				{

					if ((distance) >= (40))
					{
						std::cout << BG_RED << " " << FG_DEFAULT << BG_DEFAULT;
					}	
					else if ((distance) >= (35))
					{
						std::cout << BG_MAGENTA << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (30))
					{
						std::cout << BG_ORANGE << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (25))
					{
						std::cout << BG_YELLOW << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (20))
					{
						std::cout << BG_LIGHT_YELLOW  << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (15))
					{
						std::cout << BG_LIGHT_GREEN  << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (10))
					{
						std::cout << BG_GREEN << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (5))
					{
						std::cout << BG_CYAN << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (3))
					{
						std::cout << BG_LIGHT_BLUE  << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (2))
					{
						std::cout << BG_BLUE << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else if ((distance) >= (0))
					{
						std::cout << BG_LIGHT_BLACK << " " << FG_DEFAULT << BG_DEFAULT;
					}
					else
					{
						std::cout << BG_DARK_GRAY << " " << FG_DEFAULT << BG_DEFAULT;
					}
				}    
			} 
			std::cout << " : " << total_irrad << std::endl;
			rad_map.clear();
			
		} /* end calc_longitude */

		std::cout << "Legend:" << std::endl;	
		std::cout << "MJ/m^2 day :" << FG_RED << "* >40 \t" << FG_MAGENTA << "* >35 \t" << FG_ORANGE << "* >30 \t" << FG_YELLOW << "* >25 \t" << FG_LIGHT_YELLOW << "* >20 \t" << FG_DEFAULT << BG_DEFAULT << std::endl;
		std::cout << "MJ/m^2 day :" << FG_LIGHT_GREEN << "* >15 \t" << FG_GREEN << "* >10 \t" << FG_CYAN << "* >5 \t" << FG_LIGHT_BLUE << "* >3 \t" << FG_LIGHT_BLACK << "* >=0 \t" << FG_DEFAULT << BG_DEFAULT << std::endl;
		std::cout << std::endl;	
	} while (0) ;

Produces the following:

CC BY-SA 4.0

4 Likes

The following example routine shows how to utilize the after-mentioned classes to determine the daily light integral (DLI) for a specific location and date:

	/*************/
	/* DLI tests */
	do
	{
		/* Estimated as 33.8MJ/m^2 day */
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::April / 15 / 2018 } + std::chrono::hours(0));
		sun_element_location.set("Boston", 43.0, -71, solar_date, true);
		
		double daily_light = dli.calculate_dli(solar_date, &sun_element_location, false);
		
		std::cout << "DLI at latitude 43 north on 4/15 : " << daily_light << " kW/m^2 day" << std::endl ;
		
		daily_light = dli.calculate_dli(solar_date, &sun_element_location, true);
		std::cout << "DLI with atmosphere at latitude 43 north (" << sun_element_location.name << ") on 4/15 : " << daily_light << " kW/m^2 day" << std::endl;
		
	} while (0);
	
	do
	{
		/* Estimated as 33.8MJ/m^2 day */
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::June / 15 / 2018 } + std::chrono::hours(0));
		sun_element_location.set("Boston", 43.0, -71, solar_date, true);
		
		double daily_light = dli.calculate_dli(solar_date, &sun_element_location, false);
		
		std::cout << "DLI at latitude 43 north on 6/15 : " << daily_light << " kW/m^2 day" << std::endl;
		
		daily_light = dli.calculate_dli(solar_date, &sun_element_location, true);
		std::cout << "DLI with atmosphere at latitude 43 north (" << sun_element_location.name << ") on 6/15 : " << daily_light << " kW/m^2 day" << std::endl;
		
	} while (0);
	
	do
	{
		/* Estimated as 32.2MJ/m^2 day */
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::September / 3 / 2018 } + std::chrono::hours(0));
		sun_element_location.set("Unknown", -20.0, -112.074036, solar_date, true);
		
		double daily_light = dli.calculate_dli(solar_date, &sun_element_location, false);
		
		std::cout << "DLI at latitude -20 north on 9/3 : " << daily_light << " kW/m^2 day" << std::endl;
		
		daily_light = dli.calculate_dli(solar_date, &sun_element_location, true);
		std::cout << "DLI with atmosphere at latitude -20 (" << sun_element_location.name << ") on 9/3 : " << daily_light << " kW/m^2 day" << std::endl;
		
	} while (0);
	
	
	do
	{
		/* Estimated as 25.1MJ/m^2 day */
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		
		solar_date = date::make_zoned("Brazil/East", date::local_days{ date::May / 15 / 2018 } + std::chrono::hours(0));
		sun_element_location.set("Rio de Janeiro", -22.9032, -43.1729, solar_date, true);
		
		double daily_light = dli.calculate_dli(solar_date, &sun_element_location, false);
		std::cout << "DLI at latitude -22.9032 north (" << sun_element_location.name << ") on 5/15 : " << daily_light << " kW/m^2 day" << std::endl;
		
		daily_light = dli.calculate_dli(solar_date, &sun_element_location, true);
		std::cout << "DLI with atmosphere at latitude -22.9032 (" << sun_element_location.name << ") on 5/15 : " << daily_light << " kW/m^2 day" << std::endl;
		
	} while (0);

Generates the following output:

DLI at latitude 43 north on 4/15 : 34.2952 kW/m^2 day
DLI with atmosphere at latitude 43 north (Boston) on 4/15 : 12.8464 kW/m^2 day
DLI at latitude 43 north on 6/15 : 42.6516 kW/m^2 day
DLI with atmosphere at latitude 43 north (Boston) on 6/15 : 15.9811 kW/m^2 day
DLI at latitude -20 north on 9/3 : 32.3038 kW/m^2 day
DLI with atmosphere at latitude -20 (Unknown) on 9/3 : 12.1015 kW/m^2 day
DLI at latitude -22.9032 north (Rio de Janeiro) on 5/15 : 25.4564 kW/m^2 day
DLI with atmosphere at latitude -22.9032 (Rio de Janeiro) on 5/15 : 9.53124 kW/m^2 day

Likewise, if one wanted to estimate the solar energy at a specific hour in a day at a specific location:

	do
	{
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		boost::container::flat_set<std::pair<double, double>> the_set;
		std::pair<double, double> element;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::April / 15 / 2018 } + std::chrono::hours(10));
		sun_element_location.set("Boston", 43.0, -112.074036, solar_date, true);
		
		double hour1_light = dli.calculate_irrad(solar_date, &sun_element_location, false);
		double irrad = (hour1_light); 
		
		std::cout << "Quantum Flux (W/m^2/hr) at 10 at latitude 43 north (" << sun_element_location.name << ") on 4/15 : " << irrad << std::endl;
		
	} while (0);
	
	do
	{
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		boost::container::flat_set<std::pair<double, double>> the_set;
		std::pair<double, double> element;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::April / 15 / 2018 } + std::chrono::hours(11));
		sun_element_location.set("Boston", 43.0, -112.074036, solar_date, true);
		
		double hour1_light = dli.calculate_irrad(solar_date, &sun_element_location, false);
		double irrad = (hour1_light); 
		
		std::cout << "Quantum Flux (W/m^2/hr) at 11 at latitude 43 north (" << sun_element_location.name << ") on 4/15 : " << irrad << std::endl;
		
	} while (0);
	
	do
	{
		Daily_Light dli;
		zoned_sun_time solar_date;
		sun sun_element_location;
		boost::container::flat_set<std::pair<double, double>> the_set;
		std::pair<double, double> element;
		
		solar_date = date::make_zoned("America/New_York", date::local_days{ date::April / 15 / 2018 } + std::chrono::hours(12));
		sun_element_location.set("Boston", 43.0, -112.074036, solar_date, true);
		
		double hour1_light = dli.calculate_irrad(solar_date, &sun_element_location, false);
		double irrad = (hour1_light); 
		
		std::cout << "Quantum Flux (W/m^2/hr)  at 12 at latitude 43 north (" << sun_element_location.name << ") on 4/15 : " << irrad << std::endl;
		
	} while (0);

Produces the following:

Quantum Flux (W/m^2/hr) at 10 at latitude 43 north (Boston) on 4/15 : 542.537
Quantum Flux (W/m^2/hr) at 11 at latitude 43 north (Boston) on 4/15 : 763.418
Quantum Flux (W/m^2/hr) at 12 at latitude 43 north (Boston) on 4/15 : 943.306
4 Likes

The following class, in conjunction with the class in the OP, is the “Bird” atmospheric model for solar insolation. This model provides a greater degree of accuracy than the simple model contained within the “Daily_Light” class.

The Bird model is based on the paper “A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces”, Bird 1991.
It can be retrieved from: http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf

Additional models will be developed as time allows.

//
// bird.cpp
// Version timestamp: 9-30-2018, 7:31 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# include "bird.h"

bird::bird()
{
}

bird::bird(double solar_elevation, double ic_prime, double relative_optical_path, double air_pressure, double ozone_cm, double h2o_cm, double depth_380nm, double depth_500nm, double Ba, double tau_aerosol, double albedo)
{
	this->solar_elevation = solar_elevation;
	this->ic_prime = ic_prime;
	this->relative_optical_path = relative_optical_path;
	this->air_pressure = air_pressure;
	this->ozone_cm = ozone_cm;
	this->h2o_cm = h2o_cm;
	this->depth_380nm = depth_380nm;
	this->depth_500nm = depth_500nm;
	this->Ba = Ba;
	this->tau_aerosol = tau_aerosol;
	this->albedo = albedo;
	
	calculate_bird();
}

bird::~bird()
{
}

/* Air mass, absolute air mass at input atmospheric pressure. Return is in pascals */
double bird::absolute_air_mass(double relative_optical_path, double air_pressure)
{
	double absolute_airmass = (relative_optical_path*air_pressure) / (ATM1_AVE);
    
	return (absolute_airmass);
}

/* Air mass, absolute air mass at sea level. Return is in pascals */
double bird::absolute_air_mass(double relative_optical_path)
{
	double absolute_airmass = (relative_optical_path);
    
	return (absolute_airmass);
}

/* Incident plane angle to inbound radiation */
/* Incident radiation is on a plane that is tangent to the outer surface of the atmosphere. E.g. laying flat on the surface. */
double bird::calculate_plane_surface_angle(double solar_zenith)
{
	double angle = TO_DEGREE(cos(TO_RADIAN(solar_zenith)));
	return (angle);
}

/* Extraterrestrial radiation at surface angle (horizontal plane) to inbound radiation */
double bird::calculate_incident_surface_irradiation(double solar_elevation)
{
	double irrad;
	if (solar_elevation > 0.0)
	{
		double solar_zenith = 90.0 - solar_elevation;
		irrad = fabs(TO_RADIAN(calculate_plane_surface_angle(solar_zenith)));
	}
	else
	{
		irrad = 0;	
	}
	
	return (irrad);
}

/* Estimate broadband aerosol optical depth from depth at 380nm and 500nm. */
/* Bird and Hulstrom, "Direct Insolation Models" (1980) */
/* http://www.nrel.gov/docs/legosti/old/344.pdf */
/* Outputs range from 0.02 - 0.5 */
double bird::calculate_aerosol_taua(double depth_380nm, double depth_500nm)
{
	return (0.27583 * depth_380nm + 0.35 * depth_500nm);
}


/* Rayleigh scattering / attenuation */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
double bird::calculate_rayleigh_scattering(double relative_optical_path, double air_pressure)
{
	double rayleigh = 0.0;
	if (relative_optical_path > 0.0)
	{
		double abs_air_mass = absolute_air_mass(relative_optical_path, air_pressure) * 100.0;
		
		double a3 = powf(abs_air_mass, 1.01);
		double a4 = 1 + abs_air_mass;
		double a5 = pow(abs_air_mass, 0.84);
		double a2 = -0.0903 * a5 * (a4 - a3);
				
		rayleigh = exp(a2);
	}
    
	return (rayleigh);
}

/* Ozone Depth / attenuation */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_ozone_depth(double relative_optical_path, double ozone_cm)
{
	double ozone = 0.0;
	if (relative_optical_path > 0.0)
	{
		double o1 = ozone_cm*relative_optical_path;
		double o2 = 0.0003 * (o1) * (o1);
		double pow1 = powf((1.0 + 139.48 * o1), -0.3034);
		
		ozone = 1.0 - 0.1611 * o1 * pow1 - 0.002715 * o1 / (1.0 + 0.044 * o1 + o2);
	}
    
	return (ozone);
}

/* Mixed Gas Depth / attenuation */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_mixed_gas_depth(double relative_optical_path, double air_pressure)
{
	double mixed_gas = 0.0;
	if (relative_optical_path > 0.0)
	{
		double abs_air_mass = absolute_air_mass(relative_optical_path, air_pressure) * 100.0;
		
		mixed_gas = exp(-0.0127*powf(abs_air_mass, 0.26));
	}
    
	return (mixed_gas);
}

/* Water H2O depth / attenuation */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_h20_depth(double relative_optical_path, double h2o_cm)
{
	double h2o = 0.0;
	if (relative_optical_path > 0.0)
	{
		double abs_h2o_mass = relative_optical_path*h2o_cm;
		
		double a1 = 6.385 * abs_h2o_mass;
		double a2 = powf((1 + 79.034 *  abs_h2o_mass), 0.6828);
		h2o = 1 - ((2.4959*abs_h2o_mass) / (a2 + a1));
	}
    
	return (h2o);
}

/* Aerosol depth / attenuation. AKA turbidity */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_aerosol_depth(double relative_optical_path, double tau_aerosol)
{
	double aerosol = 0.0;
	if (relative_optical_path > 0.0)
	{
		double pow1 = -1.0*powf(tau_aerosol, 0.873);
		double pow2 = powf(tau_aerosol, 0.7088);
		double pow3 = powf(relative_optical_path, 0.9108);
		aerosol = exp((pow1)*(1 + tau_aerosol - pow2)*pow3);
	}
    
	return (aerosol);
}

double bird::calculate_aerosol_depth(double relative_optical_path, double depth_380nm, double depth_500nm)
{
	double aerosol = 0.0;
	double tau_aerosol = calculate_aerosol_taua(depth_380nm, depth_500nm);
	
	if (relative_optical_path > 0.0)
	{
		double pow1 = powf(tau_aerosol, 0.873);
		double pow2 = powf(tau_aerosol, 0.7088);
		double pow3 = powf(relative_optical_path, 0.9108);
		aerosol = exp(-1.0*(pow1)*(1 + tau_aerosol - pow2)*pow3);
	}
    
	return (aerosol);
}


/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_taa(double aerosol_tranmittance, double relative_optical_path)
{
	double taa = 0.0;
	if (relative_optical_path > 0.0)
	{
		double pow1 = powf(relative_optical_path, 1.06);
		taa = 1.0 - 0.1 * (1.0 - relative_optical_path + pow1) * (1.0 - aerosol_tranmittance);
	}
    
	return (taa);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_rs(double aerosol_transmittance, double taa, double Ba)
{
	double rs = 0.0;

	rs = 0.0685 + (1.0 - Ba) * (1.0 - aerosol_transmittance / taa);
	
	return (rs);
}

double bird::calculate_rs(double aerosol_transmittance, double taa)
{
	double rs = 0.0;
	double Ba = 0.85;
	
	rs = calculate_rs(aerosol_transmittance, taa, Ba);
	
	return (rs);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_tas(double relative_optical_path, double rayleigh_transmittance, double aerosol_transmittance, double ozone_transmittance, double h2o_transmittance, double mixed_transmittance, double taa, double Ba)
{
	double tas = 0.0;
	if (relative_optical_path > 0.0)
	{	
		double pow1 = powf(relative_optical_path, 1.02);
		tas = 0.79 * ozone_transmittance * mixed_transmittance * h2o_transmittance * taa * (0.5 * (1.0 - rayleigh_transmittance) + Ba * (1.0 - (aerosol_transmittance / taa))) / (1.0 - relative_optical_path + pow1);
	}
    
	return (tas);
}


/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_global_horizontal(double irrad_tilted, double direct_horizontal, double tas, double albedo, double rs)
{
	double global_horizontal = 0.0;
	
	global_horizontal = (direct_horizontal + (irrad_tilted * tas)) / (1.0 - albedo * rs);
				    
	return (global_horizontal);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_diffuse_horizontal(double direct_horizontal, double global_horizontal)
{
	double diffuse_horizontal = 0.0;
		
	diffuse_horizontal = global_horizontal - direct_horizontal;

	return (diffuse_horizontal);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_direct_normal_irradiation(double aerosol_tranmittance, double h20_transmittance, double mixed_gas_transmittance, double ozone_transmittance, double rayleigh_transmittance)
{
	double direct_normal = 0.0;
	
	direct_normal = 0.9662 * aerosol_tranmittance * h20_transmittance * mixed_gas_transmittance * ozone_transmittance * rayleigh_transmittance;
    
	return (direct_normal);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_bird(double relative_optical_path, double air_pressure, double ozone_cm, double h2o_cm, double depth_380nm, double depth_500nm, double Ba)
{
	double direct_normal = 0.0;
	if (relative_optical_path > 0.0)
	{	
		this->tau_aerosol = calculate_aerosol_taua(depth_380nm, depth_500nm);
		this->direct_normal = calculate_bird();
	}
    
	return (this->direct_normal);
}

/*  */
/* Bird 1991,"A Simplified Clear Sky model for Direct and Diffuse Insolation on Horizontal Surfaces" */
/* http://rredc.nrel.gov/solar/pubs/pdfs/tr-642-761.pdf */
/* Returns transmissive ratio. */
double bird::calculate_bird()
{

	if (relative_optical_path > 0.0)
	{	
		this->raleigh = calculate_rayleigh_scattering(this->relative_optical_path, this->air_pressure);/**/
		this->ozone = calculate_ozone_depth(this->relative_optical_path, this->ozone_cm);/**/
		this->mixed = calculate_mixed_gas_depth(this->relative_optical_path, this->air_pressure);/**/
		this->h2o = calculate_h20_depth(this->relative_optical_path, this->h2o_cm);/**/
		this->aerosol = calculate_aerosol_depth(this->relative_optical_path, this->tau_aerosol);/**/
		this->taa = calculate_taa(this->aerosol, this->relative_optical_path);/**/
		this->rs = calculate_rs(this->aerosol, this->taa, this->Ba);/**/
		this->direct_normal = calculate_direct_normal_irradiation(this->aerosol, this->h2o, this->mixed, this->ozone, this->raleigh);/**/
		this->tas = calculate_tas(relative_optical_path, this->raleigh, this->aerosol, this->ozone, this->h2o, this->mixed, this->taa, Ba);
		
		this->direct_beam = ic_prime * this->direct_normal;
		this->plane_angle = calculate_incident_surface_irradiation(this->solar_elevation);
		this->irrad_horizontal  = ic_prime * this->plane_angle;
		this->direct_horizontal = direct_normal * irrad_horizontal;
		this->global_horizontal = calculate_global_horizontal(irrad_horizontal, direct_horizontal, this->tas, this->albedo, this->rs);
		this->diffuse_horizontal = calculate_diffuse_horizontal(this->direct_horizontal, this->global_horizontal);
	}
	else
	{
		this->raleigh = 0.0;/**/
		this->ozone = 0.0;/**/
		this->mixed = 0.0;/**/
		this->h2o = 0.0;/**/
		this->aerosol = 0.0;/**/
		this->taa = 0.0;/**/
		this->rs = 0.0;/**/
		this->direct_normal = 0.0;/**/
		this->tas = 0.0;
		
		this->direct_beam = 0.0;
		this->plane_angle = 0.0;
		this->irrad_horizontal  = 0.0;
		this->direct_horizontal = 0.0;
		this->global_horizontal = 0.0;
		this->diffuse_horizontal = 0.0;
	}
    
	return (this->direct_normal);
}
//
// bird.h
// Version timestamp: 9-30-2018, 7:31 PM
//
// Attribution : Copyright (c) 2018 Northern_Loki (sha256::6F290BF833967127BE26C92C8F6B1C1A3949C55A7EABCEF3ECC785CD2D38D30D)
// License is granted under the Creative Commons Attribution-ShareAlike 4.0 International.  https://creativecommons.org/licenses/by-sa/4.0/
//
# pragma once
# include <math.h>
# include "Constants.h"

class bird
{
public:
	bird();
	bird(double solar_elevation, double ic_prime, double relative_optical_path, double air_pressure, double ozone_cm, double h2o_cm, double depth_380nm, double depth_500nm, double Ba, double tau_aerosol, double albedo);
	~bird();
	
	double absolute_air_mass(double relative_optical_path);
	double absolute_air_mass(double relative_optical_path, double air_pressure);
	double calculate_plane_surface_angle(double solar_zenith);
	double calculate_incident_surface_irradiation(double solar_elevation);
	
	/* Bird */
	double calculate_aerosol_taua(double depth_380nm, double depth_500nm);
	double calculate_aerosol_depth(double relative_optical_path, double depth_380nm, double depth_500nm);
	double calculate_rayleigh_scattering(double relative_optical_path, double air_pressure);
	double calculate_ozone_depth(double relative_optical_path, double ozone_cm);
	double calculate_mixed_gas_depth(double relative_optical_path, double air_pressure);
	double calculate_h20_depth(double relative_optical_path, double air_pressure);
	double calculate_aerosol_depth(double relative_optical_path, double air_pressure);
	double calculate_taa(double aerosol_tranmittance, double relative_air_mass);
	double calculate_rs(double aerosol_tranmittance, double taa, double Ba);
	double calculate_rs(double aerosol_tranmittance, double taa);
	double calculate_tas(double relative_optical_path, double rayleigh_transmittance, double aerosol_transmittance, double ozone_transmittance, double h2o_transmittance, double mixed_transmittance, double taa, double Ba);
	double calculate_global_horizontal(double irrad_tilted, double direct_horizontal, double tas, double albedo, double rs);
	double calculate_diffuse_horizontal(double direct_horizontal, double global_horizontal);
	double calculate_direct_normal_irradiation(double aerosol_tranmittance, double h20_transmittance, double mixed_gas_transmittance, double ozone_transmittance, double rayleigh_transmittance);
	double calculate_bird(double relative_optical_path, double air_pressure, double ozone_cm, double h2o_cm, double depth_380nm, double depth_500nm, double Ba);
	double calculate_bird();
    
	/* Input parameters */
	double solar_elevation;
	double ic_prime;
	double relative_optical_path;
	double air_pressure;
	double ozone_cm;
	double h2o_cm;
	double depth_380nm;
	double depth_500nm;
	double Ba;
	double tau_aerosol;
	double albedo;
	
	/* Generated results */
	double raleigh;
	double ozone;
	double mixed;
	double h2o;
	double aerosol;
	double taa;
	double rs;
	double direct_normal;
	double direct_beam;
	double tas;
	double plane_angle;
	double irrad_horizontal;
	double direct_horizontal;
	double global_horizontal;
	double diffuse_horizontal;
	
};

CC BY-SA 4.0

3 Likes

Very fast life topic > Courtesy “the search” > But i like it >

1 Like