sado_solarradiation.cpp

这是一个GPS相关的程序

	if( m_pSumDirect )
	{
		m_pSumDirect->Assign(0.0);
		m_pSumDirect->Set_Unit(_TL("Joule"));
		DataObject_Set_Colors(m_pSumDirect, c);

		if( m_bUpdateDirect )
		{
			m_TmpDirect.Create(*Get_System(), GRID_TYPE_Float);
			DataObject_Update(m_pSumDirect, true);
		}
	}

	if( m_pSumDiffus )
	{
		m_pSumDiffus->Assign(0.0);
		m_pSumDiffus->Set_Unit(_TL("Joule"));
		DataObject_Set_Colors(m_pSumDiffus, c);

		if( m_bUpdateDiffus )
		{
			m_TmpDiffus.Create(*Get_System(), GRID_TYPE_Float);
			DataObject_Update(m_pSumDiffus, true);
		}
	}

	if( m_pSumTotal )
	{
		m_pSumTotal ->Assign(0.0);
		m_pSumTotal ->Set_Unit(_TL("Joule"));
		DataObject_Set_Colors(m_pSumTotal , c);

		if( m_bUpdateTotal )
		{
			m_TmpTotal.Create(*Get_System(), GRID_TYPE_Float);
			DataObject_Update(m_pSumTotal , true);
		}
	}

	//-----------------------------------------------------
	Process_Set_Text(_TL("initialising gradient..."));

	m_Shade .Create(*Get_System(), GRID_TYPE_Byte);
	m_Slope .Create(*Get_System(), GRID_TYPE_Float);
	m_Aspect.Create(*Get_System(), GRID_TYPE_Float);

	for(y=0; y<Get_NY() && Set_Progress(y); y++)
	{
		for(x=0; x<Get_NX(); x++)
		{
			double	s, a;

			if( m_pDEM->Get_Gradient(x, y, s, a) )
			{
				m_Slope .Set_Value(x, y, s);
				m_Aspect.Set_Value(x, y, a);
			}
			else
			{
				m_Slope .Set_NoData(x, y);
				m_Aspect.Set_NoData(x, y);
			}
		}
	}

	//-----------------------------------------------------
	if( m_bBending )
	{
		Process_Set_Text(_TL("initialising planetary bending..."));

		CSG_Grid	*pLat	= Parameters("GRD_LAT")->asGrid(),
				*pLon	= Parameters("GRD_LON")->asGrid();

		m_Lat		.Create(*Get_System(), GRID_TYPE_Float);
		m_Lon		.Create(*Get_System(), GRID_TYPE_Float);
		m_Decline	.Create(*Get_System(), GRID_TYPE_Float);
		m_Azimuth	.Create(*Get_System(), GRID_TYPE_Float);

		if( pLat || pLon )
		{
			if( pLat )
			{
				m_Lat	= *pLat;
				m_Lat	*= M_DEG_TO_RAD;
			}

			if( pLon )
			{
				m_Lon	= *pLon;
				m_Lon	*= M_DEG_TO_RAD;
			}
		}

		//-------------------------------------------------
		else
		{
			double	d, dx, dy, dxA, dyA;

			d	= M_DEG_TO_RAD / (Parameters("RADIUS")->asDouble() * M_PI / 180.0);

			switch( Parameters("LON_OFFSET")->asInt() )
			{
			case 0:	dxA	= Get_System()->Get_Extent().Get_XMin();	break;	// left
			case 1:	dxA	= Get_System()->Get_Extent().Get_XCenter();	break;	// center
			case 2:	dxA	= Get_System()->Get_Extent().Get_XMax();	break;	// right
			case 3:	dxA	= Parameters("LON_REF_USER")->asDouble();	break;	// user defined coordinate
			}

			switch( Parameters("LAT_OFFSET")->asInt() )
			{
			case 0:	dyA	= Get_System()->Get_Extent().Get_YMin();	break;	// bottom
			case 1:	dyA	= Get_System()->Get_Extent().Get_YCenter();	break;	// center
			case 2:	dyA	= Get_System()->Get_Extent().Get_YMax();	break;	// top
			case 3:	dyA	= Parameters("LAT_REF_USER")->asDouble();	break;	// user defined coordinate
			}

			dxA	 = d * (Get_XMin() - dxA);
			dyA	 = d * (Get_YMin() - dyA) + m_Latitude;
			d	*= Get_Cellsize();

			for(y=0, dy=dyA; y<Get_NY() && Set_Progress(y); y++, dy+=d)
			{
				for(x=0, dx=dxA; x<Get_NX(); x++, dx+=d)
				{
					m_Lat.Set_Value(x, y, dy);
					m_Lon.Set_Value(x, y, dx);
				}
			}
		}
	}

	//-----------------------------------------------------
	return( true );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Finalise(double SumFactor)
{
	//-----------------------------------------------------
	if( m_pSumDirect )
	{
		if( m_bUpdateDirect )	m_pSumDirect->Assign(&m_TmpDirect);
		if( SumFactor != 1.0 )	*m_pSumDirect	*= SumFactor;
	}

	if( m_pSumDiffus )
	{
		if( m_bUpdateDiffus )	m_pSumDiffus->Assign(&m_TmpDiffus);
		if( SumFactor != 1.0 )	*m_pSumDiffus	*= SumFactor;
	}

	if( m_pSumTotal )
	{
		if( m_bUpdateTotal )	m_pSumTotal->Assign(&m_TmpTotal);
		if( SumFactor != 1.0 )	*m_pSumTotal	*= SumFactor;
	}

	//-----------------------------------------------------
	m_Shade		.Destroy();
	m_Slope		.Destroy();
	m_Aspect	.Destroy();

	m_Lat		.Destroy();
	m_Lon		.Destroy();
	m_Decline	.Destroy();
	m_Azimuth	.Destroy();

	m_TmpDirect	.Destroy();
	m_TmpDiffus	.Destroy();
	m_TmpTotal	.Destroy();

	//-----------------------------------------------------
	return( true );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Insolation(void)
{
	//-----------------------------------------------------
	if( Initialise() )
	{
		if( m_bMoment )
		{
			Get_Insolation(m_Day_A, m_Hour);

			Finalise();
		}

		//-------------------------------------------------
		else
		{
			for(int Day=m_Day_A; Day<=m_Day_B && Process_Get_Okay(false); Day+=m_dDays)
			{
				for(double Hour=m_Hour; Hour<24.0 && Process_Get_Okay(false); Hour+=m_dHour)
				{
					Process_Set_Text(CSG_String::Format(SG_T("%s: %d(%d-%d), %s: %f"), _TL("Day"), Day, m_Day_A, m_Day_B, _TL("Hour"), Hour));

					if( m_bUpdateDirect )	m_pSumDirect->Assign(0.0);
					if( m_bUpdateDiffus )	m_pSumDiffus->Assign(0.0);
					if( m_bUpdateTotal )	m_pSumTotal ->Assign(0.0);

					if( Get_Insolation(Day, Hour) )
					{
						if( m_bUpdateDirect )
						{
							m_TmpDirect	+= *m_pSumDirect;
							DataObject_Update(m_pSumDirect);
						}

						if( m_bUpdateDiffus )
						{
							m_TmpDiffus	+= *m_pSumDiffus;
							DataObject_Update(m_pSumDiffus);
						}

						if( m_bUpdateTotal )
						{
							m_TmpTotal	+= *m_pSumTotal;
							DataObject_Update(m_pSumTotal);
						}
					}
				}
			}

			Finalise(m_dHour / (24.0 * (1 + m_Day_B - m_Day_A)));	// *m_pSumDirect	*= m_dHour / D->size();
		}
	}

	//-----------------------------------------------------
	return( true );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Insolation(int Day, double Hour)
{
	double	Azimuth, Decline;

	if( m_bBending )
	{
		bool	bGo	= false;

		for(int y=0; y<Get_NY() && Set_Progress(y); y++)
		{
			for(int x=0; x<Get_NX(); x++)
			{
				Get_Solar_Position(Day, Hour, m_Lat.asDouble(x, y), m_Lon.asDouble(x, y), Decline, Azimuth);

				m_Azimuth.Set_Value(x, y, Azimuth);
				m_Decline.Set_Value(x, y, Decline);

				if( Decline > 0.0 )
				{
					bGo	= true;
				}
			}
		}

		if( bGo )
		{
			return( Set_Insolation(Decline, Azimuth) );
		}
	}
	else if( Get_Solar_Position(Day, Hour, m_Latitude, 0.0, Decline, Azimuth) )
	{
		return( Set_Insolation(Decline, Azimuth) );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Vapour_Exponent(double VapourPressure)
{
	return( 0.916  - (VapourPressure > 0.0 ? 0.05125 * sqrt(VapourPressure) : 0.0) );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Vapour_A(double VapourPressure)
{
	return( 0.4158 + (VapourPressure > 0.0 ? 0.0399  * sqrt(VapourPressure) : 0.0) );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Set_Insolation(double Decline, double Azimuth)
{
	int		x, y;
	double	a, e, Direct, Diffus;

	Get_Shade(Decline, Azimuth);

	e	= Get_Vapour_Exponent	(m_VP);
	a	= Get_Vapour_A			(m_VP);

	for(y=0; y<Get_NY() && Set_Progress(y); y++)
	{
		for(x=0; x<Get_NX(); x++)
		{
			if( m_pDEM->is_NoData(x, y) )
			{
				if( m_pSumDirect )	m_pSumDirect->Set_NoData(x, y);
				if( m_pSumDiffus )	m_pSumDiffus->Set_NoData(x, y);
				if( m_pSumTotal )	m_pSumTotal ->Set_NoData(x, y);
			}
			else
			{
				if( m_pVP != NULL )
				{
					e	= Get_Vapour_Exponent	(m_pVP->asDouble(x, y));
					a	= Get_Vapour_A			(m_pVP->asDouble(x, y));
				}

				if( m_bBending )
				{
					Azimuth	= m_Azimuth.asDouble(x, y);
					Decline	= m_Decline.asDouble(x, y);
				}

				if( Decline > 0.0 )
				{
					Direct	= Get_Solar_Direct(x, y, Decline, Azimuth, e);
					Diffus	= Get_Solar_Diffus(x, y, Decline, a      , e);

					if( m_pSumDirect )	m_pSumDirect->Add_Value(x, y, Direct);
					if( m_pSumDiffus )	m_pSumDiffus->Add_Value(x, y, Diffus);
					if( m_pSumTotal  )	m_pSumTotal ->Add_Value(x, y, Direct + Diffus);
				}
			}
		}
	}

	return( true );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Reduction(double Elevation, double Decline, double Reduction)
{
	return( Decline > 0.0 ? (m_Solar_Const * pow(Reduction, (1.0 - Elevation / m_Atmosphere) / sin(Decline))) : 0.0 );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Direct(int x, int y, double Decline, double Azimuth, double Exponent)
{
	if( m_Shade.asInt(x, y) == 0 )
	{
		double	Angle, Slope;

		Slope	= m_bHorizon ? 0.0 : m_Slope.asDouble(x, y);

		Angle	= cos(Slope) * cos(Decline - M_PI_090)
				+ sin(Slope) * sin(M_PI_090 - Decline)
				* cos(Azimuth - m_Aspect.asDouble(x, y));

		if( Angle > 0.0 )
		{
			return( Angle * Get_Solar_Reduction(m_pDEM->asDouble(x, y), Decline, Exponent) );
		}
	}

	return( 0.0 );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Diffus(int x, int y, double Decline, double A, double Exponent)
{
	double	Dh, Ih, d, Slope, srad;

	Ih		= Get_Solar_Reduction(m_pDEM->asDouble(x, y), Decline, Exponent) * sin(Decline);
//	Decline	*= M_RAD_TO_DEG;	// Fehler in SADO ????!!!!!