正文
java高斯正反算代码 高斯正反算软件
提示：扫一扫查出行【扫一扫了解最新限行尾号】
复制提示
《大地测量学基础》里面的高斯投影正反算公式及换带计算的VB或者C语言编程
java高斯正反算代码我有VB的java高斯正反算代码，自己很多年前写的，一直用，但是正算－反算－正算后，Y坐标与原来的差了0.5－0.7mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！代码如下：
'高斯坐标正算
Private Sub DadiZs()
Dim t As Double, Itp As Double, X0 As Double, N As Double, L0 As Double
Dim V As Double, ll As Double, W As Double, M As Double
Lat = Radian(Lat)
Lon = Radian(Lon)
L0 = Radian(Lo)
If Tq = 0 Then
a = 6378245 '54椭球参数
b = 6356863.01877305
ep = 0.006693421622966
ep1 = 0.006738525414683
f = (a - b) / a
c = a ^ 2 / b
d = b ^ 2 / a
X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7799 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6973 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)
'X0 = 111134.8611 * (Lat * 180# / Pi) - (32005.7798 * Sin(Lat) + 133.9238 * (Sin(Lat)) ^ 3 + 0.6972 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)
Else
a = 6378140 '75椭球参数
b = 6356755.28815753
ep = 0.006694384999588
ep1 = 0.006739501819473
f = (a - b) / a
c = a ^ 2 / b
d = b ^ 2 / a
X0 = 111133.0047 * (Lat * 180 / Pi) - (32009.8575 * Sin(Lat) + 133.9602 * (Sin(Lat)) ^ 3 + 0.6976 * (Sin(Lat)) ^ 5 + 0.0039 * (Sin(Lat)) ^ 7) * Cos(Lat)
End If
ll = Lon - L0
t = Tan(Lat)
Itp = ep1 * Cos(Lat) ^ 2
W = Sqr(1 - ep * Sin(Lat) ^ 2)
V = Sqr(1 + ep1 * Cos(Lat) ^ 2)
M = c / V ^ 3
N = a / W
'x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp + 4 * Itp * Itp) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp - 330 * t ^ 2 * Itp) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320
x = X0 + N * t * (Cos(Lat)) ^ 2 * ll ^ 2 / 2 + N * t * (5 - t * t + 9 * Itp ^ 2 + 4 * Itp ^ 4) * (Cos(Lat)) ^ 4 * ll ^ 4 / 24 + N * t * (61 - 58 * t ^ 2 + t ^ 4 + 270 * Itp ^ 2 - 330 * t ^ 2 * Itp ^ 2) * (Cos(Lat)) ^ 6 * ll ^ 6 / 720 + N * t * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 8 * ll ^ 8 / 40320
y = N * Cos(Lat) * ll + N * (1 - t * t + Itp) * (Cos(Lat)) ^ 3 * ll ^ 3 / 6 + N * (5 - 18 * t * t + t ^ 4 + 14 * Itp - 58 * Itp * t * t) * (Cos(Lat)) ^ 5 * ll ^ 5 / 120 + N * (61 - 479 * t ^ 2 + 179 * t ^ 4 - t ^ 6) * Cos(Lat) ^ 7 * ll ^ 7 / 5040
r = Sin(Lat) * ll + Sin(Lat) * (Cos(Lat)) ^ 2 * ll ^ 3 * (1 + 3 * Itp + 2 * Itp ^ 2) / 3 + Sin(Lat) * (Cos(Lat)) ^ 4 * ll ^ 5 * (2 - t * t) / 15
r = Degree(r)
y = y + 500000#
End Sub
'高斯反算
Private Sub DadiFs()
Dim t As Double, Itp As Double, X0 As Double, Bf As Double, N As Double
Dim v As Double, ll As Double, W As Double, M As Double, L0 As Double
L0 = Radian(Lo)
X0 = x * 0.000001
y = y - 500000#
If Tq = 0 Then
a = 6378245 '54椭球参数
b = 6356863.01877305
ep = 0.006693421622966
ep1 = 0.006738525414683
f = (a - b) / a
c = a ^ 2 / b
d = b ^ 2 / a
If X0 3 Then
Bf = 9.04353301294 * X0 - 0.00000049604 * X0 ^ 2 - 0.00075310733 * X0 ^ 3 - 0.00000084307 * X0 ^ 4 - 0.00000426055 * X0 ^ 5 - 0.00000010148 * X0 ^ 6
ElseIf X0 6 Then
Bf = 27.11115372595 + 9.02468257083 * (X0 - 3) - 0.00579740442 * (X0 - 3) ^ 2 - 0.00043532572 * (X0 - 3) ^ 3 + 0.00004857285 * (X0 - 3) ^ 4 + 0.00000215727 * (X0 - 3) ^ 5 - 0.00000019399 * (X0 - 3) ^ 6
End If
Else
a = 6378140 '75椭球参数
b = 6356755.28815753
ep = 0.006694384999588
ep1 = 0.006739501819473
f = (a - b) / a
c = a ^ 2 / b
d = b ^ 2 / a
If X0 3 Then
Bf = 9.04369066313 * X0 - 0.00000049618 * X0 ^ 2 - 0.00075325505 * X0 ^ 3 - 0.0000008433 * X0 ^ 4 - 0.00000426157 * X0 ^ 5 - 0.0000001015 * X0 ^ 6
ElseIf X0 6 Then
Bf = 27.11162289465 + 9.02483657729 * (X0 - 3) - 0.00579850656 * (X0 - 3) ^ 2 - 0.00043540029 * (X0 - 3) ^ 3 + 0.00004858357 * (X0 - 3) ^ 4 + 0.00000215769 * (X0 - 3) ^ 5 - 0.00000019404 * (X0 - 3) ^ 6
End If
End If
Bf = Bf * Pi / 180#
t = Tan(Bf)
Itp = ep1 * Cos(Bf) ^ 2
W = Sqr(1 - ep * Sin(Bf) ^ 2)
v = Sqr(1 + ep1 * Cos(Bf) ^ 2)
M = c / v ^ 3
N = a / W
Lat = Bf - 0.5 * v ^ 2 * t * ((y / N) ^ 2 - (5 + 3 * t * t + Itp - 9 * Itp * t * t) * (y / N) ^ 4 / 12 + (61 + 90 * t * t + 45 * t ^ 4) * (y / N) ^ 6 / 360)
ll = ((y / N) - (1 + 2 * t * t + Itp) * (y / N) ^ 3 / 6 + (5 + 28 * t * t + 24 * t ^ 4 + 6 * Itp + 8 * Itp * t * t) * (y / N) ^ 5 / 120) / Cos(Bf)
r = y * t / N - y ^ 3 * t * (1 + t * t - Itp) / (3 * N ^ 3) + y ^ 5 * t * (2 + 5 * t * t + 3 * t ^ 4) / (15 * N ^ 5)
Lat = Degree(Lat)
Lon = Degree(L0 + ll)
r = Degree(r)
End Sub
有了正反算，换带也就完成了！
用到的子程序：
Public Const Pi = 3.14159265358979, p = 206264.806
Public Cktq As String
'角度化弧度
Public Function Radian(a As Double) As Double
Dim Ro As Double
Dim c As Double
Dim Fs As Double
Dim Ib As Integer
Dim Ic As Integer
If a 0 Then a = -a: t = 1
Ro = Pi / 180#
Ib = Int(a)
c = (a - Ib) * 100#
Ic = Int(c + 0.000000000001)
Fs = (c - Ic) * 100#
If t = 1 Then Radian = -(Ib + Ic / 60# + Fs / 3600#) * Ro Else Radian = (Ib + Ic / 60# + Fs / 3600#) * Ro
End Function
'弧度化角度
Public Function Degree(a As Double) As Double
Dim Bo As Double
Dim Fs As Double
Dim Im As Integer
Dim Id As Integer
If a 0 Then a = -a: t = 1
Bo = a
Call DMS(Bo, Id, Im, Fs)
If t = 1 Then Degree = -(Id + Im / 100# + Fs / 10000#) Else Degree = Id + Im / 100# + Fs / 10000#
End Function
Public Sub DMS(a As Double, Id As Integer, Im As Integer, Fs As Double)
Dim Bo As Double
Dim c As Double
c = a
c = 180# / Pi * c
Id = Int(c)
Bo = (c - Id) * 60
Im = Int(Bo)
Fs = (Bo - Im) * 60
End Sub
'取位计算
Public Function Qw(a As Double, Ws As Integer) As Double
Qw = Int(a * 10 ^ Ws + 0.5) / 10 ^ Ws
End Function
相关问答
Q1: 高斯正算公式
弧长参数计算:
i_a = x1
i_a_ = x2
//if i_a = 6378245.000 and i_a_ = 298.3 then
// i_a0 = 111134.8611
// i_a1 = 32005.7799
// i_a2 = 133.9238
// i_a3 = 0.6973
// i_a4 = 0.0039
//elseif i_a = 6378140.000 and i_a_ = 298.257 then //xi80
// i_a0 = 111133.0047
// i_a1 = 32009.8575
// i_a2 = 133.9602
// i_a3 = 0.6976
// i_a4 = 0.0039
//else
// return -1
//end if
i_b = i_a - i_a/i_a_
i_c = i_a*i_a/i_b
i_e1 = (i_a + i_b)*(i_a - i_b)/i_a/i_a
i_e2 = (i_a + i_b)*(i_a - i_b)/i_b/i_b
i_a0 = 1.0 + (3.0/4.0 + (45.0/64.0 + ( 525.0/768.0 + (33075.0/49152.0 + ( 654885.0/983040.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
i_a1 = (3.0/4.0 + (15.0/16.0 + (1575.0/1536.0 + (6615.0/6144.0 + (1091475.0/983040.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
i_a2 = ((15.0/64.0 + ( 315.0/768.0 + (6615.0/12288.0 + ( 155925.0/245760.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
i_a3 = (((105.0/1536.0 + (945.0/6144.0 + (467775.0/1966080.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
i_a4 = ((((945.0/49152.0 + ( 51975.0/983040.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
i_a5 = ((((( 10395.0/1966080.0)*i_e1)*i_e1)*i_e1)*i_e1)*i_e1
double d
d = i_a*(1.0 - i_e1)
i_a0 = i_a0*d
i_a1 = -i_a1*d/2.0
i_a2 = i_a2*d/4.0
i_a3 = -i_a3*d/6.0
i_a4 = i_a4*d/8.0
i_a5 = -i_a5*d/10.0
---------------------------------------------------------------------------------------------
计算弧长:
i_a0* x +
i_a1* sin( 2*x) +
i_a2* sin( 4*x) +
i_a3* sin( 6*x) +
i_a4* sin( 8*x) +
i_a5* sin(10*x)
----------------------------------------------------------------------------
正算公式:
GX_L0 = L0
GX_L = L - L0
GX_B = B
i_tanB = tan(gx_B)
i_tanBB = i_tanB*i_tanB
i_cosB = cos(gx_B)
i_cosBBLL = i_cosB*i_cosB*GX_L*GX_L
i_nn = i_e2*i_cosB*i_cosB
i_N = i_c/sqrt(1.0 + i_nn)
GX_X =getarc(GX_B) + i_N*i_tanB*i_cosBBLL*(0.5 +
i_cosBBLL*((5.0 - i_tanBB +9.0*i_nn + 4.0*i_nn*i_nn)/24.0 +
i_cosBBLL*((61.0 -58.0*i_tanBB + i_tanBB*i_tanBB)/720.0)))
GX_Y = i_N*i_cosB*GX_L*(1.0 + i_cosBBLL*((i_nn - i_tanBB + 1.0)/6.0+
i_cosBBLL*((5.0 + (i_tanBB - 18.0)*i_tanBB+
(14.0 -58.0*i_tanBB)*i_nn)/120.0)))
;ID=55396
Q2: java生成随机数0或1，只能是0，1二者其一，最好可以使1或0出现的概率符合高斯分布（正态分布）
Java内置java高斯正反算代码了高斯正态随机数
Random.nextGaussian()获得均值为0java高斯正反算代码，标准差为1java高斯正反算代码的正态分布。
Random r=new Random();
for(int i=0;i20;i++)
System.out.println(r.nextGaussian());
=======
-0.13281412372658716
-0.01315269766483013
-0.056744105795208055
-1.0402491360478257
-0.5311391585982737
-0.3630802799191039
-0.9866340328819281
-0.8579720715397737
-0.9265915705393113
-0.21839994469893823
0.3671093015925613
-0.6560992663799461
-1.2695699508671745
0.046509998380311696
0.1952527864179513
0.9997730342084369
如果求阀值0.7java高斯正反算代码，
Random r=new Random();
int c0=0,c1=0;
for(int i=0;i10000;i++){
byte n=(byte) (r.nextInt(10)=7?1:0); //获得一个0或1
if(n==0)c0++; else c1++;
}
System.out.println("zero:"+(c0/10000.0)+" one:"+(c1/10000.0));
=========
打印概率
zero:0.7068 one:0.2932
Q3: 高斯投影正反算
//高斯投影正、反算
//////6度带宽 54年北京坐标系
//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
void GaussProjCal(double longitude, double latitude, double *X, double *Y)
{
int ProjNo=0; int ZoneWide; ////带宽
double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
double a,f, e2,ee, NN, T,C,A, M, iPI;
iPI = 0.0174532925199433; ////3.1415926535898/180.0;
ZoneWide = 6; ////6度带宽
a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
////a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int)(longitude / ZoneWide) ;
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ;
latitude0=0;
longitude1 = longitude * iPI ; //经度转换为弧度
latitude1 = latitude * iPI ; //纬度转换为弧度
e2=2*f-f*f;
ee=e2*(1.0-e2);
NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));
T=tan(latitude1)*tan(latitude1);
C=ee*cos(latitude1)*cos(latitude1);
A=(longitude1-longitude0)*cos(latitude1);
M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2
*e2/1024)*sin(2*latitude1)
+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*sin(6*l
atitude1));
xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
X0 = 1000000L*(ProjNo+1)+500000L;
Y0 = 0;
xval = xval+X0; yval = yval+Y0;
*X = xval;
*Y = yval;
}
//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)
void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) 字串9
{
int ProjNo; int ZoneWide; ////带宽
double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
iPI = 0.0174532925199433; ////3.1415926535898/180.0;
a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数
////a=6378140.0; f=1/298.257; //80年西安坐标系参数
ZoneWide = 6; ////6度带宽
ProjNo = (int)(X/1000000L) ; //查找带号
longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI ; //中央经线
X0 = ProjNo*1000000L+500000L;
Y0 = 0;
xval = X-X0; yval = Y-Y0; //带内大地坐标
e2 = 2*f-f*f;
e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));
ee = e2/(1-e2);
M = yval;
u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(
4*u)
+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);
C = ee*cos(fai)*cos(fai);
T = tan(fai)*tan(fai);
NN = a/sqrt(1.0-e2*sin(fai)*sin(fai)); 字串1
R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin
(fai)*sin(fai)));
D = xval/NN;
//计算经度(Longitude) 纬度(Latitude)
longitude1 = longitude0+(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)/cos(fai);
latitude1 = fai -(NN*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);
//转换为度 DD
*longitude = longitude1 / iPI;
*latitude = latitude1 / iPI;
}
NN卯酉圈曲率半径，测量学里面用N表示
M为子午线弧长，测量学里用大X表示 字串2
fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示 字串4
R为底点所对的曲率半径，测量学里用Nf表示
Q4: 高斯投影正反算的直接计算结果是什么
高斯投影坐标正反算
一、基本思想：
高斯投影正算公式就是由大地坐标（L ，B ）求解高斯平面坐标（x ，y ），而高斯投影反算公式则是由高斯平面坐标（x ，y)求解大地坐标（L ，B).
二、计算模型:
基本椭球参数：
椭球长半轴a
椭球扁率f
椭球短半轴:(1)b a f =-
椭球第一偏心率
第 1 页
高斯计_高精度_高效率_高品质
全面覆盖电离辐射， 非电离辐射及直流磁场的检测。低频电磁场还可在线监测， 多种功能供用户选配
点击立即咨询，了解更多详情
咨询
深圳市柯雷科技开发 广告

：e = 椭球第二偏心率

：e '=高斯投影正算公式：此公式换算的精度为0.001m
64256
442234
22)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 ''+-''+''++-''+''⋅''+
=ρηηρρ 52224255
32233
)5814185(cos 120)1(cos 6cos l t t t B N l t B
第 2 页
N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ
其中：角度都为弧度
B 为点的纬度,0l L L ''=-，L 为点的经度,0L 为中央子午线经度；
N 为子午圈曲率半径,1222
(1sin )N a e B -=-；
tan t B =； 222cos e B η'=
180
3600ρπ''=*
其中X 为子午线弧长：
2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ⎡⎤=--++-+⎢⎥⎣⎦
第 3 页
02468,,,,a a a a a 为基本常量，按如下公式计算：
2004682426844686868
83535281612815722321637816323216128m 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;;268
m a e m e m m e m m e m m e m =-====;
高斯投影反算公式：此公式换算的精度为0.0001’'。
()()()()22222432465
第 4 页
3
2235
2422250
53922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f f
f f f f f f f
t t B B y t t y
M N M N t y t t y
M N y y l t N B N B y t t t N B L l L ηηηηη=-
+++--++=-+++++++=+
其中: 0L 为中央子午线经度。
第 5 页
f B 为底点纬度，也就是当x X =时的子午线弧长所对应的纬度。按照子午线弧长公式：68240sin 2sin 4sin 6sin82468
a a a a X a B B B B B =-+-+，迭代进行计算； 初始开始时设：10f B X a =
以后每次迭代按下式计算:
10
6824(())()sin 2sin 4sin 6sin82468i
i f f i
i i i i f
f f f f B X X F B a a a a a F B B B B B +=-=-+-+
重复迭代至1i
第 6 页
i
f f B B ε+-；为止。
1222
(1sin )f f N a e B -=-;
3
2222(1)(1sin )f f M a e e B -=-- tan f f t B =；
222cos f
f e B η'=
海福特椭球（1910) 我国52年以前基准椭球 a=6378388m b=6356911。9461279m α=0.33670033670
克拉索夫斯基椭球（1940 Krassovsky)
第 7 页
北京54坐标系基准椭球 a=6378245m b=6356863。018773m α=0.33523298692
1975年I 。U 。G.G 推荐椭球（国际大地测量协会1975) 西安80坐标系基准椭球
a=6378140m b=6356755.2881575m α=0.0033528131778
WGS-84椭球（GPS 全球定位系统椭球、17届国际大地测量协会） WGS —84 GPS 基准椭球
a=6378137m b=6356752.3142451m α=0。00335281006247
三、程序代码函数:
/*＊＊**＊＊＊**＊＊高斯投影正算函数***＊＊**＊＊**＊*＊*
第 8 页
输入 : double a ，f 椭球参数，B,L 为大地坐标,L0为中央子午线的经度，单位为弧度，x,y 为高斯平面坐标,y 加上了500000常量
返回:none
＊＊＊**＊*＊*＊*＊＊＊＊*＊＊＊＊*＊*＊*＊*＊**＊＊*＊＊*＊＊*＊＊＊/
void gaosiforward (double a ,double f ，double B ，double L ，double L0，double ＆x ，double ＆y ） ｛
double b , c ，e1， e2； //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率
double l ， W ,N ， M ， daihao ；//W 为常用辅助函数，N 为子午圈曲率半径，M 为卯酉圈曲率半径
第 9 页
double X ;//子午线弧长，高斯投影的坐标
double ruo ， ita , sb , cb ，t ;
double m [5］，n ［5]；
//计算一些基本常量
{
b =a ＊(1-f )；
e1=sqrt (a ＊a -b ＊b )/a ;
e2=sqrt （a *a -b ＊b ）/b ;
c =a ＊a /b ；
m ［0］=a *(1-e1＊e1）; m ［1］=3＊（e1＊e1*m [0])/2。0;
m[2］=5*(e1＊e1*m［1]）/4。0;
第 10 页
m[3]=7*（e1*e1*m［2]）/6.0;
m[4］=9*(e1*e1*m[3]）/8。0；
n[0]=m[0］+m［1]/2+3*m［2］/8+5＊m［3]/16+35*m[4]/128;
n[1］=m［1］/2+m[2］/2+15*m[3]/32+7＊m[4]/16；
n[2］=m［2]/8+3＊m[3]/16+7*m［4]/32;
n［3]=m［3]/32+m[4]/16；
n［4］=m［4]/128; /////by kjh 2014。5。22 把改成了
｝
//由纬度计算子午线弧长
第 11 页
{
X=n［0]*B—sin（B)*cos(B）＊（（n［1]-n[2］+n[3]）+（2*n［2]-（16*n[3］/3.0))＊sin(B)＊sin(B）+16*n[3]*pow(sin(B)，4)/3。0);
}
l=L—L0；//弧度
ita=e2＊cos（B);
sb=sin(B）；
cb=cos(B）；
W=sqrt（1-e1*e1*sb*sb）；
N=a/W；
第 12 页
t=tan(B);
ruo=（180/Pi)＊3600;
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）；
y=y+500000;
}
/＊＊＊＊＊＊*＊*＊＊*＊*高斯反算函数***＊****＊＊*****
第 13 页
输入： double a ，f 椭球参数， x，y为高斯平面坐标，L0为中央子午线的经度； B,L为大地坐标，单位为弧度
＊返回:none
＊*＊＊＊＊＊＊*＊＊＊＊**＊*＊*＊**＊＊＊*＊＊＊/
void gaosibackward(double a,double f,double x，double y,double L0,double＆B,double L)
｛
double b, c,e1, e2; //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率
double Bf,itaf,tf,Nf，Mf,Wf；
double l;
第 14 页
double m[5］,n[5］;
y=y-500000;
//计算一些基本常量
｛
b=a＊（1—f）；
e1=sqrt(a*a—b*b）/a;
e2=sqrt(a*a-b＊b)/b;
c=a*a/b;
m［0］=a＊（1-e1＊e1);
m[1]=3*（e1*e1*m［0］)/2.0；
m[2]=5*（e1＊e1*m［1］)/4。0；
第 15 页
m[3］=7*(e1＊e1*m[2］）/6。0；
m[4]=9*(e1＊e1＊m［3］)/8.0；
n［0]=m［0］+m［1]/2+3＊m［2]/8+5*m［3]/16+35*m［4］/128；
n[1]=m［1］/2+m[2]/2+15＊m［3]/32+7*m[4］/16;
n[2］=m［2］/8+3*m［3］/16+7＊m[4]/32；
n[3]=m［3]/32+m[4］/16;
n［4]=m［4]/128；
｝
//计算Bf
第 16 页
｛
double Bf1,Bfi0，Bfi1，FBfi；
Bf1=x/n[0];
Bfi0=Bf1;
Bfi1=0;
FBfi=0；
int num=0；
do
｛
num=0;
FBfi=0.0-n[1］＊sin（2＊Bfi0）/2.0+n[2］*s
第 17 页
in(4*Bfi0）/4。0-n［3]*sin(6*Bfi0）/6。0；
Bfi1=（x—FBfi)/n［0];
if (fabs（Bfi1—Bfi0）(Pi*pow(10.0,-8）/（36*18）)）
{
num=1；
Bfi0=Bfi1；
｝
｝while（num==1）;
Bf=Bfi1；
}
tf=tan（Bf）；
第 18 页
Wf=sqrt（1-e1*e1＊sin（Bf）＊sin（Bf)）;
Nf=a/Wf；
Mf=a＊（1—e1*e1)/(Wf＊Wf＊Wf)；
itaf=e2*cos（Bf);
B=Bf-tf＊y＊y/（2*Mf＊Nf）+tf*(5+3*tf＊tf+itaf＊itaf—9＊itaf*itaf*tf＊tf）*pow（y，4)/(24＊Mf＊pow(Nf，3））—tf＊（61+90*tf*tf+45*pow（tf,4）)*pow（y,6）/（720*Mf*pow（Nf，5)）;
l=y/（Nf*cos(Bf））-(1+2*tf＊tf+itaf＊itaf)＊pow（y,3)/（6＊pow(Nf,3)＊cos（Bf))+(5+28＊
tf*tf+24＊pow(tf,4）+6*itaf*itaf+8*itaf＊itaf*tf*tf)*pow(y，5）/（120＊pow(Nf，5)*
第 19 页
cos(Bf）);
L=l+L0;
｝
2014-5—22
’输入: double a ，f 椭球参数，B，L为大地坐标，L0为中央子午线的经度，单位为弧度，x，y为高斯平面坐标,y加上了常量
Private Function gaosiforward（ByVal a As Double, ByVal f As Double, ByVal B As Double，ByVal L As Double，ByVal L0 As Double）As Double（)
Dim x， y， xy(2) As Double
Dim bb, c, e1， e2 As Double’短半轴，极点
第 20 页
处的子午线曲率半径，第一偏心率，第二偏心率
Dim ll, W, N， M， daihao As Double’W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径
Dim xx As Double'子午线弧长，高斯投影的坐标
Dim ruo， ita， sb, cb, t As Double
Dim mm（5), nn（5) As Double
bb = a * （1 — f）
e1 = Math。Sqrt(a ＊ a - bb ＊ bb） / a
e2 = Math.Sqrt（a ＊ a — bb ＊ bb） / bb
c = a ＊ a / bb
第 21 页
mm（0) = a ＊（1 - e1 ＊ e1）
mm(1） = 3 ＊（e1 * e1 ＊ mm(0)） / 2.0
mm(2) = 5 * （e1 ＊ e1 ＊ mm（1)) / 4.0
mm（3) = 7 ＊（e1 ＊ e1 * mm(2)） / 6。0
mm（4） = 9 * （e1 ＊ e1 * mm（3）) / 8.0
nn(0) = mm(0） + mm（1) / 2 + 3 ＊ mm（2) / 8 + 5 * mm（3) / 16 + 35 ＊ mm（4） / 128 nn(1) = mm(1) / 2 + mm(2） / 2 + 15 ＊ mm（3) / 32 + 7 ＊ mm（4） / 16
nn（2) = mm(2） / 8 + 3 ＊ mm（3) / 16 + 7 ＊ mm（4) / 32
nn(3) = mm（3) / 32 + mm(4） / 16
nn（4） = mm（4) / 128
第 22 页
xx = nn（0) ＊ B — Sin（B) * Cos（B）＊（（nn(1） - nn(2） + nn(3)) + (2 ＊ nn（2) - (16 ＊ nn(3） / 3。0）） * Sin（B）＊ Sin（B） + 16 * nn(3） * Pow(Sin（B), 4） / 3.0）
ll = L — L0 ’弧度
ita = e2 ＊ Cos（B)
sb = Sin（B)
cb = Cos(B)
W = Sqrt（1 - e1 ＊ e1 * sb * sb）
N = a / W
t = Tan(B）
ruo = (180 / PI） * 3600
第 23 页
x = (xx + N ＊ sb ＊ cb * ll ＊ ll / 2 + N ＊ sb ＊ cb ＊ cb ＊ cb * （5 - t ＊ t + 9 ＊ ita ＊ ita + 4 * ita ＊ ita * ita * ita) * ll ＊ ll * ll * ll / 24 + N * sb ＊ cb ＊ cb * cb ＊ cb ＊ cb * (61 — 58 ＊ t * t + t * t * t * t) * ll * ll ＊ ll * ll ＊ ll * ll / 720)
第 24 页
百度文库
搜索

百度文库10亿海量资料，查找管理一应俱全
打开APP
关于java高斯正反算代码和高斯正反算软件的介绍到此就结束了，不知道你从中找到你需要的信息了吗 ？如果你还想了解更多这方面的信息，记得收藏关注本站。
