%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function do_all=Welander1982_flip_flop
%% Welander's flip-flop model, from Dijkstra's book:
%%    dT/dt=A(T_a-T)+k(\rho)(T_d-T)
%%    dS/dt=B(S_a-S)+k(\rho)(S_d-S)
%%    \rho=-\alpha T+\beta S
%%    k(\rho<epsilon)=0; k(\rho>epsilon)=k_0; 
%% For EPS231; Eli Tziperman
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

%% parameters:
% convection parameters:
p.k=1.5; % original value is 5
p.epsilon=-0.01;
% eqn of stae:
p.alpha=0.2;
p.beta=1;
p.T0=0;
p.S0=0;
% deep temp and salinity:
p.T_d=0;
p.S_d=0;
p.rho_d=-p.alpha*p.T_d+p.beta*p.S_d;
p.chi_d=p.alpha*p.T_d+p.beta*p.S_d;
% atmospheric values:
p.Sa=1;
p.Ta=1;
p.rho_a=-p.alpha*p.Ta+p.beta*p.Sa;
p.chi_a=p.alpha*p.Ta+p.beta*p.Sa;
% air-sea exchange coefficients for T,S:
factor=1;
p.A=1/factor;
p.B=0.1/factor;
% parameters for time integration:
time_max=20;
ode_options=odeset('RelTol',1.e-7,'abstol',1.e-8);

% steady solutions:
p.T_nonconv=p.Ta;
p.S_nonconv=p.Sa;
p.T_conv=(p.A*p.Ta+p.k*p.T_d)/(p.A+p.k);
p.S_conv=(p.B*p.Sa+p.k*p.S_d)/(p.B+p.k);

% integrate time dependent model:
hysteresis=0;
tspan=[0 1]*time_max; 
X0=[p.T_conv+0.01 p.S_conv+0.01];
[Time,X] = ode45(@(t,X) rhs(t,X,p),tspan,X0,ode_options);

figure(1); clf
%% plot time series of different elements in the solution:
subplot(2,2,[1 2])
plot(Time,-p.alpha*X(:,1),'r-')
hold on
plot(Time,p.beta*X(:,2),'-g')
rho_timeseries=-p.alpha*X(:,1)+p.beta*X(:,2);
plot(Time,rho_timeseries,'-b')
plot(Time,rho_timeseries*0+p.epsilon,'--k')
legend('-\alpha T','\beta S','\rho','\epsilon');
h1=xlabel('time');
h2=ylabel('T, S, \rho');
h3=title('Welandar flip-flop oscillations');
set([h1,h2,h3,gca],'fontsize',14);

%% plot phase space, showing trajectory of solution and
%% non-realizable steady states:
subplot(2,2,3)
plot(p.T_conv,p.S_conv,'g*','markersize',10);
hold on
plot(p.T_nonconv,p.S_nonconv,'b*','markersize',10);
b=round(length(X(:,1))/2); e=length(X(:,1));
T=X(b:e,1);S=X(b:e,2);
for i=1:length(T)
  rho=-p.alpha*(T(i)-p.T0)+p.beta*(S(i)-p.S0);
  if rho>p.epsilon
    plot(T(i), S(i),'g.')
  else
    plot(T(i), S(i),'b.')
  end
end
axis square
%% set axis limits to include all points to be drown plus a small margin:
Tmin=min([T;p.T_nonconv;p.T_conv]);Tmax=max([T;p.T_nonconv;p.T_conv]); dT=(Tmax-Tmin)/10;
Smin=min([S;p.S_nonconv;p.S_conv]);Smax=max([S;p.S_nonconv;p.S_conv]); dS=(Smax-Smin)/10;
axis([Tmin-dT,Tmax+dT,Smin-dS,Smax+dS]);
h1=xlabel('T');
h2=ylabel('S');
h3=title('phase space');
h4=legend('Location','West','convecting','nonconvecting');
set([h1,h2,h3,gca],'fontsize',14)
set([h4],'fontsize',10)

%% plot a zoom-in in phase space on the oscillation itself, plus
%% density contours: 
subplot(2,2,4)
T=X(:,1);S=X(:,2);
b=round(length(T)/2); e=length(T);
for i=b:e
  rho=-p.alpha*(T(i)-p.T0)+p.beta*(S(i)-p.S0);
  if rho>p.epsilon
    plot(X(i,1), X(i,2),'g.')
    hold on
  else
    plot(X(i,1), X(i,2),'b.')
    hold on
  end
end
axis square
dv=0.01;
axis([0.16 0.82 0.02 0.15]);
V=axis; 
[x,y]=meshgrid(V(1):dv:V(2),V(3):dv:V(4));
density_Data=-p.alpha*x(:,:)+p.beta*y(:,:); 
VV=[-0.09 -0.06 -0.03 p.epsilon 0 0.03 0.06 0.09];
[c,h]=contour(x,y,density_Data,VV,'k');
clabel(c,h);
h1=xlabel('T');
h2=ylabel('S');
h3=title('zoom-in & \rho contours');
set([h1,h2,h3,gca],'fontsize',14)


if 0
  %% save as pdf:
  set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [7 6]);
  set(gcf, 'PaperPosition', [0 0 7 6]); % [left, bottom, width, height];
  saveas(gcf,sprintf('Figure-Welander1982_flip_flop.pdf'));
end

calculate_period(rho_timeseries,Time,p)

if 0
  %% plot nullclines, does not seem to add much insight:
  figure(2); clf
  plot_nullclines(p,T,S)

  if 0
    %% save as pdf:
    set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [7 6]);
    set(gcf, 'PaperPosition', [0 0 7 6]); % [left, bottom, width, height];
    saveas(gcf,sprintf('Figure-Welander1982_flip_flop-nullclines.pdf'));
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mixing=k(X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=X(1);S=X(2);
rho=-p.alpha*(T-p.T0)+p.beta*(S-p.S0);
% replace the switch in k(rho) with a smooth tanh
% delta_epsilon=abs(p.epsilon)/5;
% mixing=p.k*0.5*(1+tanh((rho-p.epsilon)/delta_epsilon));
if rho>p.epsilon
  mixing=p.k;
else
  mixing=0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function myrhs=rhs(t,X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=X(1);S=X(2);
mixing=k(X,p);
myrhs=[p.A*(p.Ta-T)+mixing*(p.T0-T);
     p.B*(p.Sa-S)+mixing*(p.S0-S)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_nullclines(p,T,S)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Trying to switch from T,S to 
%% rho=-alpha T+beta S, chi=alpha T+beta S,  
%% and to plot nullclines. doesn't seem to add much insight.
%   transform from T,S to new coordinates,
%   \begin{align*}
%     \rho&=-\alpha T+\beta S\\
%     \chi&=\alpha T+\beta S
%   \end{align*}
%   then, multiply by eqn of state coefficients,
%   \begin{align*}
%     \frac{d(-\alpha T)}{dt}&=A(-\alpha T_a+\alpha T)+k(\rho)(-\alpha T_d+\alpha T)\\
%     \frac{d(\beta{}S)}{dt}&=B(\beta S_a-\beta S)+k(\rho)(\beta S_d-\beta S)
%   \end{align*}
%   rewrite using,
%   \begin{align*}
%     A&=\frac12(A+B)+\frac12(A-B)\\
%     B&=\frac12(A+B)-\frac12(A-B)
%   \end{align*}
%   to find,
%   \begin{align*}
%     \frac{d(-\alpha T)}{dt}&=\frac12(A+B)(-\alpha T_a+\alpha T)
%                              +\frac12(A-B)(-\alpha T_a+\alpha T)+k(\rho)(-\alpha T_d+\alpha T)\\
%     \frac{d(\beta{}S)}{dt}&=\frac12(A+B)(\beta S_a-\beta S)
%                             -\frac12(A-B)(\beta S_a-\beta S)+k(\rho)(\beta S_d-\beta S)
%   \end{align*}
%   add these two equations and then subtract the first from the second,
%   to find,
%   \begin{align*}
%   \frac{d\rho}{dt}&=\frac12(A+B)(\rho_a-\rho)-\frac12(A-B)(\chi_a-\chi)
%                       +k(\rho)(\rho_d-\rho)\\
%     \frac{d\chi}{dt}&=\frac12(A+B)(\chi_a-\chi)-\frac12(A-B)(\rho_a-\rho)
%                       +k(\rho)(\chi_d-\chi)
%   \end{align*}
%   $\rho$-nullcline satisfies,
%   \begin{align*}
%     \frac{d\rho}{dt}&=0=\frac12(A+B)(\rho_a-\rho)-\frac12(A-B)(\chi_a-\chi)
%                       +k(\rho)(\rho_d-\rho)\\
%     \rho&=\frac{\frac12(A+B)\rho_a-\frac12(A-B)(\chi_a-\chi)
%                       +k(\rho)\rho_d}{\frac12(A+B)+k(\rho)}
%   \end{align*}
%   $\chi$-nullcline satisfies,
%   \begin{align*}
%     \frac{d\chi}{dt}&=0=\frac12(A+B)(\chi_a-\chi)-\frac12(A-B)(\rho_a-\rho)
%                       +k(\rho)(\chi_d-\chi)\\
%     \chi&=\frac{\frac12(A+B)\chi_a-\frac12(A-B)(\rho_a-\rho)
%                       +k(\rho)\chi_d}{\frac12(A+B)+k(\rho)}
%   \end{align*}
%   Where in both cases we need to check if the $\rho$ value corresponds to a convecting or non-convecting state and set $k(\rho)$ accordingly.

% set axes for plotting:
rho=-2:0.0001:2;
chi=-2:0.0001:2;

% rho nullclines:
rho_null_conv=   (0.5*(p.A+p.B)*p.rho_a-0.5*(p.A-p.B)*(p.chi_a-chi)+p.k*p.rho_d) ...
    /(0.5*(p.A+p.B)+p.k);
rho_null_nonconv=(0.5*(p.A+p.B)*p.rho_a-0.5*(p.A-p.B)*(p.chi_a-chi)) ...
    /(0.5*(p.A+p.B));
rho_null_conv(rho_null_conv<p.epsilon)=NaN;
rho_null_nonconv(rho_null_nonconv>p.epsilon)=NaN;

% chi nullclines:
chi_null_conv=(0.5*(p.A+p.B)*p.chi_a-0.5*(p.A-p.B)*(p.rho_a-rho)+p.k*p.chi_d) ...
    /(0.5*(p.A+p.B)+p.k);
chi_null_nonconv=(0.5*(p.A+p.B)*p.chi_a-0.5*(p.A-p.B)*(p.rho_a-rho)) ...
    /(0.5*(p.A+p.B));
chi_null_conv(rho<p.epsilon)=NaN;
chi_null_nonconv(rho>p.epsilon)=NaN;


% steady solutions:
rho_conv=-p.alpha*p.T_conv+p.beta*p.S_conv;
chi_conv=p.alpha*p.T_conv+p.beta*p.S_conv;
rho_nonconv=-p.alpha*p.T_nonconv+p.beta*p.S_nonconv;
chi_nonconv=p.alpha*p.T_nonconv+p.beta*p.S_nonconv;

% plot:
set(0,'defaulttextfontsize',18); set(0,'defaultaxesfontsize',18);
% rho nullclines:
plot(chi, rho_null_conv,'-r','linewidth',2);
hold on;
plot(chi, rho_null_nonconv,'-g','linewidth',2);
% chi nullclines:
plot(chi_null_conv, rho, '--r','linewidth',2);
plot(chi_null_nonconv, rho, '--g','linewidth',2);
% plot steady states:
plot(chi_conv,rho_conv,'r*','markersize',10);
plot(chi_nonconv,rho_nonconv,'g*','markersize',10);
plot(chi,p.epsilon+chi*0,'--k','linewidth',0.5);
legend('rho nullcline conv','rho nullcline nonconv' ...
       ,'chi nullcline conv','chi nullcline nonconv' ...
       ,'conv Steady State','nonconv Steady State' ...
       ,'rho=epsilon' ...
       ,'Location','northwest');

% superimpose results of time integration:
rho_timeseries=-p.alpha*T+p.beta*S;
chi_timeseries=p.alpha*T+p.beta*S;
b=round(length(T)/2); e=length(T);
plot(chi_timeseries(b:e),rho_timeseries(b:e),'k.');

%xlim([-2 2]);
%ylim([-2 2]);
xlabel('$\chi=\alpha T+\beta S$','interpreter','LaTeX');
ylabel('$\rho=-\alpha T+\beta S$','interpreter','LaTeX');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function calculate_period(rho,Time,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b=round(length(Time)/2);
e=length(Time)-1;
n=b;
% find transition from nonconvecting to convecting:
step=1;
while step==1
  n=n+1;
  if n<e
    if rho(n)<p.epsilon && rho(n+1)>p.epsilon
      time_start_conv1=Time(n);
      step=0;
    end
  else
    disp('***trouble1'),n
    return
  end
end
% find next transition from convecting to nonconvecting:
step=1;
while step==1
  n=n+1;
  if n<e
    if rho(n)>p.epsilon && rho(n+1)<p.epsilon
      time_end_conv=Time(n);
      step=0;
    end
  else
    disp('***trouble2'),n
    return
  end
end
% find next transition from nonconvecting to convecting:
step=1;
while step==1
  n=n+1;
  if n<e
    if rho(n)<p.epsilon && rho(n+1)>p.epsilon
      time_start_conv2=Time(n);
      step=0;
    end
  else
    disp('***trouble3'),n
    return
  end
end
tau_conv=time_end_conv-time_start_conv1;
tau_nonconv=time_start_conv2-time_end_conv;
period=tau_conv+tau_nonconv;

fprintf(1,'k=%g; tau_conv=%g; tau_nonconv=%g; period=%g\n', ...
        p.k,tau_conv,tau_nonconv,period);
