测量方位角计算公式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
升级会员 