/*====================================================================*/
/*  Poincare-Abbildung x -> k * x * (1 - x)                           */
/*                     (0 < x < 1, 1 <= k < 4)                        */
/*  (c't 3/87)                                                        */
/*  Ausgabe auf den Bildschirm oder auf den Drucker (NEC P6/P7/P2200) */
/*  in der Aufloesung 360*180 P/25.4mm und Doppeldruck.               */
/*  Der x-Startwert, der gewuenschte Ausschnitt (k vertikal, x hori-  */
/*  zontal), die Vorlauflaenge ('Initialisierung in c't) und die An-  */
/*  zahl der Punkte ('Wiederholungen') werden nach dem Programmstart  */
/*  abgefordert.                                                      */
/*  Bei Ausgabe auf dem Bildschirm muss mit SHIFT/HELP umgeschaltet   */
/*  werden.                                                           */
/*  Bei Ausgabe auf den Drucker wird dort eine Skala mit ausgegeben.  */
/*                                                                    */
/*      Wolfgang Windsheimer 1.4.1988   (fuer PEARL-Pool)             */
/*                                                                    */
/*  Modulname: CHAOS                                                  */
/*  Starttask: CHA                                                    */
/*====================================================================*/

S=6000; /*-M*/;/*-L*/;/*-T*/
MODULE CHAOS;

SYSTEM;
       Dialog:A1 <->;
       Druck: PP (TFU=30) ->;

PROBLEM;
       SPC Dialog DATION INOUT ALPHIC CONTROL(ALL);
       SPC Druck  DATION OUT ALPHIC CONTROL(ALL);
       SPC SCRSET ENTRY GLOBAL;
       SPC CLEAR ENTRY GLOBAL;

DCL oben  INV BIT(1) INIT ('0'B1);
DCL unten INV BIT(1) INIT ('1'B1);
DCL r8 INV FIXED(15) INIT (-8);
DCL l8 INV   BIT(16) INIT ('00FF'B4);

/*----- Groessenangabe fuer Druckseite (NEC P2200) ------------------*/

DCL prtlines INV FIXED INIT (72);                 /* Druckkopfzeilen */
DCL pixlines INV FIXED INIT (1728);   /* Pixelzeilen (24 * prtlines) */
DCL toline1 INV FIXED INIT (32);  /* Zeile 1 um ./180 inch ueber der */
                                            /* Papiereinzugsposition */
DCL pixcols INV FIXED INIT (2400);                   /* Pixelspalten */
DCL skalcols INV FIXED INIT (48);      /* Spalten fuer Skalenstriche */
DCL pcl INV FIXED INIT (49); /* skalcols+1  erste Pixelspaltennummer */
DCL pcr INV FIXED INIT (2449); /* +pixcols letzte Pixelspaltennummer */
DCL pixline (3, 2496) BIT(16);          /* 3, pixcols + 2 * skalcols */

/*-------------------------------------------------------------------*/

DCL (xmin,xmax,xrange,ymax,ymin,yrange,xwert,xpix, xanf,ywert) FLOAT;
DCL (sdiff1,sdiff2,sdiff3,snext1,snext2,snext3) FLOAT;
DCL (tiefe,punkte,tiefemax,SPALTE,ZEILE)FIXED;
DCL pchar CHAR;
DCL pb (24) CHAR;
DCL pixlno FIXED;
DCL (xskal, yskal, sline, spunkte, skalmin)FLOAT;
DCL (xrest, yrest)FIXED;
DCL teilung (4)FIXED;


/* ----------------------------------------------------------------- */
/* In FLOAT-Arithmetik groesstes modul-Vielfaches von x bestimmen.   */
/* Dabei wird der unzureichende Wertebereich von ENTIER beruecksich- */
/* tigt.                                                             */
/* ----------------------------------------------------------------- */

fmlow: PROC(x FLOAT, modul FLOAT) RETURNS (FLOAT) REENT;
DCL a FLOAT;

a = 0;
IF modul == 0 THEN RETURN(0); FIN;

IF ABS(x / 32760.0) >= ABS(modul) THEN
  a = fmlow(x, modul * 32760.0);
  x = x - a;
FIN;
RETURN (a + modul * TOFLOAT(ENTIER(x / modul)));

END;

/* ----------------------------------------------------------------- */
/* Zu gegebenen (krummen) Anfangs- und Endwerten (runde) Skalenteil- */
/* strichwerte (aehnlich wie auf einem Rechenstab) bestimmen.        */
/*                                                                   */
/*     von: <- Anfangswert                                           */
/*     bis: <- Endwert, von > bis oder von < bis moeglich            */
/*  minpkt: <- min. Anzahl von Teilstrichen                          */ 
/*    teil: -> Feld fuer die Aufteilung in die laengeren Teilstriche */
/*    smin: -> Abstand zweier Teilstriche                            */
/*      x1: -> Wert des ersten Teilstriches, von in Richtung bis     */
/*   srest: -> Rest von x1 im Raster kleinste Teilung/groesste Teilung*/
/* ----------------------------------------------------------------- */

skalgen: PROC((von,bis,minpkt) FLOAT, teil ()FIXED IDENT,
              (smin, x1) FLOAT IDENT, srest FIXED IDENT);

DCL (xdiff, xpunkt, skal2) FLOAT;
DCL (ntotal, n25) FIXED;

xdiff = bis - von;
IF xdiff == 0 THEN RETURN; FIN;

smin = SIGN(xdiff) * EXP(LN(10.0) * ENTIER(LG(ABS(xdiff) / minpkt)));
xpunkt = xdiff / smin;

ntotal = 5;
n25 = 2;
FOR i FROM 1 TO 1 UPB teil REPEAT
  teil(i) = ntotal;
  ntotal = ntotal * n25;
  n25 = 10 // n25;
END;
IF xpunkt >= 5 * minpkt THEN
  smin = smin * 5;
  FOR i FROM 1 BY 2 TO 1 UPB teil REPEAT
    teil(i) = teil(i) * 2 // 5;
  END;
  xpunkt = xpunkt / 5;
ELSE IF xpunkt >= 2 * minpkt THEN
  smin = smin * 2;
  FOR i FROM 2 BY 2 TO 1 UPB teil REPEAT
    teil(i) = teil(i) * 5 // 2;
  END;
  xpunkt = xpunkt / 2;
FIN; FIN;

x1 = fmlow(von, -smin);
skal2 = fmlow(von, smin * teil(1 UPB teilung));
srest = ROUND((x1 - skal2) / smin);

END;

/* ----------------------------------------------------------------- */
/* Pixelzeile fuer Drucker loeschen                                  */
/* ----------------------------------------------------------------- */

clrpix: PROC;

  FOR i FROM 1 TO 2496 REPEAT
    FOR j FROM 1 TO 3 REPEAT   /* Bitmuster fuer Zeile loeschen */
      pixline(j, i) = '0'B;
    END;
  END;
END;

/* ----------------------------------------------------------------- */
/* Pixelzeile auf Drucker ausgeben                                   */
/* ----------------------------------------------------------------- */

prpix: PROC (odd BIT(1));

  PUT TOCHAR(27),'$',TOCHAR(32),TOCHAR(0),
      TOCHAR(27),'*',TOCHAR(40),TOCHAR(192),TOCHAR(9) TO Druck;
                                /* Graphik 24 Nadeln 360 Punkte/Zoll */
  FOR j FROM 1 BY 8 TO 2496 REPEAT
    FOR i FROM 0 TO 7 REPEAT
      FOR k FROM 1 TO 3 REPEAT
        IF odd THEN
          pb(3*i + k) = TOCHAR(TOFIXED(pixline(k,j+i) SHIFT r8));
        ELSE
          pb(3*i + k) = TOCHAR(TOFIXED(pixline(k,j+i) AND l8));
        FIN;
      END;
    END;
    PUT pb(1),pb(2),pb(3),pb(4),pb(5),pb(6),
        pb(7),pb(8),pb(9),pb(10),pb(11),pb(12),
        pb(13),pb(14),pb(15),pb(16),pb(17),pb(18),
        pb(19),pb(20),pb(21),pb(22),pb(23),pb(24) TO Druck;
  END;
END;     /* prpix */

/* ----------------------------------------------------------------- */
/* ober bzw. untere horizontale Skala fuer Drucker bestimmen, mit    */
/* Zahlenbeschriftung und Titelzeile oben                            */
/*                                                                   */
/* uskal: <- '0'B obere Skala, '1'B untere Skala                     */
/* ----------------------------------------------------------------- */

hskala: PROC(uskal BIT(1));

DCL (w1, w2, w3, wx) BIT(16);
DCL index FIXED;
DCL xmax2 FLOAT;

IF NOT uskal THEN
  PUT TOCHAR(27), 'j', TOCHAR(toline1), TOCHAR(27), '!', TOCHAR(4),
      TOCHAR(27),'$',TOCHAR(30),TOCHAR(0),  /* 17 cpi ... */
      TOCHAR(27), 'x', TOCHAR(1),           /* NLQ ... */
      'k\x          Poincare', TOCHAR(8), TOCHAR(39),
      '-Abbildung (k,x) -> k*x*(1-x):          ' TO Druck;
  PUT 'Ausschnitt [', ymin, ', ', ymax, '] x [', xmin, ', ', xmax, ']'
      TO Druck BY A, F(7,5), A, F(7,5), A, F(6,4), A, F(6,4), A;
  PUT TOCHAR(27), 'J', TOCHAR(32) TO Druck;
ELSE
  PUT TOCHAR(27), 'J', TOCHAR(2),
      TOCHAR(27), '!', TOCHAR(0) TO Druck;  /* 10 cpi ... */
FIN;
CALL clrpix;                              /* Pixelspeicher loeschen */

CALL skalgen(xmin - xrange / 4000, xmax, 60.0, teilung, skalmin,
             xskal, xrest);             /* Skalenparameter gewinnen */
xmax2 = xmax + skalmin / 2;
IF uskal THEN
  wx = '00F0'B4;
ELSE
  wx = '000F'B4;
FIN;
WHILE xskal <= xmax2 REPEAT
  w1 = '0000'B4;
  w2 = '0000'B4;
  w3 = '00FF'B4;
  IF xrest REM teilung(1) == 0 THEN w2 = wx; FIN;
  IF xrest REM teilung(2) == 0 THEN
    w2 = '00FF'B4;
    IF NOT uskal THEN
                             /* Skala beschriften */
      index = ROUND(pixcols / 6 * (xskal - xmin) / xrange) + 25  ;
      PUT TOCHAR(27),'$',TOCHAR(index REM 256),TOCHAR(index // 256)
          TO Druck;
      PUT xskal TO Druck BY F(10,5);
    FIN;
  FIN;
  IF xrest REM teilung(3) == 0 THEN w1 = wx; FIN;
  IF xrest REM teilung(4) == 0 THEN w1 = '00FF'B4; FIN;
  IF uskal THEN
    wx = w3; w3 = w1; w1 = wx;
  FIN;
  index = ROUND(pixcols * (xskal - xmin) / xrange) + 49  ;
  pixline (1, index) = w1;
  pixline (2, index) = w2;
  pixline (3, index) = w3;

  xskal = fmlow(xskal + skalmin / 2, -skalmin);
  xrest = xrest + 1;

/*PUT xskal, xrest TO Dialog BY SKIP, F(10,5), F(3);
  GET FROM Dialog BY SKIP; */
END;

IF NOT uskal THEN
  PUT TOCHAR(27), 'J', TOCHAR(16) TO Druck;
FIN;

CALL prpix('0'B1);

PUT TOCHAR(27), 'J', TOCHAR(26) TO Druck; /* Druckkopflaenge vor */

END;

/* ----------------------------------------------------------------- */
/* Hauptprogramm                                                     */
/* ----------------------------------------------------------------- */

CHA:TASK;

DCL (dlauf, lteil) FIXED;

PUT TOCHAR(26), TOCHAR(27), '=', TOCHAR(33), TOCHAR(46),
    'Poincare-Abbildung x -> k * x*(1-x): ]0,1[ -> ]0,1[' TO Dialog;

LOOP1:
  PUT TOCHAR(27), '=', TOCHAR(38), TOCHAR(51) TO Dialog;
  PUT'xmin (>= 0): ', TOCHAR(27), 'T' TO Dialog;
  GET xmin FROM Dialog BY SKIP,F(10);
  PUT TOCHAR(27), '=', TOCHAR(39), TOCHAR(51) TO Dialog;
  PUT'xmax (<= 1): ', TOCHAR(27), 'T' TO Dialog;
  GET xmax FROM Dialog BY SKIP,F(10);
  xrange = xmax - xmin;
IF xrange < 0.001 THEN GOTO LOOP1; FIN;

LOOP2:
  PUT TOCHAR(27), '=', TOCHAR(41), TOCHAR(51) TO Dialog;
  PUT'kmin (>= 1): ', TOCHAR(27), 'T' TO Dialog;
  GET ymin FROM Dialog BY SKIP,F(10);
  PUT TOCHAR(27), '=', TOCHAR(42), TOCHAR(51) TO Dialog;
  PUT'kmax (<= 4): ', TOCHAR(27), 'T' TO Dialog;
  GET ymax FROM Dialog BY SKIP,F(10);
  yrange = ymax - ymin;
IF ymin < 1 OR ymax > 4 OR yrange < 0.001 THEN GOTO LOOP2; FIN;

LOOP3:
  PUT TOCHAR(27), '=', TOCHAR(44), TOCHAR(40) TO Dialog;
  PUT'x-Startwert (> 0, < 1): ', TOCHAR(27), 'T' TO Dialog BY A;
  GET xanf FROM Dialog BY SKIP,F(10);
IF xanf <= 0 OR xanf >= 1 THEN GOTO LOOP3; FIN;

PUT TOCHAR(27), '=', TOCHAR(45), TOCHAR(44),
   'Vorlauf   (< 2^14): ', TOCHAR(27), 'T' TO Dialog;
GET tiefemax FROM Dialog BY F(10);
PUT TOCHAR(27), '=', TOCHAR(46), TOCHAR(45),
   'Punkte   (< 2^14): ', TOCHAR(27), 'T' TO Dialog;
GET punkte FROM Dialog BY F(10);

PUT TOCHAR(27), '=', TOCHAR(48), TOCHAR(40),
   'Ausgabe auf Drucker statt Bildschirm?', TOCHAR(27), 'T' TO Dialog;
GET pchar FROM Dialog BY SKIP, A;

IF pchar == 'J' OR pchar == 'j' OR pchar == 'Y' OR pchar == 'y' THEN

                               /* Ausgabe auf den Drucker */
  CALL hskala(oben);           /* horizontale Skala oben anbringen */

  CALL skalgen(ymax - yrange / 2 / pixlines, ymin, 100.0, teilung,
               skalmin, yskal, yrest); /* Skalenparameter gewinnen */
  sline = pixlines * (yskal - ymin) / yrange;

  PUT TOCHAR(28), '3', TOCHAR(1) TO Druck; /* Zeilenabstand 1/360'' */

  FOR lineno FROM prtlines - 1 BY -1 TO 0 REPEAT
                                       /* Druckerzeile (24 Nadeln) */

    CALL clrpix;

    FOR byteno FROM 1 TO 3 REPEAT        /* Druckerbyte (8 Nadeln) */
      FOR bitno FROM 1 TO 8 REPEAT
        pixlno = 24 * lineno - 8 * byteno - bitno + 32;
        ywert = ymin + yrange * pixlno / pixlines;

        PUT TOCHAR(27), '=', TOCHAR(50), TOCHAR(40),
            pixlno, ywert TO Dialog BY (4)A, F(8), F(10,5);

                                     /* vertikale Skala anbringen: */
        IF pixlno <= ROUND(sline) THEN
          IF pixlno > sline THEN
            dlauf = 8;
              ELSE
            dlauf = 0;         /* Entscheidung fuer Druckdurchlauf */
          FIN;
          lteil = 12;             /* Laenge des Skalenteilstriches */
          FOR m FROM 1 TO 1 UPB teilung WHILE yrest REM teilung(m) == 0
          REPEAT
            lteil = lteil + 8; /* bestimmte Teilstriche sind laenger */
          END;

          IF yrest REM teilung(3) == 0 THEN
            PUT TOCHAR(27),'$',TOCHAR(8),TOCHAR(0) TO Druck;
            PUT yskal TO Druck BY F(7,5);
          FIN;
                              /* Teilstrich in die Pixelzeile setzen */
          FOR m FROM 1 TO lteil REPEAT
            pixline(byteno, skalcols-2 - m).BIT(bitno + dlauf) = '1'B;
            pixline(byteno, pcr+2 + m).BIT(bitno + dlauf) = '1'B;
          END;

          yskal = fmlow(yskal + skalmin / 2, -skalmin);
          sline = pixlines * (yskal - ymin) / yrange;
          yrest = yrest + 1;
        FIN;

        xwert=xanf;
        FOR ii FROM 1 TO tiefemax+punkte REPEAT
                               /* die eigentliche Poincare-Abbildung */
          xwert= ywert * xwert *(1-xwert);

          IF ii > tiefemax THEN
            xpix = pixcols * (xwert - xmin) / xrange + skalcols + 1;
            IF xpix >= pcl AND xpix < pcr THEN
              IF pixline(byteno, ENTIER(xpix)).BIT(bitno) THEN
                pixline(byteno, ENTIER(xpix)).BIT(bitno+8) = '1'B;
              ELSE
                pixline(byteno, ENTIER(xpix)).BIT(bitno) = '1'B;
              FIN;
            FIN;
          FIN;
        END; /* ii */

      END; /* bitno */
    END;   /* byteno */

    CALL prpix('1'B1);

    PUT TOCHAR(10) TO Druck;                  /* Vorschub 1/360 Zoll */

    CALL prpix('0'B1);

    PUT TOCHAR(10) TO Druck;                  /* Vorschub 1/360 Zoll */
    PUT TOCHAR(27), 'J', TOCHAR(23) TO Druck; /* Vorschub 46/380Zoll */

  END;     /* lineno */

  CALL hskala(unten);         /* horizontale Skala unten anbringen */
  PUT TOCHAR(27), '2', TOCHAR(12) TO Druck;
                    /* Zeilenabstand auf Normalwert, Papierauswurf */
                    /* damit ist die Druckerausgabe beendet        */
ELSE
                    /* dasselbe auf dem Bildschirm:                */
  CALL SCRSET('CHAOS');
  CALL CLEAR;
  FOR zeile FROM 399 BY -1 TO 0 REPEAT
    ywert=ymin+ (ymax-ymin)*(399-zeile)/400;
    xwert=xanf;
    FOR ii FROM 1 TO tiefemax+punkte REPEAT
                             /* die eigentliche Poincare-Abbildung */
      xwert= ywert * xwert *(1-xwert);

      IF (ii > tiefemax) AND (xwert >= xmin) AND (xwert < xmax) THEN
        CALL SETPIX(ENTIER(640*(xwert-xmin)/xrange),zeile,1);
      FIN;
    END;
  END;
FIN;

PUT TOCHAR(27), '=', TOCHAR(48), TOCHAR(84), 'fertig.' TO Dialog;

END;

MODEND;
@