以文本方式查看主题 - 中华农历论坛 (http://bbs.nongli.net/index.asp) -- 历法知识 (http://bbs.nongli.net/list.asp?boardid=2) ---- [原创]根据“瑞士星历表”计算太阳的出、没及上中天时刻 (http://bbs.nongli.net/dispbbs.asp?boardid=2&id=46825) |
||||||||||||
-- 作者:jyarmy -- 发布时间:2011/5/17 22:21:00 -- [原创]根据“瑞士星历表”计算太阳的出、没及上中天时刻 一、基本思路 由瑞士星历表直接能够得到太阳在某时刻的赤道坐标(或者黄道坐标),再依据坐标变换公式,将赤道坐标变换为某观察地的地平坐标,从而得到该时刻太阳的高度角和方位角。以此为基础,运用数值算法,计算某日太阳高度角达到-50’的时刻,即为日出和日落时刻,在日出和日落之间,运用三分法求凸性函数极值的方法,求出此段内太阳高度角的最大值时刻,即上中天时刻。 二、具体算法 (1)根据瑞士星历表计算太阳的赤道坐标值 直接使用瑞士星历表提供的swe_calc_ut()函数计算,其中,函数参数iflag要加入标志SEFLG_EQUATORIAL(表示使用赤道坐标)和SEFLG_RADIANS(返回值为弧度值而非角度值) (2)赤道坐标变换为地平坐标 设:A : 地平经度(即方位角,以正南为0度,正西为90度,范围0-360度); h : 地平纬度(范围-90度~90度); JD:要计算时刻的儒略日 H : 本地时角; φ :观察者地理纬度(北纬为正,南纬为负); L : 观察者地理经度(东经为负,西经为正); θ : 本地恒星时;θ0 : 格林尼治恒星时; δ : 赤纬;α : 赤经 首先计算H:H = θ – α 或 H =θ0 – L – α θ0 = 280.46061837 + 360.98564736629 * (JD - 2451545.0) + 0.000387933*T2 – T3/38710000; //单位为度 其中T = (JD - 2451545.0) / 36525 得到H结果后转为弧度值并调整至(-π,π)范围内 再根据: tan(A) = sin(H) / ( cos(H) * sin(φ) – tan(δ) * cos(φ) ) sin(h) = sin(φ) * sin(δ) + cos(φ)* cos(δ) * cos(H) 求出A和h。 其中,计算A使用反正切函数时可以直接使用c语言标准库提供的atan2(y,x)函数,得到(-π,π)范围内的结果。或者直接使用atan(x)得到(-π/2,π/2)内的结果,求出h,再利用公式cos(h)sin(A)=cos(δ)sin(H)修正A的值(+π或-π)。最后,再将结果调整至(0,2π)的范围。 (3)计算太阳高度角达到某个值的时刻 根据(2)的方法,可以得到一个函数calc,输入时间JD,可以输出太阳高度角h,现在我们要根据这个函数计算太阳达到某个高度角h0时JD的值,我使用了数值算法中最简单的二分法,二分法的前提是函数在计算的范围连续且单调,这个条件很容易实现,计算日出时刻时,二分法中的left设为当日当地的午夜零时,right为12小时后即可。 (4)计算太阳上中天时刻 从日出到日落时间范围内太阳到达高度角最大值时就是到达上中天时,此范围内太阳高度角的变化函数恰好是一个连续的凸函数,可直接运用三分法计算,其中left即日出时刻,right为日落时刻。 |
||||||||||||
-- 作者:hero_hacker -- 发布时间:2011/5/18 11:39:00 -- 我采用的计算方法是:根据瑞士星历表计算 方法: 1 用 swe_set_topo 来设置站心(观察点)位置 2 用 swe_rise_trans 来计算星体的上升, 下降, 中天 |
||||||||||||
-- 作者:jyarmy -- 发布时间:2011/5/18 13:02:00 -- 我也感觉这个方法太笨,看来还是没学透瑞士星历表啊。 另外,多谢楼上给的代码,还要多跟您学习啊,呵呵 |
||||||||||||
-- 作者:浪-淘-沙 -- 发布时间:2011/5/18 16:04:00 -- 期待楼主用编程语言(C或JAVA或JAVAScript)完成上述算法。 |
||||||||||||
-- 作者:jyarmy -- 发布时间:2011/5/18 22:06:00 -- 以下是引用浪-淘-沙在2011-5-18 16:04:00的发言:
已经完成啦,明天发上来
期待楼主用编程语言(C或JAVA或JAVAScript)完成上述算法。 |
||||||||||||
-- 作者:jyarmy -- 发布时间:2011/5/20 12:44:00 -- 附件为完整代码,里面已经包含了瑞士星历表的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; }
|
||||||||||||
-- 作者:xjw01 -- 发布时间:2011/5/26 8:27:00 -- 用三点插值法求极值,会比三分法快很多倍。 插值法在《天文算法》中有讲述。 一般只需要计算一至两次就可以得到很高的精度。 |