function PaleoBudyko( )
%Interactively explore PaleoBudkyo.
%
CopyrightText = ...
    '\copyright 2012 Richard McGehee, University of Minnesota';
VersionText = ...
    'Version 2.0, 21 Dec 2012';
%
% Notes about parameters:
%   Parameter   number 
%       A         1       OLR  A + BT
%       Tc        2       critical ice line temperature
%       B         3       OLR  A + BT
%       C         4       heat transport 
%       a1        5       albedo  alpha_1
%       a2        6       albedo  alpha_2
%
%   Parameter variables:
%     ParLabel{k}    Label displayed over field control.
%     ParDefault(k): Default values.
%     ParMin(k):     Minimum allowed value for user input.
%     ParMax(k):     Maximum allowed value for user input.
%     ParPlotted(k): The values of the parameters used in the plots
%                      displayed.
%     Par(k):        The values of the parameters displayed in the control
%                      fields.
%
% Notes about the data:
%   data.mat is a file in this directory containing the relevant data
%       in the form of a 15001x4 matrix:
%     data(:,1): time in kyr from -10000 to +5000 in 1 kyr increments
%     data(:,2): eccentricities
%     data(:,3): obliquities (radians)
%     data(:,4): d18O data from Lisiecki-Raymo stack, interpolated using
%                  cubic splines from -5320 to 0 to 1 kyr increments and
%                  extended to the whole interval.
%   Note for future development:  the allowed times are determined by
%     data(1,1) and data(n,1), where n is the number of rows.

% --- Color scheme
BaseColor = [0.97 0.97 0.86];
FieldColor = [0.94 0.94 0.76];
ButtonColor = [0.88 0.87 0.71];
RefreshActiveColor = [0.22 0.71 0.22];
ExitColor = [0.8 0 0];

FigTitle = 'Budyko Responds to Milankovitch';
FigTitlePosition = [0.45 0.96 0.55 0.04];
CopyrightPosition = [0.94,0.01,0.425,0.04];

MilankovitchCredits = { ...
    'Computation of the Milankovitch Cycles are due to Laskar et al,', ...
    ['"A long-term numerical solution for the insolation quantities ', ...
     'of the Earth,"'], ...
    '{\itAstronomy & Astrophysics} {\bf428} (2004), 261–285.', ...
    'The computations are posted on the Web at', ...
    'http://www.imcce.fr/Equipes/ASD/insola/earth/La2004/index.html.' };
BudykoCredits = { ...
    ['A description of Budyko''s model and the computation of its ', ...
      'stable equilibria'], ...
    'can be found in McGehee & Lehman, "A Paleoclimate Model of', ...
    'Ice-Albedo Feedback Forced by Variations in Earth''s Orbit,"', ...
    '{\itSIAM J. Applied Dynamical Systems} {\bf11} (2012), 684-707.' }; 
ClimateCredits = { ...
    'Climate data are due to Lisiecki and Raymo, "A Pliocene-Pleistocene ', ...
    'stack of 57 globally distributed benthic \delta^{18}O records",', ...
    '{\itPaleoceanography} {\bf20} (2005), PA1003.' ...
    'The data are posted on the Web at', ...
    'http://lorraine-lisiecki.com/LR04stack.txt.' };

MilanTitle = 'Milankovitch Cycles';
BudTitle = 'Budyko Model';
ClimateTitle = 'Climate Data';

% --- Get all the data. 
All = load('data');

% --- Set time label, defaults, and range.
t0Label = 'begin (kyr)';                         % <--- t0 (begin time).
t0Default = -5000;
t1Label = 'end (kyr)';                           % <--- t1 (end time).
t1Default = 0;
tmin = All.data(1,1);                        %   allowed times,
tmax = All.data(length(All.data),1);         %   as determined by the data
% tmin = -10000;  tmax = 5000;   

% --- Set parameter label, defaults, and ranges.
% number       1     2      2      2         5           6
ParLabel =   {'A' , 'Tc' , 'B'  , 'C'  , 'alpha_1' , 'alpha_2'};
ParDefault = [202 , -10  , 1.9  , 3.04 ,   0.32    ,    0.62  ];
ParMin     = [ 0  , -273 , 0.01 ,  0   ,    0      ,     0    ];
ParMax =     [400 ,  100 ,  10  ,  20  ,    1      ,     1    ];

% --- Set current values.
t0 = t0Default;  t1 = t1Default;
Par = ParDefault;

% --- Set fake plotted values to induce initial plots.
t0Plotted = 0;  t1Plotted = 0;
ParPlotted = [ 0 0 0 0 0 0 ];
PSPlotted = [true true true];
PSShouldShow = [ false false false ];
CreditsShouldShow = false;

% ---- Set up enclosing figure.
set(0,'Units','pixels');
scrnPixels = get(0,'ScreenSize');
PUF = scrnPixels(4);                  % --- screen height in pixels
FigPosPix = [0.1*PUF, 0.1*PUF, 1.1*PUF, 0.8*PUF];
hFig = figure('Position',FigPosPix, 'Color', BaseColor);

PUA = 0.8*PUF;                                 % --- Pixel units for axes

% --- Make the Title.
axes('Units','pixels', 'Position',PUA*FigTitlePosition, ...
    'Visible','off');
text(0,0,FigTitle, ... 
    'VerticalAlignment','Bottom', ...
    'FontSize',14, 'FontWeight','bold');

% --- Make the copyright declaration.
axes('Units','pixels', 'Position',PUA*CopyrightPosition, ...
    'Visible','off');
text(1,1,{VersionText,CopyrightText}, ... 
    'VerticalAlignment','top', ...
    'HorizontalAlignment','right', ...
    'FontSize',8 );

% ---- Set up three sets of axes and corresponding popups.
PUA = 0.8*PUF;                                 % --- Pixel units for axes
AxPosition =      [0.12, 0.69, 0.8,  0.2;  ...  % --- Axes Positions
                   0.12, 0.39, 0.8,  0.2;  ...
                   0.12, 0.09, 0.8,  0.2] ;
AxPopupPosition = [0.12, 0.9, 0.17, 0.03;  ... % --- Popup Menu Positions
                   0.12, 0.6, 0.17, 0.03;  ...
                   0.12, 0.3, 0.17, 0.03] ;
AxHandle = [0 0 0];                            % --- Allocate axes handles.
for k = 1:3
    AxHandle(k) = axes('Units','pixels', ...
                       'Position',PUA*AxPosition(k,:) );
    uicontrol('Style','popup', ...
          'Position',PUA*AxPopupPosition(k,:), ...
          'BackgroundColor',FieldColor, ...
          'String','time series|power spectrum', ...
          'UserData',k, ...
          'Callback',{@CBPop});
end

% ---- Compute the parameter field positions.
npar = 8;
ParPos = zeros(npar+1,4);
ParPos(1,:) = [1.06*PUA 0.65*PUA 0.1*PUA 0.03*PUA];
ParNameIncr = [0 0.032*PUA 0 -0.005*PUA];
ParDefaultIncr = [0 -0.03*PUA 0 0];
ParIncrEven = [0.12*PUA 0 0 0];
ParIncrOdd = [-0.12*PUA -0.10*PUA 0 0];
if npar > 1
    for k = 2:npar
        if mod(k,2) > 0.5                    % k is odd
            ParPos(k,:) = ParPos(k-1,:) + ParIncrOdd;
        else                                 % k is even
            ParPos(k,:) = ParPos(k-1,:) + ParIncrEven;
        end
    end
end

% ---- Make the 'Refresh' button.
RefreshHandle = uicontrol( ...
          'Position',[1.11*PUA 0.80*PUA 0.12*PUA 0.04*PUA], ... 
          'BackgroundColor', ButtonColor, ...
          'ForegroundColor', RefreshActiveColor, ...
          'String','Refresh', 'FontSize', 12, ...
          'Callback',{@CBRefresh} );
RefreshInactiveColor = get(RefreshHandle,'BackgroundColor');

% ---- Make the 'Exit' button.
uicontrol('Position',[1.12*PUA 0.25*PUA 0.1*PUA 0.03*PUA], ...
          'BackgroundColor',ButtonColor, ...
          'ForegroundColor', ExitColor, ...
          'String','Exit', 'FontWeight','bold',...
          'Callback',{@Close});

% --- Parameter field for  t0.
t0han = uicontrol('Style','edit', ...
          'Position',ParPos(1,:), ...
          'BackgroundColor',FieldColor, ...
          'String',num2str(t0), ...
          'UserData',0, ...         % <--- Indicates t0 field.
          'Callback',{@CBtime});
uicontrol('Style','text', ...
          'Position',ParPos(1,:) + ParNameIncr, ...
          'BackgroundColor', BaseColor, ...
          'String',t0Label);
      
% --- Default button for time interval.
uicontrol('Position',ParPos(1,:) + ParDefaultIncr + [0.06*PUA 0 0 0], ...
          'BackgroundColor', ButtonColor, ...
          'String','Default', ...
          'Callback',{@Defaultt01});

% --- Parameter field for  t1.
t1han = uicontrol('Style','edit', ...
          'Position',ParPos(2,:), ...
          'BackgroundColor',FieldColor, ...
          'String',num2str(t1), ...
          'UserData',1, ...         % <--- Indicates t1 field.
          'Callback',{@CBtime});
uicontrol('Style','text', ...
          'Position',ParPos(2,:) + ParNameIncr, ...
          'BackgroundColor', BaseColor, ...
          'String',t1Label);

% --- Fields, labels, and default buttons for parameters 1-6
%         (A, Tc, B, C, a1, a2)
ParHandle = zeros(1,6);
for k = 1:6
    ParHandle(k) = uicontrol('Style','edit', ...
          'Position',ParPos(k+2,:), ...
          'BackgroundColor',FieldColor, ...
          'String',num2str(ParDefault(k)), ...
          'UserData',k, ... 
          'Callback',{@CBPar});
    uicontrol('Style','text', ...
          'Position',ParPos(k+2,:) + ParNameIncr, ...
          'BackgroundColor', BaseColor, ...
          'String',ParLabel{k});
    uicontrol('Position',ParPos(k+2,:) + ParDefaultIncr, ...        
          'BackgroundColor', ButtonColor, ...
          'String','Default', ...
          'UserData',k, ...
          'Callback',{@CBDefault});
end

% ---- Make the 'Credits' button.
uicontrol('Position',[1.12*PUA 0.1*PUA 0.1*PUA 0.03*PUA], ...
          'BackgroundColor',ButtonColor, ...
          'String','Credits', ...
          'Callback',{@CBCredits});

% --- Give all these variables full scope.
t = 0;       % -- uniform times from t0 to t1, 1 kyr time step
ecc = 0;     % -- eccentricities
obl = 0;     % -- obliquitis
d18O = 0;    % -- values of d18O
eta = 0;     % -- ice line locations
GMT = 0;     % -- global mean temperatures
freq = 0;    % -- frequencies for power spectru
PSecc = 0;   % -- power spectrum for eccentricity
PSobl = 0;   % -- power spectrum for obliquity
PSd18O = 0;  % -- power spectrum for d18O
PSeta = 0;   % -- power spectrum for ice line location
PSGMT = 0;   % -- power spectrum for global mean temperature

% ======================== Get everything started. ==================

Refresh();
ItIsTimeToRefresh = false;

% ========== Main Loop =============%|
                                    %|                                         
Exit = false;                       %|
while ~Exit                         %|
    uiwait(hFig);                   %|
    if AnythingIsDifferent()        %|
        activateRefresh()           %|
    else                            %|
        cancelRefresh()             %|
    end                             %|
    if ItIsTimeToRefresh            %|
        Refresh()                   %|
        ItIsTimeToRefresh = false;  %|
    end                             %|
    if CreditsShouldShow
        ShowCredits();
        activateRefresh();
        CreditsShouldShow = false;
    end
end                                 %|
close(hFig)                         %|
return                              %|
                                    %| 
% ==================================%|


% ====================================================================
% ====================== End of controlling program ==================
% ====================================================================
% ========================= Start of subprogram ======================
% ====================================================================

% ----------------------------- Plot Refresher -----------------------

    function Refresh()
    
    % --- Check to see what needs refreshing.
        if t0 ~= t0Plotted || t1 ~= t1Plotted  % <-- Refresh everything.
            [ t, ecc, obl, d18O ] = restrictData( t0, t1, All.data );
            freq = (0:length(t)-1)'/length(t);
            PSecc = computePowerSpectrum(ecc);
            PSobl = computePowerSpectrum(obl);
            PSd18O = computePowerSpectrum(d18O);            
            [ eta, GMT ] = computeBud(ecc,obl,Par);
            PSeta = computePowerSpectrum(eta);
            PSGMT = computePowerSpectrum(GMT);
        elseif any(Par ~= ParPlotted)          % <-- Refresh only Budyko.
            [ eta, GMT ] = computeBud(ecc,obl,Par);
            PSeta = computePowerSpectrum(eta);
            PSGMT = computePowerSpectrum(GMT);     
        end
    % --- Redisplay all three axes.
            % -- Insure uniform XLim for all three axes.
        timeinterval4plots = plotLaskar(AxHandle(1),t,ecc,obl);        
        if PSShouldShow(1)
            plotPS(AxHandle(1),freq, ...
                  PSecc,         PSobl, ...
                {'eccentricity','obliquity'})
        end
        title(AxHandle(1),MilanTitle);
        if PSShouldShow(2)
            plotPS(AxHandle(2),freq, ...
                  PSGMT, PSeta, ...
                {'GMT','ice line'} )
        else
            plotBud(AxHandle(2),t,GMT,eta,timeinterval4plots)
        end
        title(AxHandle(2),BudTitle);
        if PSShouldShow(3)
            plotPS(AxHandle(3),freq,PSd18O)
        else
            plotLR(AxHandle(3),[t,d18O],timeinterval4plots)
        end
        title(AxHandle(3),ClimateTitle);
        
    % --- Reset all 'plotted' indicators.
        t0Plotted = t0;  t1Plotted = t1;
        ParPlotted = Par;
        PSPlotted = PSShouldShow;
        cancelRefresh();
        % -- Cancel Parameter Warnings
        set(t0han,'ForegroundColor','k');
        set(t1han,'ForegroundColor','k');
        for i = 1:6
            set(ParHandle(i),'ForegroundColor','k')
        end        
        return
    end

    function activateRefresh()
        set(RefreshHandle,'BackgroundColor',RefreshActiveColor);
        set(RefreshHandle,'ForegroundColor','w');
    end

    function cancelRefresh()
        set(RefreshHandle,'BackgroundColor',RefreshInactiveColor);
        set(RefreshHandle,'ForegroundColor',RefreshActiveColor);
    end

    function adviseRefresh()
        if AnythingIsDifferent()
            activateRefresh()
        else
            cancelRefresh()
        end
    end

    function trueorfalse = AnythingIsDifferent()
    % Check to see if plots need updating.
        if t0 ~= t0Plotted || t1 ~= t1Plotted ...
                || any(Par ~= ParPlotted) ...
                || any(PSShouldShow ~= PSPlotted)
            trueorfalse = true;
        else
            trueorfalse = false;
        end
    end

    function ShowCredits()
        
        % --- Milankovitch Credits
        axes(AxHandle(1));
        cla reset;
        set(AxHandle(1),'YTick',[],'XTick',[]);
        text(0.5,0.5,MilankovitchCredits,'HorizontalAlignment','center');
        title(MilanTitle);
        
        % --- Budyko Credits
        axes(AxHandle(2));
        cla reset;
        set(AxHandle(2),'YTick',[],'XTick',[]);
        text(0.5,0.5, BudykoCredits,'HorizontalAlignment','center');
        title(BudTitle);
        
        % --- Climate Data Credits
        axes(AxHandle(3));
        cla reset;
        set(AxHandle(3),'YTick',[],'XTick',[]);
        text(0.5,0.5,ClimateCredits,'HorizontalAlignment','center');
        title(ClimateTitle);
        
        return
    end

% ------------------------ Callback Functions ------------------------

% --- Callback for 'Refresh' button.
    function CBRefresh(~,~)
        ItIsTimeToRefresh = true;
        uiresume(gcbf)
        return
    end

% --- Callback for 'Exit' button.
    function Close(~,~)
        Exit = true;
        uiresume(gcbf)
    end

% --- Callback for time fields (t0 and t1).
    function CBtime(source,~)
        t0ort1 = get(source,'UserData');
        str = get(source,'String');
        ChoiceWasBad = false;
        [newt,ok] = str2num(str);
        if ok
            newt = round(newt);
            if t0ort1 == 0;        % <--- 't0' field was touched.
                if newt < tmin
                    newt = tmin;
                    ChoiceWasBad = true;
                elseif newt > t1
                    newt = t1;
                    ChoiceWasBad = true;
                end
                if newt ~= t0
                    t0 = newt;
                    set(source,'String',num2str(t0));
                    adviseRefresh();
                end                
            else                   % <--- 't1' field was touched.
                if newt > tmax
                    newt = tmax;
                    ChoiceWasBad = true;
                elseif newt < t0
                    newt = t0;
                    ChoiceWasBad = true;
                end
                if newt ~= t1;
                    t1 = newt;
                    set(source,'String',num2str(t1));
                    adviseRefresh();
                end
            end
        else
            ChoiceWasBad = true;
        end
        if ChoiceWasBad
            set(source,'ForegroundColor','r');
        else
            set(source,'ForegroundColor','k');
        end
        uiresume(gcbf)
    end

% --- Callback for 'Default' button for t0 and t1.
    function Defaultt01(~,~)
        t0 = t0Default;
        t1 = t1Default;
        set(t0han,'String',num2str(t0));
        set(t1han,'String',num2str(t1));
        set(t0han,'ForegroundColor','k');
        set(t1han,'ForegroundColor','k');
        adviseRefresh(); 
        uiresume(gcbf)
    end

% % --- Callback for 'Default' button for t0 and t1.
%     function Defaultt01(~,~)
%         if t0Default == t0 && t1Default == t1
%         else
%             t0 = t0Default;
%             t1 = t1Default;
%             set(t0han,'String',num2str(t0));
%             set(t1han,'String',num2str(t1));
%             set(t0han,'ForegroundColor','k');
%             set(t1han,'ForegroundColor','k');
%             adviseRefresh();            
%         end
%         uiresume(gcbf)
%     end

% --- Callback for parameters 1 through 6  (A,Tc,B,C,a1,a2)
    function CBPar(source,~)
        parno = get(source,'UserData');
        str = get(source,'String');
        ChoiceWasBad = false;
        [candidate,ok] = str2num(str);
        if ok
            if candidate <  ParMin(parno)
                newpar = ParMin(parno);
                ChoiceWasBad = true;
            elseif candidate > ParMax(parno)
                newpar = ParMax(parno);
                ChoiceWasBad = true;
            else
                newpar = candidate;
            end
            set(source,'String',num2str(newpar));
            Par(parno) = newpar;
            adviseRefresh()
        else
            ChoiceWasBad = true;
        end
        if ChoiceWasBad
            set(source,'ForegroundColor','r');
        else
            set(source,'ForegroundColor','k');
        end
        uiresume(gcbf)            
    end

% -- Callback for 'Default' buttons
    function CBDefault(source,~)
        parno = get(source,'UserData');
        newpar = ParDefault(parno);
        set(ParHandle(parno),'String',newpar);
        set(ParHandle(parno),'ForegroundColor','k');
        Par(parno) = newpar;
        adviseRefresh();
        uiresume(gcbf)
    end

% --- Callback for Milankovitch pop-up menu.
    function CBPop(source,~)
        plotno = get(source,'UserData');
        val = get(source,'Value');
        if val == 1
            PSShouldShow(plotno) = false;
        else
            PSShouldShow(plotno) = true;
        end
        adviseRefresh();
        uiresume(gcbf)
    end

% --- Callback for the 'Credits' button.
    function CBCredits(~,~)
        CreditsShouldShow = true;
        uiresume(gcbf)
    end

end

% ====================================================================
% ========================== End of subprograms ======================
% ====================================================================
% ============== Start of appended independent functions =============
% ====================================================================

% ============================ Plotting Functions ====================

% ----------------------- Milankovitch Cycles ------------------------

function tint = plotLaskar( axh, t, ecc, obl )
%
% In:
%   axh:  axes handle where the plot is to go
%   t:    nx1 vector of times [t0,t1], incremented by 1 kyr
%   ecc:  nx1 vector of eccentricities
%   obl:  nx1 vector of eccentricities
%
% Out: 
%   tint = XLim for this plot (t interval)
%   Eccentricity and obliquity are plotted in the axes.

[handles, ~, ~] = plotyy(axh, ...
    t, ecc, ...
    t, obl*180/pi );
tint = get(axh,'XLim');
xlabel(axh,'kyr');
set(axh,'XGrid','on');

% --- Work on eccentricity axis.
set(handles(1),'YLim',[0,0.1]);
EccLabelHandle = get(handles(1),'Ylabel');
set(EccLabelHandle,'String','Eccentricity');
EccLabelPos = get(EccLabelHandle,'Position');
EccLabelPos(2) = 0.035;
set(EccLabelHandle,'Position',EccLabelPos);

% --- Work on obliquity axis.
set(handles(2),'YLim',[20,24.5]);
OblLabelHandle = get(handles(2),'Ylabel');
set(OblLabelHandle,'String','Obliquity');
OblLabelPos = get(OblLabelHandle,'Position');
OblLabelPos(2) = 23.25;
set(OblLabelHandle,'Position',OblLabelPos);

return
end

% ------------------------ Budkyko Model -----------------------------

function plotBud( axh, t, GMT, eta , tint )
%
% In:
%   axh:  axes handle where the plot is to go
%   t:    nx1 vector of times [t0,t1], incremented by 1 kyr
%   GMT:  nx1 vector of global mean temperatures
%   eta:  nx1 vector of ice line locations
%   tint: t interval for plots (XLim)
%
% Out: GMT and ice line are plotted in the axes.

[handles, ~, ~] = plotyy(axh, ...
    t, GMT, ...
    t, eta );
% title(axh,'Budyko Model');
xlim(axh,tint);
xlabel(axh,'kyr');
set(axh,'XGrid','on');

% --- Work on GMT axis.
GMTLabelHandle = get(handles(1),'Ylabel');
set(GMTLabelHandle,'String','GMT');
GMTLabelPos = get(GMTLabelHandle,'Position');
YL = get(handles(1),'YLim');
YL(2) = YL(2) + 0.8*(YL(2) - YL(1));
set(handles(1),'YLim',YL);
set(GMTLabelHandle,'Position',GMTLabelPos);

% --- Work on eta axis.
EtaLabelHandle = get(handles(2),'Ylabel');
set(EtaLabelHandle,'String','Ice Line');
EtaLabelPos = get(EtaLabelHandle,'Position'); 
YL = get(handles(2),'YLim');
YL(1) = YL(1) - 0.8*(YL(2) - YL(1));
set(handles(2),'YLim',YL);
set(EtaLabelHandle,'Position',EtaLabelPos);

return
end

% ------------------------ Climate Data -----------------------------


function plotLR( axh, LR, tint)
plot(axh,LR(:,1),LR(:,2));
% title(axh,'Climate Data');
xlabel(axh,'kyr');
xlim(axh,tint);
set(axh,'XGrid','on');
set(axh,'YDir','reverse')
ylabel(axh,'\delta^{18}O');
end

% ---------------------------- Power Spectra ------------------------

function plotPS( axh, freq, PS1, PS2, ylabel )

if nargin < 3
    InsufficientData();
    return
elseif nargin < 4
    if any(isnan(PS1))
        InsufficientData();
        return
    end
    plot(axh, freq, PS1);
    xlabel(axh,'frequency (1/kyr)');
    xlim(axh,[0 0.06]);
    set(axh,'XGrid','on');
    ylim(axh, [0 1.1]);
    set(axh,'YTick',[]);
else    
    [handles, ~, ~] = plotyy(axh, freq, PS1, freq, PS2 );
    
    % --- Work on PS1 plot
    set(handles(1),'XLim', [0 0.06]);
    set(handles(1),'YLim',[0 2.2]);
    set(handles(1),'YTick',[]);
    PS1LabelHandle = get(handles(1),'Ylabel');
    set(PS1LabelHandle,'String',ylabel{1});
    PS1LabelPos = get(PS1LabelHandle,'Position');
    PS1LabelPos(2) = 0.5;
    set(PS1LabelHandle,'Position',PS1LabelPos);
    
    % --- Work on PS2 plot    
    set(handles(2),'XLim', [0 0.06]);
    set(handles(2),'YLim',[-1.1 1.1]);
    set(handles(2),'YTick',[]);
    PS2LabelHandle = get(handles(2),'Ylabel');
    set(PS2LabelHandle,'String',ylabel{2});
    PS2LabelPos = get(PS2LabelHandle,'Position');
    PS2LabelPos(2) = 0.5;
    set(PS2LabelHandle,'Position',PS2LabelPos);
    
    xlabel(axh,'frequency (1/kyr)');
    set(axh,'XGrid','on');
end
return

    function InsufficientData()
        axes(axh);
        cla reset;
        text(0.5,0.5,'insufficient data','HorizontalAlignment','center');
        return
    end
end

% ---------------------------- Credit Display -----------------------

% function DisplayCredits(axh,Credits)
% axes(axh);
% cla reset;
% set(axh,'YTick',[],'XTick',[]);
% text(0.5,0.5,Credits,'HorizontalAlignment','center')
% 
% % TextVertPos = 0.8;
% % VertIncr = -0.12;
% % for i = 1:length(Credits)
% %     text(0.5,TextVertPos,Credits{i}, ...
% %         'HorizontalAlignment','center');
% %     TextVertPos = TextVertPos + VertIncr;
% % end
% return
% end


% ====================================================================

% ========================= Data Restricter ==========================

function [ t, ecc, obl, d18O ] = restrictData( t0, t1, alldata )
%Restricts data to the time interval.
%   
% In: 
%   t0, t1 define the time interval [t0, t1]
%     Assumptions:  t0 and t1 are integers, t0 <= t1, 
%       and t0 and t1 are within the time range of alldata.
%   alldata: nx4 matrix
%     alldata(:,1): times in kyr, incremented by 1 kyr
%     alldata(:,2): eccentricity
%     alldata(:,3): obliquity
%     alldata(:,4): d18O
%
% Out:
%   t:    nx1 vector of times [t0,t1], incremented by 1 kyr
%   ecc:  nx1 vector of eccentricities
%   obl:  nx1 vector of eccentricities
%   d18O: nx1 vector of d18O values

% --- Convert the time interval to array indices.
N = size(alldata,1);
n1 = t0 - alldata(1,1) + 1;
n2 = t1 - alldata(1,1) + 1;

% --- Throw out ridiculous requests.
if n1 > n2 || n1 > N || n2 < 1
    t = []; ecc = []; obl = []; d18O = [];
    return
end

% --- Get the data for the time interval.
n1 = max(n1,1);
n2 = min(n2,N);
t = alldata(n1:n2,1);
ecc = alldata(n1:n2,2);
obl = alldata(n1:n2,3);
d18O = alldata(n1:n2,4);

return
end
% ====================================================================

% ====================== Power Spectrum Computer =====================

function PS = computePowerSpectrum(data)
%Computes the relative power spectrum of the data.
%   (Power spectrum is scaled so that maximum is 1.)
%
% In:
%   data: nx1 vector of data, 
%            where the time step is assumed to be 1 kyr.
% Out:
%   PS: nx1 vector giving the power at frequencies 
%       freq = (0:n-1)'/n (in units of 1/kyr)
 
n = length(data);
if n < 1
    PS = [];
    return
end

% ---- Detrend data.
ddata = detrend(data);

% ---- Compute power spectrum.
y = fft(ddata);
PS = y.*conj(y)/n;

% --- Normalize.
PS = PS/max(PS);

return
end
% ===================================================================

% ===================================================================
% =================== Budyko Model Computers ========================
% ===================================================================

% ============= Entry Point for Computation Function ================

function [eta, GMT] = computeBud( ecc, obl, Par )
% Compute the ice line and GMT for each value of eccentricity and 
%   obliquity.
%
% In:
%   ecc: mx1 vector of eccentricity values
%   obl: mx1 vector of corresponding obliquity values (radians)
%   Par(1):  A:      OLR = A + BT
%   Par(3):  B:
%   Par(4):  C:      coefficient of heat transport
%   Par(5):  alpha1: land/water albedo
%   Para(6): alpha2: ice albedo
%   Par(2):  Tc:     ice line temperature
%
% Out:
%   IL:  mx1 matrix of corresponding ice line locations
%   GMT: mx1 matrix of corresponding global mean temperatures

% --- Set Q0 so that current Q(e) is 343.
Q0 = 3.429521534956561e+002;

% --- define local variables to avoid confusion
A = Par(1);  Tc = Par(2);
B = Par(3);   C = Par(4);
alpha1 = Par(5);  alpha2 = Par(6);

% --- Construct the matrix H.
% Explanation:
%   If p is the row vector representing the polynomial h, i.e.
%      h(y) = p_1y^3 + p_2y^2 + p_3y + p_4 ,
%   then p is given by
%      p = [ 1 Q Qs_2 ] H ,
%   where Q is the GMT and s_2 is the coeff of the quadratic Legendre
%      polynomial in the expansion of s(y), the insolation distribution.

dela = alpha2 - alpha1;
    a0 = 0.5*(alpha2 + alpha1);

    c0 = (1 - a0)/B;
    c1 = C*dela/(2*B*(B+C));
    c2 = (1 - a0)/(2*(B+C));

    H = [ 0  ,  0   ,  0   , -A/B - Tc ; ...
          0  ,  0   , 2*c1 , c0 - c1 ;  ...
          c1 , 3*c2 , -c1  , -c2  ] ;

% --- Convert eccentricity to Q and obliquity to s2.
Q = Q0./sqrt( 1 - ecc.^2 );
s2 = -0.625 +  0.9375*sin(obl).^2 ;

% --- Compute the ice line.
eta = traceIL(H, Q, s2);

% --- Compute the global mean temperature.
S = eta + 0.5*s2.*eta.*(eta.^2 - 1);
GMT = ( Q.*( 1 - alpha2 + (alpha2 - alpha1)*S ) - A )/B;

return
end
% ===================================================================

% ============= Entry Point for the Tracing Function ================
% =============== (traces the stable equilibrium) ===================

function IL = traceIL( H , Q , s2 )
%Trace the ice line through various Milankovitch cycles.
%
% In: 
%   H: 3x4 matrix satisfying
%      ph = [1 Q Q*s2}*H, where ph is the polynomial representation
%         for the Widiasih h-function.
%   Q: column vector of global mean insolations
%   s2: column vector of quadratic Legendre coeffs of the insolation 
%       distribution function  s(y) .
%   Note: length(Q) = length(s2)
%
% Out:
%   IL: column vector of stable ice lines corresponding to Q and s2.
%
% External:
%   ClosestSRP( ph, eta0 )

IL = zeros(size(Q));
n = length(Q);
if n < 1
    return
end

q = [ ones(size(Q)), Q, Q.*s2 ];
eta0 = 1;
for k = 1:n
    eta = ClosestSRP( q(k,:)*H, eta0);
    IL(k) = eta;
    eta0 = eta;
end

return
end

% ============== Find the closest stable equilibrium =================

function eta = ClosestSRP( ph, eta0 )
%Find the stable restpoint for the quadratic Budyko system closest to eta0.
%
% The equation is reduced to the 1D:
%         \dot x = h(x)  on  [0,1]
% Returned is the stable rest point closest to eta0
%
% For boundary point 0:
%   0 is considered a rest point if h(0) <= 0. 
%        h'(0) = -Inf  if  h(0) < 0
%        h'(0) = h'(0) if  h(0) = 0
% For boundary point 1:
%   1 is considered a rest point if  h(1) >= 0.
%        h'(1) = +Inf  if  h(1) > 0
%        h'(1) = h'(1) if  h(1) = 0
%
% In:
%   ph is the MATLAB representation of h
%
% Out:
%   eta
%
% Note: If the equation has no stable rest point (shouldn't happen),
%         eta is empty.
%
% External:
%   RealRoots(p) (part of this file)
%   restrict(x)  (part of this file)

% --- Find the real roots of  h .
rp = RealRoots(ph);
% --- Note that the roots are in ascending order.

% --- Keep the roots in  (0,1)  and compute the derivatives.
rp = restrict(rp);
drp = polyval(polyder(ph),rp);

% --- Decide whether  0  is a restpoint.  Concatenate if so.
h0 = polyval(ph,0);
if h0 <= 0
    rp = [0, rp];
    if h0 < 0
        dh0 = -Inf;
    else
        dh0 = polyval(polyder(ph),0);
    end
    drp = [dh0, drp];
end

% --- Decide whether  1  is a restpoint.  Concatenate if so.
h1 = polyval(ph,1);
if h1 >= 0
    rp = [rp, 1];
    if h1 > 0
        dh1 = -Inf;
    else
        dh1 = polyval(polyder(ph),1);
    end
    drp = [drp, dh1];
end

% --- Now we have all the rest points in rp with the corresponding
%     eigenvalues in drp.

% --- Return if there are no rest points.  (Shouldn't happen.)
n = length(rp);
if n < 1
    eta = [];
    return
end

% --- Throw away the unstable rest points.
for k = 1:n
    i = n + 1 - k;
    if drp(i) > 0
        rp(i) = [];
    end
end

% --- Return if there are no stable rest points.  (Shouldn't happen.)
n = length(rp);
if n < 1
    eta = [];
    return
end

% --- If there is only one, it's the one.
if n == 1
    eta = rp;
    return
end

% --- Compute all the distances.
dist = abs(rp(1,:) - eta0);

% --- Find the smallest distance.
i = 1;
d = dist(1);
for k = 2:n
    if dist(k) < d;
        d = dist(k);
        i = k;
    end
end

% --- Done!
eta = rp(i);
return
end

% ===================== Miscellaneous Functions =====================

function rx = restrict(x)
%
% In:
%   x: a vector of real numbers in ascending order
%
% Out:
%   rx: a vector of those elements of x in the open interval (0,1)

n = length(x);
if n < 1
    rx = [];
    return
end

k1 = 0;
k2 = n+1;
for k = 1:n
    if x(k) <= 0
        k1 = k;
    elseif x(k) >=1
        k2 = k;
        break;
    end    
end

rx = x(k1+1:k2-1);
return

end

function rts = RealRoots( p )
%finds the real roots of a polynomial.
%
% In: polynomial p
%     p(1)x^n + p(2)x^{n-1} + ... + p(n)x + p(n+1)
%
% Out:
%   rts: real zeros of p (in increasing order)

r = roots(p);
deg = length(r);
if deg < 1
    rts = [];
    return
end

isre = zeros(1,deg);
for i = 1:deg
    isre(i) = isreal(r(i));
end

numr = sum(isre);
if numr < 1
    rts = [];
    return
end

rer = zeros(1,numr);
k = 0;
for i = 1:deg
    if isre(i)
        k = k+1;
        rer(k) = r(i);
    end
end

rts = sort(rer);

return
end

% ===================================================================
% ========================= End of Everything =======================


