The following program was used to explore the estimated catastrophe models that are
analyzed in chapters 3 and 5 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 catastrophe (input, output, Image1, Image2, Image3);
const
minwidth = 15;
maxsize = 1000;
maxbig = 32000;
f1name = 'ImageLyap';
f2name = 'ImageSum';
f3name = 'ImagePower';
ps1name = 'PScatcode';
matnn = 5;
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..maxbig, 1..matnn] of real;
HugePtr3 = ^BigArray3;
BigArray3 = array[1..10000] of real;
GraphPtr = ^GraphBig;
GraphBig = array[1..maxbig, 1..matnn] 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, lagg, again, limits: integer;
x1low, x1high, x2low, x2high, perspect, shift, mold, w1, acc, x1dlow, x1dhigh, x2dlow, x2dhigh: real;
lx1, ux1, lx2, ux2, xprim, lowp1, lowp2, highp1, highp2, rDR, timeone, expand, pi: real;
min2, max2, min3, max3, context: real;
ff1, ff2, ff3, ps1: text;
row, col, ticker, ticker2, ticker3, graphtype, traject, sasa, year, method, highlite: integer;
nn, nnn, eqnum, laggedid, count, item, anza, graphonly, series, cross, psprint: integer;
suppress, connect, condques, rotation, angle: integer;
filein: text;
y: HugePtr2;
yy: HugePtr2;
p: array[1..20] of real;
pp: 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(yy);
new(newx);
ticker2 := 0;
laggedid := 0;
pi := 3.141549;
end;
procedure condition;
begin
p[2] := p[2] + (pp[2] * w1);
p[4] := p[4] + (pp[4] * w1);
p[5] := p[5] + (pp[5] * w1);
p[6] := p[6] + (pp[6] * w1);
p[7] := p[7] + (pp[7] * w1);
p[8] := p[8] + (pp[8] * w1);
p[9] := p[9] + (pp[9] * w1);
p[19] := p[19] + (pp[19] * w1);
end;
procedure parm1;
begin
p[2] := 0;
p[4] := 0.66407;
p[5] := -0.42078;
p[6] := 0.34843;
p[7] := 0.71071;
p[8] := 0.12868;
p[9] := 0.19019;
p[19] := 0.14762;
pp[4] := 0.81578;
pp[5] := 0.09054;
pp[6] := 0.10649;
pp[7] := 0.99656;
pp[8] := 0.20632;
pp[9] := 0.48091;
pp[19] := -1.02490;
lowp1 := 0 - expand;
highp1 := 0.42 + expand;
x1low := -10.0;
x1high := 15.0;
x1dlow := 0.5;
x1dhigh := 1.0;
x2dlow := 0.0;
x2dhigh := 0.78;
x2low := 0.0 - expand;
x2high := 0.78 + expand;
condition;
end;
procedure parm2;
begin
p[2] := 0;
p[4] := 1.20418;
p[5] := -0.51405;
p[6] := 0.15501;
p[7] := 1.13798;
p[8] := 0.20683;
p[9] := 0.19554;
p[19] := -0.50866;
pp[4] := 0.12824;
pp[5] := -0.02477;
pp[6] := -0.00479;
pp[7] := -0.06923;
pp[8] := -0.49641;
pp[9] := 0.49962;
pp[19] := 0.10608;
lowp1 := 0 - expand;
highp1 := 0.7566 + expand;
x1low := -10.0;
x1high := 15.0;
x1dlow := 0.4;
x1dhigh := 1.0;
x2dlow := 0.0;
x2dhigh := 0.65;
x2low := 0.0 - expand;
x2high := 0.65 + expand;
condition;
end;
procedure parm3;
begin
p[2] := 1;
p[5] := 5.98186;
p[6] := 0.83441;
p[7] := 1.03085;
p[8] := 0.88157;
p[9] := 0.41774;
p[10] := -0.13497;
p[11] := 1.24597;
if (condques = 1) then
context := 0;
{Republican standardized}
if (condques = 2) then
begin
context := 2;
pp[5] := -0.00334;
pp[6] := -0.04221;
pp[7] := 0.00608;
pp[8] := -0.09321;
pp[9] := -0.06038;
pp[10] := -0.11694;
pp[11] := -0.01458;
end;
{Democratic standardized}
if (condques = 3) then
begin
context := 2.7;
pp[5] := 0.001897;
pp[6] := 0.008622;
pp[7] := -0.006862;
pp[8] := 0.052890;
pp[9] := 0.066099;
pp[10] := 0.078435;
pp[11] := 0.043908;
end;
{Neighbor -3R, 0M, and 3D}
if (condques = 4) then
begin
context := 3;
pp[5] := 0.20720;
pp[6] := 0.07340;
pp[7] := -0.26706;
pp[8] := 0.03561;
pp[9] := -0.04350;
pp[10] := 0.01784;
pp[11] := 0.04216;
end;
{Rae's F nonstandardized}
if (condques = 5) then
begin
context := 3;
pp[5] := -0.002982;
pp[6] := -0.000475;
pp[7] := 0.007336;
pp[8] := 0.015913;
pp[9] := 0.036962;
pp[10] := -0.018868;
pp[11] := 0.022124;
{standardized}
context := 2.4;
context := -3.7;
pp[5] := -0.00245;
pp[6] := -0.11443;
pp[7] := -0.03304;
pp[8] := -0.10406;
pp[9] := -0.02834;
pp[10] := -0.16631;
pp[11] := -0.03693;
end;
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: real; exponent: integer): real;
var
tempvar: real;
t: integer;
begin
tempvar := base;
if (exponent > 1) then
begin
for t := 1 to (exponent - 1) do
begin
tempvar := tempvar * base;
end;
end;
if (exponent < 0) then
begin
tempvar := 1 / base;
for t := 1 to (abs(exponent) - 1) do
begin
tempvar := tempvar * (1 / base);
end;
end;
if (exponent = 0) then
tempvar := 1;
power := tempvar;
end;
function powernoint (base: real; exponent: real): real;
begin
powernoint := exp(ln(base) * exponent);
end;
function eq1: real;
var
w1, w2: real;
begin
if (kitu = 1) then
eq1 := (p[2] * (power(mx[1], 3))) - (p[1] * mx[1]) + (p[3] * mx[2]);
if (kitu = 2) then
eq1 := (p[2] * (power((mx[1] - 1), 3))) + (p[1] * (mx[1] - 1)) - (p[3] * (mx[2] - 1));
if (kitu = 3) then
eq1 := ((p[2] * (power((mx[1] - 1), 3))) + (p[1] * (mx[1] - 1)) - (p[3] * (mx[2] - 1))) * (1 - mx[1]) * mx[1];
if (kitu = 4) then
eq1 := (p[19] * mx[2]) - (p[4] * mx[2] * p[1]) + (p[5] * (p[6] * (mx[1] - p[7]))) + (p[8] * (power((p[6] * (mx[1] - p[7])), 2))) + (p[9] * (power((p[6] * (mx[1] - p[7])), 3)));
{eq1 := (p[19] * mx[2]) - (p[4] * mx[2] * p[1]) + (p[5] * (p[6] * (mx[1] - p[7]))) + (p[8] * (power((p[6] * (mx[1] - p[7])), 2))) + (p[9] * (power((p[6] * (mx[1] - p[7])), 3)));}
if (kitu = 5) then
begin
{set p5=0 to get the plane, 1 otherwise.}
w1 := (1 + mx[2] - p[1]) / 2;
w2 := (mx[2] + p[1]) / 2;
eq1 := (p[5] * (p[2] * w2 * (1 - w2) * (mx[1] - p[4]) - (p[3] * power(mx[1] - p[4], 3))) + (p[6] * w1)) - mx[1];
{eq1 := (p[5] * (p[2] * (mx[1] - p[4]) - (p[3] * power(mx[1] - p[4], 3))) + w1)-mx[1];}
{eq1 := (p[5] * (p[2] * w2 * (mx[1] - p[4]) - (p[3] * power(mx[1] - p[4], 3))) + w1) - mx[1];}
{eq1 := (p[5] * (p[2] * (1 - w2) * (mx[1] - p[4]) - (p[3] * power(mx[1] - p[4], 3))) + w1)-mx[1];}
end;
if (kitu = 6) then
begin
w1 := (1 + mx[2] - p[1]) / 2;
w2 := (mx[2] + p[1]) / 2;
eq1a := (p[9] + (pp[9] * context)) * ((((p[5] + pp[5] * context) * ((p[10] + (pp[10] * context)) + w2) * (1 - (p[11] + (pp[11] * context)) * w2) * (mx[1] - (p[6] + (pp[6] * context))));
eq1 := eq1a - ((p[7] + pp[7] * context) * ((power(mx[1] - (p[6] + (pp[6] * context)), 3)))) + ((p[8] + (pp[8] * context)) * w1)) - mx[1]);
end;
end;
function eq2: real;
begin
end;
function eq3: real;
begin
end;
function eq4: real;
begin
end;
function fx (xx: real): real;
begin
mx[1] := xx;
fx := eq1;
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 initnums;
begin
count := 0;
ticker := 0;
end;
procedure post1;
begin
writeln(ps1, '150 250 translate');
writeln(ps1, '1 setlinewidth');
writeln(ps1, '/Times-Bold findfont 15 scalefont setfont');
writeln(ps1, '1 setlinejoin');
end;
procedure post2;
begin
writeln(ps1, 'stroke showpage');
end;
procedure graph;
var
k, j, width, height, transx1, transx2, switch, kwisha, fourmax, transform, diam: integer;
minx1, minx2, maxx1, maxx2, miny1, miny2, maxy1, maxy2, miny3, maxy3, axedraw1, axedraw2: real;
axeminx, axemaxx, axeminy, axemaxy, scale, labset, transx3, transx4, phi, scale1, scale2, scale3: real;
axedraw: array[1..3, 1..2] of real;
begin
width := gright - gleft;
height := gbottom - gtop;
labset := 0.05;
anza := 1;
multiplier := count + 1;
kwisha := multiplier + 7 + 6 + 3 + 1;
EraseRect(0, 0, 400, 400);
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;
{Use the next two lines to set the maximum and minimum values for the vertical axis.}
{maxy1 := 1;}
{miny1 := 0;}
for k := multiplier to (multiplier + 7) 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;
if (k = multiplier + 6) then
begin
y^[k, 1] := miny1;
y^[k, 2] := maxy2;
y^[k, 3] := maxy3;
end;
if (k = multiplier + 7) then
begin
y^[k, 1] := miny1;
y^[k, 2] := maxy2;
y^[k, 3] := miny3;
end;
end;
for k := (multiplier + 8) to (multiplier + 17) do
begin
{Y axis value labels}
if (k = multiplier + 8) then
begin
y^[k, 1] := miny1 + ((maxy1 - miny1) * 0.5 * labset);
y^[k, 2] := miny2 - ((maxy2 - miny2) * 6 * labset);
y^[k, 3] := miny3;
end;
if (k = multiplier + 9) then
begin
y^[k, 1] := maxy1 - ((maxy1 - miny1) * labset);
y^[k, 2] := miny2 - ((maxy2 - miny2) * 5 * labset);
y^[k, 3] := miny3;
end;
{X axis value labels}
if (k = multiplier + 10) then
begin
y^[k, 1] := miny1 + ((maxy1 - miny1) * 0.5 * labset);
y^[k, 2] := miny2 + ((maxy2 - miny2) * labset);
y^[k, 3] := miny3;
end;
if (k = multiplier + 11) then
begin
y^[k, 1] := miny1 - ((maxy1 - miny1) * 1.5 * labset);
y^[k, 2] := maxy2 - ((maxy2 - miny2) * 2 * labset);
y^[k, 3] := miny3;
end;
{Z axis value labels}
if (k = multiplier + 12) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2 - ((maxy2 - miny2) * 7 * labset);
y^[k, 3] := miny3 + ((maxy3 - miny3) * 3 * labset);
end;
if (k = multiplier + 13) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2 - ((maxy2 - miny2) * 6.5 * labset);
y^[k, 3] := maxy3;
end;
{The X-axis label}
if (k = multiplier + 14) then
begin
y^[k, 1] := miny1 - ((maxy1 - miny1) * 3 * labset);
y^[k, 2] := miny2 + ((maxy2 - miny2) * 0.5);
y^[k, 3] := miny3;
end;
{The Z-axis label}
if (k = multiplier + 15) then
begin
y^[k, 1] := miny1;
y^[k, 2] := miny2 - ((maxy2 - miny2) * 0.6);
y^[k, 3] := miny3 + ((maxy3 - miny3) * 0.7);
end;
{The Y-axis label}
if (k = multiplier + 16) then
begin
y^[k, 1] := miny1 + ((maxy1 - miny1) * 4 * labset);
y^[k, 2] := miny2 - ((maxy2 - miny2) * 3 * labset);
y^[k, 3] := miny3;
end;
{The TITLE}
if (k = multiplier + 17) then
begin
y^[k, 1] := maxy1 + 6 * labset;
y^[k, 2] := miny2 - ((maxy2 - miny2) * 5 * labset);
y^[k, 3] := miny3;
end;
end;
scale1 := (maxy1 - miny1);
scale2 := (maxy2 - miny2);
scale3 := (maxy3 - miny3);
phi := angle * (pi / 180);
for k := anza to kwisha do
begin
yy^[k, 1] := y^[k, 1] / scale1;
yy^[k, 2] := y^[k, 2] / scale2;
yy^[k, 3] := y^[k, 3] / scale3;
for j := 1 to 2 do
begin
if (j = 1) then
newx^[k, j] := yy^[k, j] + (yy^[k, 3] * perspect * cos(phi));
if (j = 2) then
newx^[k, j] := yy^[k, j] + (yy^[k, 3] * perspect * sin(phi));
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) and (k < kwisha)) 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(70, 20);}
{WriteDraw('xmin=', miny2 : 5 : 3, ' ymin=', miny1 : 5 : 3, ' zmin=', miny3 : 5 : 3);}
{moveto(70, 40);}
{WriteDraw('xmax=', maxy2 : 5 : 3, ' ymax=', maxy1 : 5 : 3, ' zmax=', maxy3 : 5 : 3);}
TextSize(12);
PenSize(1, 1);
if (psprint = 1) then
post1;
for k := anza to kwisha 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));
transx3 := (25 + ((newx^[k, 1] - minx1) / (maxx1 - minx1)) * (height - 50));
transx4 := (25 + ((newx^[k, 2] - minx2) / (maxx2 - minx2)) * (width - 50));
diam := 1;
if ((y^[k, 1] >= x1dlow) and (y^[k, 1] <= x1dhigh)) then
diam := 2;
if (k < multiplier) then
begin
{PaintCircle(transx2, transx1, diam);}
if (y^[k, 5] = 1) then
begin
moveto(transx2, transx1);
WriteDraw('.');
if (highlite = 1) then
begin
if ((y^[k, 1] >= x1dlow) and (y^[k, 1] <= x1dhigh)) then
begin
TextSize(13);
WriteDraw('.');
TextSize(12);
end;
{PaintCircle(transx2, transx1, diam);}
end;
end;
if (y^[k, 5] = 2) then
begin
lineto(transx2, transx1);
end;
if (psprint = 1) then
begin
if (y^[k, 5] = 1) then
begin
writeln(ps1, transx4 : 7 : 2, ' ', transx3 : 7 : 2, ' ', 'moveto ');
writeln(ps1, '(.) show');
if (highlite = 1) then
begin
if ((y^[k, 1] >= x1dlow) and (y^[k, 1] <= x1dhigh)) then
begin
writeln(ps1, '/Helvetica-Bold findfont 25 scalefont setfont');
writeln(ps1, '(.) show');
writeln(ps1, '/Times-Bold findfont 15 scalefont setfont');
end;
end;
end;
if (y^[k, 5] = 2) then
begin
writeln(ps1, transx4 : 7 : 2, ' ', transx3 : 7 : 2, ' ', 'lineto ');
end;
end;
end;
if (k = multiplier) then
begin
moveto(transx2, transx1);
if (psprint = 1) then
begin
writeln(ps1, transx4 : 7 : 2, ' ', transx3 : 7 : 2, ' ', 'moveto ');
writeln(ps1, '/Helvetica-Bold findfont');
writeln(ps1, '15 scalefont');
writeln(ps1, 'setfont');
end;
end;
if ((k > multiplier) and (k < (multiplier + 8))) then
begin
lineto(transx2, transx1);
if (psprint = 1) then
writeln(ps1, transx4 : 7 : 2, ' ', transx3 : 7 : 2, ' ', 'lineto ');
end;
if (k > (multiplier + 7)) then
begin
moveto(transx2, transx1);
if (k = (multiplier + 8)) then
WriteDraw(miny1 : 5 : 3);
if (k = (multiplier + 9)) then
WriteDraw(maxy1 : 5 : 3);
if (k = (multiplier + 10)) then
WriteDraw(miny2 : 5 : 3);
if (k = (multiplier + 11)) then
WriteDraw(maxy2 : 5 : 3);
if (k = (multiplier + 12)) then
WriteDraw(miny3 : 5 : 3);
if (k = (multiplier + 13)) then
WriteDraw(maxy3 : 5 : 3);
if (psprint = 1) then
begin
writeln(ps1, transx4 : 7 : 2, ' ', transx3 : 7 : 2, ' ', 'moveto ');
if (k = (multiplier + 8)) then
writeln(ps1, '(', miny1 : 5 : 3, ') show');
if (k = (multiplier + 9)) then
writeln(ps1, '(', maxy1 : 5 : 3, ') show');
{if (k = (multiplier + 10)) then writeln(ps1, '(' , miny2 : 5 : 3, ') show');}
{if (k = (multiplier + 11)) then writeln(ps1, '(', maxy2 : 5 : 3, ') show');}
if (k = (multiplier + 12)) then
writeln(ps1, '(', miny3 : 5 : 3, ') show');
if (k = (multiplier + 13)) then
writeln(ps1, '(', maxy3 : 5 : 3, ') show');
if (k = (multiplier + 14)) then
begin
writeln(ps1, (transx4 - 72) : 7 : 2, ' ', (transx3 - 130) : 7 : 2, ' ', 'moveto ');
writeln(ps1, '(X AXIS RANGE=[', miny2 : 5 : 3, ', ', maxy2 : 5 : 3, ']) show ');
end;
if (k = (multiplier + 15)) then
writeln(ps1, '(Z AXIS) show ');
if (k = (multiplier + 16)) then
begin
writeln(ps1, 'gsave');
writeln(ps1, '90 rotate');
writeln(ps1, '(Y AXIS) show');
writeln(ps1, 'grestore');
end;
if (k = (multiplier + 17)) then
writeln(ps1, '(TITLE) show ');
end;
end;
if (traject = 1) then
writeln('x=', y^[k, 2] : 8 : 6, ' y=', y^[k, 1] : 8 : 6, ' z=', y^[k, 3] : 8 : 6, ' count=', y^[k, 4] : 8 : 1);
end;
if (psprint = 1) then
post2;
end;
procedure newton;
var
control1, control2, i, v, control3, newtnum, x2dim: integer;
continc1, continc2, x1next, x2next, test1, test2, x1newt, x2newt, x1save, x2save: real;
begin
x2dim := dim * 3;
newtnum := dim;
acc := 0.005;
continc1 := (x1high - x1low) / newtnum;
continc2 := (x2high - x2low) / (x2dim);
mx[1] := x1low;
mx[2] := x2low;
x1newt := mx[1];
x2save := x2low;
for control2 := 1 to (x2dim + 1) do
begin
x1newt := x1low;
x1save := x1newt;
control3 := 0;
for control1 := 1 to (newtnum + 1) do
begin
mx[1] := x1newt;
for i := 1 to 9 do
begin
x1next := mx[1] - (eq1 / par1x1);
mx[1] := x1next;
end;
test1 := abs(eq1);
test2 := abs(mx[1] - x1save);
{writeln(test1 : 8 : 5);}
if ((test1 < acc) and (test2 > acc)) then
begin
count := count + 1;
y^[count, 1] := mx[1];
y^[count, 2] := mx[2];
y^[count, 3] := p[1];
y^[count, 4] := count;
y^[count, 5] := connect;
if (limits = 1) then
begin
if (mx[1] < x1low) then
y^[count, 1] := x1low;
if (mx[1] > x1high) then
y^[count, 1] := x1high;
end;
connect := 1;
{writeln(y^[count, 1], ' ', y^[count, 2], ' ', y^[count, 3], ' ', y^[count, 4], ' p2=', p[2]);}
end;
x1save := mx[1];
x1newt := x1newt + continc1;
end;
x2save := mx[2];
mx[2] := mx[2] + continc2;
{writeln(continc2, ' ', mx[2]);}
end;
writeln('Count = ', count);
end;
function rtbis (k1, k2, xacc: real): real;
label
99;
const
jmax = 40;
var
dx, f, fmid, xmid, rtb: real;
j: integer;
begin
fmid := fx(k2);
f := fx(k1);
if f * fmid >= 0.0 then
begin
writeln('Pause in rtbis: root must be bracketed for bisection');
readln;
end;
if f < 0.0 then
begin
rtb := k2;
dx := k2 - k1;
end
else
begin
rtb := k2;
dx := k1 - k2;
end;
for j := 1 to jmax do
begin
dx := dx * 0.5;
xmid := rtb + dx;
fmid := fx(xmid);
if fmid <= 0.0 then
rtb := xmid;
if (abs(dx) < xacc) or (fmid = 0.0) then
goto 99;
end;
writeln('Pause in rtbis: too many bisections');
readln;
99:
rtbis := rtb;
end;
procedure bisect;
var
control1, control2, i, v, control3, newtnum, binum: integer;
continc1, continc2, x1next, x2next, test1, test2, x1newt, x2newt, x1save, x2save: real;
begin
binum := 20;
continc1 := (x1high - x1low) / binum;
continc2 := (x2high - x2low) / dim;
mx[1] := x1low;
mx[2] := x2low;
x1newt := mx[1];
for control2 := 1 to (dim + 1) do
begin
x1newt := x1low;
control3 := 0;
for control1 := 1 to binum do
begin
mx[1] := x1newt;
count := count + 1;
y^[count, 1] := mx[1];
y^[count, 2] := mx[2];
y^[count, 3] := p[1];
y^[count, 4] := count;
if (limits = 1) then
begin
if (mx[1] < x1low) then
y^[count, 1] := x1low;
if (mx[1] > x1high) then
y^[count, 1] := x1high;
end;
{writeln(y^[count, 1], ' ', y^[count, 2], ' ', y^[count, 3], ' ', y^[count, 4], ' p2=', p[2]);}
x1newt := x1newt + continc1;
end;
mx[2] := mx[2] + continc2;
{writeln(continc2, ' ', mx[2]);}
end;
end;
procedure answer;
var
parnum, ii: integer;
begin
writeln('What model do you want? 1-fragmentation 2-frag w/shift 3-frag w/shift&lim ');
write(' 4-Nazis 5-catvoter 6-estimated cat voter ');
readln(input, kitu);
if (kitu = 1) then
nn := 3;
write('Do you want one picture or a series? 1-one 2-series ');
readln(input, series);
write('Do you want to print the graph using a postscript file? 1-yes 2-no ');
readln(input, psprint);
write('Do you want to view the variable values? 1-yes 2-no ');
readln(input, traject);
write('Do you want a cross-section? 1-yes 2-no ');
readln(input, cross);
write('Do you want Newtons method or a bisection method? 1-Newtons 2-bisection ');
readln(input, method);
if (cross = 2) then
begin
write('Do you want a 2D or a 3D phase diagram? 0=2D 1=3D ');
readln(input, perspect);
write('What rotation angle do you want? 1=120 degrees 2-other ');
readln(input, rotation);
if (rotation = 1) then
angle := 120;
if (rotation = 2) then
begin
write(' What specific angle do you want ? ');
readln(input, angle);
end;
end;
write('Do you want vertical limits on the plot? 1-yes 2-no ');
readln(input, limits);
write('What resolution do you want on your axes? ');
readln(input, dim);
if (kitu = 4) then
begin
write('What year do you want? 1:1928-30 2:1930-32 ');
readln(input, year);
write('What is the value of the conditioning variable? ');
readln(input, w1);
write('How much do you want to expand the floor axes? Any number from 0 to ?? ');
readln(input, expand);
if (year = 1) then
parm1;
if (year = 2) then
parm2;
write('Do you want to highlight the vertical data areas? 1-yes 2-no ');
readln(input, highlite);
end;
if (series = 2) then
begin
write('For a series of plots, 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);
end;
if (kitu <> 4) then
begin
if (kitu < 4) then
begin
write('What is your low value for variable x1? ');
readln(input, x1low);
write('What is your high value for variable x1? ');
readln(input, x1high);
write('What is your low value for variable x2? ');
readln(input, x2low);
write('What is your high value for variable x2? ');
readln(input, x2high);
if (cross = 2) then
begin
write('What is your lower value in the range for your z axis parameter? ');
readln(input, lowp1);
write('What is your higher value in the range for your z axis parameter? ');
readln(input, highp1);
end;
end;
if (kitu = 5) then
begin
write('How much do you want to expand the floor axes? Any number from 0 to ?? ');
readln(input, expand);
x1low := -1;
x1high := 2;
x2low := 0 - expand;
x2high := 1 + expand;
lowp1 := 0 - expand;
highp1 := 1 + expand;
end;
if (kitu = 6) then
begin
writeln('Hoiw do you want the model conditioned? 1-none 2-Republican 3-Democratic 4-Neighbors ');
write('5-Raes F ');
readln(input, condques);
write('How much do you want to expand the floor axes? Any number from 0 to ?? ');
readln(input, expand);
x1low := -1;
x1high := 2;
x2low := 0 - expand;
x2high := 1 + expand;
lowp1 := 0 - expand;
highp1 := 1 + expand;
parm3;
lowp2 := p[2];
highp2 := p[2];
end;
if (cross = 1) then
begin
write('What is the value of the z control variable? ');
readln(input, lowp1);
highp1 := lowp1;
perspect := 0;
end;
if (kitu <> 6) then
begin
if (series = 1) then
begin
write('What is the value of your second parameter? ');
readln(input, p[2]);
lowp2 := p[2];
highp2 := p[2];
end;
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;
end;
end;
end;
procedure range;
var
c1, c2, dimc2, dimc1: integer;
begin
initnums;
if (series = 1) then
dimc2 := 1;
if (series = 2) then
dimc2 := dim;
if (cross = 1) then
dimc1 := 1;
if (cross = 2) then
dimc1 := dim;
for c2 := dimc2 downto 1 do
begin
count := 0;
shift := lowp2 + ((highp2 - lowp2) * ((c2 - 1) / (dim - 1)));
p[2] := shift;
ticker := ticker + 1;
if (series = 2) then
begin
writeln('p2=', p[2] : 8 : 4, ' ticker', ticker);
end;
for c1 := 1 to dimc1 do
begin
ticker2 := ticker2 + 1;
connect := 1;
p[1] := lowp1 + ((highp1 - lowp1) * ((c1 - 1) / (dim - 1)));
if (method = 1) then
newton;
if (method = 2) then
bisect;
end;
graph;
end;
end;
begin
initialize;
showtext;
showdrawing;
answer;
if (psprint = 1) then
rewrite(ps1, ps1name);
range;
if (psprint = 1) then
close(ps1);
end.