大地测量学课程设计

大地测量学课程设计设计题目：白塞尔大地主题解算学院：矿业学院专业：测绘工程班级：********学号：120801024*学生姓名：武*指导教师：张*2014年12月31日目录1.基本原理及思想12.白塞尔法大地主题正算步骤23.白塞尔法大地主题反算步骤44.同一平行圈弧长、子午线弧长与大地线比较大小65.程序代码86.演算示例137.参考文献168.心得体会179.教师评语181白塞尔大地主题解算一：基本原理建立以椭球中心为中心，以任意长（或单位长）为半径的辅助球，按以下三个步骤计算。第一，按一定条件将椭球面元素投影到辅助球面上。第二，在球面上解算大地问题。第三，将求得的球面元素按投影关系换算到相应的椭球元素。关键：确定球面元素与椭球面元素的关系，即它们间的投影关系。二：白塞尔法解算大地主题的基本思想:以辅助球面为基础,将椭球面三角形转换为辅助球面的相应三角形,由三角形对应元素关系,将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，然后在球面上进行大地主题解算，最后再将球面上的计算结果换算到椭球面上。三：在球面上进行大地主题解算球面上的大地主题正算：球面上的大地主题反算：2球面三角元素间的相互关系：四：正反算步骤1.白塞尔法大地主题正算步骤已知a、2e、1B、1L、1A(12A)、S，计算2B、2L、2A(21A)。（1）将椭球面元素投影到球面上由1B求1u：121tan1tanBeu计算辅助量0A和1110sincossinAuA,111sectantanAu计算球面长度，将S化为)24cos(2sin)2cos(sin11S3式中系数分别为：)2561564741(642''''kkkbbA)5123784(642''''kkkAB)12812864''''kkAC（022'2coskAe256564341642kkkA51231285121516464642kkCkkkB上式右端含有，因此需要迭代计算。第一次迭代取近似值S0，第二次计算取)24cos(2sin)2cos(sin01010S以后计算用代换0代入上式迭代计算，直到所要求的精度为止。一般取''00.001|-|。（2）解算球面三角形计算2Asinsincoscoscossincostan111112uAuAuA计算2ucoscoscoscossinsin1112Auuu或)tan(costan122Au计算1111cossinsincoscossinsintanAuuA（3）将球面元素换算到椭球面上由2u求2B422222sin11sintanueuB或22'2tan1tanueB将球面经差化为椭球面经差l，求2Ll)24cos(2sin)2cos(sinsinA-1'1''0式中4'22'22642'1283)1(16)1682(kekeeeee4'22'22'32)1(16kekee4'2'256ke0222'cosAek式中'的最大值为''0.0002，故在计算时通常可以略去不计。象限的判定1sinA符号++--tan符号+-+-||||--||||12LLl1sinA符号--++2tanA符号+-+-2A||2A||2A||2A||22A其中||、||2A为锐角。2.白塞尔法大地主题反算步骤已知a、2e、1B、1L、2B、2L，计算1A(12A)、2A(21A)、S。（1）将椭球面元素投影到球面上5由B求u121tan1tanBeu，222tan1tanBeu，l=12LL211sinsinuua,212coscosuua211sincosuub,212cossinuub采用逐次趋近法，由l计算在反算中，已知椭球面上经差l，球面经差上的对应经差未知，为了由l求，由下式可知还需计算、0A、1，计算又还需量，故需要进行迭代计算。第一次趋近，取l；sincosp2u，cosq21bbsincossinsincossincostan212121uuuuuA或qpA1tan判断1A的象限p符号++--q符号+-+-1A||1A||1A||1A||21A1212112cos)coscossinsin(cossinsincossinAuuuuAucoscoscossinsincos2121uuuucossintanarc判断象限cos+-||||-110sincosusinAA61011111tansinsinutansectantanAAAul+)24cos(2sin)2cos(sin'sin1'1'0A仿照上述计算步骤迭代计算，直到|-|1ii为止。（2）将球面元素换算到椭球面上)24cos(2sin)2cos(sin-111S212112cossincossincossincostanuuuuuA或sinsincoscoscossincostan111112uAuAuA2A象限的判断与前面一致五：同一平行圈弧长、子午线弧长与大地线比较大小子午线弧长计算公式：BaBaBaBaM8cos6cos4cos2cosa86420式中：128a1632a3271638a167321522a12835165832a88866864486422864200mmmmmmmmmmmmmmm平行圈弧长公式：''1''''cosBlblNS（cosBb''1N）不同纬度对应的一些弧长的数值B子午线弧长平行圈弧长1′1″1′1″0°110576m1842.94m30.716m111321m1855.36m30.923m15°1106561844.2630.7381075521792.5429.87630°1108631847.7130.795964881608.1326.802745°1111431852.3930.873788481314.1421.90260°1114231857.0430.95155801930.0215.575°1116251860.4231.00728902481.718.02890°1116961861.631.027000利用白塞尔大地主题反算求解大地线长S纬度为30°，经差为1°的平行圈弧长S=96488m，两点间大地线长为96487.595m经度为30°，纬度差为1°的子午线弧长X=110863m，两点间大地线长为110862.869m通过比较可知，同一平行圈或同一子午线两点间大地线长度与对应的平行圈弧长或子午线弧长相等。8六：程序代码#includestdio.h#includemath.hdoublehudu(double,double,double);/*度分秒转换为弧度*/doubledu(double);/*弧度转换为度*/doublefen(double);/*弧度转换为分*/doublemiao(double);/*弧度转换为秒*/#definePI3.1415926voidmain(void){intk;printf(请选择大地主题算法，若执行正算，请输入1；若执行反算，请输入2。\n);scanf(%d,&k);/*大地主题正算*/if(k==1){doubleax,ay,az,bx,by,bz,cx,cy,cz,S,dz,ez,fz,B1,B2,L1,L2,A1,A2;intdx,dy,ex,ey,fx,fy;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;/*输入度分秒数据*/printf(请输入大地线起点纬度度分秒\n);scanf(%lf%lf%lf,&ax,&ay,&az);printf(请输入大地线起点经度度分秒\n);scanf(%lf%lf%lf,&bx,&by,&bz);printf(请输入大地方位角度分秒\n);scanf(%lf%lf%lf,&cx,&cy,&cz);printf(请输入大地线长度\n);scanf(%lf,&S);/*调用函数*/B1=hudu(ax,ay,az);L1=hudu(bx,by,bz);A1=hudu(cx,cy,cz);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));sinu1=sin(B1)*(sqrt(1-e2))/W1;cosu1=cos(B1)/W1;sinA0=cosu1*sin(A1);coto1=cosu1*cos(A1)/sinu1;sin2o1=2*coto1/(coto1*coto1+1);9cos2o1=(coto1*coto1-1)/(coto1*coto1+1);A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);o0=(S-(B+C*cos2o1)*sin2o1)/A;sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);o=o0+(B+5*C*cos2o)*sin2o/A;g=(r*o+t*(sin2o-sin2o1))*sinA0;/*求B2*/sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));/*求L2*/q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));/*判断q*/if(sin(A1)0&&tan(q)0)q=fabs(q);elseif(sin(A1)0&&tan(q)0)q=PI-fabs(q);elseif(sin(A1)0&&tan(q)0)q=-fabs(q);elseq=fabs(q)-PI;L2=L1+q-g/3600/180*PI;/*求A2*/A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o)));/*判断A2*/if(sin(A1)0&&tan(A2)0)A2=fabs(A2);elseif(sin(A1)0&&tan(A2)0)A2=PI-fabs(A2);elseif(sin(A1)0&&tan(A2)0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/dx=(int)(du(B2));dy=(int)(fen(B2));dz=miao(B2);ex=(i