GOPTIONS lfactor=10 hsize=6 in vsize=6 in horigin=1 in vorigin=2 in;
options nocenter;
TITLE1 f=swissb h=1.6 c=black 'The Lorenz 3D Chaos Model';

* Copyright (c) 2006 Courtney Brown. All rights reserved.;
* May be used freely for noncommercial purposes.;

PROC IML;

b=2.66666666;w=28.0; s=10.0;pi=3.14159;cycles=10;
x1=5.0;x2=5.0;x3=5.0;zing=0;
h=0.01;

start;
goto buildit;

RK4:
* Fourth Order Runge-Kutta;
* You can change iterations to 1024;

do kkk=1 to 2024;

x=uniform(zing);
y=sin((kkk/2024)#(cycles#2#pi));

m1=x1;m2=x2;m3=x3;
LINK EQS;
x1K1=Dx1DT;x2K1=Dx2DT;x3K1=Dx3DT;
m1=x1+(.5#h#x1K1);m2=x2+(.5#h#x2K1);m3=x3+(.5#h#x3K1);
LINK EQS;
x1K2=Dx1DT;x2K2=Dx2DT;x3K2=Dx3DT;
m1=x1+(.5#h#x1K2);m2=x2+(.5#h#x2K2);m3=x3+(.5#h#x3K2);
LINK EQS;
x1K3=Dx1DT;x2K3=Dx2DT;x3K3=Dx3DT;
m1=x1 + h#x1K3;m2=x2 + h#x2K3;m3=x3 + h#x3K3;
LINK EQS;
x1K4=Dx1DT;x2K4=Dx2DT;x3K4=Dx3DT;

x1NEW=x1+((h/6)#(x1K1+(2#x1K2)+(2#x1K3)+x1K4));
x2NEW=x2+((h/6)#(x2K1+(2#x2K2)+(2#x2K3)+x2K4));
x3NEW=x3+((h/6)#(x3K1+(2#x3K2)+(2#x3K3)+x3K4));

x1E=x1E//x1;x2E=x2E//x2;x3E=x3E//x3;xe=xe//x;ye=ye//y;
trajects=x1E||(x2E||(x3E||(xe||ye)));
x1=x1NEW;x2=x2NEW;x3=x3NEW;
end;
RETURN;

EQS:
Dx1DT = s#(m2 - m1);
Dx2DT = (w#m1) - m2 - (m1#m3);
Dx3DT = (m1#m2) - (b#m3);
RETURN;

BUILDIT:
LINK RK4;

x1max=max(x1E);x1min=min(x1E);
x2max=max(x2E);x2min=min(x2E);
x3max=max(x3E);x3min=min(x3E);
xmax=max(xe);xmin=min(xe);
ymax=max(ye);ymin=min(ye);
origin=x1min||(x2min||(x3min||(xmin||ymin)));
x1next=x1max||(x2min||(x3min||(xmin||ymin)));
x2next=x1min||(x2max||(x3min||(xmin||ymin)));
x3next=x1min||(x2min||(x3max||(xmin||ymin)));
final={99}||({99}||({99}||({99}||{99})));
axes=origin//(x1next//(origin//(x2next//(origin//(x3next//final)))));
trajects=axes//trajects;

*print trajects;
party={'X1' 'X2' 'X3' 'RANDOM' 'SINE'};
create traject from trajects (|colname=party|);
append from trajects;
close traject;

*print x1e x2e x3e;

*print trajects;

finish;run;

data traject;set traject;
perspect=1;pi=constant('PI');angle=2*pi/3;
if ((x3=99) and (x1=99)) then x1=.;
* 3-D projections into 2-D ... Two Versions;
* Version 1;
x1D3 = x1 - (perspect*x3);
x2D3 = x2 - (perspect*x3);
* Version 2;
x1D3 = x1 - ((cos(angle))*x3);
x2D3 = x2 - ((sin(angle))*x3);

sym=1;
*IF (T LT 250) THEN DELETE;

label x1D3 ='X1';
label x2D3 ='X2';
label X3 ='X3';

symbol1 color=black v=none f=simplex i=join;
symbol2 color=black f=simplex v='*';
symbol3 color=black f=simplex v='.';
symbol4 color=black f=simplex v='.';
symbol5 color=black f=simplex v='.';
proc gplot data=traject;
plot x2D3*x1D3=sym / nolegend skipmiss
noaxes;

data traject;set traject(firstobs=8);

proc spectra p s out=b(obs=2024) adjmean whitetest;
var x1 x2 x3 Random Sine;
weights 0 1 0;
data b;set b;
freq=freq/(2*3.14159);
label freq='Frequency in cycles per observation';

TITLE1 f=swissb h=1.6 c=black 'Fourier Analysis';
proc gplot data=b;
axis3 color=black
value=(h=1.5 f=swissb c=black)
label=(a=90 r=0 h=2 f=swissb c=black);
axis4 color=black
value=(h=1.5 f=swissb c=black)
label=(h=2 f=swissb c=black);
plot (p_01 p_02 p_03 p_04 p_05)*(freq)/vaxis=axis3 haxis=axis4
vminor=0 hminor=0;

proc print data=b;
var p_01 p_02 p_03 p_04 p_05 freq period;

run;
quit;