%% ================================================================
% John Kerl
% kerl.john.r@gmail.com
% 2007-02-08
%
% A Matlab script showing a curve in 3-spaces, its derivative, the image of the
% curve under a diffeomorphism (a particular polynomial map), and the
% derivative of the latter using the Jacobian.
%% ================================================================

%% ----------------------------------------------------------------
tmin =  0;
tmax =  2*pi;
nt   = 100;
dt=(tmax-tmin)/nt;
t = [ tmin : dt : tmax ];
a = pi/4;

%% ----------------------------------------------------------------
clf;         % Clear the figure in case this is a redisplay
rotate3d on; % Permit me to rotate the graph with click and drag
hold on;     % All plot3 output goes on the same graph

%% ----------------------------------------------------------------
% Draw the circle:
% gamma(t) =
%  [x]   [ cos(t) ]
%  [y] = [ sin(t) ]
%  [z]   [     0  ]

gamma_tx = cos(t);
gamma_ty = sin(t);
gamma_tz = 0.0 * ones(1, nt+1);
plot3(gamma_tx,gamma_ty,gamma_tz, 'k');

%% ----------------------------------------------------------------
% Draw the image of the circle:
%   ([x])   [ x ]
% F ([y]) = [ y ]
%   ([z])   [ z + y^2 + 1]
F_gamma_tx = x;
F_gamma_ty = y;
F_gamma_tz = z + y.^2 + 1;
plot3(F_gamma_tx,F_gamma_ty,F_gamma_tz, 'k');

%% ----------------------------------------------------------------
% xxx temp
lo = a - pi/8;
hi = a + pi/8;
t = [ lo : (hi-lo)/nt : hi];

%% ----------------------------------------------------------------
% Draw the tangent line to the circle at t=a.
% Linear approximation:
%   L(t) = gamma(a) + gamma'(a) * (t-a)
gamma_ax = cos(a);
gamma_ay = sin(a);
gamma_az = 0;

gamma_prime_ax = -sin(a);
gamma_prime_ay =  cos(a);
gamma_prime_az =  0;

Lx = gamma_ax + gamma_prime_ax * (t-a);
Ly = gamma_ay + gamma_prime_ay * (t-a);
Lz = gamma_az + gamma_prime_az * (t-a);
plot3(Lx,Ly,Lz, 'b');

%% ----------------------------------------------------------------
% Draw the tangent line to the image of the circle.
% Linear approximation:
%   FL(a) = F(gamma(a)) + (F o gamma)'(a) * (t-a).
% Using the chain rule:
%   (F o gamma(t))' = F'(gamma(t)) * gamma'(t)
% So
%   FL(t) = F(gamma(a)) + F'(gamma(a)) * gamma'(a) * (t-a).

% Since
%   ([x])   [ x ]
% F ([y]) = [ y ]
%   ([z])   [ z + y^2 + 1]
% we have F' being the 3x3 Jacobian matrix:
%      [ 1        0        0 ]
% DF = [ 0        1        0 ]
%      [ 0        2y       1 ]

F_gamma_ax = gamma_ax;
F_gamma_ay = gamma_ay;
F_gamma_az = gamma_az + gamma_ay.^2 + 1;

DF_gamma_a = [ 1, 0, 0;
               0, 1, 0;
               0, gamma_ay, 1];

gamma_prime_a = [ gamma_prime_ax; gamma_prime_ay; gamma_prime_az ];

DF_gamma_a  = DF_gamma_a * gamma_prime_a;
DF_gamma_ax = DF_gamma_a(1);
DF_gamma_ay = DF_gamma_a(2);
DF_gamma_az = DF_gamma_a(3);

FLx = F_gamma_ax + DF_gamma_ax * (t-a);
FLy = F_gamma_ay + DF_gamma_ay * (t-a);
FLz = F_gamma_az + DF_gamma_az * (t-a);
plot3(FLx,FLy,FLz, 'b');

