%% Parameter Sweep of ODEs
% This is a parameter sweep study of a 2nd order ODE system.
%
%   m*X'' + b*X' + k*X = 0
%
% We solve the ODE for a time span of 0 to 25 seconds, with initial
% conditions X(0) = 0 and X'(0) = 1. We sweep the parameters "b" and "k"
% and record the peak values of X for each condition. At the end, we plot a
% surface of the results.
%
% We speed up the time to run our simulation by the use of the parfor
% construct.
%
% Copyright 2013 The MathWorks, Inc.
% Further developed by HA 28.11.2016
%%
clear
close all

%% Parameter Sweep, parallel
% nlabs=24;  % Triton
nlabs=2; % (laptop)
parpool(nlabs)
%%
% Underdamped -> vibration
% $b^2 < 4 m k$
%
m=1;
b=0:.1:5;
k=(b.^2)/4;

%plot(b,k);grid on;xlabel('damping b');ylabel('stiffness k');
%title('k=(b^2)/4, up left: vibration')
%shg

%% Initialize parameters
m     =         1;  % mass
bVals =  .1:.1 :5;  % damping values
kVals = 0.5:.05:5;  % stiffness values
%%
% hold on
[bGrid, kGrid] = meshgrid(bVals, kVals);
BK=[bGrid(:) kGrid(:)];
% plot(BK(:,1),BK(:,2),'.')
%%
bb=BK(:,1);kk=BK(:,2);
vibr=(bb.^2 < 4*kk);
novibr=~vibr;
byes=bb(vibr);
bno=bb(novibr);
kyes=kk(vibr);
kno=kk(novibr);
%hold on
%plot(byes,kyes,'ok')
%plot(bno,kno,'y.')
%%
peakValsvibr=nan(length(byes),1);
peakValsno=nan(length(bno),1);


%% Parameter Sweep in parallel
disp('Computing Parameter Sweep...');
tic;
parfor idx = 1:numel(byes)
  % Solve ODE
  [Tyes{idx},Yyes{idx}] = ode45(@(t,y) odesystem(t, y, m, byes(idx), ...
   kyes(idx)),[0, 25], ...  % simulate for 25 seconds
   [0, 1]) ;      % initial conditions
 peakValsvibr(idx)=max(Yyes{idx}(:,1));
end
tparyes=toc  % 9.9852
             % 2.3195  numworkers=24, Triton
%%
maxindyes=find(peakValsvibr==max(peakValsvibr))
%%
tic
parfor idx = 1:numel(bno)
[Tno{idx},Yno{idx}] = ode45(@(t,y) odesystem(t, y, m, bno(idx), ...
   kno(idx)),[0, 25], ...  % simulate for 25 seconds
   [0, 1]) ;      % initial conditions

peakValsno(idx)=max(Yno{idx}(:,1));
end
tparno=toc   % 4.5696
             % 1.0864  Triton, nlabs=24
%%
maxindno=find(peakValsno==max(peakValsno))
%%

subplot(2,1,1)
plot(Tyes{1},Yyes{1})
legend('y','y''')
title(['Vibration, max amplitude, b= ' num2str(byes(1))...
    ' k= ' num2str(kyes(1)) ])
subplot(2,1,2)
plot(Tno{1},Yno{1});
legend('y','y''')
title(['No vibration, max displacement, b= ' num2str(bno(1))...
    ' k= ' num2str(kno(1)) ])
shg

%%
%{
figure
hold on
for idx=[1 10 50 500]
    plot(Tyes{idx},Yyes{idx}(:,1))
end
   %}
%%
%{
figure
hold on
for idx=[1 10 50 500]
    plot(Tno{idx},Yno{idx}(:,1))
end   
%}
%%









