% ---------------------------------------------
% Bottle silhouette (control points)
% ---------------------------------------------
u = [0 3 6 5 6   6  4  3  3 2.5 2];    % radial coordinate r(t) samples
v = [0 0 2 8 16 25 31 37 40 42 43];    % axial coordinate z(t) samples

% We set t(i) = i at the knots for simplicity (uniform parameter)
t = 1:length(u);                        % knot parameter values
tvals = t(1):0.05:t(end);               % fine parameter grid for evaluation

% Cubic-spline interpolation of (t,u) and (t,v)
uvals = spline(t, u, tvals);            % r(t) evaluated on tvals
vvals = spline(t, v, tvals);            % z(t) evaluated on tvals

% Figure cosmetics
set(gcf, 'color', 'white')
umax = 1.1*max(uvals); 
vmax = 1.1*max(vvals);

% Animate the 2D profile being drawn (in 5-point segments)
for i = 1:length(tvals)/5
  plot(uvals(1:5*i), vvals(1:5*i), 'LineWidth', 2)
  axis equal, axis([0 umax 0 vmax])
  drawnow
end

% Final static plot of the full silhouette
plot(uvals, vvals)
axis equal, axis([0 umax 0 vmax])

% ---------------------------------------------
% Surface of revolution (revolve profile around z-axis)
% ---------------------------------------------
theta = 0:pi/20:2*pi;                   % azimuth angles (row vector)
r = uvals;                              % r = u(t)
z = vvals;                              % z = v(t)

% Build surface grids by revolving (r,z)
xx = r' * cos(theta);                   % size: [num_tvals x num_theta]
yy = r' * sin(theta);                   % same size as xx
zz = repmat(z', 1, length(theta));      % replicate z along theta

% Render the surface
figure, surf(xx, yy, zz), shading interp
axis equal
colormap(winter)

% ---------------------------------------------
% Quick correctness check with a known surface
% Hemisphere of radius 1: volume = 2*pi/3
% Param: t in [-pi/2, 0], u = cos(t) >= 0, v = 1 + sin(t) increases
% ---------------------------------------------
t = -pi/2 : pi/20 : 0;
u = cos(t); 
v = 1 + sin(t);  % quarter circle shifted up -> hemisphere on revolution

% ParSpline2 returns (uvals, vvals, V) for a spline through (u,v):
[uvals, vvals, V] = ParSpline2(u, v, 500);
[ V  2*pi/3 ]                            % expect V ~ 2.094...

% ------------------------------------------------------------
% Scaling example: scale a shape to target volume (e.g., 192.228 ml)
% ------------------------------------------------------------
[uvals, vvals, V] = ParSpline2(u, v, 2000); % higher m -> finer integration
targetVol = 192.228;                        % e.g., 6.5 oz in milliliters
s_vol = abs(V) / targetVol;                 % ratio current volume / target
% NOTE: To scale the geometry to the target volume, multiply lengths by 
% s_len = (target / current)^(1/3).
s_len = (targetVol / abs(V))^(1/3);           % <-- length scale factor

% Apply uniform scaling to the control points (keep proportions)
u = s_len * u; 
v = s_len * v;

% ---------------------------------------------
% FUNCTIONS
% ---------------------------------------------

function [uvals, vvals, V] = ParSpline2(u, v, m)
% ParSpline2
% Inputs:
%   u, v : vectors of control points defining a C^2 piecewise-cubic spline
%   m    : number of evaluation points for numerical integration
% Outputs:
%   uvals, vvals : sampled spline coordinates on a chord-length parameter
%   V            : volume of revolution about the z-axis
%
% Assumptions:
%   - The profile is parameterized so that z = v(t) is (mostly) increasing.
%     (If not, the integral still works, but interpret orientation carefully.)

  % Compute chord lengths between successive points (arc-length proxy)
  dt = ((u(2:end) - u(1:end-1)).^2 + (v(2:end) - v(1:end-1)).^2).^(1/2);

  % Cumulative chord-length parameter t
  t = cumsum(dt);
  t = [0 t];                              % start at 0

  % Build cubic spline pp-forms for u(t) and v(t)
  ppu = spline(t, u);
  ppv = spline(t, v);

  % Evaluation grid in parameter t
  tvals = linspace(t(1), t(end), m);

  % Anonymous for u(t) and dv/dt (we only need dv/dt for this V-form)
  Su  = @(tau) ppval(ppu, tau);
  dSv = @(tau) spline_deriv(ppv, tau);

  % Sample the splines for plotting/output
  uvals = ppval(ppu, tvals);
  vvals = ppval(ppv, tvals);

  % Volume by the "method of cylindrical shells" in parametric form:
  % V = ∫ 2π * r * (arc-length element in z-direction)? 
  % Here, because we parameterize by chord length and use v'(t),
  % we're using: V = π ∫ [u(t)]^2 * v'(t) dt  (consistent if v increases)
  V = pi * integral(@(tau) Su(tau).^2 .* dSv(tau), t(1), t(end));
end

function dS = spline_deriv(pp, t)
% spline_deriv
%   Given a pp-form of a cubic spline, return its derivative evaluated at t.

  % Original cubic coefs (per piece): size [nPieces x 4] for cubic a3 x^3 + a2 x^2 + a1 x + a0
  A = pp.coefs;

  % Differentiate each piece: [3*a3, 2*a2, 1*a1] to get quadratic coefs
  Anew = A(:,1:3) * diag(3:-1:1);

  % Build a new pp for the derivative
  pp.coefs = Anew;
  pp.order = 3;                  % quadratic has order 3 in pp-form
  dS = ppval(pp, t);
end

function [uvals, vvals, V, A] = ParSpline3(u, v, m)
% ParSpline3
% Like ParSpline2 but also returns surface area A
% A = ∫ 2π r * sqrt((dr/dt)^2 + (dz/dt)^2) dt

  % Chord-length parameterization
  dt = ((u(2:end)-u(1:end-1)).^2 + (v(2:end)-v(1:end-1)).^2).^(1/2);
  t = cumsum(dt);
  t = [0 t];

  % Cubic splines pp-forms
  ppu = spline(t, u);
  ppv = spline(t, v);

  % Evaluation grid
  tvals = linspace(t(1), t(end), m);

  % Build anonymous functions for u, du/dt, dv/dt
  Su  = @(tau) ppval(ppu, tau);
  dSu = @(tau) spline_deriv(ppu, tau);
  dSv = @(tau) spline_deriv(ppv, tau);

  % Samples for output/plotting
  uvals = ppval(ppu, tvals);
  vvals = ppval(ppv, tvals);

  % Volume as before (consistent with v increasing)
  V = pi * integral(@(tau) Su(tau).^2 .* dSv(tau), t(1), t(end));

  % Surface area for surface of revolution
  A = 2*pi * integral(@(tau) Su(tau) .* sqrt(dSu(tau).^2 + dSv(tau).^2), t(1), t(end));
end