2015全国大学生数学建模竞赛论文.docx

《2015全国大学生数学建模竞赛论文.docx》由会员分享，可在线阅读，更多相关《2015全国大学生数学建模竞赛论文.docx（32页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、关于太阳影子定位的模型探讨摘要太阳影子定位技术是通过分析视频中物体的太阳影子改变，确定视频拍摄的地点与日期的一种方法。如何确定视频的拍摄地点与拍摄日期是视频数据分析的重要方面。而本文主要探讨的就是通过太阳影子的改变来确定物体地点与日期的问题。针对问题一，要确定直杆太阳影子长度的改变曲线，本文用太阳高度角、赤纬角、时角、经度与纬度之间的关系建立函数模型，利用MATLAB拟合出影子长度随时间的改变曲线。针对问题二，要求依据直杆影子顶点坐标求出直杆所在地的位置，依据经纬度等变量之间的关系，建立函数模型：通过MATLAB中的fsolve方法求解出直杆所在地的经纬度可知直杆所在地为海南省昌化港。并利用M

2、ATLAB拟合得到了直杆影长与时间之间的函数，进而验证了该模型的合理性。针对问题三，要求依据直杆影子顶点坐标求出直杆所在地的位置与日期，我们对问题二的模型进展反推得到附件2直杆所在经纬度，为6月1日的新疆阿克苏；附件三直杆所在经纬度，为4月24日的蒙古赛音山达旁边。针对问题四，要求确定视频的拍摄位置与日期，我们用CAD制图软件处理截取图像中相对杆长与相对影长，结合问题二中模型得到视频的拍摄地点为山西吕梁或内蒙古包头。对于拍摄时间未知的问题，利用类似方法得到直杆所在地为6月30日的内蒙古鄂尔多斯或6月29日的内蒙古包头。关键词：太阳影子定位； 最小二乘法；太阳高度角；赤纬角 一、 问题提出如何确

3、定视频的拍摄地点与拍摄日期是视频数据分析的重要方面。太阳影子定位技术就是通过分析视频中物体的太阳影子改变，确定视频拍摄的地点与日期的一种方法。依据太阳影子定位技术建立模型，解决下列问题：1.建立影子长度改变的数学模型，分析影子长度关于各个参数的改变规律，并应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬东经）3米高直杆的太阳影子长度的改变曲线。2.依据某固定直杆在程度地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将模型应用于附件1的影子顶点坐标数据中，给出若干个可能的地点。3. 依据某固定直杆在程度地面上的太阳影子顶点坐标数据，建立数学模

4、型确定直杆所处的地点与日期。将模型分别应用于附件2与附件3的影子顶点坐标数据，给出若干个可能的地点与日期。4附件4为一根高度大约2米的直杆在太阳下的影子改变的视频。建立确定视频拍摄地点的数学模型，并用模型给出若干个可能的拍摄地点。假如拍摄日期未知，能否依据视频确定出拍摄地点与日期。二、问题分析针对问题一，要求建立影子长度改变的模型，分析影子长度关于各个参数的改变规律，通过应用太阳高度角、赤纬角、时角、经度与纬度之间的关系建立模型，利用MATLAB拟合出影子长度随时间改变的曲线。针对问题二，要求依据直杆影子顶点坐标求出直杆所在地的位置。对于所求问题，我们建立关于高度角、赤纬角、经纬度等变量之间的

5、函数模型，采纳MATLAB中 fsolve的方法求解出直杆可能所在的位置。接着我们用MATLAB最小二乘法拟合曲线验证了该模型的合理性。针对问题三，要求依据直杆影子顶点坐标求出直杆所在地的位置与日期。我们在问题二模型的根底上，重新设置变量建立了新的函数模型，用MATLAB求解出直杆可能所在位置与日期。针对问题四，要求依据一根已知高度的直杆在太阳下影子改变的视频确定视频的拍摄位置与日期。我们截取了14个时刻的视频图片，用CAD制图软件测出直杆的相对影长与杆长(见附件10)，并建立函数模型，得到视频拍摄的位置与日期。三、模型假设 1.忽视地球大气层折射效应对太阳光的影响。 2.忽视视觉对杆长与影长

6、数据搜集的影响。 3.忽视视频拍摄角度对相对影长与相对杆长的影响。四、 定义与符号说明符号含义太阳高度角纬度太阳赤纬角时角时间序号（从1月1开场计时)所求地方时直杆长影长已知地方时所求地点经度时间五、 模型的建立与求解5.1 问题一模型的建立与求解 图1-1 赤纬角示意图 图1-2 高度角示意图 对于赤纬角 （1）由题意知，直杆影长改变的时间是2015年北京时间10月22日，可取 ，由 (1) 得对于地方时 （2）地球上经度每相隔，地方时相差4分钟，故有地方时与已知地方时的关系 （3）直杆所在地点是北纬39度54分26秒,东经116度23分29秒，所以可以由（2）得出从9:00到15:00各个

7、时刻的。对于太阳高度角 （4）直杆影长 、杆长 、太阳高度角 几何关系: （5） 然后，我们利用MATLAB软件得到了直杆影长在 中随时间改变的曲线，如图1-3所示（源文件见附表（1）。图1-3 影长随时间的改变曲线5.2 问题二模型的建立与求解5.21 模型的建立与求解由题意得测量日期是 2015年4月18日，可算出日期序号，由（1）得：设直杆长 ，影长 ，而 与影子顶点坐标 有： （6） 结合 (5)得： （7） 由(2)、(3)可得： （8） 结合公式（4）得到： （9） 下面我们假设函数： （10） 而后在附件1中选取六组数据分别代入（10），利用MATLAB中最小二乘法来求出杆长 、

8、经度 与纬度 ，从而确定直杆所处的地点。所得状况如图2-1(源程序见附表 (2)）。序号北京时间X坐标（米）Y坐标（米）影长（米）经/纬度（度）114:421.03650.49731.149625826108.5990/19.315015:001.24480.53111.35336404915:301.64380.58921.74620591214:421.03650.49731.149625826108.7485/19.176915:091.35680.54831.46339985315:211.5160.57151.620144515314:451.06990.50291.18219897

9、6108.6433/19.270515:151.43490.55981.54023181715:391.78010.60741.880875001414:481.10380.50851.215296955108.5369/19.358215:001.24480.53111.35336404915:241.55770.57741.661270613514:511.13830.51421.249051052108.7285/19.219215:121.39550.55411.50148162215:421.82770.61351.927918447614:541.17320.51981.28319

10、534108.3289/19.519615:181.47510.56571.57985331615:331.68820.59521.790050915图2-1每组数据得到的经纬度数值根本一样。所以附件1直杆的经纬度坐标可取，运用谷歌地球得知，此地处于海南省昌化港(图2-2)。图2-25.22 模型的验证 依据附件1中某直杆在程度地面上的太阳影子的顶点坐标，可以通过Excel软件求出直杆在在个时间段的影子长度(图2-3)。图2-3通过运用最小二乘法，由MATLAB拟合出从14:42到15:42里影长随时间的改变曲线，如图2-4（源程序见附表(3)）。图2-4从而得拟合曲线函数为： （11）而后以

11、3分钟为间隔，拟合影长从10:00到16:00的随时间的改变曲线，如图2-5（源程序见附表 (4)）。biao图2-5由图2-5得出，当时，影长最短米 在附件1中，北京时间恰是该地区正午时间，而北京时间所在经度为东经，故在公式 (3)中：，,经度差=可得所求经度：由于南北回来线相差纬度是，一年中太阳在南北回来线范围内挪动两次，所以一年太阳直射点的挪动纬度是。现取一年365天，所以太阳每天挪动纬度为。由于春分时的太阳直射点纬度，2015年春分是3月21日，而附件1的时间为2015年4月18日，可知此时的太阳直射点的纬度：依据正午太阳高度的计算公式： （12）是要测的纬度，是太阳直射点的纬度。 在

12、模型的求解中已得杆长 米，结合 (5)、(12)可得： 得纬度： 所以附件1直杆的经纬度坐标为，运用谷歌地球得知，此地处于广东省茂名市，如图2-6。图2-6验证结果与模型根本吻合，故模型的建立较为合理。 5.3 问题三模型的建立与求解设直杆长 ，影长 ，而 与影子顶点坐标有： 结合 (5)可得： 由 (2)、(3)可知时角 结合 (4)得到： （13） 结合 (1)建立模型： （14） 在附件2中选取八组数据代入函数（14），在MATLAB中用最小二乘法求出时间序号 、杆长 、经度 与纬度 ，从而确定直杆所处的位置与日期。如图 3-1(源程序见附表 (5)）。序号北京时间X坐标（米）Y坐标（米

13、）影长（米）经/纬度（度）112:41-1.23520.1731.24725680.0946/40.739812:53-1.12810.23561.1524413:08-0.9980.3081.0444613:35-0.77190.42460.880974212:44-1.20810.1891.22279580.7203/38.335013:02-1.04960.27981.08625413:20-0.89650.36190.9667913:41-0.72270.44840.850504312:41-1.23520.1731.24725681.2781/41.047912:56-1.10180

14、.25051.12991713:08-0.9980.3081.04444613:26-0.84640.38760.930928412:44-1.20810.1891.22279579.4882/39.098213:05-1.02370.2941.06508113:26-0.84640.38760.93092813:41-0.72270.44840.850504512:50-1.15460.22031.17542979.5568/39.441313:02-1.04960.27981.08625413:23-0.87140.37480.94858513:35-0.77190.42460.88097

15、4612:41-1.23520.1731.24725679.6987/39.994312:53-1.12810.23561.1524413:17-0.92170.34880.98549113:35-0.77190.42460.880974712:41-1.23520.1731.24725681.2395/41.459412:53-1.12810.23561.1524413:08-0.9980.3081.04444613:26-0.84640.38760.930928812:44-1.20810.1891.22279579.4539/39.006913:05-1.02370.2941.06508

16、113:20-0.89650.36190.9667913:41-0.72270.44840.850504图 3-1将以上八组数据得到的经纬度拟合到坐标系中，如图3-2（源程序见附表 (6)）。图 3-2所以附件2直杆的日期可由时间序号的意义求得为6月1日。由图 3-2得经纬度坐标为，运用谷歌地球得知，此地处于新疆阿克苏。如图3-3。图 3-3用与处理附件2一样的方法，在附件3中随机取出两组数据，所得状况如图 3-4(源程序见附表 (7)）。序号北京时间X坐标（米）Y坐标（米）影长（米）经/纬度（度）113:091.16373.3363.533142184110.0119/43.656613:2

17、71.51483.30483.63542598313:451.88483.28593.78808788814:092.42133.28124.077863059213:091.16373.3363.533142184109.3294/22.458213:121.22123.32993.54676802913:151.27913.32423.56179764313:181.33733.31883.578100715图 3-4所以附件3直杆的日期可由时间序号的意义求得为4月24日。经纬度坐标为，运用谷歌地球得知，该地在蒙古赛音山达旁边。如图 3-5。图 3-55.4 问题四模型的建立与求解 附件4

18、中给出2015年7月13日某地一根直杆在太阳下影子改变的视频（即）。为了确定视频中直杆的杆长与影长之比，用CAD制图截取14个时刻的视频图片（从8:55-9:34每隔三分钟截取一次），从而得到了杆长与影长之比，如图 4-1。时刻影长杆长8:5513.3711.281.1852836888:5813.2811.41.1649122819:0151.0645.111.13189989:0461.1755.831.0956475019:0758.7553.951.088971279:1071.2166.911.0642654319:1370.4166.641.0565726299:1668.6966

19、.471.0333985269:1967.3767.061.0046227269:2266.4666.680.996700669:2565.5366.530.9849691879:2847.8751.020.9382595069:3145.5849.020.9298245619:3444.2748.960.904207516图 4-1 由 (1)、(2)、(3)、(4)得： （15）忽视视频拍摄角度的影响，依据图 4-1的数据建立如下模型： （16）在图 4-1中选取三组数据，通过MATLAB求解出经度与纬度（图 4-2）(见附表(8)，从而确定直杆所在地为山西省吕梁市（图 4-3）或内蒙古包

20、头（图 4-4）。序号北京时间经度（度）纬度（度）18:551.185283688110.868637.26549:340.90420751629:011.1318998110.999637.07629:280.93825950638:581.164912281111.114241.43939:220.99670066图4-2图 4-3图 4-4对于问题四中提出的“假如拍摄日期未知，能否依据视频确定拍摄的地点与日期”的问题，我们建立以下模型： （17），其中纬度、经度以刚好间序号是未知量。在图 4-1中选取两组数据，通过MATLAB求解出纬度、经度以刚好间序号，如图 4-5(见附表(9)）：序

21、号北京时间经度（度）纬度（度）时间序号18:581.164912109.633040.4159179.78229:131.0565739:310.929825 29:041.095648110.242341.7079178.76029:191.0046239:340.904208图 4-5从而确定视频拍摄于6月30日，内蒙古鄂尔多斯市（图4-6）或6月29日内蒙古包头（图4-7）。图 4-6图 4-7六、 模型的评价6.1 模型的优点问题一中我们依据太阳高度角、赤纬角等相关量之间的关系，利用附件1中的直杆影子顶点坐标建立较为恰当的模型，解决了关于直杆影子长度改变曲线的问题，所得曲线较为合理。问

22、题二中我们依据直杆影子顶点坐标，通过合理的假设建立了函数模型，并用最小二乘法拟合验证了该模型的合理性。将问题二的模型反推即可得到问题三的结果。问题四我们采纳CAD制图软件对截取的14个时刻的数据进展处理分析，利用MATLAB中的fsolve方法得出了所求结果。6.2 模型的缺点问题一的解决中忽视了大气层折射等因素对太阳高度角的影响，所建立模型的结果与实际状况有误差。问题四中没有考虑视频拍摄角度所造成的影响，未能得到最志向的结果。七、参考文献1 郭大为，数学建模M,合肥，安徽教化出版社,2014(5)。2 陈健婷,温银婷,傅守忠等,最佳太阳方位角计算J,肇庆学院学 报,2014(2):20-21

23、。3 王国安,米鸿涛,李兰霞等,太阳高度角与日出日落时刻太阳方位 角一年改变范围的计算J,气象与环境科学,2007,35(增 刊):161-163。4 张闯,吕东辉等太阳实时位置计算及在图像光照方向的应用J,电子测量技术,2010(11):87-89。5 张文华,司德亮,徐淑通等,太阳影子倍率的计算方法及其对光伏阵列布局的影响M,技术与产品,2011(9):28-30。6 刘二根,王广超,朱旭生等,Matlab与数学试验M,北京,国防工业出版社,2014(1):180-187。七、 附录附表一（1）图 1-3程序代码t=9:0.02:15;a=23.45*sin(2*pi*579)/365);

24、 b=39.9072;c=(t-(120-116-23/60-29/3600)/15-12)*15;d=sind(a)*sind(b)+cosd(b)*cosd(a)*cosd(c);y=asin(d);z=3*cot(y);plot(t,z,k);xlabel(时刻)ylabel(影长)（2）图 2-1程序代码function f=fun(x)f(1)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-79.5)-x(3)/sqrt(1.1496262+x(3)2);f(2)=sind(x(1)*sind(10.511)+cosd

25、(x(1)*cosd(10.511)*cosd(x(2)-72.75)-x(3)/sqrt(1.46342+x(3)2);f(3)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-69.75)-x(3)/sqrt(1.6201452+x(3)2);fsolve(fun11,30,110,2)function f=fun(x)f(1)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-78.75)-x(3)/sqrt(1.1821992+x(3)2);f(2)=sind(

26、x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-71.25)-x(3)/sqrt(1.5402322+x(3)2);f(3)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-65.25)-x(3)/sqrt(1.8808752+x(3)2);fsolve(fun11,30,110,2)function f=fun(x)f(1)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-78)-x(3)/sqrt(1.2152

27、972+x(3)2);f(2)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-75)-x(3)/sqrt(1.3533642+x(3)2);f(3)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-69)-x(3)/sqrt(1.6612712+x(3)2);fsolve(fun11,30,110,2)function f=fun(x)f(1)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-77.25

28、)-x(3)/sqrt(1.2490512+x(3)2);f(2)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-72)-x(3)/sqrt(1.5014822+x(3)2);f(3)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-64.5)-x(3)/sqrt(1.9279182+x(3)2);fsolve(fun11,30,110,2)function f=fun(x)f(1)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.

29、511)*cosd(x(2)-76.5)-x(3)/sqrt(1.2831952+x(3)2);f(2)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-70.5)-x(3)/sqrt(1.5798532+x(3)2);f(3)=sind(x(1)*sind(10.511)+cosd(x(1)*cosd(10.511)*cosd(x(2)-66.75)-x(3)/sqrt(1.7900512+x(3)2);fsolve(fun11,30,110,2)(3) 图 2-4程序代码x=14.7,14.75,14.8,14.85,14.9

30、,14.95,15.00,15.05,15.1,15.15,15.2,15.25,15.3,15.35,15.4,15.45,15.5,15.55,15.6,15.65,15.7;y=1.149625826,1.182198976,1.215296955,1.249051052,1.28319534,1.317993149,1.353364049,1.389387091,1.426152856,1.463399853,1.501481622,1.540231817,1.579853316,1.620144515,1.661270613,1.703290633,1.74620591,1.7900

31、50915,1.835014272,1.880875001,1.927918447;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r);xlabel(时刻)ylabel(影长)(4) 图2-5程序代码x=10:0.05:16;y=0.1489*x.2-3.7519*x+24.1275;plot(x,y,k)z=polyval(A,x);plot(x,y,k+,x,z,r);xlabel(时刻)ylabel(影长)（5） 图 3-1程序代码function f=fun(x)f(1)=sind(x(1)*sind(23.45*sin(2*pi*(28

32、4+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-109.75)-x(3)/sqrt(1.2472562+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-103.0005)-x(3)/sqrt(1.0444462+(x(3)2);f(3)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd

33、(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-96.2505)-x(3)/sqrt(0.8809742+(x(3)2);f(4)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-106.7505)-x(3)/sqrt(1.152442+(x(3)2);fsolve(fun1,30,80,3,170)function f=fun(x)f(1)=sind(x(1)*sind(23.45*sin

34、(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-109.005)-x(3)/sqrt(1.227952+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-104.5005)-x(3)/sqrt(1.0862542+(x(3)2);f(3)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/

35、365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-100.0005)-x(3)/sqrt(0.966792+(x(3)2);f(4)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-94.7505)-x(3)/sqrt(0.8505042+(x(3)2);fsolve(fun2,30,80,3,170)function f=fun(x)f(1)=sind(x(1)*sind(

36、23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-109.75)-x(3)/sqrt(1.2472562+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-106)-x(3)/sqrt(1.1299172+(x(3)2);f(3)=sind(x(1)*sind(23.45*sin(2*pi*(284+x

37、(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-103.005)-x(3)/sqrt(1.0444462+(x(3)2);f(4)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-98.50058)-x(3)/sqrt(0.9309282+(x(3)2);fsolve(fun3,30,80,3,170)function f=fun(x)f(1)=sind(x(1)*

38、sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-109)-x(3)/sqrt(1.2227952+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-103.7505)-x(3)/sqrt(1.0650812+(x(3)2);f(3)=sind(x(1)*sind(23.45*sin(2*pi

39、*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-98.5005)-x(3)/sqrt(0.9309282+(x(3)2);f(4)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-94.7505)-x(3)/sqrt(0.8505042+(x(3)2);fsolve(fun4,30,80,3,170)function f=fun(x)f(1)=sind

40、(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-107.5)-x(3)/sqrt(1.1754292+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-104.5005)-x(3)/sqrt(1.0862542+(x(3)2);f(3)=sind(x(1)*sind(23.45*

41、sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-99.2505)-x(3)/sqrt(0.9485852+(x(3)2);f(4)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-96.2505)-x(3)/sqrt(0.8809742+(x(3)2); fsolve(fun4,30,80,3,170)function f=fun(x)f(1)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/365)*cosd(x(2)-109.75)-x(3)/sqrt(1.2472562+(x(3)2);f(2)=sind(x(1)*sind(23.45*sin(2*pi*(284+x(4)/365)+cosd(x(1)*cosd(23.45*sin(2*pi*(284+x(4)/36