Equal arc discretization method code

%% Archimedean spiral points with Equal angle discretization method
clear all
clc
% Parameters
R2 = 30; %outer radius (between this radius and R1 is constant arc increment)
b = 2; %incerement per rev, equivalent to feed
a = 0;    %inner radius
n = round((R2 - a)./(b));  %number  of revolutions and number of 
th = 2*n*pi;      %angle obtained for n number of revolution, for one revoultion 2*pi
%% parameters defined for constant angle like incremental x value and theta value
dtheta= 1;   %% incremental angle (degree)
eqangle= dtheta* pi/(180);   %% equal angle obtianed for each spiral revolution 
npt= 2*pi/(eqangle); %% number of points in each spiral
tpoints= (npt*n)+1;   %% total number of points in n number of revolutions
theta = linspace(0,th,tpoints); 
r= a + b.*theta./(2*pi);
% Convert polar coordinates to Cartesian coordinates
xr = r.* cos(theta);
yr = r.* sin(theta);
plot(xr, yr, '-o'); 

NOTE:

The above code works well for equal angle discretization method. The code consists the information about equal angle i.e, 1 degree, the total number of points, number of points in each spiral, total angle obtained, etc.

Simialry the code for equal arc length discretization method should cosnsists the following information such as for incremental arc length i.e., 1 mm, what is number of points around the edge, spiral length, total number of points etc

Please help to developed the code based on the above information.

Looking forward for the improvement in above code

[John D'Errico]

Just download my interparc utility from the file exchange. It will do the work for you. Of course, if this is just your homework assignment, then that would not be an option.

[KUNDAN PRASAD]

Thank you for the valueable information. However, i was looking for the modification in the same above code based on equal arc discritezation method.

[David Goodmanson]

Hi Kundan,

one significant problem is that the arc length is a complicated function of theta. For the simplest case r = b*theta, the arc length between theta1 and theta2 is

arclength = (b/2)*(theta*sqrt(1+theta^2) + log(theta + sqrt(1+theta^2)))

evaluated between theta1 and theta2. It would be convenient to invert this to obtain theta and r from arclength, but that seems highly unlkely. So the approach will be numerical.

[KUNDAN PRASAD]

Thank you. It seems that the arc length is a complication function of theta. So for the given incremental arc length (length between two equidistant points) is 1 mm, how would we develped the code with numerical approach. What modification should be done in the above code. Looking forward for the changes in the code. Once again Thank you

Answer 1

[Torsten]

warning("off")
syms t t1 t2
% Parameters
R2 = 12; %outer radius (between this radius and R1 is constant arc increment)
b = 2; %incerement per rev, equivalent to feed
a = 0;    %inner radius
n = round((R2 - a)./(b));  %number  of revolutions and number of 
th = 2*n*pi;      %angle obtained for n number of revolution, for one revoultion 2*pi
x = (a + b*t/(2*pi))*cos(t);
y = (a + b*t/(2*pi))*sin(t);
dx = diff(x,t);
dy = diff(y,t);
f = simplify(sqrt(dx^2+dy^2));
g = int(f,t,t1,t2)

g = 

l = 1;   % Assuming coordinates are in mm
T(1) = 0;
Tend = T(1);
i = 1;
while Tend < th
  Tend = double(solve(subs(g,t1,T(i))-l==0,t2));
  i = i+1;
  T(i) = Tend;
end
X = double(subs(x,T));
Y = double(subs(y,T));
plot(X,Y,'o')