%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function van_der_pol_periodic_forcing(gamma,tau,alpha,beta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Van-Der-Pol-Oscillator
% This iterates forward the VDP Oscillator given the forcing field F(t)
% The VDP oscillator corresponds to the equation
% x'' + \mu (x^2 - 1)x' + x = F(t)
% Defining y = x' + \mu(x^3/3 - x)
% y' = x'' + \mu(x^2 - 1)x' = F(t) - x
% x' = y - \mu(x^3/3 - x) = y + \mu(x - x^3/3)
% examples:
% van_der_pol_periodic_forcing(0,1,0.02,0) % self-sustained, weakly nonlinear
% van_der_pol_periodic_forcing(0,-2,0.02,0) % damped oscillations

if nargin < 4
    beta = .75; % Forcing Constant
end

if nargin < 3
    alpha = 30; % Forcing Constant
end

if nargin < 2
    tau = 3.6; % Years, RHS constant
end

if nargin < 1
    gamma = 0; % Forcing constant
end

% Using a simple forward step

dt = .01; % Timestep in years

Tend = 1000; % Final time in years
nt = Tend / dt; % Timesteps
T = linspace(0,Tend,nt+1); % Time vector from 0 to Tend


X = 0*T; % X vector
Y = 0*T; % Y vector
DX = X; 
DY = Y; 
X(1) = -.005;
Y(1) = .005; 


for i = 1:nt
    [DX(i),DY(i)] = calc_rhs_VDP(X(i),Y(i),T(i),gamma,tau,alpha,beta); 
    X(i+1) = X(i) + DX(i) * dt; 
    Y(i+1) = Y(i) + DY(i) * dt; 
end

%%
figure(1); clf;
plot(X,Y);
xr=get(gca,'xlim');
yr=get(gca,'ylim');
 
%%
figure(2); clf;
plot(X(100/dt+50:1/dt:end),Y(100/dt+50:1/dt:end),'ko')
xlim(xr)
ylim(yr)
 
%%
[pks,locs]=findpeaks(X);
locmods=floor(mod(locs*dt*12,12));
locmods=[locmods,[1:12]];

figure(3); clf;
histogram(locmods)

figure(4); clf
set(0,'defaulttextfontsize',20); set(0,'defaultaxesfontsize',20);
plot(T,X,'linewidth',2)
%% save as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [3 1.8]);
set(gcf, 'PaperPosition', [0 0 3 1.8]); % [left, bottom, width, height];
exportgraphics(gcf,'figure.pdf');


disp(sprintf('min,max period in month: %3.1f, %3.1f',min(locs(2:end)*dt*12-locs(1:end-1)*dt*12),max(locs(2:end)*dt*12-locs(1:end-1)*dt*12)))

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx,dy] = calc_rhs_VDP(x,y,t,gamma,tau,alpha,beta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


P1 = 1; % Years, time constant for forcing

% Forcing
F = gamma * sin((2*pi/P1)*t);

dx = (1/tau) * (F + beta - y);
dy = (alpha/tau) * (y - y^3/3 + x);

end
