function [X,Y] = quadratic_interpolant(x,y,N)
% We're given data in the vectors x and y.  The output will be a
% linear fitting where there are N-1 new points introduced between each
% x(j) and x(j+1).  

n = length(x);
if mod(n,2) == 0
    error('must have an odd number of datapoints!');
end

j=1;
j1 = 2*(j-1)+1;
j2 = 2*(j-1)+2;
j3 = 2*(j-1)+3;
dx = (x(j2)-x(j1))/N;
xx = x(j1):dx:x(j2);
dx = (x(j3)-x(j2))/N;
xx(N+1:2*N+1) = x(j2):dx:x(j3);
I = (j-1)*(2*N)+1:j*(2*N)+1;
X(I) = xx;
p1 = (xx-x(j2)).*(xx-x(j3))/((x(j1)-x(j2))*(x(j1)-x(j3)));
p2 = (xx-x(j1)).*(xx-x(j3))/((x(j2)-x(j1))*(x(j2)-x(j3)));
p3 = (xx-x(j1)).*(xx-x(j2))/((x(j3)-x(j1))*(x(j3)-x(j2)));
Y(I) = y(j1)*p1 + y(j2)*p2 + y(j3)*p3;
for j=2:(n-1)/2
  j1 = 2*(j-1)+1;
  j2 = 2*(j-1)+2;
  j3 = 2*(j-1)+3;
  dx = (x(j2)-x(j1))/N;
  xx = x(j1):dx:x(j2);
  dx = (x(j3)-x(j2))/N;
  xx(N+1:2*N+1) = x(j2):dx:x(j3);
  I = (j-1)*(2*N)+1:j*(2*N)+1;
  X(I) = xx;
  p1 = (xx-x(j2)).*(xx-x(j3))/((x(j1)-x(j2))*(x(j1)-x(j3)));
  p2 = (xx-x(j1)).*(xx-x(j3))/((x(j2)-x(j1))*(x(j2)-x(j3)));
  p3 = (xx-x(j1)).*(xx-x(j2))/((x(j3)-x(j1))*(x(j3)-x(j2)));
  Y(I) = y(j1)*p1 + y(j2)*p2 + y(j3)*p3;
end
X(N*(n-1)+1) = x(n);
Y(N*(n-1)+1) = y(n);