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

,

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

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

,

d l為正，否則為負。關于

，定義域為1,7

2 LyLsinsincoscoscost21,t01,

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

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

D 1.914180sin0.98565D2 E=720CCA A

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

Lcot y Ly1Ly,sincoscossincostarccos

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

角差值，因已經知道了向量OB,OA的坐標

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

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

j和l分別為直桿所在的緯度和經度。目標函數是兩個關于 數。其中，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須為銳角，僅僅在這種情況下，當地時間才為白天，討論才有意義。條件三給出了k為赤緯,為豎直桿所在地的緯度，豎直桿在北半球時，為正，否則為負。tk為時刻k條件五、六、七給出了、tkB的求解方式。其中d11日時，d=1，本問中d108；l為直桿所在經度，tk0k的標準時間條件八給出了i,j對應所指的時刻。最后兩個條件給出了的定義域。 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，進而以第一步的優化結果作為約束條件，即f1l，以minf2l 19109度定位（谷歌地圖析，其他可能點所在的范圍。得到最優結果之后，可以分別以f1l不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖顏色越黃，越淺即代表此的直桿吻合情況越差。f1lf1l19109度以外，還有其他的地點的0.002的點，可以得到七個中心，這些中心也是直桿可能的所直桿可能所在地緯經誤差程大致地云-尼西亞，巨越南，河海南西海海南東西沙群其次是在海南東南的北緯18，東經109度以及越南河內（北緯21度，東經106度。但因模型自身或存在系統誤差，測量坐標時亦可能出現錯誤，因此，直桿可能所在地在地圖上的位置（谷歌地圖截圖問題（三）的角差值之差盡可能為0二者為目標的多目標規劃模型。目標函數為： L'minf,l,d yi

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

j 和l分別為直桿所在的緯度和經度，d11日時，d的值為1。目標函數是兩個關于與日期的二元函數。其中，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可,為 赤緯,為豎直桿所在地的緯度，豎直桿在北半球時，為正，否則為負。tk為時刻k的時角。條件四給出了理論方位角的計算公式。條件五、六、七給出了、tk以及一個相關量B的求解方式。其中d代表天數,l為直桿所在經度，tk0為時刻k的 條件八給出了i,j對應所指的時刻。最后三個條件分別給出 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與第二問類似，采用分層求解法，利用編程遍歷求解，可得到根第三問最可能答案及在地圖上的大概位置（谷歌地圖緯經日大致地附件圖木舒克附件陜西，漢中廣元之得到最優結果之后，與第二問類似，可以分別對特定的緯度與經度l對應的最小的f1l,d不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖為了更好的描述f1ld的相對大小，可對f1ld統一設為0.001，由此在保留觀察誤差值趨勢變化的同時，便于讀者觀察。50度，南北回歸線范圍前0.000001的點的誤差值全部規定為0.000001，修正后得到的圖如下：小于10-9，得出的可能解如下：附件附件經度緯度日誤差經度緯度日誤差-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-問題（四）的信息。由于已知關系式arcsinsinsincoscoscost，且各個時刻的iarcsinsinsincoscoscosti 其中，=23.45sin360d81d是一個變量。d 當日期為11日時，d的值為1。為直桿所在的緯度，為未知數。t為時角其計算方式為：t15

l

9.87sin2B7.53cosB1.5sin

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

iarcsinsinsincoscoscosti

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

360d利用編程，考慮到天文學公式自身存在系統誤差，故轉換為非線在日期已知的情況下，以minsinsincoscoscostsin2 i目標函數，在日期未知的情況 ， i期確定與未確定的情況下，計算出的相應、日期和大致地點如下：第四問最可能答案及在地圖上的大概位置（谷歌地圖）緯經日大致地日期確內錫林郭勒盟西日期未內錫林郭勒盟不同反映的是不同的理論影長與實際測得的影長的吻合情況，圖110150度左右，誤差明顯減0.0010.001，修1104040度左右。因本題所采用的模型自身圖日期確定時修正后不同的理論影長與實際測得影長的誤差值圖給出了日期確定的情況下不同的理論影長與實際測得的影對日期不定的情況下取誤差值低于10-3日期確日期未經緯誤差經緯日誤差0-----3.73E--------------0上述的點為第四問期確認時與日期不確認時直桿可能的所在地及對應模型于一般的。用于分析同一組拍攝的地點及時間檢驗或是否被為品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