赛区评阅编号（由赛区组委会填写）：
XXXX高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号（从A/B/C/D中选择一项填写）： 我们的报名参赛队号（12位数字全国统一编号）： 参赛学校（完整的学校全称，不含院系名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)：
日期： 年 月 日
（此承诺书打印签名后作为纸质论文的封面，注意电子版论文中不得出现此页。以上内容请仔细核对，特别是参赛队号，如填写错误，论文可能被取消评奖资格。）
赛区评阅编号（由赛区组委会填写）：
XXXX高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅记录（可供赛区评阅时使用）： 评
1
一 寸 光 阴 不 可 轻
阅 人 备 注
送全国评奖统一编号（由赛区组委会填写）：
全国评阅统一编号（由全国组委会填写）：
此编号专用页仅供赛区和全国评阅使用，参赛队打印后装订到纸质论文的第二页上。注意电子版论文中不得出现此页，即电子版论文的第一页为标题和摘要
页。
基于matlab与太阳方位角的经纬度计算方法
摘要
根据影子的变化挖掘出测量地点的信息是一项有挑战性的数学工作，这一工作可能会应用到安全领域的工作之中，本文利用影子的数据挖掘出太阳高度方位信息进而求解出所测量地点的经度纬度实现了成功定位。
针对问题一：我们已知该地点位于北京，并且以北京时间计时，通过分析时角，太阳高度角，以及当天太阳直射位点的关系，我们得到了影子长度与时间的复杂关系模型，为了精确绘制函数图像，我们在这里采用了根据曲率的变化自适应采样绘图的技术，得到了较为精确的函数图像，通过分析，基本符合实际情况。
针对问题二：我们利用已知数据，挖掘出了更多有效信息，通过对影子长度以及时间累积量进行二次多项式拟合，我们找到了包括正午时间。利用正午时间与北京正午时间的差距，我们找到了当地所在的纬度。接下来我们针对x，y坐标进行散点绘图，发现它们分别呈现线性增长的特性，在这里我们利用最小二乘法找到了其中的线性关系。利用上一步求解出的正午时间，我们求解出了正午影子朝向，即正北方向。在问题一建立的数学关系模型上，我们又利用matlab求解出了相对精确的纬度信息，信息显示，这一地点大致位于我国乌鲁木齐附近。
针对问题三：大致沿用了问题二的数学模型，我们确定了几个可能的日期，求解出了三个可能的坐标：东经107.5°，北纬44.7°，拍摄日期9月30日；
2
一 寸 光 阴 不 可 轻
东经107.5°，北纬14.79°，拍摄日期11月1日；东经107.5°，北纬20.59°，拍摄日期12月1日。
针对问题四：由于需要从摄像机视频中先测量相关信息，这存在一定的误差。我们在这里一方面利用像素个数进行较为精确的计数测量，另一方面利用透视原理，对机位测量数据进行了一定的矫正，得到了较为精确的数据。继续沿用第二个，第三个模型得到了较为精确地解。其解为：拍摄时间6月23日，北纬50.5521°，东经101°，大致位于蒙古境内；拍摄时间为7月23日，北纬41.8135°，东经：101°， 大致位于内蒙古阿拉善盟；拍摄时间为8月23日， 北纬33.1815°，东经：101°， 大致位于青海省果洛藏族自治州班玛县。 关键字：最小二乘法 自适应绘图 matlab 机位矫正 数值求解
问题重述：
如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。
1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出XXXX年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。
2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。
3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。
4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？
模型的假设：
1. 太阳直射点在南北纬回归线的运动大致视为匀速运动
3
一 寸 光 阴 不 可 轻
2. 影子长度仅受太阳高度角的影响，切周围没有人工光源，玻璃幕墙的影
响。
3. 题中给出的数据是经过精确测量的。 4. 影子投影地面是光滑的，没有倾斜。
符号的说明：
A 太阳方位角 h 太阳的高度角
φ 某地的纬度 δ 太阳直射地点的纬度 t 当地在某时刻的时角 b 影子长度
前期准备：
1.某地的正午太阳高度角： H（当地）=90°-纬度差（*同一纬度
相减，异纬相加）；
2．太阳高度角随着地方时和太阳的赤纬的变化而变化。太阳赤纬（与太阳直射点纬度相等）以δ表示，观测地地理纬度用φ表示（太阳赤纬与地理纬度都是北纬为正，南纬为负），地方时(时角)以t表示，有太阳高度角的计算公式： sin h＝sin φ sin δ＋cos φ cos δ cos t
3.经过查资料，九月23日为秋分日，太阳直射赤道。太阳直射点从赤道南移的过程可大致简略为匀速运动。这一天是十月22日，这一天太阳直射纬度是某日（R）太阳直射点的地理纬度位置=0°+（R—9月23日）*（23°26′*4/365），即为0°+30*（23°26′*4/365），即为462.2460′，换算成度，即δ= 7.7041°。
4.昼夜长度的确定，某地的日出日落时间，需要根据当日太阳直射点纬度来确定。如图，晨昏线与所求解纬度的交点即为日出时间点和日落时间点，两者相隔的度数φ2*24/360即为此日的昼长。经过资料查询，日出时间T=12*arcCos（tgA*tgB）/π，B为当地纬度，A为太阳直射纬度。日出时间和
日落时间是以当地的中午12点为对称的。如果是用其他时区时间计时， 则需要计算半个白昼时间再在正午时间的基础上加减来计算日出日落时间。
4
一 寸 光 阴 不 可 轻
5．时角的确定。某地时角在当地正午为0，由于地球每小时大致自转15°，所以上午为(正午时间-当下时刻)*15°，下午为(-正午时间+当下时刻)*15°。
6通过资料查询，我们得到，一天之中的太阳方位角是一个变化的，其具体的变化方程为cosA = (sinhsinφ - sinδ) / (cosh cosφ)。
模型一：
1.1建立模型：
根据资料 sin h＝sin φ sin δ＋cos φ cos δ cos t，h为太阳高度角，φ为该地纬度，δ为这一天的太阳直射点纬度，t为时角；
又因为影子长度仅仅与物体有效高度以及此地太阳高度角有关，所以影子长度l，l=x*tan（h）；其中，h为太阳高度角，x为物体有效高度，l为影子长度。
这一天是十月22日，这一天太阳直射纬度是某日（R）太阳直射点的地理纬度位置=0°+（R—9月23日）*（23°26′*4/365），即为0°+30*（23°26′*4/365），即为462.2460′，换算成度，即δ= 7.7041°。
对于此地太阳时角t，地球自转一周为一天，即为24小时，不同的时间有不同的时角，以t表示。由于北京市以北京时间为参考计时，故北京的正午时间为北京时间12点整。地球自转一周为360°，因此每小时的时角为15°，即太阳时角表示为：t=15*(t2-12)，t2表示日照时数 。
所以，建立如下模型
（1）。sin h＝sinφsinδ＋cosφcosδcos t
5
一 寸 光 阴 不 可 轻
δ= 7.7041° =0.0214pi φ=39°54′2′′ （2）。l=x*tan（h），x=3 （3）。t=15*(t2-12)，9= 在图中，我们可以清楚地看到，一天中，影子长度是以中午为中心对称变化的，而影子长度在当地正午达到最 6 一 寸 光 阴 不 可 轻 短，这是与实际相符合的。而影子长度与时间也大致呈二次函数关系。 模型二： 2.1． 模型分析：由于原题中并没有指出如何建立的x-y坐标系，所以具体的南北正方向需要根据数据来进一步求解。通过分析我们可以知道，一天中正午时间太阳高度角达到最大，此时对于北半球的此地，太阳位于其正南方，影长达到最短，影子方向指向正北方方向。而影子方向始终背离太阳方向，也就是说我们可以根据影子方向确定每时每刻的太阳方位角。对于影子长度与时间的我们运用最小二乘法计算出它们的二次关系，进而推求出正午影长以及正午时间。对于短时间内之中x，y数据与时间的变化关系，我们通过对数据相关性的分析，可以知道x，y坐标在短时间内与累计时间呈现高度线性关系，我们在此应用最小二乘法进行线性拟合，求解出其线性函数，经过检验。而对于这一问题，这一天在春分日前后，昼长大致12小时，太阳直射赤道处。 1）。对于经度的求解，我们可以根据已经计算出来的当地正午北京时间，以北京所在的东八区为参照，求解出所在经度。 2）。对于所求的纬度，我们建立方位角与时角，太阳直射地点的关系模型，利用处理好的数据，运用二分法进行精确求解。 7 一 寸 光 阴 不 可 轻 2.2 模型的建立： 2.2.1 经度模型的建立 我们将数据导入matlab，通过计算影子长度，我们得到从14:42到15:42的影子长度相应的数据如下： 14:42 14:45 14:48 14:51 14:54 14:57 15:00 15:03 15:06 15:09 15:12 15:15 15:18 15:21 15:24 15:27 15:30 1.149626 1.182199 1.215297 1.249051 1.283195 1.317993 1.353364 1.389387 1.426153 1.4634 1.501482 1.540232 1.579853 1.6XXXX 5 1.661271 1.703291 1.746206 8 一 寸 光 阴 不 可 轻 15:33 1.790051 15:36 1.835014 15:39 1.880875 15:42 1.927918 由于这里的时间格式无法当做自变量数据，故我们在这里把时间转化为从十二点到此时此刻的分钟累计。我们以时间累积为自变量，以影子长度为因变量，用matlab拟合出它的散点图： 我们观察这一段散点图，不难发现，两者大致呈二次多项式关系。我们接下来运用matlab中的多项式拟合函数polyfit针对这两个变量进行多项式拟合，找到的时间与影子 9 一 寸 光 阴 不 可 轻 长度关系式如下： y=0.000293612207924286*x^2-0.0733391374619733*x+5.5110。 matlab中在函数绘图的时候一般情况下是先对自变量和函数值进行等间距采样，这样就会导致某些函数在曲率较大的地方函数图像信息的不准确，故我们采取一种新的函数拟合方法，即根据函数曲率随时调整采样点密度，这样可以精确绘制出函数图像，图像如下： 根据二次函数性质，二次函数在-2a/b处取得全局函数极小值。接下来我们可以找到该函数的最低点坐标：（124.8912，0.9313），也就是说影子长度在北京时间下午2点4分达到最小，此时为该地正午时分。 由于北京的坐标位于北纬39度54分26秒,东经116度23分29秒，位于东八区。而地球上每一经度对应时差为四分钟，该地按时间来说早于北京，所以，我们可以推算出该 10 一 寸 光 阴 不 可 轻 地的经度为85度，结合我国地图，如下图，该经度附近大概通过我国新疆和西藏，尼泊尔，印度。 2.3 纬度求解模型的建立 前期，我们已经分析过太阳在某地一天中高度角h变化的具体公式，sin h＝sin φ sin δ＋cos φ cos δ cos t，这一天在春分日前后，昼长大致12小时，太阳直射赤道处，所以δ=0，所以公式可以简化为sin h＝cos φ cos t，其中，t为时角，φ为当地纬度。 对于任意时刻，影子长度l与物体高度x，太阳高度角h的关系式如下l=x/tan（h）。 接下来 根据前期分析中的一天中太阳方位角变化公式：cosA = (sinhsinφ - sinδ) / (cosh cosφ)，由于δ=0，所以可以化简为cosA = (sinhsinφ) / (cosh cosφ)。 这里对于三个方程，我们有四个未知量，故我们还要找到一组模型关系。 对于太阳方位角，我们可以进一步求解。在这里我们对第二 11 一 寸 光 阴 不 可 轻 问附件1当中的x，y轴坐标进行绘图分析。得到如下散点图。 我们观察到，在相对较短时间内，随着时间的推移x轴坐标大致呈线性增长。故我们对x轴坐标和累计时间利用最小二乘法进行线性拟合，所谓的最小二乘法（generalized least squares）就是是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。 最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小。 最小二乘法通常用于。 其具体做法就是在给定的函数类中，对于给定数据 ， i=(1,2,3…….m),求出函数p(x),使误差的平方和达到最小，即： 我们利用matlab进行编程，得到以下关系模型： A=0.0131*t-1.1112 其中，A为太阳方位角，t为从北京时间十二点的累计时间。我们求得当地正午时候影子顶端的x坐标为0.5249。同理， 12 一 寸 光 阴 不 可 轻 我们针对y坐标进行绘图分析，得到如下散点图： 我们观察到，随着时间的推移，y坐标的改变呈现高度的线性增长特性。我们继续利用最小二乘法对其进行线性回归，得到以下函数： Y=0.0019*t+0.1835 Y为影子端点的纵坐标，t为时间的累积量。我们计算求得这个地方正午时影子端点的纵坐标是0.4191。 所以我们就求得了这个地方正午时分影子端点的坐标（0.5249，0.4191），这个坐标就是正北方向坐标。利用正北方向，我们可以利用余弦定理，求得这21个采样时间点的太阳方位角余弦。余弦定理是：三角形中,任意一边的平方等于另外两边的平方和减去另两边及其夹角的余弦的积的两倍.。我们求解得到这21组。 建立模型方程组如下： 1. cosA = (sinh*sinφ) / (cosh *cosφ) 13 一 寸 光 阴 不 可 轻 2. l=x/tan（h） 3. sin h＝cos φ* cos t 我们利用matlab的fzero函数进行求解，求解得到的纬度值为北纬42.7310°，所以该地的坐标是：北纬42.7310°，东经85.23°，大致在乌鲁木齐附近。 模型三： 3.1 模型的建立： 对于问题三，我们依旧有如下的方程组： sin h＝sin φ sin δ＋cos φ cos δ cos t b=x*tan（h） cosA = (sinhsinφ - sinδ) / (cosh cosφ) 分析这个方程组我们可以知道，对于问题三，我们依然可以利用模型二求出太阳方位角的余弦值cosA，影子长度b。而我们又知道，对于任意一个确定的日期，该日的太阳直射点纬度δ也是一个确定的数，因此问题求解的关键在于日期的确定。 为了简化模型，考虑到日期的任意性，我们取了3组在9月22号（秋分）之后的日期：9月30日，11月1日，12月1日，经过查资料，九月23日为秋分日，太阳直射赤道。太阳直射点从赤道南移的过程可大致简略为匀速运动。太阳直射纬度是某日（R）太阳直射点的地理纬度位置δ=0°+（R—9月23日）*（23°26′*4/365），利用matlab计算可知 14 一 寸 光 阴 不 可 轻 这三天的太阳直射点纬度分别是南纬1.8027°，9.5288°，17.2548°。 3.2 模型的求解 3.2.1经度的确定 我们拟合出了影长b关于累计时间t的关系式：[0.000702512197687887]*t^2+[-0.0707693625707730]*t+14.0424，同时利用自适应拟合得到了下图所示关系图： 由此利用二次函数的求极限值公式求出了图像的极小值点为50，影长b最小为12.2601米。由此可知当地的正午时分为北京时间12点50分，由于地球上每一经度对应时差为四分钟，而该地正午时分晚于处于东八区的北京，因此可以推算出该点的经度是东经107.5°。结合我国地图，该地在我国可能的省份有内蒙古，陕西，贵州以及广西。 3.2.2纬度的确定 （1）利用最小二乘法以及线性回归，我们确定了影子端点的坐标值与累计时间的关系，分别是： X=0.0208873*t - 0.296692 15 一 寸 光 阴 不 可 轻 Y=3.38826 - 0.000919091*t 根据模型二的分析，可以求得正午时刻，在问题三给出的坐标系下，正北方向的坐标为（ 0.7554，3.3420）。以原点，正北方向坐标点，已经影子端点构成三角形，在这个三角形中利用余弦定理可以求得太阳方位角cosA。带入21组已知的影子端点坐标我们可得到每个时刻对应的cosA。21组数据拟合成三角函数图像如下： 接下来，我们需要确定时角。根据求出的该日的正午时间由时角的定义，可以确定时角 t=15*(t2-12.83)，13.15= 16 一 寸 光 阴 不 可 轻 这里继续之前建立的关系模型： sin h＝sin φ sin δ＋cos φ cos δ cos t b=x*tan（h） cosA = (sinhsinφ - sinδ) / (cosh cosφ) t=15*(t2-12.83） 利用lingo进行数值求解。 综上所述，可能的几组拍摄日期和拍摄地点是： 1.9月30日 地理坐标东经107.5°，北纬44.7°，大致位置在蒙古国境内 2.11月1日 地理坐标东经107.5°，北纬14.79°，大致位置在越南，老挝，柬埔寨三国交界处 3.12月1日 地理坐标为东经107.5°，北纬20.59°，大致位置在越南下龙湾观澜岛 模型四： 4.1 前期分析：根据视频采集挖掘有效信息，估算影子长度，影子方向使这一问题中的不同于以上几个模型的有趣的地方。我们注意到，随着时间的推移，影子长度方向都在变化。我们还可以观察到，摄像机的视角并不是垂直朝向杆子的，而是向左偏斜了一定角度，根据透视观点的知识，在这样的情形下，我们观察到的影子是短于实际长度的。故我们需要对机位进行人工矫正。 17 一 寸 光 阴 不 可 轻 视频中时间长度为四十分钟，经过粗略估计，影子转过9°。我们在视频当中每四分钟采集一帧图像，根据比例关系量取并且计算屏幕中影子长度b1.根据最后一幅图像中的影子方向建立x轴坐标，进而建立y轴坐标。 根据中间地板线修正视角，修正角度a0=90-（地边中间的缝线与图片下沿的夹角）。影子长度为屏幕测算长度b=b1/sin(a0)。b即为实际长度。由于这四十分钟里影子大概转过9度，也就是每一次采样，转过0.9度，根据刚才建立的坐标系求解这十组影子端点坐标a 地板中间缝线就是右图： 我们通过测量得到，a0=0.8921（弧度制），经过校正后的影子长度：l2/cos(a)。 在量取影子长度的时候，由于整个过程中影子大致只是转过9度，所以我们就认为影子的旋转不会对测量造成影响，而无论是图中还是现实生活中，竿子和影子的比例都是一样的。所以，我们可以根据图中影子和杆子的长度比例关系测定并计算得到未矫正前的影子长度，影子长度如下： 经过校正后的影子长度：l2/cos(a)，这是经过矫正后的影子的长度变化图： 18 一 寸 光 阴 不 可 轻 经过校正后的影子长度如下： 在这里，我们依旧用二次多项式拟合影子长度和累计时间的关系，通过二次函数的性质，我们直接求得这个二次函数关系模型的最低点，即从0:00累计分钟时间为757.275060307610时，即中午北京时间13:04，影子长度0.6667米。 不难判断，这一点的经度为东经101度附近，大致通过我国内蒙古，甘肃，四川，云南，泰国等地。 接下来我们计算影子顶点的坐标，由于我们沿最后一张照片中的影子方向建立坐标系，而且这四十分钟里面，影子转过9度，也就是说每次采样影子都会顺时针转过1度。我们得到这几组影子端点坐标： 19 一 寸 光 阴 不 可 轻 接下来，如模型二，我们对想x，y坐标分别对时间进行先行拟合，得到以下线性方程关系模型：x=-0.0209735146747096*t+14.8514885887912 Y=-0.00797222624247121*t+4.56472399190384 我们进而可以求得正午影子顶点坐标是 (-1.5917，-1.6855)，即 (-1.5917，-1.6855)为正北方向，以及各个采样时间点的太阳方位角余弦： 0.8XXXX5874420970 0.815274000753605 0.810350671835671 0.805372038212327 0.800349618723757 0.795237983243519 0.790070169457562 0.784882768197391 20 一 寸 光 阴 不 可 轻 0.779613525981249 0.774271891675362 0.768892852513878 由以上数据可求出每一时刻的太阳高度角h，我们在上面还计算得到正午时间，所以我们还可以计算出各个时间对应的时角，我们已知太阳高度角与纬度，时角，太阳直射纬度的关系模型： sin h＝sin φ sin δ＋cos φ cos δ cos t 用采集得到的11组数据分别带入这个关系模型，求解即可。 我们得到该地点的经纬度：北纬47.2°，东经101°。 4.2问题2 在这问题之中，日期再次成为了未知，所以我们根据图片里的景色，判断当地处于夏天，所以我们找了三组日期，以一个月为步长，确定一下日期：6月23,7月23,8月23，太阳分别直射南纬 23.26，15.612876712329, 23.26°， 7.9657534246575，即δ分别等于15.612876712329°, 7.9657534246575°。 我们利用以下关系模型：这里我们用w表示纬度。 sin(h)=sin(w)*sin(δ)+cos(w)*cos(δ)*cos(t) cosA = (sinhsinφ - sinδ) / (cosh cosφ) 最终利用matlab数值求解方法解得： 6月23： 北纬50.5521° 东经：101° 大致位于蒙 21 一 寸 光 阴 不 可 轻 古国境内 7月23： 北纬41.8135° 东经：101° 大致位于内蒙古阿拉善盟 7月23： 北纬33.1815° 东经：101° 大致位于青海省果洛藏族自治州班玛县 五．总评与分析 5.1 模型的评价： 模型一，模型二，模型三中，我们分析了太阳高度角，时角，方位角，太阳直射位点的数学关系，并且建立了关系模型，在对影子端点的坐标处理的时候，利用正午影子长度最短的原理，我们找到了三个模型的正午时间，以及正北方向，这对接下来的求解产生了极大地便利。 在模型四中，我们在视频中截取了按照等间距取样原理数帧图像，手动测量的影子长度，在这里我们为了尽可能减小误差，我们采取了计算像素个数的方法。由于摄像机并不是朝向正前方，所以我们利用透视原理对机位数据进行了数据矫正，尽可能的减小了误差。 5.2模型的改进 在根据日期计算太阳直射纬度时，由于将太阳的变化简化成了匀速运动，这里可能会造成一定的误差，因为实际上太阳直射点在南北回归线之间的运动是类似三角函数的。 在进行方程组求解的过程中，由于利用插值法不断缩小 22 一 寸 光 阴 不 可 轻 精度进行变量求解，这会在一定程度上造成解得误差。故寻找更有效的函数求解方法是非常必要的。 参考文献 [1] 汪晓银，周保平，数学建模与数学实验（第二版），北京：科学出版式，XXXX.8 [2] 汪晓银，邹庭荣，周保平，数学软件与数学实验（第二版），北京：科学出版社，XXXX.8 [3] 姜启源，叶其孝，数学建模，北京：机械工业出版社，XXXX.8。 [4] 同济大学数学系，高等数学（第二版）上册，上海：同济大学出版社，XXXX.10 [5]， 解析太阳高度角 ， [6] 薛定宇，陈阳泉，高等数学问题的MATLAB求解，北京：清华大学出版社,XXXX 23 因篇幅问题不能全部显示，请点此查看更多更全内容