% plot_ptrack.m
%
% Plot particle track(s) from ptrack.txt files.
% User should configure the following variables:
%	- fname: path to ptrack.txt file
%
% Parsed data for N tracks is provided in the following variables
% 	- time: Cell array of N arrays. Each array contains the
%	  elapsed time in seconds since the start of the track. 
%	- x: Cell array of N arrays. Each array contains time series
%	  of particle x coordinates in pixels.
%	- y: Cell array of N array. Each array contains time series
%	  of particle x coordinates in pixels.
clear all; close all; clc;
fname = 'ptracks-2022-02-08-15-17-13_8rpm_cleaned.txt';

% Open ptrack file
fid = fopen(fname, 'r');

% # Create data tuples
time = {};
x = {};
y = {};

% Loop until EOF
while true

    % Find start of track while checking for EOF
    line = fgetl(fid);
    if line == -1
        break
    end
    while ~startsWith(line,'$NEW')
        line = fgetl(fid);
        if line == -1
            break
        end
    end

    % Consume header
    line = fgetl(fid);
    while ~startsWith(line, 'HR')
        line = fgetl(fid);
    end

    % Parse data
    track_time = [];
    track_x = [];
    track_y = [];
    line = fgetl(fid);
    while ~startsWith(line, '$$END')
        data = split(line);
        if length(data) == 6
            track_time(end+1) = str2double(data(1))*3600 + ...
                str2double(data(2))*60 + str2double(data(3)) + ...
                str2double(data(4))/1000;
            track_x(end+1) = str2double(data(5));
            track_y(end+1) = str2double(data(6));
        end
        line = fgetl(fid);
    end

    % Add data to tuples
    delta_time = track_time - track_time(1);
    time{end+1} = delta_time;
    x{end+1} = track_x;
    y{end+1} = track_y;
    
end

fclose(fid);

% Plot overhead view of tracks
figure; hold on;
labels = {};
for i = 1:length(time)
    plot(x{i}, y{i});
    labels{i} = sprintf('Track %d', i);
end
xlabel('x (pixels)')
ylabel('y (pixels)')
title('Overhead view')
axis equal
set(gca, 'YDir', 'reverse')
legend(labels)

% Plot time series of x coordinate
figure; hold on;
labels = {};
for i = 1:length(time)
    plot(time{i}, x{i});
    labels{i} = sprintf('Track %d', i);
end
xlabel('Elapsed time (s)')
ylabel('x (pixels)')
title('Time series')
legend(labels)

%% Example velocity calculation

% Get data for track
t = time{2};
xp = x{2};
yp = y{2};

% Plot track
figure; hold on;
plot(xp, yp);
xlabel('x (pixels)')
ylabel('y (pixels)')
title('Track 2')
axis equal

%% Convert positions to physical coordinates
xc = 322; yc = 207;   % center (pixels)
rp = 213;             % radius (pixels)
rm = 0.145;            % radius (meters)

xm = (xp - xc)*rm/rp;   % calculate x position (meters)
ym = -(yp - yc)*rm/rp;  % calculate y position (meters)

figure; hold on;
plot(xm, ym);
xlabel('x (m)')
ylabel('y (m)')
title('Track 2')
axis equal

%% Calculate radius and azimuthal velocity
r = sqrt(xm.^2 + ym.^2);            % calculate radius (meters)
vx = gradient(xm)./gradient(t);     % calculate velocity in x (m/s)
vy = gradient(ym)./gradient(t);     % calculate velocity in y (m/s)
vtheta = vy.*xm./r - vx.*ym./r;     % calculate azimuthal velocity (m/s)

figure; hold on;
plot(r, vtheta);
xlabel('r (m)')
ylabel('v_{\theta} (m s^{-1})')
title('Track 2')

%% Compare with angular-momentum-conserving prediction
T = 60/8;                           % calculate rotation period (s)
omega = 2*pi/T;                     % calculate rotation rate (rad/s)
vtheta_pred = omega*(rm.^2 - r.^2)./r;  % calculate predicted v_theta (m/s)

figure; hold on;
plot(r, vtheta, 'b');
plot(r, vtheta_pred, 'k--')
legend('data', 'theory')
xlabel('r (m)')
ylabel('v_{\theta} (m s^{-1})')
title('Track 2')

%% Compare Rossby number with prediction
Ro = abs(vtheta)./(2*omega*r);      % calculate Ro from data (nondim)
Ro_pred = 0.5*((rm./r).^2 - 1);     % calculate predicted Ro (nondim)

figure; hold on;
plot(r, Ro, 'b')
plot(r, Ro_pred, 'k--')
legend('data', 'theory')
xlabel('r (m)')
ylabel('Ro (nondim)')
title('Track 2')
set(gca, 'YScale', 'log')
ylim([0.1, 20])

