I'm coding renderer for road network, which based on the RoadXML format.
Road curves from this format has four types:
And I have problem with the last one.
Clothoid is the same with Euler spiral and Cornu spiral. In the RoadXML clotho arc is given by three parameters:
For arc triangulation I need a function like foo(t), which returns (x, y) coords for t = 0..length. I created similar methods for circle arc without problems, but I can't do it for clotho arc.
Part of the problem is that I not totally understand how to apply start and end curvature parameters in standard clothoid formulas.
For example, sample RoadXML road. RoadXML sample http://img560.imageshack.us/img560/8172/bigroandabout.png
This is clotho curve item in the red ellipse. It's parameters:
I don't know how to implement these parameters, because clothoid curvature from 0 to -0.0165 is very straight.
I will happy, if you give me a code of this function (in C++, C#, Java, Python or pseudocode) or just a formula, which I can code.
Here is my equations:
x(t) ≈ t,
y(t) ≈ (t^3) / 6,
where length = t = s = curvature.
x(-0.0165) = -0.0165,
y(-0.0165) = -7.48688E-07.
Clotho length = 0.0165,
Source length = 45.185.
Scaled coords:
x'(l) = x / clotho_length * source_length = 45.185,
y'(l) = y / clotho_length * source_length = 5.58149E-07 ≈ 0.
x'(0) = 0,
y'(0) = 0.
Thus I get (0, 0)...(45, 0) points, which is very straight.
Where is my mistake? What am I doing wrong?
Let's see. Your data is:
start curvature = 0, straight line, R=INF
end curvature = -0.0165407, circular arc, R_c = 1/k_c = 60.4569335
length = 45.185. distance along clothoid, s_c = 45.185
according to Wikipedia article,
R s = const = R_c s_c ( s ~ k = 1/R by definition of clothoid )
d(s) = R d(theta)
d(theta) = k d(s)
d(theta) / d(s) = 1 / R = k = s / R_c s_c
theta = s^2 / 2 R_c s_c = (s/a)^2 = s / 2 R = k s / 2
where ___________________
a = sqrt(2 R_c s_c) (... = 73.915445 )
~~~~~~~~~~~~~~~~~~~
and so theta_c = k_c s_c / 2 (... = 0.37369576475 = 21.411190 degrees )
( not so flat after all !! )
(note: I call a here a reciprocal of what WP article calls a). Then,
d(x) = d(s) cos(theta)
d(y) = d(s) sin(theta)
x = INT[s=0..s] cos(theta) d(s)
= INT[s=0..s] cos((s/a)^2) a d(s/a)
= a INT[u=0..(s/a)] cos(u^2) d(u) = a C( s/a )
y = a INT[u=0..(s/a)] sin(u^2) d(u) = a S( s/a )
where C(t) and S(t) are Fresnel integrals.
So that's how you do the scaling. Not just t = s, but t = s/a = sqrt(theta). Here, for the end point, t_c = sqrt( k_c s_c / 2) = sqrt( 0.0165407 * 45.185 / 2) = 0.6113066.
Now, WolframAlpha says, {73.915445 Sqrt[pi/2] FresnelC[0.6113066/Sqrt[pi/2]], 73.915445 Sqrt[pi/2] FresnelS[0.6113066/Sqrt[pi/2]]} = {44.5581, 5.57259}. (Apparently Mathematica uses a definition scaled with the additional Sqrt[pi/2] factor.)
Testing it with your functions, x ~= t --> a*(s/a) = 45.185, y ~= t^3/3 --> a*(s/a)^3/3 = 73.915445 * 0.6113066^3 / 3 = 5.628481 (sic! /3 not /6, you have an error there).
So you see, using just the first term from the Taylor series representation of Fresnel integrals is not enough - by far. You have to use more, and stop only when desired precision is reached (i.e. when the last calculated term is less than your pre-set precision value in magnitude).
Note, that if you'll just implement general Fresnel integral functions for one-off scaled clothoid calculation, you'll lose additional precision when you'll multiply the results back by a (which is on the order of 102 ... 103 normally, for roads and railways).