测量方位角计算公式VB源代码.docx

1、测量方位角计算公式VB源代码测量方位角计算公式VB源代码角度化弧度Public Function Radian(a As Double) As Double Dim Ra As Double Dim c As Double Dim FS As Double Dim Ib As Integer Dim Ic As Integer Ra = pi / 180# Ib = Int(a) c = (a - Ib) * 100# Ic = Int(c) FS = (c - Ic) * 100# Radian = (Ib + Ic / 60# + FS / 3600#) * RaEnd Function弧

2、度化角度Public Function Degree(a As Double) As Double Dim B As Double Dim Fs1 As Double Dim Im1 As Integer Dim Id1 As Integer B = a Call DMS(B, Id1, Im1, Fs1) Degree = Id1 + Im1 / 100# + Fs1 / 10000#End FunctionPublic Sub DMS(a As Double, ID As Integer, IM As Integer, FS As Double) Dim B As Double Dim c

3、 As Double c = a c = 180# / pi * c ID = Int(c + 0.0000005) B = (c - ID) * 60 + 0.0005 IM = Int(B) FS = (B - IM) * 60End Sub计算两点间的方位角Public Function azimuth(x1 As Double, y1 As Double, x2 As Double, y2 As Double) As Single Dim dx As Double Dim dy As Double Dim fwj As Double dx = x2 - x1 dy = y2 - y1

4、If dy 0 Then fwj = pi * (1 - Sgn(dy) / 2) - Atn(dx / dy) azimuth = Degree(fwj) Else If dx 0 Then azimuth = 0 Else azimuth = 180 End If End IfEnd Function5.2程序字母代表含义La起点里程， R圆曲线半径，l0两端缓和曲线长，曲线转向角，T切线长，L曲线长，E0外矢距，q切曲差，cc线路转向（cc=1时，线路向右转；cc=-1时，线路向左转）,d桩距，m边桩距，Lii点的里程，ZH直缓点，HY缓圆点，JD交点，QZ曲中点，YZ圆直点，YH圆缓点

5、，HZ缓直点，LJD交点里程，LZH直缓点里程，LHY缓圆点里程，LQZ曲中点里程，LYH圆缓点里程，LHZ缓直点里程。5.3程序设计思路及具体程序的编写程序的基本构思：本程序利用已知直、曲线要素、交点坐标，计算出曲线上逐桩坐标，根据逐桩坐标就可实地放样。此程序在不同测站测设同一曲线时，调用程序只需改变测站和后视点坐标即可。不过测站点和后视点的坐标系统与逐桩坐标系统必须一致。程序设计中主要发现的问题及解决办法：5.3.1 起点为非整桩的情况本程序采用将非整桩起点的坐标通过变换成为整桩的方法。直线起点 L1=Int(la / 10) * 10 + k)第一缓和曲线 L2=Int(Ljd T +

6、K)圆曲线 L3=Int(Ljd - T + lo + K )第二缓和曲线 L4=Int(Ljd - T + L lo + K )并计算出这些点的坐标。表5.1 起点是否整桩的判断线型整桩起算点非整桩起算点直线起终点（两端曲线的端点）L1点第一缓和曲线ZH点L2点圆曲线ZH（带缓和曲线）或ZY点L2点第二缓和曲线HZ点HZ点5.3.2 路线方位角的计算程序中采用的起算数据中不包含方位角一项，因此必须要由已知的两点坐标计算出线路的方位角。本程序中利用ZH点和JD的坐标计算。建立一个计算方位角的模块：Function ab#(Xa#, Ya#, Xb#, Yb#)ab = Atn(Abs(Yb -

7、 Ya) / Abs(Xb - Xa) If Yb - Ya = 0 And Xb - Xa = 0 Then ab = ab ElseIf Yb - Ya 0 And Xb - Xa 0 Then ab = 3.141592654 - ab ElseIf Yb - Ya = 0 And Xb - Xa 反算正算后，Y坐标与原来的差了0.50.7mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！代码如下：高斯坐标正算Private Sub DadiZs()Dim t As Double,

8、Itp As Double, X0 As Double, N As Double, L0 As DoubleDim V As Double, ll As Double, W As Double, M As DoubleLat = 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

9、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 = 63

10、78140 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 Ifll = Lon - L0t = Tan(Lat)Itp =

11、 ep1 * Cos(Lat) 2W = Sqr(1 - ep * Sin(Lat) 2)V = Sqr(1 + ep1 * Cos(Lat) 2)M = c / V 3N = a / Wx = 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

12、 * t * (1385 - 3111 * t 2 + 543 * t 4 - t 6) * Cos(Lat) 8 * ll 8 / 40320x = 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

13、* t 2 + 543 * t 4 - t 6) * Cos(Lat) 8 * ll 8 / 40320y = 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 / 5040r = Sin(Lat) * ll + Sin(

14、Lat) * (Cos(Lat) 2 * ll 3 * (1 + 3 * Itp + 2 * Itp 2) / 3 + Sin(Lat) * (Cos(Lat) 4 * ll 5 * (2 - t * t) / 15r = 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 DoubleDim v As Double, ll As Double, W As Double, M As Double, L0 As

15、 DoubleL0 = Radian(Lo)X0 = x * 0.000001y = 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

16、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 IfElse a = 6378140 75椭球参数 b = 6356755.2

17、8815753 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 *

18、(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 IfEnd IfBf = Bf * Pi / 180#t = Tan(Bf)Itp = ep1 * Cos(Bf) 2W = Sqr(1 - ep * Sin(Bf) 2)v = Sqr(1 + ep1 * Cos(Bf) 2)M = c / v 3N = a / WLat = Bf

19、 - 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

20、* (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.806Public 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