大地測量學(xué)實(shí)習(xí)報(bào)告.doc

課程編號：1300158 課程性質(zhì)：必修大地測量計(jì)算與實(shí)習(xí)實(shí) 習(xí) 報(bào) 告學(xué) 院： 測繪學(xué)院 專業(yè)方向： 測繪工程 A 方向 實(shí)習(xí)地點(diǎn)： 武漢大學(xué) 班級組號： 3 班第 1 小組 學(xué)生姓名： 田寶林 學(xué)號： 2010301610124 指導(dǎo)教師： 劉宗泉 蘇新洲 丁士俊向 東 黃海蘭 2012年 6月18日 至 2012年 7月6日 目 錄 1.概況 -4-2.已有資料情況-6-3測區(qū)自然地理況-6-4引用文件 -6-5.主要成果產(chǎn)品技術(shù)指標(biāo)-6-5.1關(guān)于控制測量及導(dǎo)線測量-6-5.2關(guān)于碎步測量 -7-5.3關(guān)于繪圖要求-8-6.具體設(shè)計(jì)方案-8-6.1圖根控制測量-8- 6.2碎步測量：-9-6.3內(nèi)業(yè)成圖-9-7. 設(shè)計(jì)評定-9- 1 實(shí)習(xí)目的與意義 本次實(shí)習(xí)的是在我們完成大地測量學(xué)基礎(chǔ)這門課程之后進(jìn)行的，通過完成水準(zhǔn)儀有關(guān)一起的檢驗(yàn)和二等水準(zhǔn)測量，使我們熟悉儀器的操作，并在實(shí)習(xí)過程中鍛煉我們的實(shí)際動手能力，提升團(tuán)隊(duì)協(xié)作能力以及鞏固我們在課堂所學(xué)的理論知識。另外，在后期的編程和外業(yè)概算過程中，對我們的發(fā)現(xiàn)問題、提出問題、解決問題的能力得到一次提升，為今后進(jìn)入社會打下良好基礎(chǔ)。2 實(shí)習(xí)任務(wù) 本次實(shí)習(xí)的任務(wù)有兩項(xiàng)，分別是：1）二等精密水準(zhǔn)測量外業(yè)觀測與概算；（約1.5周）2）大地測量計(jì)算課程設(shè)計(jì)；（約1.5周）三測區(qū)概況 本次實(shí)習(xí)的地點(diǎn)為武漢大學(xué)。校內(nèi)人員眾多，交通復(fù)雜，地勢起伏。我組測量路線為武測環(huán)和珞珈山環(huán)。其具體路線如下圖： 此圖為武測環(huán)，上為北方向 此圖為珞珈山環(huán)，上為北方向四已知高程點(diǎn)共有三個(gè)已知點(diǎn)可選用，我組所用點(diǎn)為珞珈山環(huán)的教務(wù)部點(diǎn)，已知高程為126.157m。五作業(yè)依據(jù) 國家測繪局，國家一、二等水準(zhǔn)測量規(guī)范2006-05-24 測繪出版社，2010儀器的檢驗(yàn) (1) 水準(zhǔn)儀的i角限差為15(2)標(biāo)尺的零點(diǎn)不等差為0.10mm儀器類型光學(xué)水準(zhǔn)儀數(shù)字水準(zhǔn)儀備注前后視距差1.0m1.5m任一測站前后視距累積差3.0m6.0m測站觀測順序和方法奇數(shù)站后前前后奇數(shù)站后前前后左邊表格均為往測，光學(xué)儀器返測時(shí)與往測相反，數(shù)字水準(zhǔn)儀返測與往測相同偶數(shù)站前后后前偶數(shù)站前后后前閉合差往返高差不符值 4mm ( k為測段長度，單位為km,0.1km按0,1km算)環(huán)閉合差 4mm ( F為環(huán)線長度，單位為km )六踏勘、選點(diǎn)本次實(shí)習(xí)的路線圖已經(jīng)提前下發(fā)給我們，所以選點(diǎn)比較固定。對于武測環(huán)，在武大信息學(xué)部部分，我們選擇了二食堂門口、星湖園、大學(xué)生活動中心三個(gè)點(diǎn)；在武大文理學(xué)部部分，我們選擇了明賢門門內(nèi)門外兩點(diǎn)、校醫(yī)院點(diǎn)、生科院點(diǎn)，外加教務(wù)部已知點(diǎn)點(diǎn)，共八個(gè)點(diǎn)。對于珞珈山環(huán)，我們選擇了生科院點(diǎn)、校醫(yī)院點(diǎn)、政管院點(diǎn)、楓園點(diǎn)等共七個(gè)點(diǎn)，加上教務(wù)部已知點(diǎn)共八個(gè)點(diǎn)。7 使用的儀器與儀器檢驗(yàn)我組使用的儀器是科利達(dá)的電子水準(zhǔn)儀DL2007，用于二等水準(zhǔn)測量。我們進(jìn)行了水準(zhǔn)儀i角檢驗(yàn)和水準(zhǔn)尺零點(diǎn)差檢驗(yàn)（結(jié)果另附），符合測量規(guī)范要求。八.精密二等水準(zhǔn)數(shù)據(jù)采集和外業(yè)數(shù)據(jù)概算：水準(zhǔn)線路圖見附錄，觀測日期與觀測時(shí)段在觀測記錄薄中記載詳細(xì)。9 外業(yè)概算成果概算成果見附錄。10 課程設(shè)計(jì)編程 1.編程所用語言本次編程用C+語言在VC6.0環(huán)境下編制 2.基本數(shù)學(xué)模型 （1）高斯投影正反算 正算是指：由大地坐標(biāo)（L, B）求得高斯平面坐標(biāo)（x, y）的過程。反算是指：由高斯平面坐標(biāo)（x, y）求得大地坐標(biāo)（L, B）的過程。正算：高斯投影必須滿足的三個(gè)條件：（1），中央子午線投影后為直線。（2），中央子午線投影后長度不變。（3），投影具有正性性質(zhì)，即正性投影條件。由第一個(gè)條件可知，中央子午線東西兩側(cè)的投影必然對稱于中央子午線。設(shè)在托球面上有P1 ，P2，且對稱于中央子午線。其大地坐標(biāo)為（l, B）,(-l, B)則投影后的平面坐標(biāo)一定為P1（x, y）,P2（x,-y）.由第二個(gè)條件可知，位于中央子午線上的點(diǎn)，投影后的縱坐標(biāo)x應(yīng)該等于投影前從赤道量至該點(diǎn)的子午弧長。計(jì)算公式：1.當(dāng)將克拉索夫斯基橢球帶入計(jì)算式，可得到正算公式：其中： 2.反算公式為：其中： （2）實(shí)測斜距歸算高斯平面邊長 假設(shè)1、2兩個(gè)大地點(diǎn)在橢球面上沿法線的投影點(diǎn)1和2間的大地線的長度為S，由于在橢球面上兩點(diǎn)間大地線長度與相應(yīng)法截線長度之差是極微小的，可以忽略不計(jì)，則可以將兩點(diǎn)間的法截線長度認(rèn)為是該兩點(diǎn)間的大地線長度。并且，兩點(diǎn)間的法截線的長度與半徑等于其起始點(diǎn)曲率半徑的圓弧長相差也很小，則所求的大地線長度可以認(rèn)為是半徑。其計(jì)算如下：S=D*1-(h2-h1)/D*(h2-h1)/D/(1+h1/Ra)*(1+h2/Ra)這個(gè)題目的思想是先利用題目所給的坐標(biāo)求出其平面坐標(biāo)方位角，然后計(jì)算子午線收斂角和方向改化。得出大地方位角，然后將實(shí)測距離歸算至橢球面上，最后歸算至高斯平面，具體流程圖如下。斜距D12，大地方位角A12X1和Y1、X2和Y2實(shí)測斜距D12歸算至大地線長度S12大地線S12化算至高斯平面邊長S （3）大地主題正反算 大地主題解算：知道某些大地元素推求另一些大地元素的過程。 正解是指：已知某點(diǎn)P1的大地坐標(biāo)（L2,B2），且知該點(diǎn)到另一點(diǎn)P2(L2,B2)的大地線長及其大地方位角A12，計(jì)算P2點(diǎn)的大地坐標(biāo)（L2,B2）和大地線在P2點(diǎn)的反方位角A21.的過程。 反解是指：已知P1和P2的大地坐標(biāo)（L1,B1）和P2（L2,B2）計(jì)算P1至P2的大地線長，正反方位角A12、A21的過程。 大地主題解算的基本思想：運(yùn)用高斯平均引數(shù)的方法，。正算基本思想： 在大地線中點(diǎn) M展開，收斂快，精度高； 中點(diǎn)M不好求，以兩端點(diǎn)平均緯度及平均方位角相對應(yīng)的點(diǎn)m來代替； 借助迭代法實(shí)現(xiàn)。反算基本思想3. 各個(gè)程序主要圖框與結(jié)果 （1）高斯投影正反算高斯投影正反算（采用第一組三號數(shù)據(jù)為解算數(shù)據(jù)）已知數(shù)據(jù)X=4574846.0098Y=437239.8566反算B=411819.575651L=125142.163595正算4574846.0105437239.8566（2） 距離歸化 實(shí)測斜距化算至高斯投影平面邊長（采用克拉索夫斯基橢球）：S= 578.868606;(3) 大地主題正反算大地主題解算（采用3班第1組第3序號數(shù)據(jù)）正算反算已知數(shù)據(jù)計(jì)算結(jié)果已知數(shù)據(jù)計(jì)算結(jié)果B1=41.03356874B2=41.58569529B1=41.03356874S=83999.9995L1=130.10122676L2=130.12086865L1=130.10122676A12=1.494227A12=1.4943A21=181.5103B2=41.58569529A21=181.505957S12=84000L2=130.12086865附錄5高斯投影主要代碼：void CMyDlg:OnButton1() double PI=3.14159265358979323846264338327950288; double Bhudu,L,L0,l,N,a0,a4,a5,a6,a3,x,y; UpdateData(); Bhudu=(m_du2*3600+m_fen2*60+m_miao2)*2*PI/360/3600;/化弧度 L=m_du1*3600+m_fen1*60+m_miao1;/化秒 L0=126*3600; l=(L-L0)*2*PI/360/3600;/化弧度N=6399698.902-(21562.267-(108.973-0.612*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu);a0=32140.404-(135.3302-(0.7092-0.0040*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu);a4=(0.25+0.00252*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)-0.04166;a6=(0.166*cos(Bhudu)*cos(Bhudu)-0.084)*cos(Bhudu)*cos(Bhudu);a3=(0.3333333+0.001123*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.0040*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu);m_x1=6367558.4969*Bhudu-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(Bhudu)*cos(Bhudu);m_y1=(1+(a3+a5*l*l)*l*l)*l*N*cos(Bhudu)+500000;UpdateData(FALSE);void CMyDlg:OnButton2() double PI=3.14159265358979323846264338327950288;double x,y,Bf,Be,l,L,Z,Nf,b2,b3,b4,b5,p,B;int B1,B2,L1,L2;UpdateData();x=m_x2,y=m_y2-500000;p=60*60*360/(PI)/2;Be=x/6367558.4969;Bf=Be+(50221746+(293622+(2350+22*cos(Be)*cos(Be)*cos(Be)*cos(Be)*cos(Be)*cos(Be)*1e-10*sin(Be)*cos(Be);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf);Z=y/(Nf*cos(Bf);b2=(0.5+0.003369*cos(Bf)*cos(Bf)*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf)*cos(Bf)*cos(Bf);B=(Bf-(1.0-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2)/2/(PI)*360;/度數(shù)l=(1-(b3-b5*Z*Z)*Z*Z)*Z/2/(PI)*360;L=126+l;/度數(shù)m_du3=(int)L;m_fen3=(int)(L-(int)L)*60);m_miao3=L*3600-m_du3*3600-m_fen3*60;m_du4=(int)B;m_fen4=(int)(B-(int)B)*60);m_miao4=B*3600-m_du4*3600-m_fen4*60;UpdateData(FALSE); 距離改化主要代碼void CMyDlg:OnButton1() double PI=3.14159265358979323846264338327950288; double eep=0.006738525414683;/第二偏心率的平方 double ee=0.006693421622966; double a=6378245.0000000000; double t1, n1, ym,t2,n2, q12,q21,M1,N1 ,M2, N2,R1, R2, x1, y1, x2, y2, B1, B2, Rm,h1 ,h2, D; double Bf, tf, Nf, r12, nf, Bhudu, a12,A12;/a12為坐標(biāo)方位角，r12為子午線收斂角，q12為方向改化，A12為大地方位角 UpdateData(); h1=m_h1,h2=m_h2; D=m_D; x1=m_x1,y1=m_y1-500000,x2=m_x2,y2=m_y2-500000; B1=(m_B1du+m_B1fen/60+m_B1miao/3600)*PI/180; B2=(m_B2du+m_B2fen/60+m_B2miao/3600)*PI/180; ym=0.5*(y1+y2); t1=tan(B1); n1=sqrt(eep*cos(B1)*cos(B1); t2=tan(B2); n2=sqrt(eep*cos(B2)*cos(B2); M1=a*(1-ee)/sqrt(1-ee*sin(B1)*sin(B1)*(1-ee*sin(B1)*sin(B1)*(1-ee*sin(B1)*sin(B1); N1=a/sqrt(1-ee*sin(B1)*sin(B1); M2=a*(1-ee)/sqrt(1-ee*sin(B2)*sin(B2)*(1-ee*sin(B2)*sin(B2)*(1-ee*sin(B2)*sin(B2); N2=a/sqrt(1-ee*sin(B2)*sin(B2); R1=sqrt(M1*N1); R2=sqrt(N2*M2); do Rm=R1; q12=-(x2-x1)*(2*y1+y2-ym*ym*ym/(R1*R1)/(6*R1*R1)-n1*n1*t1*(y2-y1)*ym*ym/(R1*R1*R1); q21= (x2-x1)*(2*y2+y1-ym*ym*ym/(R2*R2)/(6*R2*R2)+n2*n2*t2*(y2-y1)*ym*ym/(R2*R2*R2);a12=atan(y2-y1)/(x2-x1);Bhudu=x1/6367588.4969;Bf=Bhudu+(50221746+(293622+(2350+22*cos(Bhudu)*cos(Bhudu) )*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*cos(Bhudu)*1.0e-10*sin(Bhudu)*cos(Bhudu);tf=tan(Bf);nf=sqrt(eep*cos(Bf)*cos(Bf);Nf=a/sqrt(1-ee*sin(Bf)*sin(Bf);r12=y1*tf/Nf-y1*y1*y1*tf*(1+tf*tf-nf*nf)/(3*Nf*Nf*Nf)+y1*y1*y1*y1*y1*tf*(2+5*tf*tf+3*tf*tf*tf*tf)/(15*Nf*Nf*Nf*Nf*Nf); A12=a12+r12-q12;R1=M1*N1/(N1*cos(A12)*cos(A12)+M1*sin(A12)*sin(A12); while(fabs(Rm-R1)=1.0e-10);A12=A12*180/PI;/化度r12=r12*180/PI;a12=a12*180/PI;m_A12du=(int)A12;m_A12fen=(int)(A12-(int)A12)*60);m_A12miao=A12*3600-m_A12du*3600-m_A12fen*60;UpdateData(FALSE);void CMyDlg:OnButton2() double PI=3.14159265358979323846264338327950288; double B1,B2; double eep=0.006738525414683;/第二偏心率的平方 double ee=0.006693421622966; double a=6378245.0000000000; double P=180*3600/PI; double N1, A1 ,RA, S, D, H1, H2, h1, h2, D12 ,Bm,ym , Mm, Nm, Rm, x1, x2,y2 ,y1,f1 ,f2,f;/h1h2為正常高,D12高斯平面邊長,Bm為平均緯度，Rm平均曲率半徑UpdateData();h1=m_h1,h2=m_h2; x1=m_x1,y1=m_y1-500000,x2=m_x2,y2=m_y2-500000; B1=(m_B1du+m_B1fen/60+m_B1miao/3600)*PI/180; B2=(m_B2du+m_B2fen/60+m_B2miao/3600)*PI/180; D=m_D;/實(shí)測長度 ym=0.5*(y1+y2); A1=(m_A12du+m_A12fen/60+m_A12miao/3600)*PI/180; N1=a/sqrt(1-ee*sin(B1)*sin(B1); RA=N1/(1+eep*cos(B1)*cos(B1)*cos(A1)*cos(A1); H1=m_h1+m_h0; H2=m_h2+m_h0; f=1-(H2-H1)/D*(H2-H1)/D; f1=1+H1/RA; f2=1+H2/RA; S=D*sqrt(f/(f1*f2)+D*D*D/(24*RA*RA); Bm=0.5*(B1+B2); Mm=a*(1-ee)/sqrt(1-ee*sin(Bm)*sin(Bm)*(1-ee*sin(Bm)*sin(Bm)*(1-ee*sin(Bm)*sin(Bm); Nm=a/sqrt(1-ee*sin(Bm)*sin(Bm); Rm=sqrt(Mm*Nm); D12=(1+ym*ym/(2*Rm*Rm)+(y2-y1)*(y2-y1)/(24*Rm*Rm)+ym*ym*ym*ym/(24*Rm*Rm*Rm*Rm)*S; m_S=S; m_D12=D12; UpdateData(FALSE);大地主題正反算主要代碼：void CMyDlg:OnButton1() / TODO: Add your control notification handler code here double ee=0.006693421622966;/第一偏心率的平方 double eep=0.006738525414683;/第二偏心率的平方 double a=6378245.0000000000; double PI=3.14159265358979323846264338327950288; double P=180*3600/PI; double B1,L1,A12,S12,M1,N1,t,n,dertB0,dertL0,dertA0,Bm,Lm,Am,dertB,dertL,dertA,Mm,Nm,tm,nm,B2,L2,A21,Vm,f1,f2,f3,f4; UpdateData(); B1=(m_B1du+m_B1fen/60+m_B1miao/3600)*PI/180;/化弧度 L1=(m_L1du+m_L1fen/60+m_L1miao/3600)*PI/180; A12=(m_A12du+m_A12fen/60+m_A12miao/3600)*PI/180; S12=m_S12; f1=1-ee*sin(B1)*sin(B1); M1=a*(1-ee)/sqrt(f1*f1*f1); N1=a/sqrt(1-ee*sin(B1)*sin(B1); t=tan(B1); n=sqrt(eep*cos(B1)*cos(B1);dertB0=S12*cos(A12)/M1;dertL0=S12*sin(A12)/(N1*cos(B1);dertA0=dertL0*sin(B1);/單位弧度Bm=B1+0.5*dertB0;Lm=L1+0.5*dertL0;Am=A12+0.5*dertA0;f2=1-ee*sin(Bm)*sin(Bm);Mm=a*(1-ee)/sqrt(f2*f2*f2);Nm=a/sqrt(1-ee*sin(Bm)*sin(Bm);tm=tan(Bm);Vm=sqrt(1+eep*cos(Bm)*cos(Bm);nm=sqrt(eep*cos(Bm)*cos(Bm);f4=3*cos(Am)*cos(Am)*nm*nm*(tm*tm-1-nm*nm-4*nm*nm*tm*tm);f3=1+S12*S12/(24*Nm*Nm)*(sin(Am)*sin(Am)*(2+3*tm*tm+2*nm*nm)+f4);dertB=S12*cos(Am)/Mm*f3;f4=tm*tm*sin(Am)*sin(Am)-cos(Am)*cos(Am)*(1+nm*nm-9*nm*nm*tm*tm+nm*nm*nm*nm);f3=1+S12*S12/(24*Nm*Nm)*f4;dertL=S12*sin(Am)/Nm/cos(Bm)*f3;f4=cos(Am)*cos(Am)*(2+7*nm*nm+9*nm*nm*tm*tm+5*nm*nm*nm*nm)+sin(Am)*sin(Am)*(2+tm*tm+2*nm*nm);f3=1+S12*S12/(24*Nm*Nm)*f4;dertA=S12*sin(Am)*tm/Nm*f3;dodertB0=dertB;dertL0=dertL;dertA0=dertA;f4=3*cos(Am)*cos(Am)*nm*nm*(tm*tm-1-nm*nm-4*nm*nm*tm*tm);f3=1+S12*S12/(24*Nm*Nm)*(sin(Am)*sin(Am)*(2+3*tm*tm+2*nm*nm)+f4);dertB=S12*cos(Am)/Mm*(f3);f4=tm*tm*sin(Am)*sin(Am)-cos(Am)*cos(Am)*(1+nm*nm-9*nm*nm*tm*tm+nm*nm*nm*nm);f3=1+S12*S12/(24*Nm*Nm)*f4;dertL=S12*sin(Am)/Nm/cos(Bm)*f3;f4=cos(Am)*cos(Am)*(2+7*nm*nm+9*nm*nm*tm*tm+5*nm*nm*nm*nm)+sin(Am)*sin(Am)*(2+tm*tm+2*nm*nm);f3=1+S12*S12/(24*Nm*Nm)*f4;dertA=S12*sin(Am)*tm/Nm*f3;Bm=B1+0.5*dertB;Lm=L1+0.5*dertL;Am=A12+0.5*dertA;Mm=a*(1-ee)/sqrt(1-ee*sin(Bm)*sin(Bm)*(1-ee*sin(Bm)*sin(Bm)*(1-ee*sin(Bm)*sin(Bm);Nm=a/sqrt(1-ee*sin(Bm)*sin(Bm);tm=tan(Bm);Vm=sqrt(1+eep*cos(Bm)*cos(Bm);nm=sqrt(eep*cos(Bm)*cos(Bm);while(dertB-dertB0)1.0e-10)|(dertL-dertL0)1.0e-10)|(dertA-dertA0)1.0e-10);B2=B1+dertB;B2=B2*180/PI;L2=L1+dertL;L2=L2*180/PI;A21=A12+dertA+PI;A21=A21*180/PI;m_B2du=(int)B2;m_B2fen=(int)(B2-(int)B2)*60);m_B2miao=B2*3600-m_B2du*3600-m_B2fen*60;m_L2du=(int)L2;m_L2fen=(int)(L2-(int)L2)*60);m_L2miao=L2*3600-m_L2du*3600-m_L2fen*60;m_A21du=(int)A21;m_A21fen=(int)(A21-(int)A21)*60);m_A21miao=A21*3600-m_A21du*3600-m_A21fen*60;UpdateData(FALSE);void CMyDlg:OnButton2() / TODO: Add your control notification handler code here double ee=0.006693421622966;/第一偏心率的平方 double eep=0.006738525414683;/第二偏心率的平方 double a=6378245.0000000000; double PI=3.14159265358979323846264338327950288; double P=180*3600/PI; double B1, L1 , B2 ,L2 ,Bm, dertB,dertA, dertL, Mm, Nm ,r01, r21, r03,s10,s12, s30,t01,t21,t03,Am,U,Vm, V,t m, nm ,T,C,S,A1,A2,f1; UpdateData(); B1=(m_b1du+m_b1fen/60+m_b1miao/3600)*PI/180;/化弧度 L1=(m_l1du+m_l1fen/60+m_l1miao/3600)*PI/180; B2=(m_b2du+m_b2fen/60+m_b2miao/3600)*PI/180;/化弧度 L2=(m_l2du+m_l2fen/60+m_l2miao/3600)*PI/180; Bm=0.5*(B1+B2); dertB=B2-B1; dertL=L2-L1; f1=1-ee*sin(Bm)*sin(Bm); Mm=a*(1-ee)/sqrt(f1*f1*f1); Nm=a/sqrt(1-ee*sin(Bm)*sin(Bm); nm=sqrt(eep*cos(Bm)*cos(Bm); tm=tan(Bm); Vm=sqrt(1+eep*cos(Bm)*cos(Bm); r01=Nm*cos(Bm); r21=Nm*cos(Bm)/(24*Vm*Vm*Vm*Vm)*(1+nm*nm-9*nm*nm*tm*tm+nm*nm*nm*nm); r03=-Nm*cos(Bm)*cos(Bm)*cos(Bm)*tm*tm/24; s10=Nm/(Vm*Vm); s12=-Nm*cos(Bm)*cos(Bm)*(2+3*tm*tm+2*nm*nm)/(24*Vm*Vm); s30=Nm*(nm*nm-tm*tm*nm*nm)