function carbon_cycles()
% solves equations for albedo and CO2 partial pressure
%   da/dt = B(T) - a
%   dp/dt = alpha*(1-w*p^mu*e^T)
%   T = T(a,p)
% plots trajectories on phase plane
%
% IJH 28 Sept 2016

%% Dimensionless parameters
ap = 0.58;
am = 0.11;
c1 = 0.2;
c2 = 0.6;
mu = 0.3;
q = 1.37;
lambda = 0.25;
nu = 0.18;
alpha = .01; % 1.05 / .01
w = 1; % 1 / .15

%% Scaling parameters (for plotting dimesional quantity)
sc.p = 36; % [Pa]
sc.p = 1; % to plot as dimensionless

%% Functions
B = @(theta) ap-(ap-am)/2*(1+tanh(c1+c2*theta));
Theta = @(a,p) q/nu*(1-a)-1/nu+lambda*p; 
dadt = @(a,p) B(Theta(a,p)) - a;
dpdt = @(a,p) alpha*(1-w*p.^mu.*exp(Theta(a,p)));
dydt = @(t,y) [ dadt(y(1),y(2)) ; dpdt(y(1),y(2)) ]; % y = [a,p]

%% Nullclines
a = linspace(am,ap,100+2); a = a(2:end-1);
p = NaN*a;
for i = 1:length(a)
    p(i) = fzero(@(x) dadt(a(i),x),[-1e3 1e3]);     % root-finding for each a along a nullcline
end
a1 = a; p1 = p;
p = NaN*a;
for i = 1:length(a)
    p(i) = fzero(@(x) dpdt(a(i),x),[0 1e3]);        % root-finding for each a along p nullcline
end
a2 = a; p2 = p;

%% Plot phase plane with nullclines and setup axes for solutions
figure(1); clf; width = 2; fontsize = 14; 
set(gcf,'DefaultAxesFontSize',fontsize,'DefaultTextFontSize',fontsize);
% set(gcf,'units','centimeters','Paperpositionmode','auto','position',[2 2 20 10]);
ax(1) = axes('position',[0.1 0.2 0.35 0.75]);
    plot(sc.p*p1,a1,'k-',sc.p*p2,a2,'b-','linewidth',width);    % nullclines
    xlabel('$p / [p]$','interpreter','latex'); 
    ylabel('$a$','interpreter','latex');
    ylim([am,ap]);
    xlim([0 2*sc.p]);
    hold on;
ax(2) = axes('position',[0.6 0.6 0.35 0.35]);
    hold on;
    box on;
    set(gca,'XTickLabel',{});
    ylabel('$a$','interpreter','latex');
ax(3) = axes('position',[0.6 0.2 0.35 0.35]);
    hold on;
    box on;
    xlabel('$t / [t]$','interpreter','latex'); 
    ylabel('$p / [p]$','interpreter','latex');

%% Solve trajectories and add to figures
n_trajectories = 1;
length_trajectory = 10/alpha;
for n = 1:n_trajectories
    p0 = 2*rand; a0 = am+(ap-am)*rand;                          % random starting point
%     [p0,a0] = ginput(1); p0 = p0/sc.p;                          % graphical input for trajectory starting point (click on phase plane)
    opts = odeset('reltol',1e-6);
    [t,y] = ode45(dydt,[0 length_trajectory],[a0,p0],opts);     % use ode45 to solve equations
    a = y(:,1); p = y(:,2);
    axes(ax(1));
    plot(sc.p*p,a,'r-',sc.p*p(1),a(1),'ro','linewidth',width);   
    axes(ax(2));
    plot(t,a,'r-','linewidth',width);
    axes(ax(3));
    plot(t,sc.p*p,'r-','linewidth',width);
end

% %% Save figure
% axes(ax(1)); text(0.95,0.95,['$\alpha = $',num2str(alpha)],'units','normalized','horizontalalignment','right','verticalalignment','top','interpreter','latex');
% print(gcf,'-depsc2',[mfilename,'_alpha',num2str(alpha),'.eps'],'-loose');

end