solarradiation.cpp

计算太阳辐射的程序：(1)计算太阳在水平面上的瞬时辐照度和平均辐照度；(2)计算太阳在倾斜面上的瞬时辐照度和平均辐照度；(3)计算太阳的赤纬角、高度角、时角等。

//
//SolarRadiation.cpp
//
#include"StdAfx.h"
#include"SolarRadiation.h"

double SolarRadiation::RiJiao(double JD, double JF, double NF, double Y, double R, double S, double F)
{
	double A = NF/4;
	double N0 = 79.6764+0.2422*(NF-1985)-floor((NF-1985)/4);
	double B = A-floor(A);
	double C = 32.8;
	if(Y<=2)
		C = 30.6;
	if(B==0&&Y>2)
		C = 31.8;
	double G = floor(30.6*Y-C+0.5)+R;
	double L = (JD+JF/60)/15;                //经度修正
	double H = S-8+F/60;                     //时刻修正
	double N = G+(H-L)/24;                   //计算积日
	double rijiao = 2*PI*(N-N0)/365.2422;    //计算日角
	return rijiao;                           //返回日角（弧度值）
};

double SolarRadiation::ChiWeiJiao(double JD, double JF, double NF, double Y, double R, double S, double F)
{
	double rijiao = RiJiao(JD,JF,NF,Y,R,S,F);
	double chiweijiao = 0.3723+23.2567*sin(rijiao)+0.1149*sin(2*rijiao)-0.1712*sin(3*rijiao)-0.758*cos(rijiao)
		                +0.3656*cos(2*rijiao)+0.0201*cos(3*rijiao);
	return chiweijiao;               //返回赤纬角（角度值）
};

double SolarRadiation::ShiCha(double JD, double JF, double NF, double Y, double R, double S, double F)
{
	double rijiao = RiJiao(JD,JF,NF,Y,R,S,F);
	double shicha = 0.0028-1.9857*sin(rijiao)+9.9059*sin(2*rijiao)-7.0924*cos(rijiao)-0.6882*cos(2*rijiao);
	return shicha;                             //返回时差（以分钟为单位）
};

double SolarRadiation::ShiJiao(double JD, double JF, double NF, double Y, double R, double S, double F)
{
	double sd = S+(F-(120-(JD+JF/60))*4)/60;     //根据北京时计算出地方时
	double sr = sd+ShiCha(JD,JF,NF,Y,R,S,F)/60;  //时差修正
	double shijiao = (sr-12)*15;
	return shijiao;                              //返回时角（角度值）

};

double SolarRadiation::GaoDuJiaoZhengXian(double JD, double JF, double NF, double Y, double R, double S, double F, double WD, double WF)
{
	double chiweijiao1 = ChiWeiJiao(JD,JF,NF,Y,R,S,F);
	double chiweijiao2 = chiweijiao1*PI/180;//将赤纬角由角度值变为弧度值
	double shijiao1 = ShiJiao(JD,JF,NF,Y,R,S,F);
	double shijiao2 = shijiao1*PI/180;//将时角由角度值变为弧度值
	double weidu = (WD+WF/60)*PI/180;//将纬度由角度值变为弧度值
	double gaodujiaozhengxian = sin(chiweijiao2)*sin(weidu)+cos(chiweijiao2)*cos(weidu)*cos(shijiao2);
	return gaodujiaozhengxian;
};

double SolarRadiation::XiShouTouSheBi(double JD, double JF, double NF, double Y, double R, double S, double F, double WD, double WF,
								   double Ta, double Rh)
{
	double ludianwendu = 2371.78/(11.286-log10(Rh*611*pow(10,8.5*(Ta-0.01)/(Ta+273.15))))-273.16;
	double gaodujiaozhengxian = GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF);

	double As;                    //As为经验常数，4-6月份为-0.0229，其他月份为0.0203
	if(Y>=4&&Y<=6)
		As = -0.0229;
	else
		As = 0.0203;

	double xishoutoushebi = 1-0.271*pow(exp(0.07074*ludianwendu+As)/gaodujiaozhengxian,0.303);
	return xishoutoushebi;
};

double SolarRadiation::SanSheTouSheBi(double JD, double JF, double NF, double Y, double R, double S, double F, double WD, double WF,
									  double Pa, double Fg)
{
	double gaodujiaozhengxian = GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF);
	double F0 = 0.085-0.247*log10(Pa*gaodujiaozhengxian/1000);

	double sanshetoushebi = (1-F0)/(1-F0*Fg);
	return sanshetoushebi;
};

double SolarRadiation::ShuiPingShunShiFuZhaoDu(double JD, double JF, double NF, double Y, double R, double S, double F, double WD,
											   double WF, double Ta, double Rh, double Pa, double Fg, int YunLeiXing, double YunLiang)
{
	double gaodujiaozhengxian = GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF);
    double shuipingshunshifuzhaodu;

	if(gaodujiaozhengxian>0)
	{
		double taiyangchangshu = 1353;
		double xishoutoushebi = XiShouTouSheBi(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh);
		double sanshetoushebi = SanSheTouSheBi(JD,JF,NF,Y,R,S,F,WD,WF,Pa,Fg);
		double wuyunfuzhaodu = xishoutoushebi*0.349*taiyangchangshu*gaodujiaozhengxian+sanshetoushebi*0.651*taiyangchangshu*gaodujiaozhengxian;//无云时的水平面辐照度
		
		double Bc,Cc;//根据云的类型确定修正因子
		switch(YunLeiXing){
		case 0:
			Bc = 82.2;
			Cc = 0.079;
			break;
		case 1:
			Bc = 87.1;
			Cc = 0.148;
			break;
		case 2:
			Bc = 52.5;
			Cc = 0.112;
			break;
		case 3:
			Bc = 39.0;
			Cc = 0.063;
			break;
		case 4:
			Bc = 34.7;
			Cc = 0.104;
			break;
		case 5:
			Bc = 23.8;
			Cc = 0.159;
			break;
		case 6:
			Bc = 15.4;
			Cc = 0.028;
			break;
		default:
			break;
		}
		
		//将无云时的水平面辐照度乘以云的修正系数。若无云（云量为0），则云的修正系数为0。
		shuipingshunshifuzhaodu = wuyunfuzhaodu*(1-(1-(Bc/94.4)*exp((0.059-Cc)/gaodujiaozhengxian))*YunLiang*YunLiang);

		//即使太阳高度角正弦大于0，在太阳刚刚升起时，水平瞬时辐照度仍有可能小于0。
		if(shuipingshunshifuzhaodu<0)
			shuipingshunshifuzhaodu = 0;
	}
	else if(gaodujiaozhengxian<=0)
	{
		shuipingshunshifuzhaodu = 0;
	}
	
	return shuipingshunshifuzhaodu;
};

double SolarRadiation::TaiYangFangWeiJiao(double JD, double JF, double NF, double Y, double R, double S, double F, double WD, double WF)
{
	double gaodujiaozhengxian = GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF);
	double chiweijiao1 = ChiWeiJiao(JD,JF,NF,Y,R,S,F);
	double chiweijiao2 = chiweijiao1*PI/180;//将赤纬角由角度值变为弧度值
	double weidu = (WD+WF/60)*PI/180;//将纬度由角度值变为弧度值

	double fangweijiaoyuxian,fangweijiao;
	if(gaodujiaozhengxian==1)
		fangweijiao = 0;
	else
	{
		fangweijiaoyuxian = (gaodujiaozhengxian*sin(weidu)-sin(chiweijiao2))/(cos(weidu)*sqrt(1-gaodujiaozhengxian*gaodujiaozhengxian));
		if(fangweijiaoyuxian==1)
			fangweijiao = 0;
		else
		{
			double shijiao = ShiJiao(JD,JF,NF,Y,R,S,F);//由时角的正负判断是下午还是上午，以此确定方位角的取值
			if(shijiao>=0)
				fangweijiao = acos(fangweijiaoyuxian);
			else
				fangweijiao = 2*PI-acos(fangweijiaoyuxian);
		}
	}
	
	fangweijiao = fangweijiao*180/PI;//将方位角转换为角度值
	return fangweijiao;           //返回方位角（角度值）
};

double SolarRadiation::QingXieShunShiFuZhaoDu(double JD, double JF, double NF, double Y, double R, double S, double F, double WD, double WF, 
											  double Ta, double Rh, double Pa, double Fg, int YunLeiXing, double YunLiang, double QingXieJiao, double FangWeiJiao)
{
	double gaodujiaozhengxian = GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF);
	double qingxieshunshifuzhaodu;

	if(gaodujiaozhengxian>0)
	{
		double shuipingshunshifuzhaodu = ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang);
		double daqisanshefuzhaodu = (0.15-0.1*gaodujiaozhengxian)*shuipingshunshifuzhaodu;
		
		double taiyangfangweijiao = TaiYangFangWeiJiao(JD,JF,NF,Y,R,S,F,WD,WF);
		double SF = gaodujiaozhengxian*cos(QingXieJiao*PI/180)+sqrt(1-gaodujiaozhengxian*gaodujiaozhengxian)*sin(QingXieJiao*PI/180)*cos((taiyangfangweijiao-FangWeiJiao)*PI/180);

		qingxieshunshifuzhaodu = daqisanshefuzhaodu+(shuipingshunshifuzhaodu-daqisanshefuzhaodu)*SF/gaodujiaozhengxian;
		
		if(qingxieshunshifuzhaodu<0)
			qingxieshunshifuzhaodu = 0;
	}
	else if(gaodujiaozhengxian<=0)
	{
		qingxieshunshifuzhaodu = 0;
	}
	return qingxieshunshifuzhaodu;
};

double SolarRadiation::ShuiPingPingJunFuZhaoDu(double JD, double JF, double NF, double Y, double R, double WD, double WF, double Ta, double Rh,
											   double Pa, double Fg, int YunLeiXing, double YunLiang, double T1S, double T1F, double T2S, double T2F)
{/*
	double S = T1S, F = T1F;
	double T = (T2S-T1S)*60+(T2F-T1F);
	double sum = 0;
	for(int i=0; i<T; i++)
	{
		if(GaoDuJiaoZhengXian(JD,JF,NF,Y,R,S,F,WD,WF)>0)
		{
			if(ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang)>=0)
				sum = sum+ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang);
			F = F+1;
			if(F == 60)
			{
				S = S+1;
				F = 0;
			}
		}
		else
		{
			F = F+1;
			if(F == 60)
			{
				S = S+1;
				F = 0;
			}
		}
	}
	double shuipingpingjunfuzhaodu = sum/T;
	return shuipingpingjunfuzhaodu;
*/
	double S = T1S, F = T1F;
	double n = (T2S-T1S)*60+(T2F-T1F);
	
	double sum1 = (ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,T1S,T1F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang)
		+ ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,T2S,T2F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang)
		)*0.5;
	
	double sum2 = 0;
	if(n > 1){
	//	wavelength = wavelengthbegin;
		for(int i = 1; i <=  (n - 1); i++){
		//	wavelength = wavelength*10 + 1;
		//	wavelength = wavelength / 10;
			F = F+1;
			if(F==60)
			{
				S = S+1;
				F = 0;
			}
			sum2 += ShuiPingShunShiFuZhaoDu(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang);
		}
	}
	else{}
	
	double shuipingpingjunfuzhaodu = (sum1 + sum2)*1/n;
	return shuipingpingjunfuzhaodu;
  };

double SolarRadiation::QingXiePingJunFuZhaoDu(double JD, double JF, double NF, double Y, double R, double WD, double WF, double Ta, double Rh,
											  double Pa, double Fg, int YunLeiXing, double YunLiang, double QingXieJiao, double FangWeiJiao, 
											  double T1S, double T1F, double T2S, double T2F)
{
	double S = T1S, F = T1F;
	double n = (T2S-T1S)*60+(T2F-T1F);
	
	double sum1 = (QingXieShunShiFuZhaoDu(JD,JF,NF,Y,R,T1S,T1F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang,QingXieJiao,FangWeiJiao)
		+ QingXieShunShiFuZhaoDu(JD,JF,NF,Y,R,T2S,T2F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang,QingXieJiao,FangWeiJiao)
		)*0.5;
	
	double sum2 = 0;
	if(n > 1){
		//	wavelength = wavelengthbegin;
		for(int i = 1; i <=  (n - 1); i++){
			//	wavelength = wavelength*10 + 1;
			//	wavelength = wavelength / 10;
			F = F+1;
			if(F==60)
			{
				S = S+1;
				F = 0;
			}
			sum2 += QingXieShunShiFuZhaoDu(JD,JF,NF,Y,R,S,F,WD,WF,Ta,Rh,Pa,Fg,YunLeiXing,YunLiang,QingXieJiao,FangWeiJiao);
		}
	}
	else{}
	
	double qingxiepingjunfuzhaodu = (sum1 + sum2)*1/n;
	return qingxiepingjunfuzhaodu;
};