% This Octave program (should also run in MATLAB) calculates
% and plots Nielsen's spiral.
%
% Author: Steve Cheng
%
% Originally written for the PlanetMath article:
% http://planetmath.org/encyclopedia/CurvatureOfNielsensSpiral.html
%

%
% MATLAB has an exponential integral function that we could use
% directly to calculate the points on the spiral, but Octave does not, 
% unfortunately, so we calculate the points with our own numerical 
% method.
%
% The spiral will be plotted in two parts in different ways,
% to ensure accurate results at reasonable speed.
% The parameters below have been tuned to get a good picture.
%
% The first part is for a <= t <= b ("small" values of t), to be
% calculated by numerical quadrature on the cosine and sine integrals.
%
% The second part is for b <= t <= c ("large" values of t), to be
% calculated by an asymptotic expansion.
%

% Ranges for time parameter
a = 2.0;
b = 7.0;
c = 50.0;

% Step size for time parameter
dt = 0.05;

% Absolute tolerance for numerical integration;
% should be small enough so that the error is not detectable graphically.
tol = 0.001;            

% Degree for asymptotic expansion
N = 7;

%
% Compute first part of Nielsen's spiral using numerical quadrature
%
t = a:dt:b;
x = zeros(1, length(t));
y = zeros(1, length(t));

% Numerically integrate from 0 to a.
% The tests u==0 are there to avoid the division by zero,
% although the discontinuities in the integrand are actually removable.
x(1) = quad( @(u) (cos(u)-1)./(u+(u==0))   , 0, a, tol) + 0.57721566;
y(1) = quad( @(u) (u==0)+sin(u)./(u+(u==0)), 0, a, tol) - pi/2;

% Numerically integrate from t(i-1) to t(i) successively.
for i = 2:length(t)
  x(i) = quad( @(u) (cos(u)-1)./u, t(i-1), t(i), tol) + x(i-1);
  y(i) = quad( @(u)     sin(u)./u, t(i-1), t(i), tol) + y(i-1);
end

% Add back the "singular" part of the cosine integral
x = x + log(t);

%
% Compute second part of Nielsen's spiral using asymptotic series
%
tl = b:dt:c;

% Reciprocal of the coefficients of the Taylor polynomials 
% for cosine and sine, up to the Nth degree term
cs_coeff = factorial(0:N) .* (mod(0:N,2)==0) .* (-1).^(mod(0:N,4)==2);
ss_coeff = factorial(0:N) .* (mod(0:N,2)==1) .* (-1).^(mod(0:N,4)==3);

cr = cos(tl)./tl;
sr = sin(tl)./tl;
cs = polyval(fliplr(cs_coeff), 1./tl);
ss = polyval(fliplr(ss_coeff), 1./tl);

xl =   sr.*cs - cr.*ss;
yl = -(cr.*cs + sr.*ss);

% Plot the spiral and save it to file
plot([x, xl], [y, yl]);
%plot(x, y, xl, yl);

axis equal;
legend('off');
print('nielsen.eps', '-depsc2');