%{
  Kinematics: rotation of a point-mass with constant
  angular velocity.

Let $z\in\mathbb{C}$ represent the position of a poin
 mass on a (complex) plane. If the point describes
a circular trajectory with constant angular velocity
$\omega$, the time rate of change of its position is
given by: 
\[
\dotd{z}\,=\,\textrm{i}\,\omega\,z
\]
where $z\,\equiv\,\rho\,\textrm{e}^{\textrm{i}\,\omega\,t}$
with $\rho\,\in\,\mathbb{R}^+$ is a constant.
The particle departs from an initial position
$z(0)\,=\,z_0$.

When using the ode15s solver, the complex-values scalar problem is converted
into a real-valued system of two first-order odes.
%}
function circular_path
close all
clc

z0     = 1+1i;
omega  = 2*pi;
T      = 20*2*pi/omega;
method = 'IE';
sample = 0.12;

t            = 0;
z            = z0;
h            = 0.5/omega;
samples      = complex(nan(ceil(T/sample),2));
samples(1,:) = [t,z0];
nsamples     = 1;
nstep        = 0;
while 1

    nstep = nstep + 1;

    if (t>(T-h))
        h = T-t;
    end

    t = t+h;

    if nstep>1
        znm1 = zn;
    end
    if nstep>2
        znm2 = znm1;
    end

    zn = z;

    switch method
        case{'Heun'}
            z = heun(@velocity,t,zn,h,omega);
        case{'AB3'}
            if nstep>3
                z = AB3(@velocity,t,zn,znm1,znm2,h,omega);
            else
                z = IE(@velocity,t,zn,h,omega);
            end

        case{'CN'}
            z = CN(@velocity,t,zn,h,omega);
        case{'IE'}
            z = IE(@velocity,t,zn,h,omega);
        otherwise
            break
    end

    if (samples(nsamples,1)>(t-sample))
        nsamples = nsamples+1;
        samples(nsamples,:) = [t,z];
    end

    if t>=T
        break;
    end
end

switch method
    case{'ode15s'}
        options = odeset('RelTol',1e-6);
      [t,z] = ode15s(@(t,z)rvelocity(t,z,{omega}),unique([0:sample:T,T]),[real(z0);imag(z0)],options);
      samples = [t(:),z(:,1)+1i*z(:,2)];
end

figure
plot(real(samples(:,2)),imag(samples(:,2)),'ko');
axis equal
grid on

end

function f = velocity(t,z,parameters)

omega = parameters{1};

f = 1i*omega*z;

end

function f = rvelocity(t,z,parameters)

omega = parameters{1};

f = omega*[-z(2);z(1)];

end

function z=heun(funct,t,zn,h,omega)

fn = funct(t,zn,{omega});

z = zn + h*fn;

f = funct(t,z,{omega});

z = zn + 0.5*h*(f+fn);

end

function z=AB3(funct,t,zn,znm1,znm2,h,omega)

% Very inefficient implementation! Recalculation of
% fnm2 and fnm1 is not strictly needed as they could be stored
% and reused.

fnm2 = funct(t,znm2,{omega});

fnm1 = funct(t,znm1,{omega});

fn   = funct(t,zn  ,{omega});

z = zn + (h/12)*(23*fn-16*fnm1+5*fnm2);

end

function z=CN(funct,t,zn,h,omega)

z = zn * ((1+0.5*1i*omega*h)/(1-0.5*1i*omega*h));

end

function z=IE(funct,t,zn,h,omega)

z = zn / (1-1i*omega*h);

end