% Create a GIF of a spinning traditional sinuous antenna.
% (C) Dylan Crocker 2019/02/25

clear variables
close all
clc

%% Sinuous Antenna Design Parameters [User Input]

p = 8;
alpha = 60;
tau = .75;
delta = 22.5;
r1 = 0.05;

%% GIF Settings [User Input]

nSteps = 90; % Step size = 360/nSteps
spinRate = 30; % steps/second
gifName = 'hypnotic_sinuous.gif';

%% Build the GIF

% Build the curve (one arm)
[x0, y0] = sinuousArm(p, alpha, tau, delta, r1, 'feed', 'tri');

% Make three more arms (rotate copies of the first arm)
x1 = x0*cos(pi/2) - y0*sin(pi/2);
y1 = x0*sin(pi/2) + y0*cos(pi/2);
x2 = x1*cos(pi/2) - y1*sin(pi/2);
y2 = x1*sin(pi/2) + y1*cos(pi/2);
x3 = x2*cos(pi/2) - y2*sin(pi/2);
y3 = x2*sin(pi/2) + y2*cos(pi/2);

% Draw the sinuous in a loop
frames = cell(1, nSteps);
angles = (0:nSteps - 1)*2*pi/nSteps;
fig = figure();
for idx = 1:nSteps
    
    % Rotate the arms
    x0n = x0*cos(angles(idx)) - y0*sin(angles(idx));
    y0n = x0*sin(angles(idx)) + y0*cos(angles(idx));
    x1n = x1*cos(angles(idx)) - y1*sin(angles(idx));
    y1n = x1*sin(angles(idx)) + y1*cos(angles(idx));
    x2n = x2*cos(angles(idx)) - y2*sin(angles(idx));
    y2n = x2*sin(angles(idx)) + y2*cos(angles(idx));
    x3n = x3*cos(angles(idx)) - y3*sin(angles(idx));
    y3n = x3*sin(angles(idx)) + y3*cos(angles(idx));

    % Draw the sinuous
    clf(fig) % clear the figure before drawing the rotated version
    fill(x0n, y0n, 'k');
    hold all
    fill(x1n, y1n, 'k');
    fill(x2n, y2n, 'k');
    fill(x3n, y3n, 'k');
    xlim([-r1 r1])
    ylim([-r1 r1])
    axis equal
    axis off
    ax = gca;               % get the current axis
    ax.Clipping = 'off';    % turn clipping off
    zoom(1.15)
    
    % Save the frame
    drawnow
    frame = getframe(fig);
    frames{idx} = frame2im(frame);

end
close;

% Make the GIF
for idx = 1:nSteps
    [A, map] = rgb2ind(frames{idx}, 256);
    if idx == 1
        imwrite(A, map, gifName, 'gif', 'LoopCount', Inf, 'DelayTime', 1/spinRate);
    else
        imwrite(A, map, gifName, 'gif', 'WriteMode', 'append', 'DelayTime', 1/spinRate);
    end
end
% For more info see:
% https://www.mathworks.com/help/matlab/ref/imwrite.html#btv452g-1

%% Helper Function

function [ x, y ] = sinuousArm( p, alpha, tau, delta, r1, varargin )
%SINUOUSARM Draws a sinuous arm based on the specified parameters.
%   The return values are a set of x,y coordinates for drawing the arm.

% Verify and parse the inputs.
args = inputParser;
feeds = { 'none', 'tri', 'exp' };
addRequired(args, 'p', @isscalar);
addRequired(args, 'alpha', @(x) isscalar(x) || isvector(x));
addRequired(args, 'tau', @(x) isscalar(x) || isvector(x));
addRequired(args, 'delta', @isscalar);
addRequired(args, 'r1', @isscalar);
addParameter(args, 'feed', 'none', @(x) any(validatestring(x, feeds)));
parse(args, p, alpha, tau, delta, r1, varargin{:});

p = args.Results.p;
alpha = deg2rad(args.Results.alpha);
tau = args.Results.tau;
delta = deg2rad(args.Results.delta); 
r1 = args.Results.r1;
feed = args.Results.feed;

% Handle the vector inputs: if they are scalars create vectors for the
% correct indexing in the loop. If they are already vector make sure they
% are the correct size.
if isscalar(alpha)
    alpha = ones(p, 1)*alpha;
elseif length(alpha) ~= p
    error('Alpha must be either a sclar or vector of length p.');
end
if isscalar(tau)
    tau = ones(p, 1)*tau;
elseif length(tau) ~= p
    error('Tau must be either a sclar or vector of length p.');
end

% Define a function for converting polar coordinates to cartesian. 
pol2cart = @(r,p) deal(r.*cos(p), r.*sin(p));

% This may need to be an optional input...
pts = 31; % points per cell

% Build the curve
rEnd = r1;
rStart = r1*tau(end);
x = zeros(pts + (p - 1)*(pts - 1), 2); % allocate storage
y = zeros(pts + (p - 1)*(pts - 1), 2);
for n = p:-1:1
    
    % Get the alpha and tau values for the cell
    alpha_p = alpha(n);
    tau_p = tau(n);
    
    % Calculate the sinuous curve in polar coordinates
    r = linspace(rStart, rEnd, pts);
    phi = ((-1)^n)*alpha_p*sin(pi*log(r/rEnd)/log(tau_p));
    
    if n == 1
        
        % Calcualte the upper curve
        [xt, yt] = pol2cart(r, phi + delta);
        x(1:pts,1) = xt;
        y(1:pts,1) = yt;

        % Calculate the lower curve
        [xt, yt] = pol2cart(r, phi - delta);
        x(1:pts,2) = xt;
        y(1:pts,2) = yt;
        
    else
        
        % For all but the first arm do not keep the first point to avoid
        % overlapping nodes.
        idx = (pts + (pts - 1)*(n - 2) + 1):(pts + (pts - 1)*(n - 1));
        
        % Calcualte the upper curve
        [xt, yt] = pol2cart(r(2:end), phi(2:end) + delta);
        x(idx,1) = xt;
        y(idx,1) = yt;

        % Calculate the lower curve
        [xt, yt] = pol2cart(r(2:end), phi(2:end) - delta);
        x(idx,2) = xt;
        y(idx,2) = yt;
        
    end
    
    % Reset the radius limits for the next cell
    rEnd = rStart;
    rStart = rStart*tau_p;
    
end

% Circular curve to close the ends of the arm
xr = r1*cos(linspace(delta, -delta, 11));
yr = r1*sin(linspace(delta, -delta, 11));

% Add the feed taper
feedGap = 0.00;
if strcmpi(feed, 'none')
    xf = [];
    yf = [];
elseif strcmpi(feed, 'tri')
    xf = [feedGap/2 x(1,1)];
    yf = [0.0       y(1,1)];
elseif strcmpi(feed, 'exp')
%     xf = linspace(0, x(1,1), 21);
%     yf = (exp(log(2)*xf) - 1)*y(1,1);
end

% Create column vector outputs
x = [x(:,1); xr(2:end-1).'; flipud(x(:,2)); xf(:)];
y = [y(:,1); yr(2:end-1).'; flipud(y(:,2)); yf(:)];

end
