Rev 1221 |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
#include "solarradiation.h"
#include "hemigrid.h"
//#include "globals.h"
//#include "modell.h"
//#include "DBUtil.h"
//--------------------- SolarRad -------------
// ------ Mini-Modul zur Berechnung der Globalstrahlung
// eigentlich ähnlich wie Strahlungsmodell Picus 1.x,
// jedoch gleich selber machen ist einfacher....
void SolarRadiation::Setup(double BreiteGrad)
{
Latitude=BreiteGrad * M_PI / 180.;
for (int day=0; day<366; day++) {
calcDay(day);
}
}
//////////////////////////////////////////////////////
/// Calculate factor for optical mass (used for radiation calculation)
/// http://www.sortie-nd.org/help/manuals/developer_documentation/cplusplus/light.html
//////////////////////////////////////////////////////
double GetOpticalMass(const double elevation_rad)
{
if (elevation_rad < 5*M_PI/180.)
return 10.39;
if (elevation_rad < 15*M_PI/180.)
return 5.6;
if (elevation_rad < 25*M_PI/180.)
return 2.9;
// if >= 25°:
// cos(zenith) == cos(90° - elevation)....
// here i am not completely sure: but sec(x)=1/cos(x)
// seems to be right. sec(25°)=2.3, sec(20)=2.9 (see above!), sec(0)=1
// the concrete steps above (5,15,25°) can be observed in cumulated outputs
// e.g. global radiation over elevation - maybe a smoother variant should be used.
return 1./cos(M_PI/2. - elevation_rad);
}
//////////////////////////////////////////////////////
// calc relative radiation matrix.
//////////////////////////////////////////////////////
void SolarRadiation::calculateRadMatrix(const float Step_deg, HemiGrid &Grid)
{
// const double transmissionskoeff[12]={
// 0.4, 0.48, 0.47, 0.46, 0.47, 0.44,
// 0.48, 0.44, 0.39, 0.38, 0.32, 0.34};
// elevation angle threshold
const double min_elevation_angle = 0*M_PI/180.;
double latidude = Latitude;
Grid.setup(Step_deg);
int cntAzimuth = Grid.matrixCountAzimuth();
int cntElevation = Grid.matrixCountElevation();
int day;
double sin_declination, declination, eccentricity;
double current_transmission;
double sun_day, diff_hour;
double hour, sun_hour;
double step = 1./60.; // 1 minute
bool sun_rise, sun_down;
double total_rad_sum = 0;
double elevation, sin_elevation, cos_azimuth, azimuth;
double Iincoming;
for (day = 0; day < 366; day++) {
// skip days due to vegetation period limit:
if (day<mMinDay || day>mMaxDay)
continue;
sin_declination=0.39785*sin((278.97+0.9856*day+1.9165*sin((356.6+0.9856*day)*M_PI/180.))*M_PI/180.);
declination=asin(sin_declination); // declination of earth (rad)
// solarconstant: Gassel p. 15
SolarConstant=1356.5+48.5*cos(0.01721*(day-15));
sun_day = (day-1)*2.*M_PI/365.; // rad
// difference between "real sunhour" and hour
diff_hour = 12+0.12357*sin(sun_day)-0.004289*cos(sun_day)+0.153809*sin(2.*sun_day)+0.060783*cos(2.*sun_day);
// eccentricity: SORTIE (http://www.sortie-nd.org/help/manuals/developer_documentation/cplusplus/light.html)
eccentricity = 1.000110 + 0.034221* cos(sun_day) + 0.001280*sin(sun_day) + 0.000719 * cos(2*sun_day) + 0.000077* sin(2*sun_day);
hour = 0;
sun_rise = sun_down = false;
//current_transmission = transmissionskoeff[(day/30)%12];
current_transmission = 0.75;
while (!sun_down) {
sun_hour=(hour-diff_hour)*15./180*M_PI; // s in rad
// sin elevation: =sin(lat)*sin(decl) + cos(lat)*cos(sunhour)*cos(decl)
sin_elevation=sin(latidude)*sin_declination + cos(latidude)*cos(sun_hour)*cos(declination);
elevation =asin(sin_elevation); // Sonnenhöhe [rad]
// check for state of the sun
if (!sun_rise && elevation>min_elevation_angle)
sun_rise=true;
if (sun_rise && elevation<min_elevation_angle)
sun_down=true;
if (elevation>=min_elevation_angle) {
// calculate azimuth
cos_azimuth = ( sin_elevation * sin(latidude) - sin_declination) / ( cos(elevation)*cos(latidude) );
azimuth = acos(cos_azimuth);
// azimuth is always > 0: use sign of sun_hour: if negative, than direction is "east" (before noon)
if (sun_hour < 0.)
azimuth*=-1;
// calculate the energy (reduction due to transmissioncoeffiecients..)
// I = eccentricity * trans^optical_mass * cos(zeninthangle!)
Iincoming = eccentricity * pow(current_transmission, GetOpticalMass(elevation)) * cos(M_PI/2. - elevation);
total_rad_sum += Iincoming;
// azimuth is 0 for "south". the grid is north-oriented -> reverse (azimuth(North)=0)
if (azimuth <= 0.)
azimuth += M_PI;
else
azimuth -= M_PI;
// store in grid:
Grid.rGet(azimuth, elevation) += Iincoming;
} // elevation>min_elevation_angle
hour += step;
if (hour>24)
throw QString("tsolarrad::calcmatrix: elevation never above minimum elevation.");
} // while
} // for (each day)
// normalize the array and add the diffus radiation...
// direct: sum(all rad) should be fraction of direct radiation
double direct_rad_multiplier = (1.-mDiffusRadFraction) / ( total_rad_sum );
// diffuse: the light is evenly distributed (above the minimum angle)
// add a little to lower bound to avoid rounding down
int min_elev_index = int( (min_elevation_angle+0.0001) / (M_PI / 2.) * cntElevation);
double diffuse_rad_addon = mDiffusRadFraction / ((cntElevation - min_elev_index) * cntAzimuth);
for (int i=min_elev_index; i<cntElevation; i++) {
for (int j=0;j<cntAzimuth;j++) {
double &cell_value = Grid.rGetByIndex(j, i);
cell_value*=direct_rad_multiplier;
cell_value+=diffuse_rad_addon;
}
}
}
void SolarRadiation::calcDay(int day)
{
// Berechnungen pro Tag:
// Deklination:
//sin D=0,39785*SIN((278,97+0,9856*G8+1,9165*SIN((356,6+0,9856*G8)*M_PI()/180))*M_PI()/180)
double sinDekl=0.39785*sin((278.97+0.9856*day+1.9165*sin((356.6+0.9856*day)*M_PI/180.))*M_PI/180.);
double Dekl=asin(sinDekl); // Deklination in Grad
//S = 1356,5+48,5*cos[0,01721*(n-15)] Gassel p. 15
SolarConstant=1356.5+48.5*cos(0.01721*(day-15));
// Sonnentag
double SunDay;
SunDay=(day-1)*2.*M_PI/365.; // [rad]
// Unterschied wahre Sonnenstunde - Stunde
double Diff=12+0.12357*sin(SunDay)-0.004289*cos(SunDay)+0.153809*sin(2.*SunDay)+0.060783*cos(2.*SunDay);
// Integration über den Tag
double Sum=0., Length=0., Rad;
// Sonnenstunde
double Hour, SunHour, sinH, Altitude;
bool SunRise=false;
bool SunDown=false;
Hour=0.;
const double Step=1/60.; // Minutentakt
// berechnung der tageslänge methode 2:
// tag: wenn strahlung > threshold
double RadThreshold = -1; // not used: Picus specific: =GModell->Settings["BBDayLengthMode"].AsDouble();
const double transmissionskoeff[12]={
0.4, 0.48, 0.47, 0.46, 0.47, 0.44,
0.48, 0.44, 0.39, 0.38, 0.32, 0.34};
while (!SunDown)
{
SunHour=(Hour-Diff)*15./180*M_PI; // s in rad
// sin Sonnenhöhe: =sin(Breite)*sin(Dekl) + cos(Breite)*cos(Sonnenstunde)*cos(Deklination)
sinH=sin(Latitude)*sinDekl + cos(Latitude)*cos(SunHour)*cos(Dekl);
Altitude=asin(sinH); // Sonnenhöhe [rad]
// calculate azimuth-angle
if (!SunRise && Altitude>0)
SunRise=true;
if (SunRise && Altitude<0)
SunDown=true;
if (Altitude>0) {
// Length: Tageslänge [Stunden]
if (RadThreshold==-1)
Length+=Step; // default
else if (SolarConstant*sinH*transmissionskoeff[(day/30)%12] > RadThreshold) {
// SolarConst[W/m2]*frac*frac=[W/m2]
Length+=Step; // strahlungs-grenze-fall
}
// Rad: W * h * frac: Wh
Rad=SolarConstant * Step * sinH;
Sum+=Rad;
}
Hour+=Step;
}
RadExtraTerrestrisch[day]=Sum*3.6; // Wh-> kJ (*3600/1000)
DayLength[day]=Length*3600; // sekunden
}