附件为完整代码，里面已经包含了瑞士星历表的dll文件，可以自行下载它的数据文件，结果会精确一些。VC6.0编译通过~
下面是核心代码：
//根据赤道坐标计算地平坐标
void calculate(struct input *in, struct output *out )
{
char serr[AS_MAXCH];
double x1[6];
long iflag, iflgret;
int p;
double GMST;//格林尼治恒星时
double H;//时角
double T;
double result1;//方位角
double result2;//高度角
double c_lon; //赤经
double c_lat; //赤纬
p=SE_SUN;
iflag =0;
iflag |= SEFLG_EQUATORIAL;//赤道坐标系
iflag |= SEFLG_RADIANS; //返回值为弧度值，而非角度值
iflgret = swe_calc_ut(in->tjd_ut , p, iflag, x1, serr);
//得到太阳的赤经和赤纬
if (iflgret < 0){
printf("error: %s\n", serr);
return;
}
c_lon = x1[0];
c_lat = x1[1];
T=(in->tjd_ut - 2451545.0) / 36525 ;
GMST = 280.46061837 + 360.98564736629*(in->tjd_ut-2451545.0)+ 0.000387933*T*T - T*T*T/38710000;//单位为度
H = d2r(GMST) - in->lon - c_lon ; //时角
if(H<0){
while(H<(-1)*PI)
H += 2*PI;
}
else{
while(H>PI)
H -=2*PI;
}
result2 = asin( sin(in->lat) * sin(c_lat) +cos(in->lat)*cos(c_lat)*cos(H) );//得到高度角
result1 = atan2( sin(H) , cos(H) * sin(in->lat) - tan(c_lat) * cos(in->lat) );//得到方位角,此方位角为由南向西测量，即正南为0度，正西为90度
if(result1<0)
result1 += 2*PI;//调整为0 - 360度的范围内
out->azimuth = result1;
out->h = result2;
}
//从start的时间开始计算，返回太阳到达高度h的时间
//主要用于计算日出及日落的时刻
//计算日出时刻一般从当地当日零时开始，计算日落从当地正午开始
double WSunHight(struct input *start , double h)
{
struct output hc;//horizon coordinate地平坐标
struct input *in,in0;
double step;
int sign,multiple;
int i;
in=&in0;
in0.lat = start->lat; in0.lon = start->lon; in0.tjd_ut = start->tjd_ut;
step = 5.0/(24.0*60.0);//5分钟
calculate(in,&hc);
if(h > hc.h)//目标高度比当前高度大 用于计算日出时刻
{
for(i=0;i<12*60/step;i++){ //一般只在12小时内搜寻
in->tjd_ut += step;
calculate(in,&hc);
if(hc.h
continue;
else
{
multiple = 2;
in->tjd_ut -= step/multiple;
for(;;){
calculate(in,&hc);
if(hc.h < h ){
if( (h-hc.h) < PRECISION )
break; //get it!
else
sign=1;
}
else{
if( (hc.h-h) < PRECISION )
break; //get it
else
sign=-1;
}
multiple *=2;
in->tjd_ut += (step/multiple)*sign;
}//for
break;
}
}
}
else{//目标高度比当前高度小 用于计算日落时刻
for(i=0;i<12*60/step;i++){ //一般只在12小时内搜寻
in->tjd_ut += step;
calculate(in,&hc);
if(hc.h > h)
continue;
else
{
multiple = 2;
in->tjd_ut -= step/multiple;
for(;;){
calculate(in,&hc);
if(hc.h > h ){
if( (hc.h-h) < PRECISION )
break; //get it!
else
sign=1;
}
else{
if( (h-hc.h) < PRECISION )
break; //get it
else
sign=-1;
}
multiple *=2;
in->tjd_ut += (step/multiple)*sign;
}//for
break;
}
}
}
return in->tjd_ut;
}
//用“三分法”确定一段时间内太阳高度角极值，及此时的时间
//in标识出起始计算时刻和地理经纬度，
//length表示计算的时间段长度，在此段时间内太阳应达到最高点
//求出的最大高度角存储在toph中，返回此时的时间
double GetTop(struct input *in, double length, double *toph)
{
struct input temp;
struct output hc;
double left,right,mid,midmid;
double v_left, v_right, v_mid, v_midmid;
double x;
left=in->tjd_ut;
right= left+length;
temp.lat = in->lat;
temp.lon = in->lon;
for(;;){
temp.tjd_ut=left;
calculate(&temp,&hc); v_left=hc.h; //得到v_left
temp.tjd_ut=right;
calculate(&temp,&hc); v_right=hc.h;//得到v_right
mid=(left + right)/2.0;
temp.tjd_ut=mid;
calculate(&temp,&hc); v_mid=hc.h; //得到v_mid
midmid=(mid+ right)/2.0;
temp.tjd_ut=midmid;
calculate(&temp,&hc); v_midmid=hc.h;//得到v_midmid
x = ( v_left < v_mid )?(v_mid - v_left):( v_left - v_mid );//取v_mid与v_left的差
if( x< PRECISION ) //满足精度要求了
break;
if( v_mid > v_midmid )
right = midmid;
else
left = mid;
}
//
*toph=v_left;
return left;
}
加好友 
发短信
Post By：2011/5/17 22:21:00
好评，获得10个金币奖励
