% == sin test script == from J.S.
%
% generate ranges : 
%  1. range object, 2. matrix object, 3. linspace, 4. perfect decimation
%  5. scaled range object, 6. scaled matrix object 7. linear interp existing
%  8. linear interp alternative
N=1000;
range_test_benchmark=(str2num(sprintf('%.1f ',-N:0.1:N)));
range_test1=(-N:0.1:N);
range_test2=[-N:0.1:N];
range_test3=linspace(-N,N,20*N+1);
range_test4=(str2num(sprintf('%.12f ',-N:0.1:N)));
range_test5=(-N/0.1:1:N/0.1)*0.1;
range_test6=[-N/0.1:1:N/0.1]*0.1;
NN=20*N+1;
range_test7=zeros(1,NN);
for (i=1:1:NN)
  range_test7(i)=(i-1-N/0.1)*0.1;
endfor
range_test8=zeros(1,NN); A=-N; B=N;
for (i=1:1:NN)
  range_test8(i)=(A*(NN-i) + B*(i-1))/(NN-1);
endfor
% count sin(k*pi)<eps : what we want is 2*N+1 in case of the best decimation
sum(abs(sin(range_test1*pi))<eps(range_test1*pi))
sum(abs(sin(range_test2*pi))<eps(range_test2*pi))
sum(abs(sin(range_test3*pi))<eps(range_test3*pi))
sum(abs(sin(range_test4*pi))<eps(range_test4*pi))
sum(abs(sin(range_test5*pi))<eps(range_test5*pi))
sum(abs(sin(range_test6*pi))<eps(range_test6*pi))
sum(abs(sin(range_test7*pi))<eps(range_test7*pi))
sum(abs(sin(range_test8*pi))<eps(range_test8*pi))
% Measure how close to the IEEE float equivalent
norm(range_test1 - range_test_benchmark)
norm(range_test2 - range_test_benchmark)
norm(range_test3 - range_test_benchmark)
norm(range_test4 - range_test_benchmark)
norm(range_test5 - range_test_benchmark)
norm(range_test6 - range_test_benchmark)
norm(range_test7 - range_test_benchmark)
norm(range_test8 - range_test_benchmark)

% Check how well integerization of non-int limits works
sum((-N:0.1:N)*pi == (-N*pi:0.1*pi:N*pi))

% Illustrate how error is non-symmetric for existing linear interp algorithm.
figure(1); plot(abs(range_test1 - range_test_benchmark));
figure(2); plot(abs(range_test3 - range_test_benchmark));

% Check a commom problem spot with declaring a float as integer.
[1.8:0.05:1.9]
[1.8:0.05:1.9-1e-12]

% On my system indicates B is an integer
[1.8:0.05:1.9+1e-16]
% On my system indicates B is not an integer
[1.8:0.05:1.9+1e-15]

range_test_benchmark = str2num(sprintf('%.1f ',linspace(-2, .1, 22)));
range_test10 = linspace (-20, 1, 22) * .1;
%range_test10 = [-20:1:1]*.1;
norm(range_test10 - range_test_benchmark)
%range_test11 = linspace (-2, .1, 22);
range_test11 = (-2:.1:.1);
norm(range_test11 - range_test_benchmark)
A = -2;
B = .1;
N = 22;
for (i=1:1:N)
  range_test12(i)=((A*(N-i) + B*(i-1))/(N-1))*1.0;
endfor
norm(range_test12 - range_test_benchmark)
[range_test10' range_test11' range_test12'] - range_test_benchmark'*[1 1 1]

% Check how -0 behaves.
[-0:1:3]
-[-0:1:3]
[-3:1:-0]
[-0:pi/3:pi]
-[-0:pi/3:pi]
[-pi:pi/3:-0]

% Check addition
(0:3) + pi
(-1:3) - pi
pi + (2:-2:-2)
pi - (1:-0.5:-2)

