The following program was used to explore the highly nonlinear and potentially chaotic
models that are analyzed in chapter 6 of Serpents in the Sand (University of Michigan
Press 1995). The program is written in Think Pascal, and it was run on a Macintosh
Quadra 900. It is capable of producing postscript output that can be sent directly to
a postscript printer, or input into a word processing program capable of accepting
postscript files (such as Word). The program is copyrighted 1995 by Courtney Brown.
However, researchers may copy, paste, and share this program freely for their own
use as long as the program itself is never sold. Simply copy and paste it into your
own program editor. To understand a bit more about how the program works, I
recommend my book, Chaos and Catastrophe Theories (Sage 1995).

program chaos (input, output, Image1, Image2, Image3);


	const
		minwidth = 15;
		maxsize = 1000;
		f1name = 'ImageLyap';
{f1name = 'ImageSdev';}
		f2name = 'ImageSum';
		f3name = 'ImagePower';
{nn=# of equations, nnn=total # of equations;}
		matnn = 4;
		matnnn = 20;
		inc = 0.001;
		vv2 = 2048;
		nnpp = 9000;
		gleft = 680;
		gtop = 460;
		gright = 1080;
		gbottom = 860;

	type
		HugePtr = ^BigArray;
		BigArray = array[1..200, 1..200] of real;
		HugePtr2 = ^BigArray2;
		BigArray2 = array[1..20000, 1..matnn] of real;
		HugePtr3 = ^BigArray3;
		BigArray3 = array[1..10000] of real;
		GraphPtr = ^GraphBig;
		GraphBig = array[1..20000, 1..2] of real;
		RealArrayVV2 = array[1..vv2] of real;
		RealArrayVV3 = array[1..1025, 1..2] of real;
		RealArrayNP = array[1..nnpp] of real;

	var
		dim, c1, c2, d1, d2, iter, iternum, toss, lyap, kitu, limit, multiplier, alert, partychoice, choice, shocks: integer;
		funcdem, funcrep, N, x1, x2, z1, z2, h, thing, min, max, weigh, totvote, time, perspect, CINI: real;
		t, cumkt, upper, lower, decay, sense, limitpol, halflife, shift, policy, shock, varsum, varold, change: real;
		lx1, ux1, lx2, ux2, xprim, lowp1, lowp2, highp1, highp2, chaoslagg, volume, Et, mold, rDR, timeone: real;
		totderiv1, totderiv2, min2, max2, min3, max3, totderiv3, totpower, avefourier, varfourier, sdevfour: real;
		ff1, ff2, ff3: text;
		gizmo, demdum, repdum, row, col, ticker, ticker2, ticker3, lagg, yearslag, graphtype, traject, sasa: integer;
		politics, algorithm, nn, nnn, eqnum, laggedid, count, item, fourier, anza, graphonly, totsignal: integer;
		signal1, suppress, p1lag: integer;
		filein: text;
		y: HugePtr2;
		p: array[1..20] of real;
		party: array[1..4] of integer;
		image1: HugePtr;
		image2: HugePtr;
		image3: HugePtr;
		newx: GraphPtr;
		fourdata: RealArrayVV2;
		spectral: RealArrayVV3;
		xone: array[1..1000] of real;
		xlag: HugePtr3;
		mx: array[1..matnnn] of real;
		x: array[1..matnnn] of real;
		znorm: array[1..matnn] of real;
		gsc: array[1..matnn] of real;
		cum: array[1..matnn] of real;
		xnew: array[1..matnnn] of real;
		finame: string;
		r: Rect;
		r2: Rect;

	procedure initialize;
	begin
		r.left := 0;
		r.top := 50;
		r.right := 640;
		r.bottom := 460;
		SetTextRect(r);
		r2.left := gleft;
		r2.top := gtop;
		r2.right := gright;
		r2.bottom := gbottom;
		SetDrawingRect(r2);
		new(y);
		new(newx);
		new(xlag);
		alert := 0;
		ticker2 := 0;
		yearslag := 0;
		laggedid := 0;
	end;

	procedure parms;
	begin
		p[4] := 0.08847;
		p[2] := 0.23665;
		p[11] := 0.27780;
		p[1] := 0.11699;
		p[5] := 1.54509;
		p[3] := 0.54922;
		p[12] := 0.29619;
		p[8] := 0.87930;
		p[10] := 0.16780;
		p[9] := 0.09635;
		p[6] := 0.77792;
		p[7] := 0.50896;
		p[13] := 0.32081;
	end;

	function lag (lagvar: real): real;
		var
			iii: integer;
	begin
		if (ticker3 = 1) then
			begin
				for iii := 1 to lagg do
					begin
						xlag^[iii] := lagvar;
					end;
				lag := lagvar;
			end;
		if (ticker3 > 1) then
			begin
				for iii := (lagg - 1) downto 1 do
					begin
						xlag^[(lagg - iii)] := xlag^[(lagg - iii + 1)];
					end;
				xlag^[(lagg)] := lagvar;
				lag := xlag^[1];
			end;
	end;

	function power (base, exponent: real): real;
	begin
		if (base < 0) then
			writeln('The power exponential function cannot use a negative base argument!');
		if (base > 0) then
			power := exp(exponent * ln(base));
		if (base = 0) then
			power := 0.0;
	end;


	function eq1: real;
	begin
		if (kitu = 1) then
			eq1 := -(mx[1] / p[2]) - (sin(mx[2])) + (p[1] * cos(mx[3]));
		if (kitu = 2) then
			eq1 := p[1] * (mx[2] - mx[1]);
		if (kitu = 3) then
			begin
				limitpol := 1;
				if (politics = 2) then
					policy := p[2];
				if ((politics = 1) and (eqnum = 3)) then
					begin
						if (choice = 1) then
							policy := 0.0;
						if (choice = 0) then
							policy := p[9];
					end;
				if ((politics = 1) and (eqnum = 4)) then
					policy := mx[4];
				eq1 := ((p[1] * mx[1]) * (1 - mx[1])) - (policy * mx[2] * mx[1]) - (decay * mx[1]);
			end;
		if (kitu = 4) then
			eq1 := (p[1] * mx[1] * (1 - mold)) - (p[5] * mx[1] * mx[2]) - (p[2] * mx[1] * sqr(sin(mx[3])));
		if (kitu = 5) then
			eq1 := (((p[1] * mx[1]) + (p[2] * mx[2])) * (1 + (p[3] * CINI) + (p[13] * mx[3]))) * mx[1];
	end;

	function eq2: real;
	begin
		if (kitu = 1) then
			eq2 := mx[1];
		if (kitu = 2) then
			eq2 := (p[2] * mx[1]) - (mx[1] * mx[3]) - mx[2];
		if (kitu = 3) then
			eq2 := (p[3] * mold) * (1 - mx[2]) - (p[4] * mx[3]);
		if (kitu = 4) then
			eq2 := (p[4] * mx[1] * (1 - mx[2])) - (p[6] * mx[2] * (1 - mx[1]));
		if (kitu = 5) then
			eq2 := (((p[4] * mx[1]) - (p[5] * mx[1] * mx[3])) * (1 + (p[6] * mx[2]))) * mx[2];
	end;

	function eq3: real;
	begin
		if (kitu = 1) then
			eq3 := p[3];
		if (kitu = 2) then
			eq3 := (mx[1] * mx[2]) - (p[3] * mx[3]);
		if (kitu = 3) then
			eq3 := ((p[5] * mx[2] * mx[1]) * (1 - mx[3])) - (p[6] * mx[3]);
		if (kitu = 4) then
			eq3 := (2 * pi);
		if (kitu = 5) then
			eq3 := (p[7] + (p[8] * ln(abs((p[9] + (p[10] * mx[1]) + (p[11] * mx[2])) / mx[3]))) - (p[12] * mx[3])) * mx[3];
	end;

	function eq4: real;
		var
			partydum1, partydum2: real;
	begin
		if (kitu = 3) then
			begin
				if (choice = 1) then
					begin
						partydum1 := 1.0;
						partydum2 := 0.0;
					end;
				if (choice = 0) then
					begin
						partydum1 := 0.0;
						partydum2 := 1.0;
					end;
				eq4 := p[2] * mx[4] * ((p[7] * partydum1) + (p[8] * partydum2) - mx[4]);
			end;
	end;

	procedure special3;
		var
			iii, tround: integer;
	begin
		if (timeone > 4.0) then
			begin
				partychoice := partychoice + 1;
				party[1] := 1;
				party[2] := 1;
				party[3] := 0;
				party[4] := 0;
				choice := party[partychoice];
				if (partychoice = 4) then
					partychoice := 0;
				timeone := 0.0;
			end;
		mold := lag(x[1]);
		laggedid := 1;
		Et := 0.1 * (Random / 32768.0);
{if (Et > 0) then Et := 0;}
{writeln('Et=', Et);}
{writeln('x1=', x[1] : minwidth, '   m1old= ', m1old : minwidth, 'lag= ', lagg);}
		tround := round(time);
		if ((shocks = 1) and (tround = (toss div 2))) then
			x[1] := x[1] + shock;
	end;

	procedure special4;
	begin
		mold := lag(x[1]);
		laggedid := 1;
	end;

	procedure special5;
	begin
		if (ticker3 = 1) then
			begin
				x[1] := 0.0;
				x[2] := 0.86;
				x[3] := 0.02;
				CINI := 0.23;
			end;
		if (ticker3 = 11) then
			begin
				x[1] := 0.0;
				CINI := 0.28;
			end;
		if ((ticker3 > 20) and (x[1] > 0.25)) then
			begin
				x[1] := 0.0;
				CINI := 0.25;
			end;
	end;

	function par1x1: real;
		var
			it1, it2: real;
	begin
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
		it1 := eq1;
		mx[1] := mx[1] + (2 * inc);
		if (laggedid = 1) then
			mold := mold + (2 * inc);
		it2 := eq1;
		par1x1 := (it2 - it1) / (2 * inc);
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
	end;

	function par1x2: real;
		var
			it1, it2: real;
	begin
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
		it1 := eq1;
		mx[2] := mx[2] + (2 * inc);
		if (laggedid = 2) then
			mold := mold + (2 * inc);
		it2 := eq1;
		par1x2 := (it2 - it1) / (2 * inc);
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
	end;

	function par1x3: real;
		var
			it1, it2: real;
	begin
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
		it1 := eq1;
		mx[3] := mx[3] + (2 * inc);
		if (laggedid = 3) then
			mold := mold + (2 * inc);
		it2 := eq1;
		par1x3 := (it2 - it1) / (2 * inc);
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
	end;

	function par1x4: real;
		var
			it1, it2: real;
	begin
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
		it1 := eq1;
		mx[4] := mx[4] + (2 * inc);
		if (laggedid = 4) then
			mold := mold + (2 * inc);
		it2 := eq1;
		par1x4 := (it2 - it1) / (2 * inc);
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
	end;

	function par2x1: real;
		var
			it1, it2: real;
	begin
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
		it1 := eq2;
		mx[1] := mx[1] + (2 * inc);
		if (laggedid = 1) then
			mold := mold + (2 * inc);
		it2 := eq2;
		par2x1 := (it2 - it1) / (2 * inc);
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
	end;

	function par2x2: real;
		var
			it1, it2: real;
	begin
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
		it1 := eq2;
		mx[2] := mx[2] + (2 * inc);
		if (laggedid = 2) then
			mold := mold + (2 * inc);
		it2 := eq2;
		par2x2 := (it2 - it1) / (2 * inc);
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
	end;

	function par2x3: real;
		var
			it1, it2: real;
	begin
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
		it1 := eq2;
		mx[3] := mx[3] + (2 * inc);
		if (laggedid = 3) then
			mold := mold + (2 * inc);
		it2 := eq2;
		par2x3 := (it2 - it1) / (2 * inc);
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
	end;

	function par2x4: real;
		var
			it1, it2: real;
	begin
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
		it1 := eq2;
		mx[4] := mx[4] + (2 * inc);
		if (laggedid = 4) then
			mold := mold + (2 * inc);
		it2 := eq2;
		par2x4 := (it2 - it1) / (2 * inc);
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
	end;

	function par3x1: real;
		var
			it1, it2: real;
	begin
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
		it1 := eq3;
		mx[1] := mx[1] + (2 * inc);
		if (laggedid = 1) then
			mold := mold + (2 * inc);
		it2 := eq3;
		par3x1 := (it2 - it1) / (2 * inc);
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
	end;

	function par3x2: real;
		var
			it1, it2: real;
	begin
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
		it1 := eq3;
		mx[2] := mx[2] + (2 * inc);
		if (laggedid = 2) then
			mold := mold + (2 * inc);
		it2 := eq3;
		par3x2 := (it2 - it1) / (2 * inc);
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
	end;

	function par3x3: real;
		var
			it1, it2: real;
	begin
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
		it1 := eq3;
		mx[3] := mx[3] + (2 * inc);
		if (laggedid = 3) then
			mold := mold + (2 * inc);
		it2 := eq3;
		par3x3 := (it2 - it1) / (2 * inc);
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
	end;

	function par3x4: real;
		var
			it1, it2: real;
	begin
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
		it1 := eq3;
		mx[4] := mx[4] + (2 * inc);
		if (laggedid = 4) then
			mold := mold + (2 * inc);
		it2 := eq3;
		par3x4 := (it2 - it1) / (2 * inc);
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
	end;

	function par4x1: real;
		var
			it1, it2: real;
	begin
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
		it1 := eq4;
		mx[1] := mx[1] + (2 * inc);
		if (laggedid = 1) then
			mold := mold + (2 * inc);
		it2 := eq4;
		par4x1 := (it2 - it1) / (2 * inc);
		mx[1] := mx[1] - inc;
		if (laggedid = 1) then
			mold := mold - inc;
	end;

	function par4x2: real;
		var
			it1, it2: real;
	begin
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
		it1 := eq4;
		mx[2] := mx[2] + (2 * inc);
		if (laggedid = 2) then
			mold := mold + (2 * inc);
		it2 := eq4;
		par4x2 := (it2 - it1) / (2 * inc);
		mx[2] := mx[2] - inc;
		if (laggedid = 2) then
			mold := mold - inc;
	end;

	function par4x3: real;
		var
			it1, it2: real;
	begin
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
		it1 := eq4;
		mx[3] := mx[3] + (2 * inc);
		if (laggedid = 3) then
			mold := mold + (2 * inc);
		it2 := eq4;
		par4x3 := (it2 - it1) / (2 * inc);
		mx[3] := mx[3] - inc;
		if (laggedid = 3) then
			mold := mold - inc;
	end;

	function par4x4: real;
		var
			it1, it2: real;
	begin
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
		it1 := eq4;
		mx[4] := mx[4] + (2 * inc);
		if (laggedid = 4) then
			mold := mold + (2 * inc);
		it2 := eq4;
		par4x4 := (it2 - it1) / (2 * inc);
		mx[4] := mx[4] - inc;
		if (laggedid = 4) then
			mold := mold - inc;
	end;

	procedure equations (k: integer);
		var
			i: integer;
	begin
		if (eqnum = 3) then
			begin
				if (k = 1) then
					xprim := eq1;
				if (k = 2) then
					xprim := eq2;
				if (k = 3) then
					xprim := eq3;
				if (k > 3) then
					begin
						if (k < 7) then
							begin
								i := k - 4;
								xprim := (mx[4 + i] * par1x1) + (mx[7 + i] * par1x2) + (mx[10 + i] * par1x3);
							end;
					end;
				if (k > 6) then
					begin
						if (k < 10) then
							begin
								i := k - 7;
								xprim := (mx[4 + i] * par2x1) + (mx[7 + i] * par2x2) + (mx[10 + i] * par2x3);
							end;
					end;
				if (k > 9) then
					begin
						i := k - 10;
						xprim := (mx[4 + i] * par3x1) + (mx[7 + i] * par3x2) + (mx[10 + i] * par3x3);
					end;
			end;
		if (eqnum = 4) then
			begin
				if (k = 1) then
					xprim := eq1;
				if (k = 2) then
					xprim := eq2;
				if (k = 3) then
					xprim := eq3;
				if (k = 4) then
					xprim := eq4;
				if (k > 4) then
					begin
						if (k < 9) then
							begin
								i := k - 5;
								xprim := (mx[5 + i] * par1x1) + (mx[9 + i] * par1x2) + (mx[13 + i] * par1x3) + (mx[17] * par1x4);
							end;
					end;
				if (k > 8) then
					begin
						if (k < 13) then
							begin
								i := k - 9;
								xprim := (mx[5 + i] * par2x1) + (mx[9 + i] * par2x2) + (mx[13 + i] * par2x3) + (mx[17] * par2x4);
							end;
					end;
				if (k > 12) then
					begin
						if (k < 17) then
							begin
								i := k - 13;
								xprim := (mx[5 + i] * par3x1) + (mx[9 + i] * par3x2) + (mx[13 + i] * par3x3) + (mx[17] * par3x4);
							end;
					end;
				if (k > 16) then
					begin
						i := k - 17;
						xprim := (mx[5 + i] * par4x1) + (mx[9 + i] * par4x2) + (mx[13 + i] * par4x3) + (mx[17] * par4x4);
					end;
			end;
	end;

	procedure model;
		var
			k: integer;
			rk1, rk2, rk3, rk4: array[1..matnnn] of real;
			xx1, xx2, xx3, xx4: array[1..matnnn] of real;
	begin
{The fourth order Runge-Kutta;}
		for k := 1 to nnn do
			begin
				mx[k] := x[k];
			end;
		for k := 1 to nnn do
			begin
				equations(k);
				rk1[k] := h * xprim;
			end;
		for k := 1 to nnn do
			begin
				xx1[k] := x[k] + (rk1[k] * 0.5);
			end;
		for k := 1 to nnn do
			begin
				mx[k] := xx1[k];
			end;
		for k := 1 to nnn do
			begin
				equations(k);
				rk2[k] := h * xprim;
			end;
		for k := 1 to nnn do
			begin
				xx2[k] := x[k] + (rk2[k] * 0.5);
			end;
		for k := 1 to nnn do
			begin
				mx[k] := xx2[k];
			end;
		for k := 1 to nnn do
			begin
				equations(k);
				rk3[k] := h * xprim;
			end;
		for k := 1 to nnn do
			begin
				xx3[k] := x[k] + rk3[k];
			end;
		for k := 1 to nnn do
			begin
				mx[k] := xx3[k];
			end;
		for k := 1 to nnn do
			begin
				equations(k);
				rk4[k] := h * xprim;
			end;
		for k := 1 to nnn do
			begin
				xnew[k] := x[k] + ((1 / 6) * (rk1[k] + (2 * rk2[k]) + (2 * rk3[k]) + rk4[k]));
			end;
	end;

	procedure initcond;
		var
			ii: integer;
	begin
{Initial conditions for the nonlinear system;}
		if (kitu = 1) then
			begin
				x[1] := 0.5;
				x[2] := 0.0;
				x[3] := 0.0;
			end;
		if (kitu = 2) then
			begin
				x[1] := 5;
				x[2] := 5;
				x[3] := 5;
			end;
		if (kitu = 3) then
			begin
				x[1] := 0.1;
				x[2] := 0.1;
				x[3] := 0.1;
				if (eqnum = 4) then
					x[4] := 0.1;
			end;
		if (kitu = 4) then
			begin
				x[1] := 0.1;
				x[2] := 0.1;
				x[3] := 0.1;
			end;
		if (kitu = 5) then
			begin
				x[1] := 0.0;
				x[2] := 0.86;
				x[3] := 0.02;
			end;
{Initial conditions fo rthe orthonormal frame;}
		for ii := (nn + 1) to nnn do
			begin
				x[ii] := 0.0;
			end;
		for ii := 1 to nn do
			begin
				x[(nn + 1) * ii] := 1.0;
				cum[ii] := 0.0;
			end;
	end;

	procedure initnums;
	begin
		mold := 0.0;
		partychoice := 0;
		ticker3 := 0;
		time := 0.0;
		timeone := 0.0;
		count := 0;
		volume := 0.0;
		Getdatetime(randSeed);
		varsum := 0.0;
		fourier := 1;
		totpower := 0.0;
		totsignal := 0;
	end;

	procedure avevar (data2: RealArrayVV3; uu: integer; var ave, svar: real);
		var
			L: integer;
			s: real;
			weightsum: real;
	begin
		ave := 0.0;
		svar := 0.0;
		weightsum := 0.0;
		for L := 1 to uu do
			begin
				weightsum := weightsum + data2[L, 2];
				ave := ave + (data2[L, 1] * data2[L, 2]);
{writeln('ave=', ave, '  weightsum=', weightsum);}
			end;
		ave := ave / weightsum;
		for L := 1 to uu do
			begin
				s := data2[L, 1] - ave;
				svar := svar + (s * s * data2[L, 2] / weightsum);
			end;
	end;

	procedure four1 (var data: RealArrayVV2; vv, isign: integer);
		var
			ww, qq, v, mmax, m, q, istep, w: integer;
			wtemp, wr, wpr, wpi, wi, theta: double;
			tempr, tempi, wrs, wis: real;
	begin
		v := 2 * vv;
		q := 1;
		for ww := 1 to vv do
			begin
				w := 2 * ww - 1;
				if (q > w) then
					begin
						tempr := data[q];
						tempi := data[q + 1];
						data[q] := data[w];
						data[q + 1] := data[w + 1];
						data[w] := tempr;
						data[w + 1] := tempi;
					end;
				m := v div 2;
				while (m >= 2) and (q > m) do
					begin
						q := q - m;
						m := m div 2;
					end;
				q := q + m;
			end;
		mmax := 2;
		while (v > mmax) do
			begin
				istep := 2 * mmax;
				theta := 6.28318530717959 / (isign * mmax);
				wpr := -2.0 * sqr(sin(0.5 * theta));
				wpi := sin(theta);
				wr := 1.0;
				wi := 0.0;
				for ww := 1 to mmax div 2 do
					begin
						m := 2 * ww - 1;
						wrs := wr;
						wis := wi;
						for qq := 0 to (v - m) div istep do
							begin
								w := m + qq * istep;
								q := w + mmax;
								tempr := wrs * data[q] - wis * data[q + 1];
								tempi := wrs * data[q + 1] + wis * data[q];
								data[q] := data[w] - tempr;
								data[q + 1] := data[w + 1] - tempi;
								data[w] := data[w] + tempr;
								data[w + 1] := data[w + 1] + tempi;
							end;
						wtemp := wr;
						wr := wr * wpr - wi * wpi + wr;
						wi := wi * wpr + wtemp * wpi + wi;
					end;
				mmax := istep;
			end;
		for ww := 1 to (vv div 2) do
			begin
				spectral[ww, 2] := sqr(data[((ww * 2) - 1)]) + sqr(data[(ww * 2)]);
				spectral[ww, 1] := (ww / (vv * h));
				if (suppress = 1) then
					spectral[1, 2] := 0.0;
				if (ww < 20) then
					begin
						writeln('freq=', spectral[ww, 1] : 12 : 5, '  periodogram=  ', spectral[ww, 2] : 15 : 3);
{writeln('fourier numbers  1=', data[((ww * 2) - 1)] : 8 : 6, '    fourier numbers 2=', data[(ww * 2)] : 8 : 6);}
					end;
				totpower := totpower + spectral[ww, 2];
			end;
		totpower := totpower / vv;
		writeln('Total spectral power=', totpower : 20 : 6);
		avevar(spectral, (vv div 2), avefourier, varfourier);
		sdevfour := sqrt(varfourier);
		writeln('avefourier=', avefourier : 10 : 3, '  sdevfour=', sdevfour : 10 : 3);
	end;

	procedure graph;
		var
			k, j, width, height, transx1, transx2, switch, kwisha, fourmax, transform: integer;
			minx1, minx2, maxx1, maxx2, miny1, miny2, maxy1, maxy2, miny3, maxy3, axedraw1, axedraw2: real;
			axeminx, axemaxx, axeminy, axemaxy, scale: real;
			axedraw: array[1..3, 1..2] of real;
	begin
		width := gright - gleft;
		height := gbottom - gtop;
		if (algorithm = 1) then
			kwisha := multiplier - 1;
		if (algorithm = 2) then
			kwisha := multiplier + 5;
		EraseRect(0, 0, 400, 400);
		if (fourier = 1) then
			begin
				transform := 1;
				if ((multiplier - anza) >= 1024) then
					fourmax := 1024;
				if (((multiplier - anza) >= 512) and ((multiplier - anza) < 1024)) then
					fourmax := 512;
				if (((multiplier - anza) >= 256) and ((multiplier - anza) < 512)) then
					fourmax := 256;
				if (((multiplier - anza) >= 128) and ((multiplier - anza) < 256)) then
					fourmax := 128;
				for k := 1 to fourmax do
					begin
						fourdata[((k * 2) - 1)] := y^[((anza - 1) + k), 1];
						fourdata[(k * 2)] := 0.0;
					end;
				four1(fourdata, fourmax, transform);
			end;
		for k := anza to (multiplier - 1) do
			begin
				if (k = anza) then
					begin
						miny1 := y^[anza, 1];
						maxy1 := y^[anza, 1];
						miny2 := y^[anza, 2];
						maxy2 := y^[anza, 2];
						miny3 := y^[anza, 3];
						maxy3 := y^[anza, 3];
					end;
				if (k > anza) then
					begin
						if (y^[k, 1] < miny1) then
							miny1 := y^[k, 1];
						if (y^[k, 1] > maxy1) then
							maxy1 := y^[k, 1];
						if (y^[k, 2] < miny2) then
							miny2 := y^[k, 2];
						if (y^[k, 2] > maxy2) then
							maxy2 := y^[k, 2];
						if (y^[k, 3] < miny3) then
							miny3 := y^[k, 3];
						if (y^[k, 3] > maxy3) then
							maxy3 := y^[k, 3];
					end;
			end;
		for k := multiplier to (multiplier + 5) do
			begin
				if (k = multiplier) then
					begin
						y^[k, 1] := miny1;
						y^[k, 2] := miny2;
						y^[k, 3] := miny3;
					end;
				if (k = multiplier + 1) then
					begin
						y^[k, 1] := maxy1;
						y^[k, 2] := miny2;
						y^[k, 3] := miny3;
					end;
				if (k = multiplier + 2) then
					begin
						y^[k, 1] := miny1;
						y^[k, 2] := miny2;
						y^[k, 3] := miny3;
					end;
				if (k = multiplier + 3) then
					begin
						y^[k, 1] := miny1;
						y^[k, 2] := maxy2;
						y^[k, 3] := miny3;
					end;
				if (k = multiplier + 4) then
					begin
						y^[k, 1] := miny1;
						y^[k, 2] := miny2;
						y^[k, 3] := miny3;
					end;
				if (k = multiplier + 5) then
					begin
						y^[k, 1] := miny1;
						y^[k, 2] := miny2;
						y^[k, 3] := maxy3;
					end;
			end;
		scale := (maxy1 - miny1) / (maxy3 - miny3);
		scale := 1;
		for k := anza to kwisha do
			begin
				for j := 1 to 2 do
					begin
						if (graphtype = 2) then
							begin
								if (algorithm = 1) then
									begin
										newx^[k, j] := (y^[k, j]) * y^[k, 3] * perspect;
									end;
								if (algorithm = 2) then
									begin
										if (j = 1) then
											newx^[k, j] := (y^[k, j]) - (perspect * scale * y^[k, 3]);
										if (j = 2) then
											newx^[k, j] := (y^[k, j]) + (perspect * y^[k, 3]);
									end;
							end;
						if (graphtype = 1) then
							begin
								newx^[k, j] := y^[k, j];
{write(newx^[k, j] : 8 : 6);}
							end;
						if ((k = anza) and (j = 2)) then
							begin
								minx1 := newx^[anza, 1];
								maxx1 := newx^[anza, 1];
								minx2 := newx^[anza, 2];
								maxx2 := newx^[anza, 2];
							end;
					end;
				if (k > anza) then
					begin
						if (newx^[k, 1] < minx1) then
							minx1 := newx^[k, 1];
						if (newx^[k, 1] > maxx1) then
							maxx1 := newx^[k, 1];
						if (newx^[k, 2] < minx2) then
							minx2 := newx^[k, 2];
						if (newx^[k, 2] > maxx2) then
							maxx2 := newx^[k, 2];
					end;
			end;
		moveto(30, 20);
		WriteDraw('x1=', miny1 : 5 : 3, '    x2=', miny2 : 5 : 3, '    x3=', miny3 : 5 : 3);
		moveto(30, 40);
		WriteDraw('x1=', maxy1 : 5 : 3, '    x2=', maxy2 : 5 : 3, '    x3=', maxy3 : 5 : 3);
		TextSize(12);
		PenSize(2, 2);
		for k := anza to kwisha do
			begin
				for j := 1 to 2 do
					begin
						transx1 := round(height - (25 + ((newx^[k, 1] - minx1) / (maxx1 - minx1)) * (height - 50)));
						transx2 := round(25 + ((newx^[k, 2] - minx2) / (maxx2 - minx2)) * (width - 50));
						if ((k = anza) and (j = 1)) then
							moveto(transx2, transx1);
						if ((k = multiplier) and (j = 1)) then
							moveto(transx2, transx1);
						lineto(transx2, transx1);
{writeln(transx2, transx1);}
					end;
				if (traject = 1) then
					write('x1=', y^[k, 1] : 8 : 6, '  x2=', y^[k, 2] : 8 : 6, '  x3=', y^[k, 3] : 8 : 6, ' count=', y^[k, 4]);
				if ((eqnum = 3) and (traject = 1)) then
					writeln;
				if ((eqnum = 4) and (traject = 1)) then
					writeln('  x4=', y^[k, 4] : 8 : 6);
			end;
		if (lyap = 2) then
			writeln('p1=', p[1], '  p2=', p[2]);
		writeln(' Lowx1=', miny1 : 8 : 6, '  Highx1=', maxy1 : 8 : 6, ' Lowx2=', miny2 : 8 : 6, '  Highx2=', maxy2 : 8 : 6, '  Lowx3=', miny3 : 8 : 6, ' Highx3=', maxy3 : 8 : 6);
	end;

	procedure sumchange;
		var
			creature: real;
			k: integer;
	begin
		if (ticker3 > 1) then
			begin
				change := abs(x[item] - varold);
				varsum := varsum + change;
{writeln('var for sum=', x[item] : 8 : 6, '   varsum=', varsum : 8 : 6);}
			end;
		varold := x[item];
		for k := 1 to nn do
			begin
				mx[k] := x[k];
			end;
		if (nn = 3) then
			creature := eq1 * eq2 * eq3;
		if (nn = 4) then
			creature := eq1 * eq2 * eq3 * eq4;
		if ((creature <= 0) and (ticker3 = 1)) then
			signal1 := 1;
		if ((creature > 0) and (ticker3 = 1)) then
			signal1 := 2;
		if ((creature < 0) and (signal1 = 2)) then
			begin
				totsignal := totsignal + 1;
				signal1 := 1;
			end;
		if ((creature > 0) and (signal1 = 1)) then
			begin
				totsignal := totsignal + 1;
				signal1 := 2;
			end;
{writeln('creature=', creature, '   totsignal=', totsignal, '   signal1=', signal1, '   ticker3=', ticker3);}
	end;

	procedure engine;
		var
			k, j, l: integer;
	begin
		initcond;
		initnums;
		while (time <= (iternum + toss)) do
			begin
{Initializing the integrator;}
				timeone := timeone + h;
				ticker3 := ticker3 + 1;
				if (kitu = 3) then
					special3;
				if (kitu = 4) then
					special4;
				if (kitu = 5) then
					special5;
				model;
				multiplier := round((1 / h) * toss);
				if (multiplier > 10000) then
					multiplier := 10000;
				count := count + 1;
				sumchange;
				if (count < multiplier) then
					begin
						for j := 1 to nn do
							begin
								if (graphtype = 2) then
									y^[count, j] := x[j];
								if (graphtype = 1) then
									begin
										y^[count, 1] := x[1];
										y^[count, 2] := time;
										y^[count, 3] := choice;
									end;
							end;
						y^[count, 4] := count;
					end;
				for k := 1 to nnn do
					begin
						x[k] := xnew[k];
					end;
				if (count = multiplier) then
					graph;
{Construct a new orthonormal basis by Gram-Scmidt method;}
{Normalize first vector}
				znorm[1] := 0.0;
				for j := 1 to nn do
					begin
						znorm[1] := znorm[1] + sqr(x[nn * j + 1]);
					end;
				znorm[1] := sqrt(znorm[1]);
				for j := 1 to nn do
					begin
						x[nn * j + 1] := x[nn * j + 1] / znorm[1];
					end;
{Generate the new orthonormal set of vectors;}
				for j := 2 to nn do
					begin
{Generate j-1 coefficients}
						for k := 1 to (j - 1) do
							begin
								gsc[k] := 0.0;
								for l := 1 to nn do
									begin
										gsc[k] := gsc[k] + (x[nn * l + j] * x[nn * l + k]);
									end;
							end;
{Construct a new vector;}
						for k := 1 to nn do
							begin
								for l := 1 to j - 1 do
									begin
										x[nn * k + j] := x[nn * k + j] - (gsc[l] * x[nn * k + l]);
									end;
							end;
{Calculate the vector's norm;}
						znorm[j] := 0;
						for k := 1 to nn do
							begin
								znorm[j] := znorm[j] + sqr(x[nn * k + j]);
							end;
						znorm[j] := sqrt(znorm[j]);
{Normalize the new vector;}
						for k := 1 to nn do
							begin
								x[nn * k + j] := x[nn * k + j] / znorm[j];
							end;
					end;
{Updating running vector magnitudes;}
				if (time > toss) then
					begin
						for k := 1 to nn do
							begin
								cum[k] := cum[k] + log2(znorm[k]);
							end;
{Normalize exponent and print exponents;}
						if ((time > (iternum + toss - h)) and (time <= (iternum + toss))) then
							begin
								write('time=', time : 8 : 3, '  it #=', ticker2);
								for k := 1 to nn do
									begin
										cumkt := cum[k] / time;
										volume := volume + cumkt;
										if (graphonly = 2) then
											write('    ', cumkt : minwidth);
										if (k = 1) then
											begin
												 totderiv1 := cumkt;
{totderiv1 := sdevfour;}
												 totderiv2 := (varsum * totsignal) / time;
												 totderiv3 := (sdevfour * totpower);
											end;
										if (k > 1) then
											begin
												if ((cumkt > totderiv1) and (cumkt < -0.05)) then
												 	totderiv1 := cumkt;
											end;
									end;
								if (volume > 0) then
									alert := alert + 1;
								if (graphonly = 2) then
									writeln('    vol=', volume);
								writeln;
								writeln('turbulence 2=', totderiv2 : minwidth, '   varsum=', varsum, '  totsignal=', totsignal);
								writeln('turbulence 3=', totderiv3 : minwidth);
							end;
					end;
				time := time + h;
			end;
	end;

	procedure answer;
		var
			parnum, anza1, ii: integer;
	begin
		write('What model do you want?   1-pendulum 2-Lorenz 3-Ecosphere 4-Maharishi  5-Nazi  ');
		readln(input, kitu);
{set up number of equations}
		nn := 3;
		nnn := 12;
		eqnum := 3;
		if (kitu = 3) then
			begin
				write('Do you want the four or three equation model? 3-three 4-four  ');
				readln(input, eqnum);
			end;
		if (eqnum = 4) then
			begin
				nn := 4;
				nnn := 20;
				politics := 1;
			end;
		write('Do you want one analysis or an image?  1-One analysis 2-image  ');
		readln(input, lyap);
		write('For what variable do you want to sum the change for turbulence?  ');
		readln(input, item);
		write('Do you want to suppress the first term of the power series?  1-yes  2-no  ');
		readln(input, suppress);
		write('Do you want to view the variable values?  1-yes 2-no  ');
		readln(input, traject);
		write('Do you want Lyaponov exponents?  1-graph only  2-Lyaponov exponents   ');
		readln(input, graphonly);
		if (graphonly = 1) then
			nnn := 3;
		if ((graphonly = 1) and (kitu = 3) and (eqnum = 4)) then
			nnn := 4;
		write('What type of graph do you want?  1-Overtime  2-Phase  ');
		readln(input, graphtype);
		perspect := 1;
		if (graphtype = 2) then
			begin
				write('What type of algorithm do you want to draw the perspective?  1-multiplicative 2-additive  ');
				readln(input, algorithm);
				write('Do you want a 2D or a 3D phase diagram?  0=2D 1=3D  ');
				readln(input, perspect);
			end;
		if (kitu = 3) then
			begin
				write('How many years of lag do you want for the public pollution reaction?   ');
				readln(input, yearslag);
				write('How many years half-life do you want for the current pollution dissipation?   ');
				readln(input, halflife);
				decay := -ln(0.5) / halflife;
				writeln('The decay constant in the model is:  ', decay);
				write('Do you want shocks?  1-yes 2-no  ');
				readln(input, shocks);
				shock := 0;
				if (shocks = 1) then
					begin
						write('What size shock?  ');
						readln(input, shock);
					end;
				if (eqnum = 3) then
					begin
						write('Do you want political policy differences?  1-yes 2-no  ');
						readln(input, politics);
					end;
			end;
		if (kitu = 4) then
			begin
				write('How many time periods of lag do you want for the conflict response?   ');
				readln(input, yearslag);
			end;
		write('How many time periods do you want for the graph?  ');
		readln(input, toss);
		write('Do you want to lag the start of the graph?  1-yes  2-no   ');
		readln(input, sasa);
		if (sasa = 1) then
			begin
				write('At what time period do you want the graph to begin?   ');
				readln(input, anza1);
			end;
		iternum := 1;
		if (graphonly = 2) then
			begin
				write('How many more time periods do you want for Lyaponov analysis?  ');
				readln(input, iternum);
			end;
		write('Step size?  ');
		readln(input, h);
		if (sasa = 1) then
			anza := round(anza1 / h);
		if (sasa = 2) then
			anza := 1;
		if (kitu <> 5) then
			begin
				if (lyap = 1) then
					begin
						write('How many parameters?  ');
						readln(input, parnum);
						for ii := 1 to parnum do
							begin
								write('What is the value for parameter ', ii, '   ?');
								readln(input, p[ii]);
							end;
					end;
			end;
		if (kitu = 5) then
			parms;
		if (lyap = 2) then
			begin
				writeln('You can vary your first two parameters between ranges.');
				write('What is your lower value in the range for your first parameter?    ');
				readln(input, lowp1);
				write('What is your higher value in the range for your first parameter?    ');
				readln(input, highp1);
				write('What is your lower value in the range for your second parameter?    ');
				readln(input, lowp2);
				write('What is your higher value in the range for your second parameter?    ');
				readln(input, highp2);
				write('How many other parameters are there in the model?  ');
				readln(input, parnum);
				for ii := 3 to (parnum + 2) do
					begin
						write('What is the value for parameter ', ii, '   ?');
						readln(input, p[ii]);
					end;
				write('How large do you want the image matrix?  ');
				readln(input, dim);
			end;
		lagg := round(yearslag / h);
		if (yearslag = 0) then
			lagg := 1;
	end;

	procedure minmax;
	begin
		ticker := ticker + 1;
		if (ticker = 1) then
			begin
				min := totderiv1;
				max := totderiv1;
				min2 := totderiv2;
				max2 := totderiv2;
				min3 := totderiv3;
				max3 := totderiv3;
			end;
		if totderiv1 < min then
			min := totderiv1;
		if totderiv1 > max then
			max := totderiv1;
		if totderiv2 < min2 then
			min2 := totderiv2;
		if totderiv2 > max2 then
			max2 := totderiv2;
		if totderiv3 < min3 then
			min3 := totderiv3;
		if totderiv3 > max3 then
			max3 := totderiv3;
	end;

	procedure axes;
		var
			d1, d2: integer;
	begin
		for d2 := 1 to 2 do
			begin
				for d1 := 1 to dim do
					begin
						write(ff1, ' ', d1);
						write(ff2, ' ', d1);
						write(ff3, ' ', d1);
					end;
				writeln(ff1);
				writeln(ff2);
				writeln(ff3);
			end;
	end;

	procedure matfill;
		var
			c1, c2: integer;
	begin
		ticker := 0;
		for c2 := dim downto 1 do
			begin
				shift := lowp2 + ((highp2 - lowp2) * ((c2 - 1) / (dim - 1)));
				p[2] := shift;
				for c1 := 1 to dim do
					begin
						ticker2 := ticker2 + 1;
						p[1] := lowp1 + ((highp1 - lowp1) * ((c1 - 1) / (dim - 1)));
						engine;
						minmax;
						image1^[(dim - c2 + 1), c1] := totderiv1;
						image2^[(dim - c2 + 1), c1] := totderiv2;
						image3^[(dim - c2 + 1), c1] := totderiv3;
					end;
			end;
	end;

	procedure printit;
		var
			c1, c2: integer;
	begin
		writeln(ff1, dim, ' ', dim);
		writeln(ff2, dim, ' ', dim);
		writeln(ff3, dim, ' ', dim);
		writeln('dimension=', dim, ' by ', dim);
		writeln(ff1, max : minwidth, ' ', min : minwidth);
		writeln(ff2, max2 : minwidth, ' ', min2 : minwidth);
		writeln(ff3, max3 : minwidth, ' ', min3 : minwidth);
		writeln('max=', max : minwidth, ' min=', min : minwidth);
		writeln('max=', max2 : minwidth, ' min=', min2 : minwidth);
		writeln('max=', max3 : minwidth, ' min=', min3 : minwidth);
		writeln;
		axes;
		for c2 := dim downto 1 do
			begin
				for c1 := 1 to dim do
					begin
						write(ff1, ' ', image1^[(dim - c2 + 1), c1] : minwidth);
						write(ff2, ' ', image2^[(dim - c2 + 1), c1] : minwidth);
						write(ff3, ' ', image3^[(dim - c2 + 1), c1] : minwidth);
						if c1 < 6 then
							write(' ', image1^[(dim - c2 + 1), c1] : minwidth);
					end;
				writeln(ff1);
				writeln(ff2);
				writeln(ff3);
				writeln;
			end;
		writeln;
		for c2 := dim downto 1 do
			begin
				for c1 := 1 to dim do
					begin
						if c1 < 6 then
							write(' ', image2^[(dim - c2 + 1), c1] : minwidth);
					end;
				writeln;
			end;
		writeln;
		for c2 := dim downto 1 do
			begin
				for c1 := 1 to dim do
					begin
						if c1 < 6 then
							write(' ', image3^[(dim - c2 + 1), c1] : minwidth);
					end;
				writeln;
			end;
	end;

begin
	initialize;
	showtext;
	showdrawing;
	answer;
	if (lyap = 1) then
		engine;
	if (lyap = 2) then
		begin
			rewrite(ff1, f1name);
			rewrite(ff2, f2name);
			rewrite(ff3, f3name);
			new(image1);
			new(image2);
			new(image3);
			matfill;
			printit;
			close(ff1);
			close(ff2);
			close(ff3);
			writeln('lpwp1=', lowp1, '  highp1=', highp1, '  lowp2=', lowp2, '  highp2=', highp2, '  p3=', p[3]);
		end;
	writeln('A nondisipative system was detected', alert, '   times.');
end.

 

This entire site is Copyright © 1997-2024 by Courtney Brown. All Rights Reserved.
DISCLAIMER
URL: https://courtneybrown.com