📄 sado_solarradiation.cpp
字号:
d = m_Atmosphere / (m_Atmosphere - m_pDEM->asDouble(x, y))
* (0.0001165 * Decline*Decline - 0.0152 * Decline + A);
d = 1.0 - d;
Dh = Ih / d - Ih;
Slope = m_bHorizon ? 0.0 : m_Slope.asDouble(x, y);
srad = (Dh + Dh * cos(Slope)) / 2.0;
return( srad < 0.0 || srad > m_Solar_Const ? m_Solar_Const : srad );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Shade(double Decline, double Azimuth)
{
//-----------------------------------------------------
m_Shade.Assign(0.0);
if( !m_bBending )
{
int i, x, y;
double dx, dy, dz;
Get_Shade_Params(Decline, Azimuth, dx, dy, dz);
for(i=0; i<Get_NCells() && Set_Progress_NCells(i); i++)
{
if( m_pDEM->Get_Sorted(i, x, y) && !Get_Shade_Complete(x, y) )
{
Set_Shade(x, y, dx, dy, dz);
}
}
}
//-----------------------------------------------------
else
{
int i, x, y, iLock;
for(i=0, iLock=1; i<Get_NCells() && Set_Progress_NCells(i); i++, iLock++)
{
if( m_pDEM->Get_Sorted(i, x, y) && !Get_Shade_Complete(x, y) )
{
if( iLock >= 255 )
iLock = 1;
if( iLock == 1 )
Lock_Create();
Set_Shade_Bended(x, y, iLock);
}
}
}
return( true );
}
//---------------------------------------------------------
void CSADO_SolarRadiation::Set_Shade(int x, int y, double dx, double dy, double dz)
{
for(double ix=x+0.5, iy=y+0.5, iz=m_pDEM->asDouble(x, y); ; )
{
x = (int)(ix += dx);
y = (int)(iy += dy);
iz -= dz;
if( !is_InGrid(x, y) || m_pDEM->asDouble(x, y) > iz )
{
return;
}
m_Shade.Set_Value(x, y, 1);
}
}
//---------------------------------------------------------
void CSADO_SolarRadiation::Set_Shade_Bended(int x, int y, char iLock)
{
double dx, dy, dz;
Get_Shade_Params(m_Decline.asDouble(x, y), m_Azimuth.asDouble(x, y), dx, dy, dz);
for(double ix=x+0.5, iy=y+0.5, iz=m_pDEM->asDouble(x, y); ; )
{
x = (int)(ix += dx);
y = (int)(iy += dy);
iz -= dz;
if( !is_InGrid(x, y) || m_pDEM->asDouble(x, y) > iz || Lock_Get(x, y) == iLock )
{
return;
}
m_Shade.Set_Value(x, y, 1);
//---------------------------------------------
Lock_Set(x, y, iLock);
Get_Shade_Params(m_Decline.asDouble(x, y), m_Azimuth.asDouble(x, y), dx, dy, dz);
}
}
//---------------------------------------------------------
inline bool CSADO_SolarRadiation::Get_Shade_Complete(int x, int y)
{
if( m_Shade.asInt(x, y) == 1 )
{
for(int iy=y-1; iy<=y+1; iy++)
{
for(int ix=x-1; ix<x+1; ix++)
{
if( is_InGrid(ix, iy) && m_Shade.asInt(ix, iy) == 0 )
{
return( false );
}
}
}
return( true );
}
return( false );
}
//---------------------------------------------------------
inline void CSADO_SolarRadiation::Get_Shade_Params(double Decline, double Azimuth, double &dx, double &dy, double &dz)
{
dz = Azimuth + M_PI_180;
dx = sin(dz);
dy = cos(dz);
if( fabs(dx) > fabs(dy) )
{
dy /= fabs(dx);
dx = dx < 0 ? -1 : 1;
}
else if( fabs(dy) > fabs(dx) )
{
dx /= fabs(dy);
dy = dy < 0 ? -1 : 1;
}
else
{
dx = dx < 0 ? -1 : 1;
dy = dy < 0 ? -1 : 1;
}
dz = tan(Decline) * sqrt(dx*dx + dy*dy) * Get_Cellsize();
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
inline bool CSADO_SolarRadiation::Get_Solar_Position(int Day, double Hour, double LAT, double LON, double &Declination, double &Azimuth)
{
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 = Hour - LON * 12.0 / M_PI;
UTime = Hour;
//-----------------------------------------------------
// 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.0 * atan2(Y, (X + R));
//-----------------------------------------------------
// 4. compute sidereal time (degree) at Greenwich local sidereal time at longitude (Degree)
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) ); // previous formula gives same result but is better because of division by zero effects...
return( Declination > 0.0 );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
/*/---------------------------------------------------------
#define CHECK_SHADOW(X, Y) if( (d = SG_Get_Distance(x, y, X, Y) * Get_Cellsize()) > 0.0 )\
{\
if( (t = (m_pDEM->asDouble(X, Y) - z) / d) > tMax )\
{\
tMax = t;\
}\
\
if( (d = zMax / d) < tanDec || d < tMax )\
{\
return( false );\
}\
else if( tMax > tanDec )\
{\
return( true );\
}\
}
//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Solar_Shadow(int x, int y, double Decline, double Azimuth)
{
static double Zenith = 89.9999 * M_DEG_TO_RAD;
//-----------------------------------------------------
if( Decline >= Zenith || Decline > m_Slope.Get_ZMax() )
{
return( false );
}
if( Decline <= 0.0 )
{
return( true );
}
//-----------------------------------------------------
bool bGo;
int dx, nx, ix, iy, iy_Last, jy;
double dy, d, t, tMax, z, zMax, tanDec;
//-----------------------------------------------------
Azimuth = fmod(Azimuth, M_PI_360); if( Azimuth < 0.0 ) Azimuth += M_PI_360;
if( Azimuth < M_PI_180 )
{
dx = 1;
nx = Get_NX() - 1;
}
else
{
dx = -1;
nx = 0;
}
Azimuth = M_PI_090 - Azimuth;
if( Azimuth == M_PI_090 )
Azimuth = 89.9999 * M_DEG_TO_RAD;
else if( Azimuth == -M_PI_090 )
Azimuth = - 89.9999 * M_DEG_TO_RAD;
else if( Azimuth == -M_PI_180 )
Azimuth = -179.9999 * M_DEG_TO_RAD;
else if( Azimuth == -M_PI_270)
Azimuth = -269.9999 * M_DEG_TO_RAD;
dy = tan(Azimuth);
//-----------------------------------------------------
z = m_pDEM->asDouble(x, y);
zMax = m_pDEM->Get_ZMax();
tanDec = tan(Decline);
tMax = 0.0;
//-----------------------------------------------------
for(ix=x, iy=y, bGo=true; ix!=nx && bGo; ix+=dx)
{
iy_Last = iy;
iy = y + (int)floor(dy * (ix - x + (dx > 0 ? 0.5 : -0.5)));
if( iy < 0 )
{
bGo = false;
iy = 0;
}
else if( iy >= Get_NY() )
{
bGo = false;
iy = Get_NY() - 1;
}
//-------------------------------------------------
if( iy_Last < iy )
{
for(jy=iy_Last+1; jy<=iy; jy++)
{
CHECK_SHADOW(ix, jy);
}
}
else if( iy_Last > iy )
{
for(jy=iy_Last; jy>iy; jy--)
{
CHECK_SHADOW(ix, jy);
}
}
else if( iy_Last == iy )
{
CHECK_SHADOW(ix, iy);
}
}
return( false );
}/**/
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -