Search code examples
rmatlab

Arc length between points on Archimedean spiral using MATLAB or R

I'm needing to derive xy values and calculate arc length between each xy value, so a length value for every value in i as generated by the attached code below (excluding the origin). The points follow an Archimedean spiral path. I don't have MATLAB and am using R, but the closest I've found that I can interpret was a MATLAB example found here with credit to Jos. Below is a modified version of the MATLAB script to generate the xy data:

r = 938; %outer radius
a = 0;    %inner radius
b = 7; %increment per rev
n = (r - a)./(b); %number  of revolutions
th = 2*n*pi;      %angle  
i = linspace(0,n,n*1000); 
x = (a+b*i).* cos(2*pi*i);
y = (a+b*i).* sin(2*pi*i);

and the R equivalent:

r <- 938 # outer radius                                     
a <- 0 # inner radius    
b <- 7 # increment per revolution
n <- (r - a)/b # number  of revolutions 
th <- 2*n*pi # angle
i <- seq(0, n, length.out = n*1000) # number of points per revolution
x <- (a+b*i) * cos(2*pi*i) 
y <- (a+b*i) * sin(2*pi*i)

My assumption is that the easiest way to derive arc length between every point is to coerce i, x, and y into a MATLAB table (dataframe in R). The closest I've found for calculating arc length is this formula for calculating the total length. I'm unable to interpret math notation, so am not sure how to implement it or how to modify it to calculate arc length between every point. Using the example of the first spiral in the link above for calculating total length I tried:

sqrt((5 + 0.1289155 * 47.12389)^2 + (0.1289155)^2) * 47.12389

The link above says the result should be 378.8 but my attempt returns 521.9324. So in sum, how is the arc length between points derived in MATLAB or R?

Solution

  • The exact formula for the length, with your notations for a (start radius), r (end radius) and b (increment per revolutions) reduces to formula

    (note that in order to preserve the OP notation, there are two different meanings of the same r symbol, that might be frown upon by some)

    That formula can be implemented this way

    r <- 938 # outer radius                                     
a <- 0 # inner radius    
b <- 7 # increment per revolution
    
A <- 2 * pi / b
fa <- sqrt(1 + A^2 * a^2)
fr <- sqrt(1 + A^2 * r^2)
int_r <- (A*r*fr - log(-(A*r)+fr))/(2*A)
int_a <- (A*a*fa - log(-(A*a)+fa))/(2*A)
spiralLen <- int_r - int_a #exact formula

    394877.5

    you can also use numerical (approximative) integration in R stats integrate to evaluate the integral

    integrate(function(r){sqrt(4*pi^2*r^2/b^2+1)}, a, r)

    394877.3 with absolute error < 5.8

    Another method, that gives a rather rough approximation, but is a very good verification because it doesn't use any theoretical considerations, but just takes the data you generated - and sums up the length of the segments of all consecutive points in the data:

    dx <- x[2:length(x)] - x[1:length(x)-1]
dy <- y[2:length(x)] - y[1:length(x)-1]
len_approx = sum(sqrt(dx^2 + dy^2))

    394876.8

    As for plotting, in R, since you already have a set of points, it seems the very basic application of plot function does the job

    plot(x, y, type="l")