% Munson/Young/Okiishi's Fluid Mechanics, 8th edition
% Example 4.3
% Code to plot streakline through (0,0) at time t

% RELEASE NOTES: v2.0, 3/14/18
%  - corrected spelling errors
%  - added ability to calculate streakline for any initial point
%  - fixed time offset (see DATA CALCLATION)

% by Wayne Clark
% ENME 320, Spring 2018 (Dr. Meilin Yu)

clear
clc

% INPUT MODE

% set Mode to 2 if you want to input your own values for the constants
%  (but note that if omega is an integer multiple of pi then
%  the x-amplitude will oscillate around zero due to the sin(omega) in 
%  the formula for x(k) - see 'DATA CALCULATION') 
Mode = 1;

% USER INPUT (if selected)
if Mode == 2
    omega = input('Frequency (multiples of pi): ');
    omega = omega*pi;
    u0 = input('Initial velocity (x): ');
    v0 = input('Initial velocity (y): ');
    ti = input('Initial time: ');
    tf = input('Final time: ');
    x0 = input('Initial x-position: ');
    y0 = input('Initial y-position: ');
    n = tf;
else
    omega = pi/8;       
    u0 = 3;
    v0 = 7;
    ti = 1;         % initial time
    tf = 50;        % final time
    n = tf;         % one point per pathline, one pathline per time step
    x0 = 0;         % default initial point is the origin
    y0 = 0;
end

% DATA SETUP
T = linspace(ti,tf,n);      % T is the initial time for which each new 
                            %   pathline is calculated
t = n;                      % calculate the streakline at the final time
x = zeros(1,n);           
y = zeros(1,n);

        % DATA CALCULATION
for k=1:n
    % each pathline needs to be offset in time by one step
    %  i.e., each 'particle' will start on its (new) pathline one step
    %  after the previous particle, because the particles are
    %  being released one-by-one as the pathlines oscillate
    
    x(k) = [u0*sin(-omega*T(k))*((T(k)-1)-t) + x0];        %  equations of pathlines 
    y(k) = [v0*((t)-(T(k)-1)) + y0];                       %  through (x0,y0) at time T(k)
end             

%  - changed T(k) to T(k)-1 because the first particle is not being delayed
%  - added x0, y0 functionality instead of being at (0,0)

% DISPLAY DATA
plot(x,y,'o',x,y)
xlabel('x')
ylabel('y')
title(['Streakline at time t = ',num2str(t),])
info = sprintf('v0 = %0.2f\nu0 = %0.2f\nomega = %0.2f\nN = %0.2f points\nx0 = %0.2f\ny0 = %0.2f', v0, u0, omega, n,x0,y0);
legend(info,'location','southeast')