求由经纬度和时间计算太阳高度角和太阳方位角的公式

H当地=90°—纬度差（*同一纬度相减，异纬相加）

一般时间

太阳高度角随着地方时和太阳的赤纬的变化而变化。太阳赤纬（与太阳直射点纬度相等）以告液δ表示，观测地地理纬度用φ表示（太阳赤纬与地理纬度都是北纬为正，南纬为负），地方时(时角)以t表示，有太阳高度角的计算公式：

sin h＝sin φ sin δ＋cos φ cos δ cos t

正午时间

日升日落，同一地点一天内太阳高度角是不断变化的。日出日落时角度都为0，

正午太阳高度角

正午时太阳高度角最大，时角为0，以上公式可以简化为：

sin h＝sin φ sin δ＋cos φ cos δ

由两角和与差的三角函数公式，可得

sin h＝cos(φ－δ)

因此，

对于太阳位于天顶以北的地区而言，h＝90°－(φ－δ)；

对于太阳位于天顶以南的地区而言，h＝90°－(δ－φ)；

二者合并，因为无论是(φ－δ)还是(δ－φ)，都是为了求当地蚂皮纬度与太阳直射纬度之差袜物物，不会是负的，因此都等于它的绝对值，所以正午太阳高度角计算公式：

h=90°-|φ-δ|

{

float dY

float dX

float Julian_Day

long int liAux1

long int liAux2

float Omega

float MeanLongitude

float MeanAnomaly

float Sin_EclipticLongitude

float GreenwichMeanSiderealTime //Prime meridian

float LocalMeanSiderealTime

float LatitudeInRadians

float HourAngle

float Cos_Latitude

float Sin_Latitude

float Cos_HourAngle

t.hour=t.hour+4

HoursinDecimal=(t.hour + t.min/ 60.0)

liAux1 =(t.mon-14)/12

liAux2=(1461*(t.year + 4800 + liAux1))/4 +(367*(t.mon- 2-12*liAux1))/12- (3*((t.year + 4900+ liAux1)/100))/4+t.date-32075

Julian_Day=(float)(liAux2)-0.5+HoursinDecimal/24.0 //Calculate current Julian Day

ElapsedJulianDays =Julian_Day-2451545.0//Calculate difference between current Julian Day and January 1st noon2000 Universal Time—JD 2451545.0

Omega=2.1429-0.0010394594*ElapsedJulianDays

MeanLongitude = 4.8950630+0.017202791698*ElapsedJulianDays //Inradians

MeanAnomaly = 6.2400600+0.0172019699*ElapsedJulianDays //Inradians

EclipticLongitude = MeanLongitude +0.03341607*sin( MeanAnomaly ) + 0.00034894*sin( 2*MeanAnomaly)-0.0001134-0.0000203*sin(Omega)//Calculate //ecliptic longitude

EclipticObliquity=0.4090928-6.2140e-9*ElapsedJulianDays+0.0000396*cos(Omega)

//Calculateecliptic obliquity

Sin_EclipticLongitude= sin( EclipticLongitude)

dY = cos( EclipticObliquity ) *Sin_EclipticLongitude

dX = cos( EclipticLongitude )

RightAscension = atan2( dY,dX ) //Calculatecelestial coordinates, right ascension in radians

if( RightAscension <0.0 )

{

RightAscension = RightAscension + twopi

}

Declination = asin( sin( EclipticObliquity)*Sin_EclipticLongitude )

//Calculate celestial coordinates, declination in radians

//Calculate local coordinates (azimuth andzenith angle) in degrees

GreenwichMeanSiderealTime= 6.6974243242 +0.0657098283*ElapsedJulianDays + HoursinDecimal

LocalMeanSiderealTime =(GreenwichMeanSiderealTime*15+ Longitude)*rad

HourAngle = LocalMeanSiderealTime -RightAscension

LatitudeInRadians = Latitude*rad

Cos_Latitude = cos( LatitudeInRadians )

Sin_Latitude = sin( LatitudeInRadians )

Cos_HourAngle= cos( HourAngle )

dY = -sin( HourAngle )

dX = tan( Declination )*Cos_Latitude -Sin_Latitude*Cos_HourAngle

Azimuth = atan2( dY, dX )//Calculatelocal coordinates, azimuth in degrees

if ( Azimuth <0.0 )

{

Azimuth = Azimuth + twopi

}

Azimuth = Azimuth/rad

ZenithAngle=(acos(Cos_Latitude*Cos_HourAngle*cos(Declination)+sin(Declination)*Sin_Latitude))