数学建模标准优秀论文模板样本.doc

资料内容仅供您学习参考，如有不当或者侵权，请联系改正或者删除。
太阳影子定位
摘 要
　本文针对如何根据直杆在水平地面上的太阳影子顶点坐标, 确定其可能的地点或者日期的问题, 首先对不同时刻影子顶点坐标用Matlab软件进行分析, 然后利用太阳高度和方位角与太阳时和经纬度的天文学关系, 建立不同数学模型求解。
针对问题一, 根据常规的地理学中关于太阳高度和方位角的相关计算公式, 并引入了更加精确的年度订正、 时差订正, 建立了利用太阳方位角求解规定日期内定长直杆不同时刻影长模型, 经过利用Matlab软件进行计算, 生成了定长杆影长随时间变化的函数关系式, 并画出了该函数所对应的曲线图象, 图象表明所建模型符合与实际相符。
针对问题二, 首先规定了大地坐标系的正方向, 根据测量地点所处地理位置和当地时间的不确定性, 可将附件中每组数据( ) 分为8种可能的坐标情况。然后, 根据直杆影子坐标与方位角的关系, 又基于问题一的模型, 得出一个包含, 纬度和经度的关系式, 并利用excel软件和Matlab软件, 代入和拟合了影长随时间变化的曲线, 利用时差公式求解出经度, 然后分别代入8种可能的坐标, 最终确定了其中2种可能的当地经纬度坐标, 为( ) 即南海附近和( ) 即俄罗斯中部。
针对问题三, 其在与问题二相同数据量和条件的前提下, 多了一个未知量积日, 基于问题二中的模型, 在保证结果准确的前提, 为了减少未知量的个数使模型可解, 舍弃了年度订正。然后得出一个关于太阳方位角二元方程, 由于纬度和积日虽然为两个变量, 但在此问题中必须对应, 因此常规的穷举法和约束算法解决此类问题较为繁琐, 本模型引入二维交集点概念, 较为简洁的刻画了和的对应关系, 在Matlab软件中综合运用穷举法, 循环算法和二维交集等算法, 解得附件2对应的当地坐标为( ) 即巴基斯坦境内, 日期为7月21日, 解得附件3中当地坐标为( ) 即陕西省内, 日期为1月15日
。
针对问题四, 我们首先利用Matlab软件导入并处理了视频, 均匀截取了19帧图像, 并将其导入Digimizer 软件进行图像处理, 选定杆长为参照, 进行灰度处理, 自动测量影长长度数据, 采集到了满足预期要求的直杆影长与对应时刻( 北京时间) 数据。依据这些数据, 在已知视频拍摄日期的前提下, 利用问题一模型, 在已知影长数据反向求解测量地点坐标, 解得测量地点坐标为( ) 即印度尼西亚附近, ( ) 即新疆地区。在拟解决问题中, 假设未知视频拍摄日期, 我们利用问题三模型的主体框架, 首先得出一个关于太阳高度角 的二元方程, 再利用模型二中的算法, 在Matlab中计算得出拍摄日期为 7月15日, 地点为印度尼西亚和新疆地区, 与第四问中第一小问视频拍摄日期和地点基本吻合。

关键词: 影子定位 太阳高度角 太阳方位角 Matlab软件 二维交集点

一、 问题重述
为了确定视频拍摄的时间和点, 我们一般采用太阳影子定位技术, 也就是利用视频中物体的影子随时间的变化而计算出视频拍摄地点和时间, 这是视频数据分析的重要内容之一。
首先, 在已知当天某时间段内直杆影子顶点坐标变化, 杆长, 日期和地点的前提下, 如何确定影子长度随时间变化的规律。
然后, 在已知日期和某时间段内直杆影子顶点坐标变化的前提下, 如何确定测量地点经纬度。
如果测量日期未知, 如何再经过某时间段内直杆影子顶点坐标变化来确定日期和经纬度。
最后, 在视频拍摄日期已知的前提下, 如何根据一段视频中直杆( 杆长已知) 影子的变化来确定拍摄地点, 若拍摄日期未知, 又该如何确定拍摄日期和地点。

二、 问题分析
4.1 问题一分析
问题一给出了日期时间地点杆长所构成的数据, 让其求解影长随时间的变化曲线, 日期的给出可求得太阳的赤纬角, 全年中每天的赤纬角都在发生变化, 每一天对应一个赤纬角; 根据所给的时间段, 综合观测地时间与北京时间可求得时角, 先建立时角与影长的关系, 继而再求时间与影长的关系, 问题一即迎刃而解。
4.2问题二分析
根据问题二中给出若干个可能的地点, 可知该问至少有两个以上不同的解, 根据附件1所提供的数据, 由问题一的结果可反向推出当地影长最短时( 即当地正午时) 所处的北京时间, 由二者的时间差, 可计算出观测地与”北京时间”所在的经度圆的精度差, 继而求得观测地经度; 求其观测地纬度时引入新变量太阳方位角, 然后查阅相关文献资料获得太阳方位角, 太阳高度角, 赤纬角三者的关系, 使求解观测地纬度问题迎刃而解。问题二建立的是投影物体和阴影轨迹间的几何关系模型。
4.3问题三分析
问题三是由问题二进一步深入所得, 问题二时间已知, 即可求得太阳的赤纬角, 问题三没有给出时间, 即在问题二解题过程中新增加了一个未知量即太阳赤纬角; 可先根据附件2、 3求出各观测地所处的经度( 如上问题2中经度的求解分析) , 进而求得时角; 因为赤纬角的大小可由年积日导出, 而年积日又是存在一个范围的( 即从1至366) , 因此问题三的纬度解决方法可采取穷举法计算, 问题二中的纬度计算只是问题三中纬度计算的一部分。综上根据附件2和附件3的固定直杆在水平地面上的太阳影子顶点坐标数据, 问题三建立数学模型同样是投影物体和阴影轨迹间几何关系的数学模型, 并引入新的解题方法穷举法, 综上即可得出若干个可能的地点与日期。
4.4问题四分析
问题四的附件为一段视频, 视频中为一已知高度直杆随时间的变化引起影长的变化, 将附件四视频导入Matlab, 均匀截取19帧图像后, 导入Digimizer采用灰度处理方法, 确定出影长与时间的关系曲线, 即可求出每个时间所对应的影长( 但太阳方位角未知) ; 影长确定后又回归到问题2的求解分析, 即已知时间和影长拟合曲线求出观测地所在经度, 并求出观测地正午时刻的影长, 然后根据直杆高度, 影长, 时角( 由时间导出) , 赤纬角( 由日期) , 即可求得观测地的纬度; 问题四的另一部分在没有拍摄日期的情况下求解观测地的位置, 拍摄日期未知即新增了一个未知量赤纬角, 在将视频处理完成后, 即可采用问题3中解题方法, 采用穷举法求解最合适的日期, 进而求得观测地的地理位置。

三、 符号说明
符号
含义
单位
备注

影子顶点横坐标
米
沿正北方向为正方向

影子顶点纵坐标
米
沿正东方向为正方向

太阳方位角
度
太阳光线在地平面上的投影与当地经线的夹角, 本文近似指竖立在地面上的直线在阳光下的阴影与正南方的夹角

年份
年
公元纪年法年份

积日
—
日期在一年内的顺序号

年度订正参数
—
用于订正因每年因天数取整而忽略0.2422天给每年带来的误差

年度订正值
天
订正因每年因天数取整而忽略0.2422天给每年带来的误差所产生的量

日角
度
天文公式中的一个与太阳赤纬角等相关的变量

太阳赤纬角
度
地球赤道平面与太阳和地球中心的连线之间的夹角

平太阳时
时
平太阳时简称”平时”, 也就是我们日常生活中所使用的时间

太阳时差
时
这里指真太阳时和平太阳时的时差

时差
时
两个地区地方时之间的差别称作为时差

纬度
度
地球上某点与地球球心的连线和地球赤道面所成的线面角

经度
度
地球上某点点离本初子午线的南北方向走线以东或以西的度数

地方时
时
以一个地方太阳升到最高点的时间为正午12时, 将连续两个正午12时之间等分为24个小时, 所成的时间系统

真太阳时
时
1真太阳日分为24真太阳时的时间系统

太阳时角
度
观测点天球子午圈沿天赤道量至太阳所在时圈的角距离

太阳高度
度
太阳光的入射方向和地平面之间的夹角

四、 模型假设
( 1) 忽略因测量当地地形或者其它因素造成的直杆与水平面的夹角;
( 2) 忽略测量当地的海拔影响;
( 3) 在问题三模型中忽略平太阳时和真太阳时的时差;
( 4) 影长不受大气干扰所产生的阳光折射的影响;
( 5) 在某一天内不考虑地球的公转;
五、 模型建立与求解
5.1 问题一模型建立与求解
对于问题一, 直杆影长与原长的有如下关系:

=coth
( 1-1)
在直杆原长已知的前提下, 求得太阳高度角h即可求得影长, 从而进一步利用插值法模拟影长随时间变化曲线。
在给定具体日期、 位置坐标和特定时间段的前提下, 我们能够根据常规地理学对于太阳高度角的计算公式建立数学模型。
太阳高度角的计算公式为:

( 1-2)
联立式( 1-1) 与式( 1-2) 可得到一个能够直观分析出影子长度关于各个参数的变化规律的模型:

( 1-3)

下面分别对太阳赤纬角, 时角两个函数予以展开。
5.1.1 太阳赤纬角
图中日地中心的连线与赤道面间的夹角每天( 实际上是每一瞬间) 均处在变化之中, 这个角度称为太阳赤纬角。它在春分和秋分时刻等于零, 而在夏至和冬至时刻有极值, 分别为正负23.442°。

图1-1 地球绕太阳运行轨迹
太阳赤纬角有两个计算公式,

( 1-4)
和

( 1-5)
其中式( 1-4) 为常规计算太阳赤纬角的方法, 只引入了积日, 计算较为简便; 式( 1-5) 在式( 1-4) 的基础上又引入了订正量。
这里, 我们选用了式( 1-5) , 因为其进行了年度订正, 结果更加准确。
其中, 日角的计算公式为:

( 1-6)
年度订正为:

( 1-7)
年度订正参数为:

( 1-8)
在第一问中, 年份给定为 , 因此可求出。
5.1.2 太阳时角
太阳时角在常规算法中能够由地方时直接算出, 式为:

( 1-9)
这种算法计算比较简单, 可是不够准确, 式中应代入真太阳时, 此式未考虑真太阳时与地方平时的时差, 在本模型中引入了时差订正, 结果更为准确。
太阳时角在的正确计算公式为:

( 1-10)
其中真太阳时计算公式为:

( 1-11)
其中为地方时, 为太阳时差。
地方时与经度相关, 计算公式为:

( 1-12)
由于经度为已知, 因此可求出。
由于赤纬角和太阳时角都已求出, 因此将其代入式( 1-2) 可求出, 再将代入式( 1-1) , 可得出, 将以上模型编译后利用Matlab, 能够拟合出以下曲线:

5.2 问题二模型建立与求解
5.2.1 建立问题二模型前提声明
正东
图 1-2 规定坐标轴
正北
对于问题二, 由于附件中规定说明: 坐标系以直杆底端为原点, 水平地面为平面, 可是并未规定的正方向, 因此本模型规定的正方向如下图:

由于未知测量地点的经纬度坐标, 因此不能确定附件中测量地点是在南北哪个半球, 又不能由附件中的北京时间推算地方时间, 因此坐标数据并不知道是基于哪种正方向测量, 因此出现了多种数据可能, 这里把附件中的数据理解为影长在某两个垂直方向的投影, 认为是( ) , 把规定坐标轴上投影为( ) , 根据排列组合, 将一组( ) 用规定坐标系表示时, 会出现8组( ) , 即( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , 其中每一个坐标都代表一个不同的时间和地点组合。

5.2.2 建立问题二模型并求解
对于问题二, 我们利用问题一模型的基本公式框架, 结合用插值法拟合曲线方法, 建立了在已知测量日期、 测量时刻( 时区时间) 、 被测直杆影子( ) 的前提下, 求解测量地点经纬度坐标的数学模型。
首先, 我们按照问题一模型的基本公式, 建立问题二模型的框架。
已知有太阳方位角为:

( 2-1)
又有地理学中太阳方向角计算公式:

( 2-2)
由于附件中已经提供某时刻( 北京时间) 直杆影长( ) , 因此式( 2-1) 为已知, 因此能够联立式( 2-2) 。式( 2-2) 中包含一个函数和两个目标量和。根据问题一模型中的式( 1-2) 到式( 1-10) 能够联立得:

( 2-3)
联立式( 2-1) 与式( 2-3) 可得:

( 2-4)
易知式( 2-4) 为二元方程, 不能解出具体的经度和纬度, 因此利用真太阳时与时区时间的时差公式

( 2-5)
来求解经度。因为某地在正午太阳高度的时刻为地方时的12点, 式( 2-5) 中可由

( 2-6)
计算, 当地时间正午12点的判断特征为直杆影子最短, 因此由附件1中数据拟合的二次曲线中能够利用公式

( 2-7)
求得, 利用Matlab软件拟合的二次曲线如下:

得到=12.598, 代入式( 2-5) 得到经度=, 代入式( 2-4) 后得:

( 2-8)
得到一元一次方程( 2-8) 后, 任意从附件1中选取一组时间和坐标()后, 按照5.2.1的声明, 分别代入8组变换坐标, 在Matlab软件中计算后得出两个能够解出值的( ) 为( ) , ( ) , 对应=, ( ) , ( ) , 对应=。
因此问题二模型的解为( ) 即南海附近, ( ) 即俄罗斯中部。

5.3 问题三模型建立与求解
对于问题三, 由于其较问题二, 在相同的已知条件影子顶点坐标( ) , 北京时间S基础上, 求解更多的未知量经度, 纬度和积日。我们基于基于问题二模型的基本框架, 在Matlab软件中利用穷举算法, 循环算法, 寻求二维交集等算法, 建立了一个能够根据某时间段( 标准时间) 里直杆太阳影子顶点坐标数据, 确定直杆所处的地点和日期模型, 本模型分为两个部分: 推导二元方程和Matlab程序计算。
5.3.1 推导二元方程
根据问题一模型和问题二模型中的公式, 联立得出一个基本的求解模型:

( 3.1)
式( 3.1) 中可从影子坐标( ) 数据求得具体数值, 因此此式为一个三元方程, 根据问题二模型的插值法能够用Matlab拟合附件2和附件3中数据的二次曲线如下,

并分别根据式( 2-5) 到( 2-7) 求得影长最短时北京时间和当地经度。
附件2中=15.204, , 附件3中=12.735, 。
式( 3-1) 中代入值以后, 变为

( 3.2)
5.3.2 Matlab程序计算
因为中和相互关联, 必须一一对应, 因此在普通穷举法的基础上又必须要求解二维交集。
所谓二维交集, 是一种解决和必须对应的较为简易算法。我们把和作为一个二维点坐标值并标注下标( ) , 其中n的含义为第n个二维点, 在应用穷举法之后会产生若干二维点, 在分别用约束条件筛选时, 某个点的横坐标或纵坐标( 即) 不符合约束条件时, 整个点被排除。由此简化了对数据间复杂关联关系的简单刻画。
式( 3-2) 为一个二元方程, 依然无法用常规方法求解, 可是经分析易知如下约束关系:

且

( 3.2)

( 3.2)
这样能够利用Matlab软件利用穷举法, 计算出当时的, 构成集合与集合( 3-2) 求交集得到, 中元素与其求解其所用的一一对应, 即中元素的下角标与其对应的相等。任意一组影子顶点坐标( ) 和北京时间S数据能够得出一个对应集合。
本模型具体应用在附件2和附件3中时, 因为每一个附件都有21组影子顶点坐标( ) 和北京时间S数据, 因此会有21个, 表述方便起见, 将其标注为,
并归纳到大集合, 又因为每组( ) 在代入规定坐标系时会产生8组坐标数据, 因此事实上每个附件会产生8个, 分别对8个中所有元素求交集得出最终集合。中所有元素为最要求解的纬度值。同理, 的值也是如此筛选。同时又因为这些都带有下标, 也带有下标, 与具有复杂关系, 套用二维交集的方法, 在Matlab软件能够比较快捷得解出最终符合条件的二维点( ) 。
利用问题二模型, 解得附件2中的当地坐标为( ) 即巴基斯坦境内, 日期为····
附件3中的当地坐标为( ) 即陕西省内, 日期为····。
5.4 问题四模型建立与求解
对于问题四, 我们首先利用Matlab软件导入附件4的视频文件, 利用算法将视频均匀截取19帧图像, 再将图像导入Digimizer 软件进行图像处理, 选定杆长为参照, 进行灰度处理, 自动测量影长长度数据, 处理效果如下图
图1-3处理之后的图像

图1-4 处理之前的图像

之后将采集的19组时刻( 北京时间) 与影长的数据利用Matlab软件进行插值法拟合二次曲线

并分别根据式( 2-5) 到( 2-7) 求得影长最短时北京时间和当地经度。
视频中=13.8498, 。
然后根据是否已知拍摄时间, 即是否已知积日, 建立不同模型进行求解。
5.4.1 已知视频拍摄时间情况下模型建立与求解
在此种情况下, 能够套用问题一模型, 进行反向求解。联立式( 1-2) 和式( 1-4) 到式( 1-12) , 得到

( 4-1)
即

( 4-2)
由于为能够由采集的数据求得, 因此可求得。
我们进行了寻求多解分析, 经计算当天的太阳赤纬角, 即当天太阳直射纬度为, 因为采集视频的数据并未约束坐标问题, 因此分析得与纬度关于纬度的纬度也应该是一个解, 即

( 4-3)
因此。
因此, 在已知视频拍摄日期的情况, 解得当地坐标为( ) 即印度尼西亚附近, ( ) 即新疆地区。
5.4.2 未知视频拍摄时间情况下模型建立与求解
未知视频拍摄时间, 即未知积日, 与问题三模型非常相似, 可是有所差别, 问题三模型是应用于已知影子顶点坐标( ) 及对应北京时间S, 未知经度, 纬度和积日。而本模型是应用于已知直杆原长与影子长度及对应北京时间S, 未知经度, 纬度和积日。两个模型的公式框架基本相同, 只是问题三模型的最终函数式关于太阳方位角, 而本模型的最终函数式关于太阳高度角, 如下式( 4-1) , 整理得到以下方程:

( 4-4)
为一个二元方程, 可按照5.3的对于此类二元方程的模型求解, 解得当地坐标与日期与5.4.1中的当地坐标与日期基本符合。
六、 模型结果分析
本文中各问题的结果为:
问题序号
结果
第一问
——
第二问
( ) 即南海附近和( ) 即俄罗斯中部
第三问
附件2对应的当地坐标为( ) 即巴基斯坦境内, 日期为7月21日,
解得附件3中当地坐标为( ) 即陕西省内, 日期为1月15日。

第四问
( ) 即印度尼西亚附近, ( ) 即新疆地区

为了便于在Matlab中编程运算, 本文模型的运算过程均统一采用弧度制, 在最终结果的表示中采用了更为方便的角度来表示。
本文在问题一模型和问题二模型中引入了基于常规公式的三个订正公式, 即经度订正, 时刻订正和年度订正, 避免了真太阳时与平太阳时混淆而产生的误差, 避免了因每年真实天数为365.2422天可是计入365天( 或366天) 而产生的误差,
本文结果是经过回带检验, 不断修正后所确定, 问题三模型中为了减少未知量数而舍弃了年度订正, 与时间订正, 理论上会有一定误差, 但控制在接受范围内。
针对于问题二和问题三, 附件中有多组数据可取, 问题二模型和问题三模型在计算结果后分别更换带入数据, 再次求解模型, 结果稳定性比较好。

七、 模型评价
问题一中根据给定参数, 运用matlab拟合出直杆的太阳影子长度变化曲线, 其解法为问题二、 三、 四提供了求解经度的方法, 我们还在结果的基础上进行模型分析改进优化, 得到的结果较合理。
问题二建立了阴影轨迹几何模型, 将求解经度转化为先求太阳高度角的间接计算, 将复杂的计算分部计算简单化, 根据影子顶点坐标数据, 在经度确定的情况下给出两个可能的观测点纬度值。
问题三考虑到附件中没有给出观测日期, 即赤纬角未知, 于是在原来的基础上引入新的解题方法——穷举法, 来计算出观测地合理地日期和位置, 在原来建立的模型上不断优化。最终经过检验得到的结果比较理想。
问题四在处理视频中我们引入了背景剪切技术——深度优先搜索——主元分析三步法得出了视频中直杆的影长等数据, 继而用到了问题2、 3的解题方法, 让我们所建立的模型更趋向于应用到实际生活, 最终所得结果也基本符合实际检验。
综上我们认为我们模型的不足之处在于人为计算的工作量较大, 编程也难以做到精简, 对一些参数的求解公式采取了简化近似, 对一些影响不大的因素采取了假设或者忽略。综上可得模型不能给出最精确的解, 有较多的计算过程带有一定的主观性。

八、 模型改进方向
解决太阳影子定位问题所建立的模型为阴影轨迹几何模型, 该模型是以直杆底端为原点所建立平面和空间坐标系进行求解的, 故该模型具有局限性, 一旦题目中出现更切实际求摄像机所在位置的问题, 我们所建立的模型便束手无策。
因此模型的改进能够将建立坐标系的原点放在摄像机镜头所在位置, 从视频图像中仅经过四组阴影观测点和恢复地平线来进行相机校准, 经过建立日晷设计原理的模型来得到阴影轨迹的二次曲线,并经过此曲线来估计观测地的具体位置, 在相机校准过程中不需要投影物体的可见性, 也就是说不在需要影长做参数, 使改进模型更加趋向于我们的现实。

参考文献
[1]武琳. 基于太阳阴影轨迹的经纬度估计技术研究[D].天津大学, .
[2]王炳忠. 太阳辐射计算讲座 第一讲 太阳能中天文参数的计算[J]. 太阳能,1999,02:8-10.

附录
附录内容:
%影长关于时间的函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function l=yingchang(s)
N=295;
jd=116.3913*2*pi/360;
phi=39.9072*2*pi/360;
N0=79.6764+0.2422*( -1985)-floor(( -1985)/4);
t=N-N0;
sita=2*pi*t/365.2422;
deta=23.45*sin(2*pi*(284+N)/365)*2*pi/360;
Et=0.0028-1.9875*sin(sita)+9.9059*sin(2*sita)-7.0924*cos(sita)-0.6882*cos(2*sita);
sd=s-((2*pi/3-jd)*4)/60;%数组
s0=sd+Et/60;
tao=(s0-12)*15*2*pi/360;
h=asin(sin(deta)*sin(phi)+cos(deta)*cos(phi).*cos(tao));
l=3.*cot(h);
end

%问题二
N=172;
jd=2.025;
syms phi;
% x=1.2448;
% y=0.5331;
% s=15;
x=1.0365;
y=0.4973;
s=14.7;

N0=79.6764+0.2422*( -1985)-floor(( -1985)/4);
t=N-N0;
sita=2*pi*t/365.2422;
deta=(0.3723+23.2567*sin(sita)+0.1149*sin(2*sita)-0.1712*sin(3*sita)-0.758*cos(sita)-0.3656*cos(2*sita)+0.0201*cos(3*sita))*2*pi/360;
Et=0.0028-1.9875*sin(sita)+9.9059*sin(2*sita)-7.0924*cos(sita)-0.6882*cos(2*sita);
sd=s-((2*pi/3-jd)*4)/60;
s0=sd+Et/60;
tao=(s0-12)*15*2*pi/360;
h=asin(sin(deta)*sin(phi)+cos(deta)*cos(phi)*cos(tao));
EQ=x/sqrt(x^2+y^2)-(sin(h)*sin(phi)-sin(deta))/(cos(h)*cos(phi));
solu=solve(EQ,phi)

%问题三
digits(2)
format short;
jd=1.2555;
syms phi;
data=xlsread('2.xls','fujian2');
contain1=cell(21,1);
contain2=cell(21,1);
[m,n]=size(data);
rN=1:366;
flag=0;
for i=1:m;
s=data(i,1);
x=data(i,2);
y=data(i,3);
k=0;
p=0;
flag=flag+1;
for N=rN
N
deta=23.45*sin(2*pi*(284+N)/365)*2*pi/360;
sd=s-((2*pi/3-jd)*4)/60;
tao=(sd-12)*15*2*pi/360;
h=asin(sin(deta)*sin(phi)+cos(deta)*cos(phi).*cos(tao));
EQ=x/sqrt(x^2+y^2)-(sin(h)*sin(phi)-sin(deta))/(cos(h)*cos(phi));
solu=solve(EQ,phi)
if(solu<=pi/2&&solu>=-pi/2)
if(flag==1)
k=k+1;
rphi1(k)=solu;
rN1(k)=N;
end
if(flag>1)
fanhui=find(rphi==solu);
if(~isempty(fanhui))
p=p+1;
rphi1(p)=solu;
rN1(p)=N;
end
end
end
end
rN=rN1;
if(isempty(rN))
disp('无公共解');
return;
end
if(isempty(rN1))
disp('无公共解');
return;
end
contain{i,1}=rphi1;
contain{i,1}=rN1;
rphi=rphi1;
rphi1=[];
rN1=[];
end
rN %公共积日
rphi %公共纬度

%拟合附件1数据
data=xlsread('1.xls','fujian1');
x=data(:,1);
y=data(:,2);
p=polyfit(x,y,2)
center=-p(2)/(2*p(1))
xx=center-3:0.01:center+3;
yy=polyval(p,xx);
plot(x,y,'r*',xx,yy,'b-');
%求经度
et=center-12;
jd=(120-et*15)*pi/180

%拟合附件2

data=xlsread('1.xls','fujian2');
x=data(:,1);
y=data(:,2);
p=polyfit(x,y,2)
center=-p(2)/(2*p(1))
xx=center-3:0.01:center+3;
yy=polyval(p,xx);
plot(x,y,'r*',xx,yy,'b-');
%求经度
et=center-12;
jd=(120-et*15)*pi/180

%拟合附件3

data=xlsread('1.xls','fujian3');
x=data(:,1);
y=data(:,2);
p=polyfit(x,y,2)
center=-p(2)/(2*p(1))
xx=center-3:0.01:center+3;
yy=polyval(p,xx);
plot(x,y,'r*',xx,yy,'b-');
%求经度
et=center-12;
jd=(120-et*15)*pi/180

Digimizer 软件进行图像处理结果