[转帖]太阳位置计算程序
把它存成*.vbs运行
'定义常数圆周率π

Const pi = 3.14159265359

'定义弧度与角度的转换比例

Const radeg = 57.29577951307

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

'工具函数：把任意角度化为0-360度范围内的角度

Function in360(bigdeg)

While bigdeg < -3600
   bigdeg = bigdeg + 3600

Wend

While bigdeg >3600
   bigdeg = bigdeg - 3600

Wend

While bigdeg < 0
   bigdeg = bigdeg + 360

Wend

While bigdeg >=360
   bigdeg = bigdeg - 360

Wend

in360 = bigdeg

End Function

'工具函数：角度转弧度

Function toRad(deg)
  toRad = deg / radeg

End Function

'工具函数：弧度转角度

Function toDeg(rad)
  toDeg = rad * radeg

End Function

'工具函数：角度三角函数Sin

Function fnsind(x)
  fnsind = Sin(x/radeg)

End Function

'工具函数：角度三角函数Cos

Function fncosd(x)
  fncosd = Cos(x/radeg)

End Function

'工具函数：角度三角函数Tan

Function fntand(x)
  fntand = Tan(x/radeg)

End Function

'工具函数：角度反正弦

Function fnasind(x)
  fnasind = radeg * Atn(x/Sqr(1-x*x))

End Function

'工具函数：角度反余弦

Function fnacosd(x)
  fnacosd = 90 - fnasind(x)

End Function

'工具函数：角度反正切

Function fnatand(x)
  fnatand = radeg * Atn(x)

End Function

'工具函数：角度反正切任意象限

Function fnatan2d(y,x)
  If x = 0 Then
          fnatan2d = 90
  End If
  If x < 0 Then
    fnatan2d = radeg * Atn(y/x) - 180
  End If
  If x >0 Then
          fnatan2d = radeg * Atn(y/x)
  End If

End Function

'工具函数：弧度反正切任意象限

Function fnatan2(y,x)
  If x = 0 Then
          fnatan2 = pi/2
  End If
  If x<0 Then
    fnatan2 = Atn(y/x) - pi
  End If
  If x>0 Then
          fnatan2 = Atn(y/x)
  End If

End Function

'工具函数：立方根

Function Cbrt(x)
  If x > 0.0 Then
    Cbrt = Exp(Log(x)/3)
  Else
          If x < 0.0 Then
           Cbrt = -Cbrt(-x)
    Else
           Return 0.0
          End If
  End If

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

'把赤经化成小时、分、秒的常用形式

Function printRA(RA)

RA = in360(RA)

cHour = RA15

cMinute = Int((RA - cHour*15)/0.25)

cSecond = Round((RA - cHour*15 - cMinute*0.25)/(0.25/60),1)

printRA = cHour&"h"&cMinute&"m"&cSecond&"s"

End Function

'把赤纬化成度、分、秒的区分南北纬的常用形式

Function printDecl(Decl)

If Decl >= 360 Or Decl < 0 Then
    Decl = in360(Decl)

End If

If Decl >90 And Decl < 270 Then
        Decl = 180 - Decl

End If

If Decl >=270 And Decl < 360 Then
        Decl = Decl - 360

End If

If Decl < 0 Then
        flag = "-"

Else
        flag = "+"

End If

Decl = Abs(Decl)

arcDeg = Int(Decl)

arcMin = Int((Decl - arcDeg)*60)

arcSec = Round((Decl - arcDeg - arcMin/60)*3600,0)

printDecl = flag&arcDeg&"°"&arcMin&"'"&arcSec&Chr(34)

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

'计算一个给定的世界时日期与UT2000.0的差值

Function dayFrom2000(theYear,theMonth,theDate)

dayFrom2000 = 367*theYear-7*(theYear+(theMonth+9)12)4 + (275*theMonth)9 + theDate - 730530

End Function



'从赤经赤纬距离计算赤道平面直角坐标

Function EquatRectCoord(RA,Decl,dist)
  xequat = dist * Cos(RA) * Cos(Decl)
  yequat = dist * Sin(RA) * Cos(Decl)
  zequat = dist * Sin(Decl)
  EquatRectCoord = Array(xequat,yequat,zequat)

End Function



'从赤道平面直角坐标来计算赤道球面坐标

Function EquatSpherCoord(xequat,yequat,zequat)
  r = Sqr(xequat*xequat + yequat*yequat + zequat*zequat)
  RA = toDeg(fnatan2(yequat,xequat))'天球赤道经度
  Decl = toDeg(fnatan2(zequat,Sqr(xequat*xequat+yequat*yequat)))'天球赤道纬度
  EquatSpherCoord = Array(r,RA,Decl)

End Function



'球面赤道坐标换算为对于某地的地平坐标

'RA:赤经，Decl:赤纬,SiderealT:恒星时，geoLat:地理纬度,geoLon:地理经度,ASL:海拔高度

Function HorizonalCoord(RA,Decl,SiderealT,geoLat,geoLon,ASL)

'@todo Not completed, wait to implement

End Function



Sub Main
  Call SunPosition

End Sub



Sub SunPosition '计算太阳的位置

Dim today

today = CDate(now)

todayYear = DatePart("yyyy",today)

todayMonth = DatePart("m",today)

todayDate = DatePart("d",today)



d = dayFrom2000(todayYear,todayMonth,todayDate)'距历元的间隔

w = 282.9409 + 4.70935E-5 * d '地球的近日点日心黄经。

a = 1.0000000 '地球与太阳的平均距离，单位为天文单位au

eqe = 0.016709 - 1.151E-9 * d '地球公转轨道偏心率

M = 356.0470 + 0.9856002585 * d '太阳平均运动

oblecl = 23.4393 - 3.563E-7 * d '地球的黄赤交角

L = w + M '太阳的平均黄道经度



'解开普勒方程，求偏近点角

E = M + radeg * eqe * fnsind(M) * (1 + eqe * fncosd(M))

E = M + radeg * eqe * fnsind(E)



'黄道平面近日点直角坐标

x = fncosd(E) - eqe

y = fnsind(E) * Sqr(1 - eqe*eqe)



'地日距离

r = Sqr(x*x + y*y)



'近日点角距

v = fnatan2d(y,x)



lon = v + w'太阳的黄道经度

lon = in360(lon)'将黄道经度化归到360度角以内



xecliptic = r * fncosd(lon)'将太阳黄道经度化为黄道直角坐标

yecliptic = r * fnsind(lon)

zecliptic = 0.0



xequat = xecliptic'将黄道直角坐标旋转到赤道直角坐标

yequat = yecliptic * fncosd(oblecl) - zecliptic*fnsind(oblecl)

zequat = yecliptic*fnsind(oblecl) + zecliptic*fncosd(oblecl)



spher_Sun = EquatSpherCoord(xequat,yequat,zequat)'通过赤道直角坐标计算太阳球面坐标





MsgBox("the Day "&todayYear&"-"&todayMonth&"-"&todayDate&"'s lon is" & lon)

MsgBox("r="&r&";RA="&printRA(spher_Sun(1))&";Decl="&printDecl(spher_Sun(2)))

End Sub



Call Main