高斯正反算计算函数
[转载]高斯正反算
![[转载]高斯正反算](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 }。
高斯正反算计算函数
高斯正反算计算函数高斯正反算（Gauss forward and backward calculation）是一种用于大地测量中进行坐标计算的方法。
该方法以德国数学家高斯（Carl Friedrich Gauss）的名字命名，用于将球面坐标（经度、纬度、大地高）转换为平面直角坐标（X、Y、Z）或反之。
正算（forward calculation）是将球面坐标转换为平面直角坐标。
该过程包括以下几个步骤：1.首先，确定球面坐标的基准，通常选择一个已知的起始点。
2.然后，计算球面坐标与基准点之间的角度差值，并将其转换为弧度。
3.使用三角函数，根据球面坐标的经度、纬度以及角度差值，计算出平面直角坐标的X、Y、Z值。
反算（backward calculation）是将平面直角坐标转换为球面坐标。
该过程与正算相反，包括以下几个步骤：1.确定平面直角坐标的基准点，即已知的起始点。
2.计算平面直角坐标与基准点之间的距离差值，通常使用勾股定理计算距离。
3.使用三角函数，根据平面直角坐标的X、Y、Z值以及距离差值，计算出球面坐标的经度、纬度以及大地高。
高斯正反算的原理是基于球面三角学和球面坐标系的转换公式。
在正算中，通过球面三角学的公式计算球面坐标与基准点之间的角度差值，并使用三角函数计算平面直角坐标。
在反算中，通过勾股定理计算平面直角坐标与基准点之间的距离差值，并使用三角函数计算球面坐标。
总结起来，高斯正反算是一种用于大地测量中进行坐标计算的方法，通过球面三角学和转换公式将球面坐标和平面直角坐标进行转换。
正算将球面坐标转换为平面直角坐标，反算将平面直角坐标转换为球面坐标。
它在地理信息系统、地图制图等领域有着广泛的应用。
高斯平面直角坐标与大地坐标的相互转换——高斯投影的正算与反算.
昆明冶金高等专科学校测绘学院 （4）计算公式
3 2 2 2 4 ( 5 3 t 9 t ) y f f f f 2M f N f 2 4M f N 3 f tf 2 4 6 (6 1 9 0t f 4 5t f ) y 7 2 0M f N 5 f 1 1 2 2 3 l y (1 2t f f ) y 3 N f co s B f 6 N f co s B f 1 2 5 (5 2 8t 2 t4 2 2 f 24 f 6 f 8 f t f )y 5 1 2 0N f co s B f B Bf tf y2 tf
式中：
2 e 2 cos2 B
t 2 tan2 B l （L L0） X为B对应子午线弧长 N为卯酉圈曲率半径 20626 5
昆明冶金高等专科学校测绘学院
2
高斯投影坐标反算公式
（1）高斯投影反算：
已知某点 x, y ，求该点 L, B ，即 x, y ( L, B) 的坐标变换。 （2）投影变换必须满足的条件
昆明冶金高等专科学校测绘学院
二、高斯投影坐标正反算得实用公式及算例
1 高斯投影坐标正算公式 （1）高斯投影正算： 已知某点的 L, B ，求该点的 x, y ，即 （2）投影变换必须满足的条件： 中央子午线投影后为直线； 中央子午线投影后长度不变； 投影具有正形性质，即正形投影条件。 （3）投影过程 在椭球面上有对称于中央子午线的两点 P1 和 P2 ,它们的大地坐标 分别为 ( L1 , B1 )或（l1 , B1）及 （L2 , B2）或(l2 , B2 ) 式中 l 为椭球面上点的经 度与中央子午线 ( L0 ) 的经度差：l L L0 ,点在中央子午线之东, l 为正，在西则为负，则投影后的平面坐标一定为P1 ( x1 , y1 ) 和 P2 ( x 2 , y 2 ) 。
GAUSSLE坐标正反算fx
GAUSSLE坐标正反算fx－4850程序源程序1.正算主程序GSZSI"X0"：S"Y0"：O"K0"：G"F0"：H"KN"：P"R0"：R"RN"：Q”Q(-Z +Y)” :D=(P-R)÷(2(H-O)PR)：KL”L(-Z +Y)” :M”ANG(YJJ)”=90:（注：此处若不给M赋值，则可计算斜交点）J=Abs(K-O)：Prog"SUB1"：”FWJ=”:F=F-M:”X=”:X=X:Pause0: ”Y=”:Y=Y◢2. 反算主程序GSFSXY：Z[2]=X：Z[3] =Y：I"X0"：S"Y0"：O"K0"：G"F0"：H"KN"：P"R0"：R"RN"：Q”Q(-Z +Y)” :D=(P-R)÷(2(H-O)PR):J=Abs((Y-S)cos(G-90)－(X-I)sin(G-90))：L=0：M”M(YJJ)”=90:Lbl 0：Prog "SUB1"：L=(Z[3]-Y)cos(G-90+QJ(1÷P+JD)×180÷π)-(Z[2]-X)sin(G-90+QJ(1÷P +JD) ×180÷π)：AbsL<1E-6=>Goto1：≠>J=J+L：Goto 0Δ←┘Lbl 1：L=0：Prog "SUB1"：L=(Z[3]-Y)÷sinF:”K=”:K=O+J:Pause0:”L=”:L=L◢3. 正算子程序(SUB1)Defm 4:A=0.1184634425:B=0.2393143352:Z[4]=0.2844444444:C=0.046910077 0:E=0.2307653449:Z[1]=0.5:X=I+J(Acos(G+QCJ(1÷P+C JD)×180÷π)+Bcos(G+Q EJ(1÷P+E JD)×180÷π) +Z[4]cos(G+QZ[1]J(1÷P+Z[1]JD)×180÷π)+Bcos(G+Q(1-E)J(1÷P+(1-E)J D)×180÷π)+Acos(G+Q (1-C)J(1÷P+(1-C)JD) ×180÷π)):Y=S+J(Asin(G+QCJ(1÷P+C JD)×180÷π)+Bsin(G+Q EJ(1÷P+E JD)×180÷π) +Z[4]sin(G+QZ[1]J(1÷P+Z[1]JD)×180÷π)+Bsin(G+Q(1-E)J(1÷P+(1-E)J D)×180÷π)+Asin(G+Q (1-C)J(1÷P+(1-C)JD) ×180÷π))：F=G+QJ(1÷P+JD) ×180÷π+M：X=X+LcosF：Y=Y+LsinF4. 曲线元要素数据库：DAT-MK≥O=>K＜H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘K≥O=>K＜H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘K≥O=>K＜H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘K≥O=>K＜H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘K≥O=>K≤H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘……………………………K≥O=>K≤H=> I=**:S=**:O=**:G=**:H=**:P=**:R=**:Q=**⊿⊿←┘(注：如有多个曲线元要素继续添加入数据库DAT-M中)5、M线（坐标正算）组合程序MG-ZBProg”DAT-M”:Prog”GSZS”6、M线（坐标计算－放样）组合程序MG-FYProg”MG-ZB”:Prog”LTKZD”: Prog”FY”7、M线（坐标反算）组合程序M-GSFBProg”DAT-M”:Prog”GSFS”说明：一、程序功能及原理1.功能说明：本程序由两个主程序——正算主程序(GSZS)、反算主程序(GSFS)和两个子程——正算子程序(SUB1)、线元数据库(DAT-M)构成，可以根据曲线段——直线、圆曲线、缓和曲线（完整或非完整型）的线元要素（起点坐标、起点里程、起点切线方位角、终点里程、起点曲率半径、止点曲率半径）及里程边距或坐标，对该曲线段范围内任意里程中边桩坐标进行正反算。
高斯投影正反算公式
高斯投影坐标正反算一、基本思想：高斯投影正算公式就是由大地坐标（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 为中央子午线经度。
第四章 7高斯投影坐标正反算
——正形投影的一般条件 ——高斯投影坐标正算 ——高斯投影坐标反算 ——高斯投影几何解释
提前在黑板上写出四个m2
上一讲应掌握的内容
1、地图(数学)投影：将椭球面上元素(包括坐标，方位和 距离)按一定的数学法则投影到可展平面上。 x F1 ( L, B) 坐标投影公式： y F2 ( L, B) 2、地图投影变形几个概念： 长度比，主方向，变形椭圆 3、四种投影变形： 长度变形，方向变形，角度变形，面积变形
x m0 m 2 l 2 m 4 l 4 y m1l m3 l 3 m5 l 5
分别对l 和q 求偏导数
2) 由第三个条件正形投影条件
y x x y 和 l q l q
dm0 dm2 2 dm4 4 2 4 m1 3m3 l 5m5 l dq dq l dq l 2m l 4m l 3 dm1 l dm3 l 3 ）
将各系数代入，略去高次项，得高斯投影坐标正算公式 精度为0.001m
xX N N sin B cos Bl 2 + sin B cos 3 B(5 - t 9 2 4 4 )l 4 + 2 24
N sin B cos 5 B(61 - 58t 2 t 4 )l 6 720
dl tan Adq
2 2 2 2 E ( dq ) 2 F tan A ( dq ) G tan A ( dq ) m2 2 2 2 r2 ( dq ) tan A ( dq )
E 2 F tan A G tan 2 A = r 2 sec 2 A E cos 2 A 2 F sin A cos A G sin 2 A = r2
高斯投影坐标反算公式
大地方位角＝坐标方位 角－子午线收敛角＋方 向改化 A1 2 1 2 1 2
高斯坐标反算实用步骤
1、根据高斯坐标确定带号、计算中央子午 线经度 ①计算带号
n int( y / 1000000 )
②计算中央子午线经度 六度带 L 6n 3 0 三度带 L0 3n
11 n 24 23 n 49
0.0067 l 2 ]l 2 cos2 B}l sinB )
对于1975国际椭球
{1 (C3 C5l 2 )l 2 cos2 B}l sin B
C3 0.33332 0.00678 cos2 B C5 0.2 cos2 B 0.0667
计算子午线收敛角的意义： 1、用于大地方位角和高斯平面坐标方位角的转换； 2、高斯正反算检核坐标计算的正确性。。
B Bf t3 f 2 4M f N
3 f
tf 2M f N f
y
2
5 3t
2 f
2 9 2 t 2 y 4 f f f
5、中央子午线收敛角和经度
2 y 2 2 y tan B f [1 (1 tan B f f )] 3 Nf 3N f
1 tan sin B l (1 t 2 3 2 2 4 ) sin B cos 2 Bl 3 3 1 (2 4t 2 2t 4 ) sin B cos 4 Bl 5 15
sin B l sin B cos 2 Bl 3 (1 3 2 2 4 )
L L0 l
小结
• • • • 了解反算公式的推导思路； 掌握反算公式保符号的意义； 用反算公式会进行计算； 掌握子午线收敛角的定义及作用；
4850-5800高斯5点通用公式正反算(验证正确)适用于各种线形
Gauss-Legendre 5点通用公式坐标正反算CASIOfx-4850P程序一、(4850程序)正反算主程序(ZFS)Deg↵"1.KL⇒XY"："2.XY⇒KL"：N：I"X0"：S"Y0"：○"K0"：G"F0"：H"KN"：P" R0"：R"RN"：Q"Q(-Z+Y)":D=(P-R)÷(2(H-○)PR): M=180÷π:N=1⇒Goto1：≠⇒Goto2↵Lbl1：{KL}：KL：J=Abs(K-○)：Prog"SUB1"："X=":X◢"Y=":Y◢"FWJ=":F -90→DMS◢Goto1↵Lbl2：{XY}：XY：Z[2]=X：Z[3]=Y：Prog"SUB2"："K=":○+J◢"L=":L◢Goto2↵2.正算子程序(SUB1)A=0.1184634425:B=0.2393143352:Z[4]=0.2844444444:C=0.046910077:E= 0.2307653449:Z[1]=0.5: X=I+J(Acos(G+QCM J(1÷P+CJD))+Bcos(G+QE MJ (1÷P+EJD))+Z[4]cos(G+QZ[1]J M(1÷P+Z[1]JD))+Bco s(G+Q(1-E)JM(1÷P +(1-E)JD))+Acos(G+Q(1-C)JM(1÷P+(1-C)JD))):Y=S+J(Asin(G+QCJM(1÷P+CJD))+Bsin(G+QEJM(1÷P+EJD))+Z[4]sin(G+QZ[1]J M(1÷P+Z[1]JD))+B sin(G+Q(1-E)JM(1÷P+(1-E)JD))+Asin(G+Q(1-C)JM(1÷P+(1-C)JD)))：F= G+QJM(1÷P+JD) +90：X=X+LcosF：Y=Y+LsinF↵3. 反算子程序(SUB2)J=Abs((Y-S)cos(G-90)－(X-I)sin(G-90))：L=0：Lbl0：Prog"SUB1"：L=(Z [3]-Y)cos(G-90+QJM(1÷P+JD))-(Z[2]-X)sin(G-90+QJM(1÷P+JD))：AbsL <1E-6⇒Goto1：≠⇒J=J+L：Goto0↵Lbl1：L=0：Prog"SUB1"：L=(Z[3]-Y)÷sinF↵二、(4850程序带数据库)正反算主程序(ZFS)反算时需输入里程以判断属于哪段线元Deg↵"1.KL⇒XY"："2.XY⇒KL"：N：M=180÷π:N=1⇒Goto1：≠⇒Goto2↵Lbl1：{TL}：TL：Prog"DATA":D=(P-R)÷(2(H-○)PR):J=Abs(K-○)：Pro g"SUB1"："X=":X◢"Y=":Y◢"FWJ=":F-90→DMS◢Goto1↵Lbl2：{XYT}：XYT：Prog"DATA":D=(P-R)÷(2(H-○)PR):Z[2]=X：Z[3]=Y：Prog"SUB2"："K=":○+J◢"L=":L◢Goto2↵2.正算子程序(SUB1)----与无数据库版相同3.反算子程序(SUB2)----与无数据库版相同4.数据库(DATA)---锦赤铁路T≤172554.99⇒I=4638500.956:S=536490.9244:○=167000:G=-77°50′0. 9″:H=172554.99:P=10∧45:R=10∧45:Q=0:Goto1T≤172624.99⇒I=463967 1.6805:S=531060.7017:○=172554.99:G=-77°50′0.9″:H=172624.99:P =10∧45:R=2000:Q=-1:Goto1T≤172840.8458⇒I=4639686.033:S=530992. 19:○=172624.99:G=-78°50′10.53″:H=172840.8458:P=2000:R=2000:Q =-1:Goto1T≤172910.8458⇒I=4639716.328:S=530778.5764:○=172840. 8458:G=-85°01′12.26″:H=172910.8458:P=2000:R=10∧45:Q=-1:Goto1T≤177000⇒I=4639721.59:S=530708.7755:○=172910.8458:G=-86°01′21.89″:H=177000:P=10∧45:R=10∧45:Q=0:Goto1↵Lbl1:T<167000⇒"CHAO XIAN"◢≠⇒T>177000⇒"CHAO XIAN"◢三、(5800程序)正反算主程序(ZFS)Deg↵"1.KL⇒XY"："2.XY⇒KL"：?N："X0"?I："Y0"?S："K0"?○："F0"?G："KN"? H："R0"?P："RN"?R："Q(-Z+Y)"?Q:(P-R)÷(2(H-○)PR)→D:180÷π→M:N =2⇒Goto2↵Lbl1：?K:?L：Abs(K-○)→J：Prog"SUB1"："X=":X◢"Y=":Y◢"FWJ=":F-9 0►DMS◢Goto1↵Lbl2：?X:?Y：X→Z[2]：Y→Z[3]：Prog"SUB2"："K=":○+J→K◢"L=":L◢Goto2↵2.正算子程序(SUB1)0.1184634425→A:0.2393143352→B:0.2844444444→Z[4]:0.046910077→C:0.2307653449→E:0.5→Z[1]:I+J(Acos(G+QCJM(1÷P+CJD))+Bcos(G+QE JM(1÷P+EJD))+Z[4]cos(G+QZ[1]J M(1÷P+Z[1]JD))+Bcos(G+Q(1-E)JM(1÷P+(1-E)JD))+Acos(G+Q(1-C)JM(1÷P+(1-C)JD)))→X:S+J(Asin(G+QCJM (1÷P+CJD))+Bsin(G+QEJ M(1÷P+EJD))+Z[4]sin(G+QZ[1]J M(1÷P+Z[1]J D))+Bsin(G+Q(1-E)JM(1÷P+(1-E)JD))+Asin(G+Q(1-C)JM(1÷P+(1-C)JD)))→Y:G+QJM(1÷P+JD) +90→F：X+Lcos(F)→X：Y+Lsin(F)→Y↵3. 反算子程序(SUB2)Abs((Y-S)cos(G-90)－(X-I)sin(G-90))→J：0→L↵Lbl0：Prog"SUB1"：(Z[3]-Y)cos(G-90+QJM(1÷P+JD))-(Z[2]-X)sin(G-90 +QJM(1÷P+JD))→L：If Abs（L）<0.000001：Then0→L：Prog"SUB1"：(Z[3] -Y)÷sin(F)→L:Else J+L→J：Goto0：IfEnd↵四、(5800程序带数据库)正反算主程序(ZFS)Deg↵"1.KL⇒XY"："2.XY⇒KL"：?N：180÷π→M:N=2⇒Goto2↵Lbl1:"KP"？T:?L:Prog"DATA"：(P-R)÷(2(H-○)PR)→D:Abs(T-○)→J:Pr og"SUB1"："X=":X◢"Y=":Y◢"FWJ=":F-90►DMS◢Goto1↵Lbl2：?X:?Y:"KP"？T:Prog"DATA"：(P-R)÷(2(H-○)PR)→D：X→Z[2]：Y →Z[3]：Prog"SUB2"："K=":○+J→T◢"L=":L◢Goto2↵2.正算子程序(SUB1)---与无数据库版相同3.反算子程序(SUB2)---与无数据库版相同4.数据库(DATA)---锦赤铁路167000→○:172554.99→H:If T≤H:Then4638500.956→I:536490.9244→S:-77°50′0.9″→G:10^(45)→P:10^(45)→R:0→Q:Goto1：IfEnd↵H→○:172624.99→H:If T≤H:Then 4639671.6805→I:531060.7017→S:-77°50′0.9″→G:10^(45)→P:2000→R:-1→Q:Goto 1:IfEnd↵H→○:172840.8458→H:If T≤H:Then 4639686.033→I:530992.19→S:-7 8°50′10.53″→G:2000→P:2000→R:-1→Q:Goto 1:IfEnd↵H→○:172910.8458→H:If T≤H:Then 4639716.328→I:530778.5746→S: -85°01′12.26″→G:2000→P:10^(45)→R:-1→Q:Goto 1:IfEnd↵H→○:177000→H:If T≤H:Then 4639721.59→I:530708.7755→S:-86°0 1′21.89″→G:10^(45)→P:10^(45)→R:0→Q:Goto 1:IfEnd↵Lbl1:If T＜167000 or T＞177000:Then"CHAO XIAN"◢Stop:IfEnd:Return↵说明：一、程序功能及原理1.功能说明：本程序由一个主程序和两个子程——正算子程序(SUB1)、反算子程序(SUB2)构成，可以根据曲线段——直线、圆曲线、缓和曲线（完整或非完整型）的线元要素（起点坐标、起点里程、起点切线方位角、终点里程、起点曲率半径、止点曲率半径）及里程边距或坐标，对该线元段范围内任意里程中边桩坐标进行正反算。
python实现高斯投影正反算方式
python实现⾼斯投影正反算⽅式使⽤Python实现了⼀下我们同事的C++⾼斯投影正反算，实际跑通，可⽤。
#!/ usr/bin/python# -*- coding:utf-8 -*-import mathdef LatLon2XY(latitude, longitude):a = 6378137.0# b = 6356752.3142# c = 6399593.6258# alpha = 1 / 298.257223563e2 = 0.0066943799013# epep = 0.00673949674227#将经纬度转换为弧度latitude2Rad = (math.pi / 180.0) * latitudebeltNo = int((longitude + 1.5) / 3.0) #计算3度带投影度带号L = beltNo * 3 #计算中央经线l0 = longitude - L #经差tsin = math.sin(latitude2Rad)tcos = math.cos(latitude2Rad)t = math.tan(latitude2Rad)m = (math.pi / 180.0) * l0 * tcoset2 = e2 * pow(tcos, 2)et3 = e2 * pow(tsin, 2)X = 111132.9558 * latitude - 16038.6496 * math.sin(2 * latitude2Rad) + 16.8607 * math.sin(4 * latitude2Rad) - 0.0220 * math.sin(6 * latitude2Rad)N = a / math.sqrt(1 - et3)x = X + N * t * (0.5 * pow(m, 2) + (5.0 - pow(t, 2) + 9.0 * et2 + 4 * pow(et2, 2)) * pow(m, 4) / 24.0 + (61.0 - 58.0 * pow(t, 2) + pow(t, 4)) * pow(m, 6) / 720.0)y = 500000 + N * (m + (1.0 - pow(t, 2) + et2) * pow(m, 3) / 6.0 + (5.0 - 18.0 * pow(t, 2) + pow(t, 4) + 14.0 * et2 - 58.0 * et2 * pow(t, 2)) * pow(m, 5) / 120.0)return x, ydef XY2LatLon(X, Y, L0):iPI = 0.0174532925199433a = 6378137.0f= 0.00335281006247ZoneWide = 3 #按3度带进⾏投影ProjNo = int(X / 1000000)L0 = L0 * iPIX0 = ProjNo * 1000000 + 500000Y0 = 0xval = X - X0yval = Y - Y0e2 = 2 * f - f * f #第⼀偏⼼率平⽅e1 = (1.0 - math.sqrt(1 - e2)) / (1.0 + math.sqrt(1 - e2))ee = e2 / (1 - e2) #第⼆偏⼼率平⽅M = yvalu = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256))fai = u \+ (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * math.sin(2 * u) \+ (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * math.sin(4 * u) \+ (151 * e1 * e1 * e1 / 96) * math.sin(6 * u)\+ (1097 * e1 * e1 * e1 * e1 / 512) * math.sin(8 * u)C = ee * math.cos(fai) * math.cos(fai)T = math.tan(fai) * math.tan(fai)NN = a / math.sqrt(1.0 - e2 * math.sin(fai) * math.sin(fai))R = a * (1 - e2) / math.sqrt((1 - e2 * math.sin(fai) * math.sin(fai)) * (1 - e2 * math.sin(fai) * math.sin(fai)) * (1 - e2 * math.sin(fai) * math.sin(fai)))D = xval / NN#计算经纬度（弧度单位的经纬度）longitude1 = L0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D * D * D * D * D / 120) / math.cos(fai)latitude1 = fai - (NN * math.tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24 + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720)#换换为deglongitude = longitude1 / iPIlatitude = latitude1 / iPIreturn latitude, longitude## print LatLon2XY(40.07837722329, 116.23514827596)# print XY2LatLon(434760.7611718801, 4438512.040474475, 117.0)以上这篇python实现⾼斯投影正反算⽅式就是⼩编分享给⼤家的全部内容了，希望能给⼤家⼀个参考，也希望⼤家多多⽀持。
高斯投影坐标正反算.ppt
2
4
6
8
高斯投影坐标正算（3）
dm0 = dX
dB =M
N cos B =N cos B
,
dq dB dq
M
c
m1 = N cos B
= V
cos B
子午线曲率半径
等量纬度定义式
m2
N 2
sin B cos B
m3
N b
cos3 B(1 t 2
2)
m4
N sin B cos3 24
y 2 dn4 dx
y4
N
cos M
B
(n1
3n3
y2
5n3
y4
)
2n2 y 4n4 y3
N cos B ( dn1 y dn3 y3 dn5 y5 )
M
dx
dx
dx
由恒等式两边对应系数相等，从而得待定系数的递推公式
n1
M N cos B
• 高斯平面直角坐标系： 区分为：自然坐标；国家统一坐标。（掌握两者的换算）
§4.9.2 正形投影的一般条件
一、长度比的通用公式推导
dS 2 (MdB)2 ( N cos Bdl )2
ds2 dx2 dy2
长度比平方为：
m2
ds dS
2
dx2 dy2 (MdB)2 (N cos
dn2 dx
y2 dn4 dx
y4
N
cos M
B
(n1
(完整版)高斯投影正反算
高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程学号：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 ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩L L 经过计算可以得出232244524632235242225sin cos sin cos (594)224sin cos (6158)720cos cos (1) 6cos (5181458)120N N x X B B l B B t l N B B t t l N y N B l B t l N B t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。
「高斯投影坐标正反算公式及适合电算的高斯投影公式」
「高斯投影坐标正反算公式及适合电算的高斯投影公式」高斯投影坐标正反算公式是用于计算高斯投影坐标的数学公式。
高斯投影坐标是一种地理坐标系统，常用于测量和测绘工作中。
高斯投影坐标正算是指已知一个点的经纬度坐标，通过公式计算出该点的高斯投影坐标。
而高斯投影坐标反算是指已知一个点的高斯投影坐标，通过公式计算出该点的经纬度坐标。
一、高斯投影坐标正算公式：已知一个点的经纬度坐标(φ,λ)，其中φ为纬度，λ为经度，以及椭球体参数a、f和中央经线经度L0，可以通过以下步骤计算出该点的高斯投影坐标(X,Y)：1.计算扁率f'：f'=(a-b)/a其中，b=a*(1-f)是椭球体的短半轴。
2.计算黄赤交角ε：ε = atan(b / a)3.计算辅助量t：t = tan(π/4 - φ/2) / [(1 - f' * sin²φ)⁰.⁵ * (1 + e' *sinφ)⁰.⁵]其中，e'=f'*(2-f')是椭球体的第一偏心率。
4.计算辅助量η：η = e'^2 * cos²φ5.计算系数A、B、C和D：A = (L - L0) * cosφC = (L - L0) * cos⁵φ * (5 - tan²φ + 9e'^² + 4e'^⁴ - 24e'^² * tan²φ - 45e'^⁴ * tan²φ)D = (L - L0) * cos⁷φ * (61 - 58tan²φ + tan⁴φ + 270e'^² - 330e'^² * tan²φ)6.计算高斯坐标X和Y：X=k0*a*(A+B/2+C/4+D/6)Y=k0*a*(C/2+D/8)其中，k0是比例系数，一般情况下取1二、高斯投影坐标反算公式：已知一个点的高斯投影坐标(X,Y)，以及椭球体参数a、f、中央经线经度L0、比例系数k0和起始经度L1，可以通过以下步骤计算出该点的经纬度坐标(φ,λ)：1.计算扁率f'：f'=(a-b)/a其中，b=a*(1-f)是椭球体的短半轴。
高斯投影正反算公式83
§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。
8.3.1高斯投影坐标正算公式： B, x,yl ⇒高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。
由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为的偶函数，y 为的奇函数；，即，l l 0330'≤l 20/1/≈''''ρl 如展开为的级数，收敛。
l （8-33）+++=++++=553316644220l m l m l m y l m l m l m m x 式中是待定系数，它们都是纬度B 的函数。
,,10m m 由第三个条件知：qyl x l y q x ∂∂-=∂∂∂∂=∂∂,(8-33)式分别对和q 求偏导数并代入上式l (8-34)----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m 上两式两边相等，其必要充分条件是同次幂前的系数应相等，即l(8-35)dq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)是一种递推公式，只要确定了就可依次确定其余各系数。
0m 由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当时有：0=l(8-36)0m X x==顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2====得：(8-37,38) B Vc B N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-39)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-=依次求得并代入(8-33)式，得到高斯投影正算公式6543,,,m m m m6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ (8-42)5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ8.3.2高斯投影坐标反算公式x,y B,⇒l投影方程：(8-43)),(),(21y x l y x B ϕϕ==满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。
高斯投影正反算编程一.高斯投影正反算基本公式
高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。
二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。
（2）相关流程图1）正算2）反算三．编程的相关代码（1）正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]){FILE *r=fopen("a.txt","r");assert(r!=NULL);FILE *w=fopen("b.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系：北京54=1，西安80=2，WGS84=3：");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一，第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;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;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B))*100.0)*60.0+( g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L))*100.0)*60.0+(g [i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i].L0))*100.0)*60. 0+(g[i].L0*100-(int)(g[i].L0*100))*100.0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽：3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n*n *n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/720;printf("\n所选坐标系的转换结果：x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}（2）反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量：x,y坐标和中央子午线经度// main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L0)!=EOF)//文件为空，无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径，子午圈曲率半径//double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球：1=克拉索夫斯基椭球，2=1975国际椭球，3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球，求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B*C)-g[i] .y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28 *B+24*B*B+6*C+8*B*C);//利用公式求解经纬度//int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果：%d度%d分%lf秒%d 度%d分%lf秒\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;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;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*sin(8*B1 )/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2; }。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
//servey.h
// #ifndef SERVEY_H
// #define SERVEY_H
#include <math.h>
#include <stdio.h>
#include <iomanip.h>
const double PI = 3.149323846;
const double epsilon = 0.;
//角度（度、分、秒）化弧度(带符号）
double angle_to_radian (double alfa)
{
double alfa1,alfa2,fsign,fbeta;
if( fabs(alfa) < epsilon ) return(0.0);
fbeta=fabs(alfa); fsign=alfa/fbeta;
alfa1=floor(fbeta+epsilon)+floor((fbeta-floor(fbeta+epsilon))*100.+epsilon)/60.;
alfa2=(fbeta*100.-floor(fbeta*100.+epsilon))/36.;
alfa1+=alfa2; alfa1=fsign*alfa1*PI/180.;
return (alfa1);
}
//度分秒化为度
double angle_to_degree(double alfa)
{
double alfa_sign; //alfa的正负号
if(alfa>=0)
{
alfa_sign = 1;
}else
{
alfa_sign = -1;
}
alfa = fabs(alfa);
double alfa1,alfa2;
double A = floor(alfa+epsilon);
double B = floor((alfa-A)*100+epsilon);
alfa1 = A+B/60;
alfa2=(alfa*100.-floor(alfa*100.+epsilon))/36.;
alfa1+=alfa2;
return (alfa_sign*alfa1);
}
//弧度化度分秒
double radian_to_angle(double alfa)
{
double alfa_sign; //alfa的正负号
if(alfa>=0)
{
alfa_sign = 1;
}else
{
alfa_sign = -1;
}
alfa = fabs(alfa);
double alfa1=alfa*180./PI+epsilon;
double alfa2=floor(alfa1+epsilon)+floor((alfa1-floor(alfa1+epsilon))*60+epsilon)/100;
double alfa3=(alfa1-angle_to_degree(alfa2))*0.36;
alfa2+=alfa3;
return (alfa_sign*alfa2);
}
//高斯正算函数
//Ellipse为枚举类型
//Cent_Meridian为中央子午线经度
//x,y为返回的高斯平面坐标
int BL_xy(int ellipse, double Cent_Meridian, double B,
double L, double& x, double& y)
{
double a, b, e1, e2; //e1为第一偏心率，e2为第二偏心率
switch (ellipse) //选椭球
{
case 0 : //54椭球
{
a=6378245.; //长半轴
b=6356863.0187730473; //短半轴
e1=sqrt(pow(a,2.)-pow(b,2.))/a;
e2=sqrt(pow(a,2.)-pow(b,2.))/b;
break;。