2015年国赛高教杯数学建模
C题 月上柳梢头
月上柳梢头，人约黄昏后”是北宋学者欧阳修的名句，写的是与佳人相约的情景。请用天文学的观点赏析该名句，并进行如下的讨论：
1. 定义“月上柳梢头”时月亮在空中的角度和什么时间称为“黄昏后”。根据天文学的基本知识，在适当简化的基础上，建立数学模型，分别确定“月上柳梢头”和“人约黄昏后”发生的日期与时间。并根据已有的天文资料（如太阳和月亮在天空中的位置、日出日没时刻、月出月没时刻）验证所建模型的合理性。
2. 根据所建立的模型，分析2016年北京地区“月上柳梢头，人约黄昏后”发生的日期与时间。根据模型判断2016年在哈尔滨、上海、广州、昆明、成都、乌鲁木齐是否能发生这一情景？如果能，请给出相应的日期与时间；如果不能，请给出原因。
整体求解过程概述(摘要)
本文在认定元宵节“月上柳梢头”，“人约黄昏后”同时发生的基础上，首先建立“昏影终”模型得到元宵夜“黄昏后”的具体时间范围，然后建立“月梢头”模型，计算出元宵夜黄昏后时间所对应的“月上柳梢头”角度，在这些条件下借助模型运用MATLAB分别计算出2015年北京地区“月上柳梢头，人约黄昏后”的日期与时间，并用此结果推算2016年北京、上海等城市的日期与时间。
针对问题一，论文引用了晨昏蒙影的天文现象对黄昏时段进行定义，将太阳中心处于地平线到太阳中心处于地平线以下 -18’的时间范围定义为“黄昏后”。根据诗句情景，可知 “黄昏后”与“月上柳梢头”两个情景在元宵节同时出现，因此本文首先建立“昏影终”模型确定出2005—2015年元宵夜“黄昏后”所对应的具体时间段，然后建立“月梢头”模型确定该时间段对应的月亮在空中的角度，通过对数据的分析定义月亮高度角为 10’到 15’较为合理。最后利用昏影终模型和月梢头模型计算出2015年北京地区的“黄昏后”和“月上柳梢头”时间段（数据见表6、7）。
针对问题二，将2015年北京地区“月上柳梢头”和“人约黄昏后”的时间段进行比对，确定出2015年“月上柳梢头，人约黄昏后”发生的日期与时间，观察其分布规律，可以发现每个月的农历十三、十四、十五、十六、十七、十八这六天可能会同时出现该情景。使用问题一模型计算出2016年北京“月梢头”、“黄昏后”的时间，并对每月的上述六天进行比对，确定出2016年北京发生14次该现象（具体日期时间见表12），借助天文航海的时间计量方法以北京为基准推算哈尔滨、上海、广州等城市的日期和时间。
本文利用模型从定性、定量两个方面对“月上柳梢头”的角度和“黄昏后”的时间进行了定义，并且考虑了天体自转、公转等方面的因素。
模型假设：
假设一：太阳、月亮的高度和方位角与太阳、月亮在天球中的坐标和观测点所处的位置有关；
假设二：太阳、月球在天球中的坐标是时间的函数；
假设三：文中关于太阳、月亮的天文参数仅随时间变化；
假设四：元宵夜“月上柳梢头”与“人约黄昏后”的情景同时出现；
假设五：文中结果涉及的日期、地区天气良好利于赏月；
假设六：“昏影终”为太阳高度角等于 -18’的时刻.
问题分析：
针对问题一，题目要求分别定义“月上柳梢头”时月亮在空中的角度和“黄昏后”的时间日期与时间。由于诗句“月上柳梢头，人约黄昏后” 的背景是元宵夜，也就是说在元宵夜“月上柳梢头”和“人约黄昏后”这两个情景会同时出现，此刻的时间、角度就是问题需要的定义。因此本文首先建立“昏影终”模型确定元宵夜“黄昏后”所对应的时间段，然后建立“月梢头”模型确定该时间段对应的月亮在空中的角度，最后借助这两个模型计算出2015年“月上柳梢头”和 “人约黄昏后”分别出现的日期与时间。
针对问题二，问题一中，已经将2015年北京地区全年的“月上柳梢头，人约黄昏后”发生的日期和时间分别求出，并且验证了模型的正确性和可行性，同样根据模型将2015 年北京地区全年“月上柳梢头，人约黄昏后”发生的日期和时间分别求解出来，再对数据进行对比，进一步得出时间段的重合部分，即为“月上柳梢头，人约黄昏后”发生的日期与时间，观察其分布规律，并尝试将此规律应用于2016年，然后使用问题一模型计算出2016年相应日期的时间，进一步确认结果。最后利用天文学中的经度差公式推算出其余六个城市发生该景象的时间。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：
clc;
clear;
YY=2012;
MM=2;
mm=[31 28 31 30 31 30 31 31 30 31 30 31];
DD=6;
TT=10;
YY=YY+(mm(MM)+DD)/365;
Jt=365*(YY-1900)+(YY-1901)/4+0.5;
T=(Jt-2415020)/36525;
h=279.69668+36000.76892*T+0.0003*T^2;
e=(23+27/60+8.261/3600)+46.845/3600*T-0.0059/3600*T^2+0.00183/3600*T^3;
s=270.69668+481267.89057*T+0.00198*T^2+0.000002*T^3;
p=334.32956+4069.3403*T-0.01032*T^2+0.00001*T^3;
N=259.18328-1934.14201*T+0.00208*T^2+0.00002*T^3;
qm=s+0.109760*sind(s-p)+0.022236*sind(2-2*h+p)+0.011490*sind(2*(s-h))+0.003728*sind(2*s-2*p);
wm=0.089504*sind(s+N)+0.004897*sind(2*s-p+N)+0.004847*sind(p-N)+0.003024*sind(s-2*h+N);
xx=117;
yy=40;
cc=(TT-8)*15+h+xx-180;
CC=cosd(cc);
f=sind(yy)*(sind(e)*cosd(wm)*sind(qm)+cosd(e)*sind(wm))+cosd(yy)*(cosd(qm)*cosd(wm)*CC*(sind(qm)*cosd(wm)*cosd(e)-sind(e)*sind(wm)));
acosd(f)
clc;
clear;
YY=2015;
MM=2;
DD=12;
for TT=1:24
xx=117;
yy=40;
TG=TT-8;
T=fix(1461*(YY-1900)/4)+fix(153*MM-2)/5+DD+TG/24-36557.5;
TD=T/36525;
ET=TD+(3.17*TD+1.43)*10^(-8);
a=[0.00256 1934 235;0.00015 72002 201];
aa=0;
for i=1:2
aa=a(i,1).*cosd(a(i,2)*TD+a(i,3))+aa;
end
e=23.43928-0.01301*TD+aa;
b=load('1.txt');
ab=0;
for i=1:62
ab=b(i,1).*cosd(b(i,2)*ET+b(i,3))+ab;
end
ab1=b(63,1)*ET*cosd(b(63,2)*ET+b(63,3));
q=ab+ab1;
c=load('2.txt');
w=0;
for i=1:45
w=c(i,1).*cosd(c(i,2).*ET+c(i,3))+w;
end
aaa=atand(cosd(e)*cosd(w)*sind(q)-sind(e)*sind(w));
bbb=asind(cosd(e)*sind(w)+sind(e)*cosd(w)*sind(q));
a1=[0.0048 1934 235;-0.0004 72002 201];
d4=0;
for i=1:2
d4=a1(i,1).*sind(a1(i,2)*TD+a1(i,3))+d4;
end
sg=100.4606184+36000.77006*TD+0.00038793333*TD^2+d4+15*TG;
t=sg+xx-aaa;
hc(TT)=asind(sind(yy)*sind(bbb)+cosd(yy)*cosd(bbb)*cosd(t));
end
hc=[1:24;hc];
