
%% Exported from Jupyter Notebook
% Run each section by placing your cursor in it and pressing Ctrl+Enter

%% Markdown Cell:
% Broadly, the numerical investigation will be conducted via Matlab's `ode23t` solver.\footnote{Source Matlab files can be found on \href{http://ricopic.one/dynamic_systems_ii/source}{ricopic.one/dynamic_systems_ii/source} as \href{http://ricopic.one/dynamic_systems_ii/source/nonlinear_tank_example.m}{ricopic.one/dynamic_systems_ii/source/nonlinear_tank_example.m} and \href{http://ricopic.one/dynamic_systems_ii/source/nonlinear_fluid_state.m}{ricopic.one/dynamic_systems_ii/source/nonlinear_fluid_state.m}.}
% Of course, as with any numerical solution, specific values of the parameters must be selected.
% We begin with declaring the fluid to be water, endowing it with a density, and specify the gravitational acceleration $g$.
% Furthermore, an "anonymous" function `P_fun` is defined, accepting the height `h` of the fluid in the capacitor and returning the corresponding pressure.
% Other parameters specified include the fluid capacitance `C` and the overflow pipe inertance `L`.

%% Code Cell[116]:

save_figures = 1;

%% Code Cell[117]:

global C L % global to be used in state equation
C = 1e3; % ... fluid capacitance
L = 1e-3; % ... fluid inertance
rho = 997; % kg/m^3 ... density of water
g = 9.81; % m/s^2 ... gravitational constant
P_fun = @(h) rho*g*h; % pressure as a function of height

%% Markdown Cell:
% Next, we define the maximum height `h_max` of fluid in the capacitor, the transition height `h_t`, and the distance `dh` over which the resistor will transition from high to low impedance.

%% Code Cell[118]:

h_max = 1; % m ... maximum height of fluid
h_t = .88; % m ... transition height
dh = .05; % m ... height difference for transition

%% Markdown Cell:
% Corresponding pressures, which we prefer for computation, can be computed with `P_fun`.

%% Code Cell[119]:

P_t = P_fun(h_t); % N/m^2 ... transition pressure
dP = P_fun(dh); % N/m^2 ... pressure dif for transition

%% Markdown Cell:
% ## Nonlinear resistance
% 
% Now, let's define the variable resistance function `R_fun` ($R(P_C)$).
% We define the anonymous function via the two "limiting" resistances `R_0` ($R_0$) and `R_inf` ($R_\infty$).

%% Code Cell[120]:

R_inf = 1e-1; % N/m^2 / m^3/s ... resistance with full cap
R_0 = 1e2; % ... resistance with empty capacitor
R_fun = @(P) (R_0-R_inf)/2*(1-tanh(5/dP*(P-P_t)))+R_inf;

%% Markdown Cell:
% Let's take a moment to plot this function.
% See \autoref{fig:nonlinear_tank_example_01} for the results.
% This is a reasonable approximation of a valve that allows no flow until the capacitor fluid height reaches a threshold, then allows a significant amount of flow.

%% Code Cell[121]:

h_half = linspace(0,h_t-dh,50);
h_a = [... % heights to plot
  h_half(1:end-1),...
  linspace(h_half(end),h_max,50)...
]; 
P_a = P_fun(h_a); % N/m^2 ... pressures to plot

%% Code Cell[122]:

h = figure;
p = plot(h_a,R_fun(P_a),'linewidth',1.5);
xlabel('height (m)',...
    'interpreter','latex');
ylabel('nonlinear resistance $R$, N-s/m$^5$',...
    'interpreter','latex');
ylim([0,1.05*R_0])
grid on
hgsave(h,'figures/temp');

%% Code Cell[123]:

if save_figures
    fn = 'nonlinear_tank_example_01.tex';
    h = hgload('figures/temp');
    cleanfigure;
    matlab2tikz(['figures/',fn],...
        'figurehandle',h,...
        'width','.9\linewidth',...
        'height','.33\linewidth',...
        'extraAxisOptions',...
        [...
            'scaled ticks=false,',... 
            'tick label style={/pgf/number format/fixed},',...
            'axis lines=box'...
        ],...
        'showInfo', false...
    );
end

%% Markdown Cell:
% ### Numerical solution
% 
% The numerical ODE solver `ode23t` requires we define the first-order system of differential equations from \autoref{eq:nonlinear_fluid_example_state}.
% This is done by writing a function file `nonlinear_fluid_state.m` that the function return the time derivative of the state vector `x_a` ($\bm{x}$) at a given time.

%% Markdown Cell:
% We also pass it the nonlinear resistance function `R_fun` and the input function `Q_s_fun`.
% Let's model an offset sinusoidal input flowrate, defined as an anonymous function as follows.

%% Code Cell[124]:

Q_s_fun = @(t) 1e4*(1+sin(2*pi/1e2*t));

%% Markdown Cell:
% We're ready to simulate! The time array and zero initial conditions are specified, then simulation commences.

%% Code Cell[125]:

t_a = linspace(0,1.4e3,1e3);
x_0 = [0;0];
x_sol_struc = ode23t(...
  @(t,x) nonlinear_fluid_state(t,x,Q_s_fun,R_fun),...
  t_a,...
  x_0...
);
x_sol = deval(x_sol_struc,t_a);

%% Markdown Cell:
% We plot the results in \autoref{fig:nonlinear_tank_example_02}.
% So the overflow is relatively inactive while the capacitor fills, until $P_C$ achieves the pressure associated with a near-full capacitor.
% Then the flowrate suddenly increases rapidly due to the sudden drop in $R(P_C)$.
% Since the input is oscillating, the overflow pipe loses flowrate, then gains it again when the input flowrate increases enough to increase the capacitor pressure.
% 

%% Code Cell[126]:

h = figure;
subplot(2,1,1)
p1 = plot(...
  t_a/60,x_sol(1,:),...
  'b',...
  'linewidth',1.5...
); 
xlim([0,t_a(end)/60])
ylabel('$P_C$ response, N/m$^2$',...
    'interpreter','latex');
grid on
subplot(2,1,2)
p2 = plot(...
  t_a/60,Q_s_fun(t_a),...
  'color',.8*[1,1,1],...
  'linewidth',1.5...
);
hold on
p3 = plot(...
  t_a/60,x_sol(2,:),...
  'r',...
  'linewidth',1.5...
);
xlim([0,t_a(end)/60])
ylabel('flowrate, m$^3$/s',...
    'interpreter','latex');
grid on
xlabel('time (min)',...
    'interpreter','latex');
l = legend('$Q_s$','$Q_L$','location','northwest');
set(l,'interpreter','latex')
hgsave(h,'figures/temp');

%% Code Cell[130]:

if save_figures
    fn = 'nonlinear_tank_example_02.tex';
    h = hgload('figures/temp');
    cleanfigure;
    matlab2tikz(['figures/',fn],...
        'figurehandle',h,...
        'width','.88\linewidth',...
        'height','.88\linewidth',...
        'extraAxisOptions',...
        [...
            'scaled ticks=false,',... 
            'tick label style={/pgf/number format/fixed},',...
            'axis lines=box'...
        ],...
        'showInfo', false...
    );
end
