%%    UMD AOSC 615   SPRING 2009
%%    Lorenz 40 variable model ('95)

%% startup
clear all
close all
cpu0=cputime;
disp(sprintf('\n lorenz95: %s',datestr(now)))

%% global
global ode_prmtr

%% input
% parameter
F=8;
% number of variables
N=40;
% size of initial perturpation
sgm_x0=0.1;
% integration time
Tsu=100;   %  initial spin-up
Tot=10;       %  integration time
dTot=0.1;   %  output
% plot
ixn=[1 11];
c_x='r';
c_tx=['g','b'];
axs_tx=[0 Tot -5 15];
axs_x=[0 N -5 15];
lw_tx=2;

%% setup
% initial condition         
x0=sgm_x0*randn(N,1);
% output 
tot=[0:dTot:Tot];
KTot=length(tot);
xot=zeros(N,KTot);
% model parameter for ode
ode_prmtr=[F];    

%% initial spinup
[tsu,xsu]=ode45(@lorenz95_dxdt,[0 Tsu],x0);
xot(:,1)=xsu(end,:)'; 

%% model forward integration
for k=2:KTot;  kn1=k-1; xk0=xot(:,kn1);
    [tk,xk]=ode45(@lorenz95_dxdt,[tot(kn1) tot(k)],xk0);
    xot(:,k)=xk(end,:)'; 
end

%% plot
% t and xk
figure
for n=1:length(ixn);
    plot(tot,xot(ixn(n),:),c_tx(n),'LineWidth',lw_tx); if(n==1); hold on; end
end
legend(['x(',int2str(ixn(1)),')'],['x(',int2str(ixn(2)),')'],'Location','NorthEast')
axis(axs_tx)
xlabel('time'); ylabel('x(n)')

%% xk
figure
for k=1:KTot
    plot([0:1:N],[xot(N,k);xot(:,k)],c_x,'LineWidth',lw_tx);  hold on
    for n=1:length(ixn)
        plot(ixn(n),axs_x(3),'o','MarkerSize',8,'MarkerFaceColor',c_tx(n),'MarkerEdgeColor',c_tx(n))
    end        
    axis(axs_x)
    xlabel('n'); ylabel('x(n)')
    drawnow; hold off
end

%% finish
disp(sprintf(' ... end of lorenz95 (cpu=%s sec)\n',num2str(cputime-cpu0)))
