地理信息系统4班蔡武阳2010301130116在后面的计算过程中，我们.DOC

地理信息系统4班 蔡武阳 2010301130116 在后面的计算过程中，我们选择的是克拉索斯基椭球体，即 a=6378245，b=6356863? 【1】计算地球子午圈曲率半径、卯酉圈曲率半径和纬圈半径； 其中， 子午圈曲率半径：M= 卯酉圈曲率半径：N= 纬圈半径：r=Ncos= 第一偏心率：e= 下面是具体计算过程： 在Matlab中，取纬度变量为从0到，变化间隔为 Matlab代码实现： a=6378245; b=6356863; e=sqrt((a*a-b*b)/(a*a)); x=0:pi/180:pi/2; M=a*(1-e^2)./((1-e.^2*(sin(x).^2)).^1.5); N=a./((1-e.^2*(sin(x).^2)).^0.5); r=N.*cos(x); M=M; N=N; r=r; xlswrite(1,M,B2:B92) xlswrite(1,N,C2:C92) xlswrite(1,r,D2:D92) 【2】计算赤道到纬度为的子午线弧长； 由赤道至弧度为纬度间的子午线弧长为：s= 积分后经整理得子午线弧长的一般公式: s=a（1-e）； A=1.0050517739 B=0.00506237764 C=0.0000106245 D=0.00000002081 由于D之后的算式系数几乎为0，故只算前四项即可为从0到，变化间隔为 Matlab代码实现： a=6378245; b=6356863; e=sqrt((a*a-b*b)/(a*a)); y=0:1:60; A=1.0050517739; B=0.00506237764; C=0.0000106245; D=0.00000002081; p=180/pi; s=a*(1-e^2)*((A/p) *y-(B/2)*sind(2*y)+(C/4)*sind(4*y)-(D/6)*sind(6*y)); xlswrite(1,s, E2:E62); 【3】计算赤道到纬度为之间的1弧度经差纬差的球面梯形面积。 球面梯形面积T的计算公式为： T=(1-e) 积分并整理后得： T=K(Asin(7)+……) 式中： K=2a A=1+++ B= C= D= 令=1弧度，；；； K=2a 由于e较小，取算式前4项计算即可. 下面是具体计算过程： Matlab程序实现： a=6378245; b=6356863; e=sqrt((a*a-b*b)/(a*a)); A=1+(1/2)*e^2+(3/8)*e^4+(5/16)*e^6; B=(1/6)*e^2+(3/16)*e^4+(3/16)*e^6; C=(3/80)*e^4+(1/16)*e^6; D=(1/112)*e^6; K=2*a*a*(1-e*e); z=0:pi/180:pi/3; T=K.*(A*sin(z/2).*cos(z/2)-B*sin(3*z/2).*cos(3*z/2)+C*sin(5*z/2).*cos(5*z/2)-D*sin(7*z/2).*cos(7*z/2)); xlswrite(1,T,F2:F62); 附录： 由Matlab生成的Excel表格如下 纬度ψ（度） 子午圈曲率半径M（米） 卯酉圈曲率半径N（米） 纬圈半径r（米） 子午线弧长s（米） 球面梯形面积T（平方公里） 0 6335553 6378245 6378245 0 0 1 6335572 6378252 6377280 110576 705248 2 6335630 6378271 6374386 221153 1410286 3 6335727 6378303 6369562 331732 2114907 4 6335862 6378349 6362812 442312 2818900 5 6336036 6378407 6354135 552895 3522057 6 6336248 6378478 6343536 663482 4224171 7 6336498 6378562 6331017 774072 4925031 8 6336785 6378658 6316582 884668 5624430 9 6337110 6378767 6300234 995268 6322161 10 6337471 6378889 6281979 1105875 7018016 11 6337869 6379022 6261822 1216488 77117