% plotting of velocity, acceleration, pressure fields

clear all; % clear variables from workspace
close all; % close all figures
clc; % clear command window

t=0:0.5:24; % array of time variable
nt=length(t); % number of time steps
x=linspace(-5,5,21); %  linearly spaced vector in x
z=linspace(-5,5,21); %  linearly spaced vector in z
[X,Z]=meshgrid(x,z); % returns 2-D grid coordinates based on the coordinates 
                     % contained in vectors x and z. X is a matrix where
                     % each row is a copy of x, and Z is a matrix where 
                     % each column is a copy of z. The grid represented by 
                     % the coordinates X and Z has length(z) rows and 
                     % length(x) columns.

% constants given in the exercise                     
c1 = 5; % 1/s 
c2 = 0.3; % m
c3 = -4; % 1/s
c4 = 0.2; % m

r=1000; % kg/m3 % water density
P0 = 0; % Pa % pressure at (X0,Z0)
X0 = 0; % m
Z0 = 0; % m

% velocità
figure; % Create figure window
for i=1:nt
    U(:,:,i)=c1*X+c2*t(i); % m/s velocity in x
    W(:,:,i)=c3*Z+c4*t(i); % m/s velocity in z
    contour(X,Z,sqrt(U(:,:,i).^2+W(:,:,i).^2)) % draws a contour plot 
    hold on
    quiver(X,Z,U(:,:,i),W(:,:,i),2); % plots velocity vectors as arrows 
                                     % with components (U,W) at the points (X,Z).
    xlim([ -8,8])
    grid on;
    text(6.5,4,strcat(num2str(t(i)),' s'))  % write time on figure
    title('velocità')
    drawnow %  Update figure windows: Use this command if you modify 
            %  graphics objects and want to see the updates on the
            %  screen immediately.
    
    pause(0.1)
    hold off
end

% accelerazione
figure;
for i=1:nt
    Ax(:,:,i)=c2+(c1*X+c2*t(i))*c1; % m/s2
    Az(:,:,i)=c4+(c3*Z+c4*t(i))*c3; % m/s2
    contour(X,Z,sqrt(Ax(:,:,i).^2+Az(:,:,i).^2))
    hold on
    quiver(X,Z,Ax(:,:,i),Az(:,:,i),2);
    xlim([ -8,8])
    grid on;
    text(6.5,4,strcat(num2str(t(i)),' s'))
    title('accelerazione')
    drawnow
    pause(0.1) % Update figure windows and process callbacks
    hold off
end

% pressione
figure;
for i=1:nt
    P(:,:,i)=-r*( c1^2/2*X.^2+(c1*c2*t(i)+c2)*X-(c1^2/2*X0^2+(c1*c2*t(i)+c2)*X0) + ...
       c3^2/2*Z.^2+(c3*c4*t(i)+c4)*Z-(c3^2/2*Z0^2+(c3*c4*t(i)+c4)*Z0) ) + P0; %Pa
    contourf(X,Z,P(:,:,i),linspace(-1.4e6,3e4,21)) % Filled contour plot.
    colorbar
    caxis([-1.4e6,3e4])
    xlim([ -8,8])
    grid on;
    text(6.5,4,strcat(num2str(t(i)),' s'))
    title('pressione')
    drawnow % Update figure windows and process callbacks
    pause(0.1)
end