function do_all=linear_stability_Schopf_Suarez_1989
% solve the dispersion relation from the linearized version of the
% Schopf and Suarez delayed oscillator for a given value of alpha and
% delta.  see p 22 in 
% http://www.deas.harvard.edu/climate/eli/reprints/WH_GFD_PDFvol2001/lect_02_2001.pdf.
% 2009, 2015, Eli Tziperman for EPS231
close all; clear all

p=set_parameters;

% start with a set of random initial guesses 
fsolve_options=optimset('display','off','tolx',1.e-10,'tolfun',1.e-10);

n=0;
%% Specify the range of the random initial guesses.  This needs to
%% be fine tuned after running this once and looking at the
%% solution.  Specify the range in the real and imaginary parts
%% separately to cover the appropriate range seen in preliminary
%% solutions: 
for real_sigma=-0.1:0.01:0.3
  for imag_sigma=-5:0.5:5
    n=n+1;

    %% initial conditio for real sigma:
    X0(1)=real_sigma; %4*(rand(1,1)-0.5); 
    %% for imaginary sigma:
    X0(2)=imag_sigma; %40*(rand(1,1)-0.5);

    %% solve the dispersion relation to find real and imag sigma:
    XX = fsolve(@(X)Equations(X,p),X0,fsolve_options);

    %% make sure the equations are actually satisfied at the final solution
    %% by looking at the norm of the equations for the initial guess and
    %% final solution:
    norm_before=norm(Equations(X0,p));
    norm_after=norm(Equations(XX,p));
    if norm_after>1.e-5
      Xfsolve(n,:)=NaN*XX(:);
    else
      Xfsolve(n,:)=XX(:);
    end
  
    %fprintf('n=%g; X0=%g, %g; X=%g, %g; |Eqns| before=%g; after=%g\n' ...
    %%        ,n,X0(1),X0(2),Xfsolve(n,1),Xfsolve(n,2),norm_before,norm_after)
  end
end

figure(1); clf
%subplot(2,1,1)
plot(Xfsolve(:,1),Xfsolve(:,2),'rx');
h1=title(sprintf('Schopf-Suarez model, T=e^{\\sigma{}t}, \\alpha=%g; \\delta=%g' ...
                 ,p.alpha,p.delta));
h2=xlabel('real(\sigma)');
h3=ylabel('imag(\sigma)');
set([gca,h1,h2,h3],'fontsize',14)
grid

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eqn=Equations(X,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Return the real and imaginary parts of the dispersion relation.

sigma_R=X(1);
sigma_I=X(2);

% Full equation is
% sigma=1-3*Tbar^2-alpha*exp(-sigma*delta)

% real part:
eqn(1)=sigma_R-(1-3*p.Tbar^2-p.alpha*exp(-sigma_R*p.delta)*cos(-sigma_I*p.delta));

% imaginary part:
eqn(2)=sigma_I-(            -p.alpha*exp(-sigma_R*p.delta)*sin(-sigma_I*p.delta));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set parameters appearing in the dispersion relation.

p.alpha=0.75;
p.delta=15;
p.Tbar=sqrt(1-p.alpha);
