本人整理了线路中桩坐标和切线方位角计算的两种方法，拿出来与大家探讨，目的是对您的测量工作能有所帮助。这两种方法尚未经过实践验证，如果您发现它们有错误或缺陷，别忘了告诉我。
方法一：线路中桩坐标和切线方位角计算（CASIO-4800P程序）
适用于直线、第一缓和曲线（ZH-HY）、园曲线、卵形曲线（YH-HY，可以从大半径到小半径，也可以从小半径到大半径）、第二缓和曲线（YH-HZ）。
输入变量：A=起点桩号；B=终点桩号；D=起点切线方位角；F=起点X坐标；H=起点Y坐标；K=起点曲率；R=终点曲率（曲率为 1÷半径，直线的曲率为0）；M=判断因子（线路左转M=－1；线路右转M=1； 直线M=0）；C=桩号（起点和终点之间的任意桩号，A≦C≦B）。输出：X，Y=桩号为C点的坐标；P=桩号为C 点的切线方位角（单位：度）。
程序开始：
ABDFHKRM
Lbl 0:{C}
K-R=0=>S=C-A:Goto 5:≠=>Goto2 ⊿
Lbl 2:G=ABS((B-A)÷(R-K))
K=0=>L=0:S=C-A:Goto 1 ⊿
R=0=>L=B-A:S=L-C+A:Goto 1 ⊿
L=KG
K<R=>S=L+C-A: ≠=>S=L-C+A ⊿
Lbl1:N=S-L-(S^5-L^5)/(40G^2)+(S^9-L^9)/(3456G^4)-(S^13-L^13)/(599040G^6)+(S^17-L^17)/175472640G^8)-(S^21-L^21)/(78033715200G^10)
E=(S^3-L^3)/(6G)-(S^7-L^7)/(336G^3)+(S^11-L^11)/(42240G^5)-(S^15-L^15)/(9676800G^7)+(S^19-L^19)/(3530096640G^9)-(S^23-L^23)/(1880240947200G^11)
K>R=>N=-N:Z=-1: ≠=>N=N：Z=1 ⊿
T=D-90MZL^2÷(πG)
X=F+NcosT-MEsinT ◢
Y=H+NsinT+MEcosT ◢
P=D+90M Abs(S^2-L^2)/(πG)
P<0=>P=P+360⊿
P>360=>P=P-360⊿
P=P◢
Goto0
Lbl5:K=0=>Goto 6: ≠=>Goto 7 ⊿
Lbl6:X=F+ScosD ◢
Y=H+SsinD ◢
P=D◢
Goto0
Lbl7:X=F+2sin(90SR÷π)×cos(D+90MSR÷π) ÷R ◢
Y=H+2sin(90SR÷π)×sin(D+90MSR÷π) ÷R ◢
P=D+180SR÷π
P<0=>P=P+360⊿
P>360=>P=P-360⊿
P=P◢
Goto 0
程序结束
线路中桩坐标和切线方位角计算公式（方法一）
A=起点桩号，B=终点桩号，C=AB上任意点桩号，D=起点切线方位角，
X0=起点X坐标，Y0=起点Y坐标，M=左转为-1；右转为1；直线为0，
K=起点曲率，R=终点曲率。
直线：Ｋ＝R＝０，S=C-A
X=X0+ScosD
Y=Y0+SsinD
P=D
园曲线：Ｋ＝Ｒ≠０， S=C-A
X=X0+2sin（90SR/π）.cos（D+90MSR/π）/R
Y=Y0+2sin（90SR/π）.sin（D+90MSR/π）/R
P=D+180MSR/π
缓和曲线：曲线参数G=Abs（（B-A）/（R-K））
第一缓和曲线：Ｋ＝０，Ｒ≠０，L=0，S=C-A
第二缓和曲线：Ｋ≠０，Ｒ＝０，L=B-A，S=L-（C-A）=B-C
是卵形曲线： Ｋ≠０，Ｒ≠0，Ｋ≠Ｒ，L=KG
Ｋ＜R，S=L+(C-A), Z=1
Ｋ＞R，S=L-(C-A), Z=-1
Ｎ＝（S-L）－（S5-L5）／40G２＋（S9-L9）／3456G４-（S13-L13）／599040G6
＋（S17-L17）／175472640G8-（S21-L21）／78033715200G10
Ｅ＝（S３－L３）／６G－（S７－L７）／336G３＋（S１１－L１１）／42240G５-（S15-L15）／9676800G7＋（S19-L19）／3530096640G9-（S23-L23）／1880240947200G11
K＞R ,N=-N
T=D-90MZL2/πG
X=X0+NcosT-MEsinT
Y=Y0+NsinT+MEcosT
P=D+90M Abs（S2-L2）/πG
适用于第一缓和曲线、第二缓和曲线、卵型曲线（R1>R2或R2>R1）。该方法无须知道起点切线方位角，但要已知终点坐标。如果想用它求起点切线方位角，只要输入的桩号C比A大1毫米就可以了。
1.输入信息
A=起点桩号；B=终点桩号；XA，YA=起点坐标；XB，YB=终点坐标；R1=起点半径；R2=终点半径；K=左偏为-1，右偏为1；C=AB之间任意点桩号，A＜C≤B。直缓（缓直）点半径输1E50；Y$输“Y”，结束，回车继续下一个桩号的计算。
2.输出信息
XP，YP=桩号为C 点的坐标；AT=桩号为C 点的切线方位角（度.分秒）。
3.源程序：
10:INPUT"A=";A:INPUT"B=";B:INPUT"XA=";XA:INPUT"YA=";YA:INPUT"XB=";XB:INPUT"YB=";YB
15:INPUT "R1=";R1:INPUT"R2=";R2:INPUT "K=";K
20:G=ABS (R1*R2*(B-A)/(R1-R2)):S=SQR ((YB-YA)^2+(XB-XA)^2)
21:P0=ASN ((YB-YA)/S)
23:IF (XB-XA)<0 THEN P0=180-P0
24:IF (XB-XA)>0 AND P0<0 THENP0=P0+360
25:IF R1<R2 THEN Z=-1:M=G/R2+B-A:GOTO30
26:Z=1:M=G/R1
30:L=B-A:GOSUB 1040
31:P1=PP
35:INPUT "C=";C:L=C-A
40:GOSUB 1040
50:P2=P1-PP:P3=P0-K*P2
60:XP=XA+DP*COS P3:YP=YA+DP*SIN P3
70:P4=180*(M*L+Z*L^2/2)/G/PI:AT=P3+K*(P4-PP)
80:PRINT "XP=";XP:PRINT"YP=";YP:PRINT "AT=";DMS AT
100:INPUT "END?";Y$
110:IF Y$="Y" THEN END
120:GOTO 35
1040:X=L-L^3*(M*M/3+Z*M*L/4+L*L/20)/(2*G^2)+L^5*(M^4/5+Z*M^3*L/3+3*M^2*L^2/14+Z*M*L^3/16+L^4/144)/24/G^4
1050:Y1=L^2*(M+Z*L/3)/2/G-L^4*(M^3/2+3*Z*M^2*L/5+M*L^2/4+Z*L^3/28)/12/G^3
1060:Y2=L^6*(M^5/3+5*Z*M^4*L/7+5*M^3*L^2/8+5*Z*M^2*L^3/18+M*L^4/16+Z*L^5/176)/240/G^5
1070:Y=Y1+Y2
1080:DP=SQR (X^2+Y^2)
1090:PP=ATN (Y/X)
1095:RETURN
结束
TAG: 测量
