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

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

《测量方位角计算公式VB源代码.docx》由会员分享，可在线阅读，更多相关《测量方位角计算公式VB源代码.docx（11页珍藏版）》请在冰豆网上搜索。

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

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

角度化弧度

PublicFunctionRadian（aAsDouble）AsDouble

 DimRaAsDouble

 DimcAsDouble

 DimFSAsDouble

 DimIbAsInteger

 DimIcAsInteger

 Ra=pi/180#

 Ib=Int（a）

 c=（a-Ib）*100#

 Ic=Int（c）

 FS=（c-Ic）*100#

 Radian=（Ib+Ic/60#+FS/3600#）*Ra

EndFunction

'弧度化角度

PublicFunctionDegree（aAsDouble）AsDouble

 DimBAsDouble

 DimFs1AsDouble

 DimIm1AsInteger

 DimId1AsInteger

 B=a

 CallDMS（B,Id1,Im1,Fs1）

 Degree=Id1+Im1/100#+Fs1/10000#

EndFunction

PublicSubDMS（aAsDouble,IDAsInteger,IMAsInteger,FSAsDouble）

 DimBAsDouble

 DimcAsDouble

 c=a

 c=180#/pi*c

 ID=Int（c+0.0000005）

 B=（c-ID）*60+0.0005

 IM=Int（B）

 FS=（B-IM）*60

EndSub

'计算两点间的方位角

PublicFunctionazimuth（x1AsDouble,y1AsDouble,x2AsDouble,y2AsDouble）AsSingle

 DimdxAsDouble

 DimdyAsDouble

 DimfwjAsDouble

 dx=x2-x1

 dy=y2-y1

 Ifdy<>0Then

   fwj=pi*（1-Sgn（dy）/2）-Atn（dx/dy）

   azimuth=Degree（fwj）

 Else

   Ifdx>0Then

     azimuth=0

   Else

     azimuth=180

   EndIf

 EndIf

EndFunction

5.2程序字母代表含义

La—起点里程，R—圆曲线半径，l0—两端缓和曲线长，α—曲线转向角，T—切线长，L—曲线长，E0—外矢距，q—切曲差，cc—线路转向（cc=1时，线路向右转；cc=-1时，线路向左转）,d―桩距，m―边桩距，Li―i点的里程，ZH―直缓点，HY―缓圆点，JD―交点，QZ―曲中点，YZ―圆直点，YH―圆缓点，HZ―缓直点，LJD—交点里程，LZH―直缓点里程，LHY―缓圆点里程，LQZ―曲中点里程，LYH―圆缓点里程，LHZ―缓直点里程。

5.3程序设计思路及具体程序的编写

程序的基本构思：

本程序利用已知直、曲线要素、交点坐标，计算出曲线上逐桩坐标，根据逐桩坐标就可实地放样。

此程序在不同测站测设同一曲线时，调用程序只需改变测站和后视点坐标即可。

不过测站点和后视点的坐标系统与逐桩坐标系统必须一致。

程序设计中主要发现的问题及解决办法：

5.3.1起点为非整桩的情况

本程序采用将非整桩起点的坐标通过变换成为整桩的方法。

直线起点    L1=Int（（la/10）*10+k）

第一缓和曲线L2=Int（Ljd–T+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-Ya）/Abs（Xb-Xa））

 IfYb-Ya>=0AndXb-Xa>=0Then

  αab=αab

 ElseIfYb-Ya>0AndXb-Xa<0Then          

 αab=3.141592654-αab

 ElseIfYb-Ya<=0AndXb-Xa<=0Then

  αab=3.141592654+αab

 Else

 αab=2*3.141592654-αab

EndIf

EndFunction

αab——某点的方位角，Xa——起点X坐标（ZH或ZY点坐标），Ya——起点Y坐标（ZH或ZY点坐标），Xb——终点X坐标（JD的X坐标），Yb——终点Y坐标（JD的Y坐标）

5.3.3计算曲线边桩坐标的程序

Xi=x0+Xxi*Cos（A0）-cc*Yyi*Sin（A0） 'i点的X坐标

Yi=y0+Xxi*Sin（A0）+cc*Yyi*Cos（A0） 'i点的Y坐标

Ai=A0+cc*β

Xz=Xi+m*Cos（Ai-3.141592654/2）   '左边桩X坐标

 Yz=Yi+m*Sin（Ai-3.141592654/2）    '左边桩Y坐标

 Xy=Xi+m*Cos（Ai+3.141592654/2）   '右边桩X坐标

 Yy=Yi+m*Sin（Ai+3.141592654/2）    '右边桩Y坐标

A0计算参考起点的方位角，Ai为计算点的方位角，m为边桩距，Xi为计算点的X坐标，Yi为计算点的Y坐标，cc为转向（左转-1，右转+1），β=li^2/（2*r*lo）。

5.3.4计算数据的保存与显示

Open"C:

\DocumentsandSettings\Administrator\桌面\道路放样\计算结果\直线.txt"ForAppendAs#1        '打开文件进行写操作并保存

Print#1,"K+";li;"左边桩",Left（xil,10）,Left（yil,10）

Print#1,"K+";li;"中 桩",Left（Xi,10）,Left（Yi,10） '保存为TXT格式文档

Print#1,"K+";li;"右边桩",Left（xir,10）,Left（yir,10）

Print#1,

Close#1

Open"C:

\DocumentsandSettings\Administrator\桌面\道路放样\计算结果\直线.dat"ForAppendAs#1

Print#1,,"K+"&li;",,";Left（xil,10）;",";Left（yil,10）;",";0

Print#1,,"K+"&li;",,";Left（Yi,10）;",";Left（Xi,10）;",";0 '保存DAT格式文档

Print#1,,"K+"&li;",,";Left（xir,10）;",";Left（yir,10）;",";0

Print#1,

Close#1

“K+”表示桩号，li表示里程，xi,yi,xil,yil,xir,yir分别表示计算点的中边桩坐标。

Left（）表示取计算成果的钱10位有效数字。

在编写程序的时候，成果保存的关键词有Append,Output等，在第一次保存时用Output，目的是清除以前的旧项目数据；在同一段曲线以后计算保存中用Append，目的是在前面基础上继续添加数据，使原来的数据不被替代。

Open"C:

\DocumentsandSettings\Administrator\桌面\道路放样\计算结果\第一缓和曲线.txt"ForInputAs#1

DoWhileNotEOF

（1）                         '文本框显示计算的数据

LineInput#1,inputdata

Text5.Text=Text5.Text+inputdata+vbCrLf

Loop

Close#1

5.3.5引用CAD成图的方法

将计算结果输入到【成图数据】文件中，采用【X,Y,Z】的格式一行一行的存储在TXT文本文档中，此时Z为0。

首先引用“AutoCAD2000ObjectLibrary”，然后输入以下语句：

DimappcadAsObject      '建立Application对象

DimacadDocAsObject     '建立Document对象

DimmospaceAsObject     '建立ModelSpace对象

Setappcad=CreateObject（"autocad.application"）

IfErrThen

  Err.Clear                                   '调用CAD

  OnErrorResumeNext                         '容错

  Setappcad=CreateObject（"autocad.application"）

  IfErrThen

     Err.Clear

     MsgBox"不能运行AutoCAD,请检查是否安装！

",vbOKCancel,"警告！

"

     ExitSub

     EndIf

EndIf

appcad.Visible=True                  'CAD显示

SetacadDoc=appcad.ActiveDocument

Setmospace=acadDoc.ModelSpace

Dimstp（0To2）AsDouble

Dimetp（0To2）AsDouble

DimfnameAsString

fname="C:

\DocumentsandSettings\Administrator\桌面\道路放样\成图数据\直线中.txt"                                 'CAD成图

OpenfnameForInputAs#1            '打开数据文件

Input#1,stp（0）,stp

（1）,stp

（2）

DoWhileNotEOF

（1）

Input#1,etp（0）,etp

（1）,etp

（2）

Callmospace.AddLine（stp,etp）

stp（0）=etp（0）

stp

（1）=etp

（1）

stp

（2）=etp

（2）

Loop

Close#1

5.3.6计算的数据输出到全站仪的格式问题

软件结合公路放样的实例编制，根据实际工作需要，采用手工输入。

程序输出采用TXT文挡和DAT文档。

TXT文档便于进行打印查询及检查。

DAT文档易于将数据传输到全站仪直接进行方放样工作。

除此以外,DAT格式的文档还可直接导入到AUTOCAD中进行道路的成图,完成各种动态报表的绘制与输出。

经过学习经验和对仪器的考察，DAT文件的格式在本程序采用“点名，编码，Y,X,Z”的格式输出。

>反算－>正算后，Y坐标与原来的差了0.5－0.7mm，不知道怎么回事，这两年工作忙也没有时间再深究，但是这样的计算精度做控制足够了，如果楼主或是者是哪位同仁见此贴能顺便把这个问题解决了，咱们就一起进步了！

代码如下：

'高斯坐标正算

PrivateSubDadiZs（）

DimtAsDouble,ItpAsDouble,X0AsDouble,NAsDouble,L0AsDouble

DimVAsDouble,llAsDouble,WAsDouble,MAsDouble

Lat=Radian（Lat）

Lon=Radian（Lon）

L0=Radian（Lo）

IfTq=0Then

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）

EndIf

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#

EndSub

'高斯反算

PrivateSubDadiFs（）

DimtAsDouble,ItpAsDouble,X0AsDouble,BfAsDouble,NAsDouble

DimvAsDouble,llAsDouble,WAsDouble,MAsDouble,L0AsDouble

L0=Radian（Lo）

X0=x*0.000001

y=y-500000#

IfTq=0Then

a=6378245'54椭球参数

b=6356863.01877305

ep=0.006693421622966

ep1=0.006738525414683

f=（a-b）/a

c=a^2/b

d=b^2/a

IfX0<3Then

Bf=9.04353301294*X0-0.00000049604*X0^2-0.00075310733*X0^3-0.00000084307*X0^4-0.00000426055*X0^5-0.00000010148*X0^6

ElseIfX0<6Then

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

EndIf

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

IfX0<3Then

Bf=9.04369066313*X0-0.00000049618*X0^2-0.00075325505*X0^3-0.0000008433*X0^4-0.00000426157*X0^5-0.0000001015*X0^6

ElseIfX0<6Then

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

EndIf

EndIf

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）

EndSub

有了正反算，换带也就完成了！

用到的子程序：

PublicConstPi=3.14159265358979,p=206264.806

PublicCktqAsString

'角度化弧度

PublicFunctionRadian（aAsDouble）AsDouble

DimRoAsDouble

DimcAsDouble

DimFsAsDouble

DimIbAsInteger

DimIcAsInteger

Ifa<0Thena=-a:

t=1

Ro=Pi/180#

Ib