📄 solarradiation.cpp
字号:
m_pRadiation->Assign(m_pSum);
}
//---------------------------------------------------------
void CSolarRadiation::Execute_SumOfDays(double Latitude_RAD, double Hour_Step, double Hour_Start, double Hour_Stop, int Day_Step, int Day_Start, int Day_Stop)
{
int Day, nDays;
CSG_Grid *pSumOfDays, *pDur_Total;
//-----------------------------------------------------
nDays = Day_Stop - Day_Start;
if( nDays < 1 )
{
Execute_DailySum(Latitude_RAD, Hour_Step, Hour_Start, Hour_Stop, Day_Start);
}
//-----------------------------------------------------
else
{
pSumOfDays = SG_Create_Grid(m_pDTM);
pDur_Total = SG_Create_Grid(m_pDTM);
if( Day_Step < 1 || Day_Step > nDays )
{
Day_Step = 1;
}
//-------------------------------------------------
for(Day=Day_Start; Day<Day_Stop && Set_Progress(Day - Day_Start, nDays); Day+=Day_Step)
{
Get_DailySum(Latitude_RAD, Hour_Step, Hour_Start, Hour_Stop, Day, false);
*m_pSum *= Day_Step;
*pSumOfDays += *m_pSum;
*m_pDuration *= Day_Step;
*pDur_Total += *m_pDuration;
}
if( (Day_Step = Day_Stop - Day) > 0 ) // Rest...
{
Get_DailySum(Latitude_RAD, Hour_Step, Hour_Start, Hour_Stop, Day, false);
*m_pSum *= Day_Step;
*pSumOfDays += *m_pSum;
*m_pDuration *= Day_Step;
*pDur_Total += *m_pDuration;
}
//-------------------------------------------------
m_pRadiation ->Assign(pSumOfDays);
delete(pSumOfDays);
m_pDuration ->Assign(pDur_Total);
delete(pDur_Total);
}
}
///////////////////////////////////////////////////////////
// //
// Daily Sum //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
void CSolarRadiation::Get_DailySum(double Latitude_RAD, double Hour_Step, double Hour_Start, double Hour_Stop, int Day, bool bProgressBar)
{
bool bNight;
double time, Sol_Height, Sol_Azimuth, Solar_Angle, RDIRN, RDIFN, d;
CHillShade HillShade;
//-----------------------------------------------------
Day %= 366;
if( Day < 0 )
{
Day += 366;
}
Hour_Start += Hour_Step / 2.0;
m_pSum ->Assign(0.0);
m_pDuration ->Assign(0.0);
HillShade.Get_Parameters()->Set_Parameter(SG_T("ELEVATION") , PARAMETER_TYPE_Grid , m_pDTM);
HillShade.Get_Parameters()->Set_Parameter(SG_T("SHADE") , PARAMETER_TYPE_Grid , m_pRadiation);
HillShade.Get_Parameters()->Set_Parameter(SG_T("METHOD") , PARAMETER_TYPE_Choice , 3);
HillShade.Get_Parameters()->Set_Parameter(SG_T("EXAGGERATION") , PARAMETER_TYPE_Double , 1.0);
HillShade.Set_Show_Progress(false);
//-----------------------------------------------------
for(time=Hour_Start, bNight=false; time<Hour_Stop && is_Progress(); time+=Hour_Step)
{
Process_Set_Text(CSG_String::Format(_TL("%d. Day of Year, Local Time: %dh %02dm"), Day + 1, (int)time, (int)(60.0 * fmod(time, 1.0))));
if( bProgressBar )
{
Set_Progress(time, 24.0);
}
else
{
Process_Get_Okay(false);
}
//-------------------------------------------------
if( Get_SolarPosition(Day, time, Latitude_RAD, 0.0, Sol_Azimuth, Sol_Height, false) )
{
bNight = false;
HillShade.Get_Parameters()->Set_Parameter(SG_T("AZIMUTH") , PARAMETER_TYPE_Double, (double)M_RAD_TO_DEG * Sol_Azimuth);
HillShade.Get_Parameters()->Set_Parameter(SG_T("DECLINATION") , PARAMETER_TYPE_Double, (double)M_RAD_TO_DEG * Sol_Height);
HillShade.Execute();
for(int n=0; n<Get_NCells(); n++)
{
if( m_pRadiation->is_NoData(n) )
{
m_pSum ->Set_NoData(n);
m_pDuration ->Set_NoData(n);
}
else
{
Solar_Angle = m_pRadiation->asDouble(n);
if( Solar_Angle < M_PI_090 )
{
Get_SolarCorrection(M_PI_090 - Sol_Height, m_pDTM->asDouble(n), RDIRN, RDIFN);
d = m_SolarConstant * Hour_Step * RDIRN * cos(Solar_Angle);
// Dir[q,a] = SConst * SunDuration[q,a] * t^m(q) * cos(AngIn[q,a])
m_pSum ->Add_Value(n, d);
m_pDuration ->Add_Value(n, Hour_Step);
}
}
}
if( Parameters("UPDATE")->asBool() )
{
DataObject_Update(m_pRadiation, 20 * M_DEG_TO_RAD, 90.0 * M_DEG_TO_RAD, true);
}
}
//-------------------------------------------------
else if( !bNight )
{
bNight = true;
if( Parameters("UPDATE")->asBool() )
{
m_pRadiation->Assign(M_PI_090);
DataObject_Update(m_pRadiation, 0, 90.0 * M_DEG_TO_RAD, true);
}
}
}
}
///////////////////////////////////////////////////////////
// //
// Atmosphere //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
void CSolarRadiation::Get_SolarCorrection(double ZenithAngle, double Elevation, double &RDIRN, double &RDIFN)
{
const double AM[32] = // AM Optical air mass in 1 degree increments for zenith angles >=60 [after LIST, 1968; p. 422]
{
2.00, 2.06, 2.12, 2.19, 2.27, 2.36, 2.45, 2.55,
2.65, 2.77, 2.90, 3.05, 3.21, 3.39, 3.59, 3.82,
4.07, 4.37, 4.72, 5.12, 5.60, 6.18, 6.88, 7.77,
8.90, 10.39, 12.44, 15.36, 19.79, 26.96, 26.96, 26.96
};
const double Pressure_0 = 1013.0; // Standard atmospheric pressure = 1013 mb
int i;
double z, AMASS, AMASS2, AW, TW, TD, TDC;
//-----------------------------------------------------
z = M_RAD_TO_DEG * ZenithAngle;
if( z <= 60.0 )
{
AMASS = 1.0 / cos(ZenithAngle);
}
else
{
z -= 60.0;
i = (int)z;
AMASS = AM[i] + (z - i) * (AM[i + 1] - AM[i]);
}
z = m_Pressure / pow(10.0, Elevation * 5.4667E-05); // P Barometric pressure in mb
AMASS2 = AMASS * z / Pressure_0;
//-----------------------------------------------------
switch( m_Method )
{
case 0: default:
RDIRN = pow(m_Transmittance, AMASS2);
RDIFN = 0.271 - 0.294 * RDIRN;
break;
case 1:
AW = 1.0 - 0.077 * pow(m_Water * AMASS2, 0.3); // AW Accounts for absorption by water vapour
TW = pow(0.975, m_Water * AMASS); // TW Accounts for scattering by water vapour
TD = pow(0.950, m_Water * m_Dust / 100.0); // TD Accounts for scattering by dust
TDC = pow(0.900, AMASS2) + 0.026 * (AMASS2 - 1.0); // TDC Accounts for scattering by a dust free atmosphere
RDIRN = AW * TW * TD * TDC;
RDIFN = 0.5 * (AW - RDIRN);
break;
}
// Optical Pathlength (solar zenith angle < 80Degree) =>
// m(t) = exp(-0.000118 * Elevation - 1.638 * 10e-9 * Elevation) / cos(SolarZenithAngle)
// if( Method == 2 && Zenithangle < 80 ) // nach bla...
// RDIRN = pow(m_Transmittance, exp((-0.000118 - 1.638e-8 * Elevation) * Elevation) / cos(ZenithAngle));
}
/* The original TAPES-G source code for the optical air mass computation
C ==================================================================
SUBROUTINE SOLAR(ZA,RDIRN,RDIFN,ITEST)
COMMON/SOL1/U,D,P,TRANSM
PARAMETER (PI=3.14159265358979323846)
PARAMETER (DTOR=PI/180.)
DIMENSION AM(32)
DATA AM/2.0,2.06,2.12,2.19,2.27,2.36,2.45,2.55,2.65,2.77,2.9,
* 3.05,3.21,3.39,3.59,3.82,4.07,4.37,4.72,5.12,5.6,6.18,6.88,
* 7.77,8.9,10.39,12.44,15.36,19.79,26.96,26.96,26.96/
DATA PO/1013./
C ***************************************************************
C AM Optical air mass in 1 degree increments for zenith angles
C >=60 [LIST, 1968; p. 422]
C U Water content of a vertical slice of atmosphere in cm:
C 1.5 to 1.7, average=1.68
C D Dust factor: 1=100 ppm (standard); 2=300 ppm
C P Barometric pressure in mb
C PO Standard atmospheric pressure = 1013 mb
C TRANSM Transmittance of the atmosphere (0.6-0.8)
C AW Accounts for absorption by water vapour
C TW Accounts for scattering by water vapour
C TD Accounts for scattering by dust
C TDC Accounts for scattering by a dust free atmosphere
C **************************************************************
C Compute optical air mass
C
IF(ZA.LE.60.) THEN
AMASS=1./COS(ZA*DTOR)
ELSE
Y=ZA-59.
I=INT(Y)
AMASS=AM(I)+(Y-FLOAT(I))*(AM(I+1)-AM(I))
ENDIF
AMASS2=AMASS*P/PO
C --------------------------------------------------------------
C Account for atmospheric effects using either a lumped atmos-
C pheric transmittance approach (ITEST=1) or by calculating the
C components (ITEST=2)
C
IF(ITEST.EQ.1) THEN
RDIRN=TRANSM**AMASS2
RDIFN=0.271-0.294*RDIRN
ELSE
AW=1.0-0.077*(U*AMASS2)**0.3
TW=0.975**(U*AMASS)
TD=0.95**(U*D)
TDC=0.9**AMASS2+0.026*(AMASS2-1.0)
RDIRN=AW*TW*TD*TDC
RDIFN=0.5*(AW-RDIRN)
ENDIF
RETURN
END
C ====================================================================*/
///////////////////////////////////////////////////////////
// //
// Solar Position //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSolarRadiation::Get_SolarPosition(int Day, double Time, double LAT, double LON, double &Azimuth, double &Declination, bool bDegree)
{
static const int Day2Month[13] = { 0, 31, 49, 80, 109, 140, 170, 201, 232, 272, 303, 333, 366 };
static const double ECLIPTIC_OBL = M_DEG_TO_RAD * 23.43999; // obliquity of ecliptic
int i;
double JD, T, M, L, X, Y, Z, R, UTime,
DEC, RA, theta, tau,
Month, Year = 2000;
//-----------------------------------------------------
for(Month=1, i=0; i<=12; i++)
{
if( Day < Day2Month[i] )
{
Month = i;
Day -= Day2Month[i - 1];
break;
}
}
if( Month < 1 || Month > 12 ) // nur Tag (1 - 365) zaehlt...
{
Month = 1;
}
if( Month <= 2 )
{
Month += 12;
Year -= 1;
}
//-----------------------------------------------------
// UTime = Time - LON * 12.0 / M_PI;
UTime = Time;
//-----------------------------------------------------
// 1. Julian Date...
JD = (int)(365.25 * Year) + (int)(30.6001 * (Month + 1)) - 15 + 1720996.5 + Day + UTime / 24.0;
T = (JD - 2451545.0 ) / 36525.0; // number of Julian centuries since 2000/01/01 at 12 UT (JD = 2451545.0)
//-----------------------------------------------------
// 2. Solar Coordinates (according to: Jean Meeus: Astronomical Algorithms), accuracy of 0.01 degree
M = M_DEG_TO_RAD * (357.52910 + 35999.05030 * T - 0.0001559 * T * T - 0.00000048 * T * T * T); // mean anomaly
L = M_DEG_TO_RAD * ( (280.46645 + 36000.76983 * T + 0.0003032 * T * T) // mean longitude
+ ( (1.914600 - 0.004817 * T - 0.000014 * T * T) * sin(M)
+ (0.019993 - 0.000101 * T) * sin(2 * M) + 0.000290 * sin(3 * M) // true longitude
)
);
//-----------------------------------------------------
// 3. convert ecliptic longitude to right ascension RA and declination delta
X = cos(L);
Y = cos(ECLIPTIC_OBL) * sin(L);
Z = sin(ECLIPTIC_OBL) * sin(L);
R = sqrt(1.0 - Z*Z);
DEC = atan2(Z, R);
RA = 2 * atan2(Y, (X + R));
//-----------------------------------------------------
// 4. compute sidereal time (degree) at Greenwich local sidereal time at longitude (DegreeE)
theta = LON + M_DEG_TO_RAD * (280.46061837 + 360.98564736629 * (JD - 2451545.0) + T * T * (0.000387933 - T / 38710000.0));
//-----------------------------------------------------
// 5. compute local hour angle (degree)
tau = theta - RA;
//-----------------------------------------------------
// 6. convert (tau, delta) to horizon coordinates (h, az) of the observer
Declination = asin( sin(LAT) * sin(DEC) + cos(LAT) * cos(DEC) * cos(tau) );
Azimuth = atan2( -sin(tau) * cos(DEC), cos(LAT) * sin(DEC) - sin(LAT) * cos(DEC) * cos(tau) );
//Azimuth = atan2( -sin(Tau), cos(LAT) * tan(DEC) - sin(LAT) * cos(Tau) ); // erstere ist besser wegen division by zero effects...
if( bDegree )
{
Declination *= M_RAD_TO_DEG;
Azimuth *= M_RAD_TO_DEG;
}
return( Declination >= 0.0 );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -