課件文稿成果

學(xué)習(xí)資料，請關(guān)注淘寶：學(xué)神資料 /認(rèn)豎直桿所在的以及數(shù)據(jù)對應(yīng)的記錄日期。線，并分析長度關(guān)于各個參數(shù)的變化規(guī)律。本文利用天文學(xué)公式，建立了該豎直桿關(guān)于時間變化的計算求解模型，得出影長隨時間的變化曲線如下。此后，本文使用控制變量法，分別對、日期、桿長及當(dāng)日時間這五個參數(shù)第二問要求分析某豎直桿的頂點坐標(biāo)數(shù)據(jù)，確定其地點。本文沿用19度，東經(jīng)109度，位于我國海南省。緯經(jīng)日大致地附件525圖木舒克附件1031陜西，漢中廣元之第四問要求分析豎直桿的變化，分別在日期已知和未知的情況下求緯經(jīng)日大致地日期確713內(nèi)錫林郭勒盟西日期未625內(nèi)錫林郭勒盟問題體中的變化進(jìn)行分析是確定拍攝地點及日期的法。相關(guān)組數(shù)據(jù)皆包含21個等間距時刻（時間）直桿的x,y頂點坐標(biāo)。其地面為x,y平面，直桿垂直于地面。除三組數(shù)據(jù)以外，本題亦給出了一根直桿在下的變化，中20157138540693436秒某地高為2米的直桿的變化過程。以上信息中，各直桿的地理位置皆未知。需要解決的問建立變化的數(shù)學(xué)模型，分析長度關(guān)于各個參數(shù)的變化規(guī)律并應(yīng)20151022日時間9:00-15:00之間廣場（北395426秒，東1162329秒）3米高的直桿的長度建立數(shù)學(xué)模型，根據(jù)題目提供的附件一中的頂點坐標(biāo)數(shù)據(jù)，給出該直桿建立數(shù)學(xué)模型，根據(jù)題目提供的附件二、三種頂點坐標(biāo)數(shù)據(jù)，給出對應(yīng)分析題目提供的，確定拍攝地點的數(shù)學(xué)模型，分別討論在拍攝日期已知和未知兩種情況下，能否確定的拍攝地點與日期，并給出該可假設(shè)二：本文研究的所有對象所處地海拔為0。時刻idl時刻iL直桿在時刻i直桿在時刻i在時刻'在時刻tii時刻iE問題（一陽以定位，結(jié)合實際生驗可知，對于一個自身長度固定的豎直桿，其影問題（二用來表示角度的量應(yīng)當(dāng)能夠體現(xiàn)的方位，參考資料可知，因方位角能夠可以根據(jù)兩個時刻的頂點在同一坐標(biāo)系中的相對位置計算。建立多目標(biāo)規(guī)劃模型，找出可行的直桿對應(yīng)的值。問題（三個目標(biāo)，建立多目標(biāo)規(guī)劃模型，找出可行的直桿對應(yīng)的值與日期。問題（四第四問要求分析題目提供的，確定拍攝地點的數(shù)學(xué)模型，分別討論因第一、二、三問建模過程中，確定并不一定需要直桿長度以及影長使用了關(guān)聯(lián)高度角、日期、緯度、經(jīng)度、及時間的等式，故可利用該等式建附件四處需要對做一定的處理，找出各個時刻直桿尖端部分的變化。圖5.1.1附件4截圖5.1.1為附件4的截圖因已知豎直桿垂直于地面觀察該發(fā)現(xiàn)：此外，由于在所記錄的時間里，直桿轉(zhuǎn)動的角度十分小，而在圖示畫面中，可認(rèn)為所在的直線是和機的視線垂直的，因此可以認(rèn)為長即認(rèn)為和直桿構(gòu)成的平面垂直于機的視線。5.1.24每一幀頂部的像素坐標(biāo)C、直桿頂端的坐標(biāo)A以及上的任一其他5.1.3現(xiàn)，ACB即為高度角。問題（一）的為高度角。參考有關(guān)高度角的資料[2]，可以得到高度角的求解公式如下arcsinsinsincoscoscost其中，為太陽高度角，為太陽赤緯，的計算方式為[3=23.45sin360d L直桿所在地的緯度，t為時角。同時，由于AOBO，故sin ，L,L2 時間，tt t為本地時。本地時的計算方法為ttl，其中，l為目標(biāo)所在位置經(jīng)度， （EoT為：E9.87sin2B7.53cosB1.5sinB，B360d81t15tl9.87sin2B7.53cosB1.5sinB

,

d 經(jīng)過以上分析，可發(fā)現(xiàn)，當(dāng)目標(biāo)、日期、時間都已知時，可以求出相門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的影在給定的定義域中，畫出此直桿長度的變化曲線。16.1.1Ly1Ly

t01, 為豎直桿所在地的緯度，豎直桿在北半球時，為正，否則為負(fù)。t為時角，根6.1.1中的分析，可以得出：t15tl9.87sin2B7.53cosB1.5sinB

,

d l為正，否則為負(fù)。關(guān)于

，定義域為1,7

2 LyLsinsincoscoscost21,t01,

=23.45sin365d81 9.87sin2B7.53cosB1.5sin t15t015 12B

360d圖直桿 圖表示了2015年10月22 廣場（北緯39度54分26秒,東經(jīng)116度23分29秒3米高的直桿 長（米需要注意的是，由于本題中計算赤緯、時角誤差修正的公式并非完全精1158303.841091158直桿的 12:0（395426秒,1162329秒（上圖）395426秒，右（1222日）達(dá)到最長后又開始減小。而對于處在南半球的直桿，情況恰好圖體現(xiàn)的是6月22日時間12:00（即子午線標(biāo)準(zhǔn)時間4:00）北緯第一幅圖過任何處理，為修正前變化曲線?？梢钥闯觯谖鹘?jīng)120度左因此對其討論直桿影長毫無意義。而在此二經(jīng)度附近的區(qū)域，剛剛升起或落下，若當(dāng)?shù)卮蟮剌^為平坦，則長度極大，會影響對其他經(jīng)度區(qū)域長度的116度23分29秒,處于南緯90度到北緯90度的3米高的直桿的長度84度左右90度，時間12:00時為黑夜，因此對其討論直桿影長毫無意義。而在此南緯84度附近的區(qū)域，剛剛升起，若當(dāng)?shù)卮蟮剌^為平坦，則長度極大，會影響對其他經(jīng)度區(qū)域長度的觀察，因此本文將長度大于30米的區(qū)域統(tǒng)一設(shè)定為30米后，繪制出修正后的曲線。中心向東西延伸，長度不斷變長，直到達(dá)到南緯84度左右，。綜上所述，題目要求繪出的廣場3米高豎直桿在10月22日9至表長度關(guān)于各個參數(shù)的變化規(guī)桿長長日逐漸變短，6成側(cè)逐漸變側(cè)逐漸變于第一問的計算求解模型建立過程中，其公式為：高度角的計算公式，即：arcsinsinsincoscoscost。各項中，赤緯以及時角E都是近似式，因而本模型的主要系統(tǒng)誤差來自于此兩項。由于E皆為僅僅與日期d有關(guān)的變量，故可以選擇一個精度較高的天文計算器[4]，將之計算的d1365時的E值作為真實值，而將本文模型中計算的E值作MIDCSPACALCULATOR[5]。對比計算得出的誤6.16.1E1365時采用本文以及真實結(jié)果的比較，下方圖為赤緯的計算值和理論值的比較。其中，藍(lán)色線為本文方案的計算值，棕色線為改進(jìn)后結(jié)果，X點為SPA計算器算出的值，即真實值?？梢钥闯?，

D 1.914180sin0.98565D2 E=720CCA A

C BA1.914sinWD2AWD W 表6.2修正前后時差與赤緯誤差對比最大殘差絕對殘差平平均殘差絕對修正修正修正修正修正修正時赤可見，經(jīng)過精度更高的公式修正之后，時差和赤緯的誤差得到了較為明顯的改善。然而，由于在之后的模型中，時差和赤緯是用來計算高度角表6.3修正前后高度角與方位角誤差對平均殘差絕對值（單位為角度修正修正高度方位圖6.4修正前后高度角與方位角殘差分布圖6.4為修正前后高度角與方位角殘差分布圖。橫軸以小時為單位，從將修正后的公式代入第一問，可以得到影長在時間9至15點的最小問題（二）的本問相對于第一問，不知道直桿的長度，但因知道頂點在續(xù)時間段內(nèi)的坐標(biāo)數(shù)據(jù)，故相當(dāng)于知道了此時間段內(nèi)等間隔時刻的影長以及其角度變化（即方向。因此，合理的（直桿位置）應(yīng)不僅僅能夠擬合直桿的影長及其變化，還應(yīng)當(dāng)擬合方位角的變化。由于特定時刻，特定下的高度角和方位角皆可以計算，且同該示意圖表示的情況已過當(dāng)?shù)貢r間正午，故OB為較早時的成影，OC為較晚時的成影。與第一問相同，B，C分別為OB、OC的高度角。而第一問中，參考有關(guān)高度角的資料[1]，可以得 高度角的求解公式為arcsinsinsincoscoscost，且0,。因此，時刻t,t 2 1 Lcot 長之比即為y1 1=

Lcot y Ly1Ly,sincoscossincostarccos

2 2可以計算出兩時刻方位角理論計算值的差值的絕對值

角差值，因已經(jīng)知道了向量OB,OA的坐標(biāo)

OB。 OB。的。需要注意的是，因為附件中一共提供了21組數(shù)據(jù)，為了使一對數(shù)據(jù)的差別位角差值之差的平方和接近于0,由此建立多目標(biāo)規(guī)劃模型。差盡可能為0二者為目標(biāo)的多目標(biāo)規(guī)劃模型。目標(biāo)函數(shù)為： L'minf,l yi

L'i,j yjminf,l'' i,

j和l分別為直桿所在的緯度和經(jīng)度。目標(biāo)函數(shù)是兩個關(guān)于 數(shù)。其中，f1為每一對理論與實測的影長之比的差的平方和，Lyi,Lyj分別為時 ij測得的影長，L',L'為時刻i,j 值之差的平方和，,分別為時刻i,j測得的方位角，',分別為時刻i,j cotL'yi=cot

2k 2 sincoscossncostk k =23.45sin365d kt15k

9.87sin2B7.53cosB1.5sinB k k B360d 90,

ji條件一給出了理論影長之比的計算方式，k為k時刻的高度角ki,j。條件二給出了k的約束條件，即k須為銳角，僅僅在這種情況下，當(dāng)?shù)貢r間才為白天，討論才有意義。條件三給出了k為赤緯,為豎直桿所在地的緯度，豎直桿在北半球時，為正，否則為負(fù)。tk為時刻k條件五、六、七給出了、tkB的求解方式。其中d11日時，d=1，本問中d108；l為直桿所在經(jīng)度，tk0k的標(biāo)準(zhǔn)時間條件八給出了i,j對應(yīng)所指的時刻。最后兩個條件給出了的定義域。 L'minf,l yi

L'i,j yjminf,l'' i,

j cotL'yi=cot

2k 2 sincoscossncostk k =23.45sin365d kt15k

9.87sin2B7.53cosB1.5sinB k k B360d 90,

ji以minf1,l（ 得到多組l使得f1l，進(jìn)而以第一步的優(yōu)化結(jié)果作為約束條件，即f1l，以minf2l 19109度定位（谷歌地圖析，其他可能點所在的范圍。得到最優(yōu)結(jié)果之后，可以分別以f1l不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖顏色越黃，越淺即代表此的直桿吻合情況越差。f1lf1l19109度以外，還有其他的地點的0.002的點，可以得到七個中心，這些中心也是直桿可能的所直桿可能所在地緯經(jīng)誤差程大致地云-尼西亞，巨越南，河海南西海海南東西沙群其次是在海南東南的北緯18，東經(jīng)109度以及越南河內(nèi)（北緯21度，東經(jīng)106度。但因模型自身或存在系統(tǒng)誤差，測量坐標(biāo)時亦可能出現(xiàn)錯誤，因此，直桿可能所在地在地圖上的位置（谷歌地圖截圖問題（三）的角差值之差盡可能為0二者為目標(biāo)的多目標(biāo)規(guī)劃模型。目標(biāo)函數(shù)為： L'minf,l,d yi

L'i,j yjminf,l,d'' i,

j 和l分別為直桿所在的緯度和經(jīng)度，d11日時，d的值為1。目標(biāo)函數(shù)是兩個關(guān)于與日期的二元函數(shù)。其中，f1為每一對理論LL分別為時刻i,jL',L'為時刻i,j的理論影長。f2為方位角差值之差的平方和，i,j分別為時刻i i',j分別為時刻ij的理論方位角值。cotcotL'y1=cotLy

2k 2 sincoscossincostk k =23.45sin365d kt15k

9.87sin2B7.53cosB1.5sinB k

B

d ji條件一給出了理論影長之比的計算方式，k為k時刻的高度角，k可,為 赤緯,為豎直桿所在地的緯度，豎直桿在北半球時，為正，否則為負(fù)。tk為時刻k的時角。條件四給出了理論方位角的計算公式。條件五、六、七給出了、tk以及一個相關(guān)量B的求解方式。其中d代表天數(shù),l為直桿所在經(jīng)度，tk0為時刻k的 條件八給出了i,j對應(yīng)所指的時刻。最后三個條件分別給出 L'minf,l,d yi

L'i,j yjminf,l,d'' i,

jL'y1=cotLy

2k 2 karcsinsinsincoscoscostk incoscossincostk k s.t.=23.45sin365d t15

9.87sin2B7.53cosB1.5sinB kkk

B

d ji90,90,l與第二問類似，采用分層求解法，利用編程遍歷求解，可得到根第三問最可能答案及在地圖上的大概位置（谷歌地圖緯經(jīng)日大致地附件圖木舒克附件陜西，漢中廣元之得到最優(yōu)結(jié)果之后，與第二問類似，可以分別對特定的緯度與經(jīng)度l對應(yīng)的最小的f1l,d不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖為了更好的描述f1ld的相對大小，可對f1ld統(tǒng)一設(shè)為0.001，由此在保留觀察誤差值趨勢變化的同時，便于讀者觀察。50度，南北回歸線范圍前0.000001的點的誤差值全部規(guī)定為0.000001，修正后得到的圖如下：小于10-9，得出的可能解如下：附件附件經(jīng)度緯度日誤差經(jīng)度緯度日誤差-3.22E--7.23E--6.56E--1.96E--3.94E--7.32E--8.91E--1.19E--3.59E--9.22E--8.85E--2.81E--8.93E--8.39E--6.08E-8.75E--8.77E-8.58E-7.88E-9.94E-5.73E-6.35E-2.83E-6.88E-5.19E-5.93E-9.54E-4.67E-002.11E-1.19E-7.30E-5.94E-6.93E-9.48E-2.97E-問題（四）的信息。由于已知關(guān)系式arcsinsinsincoscoscost，且各個時刻的iarcsinsinsincoscoscosti 其中，=23.45sin360d81d是一個變量。d 當(dāng)日期為11日時，d的值為1。為直桿所在的緯度，為未知數(shù)。t為時角其計算方式為：t15

l

9.87sin2B7.53cosB1.5sin

12 k 為利用第i因而可以看出，當(dāng)日期已知時，待求解方程組中有與l知時，有與l和d

iarcsinsinsincoscoscosti

=23.45sin360d 9.87sin2B7.53cosB1.5sin tk15tk015 12B

360d利用編程，考慮到天文學(xué)公式自身存在系統(tǒng)誤差，故轉(zhuǎn)換為非線在日期已知的情況下，以minsinsincoscoscostsin2 i目標(biāo)函數(shù)，在日期未知的情況 ， i期確定與未確定的情況下，計算出的相應(yīng)、日期和大致地點如下：第四問最可能答案及在地圖上的大概位置（谷歌地圖）緯經(jīng)日大致地日期確內(nèi)錫林郭勒盟西日期未內(nèi)錫林郭勒盟不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖110150度左右，誤差明顯減0.0010.001，修1104040度左右。因本題所采用的模型自身圖日期確定時修正后不同的理論影長與實際測得影長的誤差值圖給出了日期確定的情況下不同的理論影長與實際測得的影對日期不定的情況下取誤差值低于10-3日期確日期未經(jīng)緯誤差經(jīng)緯日誤差0-----3.73E--------------0上述的點為第四問期確認(rèn)時與日期不確認(rèn)時直桿可能的所在地及對應(yīng)模型于一般的。用于分析同一組拍攝的地點及時間檢驗或是否被為品PVEducation,AzimuthAngle,of-sunlight/azimuth-angle,2015/9/12PVEducation,ElevationAngle,of-sunlight/elevation-angle,2015/9/12 ,2015/9/12RedaI,AndreasA.SolarPositionAlgorithmforSolarRadiationApplications. , , /public3/equatime.html,September13,2015 ,September13,附錄一：所用R2015aAdobeAfterEffectAdobePhotoshop附錄二 源程序代第一問解答代碼（Q1：clearall;clc;L=3;LT=9:0.001:15;DateString='22-Oct-2015';DateString2='1-Jan-2015';formatIn='dd-mmm-yyyy';d=datenum(DateString,formatIn)-datenum(DateString2,formatIn)+1;dGMT=LT-GMT;LSTM=B=(360./365).*(d-EoT=9.87.*sin(2.*B./180.*pi)-7.53.*cos(B./180.*pi)-1.5.*sin(B./180.*pi);TC=4.*(LONGITUDE-LSTM)+EoT;LST=LT+TC./60;HRA=15.*(LST-12);DEC=elevation=asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi)+azimuth=acos((sin(DEC./180.*pi).*cos(LATITUDE./180.*pi)-%azimuth=%elevation=forifLST(i)>12||azimuth(i)=2*pi-第二問解答程序data=xlsread('dataQ2','sheet1'len=sqrt(data(:,1).^2+data(:,2).^2);left=1:ceil(length(len)/2);right=ceil(length(len)/2):(length(len)%left=1:length(len)-%right=2:length(len);bleft=left;bright=rateT=len(left)./len(right)=data(bleft,=data(bright,=data(bleft,=data(bright,angleT=acos((Xl.*Xr+Yl.*Yr)./(len(bleft).*len(bright)));longRange=-180:180;latiRange=-90:90ErrorT=zeros(length(longRange),ErrorTA=zeros(length(longRange),DateString='11-Sep-2015';DateString2='1-Jan-2015';formatIn='dd-mmm-yyyy';d=datenum(DateString,formatIn)-datenum(DateString2,formatIn)+1;L=LT=14.7%LT=GMT=LT-dGMT=LT-GMT;LSTM=B=(360./365).*(d-EoT=9.87.*sin(2.*B./180.*pi)-7.53.*cos(B./180.*pi)-1.5.*sin(B./180.*pi);fori=1:forj=1:length(latiRange)LONGITUDE=longRange(i);LATITUDE=latiRange(j);TC=4.*(LONGITUDE-LSTM)+EoT;LST=LT+TC./60;HRA=15.*(LST-DEC=elevation=asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi)+azimuth=acos(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi)-azimuth(elevation<0)=angleP=abs(azimuth(bleft)-%azimuth=%elevation=

%azimuth(LST>12)=2*pi-azimuth(LSTlenP=L./tan(elevation);lenP(elevation<0)=0;rateP=lenP(left)./lenP(right);ErrorT(i,j)=sqrt(sum((rateT-rateP').^ErrorTA(i,j)=sqrt(sum((angleP-angleT').^2));Tmin=min(min(ErrorT));tol=1.01*Tmin;%ErrorT((ErrorT>tol))=tol;ErrorT((ErrorT>1))=1;ErrorTA((ErrorTA>1))=1;ErrorTA((ErrorT==tol))=1;%ErrorT=(ErrorT-min(min(ErrorT)))/(max(max(ErrorT))-min(min(ErrorT)));%ErrorTA=(ErrorTA-min(min(ErrorTA)))/(max(max(ErrorTA))-min(min(ErrorTA)));%ErrorTT=ErrorT+[X,Y]=meshgrid(longRange,%surf(X,Y,ErrorTT');contour(X,Y,ErrorT(selectX,%surf(X,Y,%contour(X,Y,%[a,b]=find(ErrorTT第三問解答程序data=xlsread('dataQ2','sheet2');xTemp=data(:,1);yTemp=data(:,T=(12+41/60):1/20:(13+41/60len=sqrt(xTemp.^2+yTemp.^%rateT=len(1:ceil(length(len)/2))./len(ceil(length(len)/2):length(len));rateT=len(1:(length(len)-1))./len(2:length(len));%angleT=acos((data(1:ceil(length(len)/2).*data(ceil(length(len)/2):(length(len)),1)+...%data(1:ceil(length(len)/2),2).*data(ceil(length(len)/2):(length(len)),)./%(len(1:ceil(length(len)/2)).*len(ceil(length(len)/2):%length(len)))%angleT=acos((data(1:length(len)-1,1).*data(2:length(len),1)+...%(data(1:length(len)-1,2).*(length(len)),2)))%(len(1:length(len)-1).*len(2:length(len))));angleT=acos((xTemp(1:length(len)-1).*xTemp(2:length(len))+...(yTemp(1:length(len)-1).*(length(len)))))(len(1:length(len)-1).*len(2:length(len))));%angleT=atan(yTemp(2:length(len))./xTemp(2:length(len)))-...%atan(yTemp(1:length(len)-1)./xTemp(1:length(len)-1));longRange=-180:latiRange=-90:90dateRange=1:365ErrorT=zeros(length(longRange),length(latiRange),length(dateRange));%ErrorTA=zeros(length(longRange),fori=1:length(longRange)forj=1:length(fork=1:length(d=kLONGITUDE=longRange(i);LATITUDE=latiRange(j);L=3;LT=GMT=LT-dGMT=LT-GMT;LSTM=B=(360./365).*(d-EoT=9.87.*sin(2.*B./180.*pi)-7.53.*cos(B./180.*pi)-1.5.*sin(B./180.*pi);TC=4.*(LONGITUDE-LSTM)+EoT;LST=LT+TC./60;HRA=15.*(LST-DEC=elevation=asin(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi)+%azimuth=acos(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi)-%azimuth(elevation<0)=%angleP=abs(azimuth(1:ceil(length(len)/2))-azimuth(ceil(length(len)/2):length(len)));%angleP=abs((azimuth(1:length(len)1)-azimuth(2:%azimuth=%elevation=%azimuth(LST>12)=2*pi-azimuth(LST

lenP=L./tan(elevation);lenP(elevation<0)=0;%rateP=lenP(1:ceil(length(len)/2))lenP(ceil(length(len)/2):length(len));rateP=lenP(1:length(len)-1)lenP(2:length(len)ErrorT(i,j,k)=(sum((rateP').^%ErrorTA(i,j)=angleT').^2));ErrorT((ErrorT>1))=1%ErrorTA((ErrorTA>0.5))=0.5%ErrorT=(ErrorT-min(min(ErrorT)))/(max(max(ErrorT))-min(min(ErrorT)));%ErrorTA=(ErrorTA-min(min(ErrorTA)))/(max(max(ErrorTA))-min(min(ErrorTA)));%ErrorTT=ErrorT+RangeX=1:361RangeY=1:[X,Y]=meshgrid(longRange(RangeX),min(ErrorT(),[],[temp,test]=min(ErrorT(),[],temp=(temp-min(min(temp)))/-%surf(X,Y,A=temp(RangeX,RangeY);A(A>0.000001)=0.000001;%contour(X,Y,%surf(X,Y,%surf(X,Y,%surf(X,Y,ErrorTA');acc=1*10^(-6);[a,b]=find(A<acc);a=a-181b=b-c=[abtest(find(A<acc))A(find(A<acc))];第四問解答程序（Q4,日期已知clc;dataraw=xlsread('keyframeData','sheet1');top=[891.2205];bottom=[891.2L=bottom(2)-top(2)len=sqrt((dataraw(:,1)-top(1)).^2(dataraw(:,2)-top(2)).^2);startTime=8+54/60+7/3600;endTime=9+34/60+46/3600;step=(endTime-startTime)/(length(dataraw)-elevationTest=L./lenT=startTime:step:endTime;longRange=-180:180;latiRange=-90:90ErrorT=zeros(length(longRange),%ErrorTA=zeros(length(longRange),DateString='13-July-2015';DateString2='1-Jan-2015';formatIn='dd-mmm-yyyy';d=datenum(DateString,formatIn)-datenum(DateString2,formatIn)+1;fori=1:length(longRange)forj=1:length(latiRange)LONGITUDE=longRange(i);LATITUDE=latiRange(j);LT=T;GMT=LT-dGMT=LT-GMT;LSTM=B=(360./365).*(d-EoT=9.87.*sin(2.*B./180.*pi)-7.53.*cos(B./180.*pi)-1.5.*sin(B./180.*pi);TC=4.*(LONGITUDE-LSTM)+EoT;LST=LT+TC./60;HRA=15.*(LST-DEC=elevation=(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi)+%azimuth=acos(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi)-%azimuth(elevation<0)=%angleP=abs(azimuth(1:ceil(length(len)/2))-azimuth(ceil(length(len)/2):length(len)));%angleP=abs(azimuth(1:length(len)-1-azimuth(2:%azimuth=%elevation=%azimuth(LST>12)=2*pi-azimuth(LST%lenP=L%lenP(elevation<0)=%rateP=lenP(1:ceil(length(len)/2))lenP(ceil(length(len)/2):length(len))%rateP=lenP(1:length(len)-1)./lenP(2:length(len));ErrorT(i,j)=sqrt(sum((elevation-elevationTest').^2));%ErrorTA(i,j)=sqrt(sum((angleP-angleT').^2));ErrorT((ErrorT>1))=1%ErrorTA((ErrorTA>0.5))=0.5ErrorT=(ErrorT-min(min(ErrorT)))/(max(max(ErrorT))-min(min(ErrorT)));%ErrorTA=(ErrorTA-min(min(ErrorTA)))/(max(max(ErrorTA))-min(min(ErrorTA)));%ErrorTT=ErrorT+[X,Y]=meshgrid(longRange,%surf(X,Y,%surf(X,Y,%surf(X,Y,ErrorTA');count=0;fori=2:length(longRange)-forj=2:length(latiRange)-if(ErrorT(i,j)-ErrorT(i,j+1)<&&ErrorT(i,j)-ErrorT(i,-1)<0ErrorT(i,j)-ErrorT(i+1j)<0&&ErrorT(i,-ErrorT(i-1j)<0

count=count+1c(:,count)=[i-181j-91ErrorT(i

contour(X,Y,%[a,b]=find(ErrorTT第四問解答程序（Q4_2,日期未知clc;dataraw=xlsread('keyframeData','sheet1');top=[891.2205];bottom=[891.2L=bottom(2)-top(2)len=sqrt((dataraw(:,1)-top(1)).^2(dataraw(:,2)-top(2)).^2);startTime=8+54/60+7/3600;endTime=9+34/60+46/3600;step=(endTime-startTime)/(length(dataraw)-elevationTest=L./lenT=startTime:step:endTime;longRange=-180:180;latiRange=-90:90dateRange=1:ErrorT=zeros(length(longRange),length(latiRange),length(dateRange));%ErrorTA=zeros(length(longRange),fori=1:length(longRange)forj=1:length(fork=1:length(dateRange)d=k;LONGITUDE=longRange(i);LATITUDE=latiRange(j);LT=T;GMT=LT-dGMT=LT-GMT;LSTM=B=(360./365).*(d-EoT=9.87.*sin(2.*B./180.*pi)-7.53.*cos(B./180.*pi)-1.5.*sin(B./180.*pi);TC=4.*(LONGITUDE-LSTM)+EoT;LST=LT+TC./60;HRA=15.*(LST-DEC=elevation=(sin(DEC./180.*pi).*sin(LATITUDE./180.*pi)+%azimuth=acos(sin(DEC./180.*pi).*cos(LATITUDE./180.*pi)-%azimuth(elevation<0)=%angleP=abs(azimuth(1:ceil(length(len)/2))-azimuth(ceil(length(len)/2):length(len)));%angleP=abs(azimuth(1:length(len)-1-azimuth(2:%azimuth=%elevation=%azimuth(LST>12)=2*pi-azimuth(LST

%lenP=L%lenP(elevation<0)=%rateP=lenP(1:ceil(length(len)/2))lenP(ceil(length(len)/2):length(len))%rateP=lenP(1:length(len)-1)./lenP(2:length(len));ErrorT(i,j,k)=sqrt(sum((elevation-elevationTest').^2));%ErrorTA(i,j)=sqrt(sum((angleP-angleT').^2));ErrorT((ErrorT>30))=30RangeX=1:361RangeY=1:[X,Y]=meshgrid(longRange(RangeX),[temp,test]=min(ErrorT(),[],temp=(temp-min(min(temp)))/-A=temp(RangeX,RangeY);A(A>0.1)=0.1;surf(X,Y,acc=1*10^(-6);[a,b]=find(A<acc);a=a-181b=b-c=[abtest(find(A<acc))A(find(A<acc))];赤緯誤差分析d=DEC2=-asin(0.39779*cosd(0.98565*(d+10)1.91