%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .% 
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.%
%.*.                                                                   .*.% 
%.*.                Program to Simulate the Lorenz-Attraktor           .*.% 
%.*.                   (Patrick Scholz, AWI, 06.01.2010)               .*.%
%.*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .*.% 
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.%


%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
%PARAMETER LORENTZ-ATRAKTOR
s=10;
r=28; %28.0;
b=8/3;

%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
% TIME STEP PARAMETER
deltat=0.001; %deltat
t_step=10000; % number of time steps
%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
% INITIAL POINTS X0;Y0;Z0
x0=10;
y0=20;
z0=0.3;

%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
% EULER METHOD INITIAL-CONDITION
x_eul=zeros(t_step+1,1)*NaN;         %allocate x vector euler
y_eul=zeros(t_step+1,1)*NaN;         %allocate y vector euler
z_eul=zeros(t_step+1,1)*NaN;         %allocate z vector euler
x_eul(1)=x0;y_eul(1)=y0;z_eul(1)=z0; %set initial condition euler

%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
% RUNGE KUTTA METHOD INITIAL-CONDITION
x_rkm=zeros(t_step+1,1)*NaN;         %allocate x vector runge-kutta
y_rkm=zeros(t_step+1,1)*NaN;         %allocate y vector runge-kutta
z_rkm=zeros(t_step+1,1)*NaN;         %allocate z vector runge-kutta
x_rkm(1)=x0;y_rkm(1)=y0;z_rkm(1)=z0; %set initial condition runge-kutta

%---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
% SETUP DIFFERENTIAL-EQUATION-SYSTEM
DGL_RHS_x= @(x,y,z,s,r,b) s*( -x + y);
DGL_RHS_y= @(x,y,z,s,r,b) ( r*x - y - x*z);
DGL_RHS_z= @(x,y,z,s,r,b) ( x*y - b*z );
 


for tt=2:t_step+1
    
    if mod(tt,500)==0;fprintf('timestep = %i\n',tt);end
    %---------------------------------------------------------------------%
    % EULER METHOD
    x_eul(tt)=x_eul(tt-1) + deltat * DGL_RHS_x(x_eul(tt-1),y_eul(tt-1),z_eul(tt-1),s,r,b);
    y_eul(tt)=y_eul(tt-1) + deltat * DGL_RHS_y(x_eul(tt-1),y_eul(tt-1),z_eul(tt-1),s,r,b);
    z_eul(tt)=z_eul(tt-1) + deltat * DGL_RHS_z(x_eul(tt-1),y_eul(tt-1),z_eul(tt-1),s,r,b);
    
    %break off simulation when values of x,y,z run towards infinity
    if (x_eul(tt)>1e20) || (y_eul(tt)>1e20) || (z_eul(tt)>1e20)
        break
    end
    %---------------------------------------------------------------------%
    % RUNGE KUTTA METHOD
   
    k1_x = DGL_RHS_x(x_rkm(tt-1)               , y_rkm(tt-1)               , z_rkm(tt-1)              ,s,r,b);
    k2_x = DGL_RHS_x(x_rkm(tt-1)+deltat*k1_x/2 , y_rkm(tt-1)+deltat*k1_x/2 , z_rkm(tt-1)+deltat*k1_x/2,s,r,b);
    k3_x = DGL_RHS_x(x_rkm(tt-1)+deltat*k2_x/2 , y_rkm(tt-1)+deltat*k2_x/2 , z_rkm(tt-1)+deltat*k2_x/2,s,r,b);
    k4_x = DGL_RHS_x(x_rkm(tt-1)+deltat*k3_x   , y_rkm(tt-1)+deltat*k3_x   , z_rkm(tt-1)+deltat*k3_x  ,s,r,b);
    x_rkm(tt)=x_rkm(tt-1) + deltat*( 1/6*k1_x + 1/3*k2_x + 1/3*k3_x + 1/6*k4_x);
    
    k1_y = DGL_RHS_y(x_rkm(tt-1)               ,y_rkm(tt-1)                , z_rkm(tt-1)              ,s,r,b);
    k2_y = DGL_RHS_y(x_rkm(tt-1)+deltat*k1_y/2 , y_rkm(tt-1)+deltat*k1_y/2 , z_rkm(tt-1)+deltat*k1_y/2,s,r,b);
    k3_y = DGL_RHS_y(x_rkm(tt-1)+deltat*k2_y/2 , y_rkm(tt-1)+deltat*k2_y/2 , z_rkm(tt-1)+deltat*k2_y/2,s,r,b);
    k4_y = DGL_RHS_y(x_rkm(tt-1)+deltat*k3_y   , y_rkm(tt-1)+deltat*k3_y   , z_rkm(tt-1)+deltat*k3_y  ,s,r,b);
    y_rkm(tt)=y_rkm(tt-1) + deltat*( 1/6*k1_y + 1/3*k2_y + 1/3*k3_y + 1/6*k4_y);
    
    k1_z = DGL_RHS_z(x_rkm(tt-1)               , y_rkm(tt-1)               , z_rkm(tt-1)              ,s,r,b);
    k2_z = DGL_RHS_z(x_rkm(tt-1)+deltat*k1_z/2 , y_rkm(tt-1)+deltat*k1_z/2 , z_rkm(tt-1)+deltat*k1_z/2,s,r,b);
    k3_z = DGL_RHS_z(x_rkm(tt-1)+deltat*k2_z/2 , y_rkm(tt-1)+deltat*k2_z/2 , z_rkm(tt-1)+deltat*k2_z/2,s,r,b);
    k4_z = DGL_RHS_z(x_rkm(tt-1)+deltat*k3_z   , y_rkm(tt-1)+deltat*k3_z   , z_rkm(tt-1)+deltat*k3_z  ,s,r,b);
    z_rkm(tt)=z_rkm(tt-1) + deltat*( 1/6*k1_z + 1/3*k2_z + 1/3*k3_z + 1/6*k4_z);
    
    %break off simulation when values of x,y,z run towards infinity
    if (x_rkm(tt)>1e20) || (y_rkm(tt)>1e20) || (z_rkm(tt)>1e20)
        break
    end
end

hcp=[];hcm=[];h_rkm=[];h_eul=[];h_end=[];h_start=[];
hf=figure;
    %---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
    ax1=subplot(1,2,1);
    hold on
        h_eul=plot3(x_eul,y_eul,z_eul,'b') ;
        h_end=plot3(x_eul(end),y_eul(end),z_eul(end),'ob','MarkerFaceColor','c','MarkerSize',8); %Endpunkte
        h_start=plot3(x_eul(1),y_eul(1),z_eul(1),'oc','MarkerFaceColor','b','MarkerSize',7); %Startpunkte
        %---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
        % 1) Fixpoint (Gleichgewichtspunkte) --> (x,y,z)=(0,0,0)
        hcr=plot3(0,0,0,'or','MarkerFaceColor','r','MarkerSize',6); %Ruhelage
        if r>1
            % 2) & 3) Analytical calculated Fixpoint (Gleichgewichtspunkte)
            % C_+  = (x_+ , y_+ , z_+) = ( sqrt(b*(r-1)) ,  sqrt(b*(r-1)) , r-1) 
            % C_+  = (x_- , y_- , z_-) = (-sqrt(b*(r-1)) , -sqrt(b*(r-1)) , r-1) 
            hcp=plot3(sqrt(b*(r-1)),sqrt(b*(r-1)),r-1,'vr','MarkerFaceColor',[255,153,0]/255,'MarkerSize',8); %Fixpunkt C+
            hcm=plot3(-sqrt(b*(r-1)),-sqrt(b*(r-1)),r-1,'^r','MarkerFaceColor',[255,153,0]/255,'MarkerSize',8); %Fixpunkt C-
        end
    hold off 
    view(3)    
    xlabel('X')
    ylabel('Y')
    zlabel('Z')
    title('EULER-METHODE')
    box on 
    grid on 
    if isempty(hcp)==1 || isempty(hcm)==1
         leg_eul=legend([h_eul h_end h_start hcr],'Euler','Endpunkt','Startpunkt','Ruhelage (0,0,0)',...
                'Location','North');...        'Orientation','horizontal')
    else
         leg_eul=legend([h_eul h_end h_start hcr hcp hcm],'Euler','Endpunkt','Startpunkt','Ruhelage (0,0,0)','Fixpunkt C^{+}','Fixpunkt C^{-}',...
                'Location','North');...        'Orientation','horizontal')   
    end
    %---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
    ax2=subplot(1,2,2);
    hold on
        h_rkm=plot3(x_rkm,y_rkm,z_rkm,'g') ;
        h_end=plot3(x_rkm(end),y_rkm(end),z_rkm(end),'og','MarkerFaceColor','y','MarkerSize',8); %Endpunkte
        h_start=plot3(x_rkm(1),y_rkm(1),z_rkm(1),'oy','MarkerFaceColor','g','MarkerSize',7); %Startpunkte
        %---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
        % 1) Fixpoint (Gleichgewichtspunkte) --> (x,y,z)=(0,0,0)
        hcr=plot3(0,0,0,'or','MarkerFaceColor','r','MarkerSize',6); %Ruhelage
        if r>1
            % 2) & 3) Analytical calculated Fixpoint (Gleichgewichtspunkte)
            % C_+  = (x_+ , y_+ , z_+) = ( sqrt(b*(r-1)) ,  sqrt(b*(r-1)) , r-1) 
            % C_+  = (x_- , y_- , z_-) = (-sqrt(b*(r-1)) , -sqrt(b*(r-1)) , r-1) 
            hcp=plot3(sqrt(b*(r-1)),sqrt(b*(r-1)),r-1,'vr','MarkerFaceColor',[255,153,0]/255,'MarkerSize',8); %Fixpunkt C+
            hcm=plot3(-sqrt(b*(r-1)),-sqrt(b*(r-1)),r-1,'^r','MarkerFaceColor',[255,153,0]/255,'MarkerSize',8); %Fixpunkt C-
        end
    hold off 
    view(3)    
    xlabel('X')
    ylabel('Y')
    zlabel('Z')
    title('RUNGE-KUTTA-METHODE')
    box on 
    grid on
    if isempty(hcp)==1 || isempty(hcm)==1
         leg_rkm=legend(ax2,[h_rkm h_end h_start hcr],'Runge-Kutta','Endpunkt','Startpunkt','Ruhelage (0,0,0)',...
                'Location','North');...        'Orientation','horizontal')
    else
         leg_rkm=legend(ax2,[h_rkm h_end h_start hcr hcp hcm],'Runge-Kutta','Endpunkt','Startpunkt','Ruhelage (0,0,0)','Fixpunkt C^{+}','Fixpunkt C^{-}',...
                'Location','North');...        'Orientation','horizontal')   
    end
    %---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|--
set(leg_eul,'FontSize',8);
set(leg_rkm,'FontSize',8);

disp('FINISH')