测量专业基本计算与编程(6)

参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 40.02356784 130.10122627 1.49432980 80000.0045 40.45479027 130.12011040 181.50538463

3 4 40.02356784 115.10000000 36.12010267 414306.5362 43.00558784 118.10030000 218.11267964

5 6 68.58000000 33.05000000 339.49563867 7999606.3891 37.44599755 -122.26000057 9.01078088

7 8 35.00002200 90.00001100 100.00003279 15000000.3380 -30.29209640 215.59043380 290.32533887

3） 贝塞尔算法计算程序

Sub DDZT_Bessel (B1 As Double, L1 As Double, A1 As Double, S12 As Double, _ B2 As Double,L2 As Double, A2 As Double, k As Integer) eps = e2 / (1# - e2) b = a / Sqr(1# + eps) Select Case k Case 1

U1 = Atn(Sqr(1 - e2) * Tan(B1))

sinA0 = Cos(U1) * Sin(A1): cosA0 = Sqr(1 - sinA0 * sinA0) sigma1 = Atn(Tan(U1) / Cos(A1)) xk2 = eps * cosA0 * cosA0 xk4 = xk2 *xk2 xk6 = xk4 * xk2

alpha = (1 - xk2 / 4 + 7 * xk4 / 64 - 15 * xk6 / 256) / b beta = xk2 / 4 - xk4 / 8 + 37 * xk6 / 512 gamma = xk4 / 128 - xk6 / 128 sigma = alpha * S12 Do

sigma0 = sigma

sigma = alpha * S12 + beta * Sin(sigma0) * Cos(2 * sigma1 + sigma0) sigma = sigma + gamma * Sin(2 * sigma0) * Cos(4 * sigma1 + 2 * sigma0) Dsigma = Abs(sigma - sigma0) * P0 '常数P0=206265 Loop While Dsigma > 0.0001 '计算反方位角A2

sinA2 = Cos(U1) * Sin(A1)

cosA2 = Cos(U1) * Cos(sigma) * Cos(A1) - Sin(U1) * Sin(sigma) tanA2 = sinA2 / cosA2: A2 = Abs(Atn(sinA2 / cosA2)) SinA1 = Sin(A1)

If SinA1 < 0 And tanA2 > 0 Then A2 = A2 If SinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If SinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If SinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 '计算大地纬度B2

sinU2 = Sin(U1) * Cos(sigma) + Cos(U1) * Cos(A1) * Sin(sigma) B2 = Atn(sinU2 / Sqr(1 - e2) / Sqr(1 - sinU2 * sinU2)) '计算大地经度L2

sinl = Sin(A1) * Sin(sigma)

cosl = Cos(U1) * Cos(sigma) - Sin(U1) * Sin(sigma) * Cos(A1) tanlambda = sinl / cosl: lambda = Abs(Atn(sinl / cosl)) If tanlambda > 0 And SinA1 > 0 Then lambda = lambda If tanlambda < 0 And SinA1 > 0 Then lambda = pi - lambda If tanlambda < 0 And SinA1 < 0 Then lambda = -lambda If tanlambda > 0 And SinA1 < 0 Then lambda = lambda - pi e4 = e2 * e2 e6 = e4 * e2

xk2 = e2 * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk4 * xk2

alpha1 = (e2 / 2 + e4 / 8 + e6 / 16) - e2 * (1 + e2) * xk2 / 16 + 3 * xk4 * e2 / 128 beta1 = e2 * (1 + e2) * xk2 / 16 - e2 * xk4 / 32 gamma1 = e2 * xk4 / 256

xx = alpha1 * sigma + beta1 * Sin(sigma) * Cos(2 * sigma1 + sigma) xx = xx + gamma1 * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma) l = lambda - sinA0 * xx L2 = L1 + l Case 2

U1 = Atn(Sqr(1 - e2) * Tan(B1)): U2 = Atn(Sqr(1 - e2) * Tan(B2)) DL = L2 - L1

sa1 = Sin(U1) * Sin(U2): sa2 = Cos(U1) * Cos(U2) cb1 = Cos(U1) * Sin(U2): cb2 = Sin(U1) * Cos(U2) lambda = DL Do

lambda0 = lambda

p = Cos(U2) * Sin(lambda0): q = cb1 - cb2 * Cos(lambda0) '计算方位角A1 A1 = Abs(Atn(p / q))

If p > 0 And q > 0 Then A1 = A1 If p > 0 And q < 0 Then A1 = pi - A1 If p < 0 And q < 0 Then A1 = pi + A1 If p < 0 And q > 0 Then A1 = 2 * pi - A1 '计算sigma

Ssigma = p * Sin(A1) + q * Cos(A1) Csigma = sa1 + sa2 * Cos(lambda0) sigma = Abs(Atn(Ssigma / Csigma)) If Csigma > 0 Then sigma = sigma If Csigma < 0 Then sigma = pi - sigma '计算A0 与 sigma1

sinA0 = Cos(U1) * Sin(A1)

sigma1 = Atn(Tan(U1) / Cos(A1)) '计算椭球面经差lambda

cosA0 = Sqr(1 - sinA0 * sinA0) e4 = e2 * e2 e6 = e4 * e2

xk2 = e2 * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk4 * xk2

alpha1 = (e2 / 2 + e4 / 8 + e6 / 16) - e2 * (1 + e2) * xk2 / 16 + _ 3 * xk4 * e2 / 128

beta1 = e2 * (1 + e2) * xk2 / 16 - e2 * xk4 / 32 gamma1 = e2 * xk4 / 256

xx = alpha1 * sigma + beta1 * Sin(sigma) * Cos(2 * sigma1 + sigma) xx = xx + gamma1 * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma) lambda = DL + sinA0 * xx

Dlambda = Abs(lambda - lambda0) * P0 'P0=206265 Loop While Dlambda > 0.0001 '计算椭球面距离S12

cosA0 = Sqr(1 - sinA0 * sinA0) xk2 = eps * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk2 * xk4

alpha = (1 - xk2 / 4 + 7 * xk4 / 64 - 15 * xk6 / 256) / b beta = xk2 / 4 - xk4 / 8 + 37 * xk6 / 512 gamma = xk4 / 128 - xk6 / 128

xs12 = gamma * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma)

S12 = (sigma - beta * Sin(sigma) * Cos(2 * sigma1 + sigma) - xs12) / alpha '计算椭球面反方位角A21 sinA2 = Cos(U1) * Sin(A1)

cosA2 = Cos(U1) * Cos(sigma) * Cos(A1) - Sin(U1) * Sin(sigma) tanA2 = sinA2 / cosA2: A2 = Abs(Atn(sinA2 / cosA2)) SinA1 = Sin(A1)

If SinA1 < 0 And tanA2 > 0 Then A2 = A2 If SinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If SinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If SinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 End Select End Sub

4.3 嵌套系数法解算任意距离的大地主题

目前大地主题解算方法虽然比较完善，但大多数解算方法（特别是反算）仍然存在一些缺陷，主要表现在应用范围受到一定的限制。受其计算精度的影响，它们几乎都只适合一定范围而不能解任意大地距离的主题问题。另外产生奇异性，当两大地点位于同一子午圈或同

在赤道时算法不能正常计算。电脑的广泛应用为改进大地主 …… 此处隐藏：3259字，全部文档内容请下载后查看。喜欢就下载吧 ……