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.