求C语言计算道路缓和曲线坐标程序代码

说明：该程序适用于计算器 CASIO fx-4800，4850P，可计算线路中心的缓和曲线、圆曲线、直线段，中、边桩

坐标及切线方位角。

1、A? 输入转角:左转为负，右转为正

2、R? 输入圆曲线半径

3、LS? 输入缓和曲线长度

4、JD(DK)? 输入交点里程桩号

5、X(JD)? 输入本交点X坐标

6、Y(JD)? 输入本交点Y坐标

7、FWJ? 输入待求点切线方位角

9、J? 输入0程序计算中桩，输入1程序计算边桩

10、Z? 输入里程桩号

1 A：R：C“LS”：D“JD(DK)”

2 P=C∧2/24/R-C∧4/2688/R∧3

3 Q=C/2-C∧3/240/R∧2

4 B=90C/兀/R

5 T=(R+P)tan(AbsA/2)+Q◢

6 W=(R+P)/cos(A/2)-R◢

7 L=((AbsA)-2B)兀R/180+2C◢

8 G“ZH”=D-T◢

9 H“HY”=G+C◢

10 I“QZ”=G+L/2◢

11 K“YH”=G+L-C◢

12 M“HZ”=G+L◢

13 N”X(JD)”：E”Y(JD)”：F”FWJ”: J

14 A＜0=>S=-1:≠=>S=1⊿ (提示：0为数字“0”)

15 U=F+A/2+90S

16 V=W+R

17 B=N+VcosU

18 O=E+VsinU (提示：O为字母“O”)

19 Lbl 1

20 {Z}

21 Z≤G=>L=T+G-Z

22 V=F+180

23 U=F

24 Goto 2⊿

25 Z≤H=>L=Z-G

26 V=L-L∧5/(90R∧2C∧2)

27 L=30L∧2S/(兀RC)

28 P=F+180

29 Q=F+L

30 U“FWJ”=F+3L◢

31 Goto 4⊿

32 Z≤K=>L=F+A/2+90S+180+180(Z-I)S/R/兀

33 U“FWJ”=L+90S◢

34 Goto 5⊿

35 Z≤M=>L=M-Z

36 V=L-L∧5/(90R∧2C∧2)

37 L=30SL∧2/(兀RC)

38 P=F+A

39 Q=F+A+180-L

40 U=F-3L+A◢

41 Goto4⊿

42 Z>M=>L=Z-M+T

43 U=F+A

44 V=U

45 Goto 2

46 Lbl 2

47 X=N+LcosV◢

48 Y=E+LsinV◢

49 Goto 6⊿

50 Lbl 3

51 {W}

52 P“XL”=X+Wcos(U-90) ◢

53 Q“YL”=Y+Wsin(U-90) ◢

54 P“XR”=X+Wcos(U+90) ◢

55 Q“YR”=Y+Wsin(U+90) ◢

56 Goto 1

57 Lbl 4

58 X=N+TcosP+VcosQ◢

59 Y=E+TsinP+VsinQ◢

60 Goto 6

61 Lbl 5

62 X=B+R cosL◢

63 Y=O+RsinL◢ (提示：O为字母“O”)

64 Goto 6

65 Lbl 6

66 J=1=>Goto 3⊿

67 Goto 1

注:

1、◢ 为输出指令，若在后面加上，即可显示前面的计算结果输出在屏幕上。

2、括号（）为说明，请不要输入。

计算参数:

A?=53°12′46.1〃

R?=4500

LS?=360

JD(DK)?=10021.359

T=2434.65260

W=534.31251

L=4539.32398

ZH=7586.70640

HY=7946.70640

QZ=9856.36839

YH=11766.03038

HZ=12126.03038

X(JD)?=3378226.731

Y(JD)?=456053.721

FWJ?=98°56′55.62〃

J?=1

Z?=10000

核对结果:

FWJ=487°23′2.26〃

X=3377706.668

Y=455858.525

W?=10（……偏距）

XL=3377714.613（……偏左坐标）

YL=455864.5966（……偏左坐标）

XR=3377698.722（……偏右坐标）

YR=455852.4535（……偏右坐标）

你好，计算圆曲线上的逐点坐标，计算起点为ZY（直圆）点，并非需要已知交点坐标，不过一些程序为了表达输入形式美观便捷，所以只允许输入交点坐标。

已知ZY（直圆）点坐标推算交点坐标公式：

JDx=ZYx+T*cos（F）

JDy=ZYy+T*sin（F）

式中：

JDx为交点坐标X、JDy同理。

ZYx为直圆点坐标X，ZYy同理。

T为切线长。

F为直线方位角（ZY到JD直线方位角，或前交点为本交点直线方位角）。

关于程序可以使用“[Excel]曲线坐标计算程序VBA 4.9”，这个程序功能强大，可选择交点法、线元法等计算曲线坐标。

提供几种不同的做法，供参考。

方法1：

=====

直接从绘图数据插值（经检验z数据是单调增加的），代码如下：

syms x y z

eq1=-2.*pi.*0.05415.*0.0000002.*sin(x).*sin((5.*pi./6)+x)-4./3.*pi.* ...

    0.0000002.^3.*2000.*z.*9.8+2./3.*pi.*z.*9.8.*0.0000002.^3.* ...

    (1820-1000).*(cos(x)).^3-2./3.*pi.*z.*9.8.*0.0000002.^3.* ...

    (1820+1000)-pi.*z.*9.8.*0.0000002.^2.*(y+0.0000002.*cos(x)).* ...

    (1820-1000).*(sin(x)).^2 

eq2=-sin((5.*pi./6)+x).*besselk(0,(z.*9.8.*(1820-1000)./0.05415).^ ...

    0.5.*0.0000002.*sin(x))+(z.*9.8.*(1820-1000)./0.05415).^0.5.* ...

    y.*besselk(1,(z.*9.8.*(1820-1000)./0.05415).^0.5.*0.0000002.*sin(x))

Z=solve(eq1,z)

eq=subs(eq2,z,Z)

h=ezplot(eq,[0.52 0.5245],[-1e-9 0])

X=get(h,'XData')

Y=get(h,'YData')

view(-10,35)

grid on

zz=subs(Z,{x y},{X Y})

X(zz>20000)=[]

Y(zz>20000)=[]

zz(zz>20000)=[]

set(h, 'XData', X, 'YData', Y, 'ZData', zz, 'linewidth', 2)

title('')

zlabel('z')

axis tight

box off

zi = [1 10 100 1000 5000 10000 19000 20000]

xi = interp1(zz, X, zi, 's')

yi = interp1(zz, Y, zi, 's')

vpa([xi yi].',10)

举例：

------

上面的代码中z分别取[1 10 100 1000 5000 10000 19000 20000]，得到的结果如下（每行代表一组对应的x、y）：

[      .5235988072,  .1770687259e-12]

[      .5235990919, -.3682753278e-13]

[      .5236019385, -.2134703182e-11]

[      .5236304028, -.2157265249e-10]

[      .5237568792, -.9522223488e-10]

[      .5239149030,  -.1794854554e-9]

[      .5241991463,  -.3217413726e-9]

[      .5242307131,  -.3370566367e-9]

把数据代回方程，看一下精度：

>> eq = subs([eq1 eq2], {x y z}, {xi yi zi})

eq =

  1.0e-005 *

    0.0000   -0.0000   -0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0000

    0.2096    0.2519    0.3885    0.0182    0.0055    0.0072    0.0034   -0.0010

>> norm(eq)

ans =

  5.0868e-006

说明：

------

（1）你所贴代码的前半部分本来用于说明两个曲面相交的，如果只画线，该部分属于多余；

（2）顺便修正一点小错误——原代码的这两句有点小问题：

x(z>20000)=NaN

y(z>20000)=NaN

我上次还有点奇怪，为什么坐标范围看上去不太对劲，但没有深究。刚才发现是这两句写错了，x、y应该是大写的。另外，现在为符合插值的需要应把右侧NaN改为空矩阵（[]）。

（3）有多种插值方法可用，我使用了最有利于平滑曲线的样条插值，楼主也可以试试其它做法。

（4）注意z0取值的范围最好不要超过[min(z) max(z)]=[8.45 19943]，如果超出该范围，部分方法需要通过参数指定允许interp1外插（默认情况下，样条插值允许外插，但线性插值不允许）。

（5）由于方法本身是基于部分离散点信息得到未知点的，所以很难说是否能保证精度。

（6）兼容性：这里的绘图代码在2007b上测试没问题，但在6.5上不行，更高版本没试，也可能会有问题。

方法2：

=====

从原始问题出发，直接解方程，代码如下：

syms x y z

eq1=-2.*pi.*0.05415.*0.0000002.*sin(x).*sin((5.*pi./6)+x)-4./3.*pi.* ...

    0.0000002.^3.*2000.*z.*9.8+2./3.*pi.*z.*9.8.*0.0000002.^3.* ...

    (1820-1000).*(cos(x)).^3-2./3.*pi.*z.*9.8.*0.0000002.^3.* ...

    (1820+1000)-pi.*z.*9.8.*0.0000002.^2.*(y+0.0000002.*cos(x)).* ...

    (1820-1000).*(sin(x)).^2 

eq2=-sin((5.*pi./6)+x).*besselk(0,(z.*9.8.*(1820-1000)./0.05415).^ ...

    0.5.*0.0000002.*sin(x))+(z.*9.8.*(1820-1000)./0.05415).^0.5.* ...

    y.*besselk(1,(z.*9.8.*(1820-1000)./0.05415).^0.5.*0.0000002.*sin(x))

zi = [1 10 100 1000 5000 10000 19000 20000]

xi = sym( zi*0 )

yi = xi

for i = 1 : length(zi)

    z0 = zi(i)

    [xi(i), yi(i)] = solve(subs(eq1,z,z0), subs(eq2,z,z0))

end

vpa([xi-2*pi yi].',10)

举例：

------

与方法1 的测试用例相同，得到的结果如下：

[       .523598806, -.3251499235e-13]

[       .523599091, -.2887359289e-12]

[       .523601938, -.2523219710e-11]

[       .523630402, -.2159081181e-10]

[       .523756878, -.9522774887e-10]

[       .523914902,  -.1794926954e-9]

[       .524199145,  -.3217447875e-9]

[       .524230712,  -.3370556700e-9]

同样，检验一下精度：

>> eq = subs([eq1 eq2], {x y z}, {xi-2*pi yi zi})

>> eq=double(eq)

eq =

  1.0e-018 *

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

    0.0000    0.0002    0.0021    0.0182    0.0812    0.1539    0.2774    0.2907

>> norm(eq)

ans =

  4.3823e-019

说明：

------

（1）由于求解方程组得到x总是位于要求的坐标范围之外，所以把它减小2*pi。

（2）这种方法需要符号数学工具箱的支持，精度远高于前一种方法。

（3）兼容性：用solve解方程的结果可能和符号数学工具箱版本相关，我用Maple内核的6.5和2007b做了测试，可以求出想要的结果，但在MuPad内核的版本上可能会有问题。

方法3：

=====

使用数值方法解方程。本来以为第2种方法有些条件下会失效，所以考虑这种做法，但从实际情况看，方法2的效果很好，所以，这个就不做了。

====================================

最后，八卦几句：

已经多次回答楼主的问题了，有点好奇您是从哪里找到这么多复杂的表达式要求解的，和您的工作有关吗？而且难得的是，问题看上去虽然比较复杂，但根据楼主提供的信息，刚好都是可以求解的，简直就像是考试的试题一般。

顺便说一下，我注意到，这次的方程与上次我回答的问题相比，有7处由9.78换成了9.8（想来应该是重力加速度），应该是有意为之的吧？

---