OVERLAY PROCEDURE FFT_MENU;

    { This procedure displays the driving menu for the FFT algorithm.
      The FFT_PROCEDURE procedure is called from this menu.
      The entire FFT program is a subprocedure to the FFT_MENU.
      This facilitates the concept of overlaying the programs }

  PROCEDURE FFT_PROCEDURE;

    CONST

    { "Tuning constants" to adjust accuracy, space requirements and
       speed of execution. }

      FRACTIONBIGGER =  0.2;
      WORLDX         =   31;      { SIZE OF THE WORLD }
      WORLDY         =   31;
      WORLDMAX       =   31;      { UPPER LIMIT TO EITHER WORLDX OR WORLDY }
      WORLDSIZE      = 1024;      { (WORLDX+1) * (WORLDY+1) }

    { The following unusual type declaration allows WORLD to be
      treated as either two-dimensional (D2) or one-dimensional (D1) }

    TYPE
      ARRAYD1 = ARRAY[1..WORLDSIZE] OF REAL;
      ARRAYD2 = ARRAY[0..WORLDX, 0..WORLDY] OF REAL;
      WORLDARRAY = RECORD
                     CASE BOOLEAN OF
                       TRUE : (D1 : ARRAYD1);
                       FALSE: (D2 : ARRAYD2)
                   END;

    { Global data known to main procedure and all sub procedures }

    VAR
      NPOINTS : INTEGER;          { NUMBER OF POINTS TO PROCESS }

      SUMOFWORLD  : REAL;
      LOWXWORLD   : REAL;
      HIGHXWORLD  : REAL;
      LOWYWORLD   : REAL;
      HIGHYWORLD  : REAL;

      WORLD       : WORLDARRAY;

    PROCEDURE OUTMATRIX (VAR OUT:ARRAYD2);

      { Writes the utilization distribution matrix out to the work
        file. This file is used later to create the maps. }

      VAR
        C     : INTEGER;
        I     : INTEGER;
        J     : INTEGER;

      BEGIN  { OUTMATRIX PROCEDURE }

        WRITELN;
        WRITELN ('OUTMATRIX');

        FOR I := 0 TO WORLDX DO
          BEGIN
            C := 0;
            FOR J := 0 TO WORLDY DO
              BEGIN
                C := C + 1;
                WRITE (FILEVAR5,' ',OUT[I,J]:10:3);
                IF C = 8 THEN
                  BEGIN
                    WRITELN (FILEVAR5);
                    C := 0
                  END;
              END;  { FOR J = 0 TO WORLDY }

            IF C <> 0 THEN WRITELN (FILEVAR5);

          END; { FOR I = 0 TO WORLDX }

      END;  { OUTMATRIX PROCEDURE }

    PROCEDURE MAP;

      { Calculates the map index from a UD which has first been sorted
        by SHELLSORTUD.
        2) Extracts the sorted value at the map index for isoline plotting. }

      CONST
        HEADING1 = '    MAP %      MAP         WORLD      PERCENT OF';
        HEADING2 = '      TOTAL       NUMBER OF';
        HEADING3 = '              VALUE        VALUE      TOTAL AREA';
        HEADING4 = '       AREA        SQUARES';

      VAR
        I     : INTEGER;
        J     : INTEGER;
        NSQRS : INTEGER;

        PROP         : REAL;
        AREAOFWORLD  : REAL;
        LIMIT        : REAL;
        SUM          : REAL;
        PVALUE       : REAL;
        MAP          : REAL;

      BEGIN  { MAP PROCEDURE }


        AREAOFWORLD := (HIGHXWORLD-LOWXWORLD) * (HIGHYWORLD-LOWYWORLD);

        { PRINT THE MAP HEADINGS }

        IF UPCASE(PRINTCHAR) IN ['D','B'] THEN
          BEGIN
            CLRSCR;
            GOTOXY(1,1);
            WRITELN (HEADING1,HEADING2);
            WRITELN (HEADING3,HEADING4);
          END;  { IF PRINTCHAR = D OR B }

        IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
          BEGIN
            WRITELN (LST,HEADING1,HEADING2);
            WRITELN (LST,HEADING3,HEADING4);
          END;  { IF PRINTCHAR = P OR B }

        J := 5;
        WHILE J < 100 DO
          BEGIN
            PVALUE:= J/100;
            J := J + 5;
            LIMIT := (1.0-PVALUE) * SUMOFWORLD;
            I := 0;
            SUM := 0.0;

            REPEAT
              I := I + 1;
              SUM := SUM + WORLD.D1[I];
            UNTIL ( SUM > LIMIT ) OR ( I=WORLDSIZE );


            IF SUM > LIMIT THEN
              BEGIN
                I := I - 1;
                MAP := (WORLDSIZE-I) * AREAOFWORLD / WORLDSIZE;
                PROP:= (WORLDSIZE-I)/WORLDSIZE;
                NSQRS:= WORLDSIZE-I;

                { WRITE / DISPLAY THE TABLE INFORMATION }

                IF UPCASE(PRINTCHAR) IN ['D','B'] THEN
                  BEGIN
                    WRITE ('   (',PVALUE:3:2,')   ',MAP:8:5,'    ');
                    WRITE (WORLD.D1[I]:10:6,'     ',PROP:6:4,'      ');
                    WRITE (AREAOFWORLD:10:6,'       ',NSQRS:3);
                    WRITELN;
                  END;  { IF PRINTCHAR = D OR B }
                IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
                  BEGIN
                    WRITE (LST,'   (',PVALUE:3:2,')   ',MAP:8:5,'    ');
                    WRITE (LST,WORLD.D1[I]:10:6,'     ',PROP:6:4,'      ');
                    WRITE (LST,AREAOFWORLD:10:6,'       ',NSQRS:3);
                    WRITELN (LST);
                  END;  { IF PRINTCHAR = P OR B }

              END;  { IF SUM > LIMIT ... }

          END;  { WHILE J < 100 }

        KEYBOARD_PAUSE

      END;  { MAP PROCEDURE }

  PROCEDURE SHELLSORTUD (VAR WORLD:ARRAYD1);

    { SHELLSORT.  Treats the UD stored in world as a one dimensional
     array and sorts it in increasing order. This is necessary before
     calling map. Adapted from KNUTH xxxx        }

    CONST
      T = 4;
      N = WORLDSIZE;

    TYPE
      INDEX =-15..WORLDSIZE;
      ITEM = REAL;

    VAR
      I       : INDEX;
      J       : INDEX;
      K       : INDEX;
      S       : INDEX;

      X       : ITEM;

      M       : 1..T;

      ARAY    : ARRAY[INDEX] OF REAL;

      H       : ARRAY [1..T] OF INTEGER;

    BEGIN

      WRITELN('SORTING');

      FOR I:= 1 TO WORLDSIZE DO
        ARAY[I] := WORLD[I];
      H[1]:=15;
      H[2]:=7;
      H[3]:=3;
      H[4]:=1;

      { SORT }

      FOR M := 1 TO T DO
        BEGIN
          K := H[M];
          S := -K;
          FOR I := K+1 TO N DO
            BEGIN
              X := ARAY[I];
              J := I-K;
              IF S=0 THEN S := -K;
              S := S+1;
              ARAY[S] :=X;
              WHILE X < ARAY[J] DO
                BEGIN
                  ARAY[J+K] := ARAY[J];
                  J := J-K;
                END;
              ARAY[J+K] := X;
            END;  { FOR I = K+1 TO N }

        END; { FOR M = 1 TO T }

      FOR I := 1 TO WORLDSIZE DO
        WORLD[I] := ARAY[I];
      SUMOFWORLD := 0.0;
      FOR I := 1 TO WORLDSIZE DO
        SUMOFWORLD := SUMOFWORLD + WORLD[I];

    END;  { SHELLSORTUD }

  PROCEDURE DENSITY;

    { Estimates the UD from the ovservations stored in variable
      WORLD using the FOURIER TRANSFROM METHOD (ANDERSON, 1980).
      WORLD and NPOINTS are used. }

    CONST
      { CUTOFFCOUNT should normally be set to 2.  When a CUTOFFCOUNT
     by CUTOFFCOUNT square of FOURIER COEFFICIENTS drops below
     cutoff, they are set to 0.  Detailed explanation is found at the
     end of APPENDIX A (ANDERSON, 1981) }

      CUTOFFCOUNT = 2;

    TYPE
      COMPLEX = RECORD
                  RE : REAL;
                  IM : REAL;
                END;

      { DIRECT = FORWARD FFT, INVERSE = INVERSE FFT }

       DIRECTINVERSE = (DIRECT, INVERSE);
       XROWYCOLUMN = (XROW, YCOLUMN);
       INDEX = RECORD
                 L : INTEGER;
                 N : INTEGER;
               END;

    VAR
      X       : INTEGER;
      Y       : INTEGER;
      NSWAP   : INTEGER;

      NSQUARED : REAL;

      IMAGWORLD : ARRAYD2;
      COEFFTABLE : ARRAY[0..WORLDMAX] OF COMPLEX;
      SWAPINDEX : ARRAY[0..WORLDMAX] OF INDEX;

    PROCEDURE FFT (FFTDIRECTION:XROWYCOLUMN;
                   ROWORCOLUMN:INTEGER);

      { RECURSIVE FAST FOURIER TRANSFORM.  SEE FIG. 6.10 P297
        OPPENHEIM & SCHAFER (1975) }

    VAR
      LENGTH      : INTEGER;
      LENGTHOVER2 : INTEGER;

      PROCEDURE SUBDIVIDE (START, HALFSIZE:INTEGER);

        (*RECURSIVELY SUBDIVIDE DATA*)

        PROCEDURE SWARM;

          { CALCULATE A SWARM OF BUTTERFLIES }

        VAR
          RINC    : INTEGER;
          R       : INTEGER;
          P       : INTEGER;
          I1      : INTEGER;
          I2      : INTEGER;
          J1      : INTEGER;
          J2      : INTEGER;

          XP          : COMPLEX;
          XQ          : COMPLEX;
          TERM        : COMPLEX;
          COEFFICIENT : COMPLEX;

         BEGIN  { SWARM }

           RINC := LENGTHOVER2 DIV HALFSIZE;
           R := 0;
           FOR P := START TO (START + HALFSIZE - 1) DO
             BEGIN
               COEFFICIENT := COEFFTABLE[R];
               CASE FFTDIRECTION OF
                 XROW : BEGIN
                        I1 := P;
                        J1 := ROWORCOLUMN;
                        I2 := P + HALFSIZE;
                        J2 := ROWORCOLUMN
                        END;
                YCOLUMN : BEGIN
                        I1 := ROWORCOLUMN;
                        J1 := P;
                        I2 := ROWORCOLUMN;
                        J2 := P + HALFSIZE
                        END;
                 END; { CASE FFTDIRECTION }

               XP.RE := WORLD.D2[I1,J1];
               XP.IM := IMAGWORLD[I1,J1];
               XQ.RE := WORLD.D2[I2,J2];
               XQ.IM := IMAGWORLD[I2,J2];
               TERM.RE := COEFFICIENT.RE*XQ.RE - COEFFICIENT.IM*XQ.IM;
               TERM.IM := COEFFICIENT.RE*XQ.IM + COEFFICIENT.IM*XQ.RE;
               WORLD.D2[I1,J1] := XP.RE + TERM.RE;
               IMAGWORLD[I1,J1] := XP.IM + TERM.IM;
               WORLD.D2[I2,J2] := XP.RE - TERM.RE;
               IMAGWORLD[I2,J2] := XP.IM - TERM.IM;
               R := R + RINC

             END;  { FOR P = START TO ... }

         END;  { SWARM }

       BEGIN  { SUBDIVIDE PROCEDURE }

         IF HALFSIZE > 1 THEN
           BEGIN
             SUBDIVIDE (START, HALFSIZE DIV 2);
             SUBDIVIDE (START+HALFSIZE, HALFSIZE DIV 2)
           END;

        SWARM;

      END;  { SUBDIVIDE PROCEDURE }

      PROCEDURE BITREVERSE;
      VAR
        N, I1, I2, J1, J2 : INTEGER;
        SWAP : REAL;

       BEGIN  { BITREVERSE PROCEDURE }

         FOR N := 0 TO NSWAP DO
           BEGIN
             CASE FFTDIRECTION OF
               XROW : BEGIN
                      I1 := SWAPINDEX[N].L;
                      J1 := ROWORCOLUMN;
                      I2 := SWAPINDEX[N].N;
                      J2 := ROWORCOLUMN;
                      END;
              YCOLUMN : BEGIN
                      I1 := ROWORCOLUMN;
                      J1 := SWAPINDEX[N].L;
                      I2 := ROWORCOLUMN;
                      J2 := SWAPINDEX[N].N;
                      END;
             END;  { CASE FFTDIRECTION }

             SWAP := WORLD.D2[I1,J1];
             WORLD.D2[I1,J1] := WORLD.D2[I2,J2];
             WORLD.D2[I2,J2] := SWAP;
             SWAP := IMAGWORLD[I1,J1];
             IMAGWORLD[I1,J1] := IMAGWORLD[I2,J2];
             IMAGWORLD[I2,J2] := SWAP ;

           END;  { FOR I = 0 TO NSWAP }

       END;  { BITREVERSE PROCEDURE }

     BEGIN  { FFT PROCEDURE }

       WRITE ('f.');
       CASE FFTDIRECTION OF
         XROW : LENGTH := WORLDX +1;
         YCOLUMN : LENGTH := WORLDY +1
       END;  { CASE FFTDIRECTION }

       BITREVERSE;

       LENGTHOVER2 := LENGTH DIV 2;

       SUBDIVIDE (0, LENGTHOVER2);

     END;  { FFT PROCEDURE }

    PROCEDURE SETUP (KIND : DIRECTINVERSE;
                     SIZE : INTEGER);

      CONST
        TWOPI = 6.2831853;

      VAR
        R        : INTEGER;
        N        : INTEGER;
        L        : INTEGER;
        LOG2SIZE : INTEGER;

        TWOPIOVERSIZE : REAL;
        THETA         : REAL;

      FUNCTION REVERSE (N:INTEGER):INTEGER;

        { BIT REVERSE N }

        VAR
          I      : INTEGER;
          RESULT : INTEGER;

       BEGIN  { REVERSE FUNCTION }

         WRITE ('r');
         RESULT := 0;
         FOR I := 1 TO LOG2SIZE DO
           BEGIN
             RESULT := RESULT*2;
             IF (N MOD 2) <> 0  THEN RESULT := RESULT +1;
             N := N DIV 2;
           END;

         REVERSE := RESULT;

       END;  { REVERSE FUNCTION }

     BEGIN  { SETUP PROCEDURE }

       { CALCULATE SINE AND COSINE COEFFICIENTS }

       WRITELN ('BEGINNING SETUP');

       TWOPIOVERSIZE := TWOPI / SIZE;
       FOR R := 0 TO (SIZE - 1) DO
         BEGIN
           THETA := R * TWOPIOVERSIZE;
           COEFFTABLE[R].RE := COS (THETA);
           CASE KIND OF
             DIRECT  : COEFFTABLE[R].IM :=  SIN(THETA);
             INVERSE : COEFFTABLE[R].IM := -SIN (THETA)
           END; { CASE KIND OF }

         END; { FOR R = O TO SIZE-1 }

       { CALCULATE LOG BASE 2 OF SIZE }

       N := SIZE;
       LOG2SIZE := 0;
         REPEAT
           N := N DIV 2;
           LOG2SIZE := LOG2SIZE +1
         UNTIL N = 1;

       { CALCULATE COEFFICIENTS FOR BIT REVERSAL }

       NSWAP := -1;
       FOR N := 0 TO (SIZE-1) DO
         BEGIN
           L := REVERSE(N);
           IF L > N THEN
             BEGIN
               NSWAP := NSWAP + 1;
               SWAPINDEX[NSWAP].L := L;
               SWAPINDEX[NSWAP].N := N
             END; { IF L > N }

         END;  { FOR N = 0 TO SIZE - 1 }

         WRITELN;

     END;  { SETUP PROCEDURE }

    PROCEDURE EXCISE (XBOTTOM,XTOP,YBOTTOM,YTOP:INTEGER);

      { Selectively excises FOURIER COEFFICIENTS within the rectangular
       region xbottom to xtop and ybottom to ytop.  If any
       CUTOFFCOUNT by CUTOFFCOUNT sized square of FOURIER COEFFICIENTS
       within the region drops below cutoff then all higher frequency
       coefficients in the region are set to zero.

       I'M NOT HAPPY WITH THIS PROCEDURE.  IT IS TOO LONG AND
       COMPLICATED CONSIDERING WHAT IT DOES.  IT COULD PROBABLY
       BE REWRITTEN TO IMPROVE CLARITY.  ALSO, IT MIGHT BE POSSIBLE
       TO MAKE USE OF THE FACT THAT OF THE 4 CALLS TO EXCISE, 2
       CARRY OUT IDENTICAL CALCULATIONS IN DIFFERENT REGIONS. }

    TYPE
      XY = RECORD
             X : INTEGER;
             Y : INTEGER;
           END;

    VAR
      I       : INTEGER;
      X       : INTEGER;
      Y       : INTEGER;
      LAST    : INTEGER;
      XINC    : INTEGER;
      YINC    : INTEGER;
      OLDXTOP : INTEGER;
      OLDYTOP : INTEGER;

      CUTOFF  : REAL;

      { STORAGE FOR THE X,Y LOCATIONS ABOVE WHICH THE FOURIER
       COEFFICIENTS MUST BE SET TO ZERO. }

      STOREXY : ARRAY[1..WORLDMAX] OF XY;

      FUNCTION ZEROSQUARE (X,Y:INTEGER):BOOLEAN;

        { Checks to see if the square X to X+CUTOFFCOUNT,
         Y TO Y+CUTOFFCOUNT has dropped below cutoff and needs
         to be set to zero. }

      VAR
        XEND  : INTEGER;
        YEND  : INTEGER;
        OLDX  : INTEGER;

      BEGIN  { ZEROSQUARE FUNCTION }

        ZEROSQUARE := TRUE;
        XEND := X + XINC*CUTOFFCOUNT;
        YEND := Y + YINC*CUTOFFCOUNT;
        OLDX := X;
        REPEAT
          X := OLDX;
          REPEAT
            IF ( SQR(WORLD.D2[X,Y]) + SQR(IMAGWORLD[X,Y]) ) >= CUTOFF
              THEN ZEROSQUARE := FALSE;
            X := X+XINC;
          UNTIL X = XEND;
          Y := Y + YINC;
        UNTIL Y = YEND;

      END; { ZEROSQUARE FUNCTION }

    BEGIN  { EXCISE }

      WRITELN;
      WRITELN ('EXCISE');

      CUTOFF := ( 2 / (NPOINTS+1) ) * NPOINTS * NPOINTS;

      { NOTE THAT CUTOFF IS ADJUSTED BY SQR(NPOINTS) BECAUSE
        THE FOURIER TRANSFORM WAS NOT DIVIDED BY NPOINTS }

      IF XTOP > XBOTTOM
        THEN XINC := 1
        ELSE XINC := -1;
      IF YTOP > YBOTTOM
        THEN YINC := 1
        ELSE YINC := -1;
      LAST := 0;

      { SAVE OLD TOP FOR LATER }

      OLDXTOP := XTOP;
      OLDYTOP := YTOP;

      { ADJUST TOP BECAUSE THE HIGHEST CUTOFFCOUNT BY CUTOFFCOUNT
        SQUARE SHOULD NOT BE ABOVE ORIGINAL XTOP AND YTOP. }

      XTOP := XTOP + XINC*(2-CUTOFFCOUNT);
      YTOP := YTOP + YINC*(2-CUTOFFCOUNT);
      Y := YBOTTOM;

      { SCAN THROUGH THE WORLD CALCULATING THE X,Y COORDINATES
        ABOVE WHICH THE FOURIER COEFFICIENTS MUST BE SET TO ZERO. }

      WHILE Y <> YTOP DO
        BEGIN
          X := XBOTTOM;
          WHILE X <> XTOP DO
            BEGIN
              IF ( SQR(WORLD.D2[X,Y]) + SQR(IMAGWORLD[X,Y]) ) < CUTOFF THEN
                BEGIN
                  WORLD.D2[X,Y] := 0.0;
                  IMAGWORLD[X,Y] := 0.0;
                  IF ZEROSQUARE(X,Y) THEN
                    BEGIN
                      LAST := LAST + 1;
                      STOREXY[LAST].X := X;
                      STOREXY[LAST].Y := Y;
                      XTOP := X
                    END
                  ELSE
                    X := X + XINC;
                END
              ELSE
                X := X + XINC;
            END;

          Y := Y + YINC
        END;

      { RESTORE OLD TOP }

      XTOP := OLDXTOP + XINC;
      YTOP := OLDYTOP + YINC;

      { SET COEFFICIENTS TO ZERO *}

      FOR I := 1 TO LAST DO
        BEGIN
          Y := STOREXY[I].Y;
          WHILE Y <> YTOP DO
            BEGIN
              X := STOREXY[I].X;
              WHILE (X*XINC) < (XTOP*XINC) DO
                BEGIN
                  WORLD.D2[X,Y] := 0.0;
                  IMAGWORLD[X,Y] := 0.0;
                  X := X + XINC
                END;
              Y := Y + YINC
            END; { WHILE Y <> YTOP }

          XTOP := STOREXY[I].X;

        END; { FOR I = 1 TO LAST }

      END; { EXCISE PROCEDURE }

    BEGIN  { DENSITY PROCEDURE }

      WRITELN ('BEGINNING DENSITY ');

      { ZERO OUT IMIGINARY DATA }

       FOR X :=0 TO WORLDX DO
         FOR Y :=0 TO WORLDY DO
           IMAGWORLD[X,Y] := 0;

      { DO 2 DIMENSIONAL DIRECT (OR FORWARD) FFT }
      { DO ROWS }

      SETUP (DIRECT, WORLDX+1 );

      FOR X :=0 TO WORLDX DO
        FFT (XROW, X);

      { THEN COLUMNS }

      IF WORLDY <> WORLDX THEN
        SETUP (DIRECT, WORLDY+1 );

      FOR Y := 0 TO WORLDY DO
        FFT (YCOLUMN, Y);

      { FILTER BY EXCISING COEFFICIENTS WITHIN SPECIFIED REGION }

      EXCISE ( 0,      (WORLDX DIV 2),      0,      (WORLDY DIV 2)     );
      EXCISE ( WORLDX, ((WORLDX+1) DIV 2),  0,      (WORLDY DIV 2)     );
      EXCISE ( 0,      (WORLDX DIV 2),      WORLDY, ((WORLDY+1) DIV 2) );
      EXCISE ( WORLDX, ((WORLDX+1) DIV 2),  WORLDY, ((WORLDY+1) DIV 2) );

      { THE DIRECT FFT DOES NOT DIVIDE BY N SO DO IT NOW }

      NSQUARED := (WORLDX+1) * (WORLDY+1);
      FOR X := 0 TO WORLDX DO
        FOR Y := 0 TO WORLDY DO
          BEGIN
            WORLD.D2[X,Y] := WORLD.D2[X,Y] / NSQUARED;
            IMAGWORLD[X,Y] := IMAGWORLD[X,Y] / NSQUARED
          END;

      { DO TWO DIMENSTIONAL INVERSE FFT }

     SETUP (INVERSE, WORLDY+1 ); (*DO COLUMNS*)
     FOR Y := 0 TO WORLDY DO
     FFT (YCOLUMN, Y);
     IF WORLDX<>WORLDY  (*THEN ROWS*)
      THEN
      SETUP (INVERSE, WORLDX+1 );
     FOR X := 0 TO WORLDX DO
     FFT (XROW, X);
   END (*DENSITY*);

  PROCEDURE READDATA;

    { READS DATA (X,Y LOCATIONS) FROM INPUT DATA FILE. }

    CONST
      MAXWARN = 5;                { MAXIMUM NUMBER OF WARNINGS ISSUED }

    VAR
      XINDEX    : INTEGER;
      YINDEX    : INTEGER;
      I         : INTEGER;
      ERRORCOUNT : INTEGER;       { NUMBER OF ERRORS THAT HAVE OCCURRED }
      NPOINTS   : INTEGER;        { NUMBER OF OBSERVATIONS READ }

    BEGIN  { READDATA PROCEDURE }

      { Load the min and max values from second and third records of data }

      WRITELN;
      WRITELN('Loading min and max values.');

     (*INITIALIZE VARIABLES*)

      SEEK (FILEVAR2,1);
      READ (FILEVAR2,FILE2REC);

      XLOW  := FILE2REC.X;
      YLOW  := FILE2REC.Y;
      ZLOW  := FILE2REC.Z;
      NPOINTS := FILE2REC.RECORD_NO;

      READ (FILEVAR2,FILE2REC);        { GET THE MAXIMUM VALUES }

      XHIGH := FILE2REC.X;
      YHIGH := FILE2REC.Y;
      ZHIGH := FILE2REC.Z;

      IF XLOW = XHIGH THEN             { HOME RANGE IS A VERTICAL LINE }
        BEGIN
          XLOW := XLOW - 1.0;          { OFFSET THE GRID BOUNDRIES SO PROGRAM }
          XHIGH := XHIGH + 1.0;        { WILL WORK. }
          LINE_OUT(1,24,SPACE79);
          LINE_OUT(1,25,SPACE79);
          LINE_OUT(1,24,'WARNING: All points have same X location');
          GOTOXY(1,25);
          KEYBOARD_PAUSE;
        END;

      IF YLOW = YHIGH THEN             { HOME RANGE IS A HORIZONTAL LINE }
        BEGIN
          YLOW := YLOW - 1.0;
          YHIGH := YHIGH + 1.0;
          LINE_OUT(1,24,SPACE79);
          LINE_OUT(1,25,SPACE79);
          LINE_OUT(1,24,'WARNING: All points have the same Y location');
          GOTOXY(1,25);
          KEYBOARD_PAUSE;
        END;  { YLOW = YHIGH }

      { ASK USER IF THEY WANT TO USE THE STUDY AREA OVERLAY TO DEFINE GRID }

      IF STUDYFILE <> '' THEN
       BEGIN
         LINE_OUT(1,25,SPACE79);
         LINE_OUT(1,25,'Do you want to use the STUDY FILE for grid boundry definition? (Y|N) ');
         KEYVALUE := ' ';
         WHILE NOT (UPCASE(KEYVALUE) IN ['Y','N']) DO
           READ    (KBD,KEYVALUE);

         KEYVALUE := UPCASE(KEYVALUE);

         IF KEYVALUE = 'Y' THEN
           BEGIN
             ASSIGN (FILEVAR3,STUDYFILE);
             RESET  (FILEVAR3);
             SEEK   (FILEVAR3,1);
             READ   (FILEVAR3,FILE3REC);

             XLOW := MIN_VAL(FILE3REC.X1,XLOW);
             YLOW := MIN_VAL(FILE3REC.Y1,YLOW);
             ZLOW := MIN_VAL(FILE3REC.Z1,ZLOW);
             XHIGH := MAX_VAL(FILE3REC.X2,XHIGH);
             YHIGH := MAX_VAL(FILE3REC.Y2,YHIGH);
             ZHIGH:= MAX_VAL(FILE3REC.Z2,ZHIGH);

           END;  { IF KEYVALUE = Y }

       END;  { IF STUDYFILE <> '' };

      { If a study area is used as input to the user definition of the boundry
        then we need to determine the four extreme points of the study area. }

      { <<<<<<<<<<<< Later, install code to find four extremes of the study
        area file. >>>>>>>>>>>>>>>>>>>>}

      { Determine the delta x and delta y for the grid }

      CLRSCR;
      WRITELN;
      WRITELN('Outer limits of the grid are: ');
      WRITELN('   X: ',XLOW:10:5,' TO ',XHIGH:10:5);
      WRITELN('   Y: ',YLOW:10:5,' TO ',YHIGH:10:5);
      WRITELN('There are ',NPOINTS,' data points.');

 {    CLRSCR;   }
      LOWXWORLD  := XLOW  - (XHIGH-XLOW) * FRACTIONBIGGER;
      HIGHXWORLD := XHIGH + (XHIGH-XLOW) * FRACTIONBIGGER;
      LOWYWORLD  := YLOW  - (YHIGH-YLOW) * FRACTIONBIGGER;
      HIGHYWORLD := YHIGH + (YHIGH-YLOW) * FRACTIONBIGGER;

      { WORLDX AND WORLDY ARE THE SAME AS GRIDSIZE IN HISTOGRAM }

      XDELTA := (HIGHXWORLD - LOWXWORLD) / WORLDX;
      YDELTA := (HIGHYWORLD - LOWYWORLD) / WORLDY;

      NPOINTS := 0;
      ERRORCOUNT := 1;

      { SET WORLD GRID VALUES TO ALL 0.0 };

      FOR I := 1 TO WORLDSIZE DO
        WORLD.D1[I] := 0.0;

      { Second pass is to put the elements into the grid cells. }

      WRITELN;
      WRITELN ('SECOND PASS, GETTING FREQUENCY OF CELL COUNTS' );

      SEEK (FILEVAR2,3);                { SET PAST THE FILE HEADER }

      WHILE NOT EOF(FILEVAR2) DO
        BEGIN
          READ  (FILEVAR2,FILE2REC);
          WITH FILE2REC DO
            BEGIN
              XINDEX := ROUND ((X-LOWXWORLD) / XDELTA);
              YINDEX := ROUND ((Y-LOWYWORLD) / YDELTA);

              { IF POINT IS WITH IN THE GRID, ADD IT TO THE GRID FREQUENCIES }

              IF (XINDEX IN [0..WORLDX]) AND (YINDEX IN [0..WORLDX]) THEN
                BEGIN
                  NPOINTS := NPOINTS + 1;
                  IF ( NPOINTS = 1 ) OR ( X < XLOW ) THEN XLOW  := X;
                  IF ( NPOINTS = 1 ) OR ( X > XHIGH) THEN XHIGH := X;
                  IF ( NPOINTS = 1 ) OR ( Y < YLOW ) THEN YLOW  := Y;
                  IF ( NPOINTS = 1 ) OR ( Y > YHIGH) THEN YHIGH := Y;
                  WORLD.D2[XINDEX,YINDEX] := WORLD.D2[XINDEX,YINDEX] + 1.0
                END  { IF INDEXES WITHIN THE WORLD }
              ELSE
                IF ERRORCOUNT <= MAXWARN THEN
                  BEGIN

                    { DISPLAY/LIST ERROR MESSAGES }

                    IF UPCASE(PRINTCHAR) IN ['D','B'] THEN
                      WRITELN (OUTPUT,' ?HOMERANGE-W-LOCATION (', X:1:3, ',',
                          Y:1:3, ') IGNORED; IT IS OUTSIDE OF WORLD');
                    IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
                      WRITELN (LST,' ?HOMERANGE-W-LOCATION (', X:1:3, ',',
                          Y:1:3, ') IGNORED; IT IS OUTSIDE OF WORLD');

                    ERRORCOUNT := ERRORCOUNT + 1;

                    IF ERRORCOUNT > MAXWARN THEN
                      BEGIN

                        { DISPLAY/LIST FINAL WARNING }

                        IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
                          WRITELN (OUTPUT,' ?HOMERANGE-W-TOO MANY WARNINGS; ',
                           'FUTURE WARNINGS WON''T BE PRINTED');
                        IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
                          WRITELN (LST,' ?HOMERANGE-W-TOO MANY WARNINGS; ',
                           'FUTURE WARNINGS WON''T BE PRINTED');
                      END;

                  END;  { IF ERRORCOUNT <= MAXWARN }

            END; { WITH FILE2REC DO }

        END;  { WHILE NOT EOF(FILEVAR2) }

      IF XLOW = XHIGH THEN
        BEGIN
          XLOW  := XLOW  - 1.0;
          XHIGH := XHIGH + 1.0;
        END;
      IF YLOW = YHIGH THEN
        BEGIN
          YLOW  := YLOW  - 1.0;
          YHIGH := YHIGH + 1.0;
        END;
      AREA := (XHIGH-XLOW) * (YHIGH-YLOW);
      AREA := AREA / (( HIGHXWORLD - LOWXWORLD ) * ( HIGHYWORLD - LOWYWORLD ));

      IF (AREA < 0.2) THEN
        BEGIN

          { DISPLAY / LIST MESSAGES }

          IF UPCASE(PRINTCHAR) IN ['D','B'] THEN
            WRITELN (OUTPUT,' ?HOMERANGE-W-OBSERVATIONS FELL WITHIN ONLY ',
                   (AREA*100.0):1:1,'& OF THE WORLD');
          IF UPCASE(PRINTCHAR) IN ['P','B'] THEN
            WRITELN (LST,' ?HOMERANGE-W-OBSERVATIONS FELL WITHIN ONLY ',
                   (AREA*100.0):1:1,'& OF THE WORLD');

        END;  { IF AREA < .2 }

     END;  { READDATA PROCEDURE }


  BEGIN  { FFT_PROGRAM PROCEDURE }

    { OPEN AND SET THE FILES }

    ASSIGN  (FILEVAR2,FILEOUT);
    RESET   (FILEVAR2);

    { READ THE DATA }

    READDATA;
    WRITELN ('END OF READDATA');

    { ESTIMATE THE UD }

    DENSITY;

    { The UD is now available in WORLD.D2 to be output if desired  }

    IF UDRESTART <> '' THEN
      BEGIN
        ASSIGN (FILEVAR5,UDRESTART);
        REWRITE(FILEVAR5);
        UDOPEN := TRUE;
        OUTMATRIX (WORLD.D2);
      END;

    SHELLSORTUD (WORLD.D1);

    { MAP(P) may be called for other values of P if desired  }

    MAP;

    KEYVALUE := 'R';

    CLOSE (FILEVAR2);
    IF UDOPEN THEN CLOSE (FILEVAR5);

  END;  { FFT_PROCEDURE PROCEDURE }

  BEGIN  { FFT_MENU PROCEDURE }

    { Display the options menu for FFT_PROCEDURE }

    KEYVALUE := 'R';

    WHILE KEYVALUE = 'R' DO
    BEGIN

    CLRSCR;
    DISPLAY_COLOR(CYAN);

    LINE_OUT(20,2,'(9) Fast Fourier Transformation (no graphics)');
    LINE_OUT(10,4,'1  Data file:');
    LINE_OUT(10,6,'2  Study area file:');
    LINE_OUT(10,8,'3  ');
    LINE_OUT(10,10,'4  Text output displayed on:');
    LINE_OUT(10,12,'5  ');
    LINE_OUT(10,14,'6  Restart file:');
    LINE_OUT(10,16,'7  ');
    LINE_OUT(10,18,'8  ');

    TEXTCOLOR(YELLOW);
    LINE_OUT(24,4,FILEOUT);
    LINE_OUT(30,6,STUDYFILE);
    LINE_OUT(27,14,UDRESTART);

    CASE UPCASE(PRINTCHAR) OF
      'D' : PRINTSTR := 'DISPLAY ONLY';
      'P' : PRINTSTR := 'PRINTER ONLY';
      'B' : PRINTSTR := 'BOTH DISPLAY AND PRINTER';
    END;
    LINE_OUT(39,10,PRINTSTR);

    DISPLAY_COLOR(RED);

    LINE_OUT(7,21,'P  Process Data');
    LINE_OUT(37,21,'E  Exit to Main Menu');
    TEXTCOLOR(YELLOW);
    GOTOXY(1,24);

    KEYVALUE := ' ';

    WHILE NOT (UPCASE(KEYVALUE) IN ['E','R']) DO
      BEGIN

        KEYVALUE := ' ';
        WHILE NOT (UPCASE(KEYVALUE) IN ['1','2','4','6','E','P']) DO
          READ (KBD,KEYVALUE);
        KEYVALUE := UPCASE(KEYVALUE);

        GOTOXY(1,23);
        DELLINE;
        DELLINE;

        CASE KEYVALUE OF
          '1': MENU_CALL(1);
          '2': MENU_CALL(2);
          '4': MENU_CALL(4);
          '6': BEGIN
               LINE_OUT(1,23,'Enter the name of the FFT RESTART file ');
               READLN  (INPUT,UDRESTART);
               GOTOXY(1,23);
               DELLINE;
               CAPIT(UDRESTART);
               LINE_OUT(27,14,UDRESTART);
               END;
          'P': BEGIN
               KEEP_GOING := TRUE;

               IF FILEOUT = '' THEN
                 BEGIN
                   LINE_OUT(1,23,'The input data set must be defined.');
                   KEEP_GOING := FALSE;
                 END;

               IF (KEEP_GOING) AND (NOT EXIST(FILEOUT)) THEN
                 BEGIN
                   GOTOXY(1,23);
                   WRITE('Formatted data file ',FILEOUT,' does not exist.');
                   KEEP_GOING := FALSE;
                 END;

               IF KEEP_GOING THEN
                 BEGIN
                   VERIFY_DATA_FILE;
                   IF NOT KEEP_GOING THEN
                     BEGIN
                       GOTOXY(1,23);
                       WRITE('The input file ',FILEOUT,' is not a formatted data file.');
                       IF (XLOW = -5) AND (YLOW = -3) AND (ZLOW = -7)
                         THEN LINE_OUT(1,24,'It is a formatted STUDY AREA file.')
                         ELSE LINE_OUT(1,24,'It may be an unformatted data file or a version 1.0 file.');
                     END;
                 END;

               IF (KEEP_GOING) AND (STUDYFILE <> '') THEN
                 IF NOT EXIST(STUDYFILE) THEN
                   BEGIN
                     GOTOXY(1,23);
                     WRITE ('Formatted study area file ',STUDYFILE,' does not exist.');
                     KEEP_GOING := FALSE;
                   END;

               IF (KEEP_GOING) AND (STUDYFILE <> '') THEN
                 BEGIN
                   VERIFY_STUDY_AREA_FILE;
                   IF NOT KEEP_GOING THEN
                     BEGIN
                       GOTOXY(1,23);
                       WRITE ('The study area file ',STUDYFILE,'  is not a formatted study area file.');
                       IF (SXLOW = -1) AND (SYLOW = -2) AND (SXHIGH = -3)
                         THEN LINE_OUT(1,24,'It is a formatted INPUT DATA file.')
                         ELSE LINE_OUT(1,24,'It may be an unformatted data file or unformatted study area file.');
                       END;
                     END;

               IF KEEP_GOING THEN
                 BEGIN
                   CLRSCR;
                   FFT_PROCEDURE;
                 END;

               END;  { CASE 'P' }

        END;  { CASE KEYVALUE OF }

        GOTOXY(1,25);

      END;  { WHILE NOT KEYVALUE IN [E,R] }

    END;  { WHILE KEYVALUE = 'R' }

    KEYVALUE := ' ';

  END;  { FFT_MENU PROCEDURE }