課件文稿成果

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

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

,

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

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

D 1.914180sin0.98565D2 E=720CCA A

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

Lcot y Ly1Ly,sincoscossincostarccos

2 2可以計(jì)算出兩時(shí)刻方位角理論計(jì)算值的差值的絕對(duì)值

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

OB。 OB。的。需要注意的是，因?yàn)楦郊幸还蔡峁┝?1組數(shù)據(jù)，為了使一對(duì)數(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ù)是兩個(gè)關(guān)于 數(shù)。其中，f1為每一對(duì)理論與實(shí)測的影長之比的差的平方和，Lyi,Lyj分別為時(shí) ij測得的影長，L',L'為時(shí)刻i,j 值之差的平方和，,分別為時(shí)刻i,j測得的方位角，',分別為時(shí)刻i,j cotL'yi=cot

2k 2 sincoscossncostk k =23.45sin365d kt15k

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

ji條件一給出了理論影長之比的計(jì)算方式，k為k時(shí)刻的高度角ki,j。條件二給出了k的約束條件，即k須為銳角，僅僅在這種情況下，當(dāng)?shù)貢r(shí)間才為白天，討論才有意義。條件三給出了k為赤緯,為豎直桿所在地的緯度，豎直桿在北半球時(shí)，為正，否則為負(fù)。tk為時(shí)刻k條件五、六、七給出了、tkB的求解方式。其中d11日時(shí)，d=1，本問中d108；l為直桿所在經(jīng)度，tk0k的標(biāo)準(zhǔn)時(shí)間條件八給出了i,j對(duì)應(yīng)所指的時(shí)刻。最后兩個(gè)條件給出了的定義域。 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度定位（谷歌地圖析，其他可能點(diǎn)所在的范圍。得到最優(yōu)結(jié)果之后，可以分別以f1l不同反映的是不同的理論影長與實(shí)際測得的影長的吻合情況，圖顏色越黃，越淺即代表此的直桿吻合情況越差。f1lf1l19109度以外，還有其他的地點(diǎn)的0.002的點(diǎn)，可以得到七個(gè)中心，這些中心也是直桿可能的所直桿可能所在地緯經(jīng)誤差程大致地云-尼西亞，巨越南，河海南西海海南東西沙群其次是在海南東南的北緯18，東經(jīng)109度以及越南河內(nèi)（北緯21度，東經(jīng)106度。但因模型自身或存在系統(tǒng)誤差，測量坐標(biāo)時(shí)亦可能出現(xiàn)錯(cuò)誤，因此，直桿可能所在地在地圖上的位置（谷歌地圖截圖問題（三）的角差值之差盡可能為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日時(shí)，d的值為1。目標(biāo)函數(shù)是兩個(gè)關(guān)于與日期的二元函數(shù)。其中，f1為每一對(duì)理論LL分別為時(shí)刻i,jL',L'為時(shí)刻i,j的理論影長。f2為方位角差值之差的平方和，i,j分別為時(shí)刻i i',j分別為時(shí)刻ij的理論方位角值。cotcotL'y1=cotLy

2k 2 sincoscossincostk k =23.45sin365d kt15k

9.87sin2B7.53cosB1.5sinB k

B

d ji條件一給出了理論影長之比的計(jì)算方式，k為k時(shí)刻的高度角，k可,為 赤緯,為豎直桿所在地的緯度，豎直桿在北半球時(shí)，為正，否則為負(fù)。tk為時(shí)刻k的時(shí)角。條件四給出了理論方位角的計(jì)算公式。條件五、六、七給出了、tk以及一個(gè)相關(guān)量B的求解方式。其中d代表天數(shù),l為直桿所在經(jīng)度，tk0為時(shí)刻k的 條件八給出了i,j對(duì)應(yīng)所指的時(shí)刻。最后三個(gè)條件分別給出 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é)果之后，與第二問類似，可以分別對(duì)特定的緯度與經(jīng)度l對(duì)應(yīng)的最小的f1l,d不同反映的是不同的理論影長與實(shí)際測得的影長的吻合情況，圖為了更好的描述f1ld的相對(duì)大小，可對(duì)f1ld統(tǒng)一設(shè)為0.001，由此在保留觀察誤差值趨勢變化的同時(shí)，便于讀者觀察。50度，南北回歸線范圍前0.000001的點(diǎn)的誤差值全部規(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，且各個(gè)時(shí)刻的iarcsinsinsincoscoscosti 其中，=23.45sin360d81d是一個(gè)變量。d 當(dāng)日期為11日時(shí)，d的值為1。為直桿所在的緯度，為未知數(shù)。t為時(shí)角其計(jì)算方式為：t15

l

9.87sin2B7.53cosB1.5sin

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

iarcsinsinsincoscoscosti

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

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