月上柳梢头人约黄昏后”的天文解释.docx

第 3 页 共 51 页参赛报名号2015075 “月上柳梢头，人约黄昏后”的天文解释摘要 月球，是环绕地球运行的一颗卫星，是被人们研究得最彻底的天体。在古时，诗人常常借着月亮表达自己的情感。“月上柳梢头，人约黄昏后”正是诗人欧阳修在元宵之时与佳人相约之景的写照，富有诗情画意。在天文学中，月亮的升起和太阳的落下有着不同的运动规律，而欧阳修当时所遇之景并非在何时何地均可发生。当太阳缓缓地从地平线上落下，月亮已升起至柳梢之上正是这句诗的表面解释，而要想更精确地确定月亮和太阳的位置，我们可用月亮与太阳的高度角来对这一情景进行诠释。对于问题一，在本文中，我们先对柳树和人的高度与人和柳树间的距离大致确定月亮在天空中的角度，然后考虑到黄昏是指日落以后到天还没有完全黑的这段时间，我们将运用天文学的计算理论计算太阳高度角以确定日落时间从而定义黄昏的时间。在模型的假设条件完成之后，我们查阅相关天文理论的资料，用ELP-2000/82 月球理论计算月球位置，通过球面三角知识与月球位置求得的月球高度角，利用c语言的for语句可得特定日期“月上柳梢头”的时间范围，而通过太阳高度角计算的相关公式又可以算得某地点的日落时间与天黑时间从而确定“人约黄昏后”的时间范围。在这些计算中，我们将使用c语言进行数据的运算。对于问题二，我们根据第一问所建立的模型用编写的程序求解2016某个地点每天的日落时间和天黑时间确定黄昏的时间段，然后先将黄昏时间每隔10分钟取一个时间点，利用c语言的for语句将每个时间点所对应的月亮高度角计算出来，如果月亮的高度角符合“月上柳梢头”的角度范围，且月亮处于上升状态，则会输出此时的日期和时间。最后，我们通过此程序得知一年中各地发生此现象的所有日期和时间。关键词：c语言 太阳高度角 月亮高度角 一、问题重述“月上柳梢头，人约黄昏后”是北宋学者欧阳修的名句，写的是与佳人相约的情景。请用天文学的观点赏析该名句，并进行如下的讨论：1、定义“月上柳梢头”时月亮在空中的角度和什么时间称为“黄昏后”。根据天文学的基本知识，在适当简化的基础上，建立数学模型，分别确定“月上柳梢头”和“人约黄昏后”发生的日期与时间，由于大部分地点（除了极昼和极夜的地区）都能分别发生以上两个现象，因此我们计算某个特定的日期和地点分别发生这两个现象的时间，并根据已有的天文资料（如太阳和月亮在天空中的位置、日出日没时刻、月出月没时刻）验证所建模型的合理性。2、根据所建立的模型，分析2016年北京地区“月上柳梢头，人约黄昏后”发生的日期与时间。根据模型判断2016年在哈尔滨、上海、广州、昆明、成都、乌鲁木齐是否能发生这一情景？如果能，请给出相应的日期与时间；如果不能，请给出原因。二、模型假设1. 由于人和柳树的高度及人与柳树间的距离具有不确定性，不妨假设人的高度为h1 =1.7米，柳梢头的高度为h2=8米，人与柳树的距离在25-50米。2. 假设黄昏从太阳高度角降至0开始。3. 假设人眼的水平视线与柳梢头的距离差即为人和柳树高度的差值。4. 忽略地球章动与光行差。5. 不考虑月食的出现。6. 忽略地形的影响且观测条件良好。三、符号说明符号意义黄经黄纬JD儒略日观测者纬度观测者经度H本地时角黄赤交角本地恒星时赤经赤纬UT北京时间h地平纬度（即高度角）T儒略世纪数N积日四、问题分析4.1诗文赏析“月上柳梢头，人约黄昏后”这句诗源自著名诗人欧阳修脍炙人口的诗篇生查子元夕，此诗情景交融，展现了美的意境。而诗中的“去年元月时”告诉了我们当时月上柳梢头的时间。其中“月上柳梢头，人约黄昏后”这句诗蕴含着天文知识。4.2问题一的分析题目首先要求定义“月上柳梢头”时月亮在空中的角度和黄昏后的时间范围，而由模型假设可以知道月上柳梢头时月亮所对应的最小高度角的正切值tan(t1)= ，最大高度角的正切值tan(t2)= ，由此可得t1 ，t2因此定义月亮在空中的高度角为918，即月亮的高度角为918。因为在太阳高度角降为0之后，天空渐渐逐步变暗，而天文学上定义日落与天黑之间的时间为蒙影，且蒙影的民用太阳高度角为-60。因此我们不妨定义太阳高度角在-60时为黄昏以确定黄昏的时间段。在以上两个范围确定之后，我们将利用天文学中的理论建立数学模型，进而能够确定某经纬度的月亮高度角和太阳高度角与日期和时间的关系从而分别确定“月上柳梢头”和“人约黄昏后”对应的时间。在此计算过程中，我们将根据模型编写c语言程序进行第 5 页 共 51 页参赛报名号2015075 第3页共14页较为精确的计算。4.3问题二的分析题目要求我们分析2016年北京发生“月上柳梢头，人约黄昏后”的日期和时间，这时由第一问中建立的模型可以先通过计算黄昏后中某些时间点时月亮所对应的月亮高度角进而确定月亮的高度角是否符合之前定义的范围。如果符合，且月亮处于上升状态，我们即找到了所要求得日期和时间。同样的道理，对于哈尔滨等地用相同的方法可以使用程序得知可否输出日期和时间，进而判断是否可发生此现象并给出最后的结果。五、模型建立和解决5.1问题的模型建立和解决5.1.1模型的准备查阅相关的天文知识，我们在计算过程中将用到儒略日，恒星时，月球黄经黄纬，因此在此引用相关资料先对三者进行说明。【2】儒略日： 儒略日是指从公元-4712年开始连续计算日数得出的天数及不满一日的小数，通常记为JD(*)。传统上儒略日的计数是从格林尼治平午，即世界时12点开始的。若以力学时（或历书时）为标尺，这种计数通常表达为“儒略历书日”，即JDE，中E只是一种表征，即按每天86400个标准秒长严格地计。Y为给定年份，M为月份，D为该月日期（可以带小数）。若M2，Y和M不变，若M=1或2，以Y1代Y，以M+12代M，换句话说，如果日期在1月或2月，则被看作是在前一年的13月或14月。A=INT（Y/100）B=2-A+INT(A/4)JD=INT(365.25(Y+4716)+INT(30.6001(M+1)+D+B-1524.5 INT为取整函数【2】恒星时：恒星时是根据地球自转来计算的，它的基础是恒星日。本地恒星时的定义是一个地方的子午圈与天球的春分点之间的时角。本文所述的时间指世界时UT，0hUT表示世界时0点。对于一个给定的日期，格林尼治子午圈的恒星时可按如下计算。先计算当天（0hUT）对应的儒略日JD，它是以.5结束的数字。那么就有：T=(JD-2451545.0)/36525那么(0hUT)的格林尼治平恒星时使用以下表达式计算为了找出任意UT时刻的恒星时，先计算UT*1.00273790935，再把计算结果加入(0hUT)时刻对应的平恒星时，即得到UT时刻的平恒星时。对于北京时间，也就是东经120度时间UT，任意地点的恒星时有恒星时不是时间单位，这里单位是度。为当地经度，东经为正，西经为负。第 51 页 共 51 页参赛报名号2015075 【1】【2】月球位置为了准确计算出某时刻月球的准确位置，须计算月球黄经黄纬及距离的数百个周期项，这里仅考虑主要的周期项。利用本算法，可得到月心位置坐标：黄经()、黄纬()。一、计算方法:本文的周期项是基于ELP-2000/82月球理论。T表达为J2000起算的儒略世纪数。使用以下表达式计算角度L，D，M，M，F，角度单位是度。月球平黄经：L=218.3164591+481267.88134236T-0.0013268T2+T3/538841-T4/65194000月日距角：D=297.8502042+445267.1115168T-0.0016300T2+T3/545868-T4/113065000太阳平近点角：M=357.5291092+35999.0502909T-0.0001536T2+T3/24490000月亮平近点角：M=134.9634114+477198.8676313T+0.0089970T2+T3/23569699-T4/14712000月球经度参数：F=93.2720993+483202.0175273T-0.0034029T2-T3/3526000+T4/863310000三个必要的参数：A1=119.75+131.849TA2=53.09+479264.290TA3=313.45+481266.484T取和计算45.A表中各项(I)，取和计算45.B表中各项(b)。这里I与b是正弦项取和，A*sin()，式中的是表中D、M、M、F的线性组合，组合系数在表45.A及45.B相应的列中，A是振幅。以表45.A第8行为例：同理可计算第1、2、3、4.各行,得到I1、I2、I3.最后I=I1+I2+I3+.。然而,表中的这些项包含了了M(太阳平近点角)，它与地球公转轨道的离心率有关，就目前而言离心率随时间不断减小。由于这个原因，振幅A实际上是个变量(并不是表中的常数)，角度中含M或-M时，还须乘上E，含2M或-2M时须乘以E的平方进行修正。E的表达式如下：E=1-0.002516T-0.0000074T2此外，还要处理主要的行星摄动问题（A1与金星摄动相关，A2与木星摄动相关，L与地球扁率摄动相关）：I+=+3958*sin(A1)+1962*sin(L-F)+318*sin(A2)b+=-2235*sin(L)+382*sin(A3)+175*sin(A1-F)+175*sin(A1+F)+127*sin(L-M)-115236*sin(L+M)最后得到月球的坐标如下：=L+I/1000000(黄经单位：度)=b/1000000(黄纬单位：度)因45.A及45.B表中的振幅系数的单位是10-6度，所以上式计算时除以1000000 二、两个计算用的表:表45.A月球黄经周期项(I)黄经单位：0.000001度，角度的组合系数I的各项振幅A(正弦振幅)DMMF0010628874420-101274027200065831400202136180100-1851160002-11433220-20587932-1-10570662010533222-1004575801-10-409231000-347200110-30383200-2153270012-12528001-21098040-101067500301003440-20854821-10-78882100-676610-10-5163110049872-1104036202039944000386120-30366501-20-268920-12-26022-1-2023901010-23482-20022360120-21200200-20692-2-102048201-2-17732002-15954-1-1012150022-111030-10-8922110-8104-1-2075902-10-71322-10-70021-206912-10-2596401054900405374-10052010-20-487210-2-399002-2-381111035130-20-34040-303302-1203270210-32311-10299203029420-1-20 表 45.B 月球黄纬周期项(b)，单位：0.000001 度 - 角度的组合系数b的各项振幅A(正弦振幅)DMMF000151281220011280602001-1277693200-117323720-115541320-1-146271200132573002117198201-19266002-188222-10-1821620-2-1432420114200210-1-33592-1-1124632-10122112-1-1-1206501-1-1-187040-1-118280101-17940003-174901-11-15651001-14910111-1475011-1-1410010-1-1344100-1-133500311107400-1102140-11833001-377740-21671200-3607202-15962-11-149120-21-451003-1439202142220-3-142121-11-3662101-35140013312-1113152-20-13020013-283211-1-229110-1223110122301-2-1-22021-1-1-2201011-1852-1-2-11810121-17740-2-11764-1-1-1166101-1-164401-113210-1-1-1194-10-11152-2011075.1.2模型的建立 为了帮助理解，我们先将月球高度角求解的流程图表示出来 月亮高度角计算首先，计算月亮高度角需要先求出月球的黄经和黄纬，根据上面的准备， (黄经单位：度) (黄纬单位：度) 接着，通过坐标转换的知识将黄经、黄纬转换为赤经、赤纬 然后，再利用时间和经纬度的数据可得月亮高度角h为本地恒星时，为赤经【3】太阳高度角的计算N0= 79.6764+0.2422*(年份-1985)-INT(年份-1985)/4M=N-N0 N为积日，即这天在一年中的天数，比如2月1日，N=32。日角的计算公式太阳赤纬表达式 太阳高度角的表达式t为当地太阳时，即当地时间。计算日落时间时，令h=0。计算天黑时间时，令h=-6，可求出相应时间t。北京时间与当地时间关系，为当地经度。 5.1.3模型的解决 第一问：利用建立的模型编写程序，用for语句搜索当天各时间的太阳和月亮高度角，筛选符合要求的时间。下面计算北京（116E，40N）2016年1月22日的“月上柳梢头”和“人约黄昏后”的时间。 “月上柳梢头”时间：17时-18时 “人约黄昏后”开始时间：17时23分第二问：先计算2016年每一天的黄昏时间，再由第一问中建立的模型可以先通过计算黄昏内某些时间点月亮所对应的月亮高度角进而确定月亮的高度角是否符合之前定义的范围，最后判断月亮是否处于上升状态，如果符合，我们即找到了所要求得日期和时间。判断月亮是否为上升状态的方法：计算某时刻t的月亮高度角h1与此时刻一小时后的时刻t+1的月亮高度角h2，若h2h1，则月亮上升状态。因为1小时内月亮无法从18运行到月中天，所以月亮先上升后下降的情况不存在，任何情况此方法有效。，通过程序的运算得到北京（116E，40N）2016年发生“月上柳梢头，人约黄昏后”的时间为2016年2月21日17时53分,10.888885度2016年2月22日18时4分,11.725196度2016年3月22日18时24分,14.006413度2016年4月20日18时51分,16.229372度2016年4月21日19时2分,17.129314度2016年5月21日19时14分,9.402243度2016年6月19日19时28分,11.637392度2016年7月18日19时27分,14.651909度2016年8月17日19时4分,13.086981度2016年9月14日18时24分,17.040792度2016年9月15日18时32分,10.973973度2016年10月13日17时37分,15.775666度2016年10月14日17时45分,9.556035度2016年11月13日17时4分,12.350664度2016年12月11日16时59分,17.109592度2016年12月13日17时9分,11.384358度查阅日历后，我们发现2016年2月22日是元宵节，且其他日期都在农历十三至十六，这与欧阳修诗中提及的元宵节相符合，因此我们的模型还是基本正确的。继续运行程序计算其他六个地点2016年发生“月上柳梢头，人约黄昏后”的时间。哈尔滨（126E，45N）2016年1月22日16时31分,9.239304度2016年2月21日17时6分,11.653321度2016年3月21日17时43分,12.957748度2016年3月22日18时4分,13.519439度2016年4月20日18时19分,14.803655度2016年5月19日18时45分,15.843122度2016年5月20日19时6分,15.968928度2016年6月18日19时2分,16.635395度2016年7月17日19时1分,17.840364度2016年7月18日19时10分,11.376996度2016年8月15日18时36分,13.608990度2016年9月12日17时50分,16.100756度2016年9月14日18时6分,14.059811度2016年10月11日16时55分,14.197354度2016年10月12日17时4分,17.999767度2016年10月13日17时2分,13.525196度2016年11月11日16时14分,13.309301度2016年12月11日16时5分,16.904108度2016年12月12日16时15分,9.858195度上海（121E，31N）2016年1月22日17时21分,17.448652度2016年2月21日17时43分,9.311610度2016年2月22日18时4分,10.825257度2016年3月22日18时4分,14.608612度2016年5月20日18时36分,11.731987度2016年5月21日19时6分,13.350939度2016年6月18日18时47分,15.249158度2016年6月19日19时7分,16.624075度2016年7月18日18时48分,9.599825度2016年7月19日19时8分,11.712987度2016年8月16日18时33分,15.078372度2016年9月15日18时9分,14.541442度2016年10月14日17时23分,11.612833度2016年11月13日17时0分,12.299422度2016年12月11日16时59分,17.141016度2016年12月13日17时0分,9.362485度广州（113E，23N）2016年1月23日18时7分,17.243105度2016年2月22日18时23分,9.800927度2016年3月22日18时36分,14.836303度2016年3月23日19时6分,17.572811度2016年4月21日18时45分,9.247188度2016年4月22日19时6分,11.952693度2016年5月20日18时55分,14.397724度2016年5月21日19时5分,16.632683度2016年6月20日19时4分,10.518424度2016年7月19日19时6分,15.593328度2016年8月17日18时54分,9.478863度2016年8月18日19时4分,12.672923度2016年9月15日18时29分,17.451839度2016年10月15日18时9分,16.997084度2016年11月13日17时44分,11.995618度2016年11月14日18时3分,14.646480度昆明（103E，25N）2016年1月22日18时43分,16.181719度2016年1月23日19时4分,17.653614度2016年2月22日19时2分,10.075622度2016年3月22日19时16分,14.806972度2016年5月20日19时38分,13.753357度2016年5月21日20时8分,15.836126度2016年6月20日20时8分,9.660623度2016年7月18日19时50分,12.424285度2016年7月19日20时0分,14.643270度2016年8月18日20时6分,11.982417度2016年9月15日19时9分,16.750822度2016年10月14日18时39分,12.833344度2016年10月15日19时8分,16.738533度2016年11月13日18时21分,12.094069度成都（104，31N）2016年1月22日18时29分,17.448652度2016年2月21日18时51分,9.311610度2016年2月22日19时2分,10.825257度2016年3月22日19时12分,14.608612度2016年5月20日19时44分,11.731987度2016年5月21日20时4分,13.350939度2016年6月18日19时55分,15.249158度2016年6月19日20时5分,16.624075度2016年7月18日19时56分,9.599825度2016年7月19日20时6分,11.712987度2016年8月16日19时41分,15.078372度2016年8月17日20时0分,17.952955度2016年8月18日20时9分,9.836185度2016年9月15日19时7分,14.541442度2016年10月14日18时31分,11.612833度2016年11月13日18时8分,12.299422度2016年12月13日18时8分,9.362485度乌鲁木齐（88E，44N）2016年1月23日19时6分,9.770240度2016年2月21日19时39分,11.507284度2016年2月22日20时0分,12.033051度2016年3月22日20时16分,13.625275度2016年4月20日20时49分,15.097610度2016年4月21日21时0分,15.590712度2016年5月20日21时16分,16.527796度2016年6月18日21时31分,17.322170度2016年6月19日21时41分,9.351733度2016年7月18日21时30分,12.036457度2016年8月16日21时5分,16.322237度2016年8月17日21时13分,10.845237度2016年9月14日20时18分,14.662129度2016年9月15日20时26分,9.310205度2016年10月13日19时25分,13.982441度2016年11月11日18时48分,13.592419度2016年11月13日19时7分,12.274093度2016年12月11日18时40分,16.955854度2016年12月13日19时0分,12.198712度六、模型的评价6.1问题一的评价： 问题一的主要难题是建立时间地点与太阳月亮高度角的关系。我们采用了天文资料上比较精确的算法，所输出结果与“寿星万年历”软件的结果相比，误差很小，因此判定我们的模型基本正确且精度较高。6.2问题二的评价：问题二的主要问题是设置循环算法，利用计算机搜索符合题目要求的日期和时间。一方面，我们利用“stellarium”星图软件对输出结果进行验证，发现模拟的场景基本符合题目要求。另一方面，输出日期皆在农历十五前后三天内，与诗中提及的元宵佳节相吻合。由此可见，本模型基本正确且精度较高。七、参考文献 1 Jean.Meeus,Astronomical.Algorithms,America:Willmann Bell,1991 2许建伟，天文算法，/p-309119695.html，2015年10月25日 3未知作者，太阳高度角和太阳方位角，/link?url=1fFyf9wBQPosyhtej0Rw1kNKr7xjNulnskbVlv10LOaoUATR-OqWy6PoL7XGRMkdOdqGc93qg0rwbpN8XXY0vjSO-8wVYJoRPHobQVnHzP3，2015年10月25日附录: （注意，程序中经纬度输入格式为“经度，纬度”，东经为正西经为负，北纬为正南纬为负，不带E W N S标识。）附录计算1月22日日落时间源代码#include#include#define PI 3.1415926#define DE 3.141592654/180 double osin(double can)float can2,shuchu1;/can2=can/180*PI;shuchu1=sin(can2);return shuchu1;double otan(double can7)double can8,shuchu4;can8=can7/180*PI;shuchu4=tan(can8);return shuchu4;double ocos(double can5)double can6,shuchu3;can6=can5/180*PI;shuchu3=cos(can6);return shuchu3;double oasin(double can3)double can4,shuchu2;can4=asin(can3);shuchu2=can4/PI*180;return shuchu2;double oacos(double can12)double can14,shuchu8;can14=acos(can12);shuchu8=can14/PI*180;return shuchu8;double oatan(double can9)double can10,shuchu5;can10=atan(can9);shuchu5=can10/PI*180;return shuchu5;/C中三角函数使用弧度，所以先化弧度为角度。int getyear(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return year;/儒略日化年份int getmonth(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return month;/儒略日化月份int getday(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return day;/儒略日化日期float suntime(float x,float y,int days)float t,t2,t3,ED,z,N,ET;#define PI 3.1415926N=79.6764+0.2422*(getyear(days)-1985)-floor(getyear(days)-1985)/4);t=days-2457389-N;z=360*t/365.2422;ED=0.3723+23.2567*osin(z)+0.1149*osin(2*z)-0.1712*sin(3*z)-0.758*ocos(z)+0.3656*ocos(2*z)+0.0201*ocos(3*z);/太阳赤纬ET=0.0028-1.9857*osin(z)+9.9059*osin(2*z)-7.0924*ocos(z)-0.6882*ocos(2*z);t2=(oacos(-osin(y)*osin(ED)/ocos(y)*ocos(ED)/15+(120.0-x)/15+12-ET/60;/落日时间return t2;void main()float time,x,y;int Month,Day,Hour,Min,Year;scanf(%f,%f,&x,&y);for(int days=2457389;days=2457754;)time=suntime(x,y,days);Year=getyear(days);Month=getmonth(days);Day=getday(days);Hour=floor(time);Min=floor(time-floor(time)*60);printf(日期是%d年%d月%d日%d时%d分n,Year,Month,Day,Hour,Min);days=days+1;附录问题一：计算1月22日月出月落时间源代码#include#include#define PI 3.1415926#define DE 3.141592654/180 double osin(double can)float can2,shuchu1;/can2=can/180*PI;shuchu1=sin(can2);return shuchu1;double otan(double can7)double can8,shuchu4;can8=can7/180*PI;shuchu4=tan(can8);return shuchu4;double ocos(double can5)double can6,shuchu3;can6=can5/180*PI;shuchu3=cos(can6);return shuchu3;double oasin(double can3)double can4,shuchu2;can4=asin(can3);shuchu2=can4/PI*180;return shuchu2;double oacos(double can12)double can14,shuchu8;can14=acos(can12);shuchu8=can14/PI*180;return shuchu8;double oatan(double can9)double can10,shuchu5;can10=atan(can9);shuchu5=can10/PI*180;return shuchu5;/化弧度为角度float JD(int year,int month,int day)int A,B,days;A=floor(year/100);B=2-A+floor(A/4);if(month2)days=floor(365.25*(year+4716)+floor(30.6001*(month+1)+day+B-1524.5;elsedays=floor(365.25*(year+4715)+floor(30.6001*(month+13)+day+B-1524.5;return days;/儒略日double simplify(double jiaodu)double jiaodu2;if(jiaodu0)jiaodu2=jiaodu-floor(jiaodu/360)*360;else if(jiaodu0)jiaodu2=jiaodu+floor(fabs(jiaodu/360)+1)*360;else jiaodu2=0;return jiaodu2;/角度化简int getyear(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return year;/儒略日化年份int getmonth(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return month;/儒略日化月份int getday(int days)int a,b,c,d,e,z,k,year,month,day,f;z=floor(days+0.5);f=days-z;if(z2299161)a=z;elsek=floor(z-1867216.25)/36524.25;a=z+1+k-floor(k/4);b=a+1524;c=floor(b-122.1)/365.25);d=floor(365.25*c);e=floor(b-d)/30.6001);day=b-d-floor(30.60001*e)+f;if(e2)year=c-4716;else year=c-4715;return day;/儒略日化日期float jde(int Y,int M,int D,int hour,int min,int sec) int f,g; double mid1,mid2,J,JDE,A; if(M=3) f=Y; g=M; if(M=1|M=2) f=Y-1; g=M+12; ; mid1=floor(365.25*f); mid2=floor(30.6001*(g+1); A=2-floor(f/100)+floor(f/400); J=mid1+mid2+D+A+1720994.5; JDE=J+hour/24+min/1440+sec/86400; return JDE; /历书儒略日/ ; double moonangel(int X,int Y,int year,int month,int day,int hour,int min,int sec) / int i,c; float JDE,T,L1