全國大學生數(shù)學建模競賽論文-關于太陽影子定位的模型研究

關于太陽影子定位的模型研究摘要太陽影子定位技術是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點與日期的一種方法。如何確定視頻的拍攝地點和拍攝日期是視頻數(shù)據(jù)分析的重要方面。而本文主要研究的就是通過太陽影子的變化來確定物體地點和日期的問題。針對問題一，要確定直桿太陽影子長度的變化曲線，本文用太陽高度角、赤緯角、時角、經(jīng)度和緯度之間的關系建立函數(shù)模型，利用MATLAB擬合出影子長度隨時間的變化曲線。針對問題二，要求根據(jù)直桿影子頂點坐標求出直桿所在地的位置，根據(jù)經(jīng)緯度等變量之間的關系，建立函數(shù)模型：通過MATLAB中的fsolve方法求解出直桿所在地的經(jīng)緯度可知直桿所在地為海南省昌化港。并利用MATLAB擬合得到了直桿影長與時間之間的函數(shù)，進而驗證了該模型的合理性。針對問題三，要求根據(jù)直桿影子頂點坐標求出直桿所在地的位置與日期，我們對問題二的模型進行反推得到附件2直桿所在經(jīng)緯度，為6月1日的新疆阿克蘇；附件三直桿所在經(jīng)緯度，為4月24日的蒙古賽音山達附近。針對問題四，要求確定視頻的拍攝位置與日期，我們用CAD制圖軟件處理截取圖像中相對桿長與相對影長，結合問題二中模型得到視頻的拍攝地點為山西呂梁或內蒙古包頭。對于拍攝時間未知的問題，利用類似方法得到直桿所在地為6月30日的內蒙古鄂爾多斯或6月29日的內蒙古包頭。關鍵詞：太陽影子定位；最小二乘法；太陽高度角；赤緯角問題提出如何確定視頻的拍攝地點和拍攝日期是視頻數(shù)據(jù)分析的重要方面。太陽影子定位技術就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法。根據(jù)太陽影子定位技術建立模型，解決下列問題：1.建立影子長度變化的數(shù)學模型，分析影子長度關于各個參數(shù)的變化規(guī)律，并應用建立的模型畫出2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯東經(jīng)）3米高直桿的太陽影子長度的變化曲線。2.根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點。將模型應用于附件1的影子頂點坐標數(shù)據(jù)中，給出若干個可能的地點。3.根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點和日期。將模型分別應用于附件2和附件3的影子頂點坐標數(shù)據(jù)，給出若干個可能的地點與日期。4．附件4為一根高度大約2米的直桿在太陽下的影子變化的視頻。建立確定視頻拍攝地點的數(shù)學模型，并用模型給出若干個可能的拍攝地點。如果拍攝日期未知，能否根據(jù)視頻確定出拍攝地點與日期。二、問題分析針對問題一，要求建立影子長度變化的模型，分析影子長度關于各個參數(shù)的變化規(guī)律，通過應用太陽高度角、赤緯角、時角、經(jīng)度和緯度之間的關系建立模型，利用MATLAB擬合出影子長度隨時間變化的曲線。針對問題二，要求根據(jù)直桿影子頂點坐標求出直桿所在地的位置。對于所求問題，我們建立關于高度角、赤緯角、經(jīng)緯度等變量之間的函數(shù)模型，采用MATLAB中fsolve的方法求解出直桿可能所在的位置。接著我們用MATLAB最小二乘法擬合曲線驗證了該模型的合理性。針對問題三，要求根據(jù)直桿影子頂點坐標求出直桿所在地的位置與日期。我們在問題二模型的基礎上，重新設置變量建立了新的函數(shù)模型，用MATLAB求解出直桿可能所在位置與日期。針對問題四，要求根據(jù)一根已知高度的直桿在太陽下影子變化的視頻確定視頻的拍攝位置與日期。我們截取了14個時刻的視頻圖片，用CAD制圖軟件測出直桿的相對影長與桿長(見附件10)，并建立函數(shù)模型，得到視頻拍攝的位置與日期。三、模型假設1.忽略地球大氣層折射效應對太陽光的影響。2.忽略視覺對桿長與影長數(shù)據(jù)收集的影響。3.忽略視頻拍攝角度對相對影長與相對桿長的影響。定義與符號說明符號含義太陽高度角緯度太陽赤緯角時角時間序號（從1月1開始計時)所求地方時直桿長影長已知地方時所求地點經(jīng)度時間模型的建立與求解5.1問題一模型的建立與求解圖1-1赤緯角示意圖圖1-2高度角示意圖對于赤緯角（1）由題意知，直桿影長變化的時間是2015年北京時間10月22日，可取，由(1)得對于地方時（2）地球上經(jīng)度每相隔，地方時相差4分鐘，故有地方時與已知地方時的關系（3）直桿所在地點是北緯39度54分26秒,東經(jīng)116度23分29秒，所以可以由（2）得出從9:00到15:00各個時刻的。對于太陽高度角（4）直桿影長、桿長、太陽高度角幾何關系:（5）然后，我們利用MATLAB軟件得到了直桿影長在中隨時間變化的曲線，如圖1-3所示（源文件見附表（1））。圖1-3影長隨時間的變化曲線5.2問題二模型的建立與求解5.21模型的建立與求解由題意得測量日期是2015年4月18日，可算出日期序號，由（1）得：設直桿長，影長，而與影子頂點坐標有：（6）結合(5)得：（7）由(2)、(3)可得：（8）結合公式（4）得到：（9）下面我們假設函數(shù)：（10）而后在附件1中選取六組數(shù)據(jù)分別代入（10），利用MATLAB中最小二乘法來求出桿長、經(jīng)度與緯度，從而確定直桿所處的地點。所得情況如圖2-1(源程序見附表(2)）。序號北京時間X坐標（米）Y坐標（米）影長（米）經(jīng)/緯度（度）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.182198976108.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.28319534108.3289/19.519615:181.47510.56571.57985331615:331.68820.59521.790050915圖2-1每組數(shù)據(jù)得到的經(jīng)緯度數(shù)值基本一致。所以附件1直桿的經(jīng)緯度坐標可取，使用谷歌地球得知，此地處于海南省昌化港(圖2-2)。圖2-25.22模型的驗證根據(jù)附件1中某直桿在水平地面上的太陽影子的頂點坐標，可以通過Excel軟件求出直桿在在個時間段的影子長度(圖2-3)。圖2-3通過使用最小二乘法，由MATLAB擬合出從14:42到15:42里影長隨時間的變化曲線，如圖2-4（源程序見附表(3)）。圖2-4從而得擬合曲線函數(shù)為：（11）而后以3分鐘為間隔，擬合影長從10:00到16:00的隨時間的變化曲線，如圖2-5（源程序見附表(4)）。biao圖2-5由圖2-5得出，當時，影長最短米在附件1中，北京時間恰是該地區(qū)正午時間，而北京時間所在經(jīng)度為東經(jīng)，故在公式(3)中：，,經(jīng)度差=可得所求經(jīng)度：由于南北回歸線相差緯度是，一年中太陽在南北回歸線范圍內移動兩次，所以一年太陽直射點的移動緯度是?，F(xiàn)取一年365天，所以太陽每天移動緯度為。由于春分時的太陽直射點緯度，2015年春分是3月21日，而附件1的時間為2015年4月18日，可知此時的太陽直射點的緯度：根據(jù)正午太陽高度的計算公式：（12）是要測的緯度，是太陽直射點的緯度。在模型的求解中已得桿長米，結合(5)、(12)可得：得緯度：所以附件1直桿的經(jīng)緯度坐標為，使用谷歌地球得知，此地處于廣東省茂名市，如圖2-6。圖2-6驗證結果與模型基本吻合，故模型的建立較為合理。5.3問題三模型的建立與求解設直桿長，影長，而與影子頂點坐標有：結合(5)可得：由(2)、(3)可知時角結合(4)得到：（13）結合(1)建立模型：（14）在附件2中選取八組數(shù)據(jù)代入函數(shù)（14），在MATLAB中用最小二乘法求出時間序號、桿長、經(jīng)度與緯度，從而確定直桿所處的位置與日期。如圖3-1(源程序見附表(5)）。序號北京時間X坐標（米）Y坐標（米）影長（米）經(jīng)/緯度（度）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.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.880974612: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.06508113:20-0.89650.36190.9667913:41-0.72270.44840.850504圖3-1將以上八組數(shù)據(jù)得到的經(jīng)緯度擬合到坐標系中，如圖3-2（源程序見附表(6)）。圖3-2所以附件2直桿的日期可由時間序號的意義求得為6月1日。由圖3-2得經(jīng)緯度坐標為，使用谷歌地球得知，此地處于新疆阿克蘇。如圖3-3。圖3-3用與處理附件2相同的方法，在附件3中隨機取出兩組數(shù)據(jù)，所得情況如圖3-4(源程序見附表(7)）。序號北京時間X坐標（米）Y坐標（米）影長（米）經(jīng)/緯度（度）113:091.16373.3363.533142184110.0119/43.656613:271.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日。經(jīng)緯度坐標為，使用谷歌地球得知，該地在蒙古賽音山達附近。如圖3-5。圖3-55.4問題四模型的建立與求解附件4中給出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.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）忽略視頻拍攝角度的影響，根據(jù)圖4-1的數(shù)據(jù)建立如下模型：（16）在圖4-1中選取三組數(shù)據(jù)，通過MATLAB求解出經(jīng)度與緯度（圖4-2）(見附表(8))，從而確定直桿所在地為山西省呂梁市（圖4-3）或內蒙古包頭（圖4-4）。序號北京時間經(jīng)度（度）緯度（度）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對于問題四中提出的“如果拍攝日期未知，能否根據(jù)視頻確定拍攝的地點和日期”的問題，我們建立以下模型：（17），其中緯度、經(jīng)度以及時間序號是未知量。在圖4-1中選取兩組數(shù)據(jù)，通過MATLAB求解出緯度、經(jīng)度以及時間序號，如圖4-5(見附表(9)）：序號北京時間經(jīng)度（度）緯度（度）時間序號18:581.164912109.633040.4159179.78229:131.0565739:310.92982529: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模型的優(yōu)點問題一中我們根據(jù)太陽高度角、赤緯角等相關量之間的關系，利用附件1中的直桿影子頂點坐標建立較為恰當?shù)哪Ｐ?，解決了關于直桿影子長度變化曲線的問題，所得曲線較為合理。問題二中我們根據(jù)直桿影子頂點坐標，通過合理的假設建立了函數(shù)模型，并用最小二乘法擬合驗證了該模型的合理性。將問題二的模型反推即可得到問題三的結果。問題四我們采用CAD制圖軟件對截取的14個時刻的數(shù)據(jù)進行處理分析，利用MATLAB中的fsolve方法得出了所求結果。6.2模型的缺點問題一的解決中忽略了大氣層折射等因素對太陽高度角的影響，所建立模型的結果與實際情況有誤差。問題四中沒有考慮視頻拍攝角度所造成的影響，未能得到最理想的結果。七、參考文獻[1]郭大為，數(shù)學建模[M],合肥，安徽教育出版社,2014(5)。[2]陳健婷,溫銀婷,傅守忠等,最佳太陽方位角計算[J],肇慶學院學報,2014(2):20-21。[3]王國安,米鴻濤,李蘭霞等,太陽高度角和日出日落時刻太陽方位角一年變化范圍的計算[J],氣象與環(huán)境科學,2007,35(增刊):161-163。[4]張闖,呂東輝等太陽實時位置計算及在圖像光照方向的應用[J],電子測量技術,2010(11):87-89。[5]張文華,司德亮,徐淑通等,太陽影子倍率的計算方法及其對光伏陣列布局的影響[M],技術與產(chǎn)品,2011(9):28-30。[6]劉二根,王廣超,朱旭生等,Matlab與數(shù)學實驗[M],北京,國防工業(yè)出版社,2014(1):180-187。附錄附表一（1）圖1-3程序代碼t=9:0.02:15;a=23.45*sin((2*pi*579)/365);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程序代碼functionf=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.149626^2+x(3)^2);f(2)=sind(x(1))*sind(10.511)+cosd(x(1))*cosd(10.511)*cosd(x(2)-72.75)-x(3)/sqrt(1.4634^2+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.620145^2+x(3)^2);fsolve(@fun11,[30,110,2])functionf=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.182199^2+x(3)^2);f(2)=sind(x(1))*sind(10.511)+cosd(x(1))*cosd(10.511)*cosd(x(2)-71.25)-x(3)/sqrt(1.540232^2+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.880875^2+x(3)^2);fsolve(@fun11,[30,110,2])functionf=fun(x)f(1)=sind(x(1))*sind(10.511)+cosd(x(1))*cosd(10.511)*cosd(x(2)-78)-x(3)/sqrt(1.215297^2+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.353364^2+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.661271^2+x(3)^2);fsolve(@fun11,[30,110,2])functionf=fun(x)f(1)=sind(x(1))*sind(10.511)+cosd(x(1))*cosd(10.511)*cosd(x(2)-77.25)-x(3)/sqrt(1.249051^2+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.501482^2+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.927918^2+x(3)^2);fsolve(@fun11,[30,110,2])functionf=fun(x)f(1)=sind(x(1))*sind(10.511)+cosd(x(1))*cosd(10.511)*cosd(x(2)-76.5)-x(3)/sqrt(1.283195^2+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.579853^2+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.790051^2+x(3)^2);fsolve(@fun11,[30,110,2])(3)圖2-4程序代碼x=[14.7,14.75,14.8,14.85,14.9,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.790050915,1.835014272,1.880875001,1.927918447];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r');xlabel('時刻')ylabel('影長')圖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('影長')圖3-1程序代碼functionf=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.247256^2+(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.044446^2+(x(3)^2));f(3)=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.880974^2+(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.15244^2+(x(3)^2));fsolve(@fun1,[30,80,3,170])functionf=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.005)-x(3)/sqrt(1.22795^2+(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.086254^2+(x(3)^2));f(3)=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)-100.0005)-x(3)/sqrt(0.96679^2+(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.850504^2+(x(3)^2));fsolve(@fun2,[30,80,3,170])functionf=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.247256^2+(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.129917^2+(x(3)^2));f(3)=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.005)-x(3)/sqrt(1.044446^2+(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.930928^2+(x(3)^2));fsolve(@fun3,[30,80,3,170])functionf=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)-x(3)/sqrt(1.222795^2+(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.065081^2+(x(3)^2));f(3)=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.5005)-x(3)/sqrt(0.930928^2+(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.850504^2+(x(3)^2));fsolve(@fun4,[30,80,3,170])functionf=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)-107.5)-x(3)/sqrt(1.175429^2+(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.086254^2+(x(3)^2));f(3)=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)-99.2505)-x(3)/sqrt(0.948585^2+(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.880974^2+(x(3)^2));fsolve(@fun4,[30,80,3,170])functionf=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.247256^2+(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.75)-x(3)/sqrt(1.15244^2+(x(3)^2));f(3)=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)-100.7505)-x(3)/sqrt(0.985491^2+(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.880974^2+(x(3)^2));fsolve(@fun6,[30,80,3,170])functionf=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.247256^2+(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.75)-x(3)/sqrt(1.15244^2+(x(3)^2));f(3)=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.005)-x(3)/sqrt(1.044446^2+(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.5005)-x(3)/sqrt(0.930928^2+(x(3)^2));fsolve(@fun7,[30,80,3,170])functionf=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)-x(3)/sqrt(1.222795^2+(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.065081^2+(x(3)^2));f(3)=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)-100.0005)-x(3)/sqrt(0.96679^2+(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.850504^2+(x(3)^2));fsolve(@fun8,[30,80,3,170])(6)圖3-2程序代碼x=[40.739838.335041.047939.098239.441339.994341.459439.0069]y=[80.094680.720381.278179.488279.556879.698781.239579.4539]plot(y,x,'k+');xlabel('緯度X')ylabel('經(jīng)度Y')title('經(jīng)緯度坐標示意圖')（7）圖3-4程序代碼functionf=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)-102.75)-x(3)/sqrt(3.533142^2+(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)-102)-x(3)/sqrt(3.546768^2+(x(3)^2));f(3)=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)-101.25)-x(3)/sqrt(3.561798^2+(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)-100.5)-x(3)/sqrt(3.578101^2+(x(3)^2));fsolve(@fun1,[40,110,4,100])functionf=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)-102.75)-x(3)/sqrt(3.53142^2+(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)-98.25)-x(3)/sqrt(3.635426^2+(x(3)^2));f(3)=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)-93.75)-x(3)/sqrt(3.788088^2+(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)-87.75)-x(3)/sqrt(4.077863^2+(x(3)^2));fsolve(@fun2,[40,110,4,100])(8)圖4-2程序代碼functionh=fun(t)h