%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function van_der_pol_periodic_forcing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% 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)

% 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) = -.5;
Y(1) = .5; 

gamma =5; % Forcing Amplitude


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

%%
figure
 plot(X,Y);
 xr=get(gca,'xlim');
 yr=get(gca,'ylim');
 
 %%
 figure
 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
histogram(locmods)

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)))

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

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

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

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

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


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
