高斯正反算(零误差)
[转载]高斯正反算
![[转载]高斯正反算](https://img.taocdn.com/s3/m/d13bd7fbe109581b6bd97f19227916888486b9df.png)
[转载]⾼斯正反算⼤地坐标向笛卡尔坐标转换⾼斯正反算采⽤不同椭球实现⾼斯克⾥格投影，将经纬度坐标转换为⾼斯平⾯坐标：正算⾼斯平⾯坐标转换为不同椭球下的经纬度坐标：反算1void GaussProjectDirect(double a,double efang,double B,double L,double L0,double& x,double &y,double& R)//⾼斯投影正算克⽒2 {34double b=aefangtob(a,efang);5double e2=seconde(a,b);6double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f",W);7double N=a/W;printf("N=%f",N);8double M=a*(1-efang)/pow(W,3);printf("M=%f",M);9double t=tee(B);10double eitef=eitefang(a,b,B);11double l=L-L0;12//主曲率半径计算13double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;14 m0=a*(1-efang); n0=a;15 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;16 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;17 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;18 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;19//⼦午线曲率半径20double a0,a2,a4,a6,a8;21 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;22 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;23 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;24 a6=m6/32+m8/16;25 a8=m8/128;2627double X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8*sin(8*B);28 x=X+N/2*t*cos(B)*cos(B)*l*l+N/24*t*(5-t*t+9*eitef+4*pow(eitef,2))*pow(cos(B),4)*pow(l,4)+N/720*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6);29 y=N*cos(B)*l+N/6*(1-t*t+eitef)*pow(cos(B),3)*pow(l,3)+N/120*(5-18*t*t+pow(t,4)+14*eitef-58*eitef*t*t)*pow(cos(B),5)*pow(l,5);30 R=sqrt(M*N);31 }323334//⾼斯投影反算353637void GaussProjectInvert(double a,double efang,double x,double y,double L0,double &B,double& L,double& R)38 {39double b=aefangtob(a,efang);404142double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;43 m0=a*(1-efang); n0=a;44 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;45 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;46 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;47 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;484950//⼦午线曲率半径51double a0,a2,a4,a6,a8;52 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;53 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;54 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;55 a6=m6/32+m8/16;56 a8=m8/128;575859double X=x;60double FBf=0;61double Bf0=X/a0,Bf1=0;62while((Bf0-Bf1)>=0.0001)63 { Bf1=Bf0;64 FBf=a0*Bf0-a2/2*sin(2*Bf0)+a4/4*sin(4*Bf0)-a6/6*sin(6*Bf0)+a8/8*sin(8*Bf0);65 Bf0=(X-FBf)/a0;66 }67double Bf=Bf0;68double Vf=bigv(a,b,Bf);69double tf=tee(Bf);70double Nf=bign(a,b,Bf);71double eiteffang=eitefang(a,b,Bf);72double Bdu=rad_deg(Bf)-1/2.0*Vf*Vf*tf*(pow((y/Nf),2)-1.0/12*(5+3*tf*tf+eiteffang-9*eiteffang*tf*tf)*pow((y/Nf),4)+1.0/360.0*(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6))*180/PI; 73double ldu=1.0/cos(Bf)*(y/Nf+1.0/6.0*(1+2*tf*tf+eiteffang)*pow((y/Nf),3)+1.0/120.0*(5+28*tf*tf+24*tf*tf+6*eiteffang+8*eiteffang*tf*tf)*pow((y/Nf),5))*180.0/PI;747576 B=deg_int(Bdu);77 L=L0+deg_int(ldu);78double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f\n",W); 79double N=a/W;printf("N=%f\n",N);80double M=a*(1-efang)/pow(W,3);printf("M=%f\n",M);81 R=sqrt(M*N);828384 }。
高斯投影正反算
高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程 学号：X51414012：超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差围的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。
高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。
由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。
高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。
二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件 1）中央子午线投影后为直线 2）中央子午线投影后长度不变 3）投影具有正形性质，即正形投影条件2.高斯正算公式推导1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。
2）由于高斯投影是换带投影，在每带经差ｌ是不大的，ｌρ是一个微小量，所以可以将X=X （ｌ，q ），Y=Y （l ，q ）展开为经差为l 的幂级数，它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数，他们都是纬度B 的函数。
3）由第三个条件：∂y ∂l =∂x ∂q 和∂x ∂l =-∂y∂q ，将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩经过计算可以得出232244524632235242225sin cos sin cos (594)224 sin cos (6158)720cos cos (1)6cos (5181458)120N N x X B B l B B t l NB B t t l Ny N B l B t l NB t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。
高斯投影正反算公式
⾼斯投影正反算公式⾼斯投影坐标正反算⼀、基本思想：⾼斯投影正算公式就是由⼤地坐标（L ，B ）求解⾼斯平⾯坐标（x ，y ），⽽⾼斯投影反算公式则是由⾼斯平⾯坐标（x ，y ）求解⼤地坐标（L ，B ）。
⼆、计算模型：基本椭球参数：椭球长半轴a椭球扁率f椭球短半轴：(1)b a f =-椭球第⼀偏⼼率：e a= 椭球第⼆偏⼼率：e b'=⾼斯投影正算公式：此公式换算的精度为0.001m6425644223422)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+''++-''+''?''+=ρηηρρ 5222425532233)5814185(cos 120)1(cos 6cos l t t t B N l t B N l B N y ''-++-''+''+-''+''?''=ηηρηρρ其中：⾓度都为弧度B 为点的纬度，0l L L ''=-，L 为点的经度，0L 为中央⼦午线经度； N 为⼦午圈曲率半径，1222(1sin )N a e B -=-；tan t B =； 222cos e B η'=1803600ρπ''=*其中X 为⼦午线弧长：2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ??=--++-+02468,,,,a a a a a 为基本常量，按如下公式计算：200468242684468686883535281612815722321637816323216128m a m m m m m m a m m m a m m m m a m a ?=++++=+++=++=+ =??02468,,,,m m m m m 为基本常量，按如下公式计算：22222020426486379(1);;5;;268m a e m e m m e m m e m m e m =-====；⾼斯投影反算公式：此公式换算的精度为0.0001’’.()()()()2222243246532235242225053922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f ff f f f f f ft t B B y t t yM N M N t y t t yM N y y l t N B N B y t t t N B L l L ηηηηη=-+++--++=-+++++++=+其中： 0L 为中央⼦午线经度。
高斯坐标正反算
正形投影的一般条件基本出发点：在正形投影中，长度比与方向无关。
1、长度比的通用公式如图4-42，在微分直角三角形P1P2P3及P1′P2′P3′中有：其中l=L-L0，L0通常是中央子午线的经度，L是P点的经度令：()()222222d=d cos dd=d dS M B N B ls x y++(1)m平方可为：()()()22222222222d d d d dd d cos d dcos dcoss x y x ymS M B N B l M BN B lN B++⎛⎫===⎪⎡⎤⎝⎭+⎛⎫+⎢⎥⎪⎝⎭⎢⎥⎣⎦(2)为简化公式，令：ddcosM BqN B=dc o sB M BqN B=⎰(3) q称为等量纬度，因为它只与纬度B有关。
这样，式(2)可表示为：()()222222d dd dx ymr q l+=⎡⎤+⎣⎦(4)我们投影的目的是：建立平面坐标xy和大地坐标BL之间的函数关系，由式(3)可知，即建立xy和bl的函数关系。
令()(),,x x l q y y l q==(5) 对上式进行全微分可得：d d dd d dx xx q lq ly yy q lq l∂∂⎧=+⎪∂∂⎪⎨∂∂⎪=+⎪∂∂⎩(6)将上式代入式(1)中第二项，并令：2222x yEq qx x y yFq l q lx yGl l⎧⎛⎫⎛⎫∂∂=+⎪ ⎪ ⎪∂∂⎪⎝⎭⎝⎭⎪∂∂∂∂⎪=+⎨∂∂∂∂⎪⎪∂∂⎛⎫⎛⎫⎪=+⎪ ⎪∂∂⎪⎝⎭⎝⎭⎩(7)可得： ()()()()222d d 2d d d s E q F q l G l =++ (8) 则式(4)可写为： ()()()()()()222222d 2d d d d d E q F q l G l m r q l ++=⎡⎤+⎣⎦ (9)2 柯西-黎曼条件在上式引入方向，如图4-42所示：2313d d cot d d P P M B q A PP r l l === (10) 即： d tan d l A q = (11)将式(11)代入式(9)可得：注意sec 1cos A A =()()()()()222222222222222d 2tan d tan d d tan d 2tan tan sec cos 2sin cos sin E q F A q G A q m r q A q E F A G A r AE AF A AG A r ++=⎡⎤+⎣⎦++=++=(12)要想让m 和A 无关，必须使F=0，E=G ，即22220x x y y q l q l x y x y q q l l ∂∂∂∂⎧+=⎪∂∂∂∂⎪⎨⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫⎪+=+⎪ ⎪ ⎪ ⎪⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎩ (13) 由上式第一式可得：y y x q l x lq∂∂∂∂∂=-∂∂∂(14)代入第二式可得： 222222y x y y x lq q q q x q ∂⎛⎫⎪⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫∂∂∂∂∂⎝⎭+=+⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎛⎫∂⎣⎦⎪∂⎝⎭(15) 消去公共项可得： 22x y q l ⎛⎫∂∂⎛⎫= ⎪ ⎪∂∂⎝⎭⎝⎭ (16)开方并代入式(13)的第一项：x y q l x y lq ∂∂⎧=⎪∂∂⎪⎨∂∂⎪=-⎪∂∂⎩ (17)高斯投影坐标正算高斯投影三条件：L0为直线；L0长度不变；正形投影 1、幂级数展开公式（x 偶y 奇）l /ρ微小量(ρ''=206265)，可进行级数展开，可得：2402435135x m m l m l y m l m l m l ⎧=+++⎪⎨=+++⎪⎩(18) 式中mi 为待定系数，是q 、B 的函数。
高斯正算和高斯反算精度
高斯正算和高斯反算精度高斯正算和高斯反算是地理测量中常用的计算方法，用于求解地球上两点之间的距离、方位角和坐标等问题。
高斯正算是根据已知起点坐标、方位角和距离，计算出终点的坐标；而高斯反算则是根据已知起点和终点的坐标，计算出两点之间的距离和方位角。
高斯正算是地理测量中常用的一种计算方法，通过已知起点的坐标、方位角和距离，计算出终点的坐标。
这种方法主要基于高斯投影的理论，通过将地球表面的曲面投影到平面上，使得计算更加简便。
在进行高斯正算时，首先需要确定起点的坐标、方位角和距离，然后通过一系列的计算步骤，得到终点的坐标。
高斯正算可以应用于各种地理测量问题，例如导航定位、地图制作等。
在进行高斯反算时，已知起点和终点的坐标，需要计算出两点之间的距离和方位角。
与高斯正算相反，高斯反算是从已知坐标求解未知距离和方位角的过程。
通过一系列的计算步骤，可以得到两点之间的距离和方位角。
高斯反算也是地理测量中常用的计算方法，可以应用于测量控制点、确定地理位置等方面。
高斯正算和高斯反算的精度是地理测量中非常重要的指标。
精度是指测量结果与真实值之间的偏差程度，通常用误差来表示。
在进行高斯正算和高斯反算时，精度的提高可以通过多种方式实现。
首先，使用更加精确的测量设备和仪器可以提高计算结果的精度。
其次，选择合适的计算方法和算法也可以提高精度。
此外，对输入数据的准确性进行检查和验证，以及进行误差分析和修正，也是提高精度的关键步骤。
在实际应用中，高斯正算和高斯反算的精度要求因任务而异。
对于一般的地理测量任务，如导航定位、地图制作等，一般要求精度在几米到十几米之间。
而对于精密测量任务，如大地测量、工程测量等，精度要求可能在毫米甚至更高。
为了满足不同任务的精度要求，需要根据具体情况选择合适的测量方法和技术，以及进行精确的数据处理和分析。
高斯正算和高斯反算是地理测量中常用的计算方法，用于求解地球上两点之间的距离、方位角和坐标等问题。
在进行高斯正算和高斯反算时，需要选择合适的计算方法和算法，并注意精度的控制和提高。
(整理)高斯投影坐标正反算公式
高斯投影坐标正反算公式未知2010-04-03 10:47:15 本站§高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（ C-R 偏微分方程），还有它本身的特殊条件。
1.1 高斯投影坐标正算公式： B, x,y高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。
由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即 (8-10) 式中， x 为的偶函数， y 为的奇函数；，即，如展开为的级数，收敛。
（ 8-33 ）式中是待定系数，它们都是纬度 B 的函数。
由第三个条件知：(8-33) 式分别对和 q 求偏导数并代入上式(8-34)上两式两边相等，其必要充分条件是同次幂前的系数应相等，即(8-35)(8-35) 是一种递推公式，只要确定了就可依次确定其余各系数。
由第二条件知 : 位于中央子午线上的点，投影后的纵坐标 x 应等于投影前从赤道量至该点的子午线弧长 X ，即 (8-33) 式第一式中，当时有：(8-36)顾及 ( 对于中央子午线 )得：(8-37,38)(8-39)依次求得并代入 (8-33) 式，得到高斯投影正算公式(8-42)1.2 高斯投影坐标反算公式x,y B,投影方程：(8-43)满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。
高斯投影坐标反算公式推导要复杂些。
①由 x 求底点纬度 ( 垂足纬度 ), 对应的有底点处的等量纬度，求 x,y 与的关系式，仿照 (8-10) 式有，由于 y 和椭球半径相比较小 (1/16.37) ，可将展开为 y 的幂级数；又由于是对称投影， q必是 y 的偶函数，必是 y 的奇函数。
(8-45)是待定系数，它们都是 x 的函数 .由第三条件知：，， (8-21)(8-45) 式分别对 x 和 y 求偏导数并代入上式上式相等必要充分条件，是同次幂 y 前的系数相等，第二条件，当 y=0 时，点在中央子午线上，即 x=X ，对应的点称为底点，其纬度为底点纬度，也就是 x=X 时的子午线弧长所对应的纬度，设所对应的等量纬度为。
高斯公式正负
高斯公式正负
高斯公式是数学和物理中非常重要的一个公式，它可以用来计算积分，特别是在计算一些难以直接求解的积分时。
正负号的判断在高斯公式中是一个关键问题，因为它决定了积分的方向和曲面的闭合性。
如果所积分的曲面是闭合的，那么方向就是负号。
这是因为闭合曲面内的积分值是有限的，而闭合曲面外的积分值是无限的，所以方向向内是负号。
相反，如果所给的曲面不是闭合的，这时需要作辅助面使其成为闭合的曲面。
这时，方向向里为负号，外为正号。
这是因为辅助面内的积分值是有限的，而辅助面外的积分值也是有限的，所以方向向外是正号。
高斯公式的正负号判断是一个相对复杂的问题，需要综合考虑曲面的形状、积分的方向以及积分的范围等因素。
在实际应用中，需要根据具体情况进行判断，而且有时候可能需要借助一些辅助工具或方法来帮助判断正负号。
总之，正确理解高斯公式的正负号判断对于正确应用高斯公式非常重要。
高斯投影正反算
高斯投影正反算LT2340123423401234...........q m m y m y m y m y l n n y n y n y n y '''''=+++++'''''=+++++3.引入高斯投影条件之一：正形条件4.由于可得到00n '=，带入上式可得到 2402435135...........q m m y m y l n y n y n y '''=+++'''=+++5.引入高斯投影条件之三：中央子午线投影后长度不变001122223322442255sec sec 122sec (12)6sec (56)24sec (5286)120f f f f f f f f f f f f f f f f f f f m q B dm n dx N t B dn m dx N B n t N t B m t N B n t N ηηη'=⎧⎪'⎪'==⎪⎪⎪''=-=-⎪⎪⎪⎨'=-++⎪⎪⎪'=++⎪⎪⎪⎪'=++⎪⎩四、高斯投影的特点1.当l 等于常数时，随着B 的增加x 的值增大，y 的值减小，无论B 值为正或为负，y 值不变。
这就是说，椭球面上除中央子午线外，其它子午线投影后，均向中央子午线弯曲，并向两极收敛，同时还对称于中央子午线和赤道。
2.当B 等于常数时，随着l 的增加，x 值和y 值都增大。
所以在椭球面上对称于赤道的纬圈，投影后仍成为对称的曲线，同时与子午线的投影曲线互相垂直凹向两极。
3.据中央子午线越远的子午线，投影后弯曲越厉害，长度变形越大。
五、MATLAB编程实现坐标正反算1.编写main函数function maindisp('欢迎使用高斯投影正反算及相邻带的坐标换算程序');disp('1：高斯正算 2：高斯反算 3：换带计算');K=0;while (K<1||K>3)K=input('请根据上列选择计算类型 K=');switch Kcase 1GSZS;case 2GSFS;case 3HDJS;otherwisedisp('K 值无效（1-3）');enddisp('程序作者：亚里士多墩');disp('指导老师：亚里士多德');end2.编写高斯正算GSZS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3-0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效（1-3）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24; xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End3.编写高斯反算公式GSFS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3 -0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效（1-3）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24;xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4); r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End六、代码测试1.L=111°47'24〞.8974，B=31°04'41〞.6832，L0=111°2.L=111.47532575，B=31.23484275，L0=1113.L=114.20，B=30.30，L0=1114.L=118.54152206，B=32.24576522，L0=117。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
//84的椭球final double a = 6378137;final double Alfa = 1.0 / 298.257223563;double centreL, x, y, b, e1, ee;double a0, a2, a4, a6, a8, Bf0;double[] Coeficient_a0 = new double[5];double sinBf, cosBf;double FBf, Bf1, dB, bf;double c, v, Nf, Mf, tf;double itaf, dietaB, dietaL;double B, L;double dmsB, dmsL, dmsCentreL1;double radlat, radlon, radl0, l;double sb, cb, t, ita, X, N;public String MakeProject(double L, double B, double CentreLon) //高斯正算{/*输入已知数据:经度\纬度\ 中央子午线*/dmsB = B;dmsL = L;dmsCentreL1 = CentreLon;radlat = DMSTORAD(dmsB);radlon = DMSTORAD(dmsL);radl0 = DMSTORAD(dmsCentreL1);l = radlon - radl0;b = a * (1 - Alfa);sb = Math.sin(radlat);cb = Math.cos(radlat);t = sb / cb;e1 = Math.sqrt((a / b) * (a / b) - 1);ee = Math.sqrt(1 - (b / a) * (b / a));ita = e1 * cb;a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a6 = Coeficient_a0[3];a8 = Coeficient_a0[4];X = a0 * radlat - sb * cb * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sb * sb + 16 * a6 * sb * sb * sb * sb / 3.0);//计算卯酉圈半径Nc = a * a / b;v = Math.sqrt(1 + e1 * e1 * cb * cb);N = c / v;//计算未知点的坐标x = X + N * sb * cb * l * l / 2 + N * sb * cb * cb * cb * (5 - t * t + 9 * ita * ita + 4 * ita * ita * ita * ita) * l * l * l * l / 24 + N * sb * cb * cb * cb * cb * cb * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / 720;y = N * cb * l + N * cb * cb * cb * (1 - t * t + ita * ita) * l * l * l / 6 + N * cb * cb * cb * cb * cb * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * l * l * l * l * l / 120;//输出未知点坐标return (y + 500000)+","+x;}public Point MoveProject(double CentreLon, double X1, double Y1){centreL = CentreLon * Math.PI / 180.0;x = X1;while(Y1>1000000)Y1=Y1-1000000;y = Y1 - 500000.0;b = a * (1 - Alfa);e1 = Math.sqrt((a / b) * (a / b) - 1);//第二偏心率ee = Math.sqrt(1 - (b / a) * (b / a));a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a4 = Coeficient_a0[2];a8 = Coeficient_a0[4];Bf0 = x / a0;do{sinBf = Math.sin(Bf0); cosBf = Math.cos(Bf0);FBf = -sinBf * cosBf * ((a2 - a4 + a6) + (2.0 * a4 - 16.0 * a6 / 3.0) * sinBf * sinBf +(16.0 / 3.0) * a6 * (sinBf * sinBf * sinBf * sinBf));Bf1 = (x - FBf) / a0;dB = Bf1 - Bf0;Bf0 = Bf1;} while (Math.abs(dB * 180.0 / Math.PI * 3600.0) > 0.00001);bf = Bf1;c = a * a / b;v = Math.sqrt(1 + e1 * e1 * Math.cos(bf) * Math.cos(bf));Nf = c / v;Mf = c / (v * v * v);tf = Math.sin(bf) / Math.cos(bf);itaf = e1 * Math.cos(bf);dietaB = tf * Math.pow(y, 2) / (2 * Mf * Nf) - tf * (5 + 3 * Math.pow(tf, 2) + Math.pow(itaf, 2) - 9 * Math.pow(tf, 2) * Math.pow(itaf, 2)) * Math.pow(y, 4) / (24 * Mf * Math.pow(Nf, 3)) + (61 + 90 * Math.pow(tf, 2) + 45 * Math.pow(tf, 4)) * Math.pow(y, 6) / (720 * Mf * Math.pow(Nf, 5));dietaL = y / (Nf * Math.cos(bf) + (1 + 2 * tf * tf + itaf * itaf) * Math.cos(bf) * y * y / (6 * Nf)) + (5 + 44 * tf * tf + 32 * tf * tf * tf * tf - 2 * itaf * itaf - 16 * itaf * itaf * tf * tf) / (360 * Nf * Nf * Nf * Mf * Mf * Math.cos(bf));B = bf - dietaB;L = centreL + dietaL;dmsB = RADTODMS(B);dmsL = RADTODMS(L);return new Point(dmsL , dmsB);}public void a0a2a4a6a8(double a, double e, double[] Coeficient_a0){double m0, m2, m4, m6, m8;m0 = a * (1 - e * e);m2 = 3 * e * e * m0 / 2;m4 = 5 * e * e * m2 / 4;m6 = 7 * e * e * m4 / 6;m8 = 9 * e * e * m6 / 8;Coeficient_a0[0] = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;Coeficient_a0[1] = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;Coeficient_a0[2] = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;Coeficient_a0[3] = m6 / 32 + m8 / 16;Coeficient_a0[4] = m8 / 128;}static int intSignOfDms, intSignOfRad;static double radAngle, dmsAngle;public static double DMSTORAD(double dmsAngle1){intSignOfDms = 1;if (dmsAngle1 < 0) intSignOfDms = -1;dmsAngle1 = Math.abs(dmsAngle1);radAngle = dmsAngle1 * Math.PI / 180.0;radAngle = radAngle * intSignOfDms;return radAngle;}public static double RADTODMS(double radAngle){intSignOfRad = 1;if (radAngle < 0) intSignOfRad = -1;radAngle = Math.abs(radAngle);dmsAngle = radAngle * 180 / Math.PI;dmsAngle = dmsAngle * intSignOfRad;return dmsAngle;}。