(* 
   Mathematica program that uses a 4  th - 
          order Runge - 
  Kutta method to numerically approximate points on 
the Lorenz attractor . 
                by Lewis Dartnell
*)

<< 
          Graphics`Graphics3D`;


      Sigma  =  10;

      r  =  28;

      b  =  8/3;

      h = 0.01;


      CurrentVector = {-2.5760498013251096`, -4.369978021478614`, 
            13.840153578911828`};

      SavedData  =  {};


      CalculateA[{x_, y_, 
              z_}] := 
(
VectorA = {Sigma (y - x), r x  -  y  -  x z, 
                  x y  -  
                    b z};
);


      CalculateB[{x_, y_, 
              z_}] := 
(
TempBx = 
              x +  0.5 h VectorA[[1]]; 
TempBy = 
              y +  0.5 h VectorA[[2]]; 
TempBz = 
              z +  0.5 h VectorA[[3]]; 
VectorB = {Sigma (TempBy - TempBx),  
                r TempBx  -  TempBy  -  TempBx TempBz,  
                TempBx TempBy  -  
                  b TempBz};
);


      CalculateC[{x_, y_, 
              z_}] := 
(
TempCx = 
              x +  0.5 h VectorB[[1]]; 
TempCy = 
              y +  0.5 h VectorB[[2]]; 
TempCz = 
              z +  0.5 h VectorB[[3]]; 
VectorC = {Sigma (TempCy - TempCx),  
                r TempCx  -  TempCy  -  TempCx TempCz,  
                TempCx TempCy  -  
                  b TempCz};
);


      CalculateD[{x_, y_, 
              z_}] := 
(
TempDx = 
              x +   h VectorC[[1]]; 
TempDy = 
              y +   h VectorB[[2]]; 
TempDz = 
              z +   h VectorB[[3]]; 
VectorD = {Sigma (TempDy - TempDx),  
                r TempDx  -  TempDy  -  TempDx TempDz,  
                TempDx TempDy  -  
                  b TempDz};
);


 (* 
        Loop to reiterate the process on each new point,  
        and dump data to an array ...  *) 


      Do[

SavedData = 
            Append[SavedData, 
              CurrentVector]; 

CalculateA[{CurrentVector[[1]], 
              CurrentVector[[2]], 
              CurrentVector[[3]]}]; 
CalculateB[{CurrentVector[[1]], 
              CurrentVector[[2]], 
              CurrentVector[[3]]}]; 
CalculateC[{CurrentVector[[1]], 
              CurrentVector[[2]], 
              CurrentVector[[3]]}]; 
CalculateD[{CurrentVector[[1]], 
              CurrentVector[[2]], 
              CurrentVector[[3]]}]; 

NextVectorx = 
            CurrentVector[[1]] +  
              1/6 h (VectorA[[1]]  +  2 VectorB[[1]]  +  
                    2 VectorC[[1]]  +  
                    VectorD[[1]]); 
NextVectory = 
            CurrentVector[[2]] +  
              1/6 h (VectorA[[2]]  +  2 VectorB[[2]]  +  
                    2 VectorC[[2]]  +  
                    VectorD[[2]]); 
NextVectorz = 
            CurrentVector[[3]] +  
              1/6 h (VectorA[[3]]  +  2 VectorB[[3]]  +  
                    2 VectorC[[3]]  +  
                    VectorD[[3]]); 

CurrentVector = {NextVectorx, NextVectory, 
              NextVectorz};


, 
{5000}
] // 
        Timing

\(ScatterPlot3D[
          SavedData, \ 
PlotJoined -> 
            False, \ 
PlotStyle -> {GrayLevel[0.2], \ 
              Thickness[0.001], 
              PointSize[0.001]}, 
Boxed -> 
            False, 
Axes -> 
            True, 
AxesEdge -> {{\(-1\), \(-1\)}, {\(-1\), \(-1\)}, {\(-1\), \(-1\)}}, 
AxesLabel -> {x, 
              y, z}, \ 
Ticks -> 
            True, 
AxesStyle -> {RGBColor[0.2, 0, 
                0.8], Thickness[
                0.001]}, 
TextStyle -> {FontFamily -> "Arial", FontSize -> 12, \ 
              FontColor -> 
                RGBColor[0.2, 0, 
                  0.8]}, 
ViewPoint -> {20`, \(-20`\), 
              10.`}, 
ViewCenter -> {1/2, 1/2, 
              1/2}, 
SphericalRegion -> 
            True, 
ImageSize -> 800] // 
        Timing\)

\(\(Loop = 1;\)

      \(Radius = 10;\)


      \(Do[
LewisFrame = 
            ScatterPlot3D[
              SavedData, \ 
PlotJoined -> 
                False, \ 
PlotStyle -> {Thickness[
                    0.001], PointSize[0.001]}, 
Boxed -> 
                False, 
Axes -> 
                True, 
AxesEdge -> {{\(-1\), \(-1\)}, {\(-1\), \(-1\)}, {\(-1\), \(-1\)}}, 
 (*\ \(AxesLabel -> {x, y, z}\)\(,\)\ *) 
AxesLabel -> 
                None, 
Ticks -> 
                False, 
AxesStyle -> {RGBColor[0.2, 0, 
                    0.8], Thickness[
                    0.001]}, 
TextStyle -> {FontFamily -> "Arial", FontSize -> 14, \ 
                  FontColor -> 
                    RGBColor[0.2, 0, 
                      0.8]}, 
ViewPoint -> {Radius\ Cos[Theta], Radius\ Sin[Theta], 
                  10\ Sin[Theta/2]\  + \ 
                    10}, 
ViewCenter -> {1/2, 1/2, 
                  1/2}, 
SphericalRegion -> 
                True, 
ImageSize -> 
                300
]; 

Export["Lorenz" <> ToString[Loop] <> ".bmp", 
            LewisFrame, "BMP"]; 
Loop = 
            Loop + 1;

, 
{Theta, 0, 
            4\ Pi, \ \((2\ Pi)\)/
              40}
];\)

      \)

\(\(MaxLoop = 
          1000;\)


      \(Do[
LewisFrame = 
            ScatterPlot3D[
              Take[SavedData, Loop], \ 
PlotJoined -> 
                True, \ 
PlotStyle -> {Thickness[
                    0.001], PointSize[0.001]}, 
Boxed -> 
                False, 
Axes -> 
                True, 
AxesEdge -> {{\(-1\), \(-1\)}, {\(-1\), \(-1\)}, {\(-1\), \(-1\)}}, 
AxesLabel -> {x, 
                  y, z}, 
Ticks -> 
                False, 
AxesStyle -> {RGBColor[0.2, 0, 
                    0.8], Thickness[
                    0.001]}, 
TextStyle -> {FontFamily -> "Arial", FontSize -> 14, \ 
                  FontColor -> 
                    RGBColor[0.2, 0, 
                      0.8]}, 
ViewPoint -> {20`, \(-20`\), 10.`}, 
ViewCenter -> {1/2, 1/2, 
                  1/2}, 
SphericalRegion -> 
                True, 
ImageSize -> 
                300, 
PlotRange -> {{\(-20\), 
                    20}, {\(-25\), 25}, {0, 
                    45}}
]; 

If[
            Mod[Loop, 20] == 0\ \  || \ 
              Loop == 0\ , 
\(Export["Lorenz" <> 
                  ToString[Loop] <> ".bmp", 
                LewisFrame, "BMP"];\)
, 
]

, 
{Loop, \ 1, 
            MaxLoop}
];\)


      \)

ScatterPlot3D // Options