Personal tools
You are here: Home Projects ALGOL Source Code Hemker et al. Multigrid libraries. P. W. Hemker. A library of Euler multigrid routines.
Document Actions

P. W. Hemker. A library of Euler multigrid routines.

by Paul McJones last modified 2010-07-17 15:19

P. W. Hemker. A library of Euler multigrid routines. Centrum voor Wiskunde en Informatica, Numerieka Wiskunde, 1986.

Click here to get the file

Size 167.4 kB - File type text/plain

File contents

 #********************************************************************#
 #                                                                    #
 #                                                     PWH/860414     #
 #                                                                    #
 #                                                                    #
 #           A LIBRARY  OF EULER  MULTIGRID ROUTINES                  #
 #                                                                    #
 #           A PROGRAMME TO STUDY THE NUMERICAL BEHAVIOUR             #                                            
 #                                                                    #
 #           OF A MULTIGRID TECHNIQUE FOR THE SOLUTION OF             #
 #           THE EULER EQUATIONS FOR COMPRESSIBLE FLUID FLOW.         #
 #                                                                    #
 #           P.W. Hemker                                              #
 #                                                                    #
 #                                                                    #
 #           REFERENCES:                                              #
 #                                                                    #
 #           -- Hemker, P.W. & Spekreijse, S.P.,                      #
 #           Multigrid solution of the steady Euler equations. In:    #
 #           Advances in Multi-Grid methods, (D.Braess, W.Hackbusch & #
 #           U.Trottenberg eds) pp. 33-44, Proceedings of the         #
 #           conference held in Oberwolfach, Dec. 1984, Vieweg Publ.  #
 #           Comp., Braunschweig, 1985                                #
 #                                                                    #
 #           -- Hemker, P.W.,                                         #
 #           Multigrid and improved accuracy for the Euler equations. #
 #           In: Discretization in differential equations and         #
 #           enclosures, (E.Adams, R.Ansorge, Ch.Groszmann and        #
 #           H.-G.Roos eds) pp. 123-136, Akademie-Verlag Berlin, 1987 #
 #                                                                    #
 #           -- Hemker, P.W. & Spekreijse, S.,                        #
 #           Multiple Grid and Osher's Scheme for the Efficient       #
 #           Solution of the Steady Euler Equations. J. Appl. Num.    #
 #           Math. 2 (1986) 475-493.                                  #
 #                                                                    #
 #           -- Hemker, P.W.,                                         #
 #           Defect correction and higher order schemes for the       #
 #           multigrid solution of the steady Euler equations. In:    #
 #           Multi-Grid Methods II (W. Hackbusch and U. Trottenberg   #
 #           eds) pp. 149-165, Lecture Notes in Mathematics 1228,     #
 #           Springer Verlag, 1986.                                   #
 #                                                                    #          
 #           --  Hemker, P.W. & Johnson, G.M.,                        #
 #           Multigrid approaches to the Euler equations. Chapter 3 in:
 #           Multigrid Methods, (S.F. McCormick ed.) pp. 57-72, SIAM  #
 #           Frontiers in Appl. Math. Vol.3, 1987.                    #
 #                                                                    #
 #                                                                    #
 #                                                                    #
 #           HISTORY:                                                 #
 #                                                                    #           
 #           PF=EUMG001                   VSN.840222                  #
 #           PF=EUMG002                   VSN.840306                  #
 #           PF=EUMG003                   VSN.840308                  #
 #           PF=EUMG004                   VSN.840316                  #
 #           PF=EUMG005 (BSTATE)          VSN.840327                  #
 #           PF=EUMG006 (2ND ORDER)       VSN.840517                  #
 #           PF=EUMG007 (SEIDEL)(S)       VSN.840604                  #
 #           PF=EUMG008 (ILLU   )         VSN.840608                  #
 #           PF=EUMG009 (LIBRARY)         VSN.840613                  #
 #           PF=EUMG010 (STORAGE)         VSN.840615                  #
 #           PF=EUMG011 (U,V,C,Z)         VSN.840702                  #
 #           PF=EUMG012 (LINSYS)(FTN,P)   VSN.840709                  #
 #           PF=EUMG013 (MAPPING)         VSN.840724                  #
 #           PF=EUMG014 (NEW FIELDS)(S)   VSN.840823                  #
 #                      (WALLS,FLUXT)                                 #
 #           PF=EUMG015 (OSH-NAT)         VSN.840927                  #
 #                      (MGCS)                                        #
 #           PF=EUMG016 (RB4/2 RELAX)     VSN.841114                  #
 #                      (FMG,PROLON2)                                 #
 #           PF=EUMG017 (RAP)(P)          VSN.850102                  #
 #                      (NEWTON)(P)                                   #
 #                      (SEIDEL:LEX/RB)                               #
 #           PF=EUMG018 (NEWTON IT)       VSN.850322                  #
 #                      (RESIDU2ND)                                   #
 #           PF=EUMG019 (ST.PROBLEMS)     VSN.850515                  #
 #                      (SMOOTH)                                      #
 #           PF=EUMG020 (VL 2ND ORDER)    VSN.850731                  #
 #                      (JACOBI)                                      #
 #           PF=EUMG021 (STANDARD PROBLS) VSN.850930                  #
 #           PF=EUMG022 (STRUCTS)  (JR)   VSN.860110                  #
 #                      (ST/SF EXTENDED)(JR)                          #
 #           PF=EUMG023 (TOTVAR/SETFIELD) VSN.860224                  #
 #                      (PROBLEM 3B)                                  #
 #                      (STAR ROUTINES )                              #
 #           PF=EUMG024 (MDCP)            VSN.860414                  #
 #                      (2ND INTERP)(P)                               #
 #                                                                    #
 #********************************************************************#
 
 EULIB:                                               'PR' EJECT 'PR'
 'BEGIN' # TWO-DIMENSIONAL NON-ISENTROPIC EULERIAN FLOW #
 
   'PRIO'  'MINRESIDU' = 1, 'PLUSRESIDU' = 1,  'CORRECT' = 1,
           'PRINT'     = 1,          -*  = 4;
 
                                                  # TYPE            #
   'MODE' 'FLUX'  = 'STRUCT'('REAL' C1,C2,C3,C4); # 0 :             #
   'MODE' 'STATE' = 'STRUCT'('REAL' C1,C2,C3,C4); # 1 :(R,RU,RV,RE) #
                                                  # 2 :(U, V, C, Z) #
   'MODE' 'FIELD' = 'STRUCT'('REAL' DDX,DDY,
                             'REF''INT' TYPE,   'REF'[,]'STATE' Q);
   'MODE' 'FLUXAC'= 'STRUCT'('REAL'C11,C12,C13,C14,C21,C22,C23,C24,
                                   C31,C32,C33,C34,C41,C42,C43,C44);
   'MODE' 'STATAC'= 'STRUCT'('REAL'C11,C12,C13,C14,C21,C22,C23,C24,
                                   C31,C32,C33,C34,C41,C42,C43,C44);
 
   'MODE' 'LINMAT'= 'STRUCT'('REF''REF'[,]'FLUXAC' ARRAY,
                             'INT' JL,JU,KL,KU,P);
   'MODE' 'NETMAT'= 'REF'[]'LINMAT';
 
   'MODE' 'POINT' = 'STRUCT'('REAL'  X, Y);
   'MODE' 'WALL'  = 'STRUCT'('REAL'  C, S,DS);
 
 #     .-------------------------------->  YY                    #
 #     !                          ----->                         #
 #     !                            J                            #
 #     !  !       JL . . . . . JU                                #
 #     !  !                                        NW      NO    #
 #     !  !   IL   * * * * * * *                                 #
 #     !  !    .   * * * * * * *       MAPPING                   #
 #     !  !    .   * * * * * * *      ---------->    (X,Y)       #
 #     !  !    .   * * * * * * *                                 #
 #     !  !   IU   * * * * * * *                                 #
 #     !  !                                        ZW      ZO    #
 #  XX !  V                                                      #
 #     V    I                                                    #
 #                                                               #
 #         COMPUTATIONAL DOMAIN                  PHYSICAL DOMAIN #
 
 
   #$H#
   'PROC'  COMP FIELD = ('REAL' DDX,DDY, 'INT' XXL,XXU,YYL,YYU)'FIELD':
   #$B CMP.FLD B$#
   'BEGIN' XXN:= XXL*DDX; XXZ:= XXU*DDX; YYW:= YYL*DDY; YYO:= YYU*DDY;
           (DDX,DDY,'HEAP''INT',
            'HEAP'[(XXL+1):XXU,(YYL+1):YYU]'STATE')
   'END'; #$E#
          #E$#
 
   'REAL'  XXN,XXZ,YYW,YYO;  # GRENZEN VAN HET REKENGEBIED #
   'PROC'  MAPPING := ('REAL' XX,YY)'POINT': (XX,YY);
 
 # BACK-STORAGE SYSTEM            #                   'PR' EJECT 'PR'
 
   'PROC' BACKSTORE = 'VOID': # PREPARES THE BACK-FILE #
                'PR' XREF A68FTN,START'PR'  'SKIP';
   'PROC' PPUT = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID':
                'PR' XREF A68FTN,PUT  'PR'  'SKIP';
   'PROC' GGET = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID':
                'PR' XREF A68FTN,GET  'PR'  'SKIP';
   'PROC' GGEN = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID':
                'PR' XREF A68FTN,GEN  'PR'  'SKIP';
 
   'INT' LONGLINE:= 5;
 
   #$H#
   'PROC' PET = ('BOOL' CHANGE, 'LINMAT' L)'VOID':
   #$B PET.68  B$#
   'IF'   'INT'P:= P'OF'L; P>0
   'THEN' 'IF'   CHANGE
          'THEN' 'INT' JL = JL'OF'L, KL = KL'OF'L;
                 'INT'S:=(KU'OF'L-KL+1)*(JU'OF'L-JL+1)*16;
                 'REF'[,]'FLUXAC' ARRAY = ARRAY'OF'L;
                 PPUT(C11'OF'ARRAY[JL,KL],S,P)
          'FI' ;
          ARRAY'OF'L:= 'REF'[,]'FLUXAC'('NIL')
   'FI' ; #$E#
          #E$#
 
   #$H#
   'PROC' GET = ('LINMAT' L)'REF'[,]'FLUXAC':
   #$B GET.68  B$#
   'IF'   'REF'[,]'FLUXAC'(ARRAY'OF'L) 'ISNT' 'REF'[,]'FLUXAC'('NIL')
   'THEN' ARRAY'OF'L
   'ELSE' 'INT' JL=JL'OF'L, JU=JU'OF'L, KL=KL'OF'L, KU=KU'OF'L;
          'INT' P:= P'OF'L,
                S:= (KU-KL+1)*(JU-JL+1)*16;
          'REF'[,]'FLUXAC' ARRAY = 'HEAP'[JL:JU,KL:KU]'FLUXAC';
          ARRAY'OF'L:= ARRAY;
          GGET(C11'OF'ARRAY[JL,KL],S,P); ARRAY
   'FI' ; #$E#
          #E$#
 
   #$H#
   'PROC' GEN = ('REF''LINMAT' L, 'INT' JL,JU,KL,KU)'REF'[,]'FLUXAC':
   #$B GEN.68  B$#
   'BEGIN'
    'INT' S,P:= 0, 'INT' LJ = JU-JL+1,
    'REF'[,]'FLUXAC' ARRAY = 'HEAP'[JL:JU,KL:KU]'FLUXAC';
    'IF'   LJ>LONGLINE
    'THEN' S:= (KU-KL+1)*LJ*16;
           GGEN(C11'OF'ARRAY[JL,KL],S,P)
    'FI'   ;
    'REF''REF'[,]'FLUXAC' ARR = 'HEAP''REF'[,]'FLUXAC':= ARRAY;
    L:= (ARR,JL,JU,KL,KU,P); ARRAY
   'END' ; #$E#
           #E$#
 
 
 #$H#
 'PROC' ATTACH = ('STRING' LFN, 'INT' NR FIELDS )'VOID':
 #$B ATTACH B$#
 'BEGIN' OUTCHECK:= 0; OUTTOP:= NR FIELDS; WRITE:= NR FIELDS=0;
         ESTABLISH(FILE,LFN,BINARYSEQUENTIALCHANNEL,1,5000,3500)
 'END'; #$E#
        #E$#
 
 'FILE' FILE; 'BOOL' WRITE; 'INT' OUTCHECK, OUTTOP;
 
  #$H#
  'PROC' OUT = ('INT' NO, 'FIELD' F)'VOID':
  #$B FLD.OUT B$#
  'BEGIN' PRINT((NEWLINE, "OUT: " ,WHOLE(NO,0) ));
          'REF' [,]'STATE' Q = Q'OF'F;
          'INT' L1=1'LWB'Q, U1=1'UPB'Q, L2=2'LWB'Q, U2=2'UPB'Q,
                TYPE = TYPE'OF'F;
          'REAL' DDX = DDX'OF'F, DDY = DDY'OF'F;
          PUTBIN(FILE, (NO,L1,U1,L2,U2,TYPE,DDX,DDY,NEWLINE) );
          PUTBIN(FILE, (Q, 221141,NEWLINE) )
  'END'; #$E#
         #E$#
 
  #$H#
  'PROC' IN = ('INT' NO, 'REF''FIELD' F)'VOID':
  #$B FLD.IN B$#
  'BEGIN' PRINT((" IN : " ,WHOLE(NO,0) ));
          'INT' L1,U1,L2,U2,TYPE,NOC;
          'REAL' DDX, DDY;
          GETBIN(FILE, (NOC,L1,U1,L2,U2,TYPE,DDX,DDY,NEWLINE) );
          ( NOC /= NO ! PRINT((NEWLINE,"IN-FOUT: ", NO, NOC)) );
          'HEAP' [L1:U1,L2:U2] 'STATE' Q;
          GETBIN(FILE, (Q, NOC,NEWLINE) );
          ( NOC /= 221141 ! PRINT((NEWLINE,"IN-FOUT SLOT:", NO)) );
          F:= 'FIELD'(DDX,DDY,'HEAP''INT':=TYPE,Q)
  'END'; #$E#
         #E$#
 
  #$H#
  'PROC' REWIND = 'VOID':
  #$B REWIND B$#
  'BEGIN' ( WRITE ! OUTTOP:= OUTCHECK; PUTBIN(FILE, (NEWLINE)) );
          WRITE:='FALSE'; RESET(FILE); OUTCHECK:=0
  'END'; #$E#
         #E$#
 
  #$H#
  'PROC' SKIPF = ('INT' N)'VOID':
  #$B SKIPF B$#
  'BEGIN' 'FIELD' SF;
          'IF' OUTCHECK+N <= OUTTOP
          'THEN' 'TO' N 'DO' IN(OUTCHECK+:=1,SF) 'OD'
          'FI'
  'END' ; #$E#
          #E$#
 
  #$H#
  'PROC' IN FRONT = ('INT' N)'VOID':
  #$B IN.FRNT B$#
  'BEGIN' PRINT((NEWLINE, "FRONT : " ,WHOLE(N,0) ));
          ( N > OUTTOP ! ERROR );
          ( OUTCHECK>=N ! REWIND );  SKIPF(N-1-OUTCHECK)
  'END'; #$E#
         #E$#
 
  #$H#
  'PROC' TAKE = ('INT' N, 'REF''FIELD' SF)'VOID':
  #$B FL.TAKE B$#
  'BEGIN' IN FRONT (N);
          IN(OUTCHECK+:=1,SF);
          PRINT((" TAKEN : " ,WHOLE(OUTCHECK,0) ))
  'END'; #$E#
         #E$#
 
 # PICTURE PREPARATION            #                   'PR' EJECT 'PR'
 
 'FILE' TRANSPORT;
 
 'PROC' TRANSA = 'VOID':
        ESTABLISH(TRANSPORT,"TAPE10",ZTYPE CHANNEL,1,100000,40);
 
   #$H#
 'PROC' TRANSB = 'VOID':
   #$B TRANSB  B$#
 'BEGIN' 'REAL' ENDOFFILE = 77777;
         PUT ( TRANSPORT,(FLOAT(ENDOFFILE,12,6,2),NEWLINE) ) ;
         CLOSE(T RANSPORT)
 'END' ; #$E#
         #E$#
 # INFOPL    MEANING                       #
 #        IF X-DIM AND Y-DIM > 1           #
 #   +/-3    3D-PLOT                       #
 #   +/-2    CONTOUR-PLOT                  #
 #        IF X-DIM OR  Y-DIM = 1           #
 #     -1    CURVE  -PLOT, NO LEGEND       #
 #      0    CURVE  -PLOT, 1ST ORDER       #
 #      1    CURVE  -PLOT, 1  DCP          #
 #      2    CURVE  -PLOT, 2  DCP          #
 #      3    CURVE  -PLOT, 3  DCP          #
 #      4    CURVE  -PLOT, 4  DCP          #
 #      5    CURVE  -PLOT, EXACT           #
 
   #$H#
  'PROC' TRANSV = ('INT' INFOPL, 'VAR' VARIA) 'VOID':
   #$B TRANSV  B$#
  ([1:20]'REAL' DUMP,
   'REF'[,]'REAL' Q = Q'OF'VARIA,
   'INT' INT = (INTERN'OF'VARIA!0!1);
   'REAL' R, 'INT' KEY,
   'INT' L1 = 1'LWB'Q, U1 = 1'UPB'Q, L2 = 2'LWB'Q, U2 = 2'UPB'Q;
   'FOR' K 'FROM' 1 'TO' 20 'DO' DUMP[K]:= 0 'OD';
   'IF'   (U1-L1+1)>1 'AND' (U2-L2+1)>1
   'THEN' KEY:= 54321;
          DUMP[1:12]:= (KEY,(U1-L1+1),(U2-L2+1),L1,U1,L2,U2,
                      DDX'OF'VARIA,DDY'OF'VARIA,INFOPL,VAR'OF'VARIA,INT)
   'ELSE' KEY:= 76543;
          DUMP[1:8]:=
          ((U1-L1+1)>1 !(KEY,(U1-L1+1),L1,U1,DDX'OF'VARIA,INFOPL,
                         VAR'OF'VARIA,INT)
                       !(KEY,(U2-L2+1),L2,U2,DDY'OF'VARIA,INFOPL,
                         VAR'OF'VARIA,INT) )
   'FI' ;
   'FOR' K 'TO' (KEY=54321!20!10)
   'DO' ( 'ABS'(R:= DUMP[K]) >= 1.0E10 ! ERROR
                      ! PUT(TRANSPORT,(FLOAT(R,12,6,2),NEWLINE)) )
   'OD';
   'FOR' I 'FROM' L1 'TO' U1
   'DO'  'FOR' J 'FROM' L2 'TO' U2
         'DO' R:= Q[I,J];
              ('ABS'R <= 1.0E-10 ! R:=0 );
              ('ABS'R >= 1.0E+10 ! ERROR);
              PUT(TRANSPORT,(FLOAT(R,12,6,2),NEWLINE))
         'OD'
   'OD'); #$E#
          #E$#
 
 # OPERATORS                      #                   'PR' EJECT 'PR'
 
   'OP' ** = ('REAL' A,B)'REAL': EXP( B*LN(A) );
 
   'OP' - = ('FLUX' F)'FLUX':(-C1'OF'F,-C2'OF'F,-C3'OF'F,-C4'OF'F);
 
   'OP' + = ('FLUX' A,B)'FLUX':
   (C1'OF'A+C1'OF'B,C2'OF'A+C2'OF'B,C3'OF'A+C3'OF'B,C4'OF'A+C4'OF'B);
 
   'OP' - = ('FLUX' A,B)'FLUX':
   (C1'OF'A-C1'OF'B,C2'OF'A-C2'OF'B,C3'OF'A-C3'OF'B,C4'OF'A-C4'OF'B);
 
   'OP' * = ('FLUX' A,B)'REAL':
    C1'OF'A*C1'OF'B+C2'OF'A*C2'OF'B+C3'OF'A*C3'OF'B+C4'OF'A*C4'OF'B ;
 
   'OP' * = ('REAL' B, 'FLUX' A)'FLUX':
   ( B*C1'OF'A , B*C2'OF'A , B*C3'OF'A , B*C4'OF'A );
 
   'OP' / = ('FLUX' A, 'REAL' B)'FLUX':
   ( C1'OF'A/B , C2'OF'A/B , C3'OF'A/B , C4'OF'A/B ) ;
 
   'OP' +:= = ('REF''FLUX' A, 'FLUX' B)'REF''FLUX':
        ((C1'OF'A +:= C1'OF'B , C2'OF'A +:= C2'OF'B ,
          C3'OF'A +:= C3'OF'B , C4'OF'A +:= C4'OF'B );A ) ;
 
   'OP' -:= = ('REF''FLUX' A, 'FLUX' B)'REF''FLUX':
        ((C1'OF'A -:= C1'OF'B , C2'OF'A -:= C2'OF'B ,
          C3'OF'A -:= C3'OF'B , C4'OF'A -:= C4'OF'B );A ) ;
 
   'OP' *:= = ('REF''FLUX' A, 'REAL' B)'REF''FLUX':
    ((C1'OF'A*:=B , C2'OF'A*:=B , C3'OF'A*:=B , C4'OF'A*:=B);A);
 
   'OP' /:= = ('REF''FLUX' A, 'REAL' B)'REF''FLUX':
    ((C1'OF'A/:=B , C2'OF'A/:=B , C3'OF'A/:=B , C4'OF'A/:=B);A);
 
   'OP' 'ABS' = ('FLUX' A)'REAL':
     'ABS'C1'OF'A+ 'ABS'C2'OF'A+ 'ABS'C3'OF'A+ 'ABS'C4'OF'A ;
 
   #$H#
   'OP' - = ('FLUXAC' A)'FLUXAC' :
   #$B MINFA   B$#
   'FLUXAC'( -C11'OF'A,-C12'OF'A,-C13'OF'A,-C14'OF'A ,
             -C21'OF'A,-C22'OF'A,-C23'OF'A,-C24'OF'A ,
             -C31'OF'A,-C32'OF'A,-C33'OF'A,-C34'OF'A ,
             -C41'OF'A,-C42'OF'A,-C43'OF'A,-C44'OF'A
           ) ; #$E#
               #E$#
 
   #$H#
   'OP' + = ('FLUXAC' A,B)'FLUXAC':
   #$B FAPLFA  B$#
     ( C11'OF'A + C11'OF'B , C12'OF'A + C12'OF'B ,
       C13'OF'A + C13'OF'B , C14'OF'A + C14'OF'B ,
       C21'OF'A + C21'OF'B , C22'OF'A + C22'OF'B ,
       C23'OF'A + C23'OF'B , C24'OF'A + C24'OF'B ,
       C31'OF'A + C31'OF'B , C32'OF'A + C32'OF'B ,
       C33'OF'A + C33'OF'B , C34'OF'A + C34'OF'B ,
       C41'OF'A + C41'OF'B , C42'OF'A + C42'OF'B ,
       C43'OF'A + C43'OF'B , C44'OF'A + C44'OF'B
      ) ; #$E#
          #E$#
 
   #$H#
   'OP' - = ('FLUXAC' A,B)'FLUXAC':
   #$B FAMNFA  B$#
     ( C11'OF'A - C11'OF'B , C12'OF'A - C12'OF'B ,
       C13'OF'A - C13'OF'B , C14'OF'A - C14'OF'B ,
       C21'OF'A - C21'OF'B , C22'OF'A - C22'OF'B ,
       C23'OF'A - C23'OF'B , C24'OF'A - C24'OF'B ,
       C31'OF'A - C31'OF'B , C32'OF'A - C32'OF'B ,
       C33'OF'A - C33'OF'B , C34'OF'A - C34'OF'B ,
       C41'OF'A - C41'OF'B , C42'OF'A - C42'OF'B ,
       C43'OF'A - C43'OF'B , C44'OF'A - C44'OF'B
     ) ; #$E#
         #E$#
 
   #$H#
   'OP' * = ('REAL' B, 'FLUXAC' A)'FLUXAC':
   #$B RTIMFA  B$#
     ( B*C11'OF'A,B*C12'OF'A,B*C13'OF'A,B*C14'OF'A ,
       B*C21'OF'A,B*C22'OF'A,B*C23'OF'A,B*C24'OF'A ,
       B*C31'OF'A,B*C32'OF'A,B*C33'OF'A,B*C34'OF'A ,
       B*C41'OF'A,B*C42'OF'A,B*C43'OF'A,B*C44'OF'A
     ) ; #$E#
         #E$#
 
   #$H#
   'OP' / = ('FLUXAC' A,'REAL' B )'FLUXAC':
   #$B FADIVR  B$#
     ( C11'OF'A/B,C12'OF'A/B,C13'OF'A/B,C14'OF'A/B ,
       C21'OF'A/B,C22'OF'A/B,C23'OF'A/B,C24'OF'A/B ,
       C31'OF'A/B,C32'OF'A/B,C33'OF'A/B,C34'OF'A/B ,
       C41'OF'A/B,C42'OF'A/B,C43'OF'A/B,C44'OF'A/B
     ) ; #$E#
         #E$#
 
   #$H#
   'OP' +:= = ('REF''FLUXAC' A, 'FLUXAC' B)'REF''FLUXAC':
   #$B FAPBFA  B$#
   'BEGIN' (C11'OF'A +:= C11'OF'B , C12'OF'A +:= C12'OF'B ,
            C13'OF'A +:= C13'OF'B , C14'OF'A +:= C14'OF'B ,
            C21'OF'A +:= C21'OF'B , C22'OF'A +:= C22'OF'B ,
            C23'OF'A +:= C23'OF'B , C24'OF'A +:= C24'OF'B ,
            C31'OF'A +:= C31'OF'B , C32'OF'A +:= C32'OF'B ,
            C33'OF'A +:= C33'OF'B , C34'OF'A +:= C34'OF'B ,
            C41'OF'A +:= C41'OF'B , C42'OF'A +:= C42'OF'B ,
            C43'OF'A +:= C43'OF'B , C44'OF'A +:= C44'OF'B );A
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' -:= = ('REF''FLUXAC' A, 'FLUXAC' B)'REF''FLUXAC':
   #$B FAMBFA  B$#
   'BEGIN' (C11'OF'A -:= C11'OF'B , C12'OF'A -:= C12'OF'B ,
            C13'OF'A -:= C13'OF'B , C14'OF'A -:= C14'OF'B ,
            C21'OF'A -:= C21'OF'B , C22'OF'A -:= C22'OF'B ,
            C23'OF'A -:= C23'OF'B , C24'OF'A -:= C24'OF'B ,
            C31'OF'A -:= C31'OF'B , C32'OF'A -:= C32'OF'B ,
            C33'OF'A -:= C33'OF'B , C34'OF'A -:= C34'OF'B ,
            C41'OF'A -:= C41'OF'B , C42'OF'A -:= C42'OF'B ,
            C43'OF'A -:= C43'OF'B , C44'OF'A -:= C44'OF'B );A
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' *:= = ('REF''FLUXAC' A, 'REAL' B)'REF''FLUXAC':
   #$B FATBRL  B$#
   'BEGIN'
     (C11'OF'A*:=B,C12'OF'A*:=B,C13'OF'A*:=B,C14'OF'A*:=B,
      C21'OF'A*:=B,C22'OF'A*:=B,C23'OF'A*:=B,C24'OF'A*:=B,
      C31'OF'A*:=B,C32'OF'A*:=B,C33'OF'A*:=B,C34'OF'A*:=B,
      C41'OF'A*:=B,C42'OF'A*:=B,C43'OF'A*:=B,C44'OF'A*:=B);A
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' /:= = ('REF''FLUXAC' A, 'REAL' B)'REF''FLUXAC':
   #$B FADBRL  B$#
   'BEGIN'
     (C11'OF'A/:=B,C12'OF'A/:=B,C13'OF'A/:=B,C14'OF'A/:=B,
      C21'OF'A/:=B,C22'OF'A/:=B,C23'OF'A/:=B,C24'OF'A/:=B,
      C31'OF'A/:=B,C32'OF'A/:=B,C33'OF'A/:=B,C34'OF'A/:=B,
      C41'OF'A/:=B,C42'OF'A/:=B,C43'OF'A/:=B,C44'OF'A/:=B);A
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' * = ('FLUXAC' A, 'FLUX' B)'FLUX':
   #$B FATIMF  B$#
   'BEGIN'
     'REAL' B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B ;
     'FLUX'(C11'OF'A*B1+C12'OF'A*B2+C13'OF'A*B3+C14'OF'A*B4,
            C21'OF'A*B1+C22'OF'A*B2+C23'OF'A*B3+C24'OF'A*B4,
            C31'OF'A*B1+C32'OF'A*B2+C33'OF'A*B3+C34'OF'A*B4,
            C41'OF'A*B1+C42'OF'A*B2+C43'OF'A*B3+C44'OF'A*B4)
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' * = ('FLUX' B, 'FLUXAC' A)'FLUX':
   #$B FTIMFA  B$#
   'BEGIN'
     'REAL' B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B ;
     'FLUX'(B1*C11'OF'A+B2*C21'OF'A+B3*C31'OF'A+B4*C41'OF'A,
            B1*C12'OF'A+B2*C22'OF'A+B3*C32'OF'A+B4*C42'OF'A,
            B1*C13'OF'A+B2*C23'OF'A+B3*C33'OF'A+B4*C43'OF'A,
            B1*C14'OF'A+B2*C24'OF'A+B3*C34'OF'A+B4*C44'OF'A)
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' * = ('FLUXAC' A,B)'FLUXAC' :
   #$B FATIFA  B$#
   'FLUXAC' (
C11'OF'A*C11'OF'B+C12'OF'A*C21'OF'B+C13'OF'A*C31'OF'B+C14'OF'A*C41'OF'B,
C11'OF'A*C12'OF'B+C12'OF'A*C22'OF'B+C13'OF'A*C32'OF'B+C14'OF'A*C42'OF'B,
C11'OF'A*C13'OF'B+C12'OF'A*C23'OF'B+C13'OF'A*C33'OF'B+C14'OF'A*C43'OF'B,
C11'OF'A*C14'OF'B+C12'OF'A*C24'OF'B+C13'OF'A*C34'OF'B+C14'OF'A*C44'OF'B,
C21'OF'A*C11'OF'B+C22'OF'A*C21'OF'B+C23'OF'A*C31'OF'B+C24'OF'A*C41'OF'B,
C21'OF'A*C12'OF'B+C22'OF'A*C22'OF'B+C23'OF'A*C32'OF'B+C24'OF'A*C42'OF'B,
C21'OF'A*C13'OF'B+C22'OF'A*C23'OF'B+C23'OF'A*C33'OF'B+C24'OF'A*C43'OF'B,
C21'OF'A*C14'OF'B+C22'OF'A*C24'OF'B+C23'OF'A*C34'OF'B+C24'OF'A*C44'OF'B,
C31'OF'A*C11'OF'B+C32'OF'A*C21'OF'B+C33'OF'A*C31'OF'B+C34'OF'A*C41'OF'B,
C31'OF'A*C12'OF'B+C32'OF'A*C22'OF'B+C33'OF'A*C32'OF'B+C34'OF'A*C42'OF'B,
C31'OF'A*C13'OF'B+C32'OF'A*C23'OF'B+C33'OF'A*C33'OF'B+C34'OF'A*C43'OF'B,
C31'OF'A*C14'OF'B+C32'OF'A*C24'OF'B+C33'OF'A*C34'OF'B+C34'OF'A*C44'OF'B,
C41'OF'A*C11'OF'B+C42'OF'A*C21'OF'B+C43'OF'A*C31'OF'B+C44'OF'A*C41'OF'B,
C41'OF'A*C12'OF'B+C42'OF'A*C22'OF'B+C43'OF'A*C32'OF'B+C44'OF'A*C42'OF'B,
C41'OF'A*C13'OF'B+C42'OF'A*C23'OF'B+C43'OF'A*C33'OF'B+C44'OF'A*C43'OF'B,
C41'OF'A*C14'OF'B+C42'OF'A*C24'OF'B+C43'OF'A*C34'OF'B+C44'OF'A*C44'OF'B
     ) ; #$E#
         #E$#
 
   #$H#
   'OP' -* = ('FLUXAC' A,B)'FLUXAC':
   #$B FAMTFA  B$#
   'BEGIN'         -(A*B)
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' /   = ('INT' ONE, 'FLUXAC' A)'FLUXAC':
   #$B GAUSIN  B$#
     'BEGIN'
     # 4*4 MATRIX INVERSION #
       (ONE /= 1 ! ERROR );
       'PROC' INV = ('REF''REAL' A, B, PIVOT)'VOID':
         'PR' XREF A68FTN,INVERT 'PR' 'SKIP';
       [1:4,1:4]'REAL' AA,B ;'REAL' PIVOT ;
       AA := ( (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A) ,
               (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A) ,
               (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A) ,
               (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A) );
       INV(AA[1,1],B[1,1],PIVOT);
       (PIVOT<1.0E-5 ! PRINT((NEWLINE," MAYBE SINGULAR MATRIX,",
                              " PIVOT = ",PIVOT)); (PIVOT=0!ERROR) );
       'FLUXAC' ( B[1,1],B[2,1],B[3,1],B[4,1],
                  B[1,2],B[2,2],B[3,2],B[4,2],
                  B[1,3],B[2,3],B[3,3],B[4,3],
                  B[1,4],B[2,4],B[3,4],B[4,4] )
     'END';#$E#
           #E$#
 
   #$H#
   'OP' / = ('FLUX' Y, 'FLUXAC' A)'STATE':
   #$B GAUSEL  B$#
    'BEGIN' 'PROC' GAUSEL = ('REF''REAL' C,X, PIVOT)'VOID':
               'PR' XREF A68FTN,GAUSS 'PR' 'SKIP';
       [1:4] 'REAL' X, [1:5,1:4] 'REAL' C, 'REAL' PIVOT;
       # A68 <-> FORTRAN TRANPOSITION OF ROWS AND COLUMNS. #
       C := ( (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A),
              (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A),
              (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A),
              (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A),
              (C1 'OF'Y,C2 'OF'Y,C3 'OF'Y,C4 'OF'Y) );
        GAUSEL(C[1,1],X[1],PIVOT);
        (PIVOT<1.0E-5 ! PRINT((NEWLINE," MAYBE SINGULAR MATRIX,",
                               " PIVOT = ",PIVOT)); (PIVOT=0!ERROR));
        'STATE' ( X[1],X[2],X[3],X[4] )
   'END'; #$E#
          #E$#
 
   #$H#
  'OP' / = ('REF'[]'FLUX' B, 'REF'[,]'FLUXAC' A)'REF'[]'FLUX':
   #$B DIV4    B$#
  'BEGIN' [1:4,1:4]'FLUXAC' MAT;
          'HEAP'[1:4]'FLUX' VEC;
          'PROC' SOL4 = ('REF'[,]'FLUXAC' A, 'REF'[]'FLUX' R)'VOID':
        # NB: OVERSCHRIJFT A EN R !!!!!!!!!!!!!!! #
        'BEGIN' 'FOR' I 'TO' 4
                'DO' A[I,I]:= 1/A[I,I];
                     'FOR' J 'FROM' I+1 'TO' 4
                     'DO'      A[I,J] := A[I,I]*A[I,J];
                          'FOR' K 'FROM' I+1 'TO' 4
                          'DO' A[K,J]-:= A[K,I]*A[I,J]
                     'OD' 'OD';
                               R[I  ] := A[I,I]*R[I  ];
                          'FOR' K 'FROM' I+1 'TO' 4
                          'DO' R[K  ]-:= A[K,I]*R[I  ]
                'OD'      'OD';
 
                'FOR' I 'FROM' 3 'BY' -1 'TO' 1
                'DO''FOR' J 'FROM' I+1 'TO' 4
                    'DO' R[I  ]-:= A[I,J]*R[J  ]
                'OD''OD'
          'END';
          'FOR' I 'TO' 4
          'DO' VEC[I]:= B[I];
               'FOR' J 'TO' 4
               'DO' MAT[I,J]:= A[I,J]
          'OD' 'OD';
          SOL4(MAT,VEC);
          VEC
   'END'; #$E#
          #E$#
 
 # COLUMN HANDLING PROCEDURES     #                   'PR' EJECT 'PR'
 
   'PRIO' 'FC'=9 ;
   #$H#
   'OP' 'FC' = ('INT' I,'FLUXAC'A)'FLUX':   # FREE A COLUMN #
   #$B FRE.COL B$#
   (I! 'FLUX' (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A),
       'FLUX' (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A),
       'FLUX' (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A),
       'FLUX' (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A)
     ! NUL ) ; #$E#
               #E$#
 
 # ASSIGNMENT     OF FLUX B TO COLUMN I OF A     #
   #$H#
   'PROC' CASS = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B COL.ASS B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A:=B1,C21'OF'A:=B2,C31'OF'A:=B3,C41'OF'A:=B4)),
       ((C12'OF'A:=B1,C22'OF'A:=B2,C32'OF'A:=B3,C42'OF'A:=B4)),
       ((C13'OF'A:=B1,C23'OF'A:=B2,C33'OF'A:=B3,C43'OF'A:=B4)),
       ((C14'OF'A:=B1,C24'OF'A:=B2,C34'OF'A:=B3,C44'OF'A:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -CASS- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # ADDITION       OF FLUX B TO COLUMN I OF A     #
   #$H#
   'PROC' CADD = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B COL.ADD B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A+:=B1,C21'OF'A+:=B2,C31'OF'A+:=B3,C41'OF'A+:=B4)),
       ((C12'OF'A+:=B1,C22'OF'A+:=B2,C32'OF'A+:=B3,C42'OF'A+:=B4)),
       ((C13'OF'A+:=B1,C23'OF'A+:=B2,C33'OF'A+:=B3,C43'OF'A+:=B4)),
       ((C14'OF'A+:=B1,C24'OF'A+:=B2,C34'OF'A+:=B3,C44'OF'A+:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -CADD- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # SUBSTRACTION   OF FLUX B TO COLUMN I OF A     #
   #$H#
   'PROC' CSUB = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B COL.SUB B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A-:=B1,C21'OF'A-:=B2,C31'OF'A-:=B3,C41'OF'A-:=B4)),
       ((C12'OF'A-:=B1,C22'OF'A-:=B2,C32'OF'A-:=B3,C42'OF'A-:=B4)),
       ((C13'OF'A-:=B1,C23'OF'A-:=B2,C33'OF'A-:=B3,C43'OF'A-:=B4)),
       ((C14'OF'A-:=B1,C24'OF'A-:=B2,C34'OF'A-:=B3,C44'OF'A-:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -CSUB- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # MULTIPLICATION OF COLUMN I OF A BY A FACTOR B #
   #$H#
   'PROC' CMUL = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC':
   #$B COL.MUL B$#
   ( 'IF' (A'ISNT''NIL')'AND'(B/=1.0) 'THEN'
    (I!((C11'OF'A*:=B,C21'OF'A*:=B,C31'OF'A*:=B,C41'OF'A*:=B)),
       ((C12'OF'A*:=B,C22'OF'A*:=B,C32'OF'A*:=B,C42'OF'A*:=B)),
       ((C13'OF'A*:=B,C23'OF'A*:=B,C33'OF'A*:=B,C43'OF'A*:=B)),
       ((C14'OF'A*:=B,C24'OF'A*:=B,C34'OF'A*:=B,C44'OF'A*:=B))
      ! PRINT((NEWLINE,"WRONG INDEX IN -CMUL- CALL"));ERROR )
     'FI';A) ; #$E#
               #E$#
 
 # DIVISION       OF COLUMN I OF A BY A FACTOR B #
   #$H#
   'PROC' CDIV = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC':
   #$B COL.DIV B$#
   ('IF' (A'ISNT''NIL')'AND'(B/=0.0)'AND'(B/=1.0) 'THEN'
     (I!((C11'OF'A/:=B,C21'OF'A/:=B,C31'OF'A/:=B,C41'OF'A/:=B)),
        ((C12'OF'A/:=B,C22'OF'A/:=B,C32'OF'A/:=B,C42'OF'A/:=B)),
        ((C13'OF'A/:=B,C23'OF'A/:=B,C33'OF'A/:=B,C43'OF'A/:=B)),
        ((C14'OF'A/:=B,C24'OF'A/:=B,C34'OF'A/:=B,C44'OF'A/:=B))
       ! PRINT((NEWLINE,"WRONG INDEX IN -CDIV- CALL"));ERROR )
     'FI';A) ; #$E#
               #E$#
 
 # ROW    HANDLING PROCEDURES     #                   'PR' EJECT 'PR'
 
   'PRIO' 'FR'=9 ;
   #$H#
   'OP' 'FR' = ('INT' I,'FLUXAC'A)'FLUX':   # FREE A ROW    #
   #$B FRE.ROW B$#
   (I! 'FLUX' (C11'OF'A,C12'OF'A,C13'OF'A,C14'OF'A),
       'FLUX' (C21'OF'A,C22'OF'A,C23'OF'A,C24'OF'A),
       'FLUX' (C31'OF'A,C32'OF'A,C33'OF'A,C34'OF'A),
       'FLUX' (C41'OF'A,C42'OF'A,C43'OF'A,C44'OF'A)
     ! NUL ) ; #$E#
               #E$#
 
 # ASSIGNMENT     OF FLUX B TO ROW    I OF A     #
   #$H#
   'PROC' RASS = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B ROW.ASS B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A:=B1,C12'OF'A:=B2,C13'OF'A:=B3,C14'OF'A:=B4)),
       ((C21'OF'A:=B1,C22'OF'A:=B2,C23'OF'A:=B3,C24'OF'A:=B4)),
       ((C31'OF'A:=B1,C32'OF'A:=B2,C33'OF'A:=B3,C34'OF'A:=B4)),
       ((C41'OF'A:=B1,C42'OF'A:=B2,C43'OF'A:=B3,C44'OF'A:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -RASS- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # ADDITION       OF FLUX B TO ROW    I OF A     #
   #$H#
   'PROC' RADD = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B ROW.ADD B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A+:=B1,C12'OF'A+:=B2,C13'OF'A+:=B3,C14'OF'A+:=B4)),
       ((C21'OF'A+:=B1,C22'OF'A+:=B2,C23'OF'A+:=B3,C24'OF'A+:=B4)),
       ((C31'OF'A+:=B1,C32'OF'A+:=B2,C33'OF'A+:=B3,C34'OF'A+:=B4)),
       ((C41'OF'A+:=B1,C42'OF'A+:=B2,C43'OF'A+:=B3,C44'OF'A+:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -RADD- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # SUBSTRACTION   OF FLUX B TO ROW    I OF A     #
   #$H#
   'PROC' RSUB = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC':
   #$B ROW.SUB B$#
   ('IF' A'ISNT''NIL'
    'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B;
    (I!((C11'OF'A-:=B1,C12'OF'A-:=B2,C13'OF'A-:=B3,C14'OF'A-:=B4)),
       ((C21'OF'A-:=B1,C22'OF'A-:=B2,C23'OF'A-:=B3,C24'OF'A-:=B4)),
       ((C31'OF'A-:=B1,C32'OF'A-:=B2,C33'OF'A-:=B3,C34'OF'A-:=B4)),
       ((C41'OF'A-:=B1,C42'OF'A-:=B2,C43'OF'A-:=B3,C44'OF'A-:=B4))
      ! PRINT((NEWLINE,"WRONG INDEX IN -RSUB- CALL"));ERROR )
    'FI';A) ; #$E#
              #E$#
 
 # MULTIPLICATION OF ROW    I OF A BY A FACTOR B #
   #$H#
   'PROC' RMUL = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC':
   #$B ROW.MUL B$#
   ( 'IF' (A'ISNT''NIL')'AND'(B/=1.0) 'THEN'
    (I!((C11'OF'A*:=B,C12'OF'A*:=B,C13'OF'A*:=B,C14'OF'A*:=B)),
       ((C21'OF'A*:=B,C22'OF'A*:=B,C23'OF'A*:=B,C24'OF'A*:=B)),
       ((C31'OF'A*:=B,C32'OF'A*:=B,C33'OF'A*:=B,C34'OF'A*:=B)),
       ((C41'OF'A*:=B,C42'OF'A*:=B,C43'OF'A*:=B,C44'OF'A*:=B))
      ! PRINT((NEWLINE,"WRONG INDEX IN -RMUL- CALL"));ERROR )
     'FI';A) ; #$E#
               #E$#
 
 # DIVISION       OF ROW    I OF A BY A FACTOR B #
   #$H#
   'PROC' RDIV = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC':
   #$B ROW.DIV B$#
   ( 'IF' (A'ISNT''NIL')'AND'(B/=0.0)'AND'(B/=1.0) 'THEN'
    (I!((C11'OF'A/:=B,C12'OF'A/:=B,C13'OF'A/:=B,C14'OF'A/:=B)),
       ((C21'OF'A/:=B,C22'OF'A/:=B,C23'OF'A/:=B,C24'OF'A/:=B)),
       ((C31'OF'A/:=B,C32'OF'A/:=B,C33'OF'A/:=B,C34'OF'A/:=B)),
       ((C41'OF'A/:=B,C42'OF'A/:=B,C43'OF'A/:=B,C44'OF'A/:=B))
      ! PRINT((NEWLINE,"WRONG INDEX IN -RDIV- CALL"));ERROR )
     'FI';A) ; #$E#
               #E$#
 
 # ILLU RELAX                     #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' ILLU RELAX = ('NETMAT' DEC,JAC, 'FIELD' QF,RF )'VOID':
   #$B ILLUR   B$#
   'BEGIN' 'REF'[,]'STATE'QQ = Q'OF'QF, RH=Q'OF'RF;
           'INT' L1= 1'LWB'QQ, U1= 1'UPB'QQ,
                 L2= 2'LWB'QQ, U2= 2'UPB'QQ;
           'LOC'[L2:U2]'STATE' DUI;
 
     'PROC' SOLL = ('REF'[]'FLUXAC' L,D,U, 'REF'[]'STATE' Z )'VOID':
     ('FOR'J'FROM' L2+1         'TO' U2 'DO' Z[J]+:= L[J]*Z[J-1] 'OD';
      'FOR'J'FROM' L2           'TO' U2 'DO' Z[J] := D[J]*Z[J  ] 'OD';
      'FOR'J'FROM' U2-1 'BY' -1 'TO' L2 'DO' Z[J]+:= U[J]*Z[J+1] 'OD'
     );
 
           ( 'REF'[,]'FLUXAC' DECL1 = GET(DEC[L1]);
                SOLL(DECL1[,-1],DECL1[,0],DECL1[,1],RH[L1,]);
                PET('FALSE',DEC[L1])
           );
           'FOR' I 'FROM' L1+1 'TO' U1
           'DO' 'REF'[,]'FLUXAC' DECI = GET(DEC[I]),
                                 JACI = GET(JAC[I]);
                'FOR' J 'FROM' L2 'TO' U2
                'DO' RH[I,J]-:= JACI[J,-2]*RH[I-1,J  ] 'OD';
                SOLL(DECI[,-1],DECI[,0],DECI[,1],RH[I,]);
                PET('FALSE',DEC[I]); PET('FALSE',JAC[I])
           'OD';
                'FOR' J 'FROM' L2 'TO' U2
                'DO' QQ[U1,J]+:= ( DUI[J] := RH[U1,J] ) 'OD';
           'FOR' I 'FROM' U1-1 'BY' -1 'TO' L1
           'DO' 'REF'[,]'FLUXAC' DECI = GET(DEC[I]),
                                 JACI = GET(JAC[I]);
                'FOR' J 'FROM' L2 'TO' U2
                'DO' DUI[J] := JACI[J, 2]*DUI[J] 'OD';
                SOLL(DECI[,-1],DECI[,0],DECI[,1],DUI);
                PET('FALSE',DEC[I]); PET('FALSE',JAC[I]);
                'FOR' J 'FROM' L2 'TO' U2
                'DO' QQ[I,J]+:= ( DUI[J] := RH[I,J]-DUI[J] )'OD'
           'OD'
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' ILLUDEC = ('NETMAT' JAC, 'REF''NETMAT' DECN )'VOID':
   # [L1:U1,L2:U2,-1:+1]'FLUXAC' DEC           #
   # [L1:U1]'LINMAT'             DEC           #
   #$B ILLUD   B$#
   'BEGIN' 'INT' L1 =  'LWB'JAC,     U1 =  'UPB'JAC;
           'INT' L2 = JL'OF'JAC[L1], U2 = JU'OF'JAC[L1];
           DECN:= 'HEAP'[L1:U1]'LINMAT';
           [L2:U2,-1:+1]'FLUXAC' D, DINV;
           'NETMAT' DEC = DECN;
           (L1/='LWB'DEC 'OR' U1/='UPB'DEC ! ERROR);
 
           'PROC' DECC = ('REF'[,]'FLUXAC' DEC, D)'VOID':
           (     DEC[L2, 0]:= 1/ D[L2, 0];
            'FOR' J 'FROM' L2+1 'TO' U2
            'DO' DEC[J ,-1]:=    D[J ,-1] -*DEC[J-1, 0];
                 DEC[J , 0]:= 1/(D[J , 0] + DEC[J  ,-1]*D[J-1,1])
            'OD' ;
            'FOR' J 'FROM' L2   'TO' U2-1
            'DO' DEC[J ,+1]:=    DEC[J,0] -*D[J,1] 'OD'
           );
 
           ('REF'[,]'FLUXAC' JACL1 = GET(JAC[L1]),
                             DECL1 = GEN(DEC[L1],L2,U2,-1,1);
            DECC(DECL1,JACL1)
           );
 
           'FOR' I 'FROM' L1 'TO' U1-1
           'DO'  'REF'[,]'FLUXAC' JACI  = GET(JAC[I  ]),
                                  JACIP = GET(JAC[I+1]),
                                  DECI  = GET(DEC[I  ]),
                                  DECIP = GEN(DEC[I+1],L2,U2,-1,1);
 
                 DINV[U2,0]    := DECI[U2,0];
                'FOR' J 'FROM' U2-1 'BY' -1 'TO' L2
                'DO' DINV[ J,0]:= DECI[J,0]
                                 +DECI[J,1]*DINV[J+1,0]*DECI[J+1,-1]
                'OD';
 
                'FOR' J 'FROM' L2+1 'TO' U2
                'DO'  DINV[J  ,-1]:= DINV[J  ,0]*DECI[J,-1];
                      DINV[J-1, 1]:= DECI[J-1,1]*DINV[J, 0]
                'OD';
 
                'FOR' K 'FROM'           -1   'TO'          1    'DO'
                'FOR' J 'FROM' (K=-1!L2+1!L2) 'TO' (K=1!U2-1!U2)
                'DO'  D[J,K] := JACIP[J, K]
                               -JACIP[J,-2]*DINV[J,K]*JACI[J+K,2]
                'OD'                                             'OD';
 
                DECC(DECIP,D); PET('FALSE',JAC[ I]); PET('TRUE',DEC[ I])
           'OD'              ; PET('FALSE',JAC[U1]); PET('TRUE',DEC[U1])
   'END'; #$E#
          #E$#
 
 # CONSTANTS AND 'NUL'            #                   'PR' EJECT 'PR'
 
   'REAL' GAM   = 1.4 ;
 
   'REAL' OGAM       =  1/GAM          ,  OTGAM   = 0.5/ GAM    ,
          GAMM1      =    GAM-1        ,  OGAMM1  =   1/(GAM-1) ,
          HGAMM1     =   (GAM-1)/2     , TOGAMM1  =   2/(GAM-1) ,
          TOGAMP1    = 2/(GAM+1)       ,
          OGAMGAMM1  = 1/(GAM*(GAM-1)) ,
          GAMM1OGAMP1= (GAM-1)/(GAM+1) ;
 
   'FLUX'   NUL       = (0,0,0,0);
   'FLUXAC' NULFLUXAC = (0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 );
 
   'OP' 'NUL' = ('REF''FLUXAC' A)'VOID' :  A:=NULFLUXAC ;
 
   #$H#
   'OP' 'NUL' = ('REF'[,]'STATE' A )'VOID':
   #$B NL.RRS  B$#
   'BEGIN'
     'INT' IL=1'LWB'A,IU=1'UPB'A,JL=2'LWB'A,JU=2'UPB'A ;
     'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE'AI=A[I,];
     'FOR' J 'FROM' JL 'TO' JU 'DO' AI[J] := NUL
     'OD' 'OD'
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'NUL' = ('REF'[]'FLUXAC' A )'VOID':
   #$B NL.RFA  B$#
    'FOR' I 'FROM' 'LWB'A 'TO' 'UPB'A
    'DO' A[I] := NULFLUXAC
    'OD' ; #$E#
           #E$#
 
   #$H#
   'OP' 'NUL' = ('REF'[,]'FLUXAC' A )'VOID':
   #$B NL.RRFA B$#
   'BEGIN'
     'INT' IL=1'LWB'A,IU=1'UPB'A,JL=2'LWB'A,JU=2'UPB'A ;
     'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'FLUXAC'AI=A[I,];
     'FOR' J 'FROM' JL 'TO' JU 'DO' AI[J] := NULFLUXAC
     'OD' 'OD'
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'NUL'    = ('FIELD' A)'FIELD':
   #$B NL.FLD  B$#
   'BEGIN' 'NUL'(Q'OF'A) ; A
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'NULRHS' = ('FIELD' Q)'FIELD':
   #$B NL.RHS  B$#
   'BEGIN' 'FIELD' F = 'NUL''COPY'Q ; TYPE'OF'F := 0 ; F
   'END' ; #$E#
           #E$#
 
                                                      'PR' EJECT 'PR'
   #$H#
   'OP' 'COPY' = ('FIELD' F)'FIELD':
   #$B COPY.FL B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'F;
           (DDX'OF'F, DDY'OF'F, 'HEAP''INT':= TYPE'OF'F,
            'HEAP'[1'LWB'Q:1'UPB'Q,2'LWB'Q:2'UPB'Q]'STATE':= Q )
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'NORM' = ('FIELD'F)[,]'REAL':
   #$B NORM B$#
   'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
           'INT' N1=1'UPB'Q, N2=2'UPB'Q;
           [1:4,-1:2]'REAL' NRM;
           'FOR' K 'TO' 4
           'DO' 'REF'[,]'REAL' QK=(K!C1'OF'Q,C2'OF'Q,C3'OF'Q,C4'OF'Q);
                'REAL' T, S1:= 0, S2:= 0, MAX:= 0;
                'FOR' I 'TO' N1 'DO' 'FOR'J 'TO' N2 'DO'
                      T:= 'ABS' QK[I,J]; S1+:=T; S2+:=T*T;
                      ( T>MAX ! MAX:=T )
                'OD' 'OD';
                NRM[K,'AT'1]:= (S1,MAX, S1/(N1*N2), SQRT(S2/(N1*N2)) )
           'OD';NRM
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'NORM1' = ('FIELD'F)[]'REAL':
   #$B NORM1 B$#
   'BEGIN' [,]'REAL'N= 'NORM'F; N [ , (TYPE'OF'F=0!-1!1) ]
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'PRINTNORM' = ('FIELD'F)'VOID':
   #$B NRM.PRT B$#
   'BEGIN' 'INT' TYPE=TYPE'OF'F;
           [,]'REAL' N = 'NORM' F;
           PRINT((NEWLINE,"INF-NORM:   ",N[, 0],
                  NEWLINE,"  2-NORM:   ",N[, 2],
                  NEWLINE,"  1-NORM:   ",N[,(TYPE=0!-1!1)],NEWLINE))
   'END' ; #$E#
           #E$#
 
 
   #$H#
   'OP' + = ('FIELD' A,B)'FIELD':
   #$B FLPLFL B$#
   'BEGIN' 'COPY' A +:= B
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' - = ('FIELD' A,B)'FIELD':
   #$B FLMNFL B$#
   'BEGIN' 'COPY' A -:= B
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' +:= = ('FIELD' A, 'FIELD' B)'FIELD':
   #$B PL.AFL  B$#
   'BEGIN' (TYPE'OF'A /= TYPE'OF'B ! ERROR );
     'REF'[,]'STATE' QA = Q'OF'A, QB = Q'OF'B;
     'INT'  L1=1'LWB'QB, U1=1'UPB'QB, L2=2'LWB'QB, U2=2'UPB'QB;
     'FOR' I 'FROM' L1 'TO' U1 'DO' 'REF'[]'STATE'AI=QA[I,],BI=QB[I,];
     'FOR' J 'FROM' L2 'TO' U2 'DO' 'REF''STATE' AIJ=AI[J],BIJ=BI[J];
       ( C1'OF'AIJ +:= C1'OF'BIJ , C2'OF'AIJ +:= C2'OF'BIJ ,
         C3'OF'AIJ +:= C3'OF'BIJ , C4'OF'AIJ +:= C4'OF'BIJ )
     'OD' 'OD';  A
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' -:= = ('FIELD' A, 'FIELD' B)'FIELD':
   #$B MN.AFL  B$#
   'BEGIN' (TYPE'OF'A /= TYPE'OF'B ! ERROR );
     'REF'[,]'STATE' QA = Q'OF'A, QB = Q'OF'B;
     'INT'  L1=1'LWB'QB, U1=1'UPB'QB, L2=2'LWB'QB, U2=2'UPB'QB;
     'FOR' I 'FROM' L1 'TO' U1 'DO' 'REF'[]'STATE'AI=QA[I,],BI=QB[I,];
     'FOR' J 'FROM' L2 'TO' U2 'DO' 'REF''STATE' AIJ=AI[J],BIJ=BI[J];
       ( C1'OF'AIJ -:= C1'OF'BIJ , C2'OF'AIJ -:= C2'OF'BIJ ,
         C3'OF'AIJ -:= C3'OF'BIJ , C4'OF'AIJ -:= C4'OF'BIJ )
     'OD' 'OD';  A
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' *:= = ('FIELD' A, 'REAL' B)'FIELD':
   #$B TM.AFL  B$#
    'IF' B = 1.0 'THEN' A
    'ELSE'
      'REF'[,]'STATE' QA = (Q'OF'A)[@1,@1];
      'INT'  N1 = 1'UPB'QA, N2 = 2'UPB'QA;
      'FOR' I 'TO' N1 'DO' 'REF'[]'STATE' AI=QA[I,];
      'FOR' J 'TO' N2 'DO' 'REF''STATE' AIJ =AI[J] ;
        ( C1'OF'AIJ*:=B,C2'OF'AIJ*:=B,C3'OF'AIJ*:=B,C4'OF'AIJ*:=B )
      'OD' 'OD' ; A
    'FI' ; #$E#
           #E$#
 
   #$H#
   'OP' * = ('NETMAT' JAC, 'FIELD' F) 'FIELD':
   #$B TIM     B$#
   'BEGIN'# PDZ #
          'REF'[,]'STATE' QF = Q'OF'F;
    'INT' L1 = 1'LWB'QF, U1 = 1'UPB'QF,
          L2 = 2'LWB'QF, U2 = 2'UPB'QF;
    (L1/='LWB'JAC 'OR' U1/='UPB'JAC ! ERROR);
    'HEAP'[L1:U1,L2:U2]'STATE' QG;
    'FIELD' G = (DDX'OF'F,DDY'OF'F, 'HEAP''INT':= 0, QG);
    'FOR' I 'FROM' L1 'TO' U1
    'DO'  'BOOL' LWB = (I=L1), UPB = (I=U1);
          'REF'[,]'FLUXAC' JACI = GET(JAC[I]),
          'REF' []'STATE'  QFM  = QF[(LWB!L1!I-1),],
                           QFI  = QF[        I   ,],
                           QFP  = QF[(UPB!U1!I+1),],
                           QGI  = QG[        I   ,];
          'INT' L3 = ('INT'L=2'LWB'JACI; (LWB!(L<-1!-1!L)!L)),
                U3 = ('INT'U=2'UPB'JACI; (UPB!(U> 1! 1!U)!U));
          'FOR' J 'FROM' L2 'TO' U2
          'DO'  'BOOL' NOTL = (J>L2), NOTR = (J<U2),
                'REF'[] 'FLUXAC' JACIJ = JACI[J,],
                'REF'   'STATE'  QGIJ  =  QGI[J]; QGIJ:=NUL ;
                'FOR' K 'FROM' L3 'TO' U3
                'DO' 'CASE' K+3
                     'IN'           QGIJ+:= JACIJ[K]*QFM[J  ] ,
                            (NOTL ! QGIJ+:= JACIJ[K]*QFI[J-1]),
                                    QGIJ+:= JACIJ[K]*QFI[J  ] ,
                            (NOTR ! QGIJ+:= JACIJ[K]*QFI[J+1]),
                                    QGIJ+:= JACIJ[K]*QFP[J  ]
                     'OUT'  ERROR
                     'ESAC'
                'OD'
          'OD'  ;
          PET('FALSE',JAC[I])
    'OD' ; G
   'END'; #$E#
          #E$#
 
 # FIELD MANIPULATION             #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' INITIATE = ('FIELD' F, 'STATE' A, 'INT' TYPE )'FIELD':
   #$B INI.FLD B$#
   'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
           TYPE'OF'F:=TYPE ; 'INT' N1=1'UPB'Q,N2=2'UPB'Q;
           'FOR' I 'TO' N1 'DO' 'REF'[]'STATE' FI=Q[I,];
           'FOR' J 'TO' N2 'DO' FI[J] := A
           'OD' 'OD'; F
   'END'; #$E#
          #E$#
 
   'PROC'('REAL','REAL')'REAL' F0,F1,F2,F3,F4;
 
   #$H#
   'PROC' SETFIELD = ('FIELD' QQ, 'INT' TYPE, 'BOOL' STAR,
                        'PROC'('REAL','REAL')'REAL' F1,F2,F3,F4)'FIELD':
   #$B SET.FLD B$#
   'BEGIN'
     'REAL' DDX = DDX'OF'QQ, DDY = DDY'OF'QQ   ; TYPE'OF'QQ:=TYPE ;
     'POINT' P; 'REF''REAL' X=X'OF'P, Y=Y'OF'P ; 'REAL' XX, YY    ;
     'PROC' HERE = ('REAL' XX,YY)'STATE':
         ( P:= MAPPING (XX,YY);
           'STATE'(F1(X,Y),F2(X,Y),F3(X,Y),F4(X,Y))  );
     'REF'[,]'STATE' Q  = Q'OF'QQ;
     'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
     'INT'  ILH = (IL+1)'OVER'2, IUH = IU'OVER'2,
            JLH = (JL+1)'OVER'2, JUH = JU'OVER'2;
     (STAR
     ! 'IF' 'NOT' ( ILH*2 = IL+1 'AND' IUH*2 = IU 'AND'
                    JLH*2 = JL+1 'AND' JUH*2 = JU )
       'THEN' ERROR; QQ
       'ELSE' 'REAL' X23 = 2*DDX/3, Y23 = 2*DDY/3;
              'FOR' IH 'FROM' ILH 'TO' IUH
              'DO' 'INT'I = 2*IH-1; XX:= I*DDX ;
                   'REF'[]'STATE' QI=Q[I,],QIP=Q[I+1,];
                   'FOR' JH 'FROM' JLH 'TO' JUH
                   'DO' 'INT'J = 2*JH-1; YY:= J*DDY; 'INT'JP=J+1;
                        QI  [J ] := HERE ( XX     , YY-Y23 ) ;
                        QI  [JP] := HERE ( XX-X23 , YY     ) ;
                        QIP [JP] := HERE ( XX     , YY+Y23 ) ;
                        QIP [J ] := HERE ( XX+X23 , YY     )
              'OD''OD'; QQ
       'FI'
     !  'FOR' I 'FROM' IL  'TO' IU
        'DO' XX:= (I-0.5)*DDX; 'REF'[]'STATE'QI=Q[I,] ;
             'FOR' J 'FROM' JL 'TO' JU
             'DO'  QI[J]:= HERE( XX,(J-0.5)*DDY )
        'OD' 'OD'; QQ
     ) 'END' ; #$E#
               #E$#
 
   #$H#
   'OP' 'RESTRICT'    = ('FIELD' A)'FIELD':
   #$B RESTR   B$#
   'BEGIN' 'REAL' W = ( TYPE'OF'A ! 0.25,0.25 ! 1.0);
           'REF'[,]'STATE' Q  = Q'OF'A;
           'INT' L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
           'INT' LL2:= (L2+1)'OVER'2, UU2:= U2'OVER'2, TI, TJ,
                 LL3:= (L3+1)'OVER'2, UU3:= U3'OVER'2;
           'HEAP'[LL2:UU2,LL3:UU3]'STATE' QN;
           'FOR' I 'FROM' LL2 'TO' UU2 'DO' TI:= I+I;
                   'REF'[]'STATE' QNI=QN[I,],QTIM=Q[TI-1,],QTI=Q[TI,];
           'FOR' J 'FROM' LL3 'TO' UU3 'DO' TJ:= J+J;
               QNI[J]:= W* (QTIM[TJ-1]+QTIM[TJ]+QTI[TJ-1]+QTI[TJ])
           'OD' 'OD';
           'FIELD'(2*(DDX'OF'A),2*(DDY'OF'A),'HEAP''INT':=TYPE'OF'A,QN)
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'PROLON' = ('FIELD'A)'FIELD':
   #$B PROL.1  B$#
   'BEGIN' 'REF'[,]'STATE' Q  = Q'OF'A;
           'INT' L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
           'INT' TI, TJ;
           'HEAP'[2*L2-1:2*U2,2*L3-1:2*U3]'STATE' QN;
           'FOR' I 'FROM' L2 'TO' U2 'DO' TI:= I+I;
               'REF'[]'STATE' QNTIM=QN[TI-1,],QNTI=QN[TI,],QI=Q[I,];
           'FOR' J 'FROM' L3 'TO' U3 'DO' TJ:= J+J;
              QNTIM[TJ-1]:=QNTIM[TJ]:=QNTI[TJ-1]:=QNTI[TJ]:=QI[J]
           'OD' 'OD';
           'FIELD'((DDX'OF'A)/2,(DDY'OF'A)/2,'HEAP''INT':=TYPE'OF'A,QN)
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'OP' 'PROLON2' = ('FIELD' F)'FIELD':
   #$B PROLO2  B$#
   'BEGIN' 'REF'[,]'STATE' QQ  = (Q'OF' F);
           'INT'  L2= 1'LWB'QQ, U2= 1'UPB'QQ,
                  L3= 2'LWB'QQ, U3= 2'UPB'QQ;
    'IF'   'ODD'(U2-L2+1) 'OR' 'ODD'(U3-L3+1)
    'THEN' 'PROLON' F
    'ELSE' 'HEAP'[2*L2-1:2*U2,2*L3-1:2*U3]'STATE'QQQ;
           'INT'  TI, TJ, PI, PJ;
           'REAL' X,Y;
           'STATE' A,B,C,D,P,Q,R,S;
           'FOR' I 'FROM' L2+1 'BY' 2 'TO' U2 'DO' TI:= I+I;
                   'REF'[]'STATE' QQIM=QQ[I-1,],QQI=QQ[I,] ;
           'FOR' J 'FROM' L3+1 'BY' 2 'TO' U3 'DO' TJ:= J+J;
                    A:= QQIM[J-1]    ; B := QQI[J-1]; C := QQIM[J];
                    D:= QQI [J]-B-C+A; B-:= A ; C-:= A ; X:=-0.75 ;
           'FOR' K 'TO' 4 'DO' P:=C+(X+:=0.5)*D ; Q:=A+X*B ;
                          Y:=-0.75;'REF'[]'STATE'QQQK=QQQ[TI-4+K,];
           'FOR' L 'TO' 4 'DO' QQQK[TJ-4+L] := Q+(Y+:=0.5)*P
           'OD' 'OD'
           'OD' 'OD';
           'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
                    QQQ)
    'FI'
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'PROLIN2' = ('FIELD' F)'FIELD':
   #$B PROLI2  B$#
   'BEGIN' 'REF'[,]'STATE' QQ  = (Q'OF' F);
           'INT'  L1= 1'LWB'QQ, U1= 1'UPB'QQ,
                  L2= 2'LWB'QQ, U2= 2'UPB'QQ;
    'IF'   (L1=U1)'OR'(L2=U2)
    'THEN' 'PROLON'F
    'ELSE' 'HEAP'[2*L1-1:2*U1,2*L2-1:2*U2]'STATE'QQQ;
           'INT'  TI, TJ, IP, JP,
           'STATE' P,Q,R,S,
           'BOOL' LWBI, UPBI;
           #  25  15   5  -5  #
           #  15   9   3  -3  #   # /16 #
           #   5   3   1  -1  #
           #  -5  -3  -1   1  #
           'REAL' Z25 = 25/16, Z15 = 15/16, Z9 = 9/16,
                   Z5 =  5/16,  Z3 =  3/16, Z1 = 1/16;
           'FOR' I 'FROM' L1+1  'TO' U1
           'DO' IP:=2*I-2;
                LWBI:=(I=L1+1);
                UPBI:=(I=U1  );
                'FOR' J 'FROM' L2+1  'TO' U2
                'DO' JP:= 2*J-2;
                     P:= QQ[I-1,J-1];    Q:= QQ[I  ,J-1];
                     R:= QQ[I-1,J  ];    S:= QQ[I  ,J  ];
                     QQQ[IP  ,JP  ]:=Z9*P+Z3*(Q+R)+Z1*S;
                     QQQ[IP+1,JP  ]:=Z9*Q+Z3*(P+S)+Z1*R;
                     QQQ[IP  ,JP+1]:=Z9*R+Z3*(P+S)+Z1*Q;
                     QQQ[IP+1,JP+1]:=Z9*S+Z3*(Q+R)+Z1*P;
                     (J=L2+1!(LWBI ! QQQ[IP-1,JP-1]:=
                                     Z25*P-Z5*(Q+R)+Z1*S);
                             (UPBI ! QQQ[IP+2,JP-1]:=
                                     Z25*Q-Z5*(P+S)+Z1*R);
                              QQQ[IP  ,JP-1]:=
                              Z15*P+Z5*Q-Z3*R-Z1*S;
                              QQQ[IP+1,JP-1]:=
                              Z15*Q+Z5*P-Z3*S-Z1*R    );
                     (J=U2 ! (LWBI ! QQQ[IP-1,JP+2]:=
                                     Z25*R-Z5*(P+S)+Z1*Q);
                             (UPBI ! QQQ[IP+2,JP+2]:=
                                     Z25*S-Z5*(Q+R)+Z1*P);
                              QQQ[IP  ,JP+2]:=
                              Z15*R+Z5*S-Z3*P-Z1*Q;
                              QQQ[IP+1,JP+2]:=
                              Z15*S+Z5*R-Z3*Q-Z1*P    );
                     (LWBI ! QQQ[IP-1,JP  ]:=
                              Z15*P+Z5*R-Z3*Q-Z1*S;
                             QQQ[IP-1,JP+1]:=
                              Z15*R+Z5*P-Z3*S-Z1*Q    );
                     (UPBI ! QQQ[IP+2,JP  ]:=
                              Z15*Q+Z5*S-Z3*P-Z1*R;
                             QQQ[IP+2,JP+1]:=
                              Z15*S+Z5*Q-Z3*R-Z1*P    )
           'OD' 'OD';
           'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
                    QQQ)
    'FI'
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'CLLPNT' = ('FIELD' F)'FIELD':
   #$B CLLPNT B$#
   'BEGIN' 'FIELD' WF:= ADD WALL('TRUE',F);
      'REF'[,]'STATE' QF  = (Q'OF' WF);
      'INT'  L1= 1'LWB'QF, U1= 1'UPB'QF,
             L2= 2'LWB'QF, U2= 2'UPB'QF;
      'INT' U1M = U1-1, U2M = U2-1;
      'HEAP'[L1:U1M,L2:U2M]'STATE' QPF;
 
      'FOR' I 'FROM' L1 'TO' U1M
      'DO'  'FOR' J 'FROM' L2 'TO' U2M
            'DO'  QPF[I,J]:= 0.25*( QF[I  ,J  ]+QF[I  ,J+1]+
                                    QF[I+1,J  ]+QF[I+1,J+1] )
            'OD'
      'OD'  ;
 
      Q'OF'WF:='NIL';
      'FIELD'( DDX'OF'F, DDY'OF'F, 'HEAP''INT':=TYPE'OF'F, QPF)
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'PROSCND' = ('FIELD' C) 'FIELD':
   #$B PROSCND B$#
   'IF'   'REF'[,]'STATE' QC  = (Q'OF' C);
          'ODD'(1'UPB'QC-1'LWB'QC+1)'OR''ODD'(2'UPB'QC-2'LWB'QC+1)
   'THEN' 'PROLON'C
   'ELSE' 'FIELD' PC:= 'CLLPNT'C;
          'REF'[,]'STATE' QP  = (Q'OF'PC);
          'INT'  L1= 1'LWB'QP, U1= 1'UPB'QP,
                 L2= 2'LWB'QP, U2= 2'UPB'QP;
          'HEAP'[2*(L1+1)-1:2*U1,2*(L2+1)-1:2*U2]'STATE'QF;
          'STATE' S11, S12, S13, S21, S22, S23, S31, S32, S33;
          'REAL' Z1:=1/16, Z2:=1/8,  Z3 :=3/16, Z4 := 1/4,  Z5 := 5/16,
                 Z6:=3/8,  Z9:=9/16, Z10:=5/8,  Z15:=15/16, Z25:=25/16;
          ( TYPE'OF'C = 0 ! Z1 *:=0.25; Z2 *:=0.25; Z3 *:=0.25;
                            Z4 *:=0.25; Z5 *:=0.25; Z6 *:=0.25;
                            Z9 *:=0.25; Z10*:=0.25; Z15*:=0.25;
                            Z25*:=0.25 );
          'INT' I1:= 2*(L1+1)-1;
          'FOR' I 'FROM' L1 'BY' 2 'TO' (U1-2)
          'DO'  'INT' I2 = I1+1, I3 = I1+2, I4 = I1+3;
                S13:= QP[I  ,L2]; S23:= QP[I+1,L2]; S33:= QP[I+2,L2];
                'INT' J1:= 2*(L2+1)-1;
                'FOR' J 'FROM' L2 'BY' 2 'TO' (U2-2)
                'DO'  'INT' J2 = J1+1, J3 = J1+2, J4 = J1+3;
                      S11:= S13; S12:= QP[I  ,J+1]; S13:= QP[I  ,J+2];
                      S21:= S23; S22:= QP[I+1,J+1]; S23:= QP[I+1,J+2];
                      S31:= S33; S32:= QP[I+2,J+1]; S33:= QP[I+2,J+2];
                      QF[I1,J1]:=   Z4 *S11 + Z6 *S12 - Z2 *S13
                                  + Z6 *S21 + Z9 *S22 - Z3 *S23
                                  - Z2 *S31 - Z3 *S32 + Z1 *S33;
                      QF[I1,J2]:=             Z10*S12 - Z2 *S13
                                            + Z15*S22 - Z3 *S23
                                            - Z5 *S32 + Z1 *S33;
                      QF[I1,J3]:= - Z2 *S11 + Z10*S12
                                  - Z3 *S21 + Z15*S22
                                  + Z1 *S31 - Z5 *S32;
                      QF[I1,J4]:= - Z2 *S11 + Z6 *S12 + Z4 *S13
                                  - Z3 *S21 + Z9 *S22 + Z6 *S23
                                  + Z1 *S31 - Z3 *S32 - Z2 *S33;
                      QF[I2,J1]:=   Z10*S21 + Z15*S22 - Z5 *S23
                                  - Z2 *S31 - Z3 *S32 + Z1 *S33;
                      QF[I2,J2]:=             Z25*S22 - Z5 *S23
                                            - Z5 *S32 + Z1 *S33;
                      QF[I2,J3]:= - Z5 *S21 + Z25*S22
                                  + Z1 *S31 - Z5 *S32;
                      QF[I2,J4]:= - Z5 *S21 + Z15*S22 + Z10*S23
                                  + Z1 *S31 - Z3 *S32 - Z2 *S33;
                      QF[I3,J1]:= - Z2 *S11 - Z3 *S12 + Z1 *S13
                                  + Z10*S21 + Z15*S22 - Z5 *S23;
                      QF[I3,J2]:=           - Z5 *S12 + Z1 *S13
                                            + Z25*S22 - Z5 *S23;
                      QF[I3,J3]:=   Z1 *S11 - Z5 *S12
                                  - Z5 *S21 + Z25*S22;
                      QF[I3,J4]:=   Z1 *S11 - Z3 *S12 - Z2 *S13
                                  - Z5 *S21 + Z15*S22 + Z10*S23;
                      QF[I4,J1]:= - Z2 *S11 - Z3 *S12 + Z1 *S13
                                  + Z6 *S21 + Z9 *S22 - Z3 *S23
                                  + Z4 *S31 + Z6 *S32 - Z2 *S33;
                      QF[I4,J2]:=           - Z5 *S12 + Z1 *S13
                                            + Z15*S22 - Z3 *S23
                                            + Z10*S32 - Z2 *S33;
                      QF[I4,J3]:=   Z1 *S11 - Z5 *S12
                                  - Z3 *S21 + Z15*S22
                                  - Z2 *S31 + Z10*S32;
                      QF[I4,J4]:=   Z1 *S11 - Z3 *S12 - Z2 *S13
                                  - Z3 *S21 + Z9 *S22 + Z6 *S23
                                  - Z2 *S31 + Z6 *S32 + Z4 *S33;
                      J1+:= 4
                'OD' ;
                I1+:= 4
          'OD'  ;
           Q'OF'PC:= 'NIL';
          'FIELD'((DDX'OF'C)/2,(DDY'OF'C)/2,'HEAP''INT':=TYPE'OF'C,QF)
   'FI' ; #$E#
          #E$#
 # FIRST-ORDER STAR-INTERPOLATION #                   'PR' EJECT 'PR'
 
   #$H#
   'OP' 'STARPROL' = ('FIELD' F)'FIELD':
   #$B STRPRL  B$#
   'IF' 'REF'[,]'STATE' QC  = Q'OF'F;
        'INT' IL= 1'LWB'QC, IU= 1'UPB'QC,
              JL= 2'LWB'QC, JU= 2'UPB'QC;
        'ODD'(IU-IL+1) 'OR' 'ODD'(JU-JL+1)
   'THEN' 'PROLON' F
   'ELSE' 'HEAP'[2*IL-1:2*IU,2*JL-1:2*JU]'STATE'QF;
          'INT'  TI, TJ;
          'FOR' I 'FROM' IL+1 'BY' 2 'TO' IU 'DO' TI:= I+I;
          'FOR' J 'FROM' JL+1 'BY' 2 'TO' JU 'DO' TJ:= J+J;
                QF[TI-3,TJ-2]:= QF[TI-2,TJ-2]:=
                QF[TI-3,TJ-1]:= QF[TI-3,TJ  ]:= QC[I-1,J  ];
                QF[TI-2,TJ  ]:= QF[TI-2,TJ-1]:=
                QF[TI-1,TJ  ]:= QF[TI  ,TJ  ]:= QC[I  ,J  ];
                QF[TI  ,TJ-1]:= QF[TI-1,TJ-1]:=
                QF[TI  ,TJ-2]:= QF[TI  ,TJ-3]:= QC[I  ,J-1];
                QF[TI-1,TJ-3]:= QF[TI-1,TJ-2]:=
                QF[TI-2,TJ-3]:= QF[TI-3,TJ-3]:= QC[I-1,J-1]
          'OD' 'OD';
          'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
                   QF)
   'FI'; #$E#
         #E$#
 # STATE FUNCTIONS                #                   'PR' EJECT 'PR'
 
   # THESE FUNCTIONS ASSUME STATE (1) #
 
   'PROC' DENSITY     = ('STATE' Q)'REAL' : C1'OF'Q ;
   'PROC' TOTENERGY   = ('STATE' Q)'REAL' : C4'OF'Q ;
   'PROC' XVELOCITY   = ('STATE' Q)'REAL' : C2'OF'Q/C1'OF'Q ;
   'PROC' YVELOCITY   = ('STATE' Q)'REAL' : C3'OF'Q/C1'OF'Q ;
   'PROC' XMOMENT     = ('STATE' Q)'REAL' : C2'OF'Q ;
   'PROC' YMOMENT     = ('STATE' Q)'REAL' : C3'OF'Q ;
   'PROC' LOCMACH     = ('STATE' Q)'REAL': SPEED(Q)/SOUND(Q);
   'PROC' PRESSURE    = ('STATE' Q)'REAL':
                    ('REAL' R=C1'OF'Q; 'REAL' RU=C2'OF'Q, RV=C3'OF'Q;
                            ((RU*RU+RV*RV)/(-2*R) + C4'OF'Q)*GAMM1 );
   'PROC' DIREC       = ('STATE' Q)'REAL':
                    ('REAL' U=C2'OF'Q, V= C3'OF'Q;
                     'REAL' D=(U=0!'SIGN'V!2*ARCTAN(V/U)/PI);
                               U<0! ( V<0 ! D-2 ! D+2) ! D );
   'PROC' ENTROPY = ('STATE' Q)'REAL':
          ('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
            ((RU*RU+RV*RV)/(-2*R)+ C4'OF'Q)*GAMM1* R**-GAM);
   'PROC' TOTENTHALPY = ('STATE' Q)'REAL':
          ('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
            ((RU*RU+RV*RV)*GAMM1/(-2*R)+ GAM*C4'OF'Q) /R  );
   'PROC' SOUND   = ('STATE' Q)'REAL': SQRT(GAM*GAMM1*TEMPERAT(Q));
   'PROC' SPEED   = ('STATE' Q)'REAL':
               SQRT(C2'OF'Q*C2'OF'Q+C3'OF'Q*C3'OF'Q)/C1'OF'Q;
   'PROC' TEMPERAT= ('STATE' Q)'REAL':
          ('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
          # I.E. P/(R*GAMM1) #   ((RU*RU+RV*RV)/(-2*R) + C4'OF'Q)/R );
 
   []'PROC'('STATE')'REAL' VAR =
                        (DENSITY, XMOMENT, YMOMENT, TOTENERGY,
                         XVELOCITY, YVELOCITY, SOUND, ENTROPY,
                         PRESSURE, SPEED, LOCMACH, DIREC,
                         TOTENTHALPY, TEMPERAT);
 
   'MODE' 'VAR'   = 'STRUCT'('BOOL' INTERN, 'REAL' DDX,DDY,
                                  'INT' VAR, 'REF'[,]'REAL'Q);
 
   #$H#
   'PROC' VARIABLE = ('STRING' NAME, 'BOOL' INTERN, 'FIELD' F) 'VAR':
   #$B VAR     B$#
   'BEGIN' 'INT' VARNO:= 0;
                (  NAME >  "YVELOCITY"   !  ERROR; VARNO:= 0 );
                (  NAME <= "YVELOCITY"   !  VARNO:=  6);
                (  NAME <= "YMOMENT"     !  VARNO:=  3);
                (  NAME <= "XVELOCITY"   !  VARNO:=  5);
                (  NAME <= "XMOMENT"     !  VARNO:=  2);
                (  NAME <= "TOTENTHALPY" !  VARNO:= 13);
                (  NAME <= "TOTENERGY"   !  VARNO:=  4);
                (  NAME <= "TEMPERATURE" !  VARNO:= 14);
                (  NAME <= "SPEED"       !  VARNO:= 10);
                (  NAME <= "SOUND"       !  VARNO:=  7);
                (  NAME <= "PRESSURE"    !  VARNO:=  9);
                (  NAME <= "LOCMACH"     !  VARNO:= 11);
                (  NAME <= "ENTROPY"     !  VARNO:=  8);
                (  NAME <= "DIRECTION"   !  VARNO:= 12);
                (  NAME <= "DENSITY"     !  VARNO:=  1);
                (  NAME  = ""            !  ERROR; VARNO:=  0);
 
           ( TYPE'OF'F /=0 ! STATE E(F) );
           [,]'STATE' Q = (Q'OF'F)[@1,@1];
           'INT' L1Q = 1'LWB'Q, U1Q= 1'UPB'Q,
                 L2Q = 2'LWB'Q, U2Q= 2'UPB'Q;
           'HEAP'[L1Q:U1Q,L2Q:U2Q]'REAL' V;
                  'PROC'('STATE')'REAL' COMP = VAR[VARNO];
                  'FOR' I 'FROM' L1Q 'TO' U1Q 'DO'
                  'FOR' J 'FROM' L2Q 'TO' U2Q 'DO' V[I,J]:= COMP(Q[I,J])
                  'OD' 'OD';
           'VAR'(INTERN,DDX'OF'F, DDY'OF'F, VARNO, V)
   'END'; #$E#
          #E$#
 
 # FLUX FUNCTION                  #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' ZFLUX = ('STATE' Q, 'REF''FLUXAC' FAC )'FLUX':
   #$B ZFLUX   B$#
   'BEGIN' 'REAL' U:=C1'OF'Q,V:=C2'OF'Q,C:=C3'OF'Q,Z:=C4'OF'Q,
                  C2,U2,POR,R,RU,RUV,RUUP,W2,EPR,RUEPR,OC,OCU,OS;
                  C2   := C*C;
                  U2   := U*U;
                  POR  := C2*OGAM;
                  R    := EXP((LN(POR)-Z)*OGAMM1);
                  RU   := R*U;
                  RUV  := RU*V;
                  RUUP := R*(U2+POR);
                  W2   := (V*V+U2)*0.5;
                  EPR  := C2*OGAMM1 + W2;
                  RUEPR:= RU*EPR;
     ( FAC 'ISNT''NIL'
     ! OS:= -OGAMM1; OC:= TOGAMM1*R/C; OCU:= OC*U;
       FAC := 'FLUXAC'
               (          R,  0, OCU         , OS*RU    ,
                       2*RU,  0, OC * (U2+C2), OS*RUUP  ,
                        R*V, RU, OCU*      V , OS*RUV   ,
                 (U2+EPR)*R,RUV, OCU*(EPR+C2), OS*RUEPR )
     );
     'FLUX'( (RU, RUUP, RUV, RUEPR) )
   'END';#$E#
         #E$#
 
 # SIMULTANEOUS STATE FUNCTIONS   #                   'PR' EJECT 'PR'
 
   'BOOL' WARNING:= 'TRUE';
 
   'PROC' REVEAL = ('REAL' C,Z, 'REF''REAL' P,R)'VOID':
          ('REAL'POR:= C*C/GAM; R:= EXP((LN(POR)-Z)*OGAMM1); P:= POR*R);
 
   #$H#
   'OP' 'EZ' = ('STATE' A)'STATE':
   #$B EZ      B$#
   ('REAL' R ,U:=C1'OF'A, V:= C2'OF'A, P,
           POR:= C3'OF'A*C3'OF'A/GAM;
           R:= EXP((LN(POR)-C4'OF'A)*OGAMM1); P:= POR*R;
    'STATE' (R,R*U,R*V,(U*U+V*V)*R/2 + P/GAMM1)
   ); #$E#
      #E$#
 
   #$H#
   'OP' 'ZE' = ('STATE' A)'STATE':
   #$B ZE      B$#
   ('REAL' R:=  C1'OF'A, M:= C2'OF'A, N:= C3'OF'A, P;
           P:= ((M*M+N*N)*-0.5/R + C4'OF'A)*GAMM1;
    'STATE' (M/R,N/R,SQRT(GAM*P/R),LN(P)-GAM*LN(R))
   ); #$E#
      #E$#
 
   #$H#
   'OP' 'EP' = ([]'REAL' A)'STATE':
   #$B EP      B$#
   ('REAL' R = A[1] ; 'REAL' RU = R*A[2] , RV = R*A[3] ;
    'STATE'(R,RU,RV, (RU*RU+RV*RV)*0.5/R + A[4]*OGAMM1)
   ); #$E#
      #E$#
 
   #$H#
   'OP' 'PE' = ('STATE' A)[]'REAL':
   #$B PE      B$#
   ('REAL' R:=  C1'OF'A; 'REAL' U = C2'OF'A/R, V = C3'OF'A/R;
    'REAL' P:= ((U*U+V*V)*-0.5*R + C4'OF'A)*GAMM1;
    ( P<=0  'OR' R<=0 ! R:= P:= 1.0E-10;
    ( WARNING         ! PRINT((NEWLINE,"NEGATIVE RHO OR P "+50*"!"));
                        WARNING := 'FALSE' ));
    []'REAL'(R,U,V,P, SQRT(GAM*P/R), LN(P)-GAM*LN(R) )
   ); #$E#
      #E$#
 
   #$H#
   'OP' 'PZ' = ('STATE' A)[]'REAL':
   #$B PZ      B$#
    ( 'REAL' P,R; REVEAL(C3'OF'A,C4'OF'A,P,R);
      []'REAL'(R,C1'OF'A,C2'OF'A,P,C3'OF'A,C4'OF'A)
    ); #$E#
       #E$#
 
   #$H#
   'OP' 'ZP' = ([]'REAL' A)'STATE':
   #$B ZP      B$#
   ( 'ZE''EP' A
   ); #$E#
      #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' STATE Z = ('FIELD' F)'FIELD':
   #$B STT.Z   B$#
   'IF'   TYPE'OF'F = 0  'THEN' ERROR; 'SKIP'
   'ELIF' TYPE'OF'F = 1
   'THEN' TYPE'OF'F:= 2;
     'REF'[,]'STATE'Q = Q'OF'F;
     'FOR' I 'FROM' 1'LWB'Q 'TO' 1'UPB'Q 'DO' 'REF'[]'STATE' QI=Q[I,];
     'FOR' J 'FROM' 2'LWB'Q 'TO' 2'UPB'Q 'DO' QI[J]:='ZE'QI[J]
     'OD''OD'; F
   'ELSE' F
   'FI';#$E#
        #E$#
 
   #$H#
   'PROC' STATE E = ('FIELD' F)'FIELD':
   #$B STT.E   B$#
   'IF'   TYPE'OF'F = 0  'THEN' ERROR; 'SKIP'
   'ELIF' TYPE'OF'F = 2
   'THEN' TYPE'OF'F:= 1;
          'REF'[,]'STATE'Q = Q'OF'F;
     'FOR' I 'FROM' 1'LWB'Q 'TO' 1'UPB'Q 'DO' 'REF'[]'STATE'QI=Q[I,];
     'FOR' J 'FROM' 2'LWB'Q 'TO' 2'UPB'Q 'DO' QI[J]:='EZ'QI[J]
     'OD''OD'; F
   'ELSE' F
   'FI';#$E#
        #E$#
 
 # NUMERICAL FLUX                 #                   'PR' EJECT 'PR'
 
   'PROC' NFLUX:= ('STATE' QO,Q1, 'REF''FLUXAC' FAC0,FAC1)'FLUX':
                  ( ERROR; NUL );
   'INT' SOS:= 1;
   'BOOL' CAVITATION:= 'FALSE';
 
   #$H#
   'PROC' OSHER = ('STATE' Q0,Q1,'REF' 'FLUXAC' FAC0,FAC1)'FLUX':
   #$B OSHER   B$#
   'BEGIN' 'REAL' U0=C1'OF'Q0,V0=C2'OF'Q0,C0=C3'OF'Q0,Z0=C4'OF'Q0,
                  U1=C1'OF'Q1,V1=C2'OF'Q1,C1=C3'OF'Q1,Z1=C4'OF'Q1,
                  SGN  = SOS,    SHGAMM1 = SOS*HGAMM1;
           'BOOL' FACP = FAC1 'ISNT''NIL',
                  FACM = FAC0 'ISNT''NIL';
           'REAL' PSI0 = U0-C0/SHGAMM1,
                  PSI1 = U1+C1/SHGAMM1,
                  ALFA = EXP((Z1-Z0)*OTGAM);
           'REAL' OALFAP1=1/(1+ALFA);
           'REAL' UH   = (PSI1+ALFA*PSI0)*OALFAP1 +1.0E-200,
                  CL   = SHGAMM1*(PSI1-PSI0)*OALFAP1;
                  (CL<0! # CAVITATION # ERROR );
           'REAL' CR   = ALFA*CL,
                  LAM0 = U0+SGN*C0,  LAML = UH+SGN*CL,
                  LAM1 = U1-SGN*C1,  LAMR = UH-SGN*CR;
 
      'FLUX'        FLUC ,          F   := NUL;
      'REF''FLUXAC' FLAC = (FACM 'OR' FACP
                           ! (FACP !FAC1:=NULFLUXAC);
                             (FACM !FAC0:=NULFLUXAC);
                             'LOC''FLUXAC'
                           ! 'NIL'                  );
 
      (UH*LAM0*LAML*LAMR*LAM1=0 ! ERROR);
      (CL<0 ! CAVITATION := 'TRUE' ; EXCEPTION ) ;
 
      'IF' LAM0>0
      'THEN' (FACM ! F+:=ZFLUX(Q0, FAC0)
                   ! F+:=ZFLUX(Q0,'NIL') )
      'FI';
 
      'IF' LAM1<0
      'THEN' (FACP ! F+:=ZFLUX(Q1, FAC1)
                   ! F+:=ZFLUX(Q1,'NIL') )
      'FI';
 
      'IF' LAM0*LAML < 0
      'THEN' 'REAL' U= GAMM1OGAMP1*PSI0;
             FLUC:= ZFLUX((U,V0,-SGN*U,Z0),(FACM ! FLAC ! 'NIL'));
             (FACM!'FLUX' K = (1'FC'FLAC)-SGN*(3'FC'FLAC) ;
                   CASS(FLAC,1,GAMM1OGAMP1   *K ) ;
                   CASS(FLAC,3,(-SGN*TOGAMP1)*K ) ;
                   (LAM0<0 ! FAC0+:=FLAC ! FAC0-:=FLAC ));
             (LAM0<0 ! F+:=FLUC ! F-:=FLUC )
      'FI';
 
      'IF' LAMR*LAM1 < 0
      'THEN' 'REAL' U= GAMM1OGAMP1*PSI1;
             FLUC:= ZFLUX((U,V1,SGN*U,Z1),(FACP ! FLAC ! 'NIL' ));
             (FACP!'FLUX' K= 1'FC'FLAC+SGN*3'FC'FLAC;
                   CASS(FLAC,1,GAMM1OGAMP1*K);
                   CASS(FLAC,3,(SGN*TOGAMP1)*K);
                   (LAM1>0 ! FAC1+:=FLAC ! FAC1-:=FLAC ));
             (LAM1>0 ! F+:=FLUC ! F-:=FLUC )
      'FI';
 
      'IF' LAML*UH<0
      'THEN' FLUC:= ZFLUX((UH,V0,CL,Z0),FLAC);
             (FACM!'REAL' DUM= CR*OALFAP1*OGAMGAMM1;
                   'FLUX' K1= (ALFA*OALFAP1)*(1'FC'FLAC)
                             -(SHGAMM1*OALFAP1)*(3'FC'FLAC),
                          K2= DUM*(1'FC'FLAC+SHGAMM1*(3'FC'FLAC));
                   CSUB(FAC0,1,SGN*K1            );
                   CSUB(FAC0,2,SGN*(2'FC'FLAC)   );
                   CSUB(FAC0,3,-TOGAMM1*K1       );
                   CSUB(FAC0,4,SGN*(4'FC'FLAC)+K2) );
             (FACP!'REAL' DUM=-CR*OGAMGAMM1;
                   'FLUX' K = OALFAP1*
                              (1'FC'FLAC+ SHGAMM1 *(3'FC'FLAC));
                   CSUB(FAC1,1, SGN    *K);
                   CSUB(FAC1,3, TOGAMM1*K);
                   CSUB(FAC1,4, DUM    *K) );
             (SGN>0! F-:= FLUC ! F+:= FLUC)
      'ELIF' UH*LAMR<0
      'THEN' FLUC:= ZFLUX((UH,V1,CR,Z1),FLAC);
             (FACM!'REAL' DUM= CL*OGAMGAMM1;
                   'FLUX'K= ALFA*OALFAP1*
                              (1'FC'FLAC- SHGAMM1 *(3'FC'FLAC));
                   CSUB(FAC0,1,        SGN*K);
                   CSUB(FAC0,3,   -TOGAMM1*K);
                   CSUB(FAC0,4,        DUM*K) );
             (FACP!'REAL' DUM=-CR*OALFAP1*OGAMGAMM1;
                   'FLUX' K1= OALFAP1*
                              (1'FC'FLAC+(SHGAMM1*ALFA)*3'FC'FLAC),
                          K2= DUM*(1'FC'FLAC- SHGAMM1  *3'FC'FLAC);
                   CSUB(FAC1,1,            SGN*K1);
                   CSUB(FAC1,2,     SGN*2'FC'FLAC);
                   CSUB(FAC1,3,        TOGAMM1*K1);
                   CSUB(FAC1,4,  SGN*4'FC'FLAC+K2) );
             (SGN>0! F-:= FLUC ! F+:= FLUC)
      'FI';
 
      EXCEPTION : 'SKIP' ;
                  'FLUX'(F)
   'END';#$E#
         #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' POSINT = ('STATE' Q0,Q1,
                    'REF''FLUXAC' FAC0,FAC1)'FLUX':
   #$B POSINT  B$#
   'BEGIN' 'FLUX' F:= OSHER(Q0,Q1,FAC0,FAC1);
           (FAC1 'IS' 'NIL'
                    ! F:=ZFLUX(Q1,'NIL')-F
                    ! 'FLUXAC' FLAC;
                      F:=ZFLUX(Q1,FLAC)-F;
                      FAC1:=FLAC-FAC1     );
           (FAC0 'ISNT''NIL'
                    ! FAC0:= -FAC0);
           'FLUX'(F)
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' NEGINT = ('STATE' Q0,Q1,
                    'REF' 'FLUXAC' FAC0,FAC1)'FLUX':
   #$B NEGINT  B$#
   'BEGIN' 'FLUX' F:= OSHER(Q0,Q1,FAC0,FAC1);
           (FAC0  'IS' 'NIL'
                    ! F-:= ZFLUX(Q0,'NIL')
                    ! 'FLUXAC' FLAC;
                      F-:= ZFLUX(Q0, FLAC);
                      FAC0-:=FLAC);
           'FLUX'(F)
   'END'; #$E#
          #E$#
 
 # BOUNDARY STATES                #                   'PR' EJECT 'PR'
 
 'PROC'('WALL','REAL','REAL','STATE','REF''STATAC')'STATE'
        NOORD, ZUID, OOST, WEST;
 
   #$H#
 'PROC' FIXED UVCZ    = ('WALL' W, 'REAL' U,V,C,Z,'STATE' QIN,
                         'REF''STATAC' QAC)'STATE':
   #$B BC.FIX  B$#
 (      (QAC 'ISNT''NIL' ! QAC:=NULFLUXAC); 'STATE'(U,V,C,Z)
 ); #$E#
    #E$#
 
   #$H#
 'PROC' SUB OUT PRESS = ('WALL' W, 'BOOL' LEFT, 'REAL' P1, 'STATE' QIN,
                         'REF''STATAC' QAC)'STATE':
   #$B BC.RSO  B$#
 'BEGIN' # NOT CHECKED: LEFT='TRUE' #
     'REAL' C =C'OF'W, S =S'OF'W, C0=C3'OF'QIN, Z=C4'OF'QIN;
     'REAL' U0=C*C1'OF'QIN+S*C2'OF'QIN,V=-S*C1'OF'QIN+C*C2'OF'QIN;
     'REAL' R1,C1,U1, STOGM:= ( LEFT ! -TOGAMM1 ! TOGAMM1 );
            R1:=EXP((LN(P1)-Z)*OGAM);
            C1:=SQRT(GAM*P1/R1);
            U1:=U0+STOGM*(C0-C1);
     'IF'   QAC 'ISNT''NIL'
     'THEN' 'REAL' DUM=C1*OGAM;
            QAC:=( 1 ,0 ,C*STOGM ,-C*OGAMM1*DUM ,
                   0 ,1 ,S*STOGM ,-S*OGAMM1*DUM ,
                   0 ,0 ,0       ,      0.5*DUM ,
                   0 ,0 ,0       ,        1     )
     'FI';  'STATE'(C*U1-S*V,S*U1+C*V,C1,Z)
   'END'; #$E#
          #E$#
 
   #$H#
 'PROC' WALL          = ('WALL' W, 'BOOL' LEFT, 'REAL' NOR,'STATE' QIN,
                         'REF''STATAC' QAC)'STATE':
   #$B BC.LW   B$#
 'BEGIN'
     'REAL' C   = C'OF'W, S= S'OF'W, C1= C3'OF'QIN, Z = C4'OF'QIN;
     'REAL' U1 =C*C1'OF'QIN+S*C2'OF'QIN,V= -S*C1'OF'QIN+C*C2'OF'QIN,
            SHG= (LEFT!-HGAMM1!+HGAMM1);
     'REAL' U0:=NOR*V,C0;
            C0:=C1-SHG*(U0-U1);
     'IF'   QAC 'ISNT''NIL'
     'THEN' QAC:=( S*( -C*NOR+S)   ,   C*(  C*NOR-S ) ,0   ,0 ,
                   S*( -S*NOR-C)   ,   C*(  S*NOR+C ) ,0   ,0 ,
                   (C+S*NOR)*SHG   ,   (S-C*NOR)*SHG  ,1   ,0 ,
                             0     ,           0      ,0   ,1 )
     'FI';  'STATE'(C*U0-S*V,S*U0+C*V,C0,Z)
   'END'; #$E#
          #E$#
 
                                                      'PR' EJECT 'PR'
   #$H#
 'PROC' SUB IN UVZ = ('WALL' W,'BOOL' LEFT,'REAL' UX,VY,Z0,'STATE' QIN,
                      'REF''STATAC' QAC)'STATE':
   #$B BC.LSI  B$#
 'BEGIN' # NOT CHECKED: LEFT='FALSE' #
         'REAL' C =C'OF'W, S =S'OF'W, C1=C3'OF'QIN,  Z1=C4'OF'QIN;
         'REAL' U0=C*UX    +S*VY    , V0=-S*UX    +C*VY    ,
                U1=C*C1'OF'QIN+S*C2'OF'QIN,V=-S*C1'OF'QIN+C*C2'OF'QIN,
                SHG= (LEFT !HGAMM1 !-HGAMM1);
         'REAL' C0, DD, D1, D2;
       ( SOS>0  #I.E. OSHERS SCHEME #
       !'REAL' P1,R1,RH,CH;
               REVEAL(C1,Z1,P1,R1);
               RH:=EXP(OGAM*(LN(P1)-Z0));
               CH:=SQRT(GAM*P1/RH);
               C0:=CH+SHG*(U0-U1);
               DD:= R1*C1/(RH*CH);
               D1:= -SHG;
               D2:= -CH*OTGAM
       !'REAL' PH,RH,CH,R0;
               CH:= C1+SHG*(U0-U1);
               REVEAL(CH,Z1,PH,RH);
               R0:= EXP(OGAM*(LN(PH)-Z0));
               C0:= SQRT(GAM*PH/R0);
               DD:= RH*CH/(R0*C0);
               D1:= -SHG*DD;
               D2:= -C0*OTGAM
       );
       ( QAC 'ISNT''NIL'
       ! QAC:= ( 0  ,  0  , 0  ,0  ,
                 0  ,  0  , 0  ,0  ,
                D1*C, D1*S, DD ,D2 ,
                 0  ,  0  , 0  ,0  )
       );'STATE'(UX,VY,C0,Z0)
   'END'; #$E#
          #E$#
 
 # EXTEND FIELD WITH WALL STATES  #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' ADD WALL STATES = ('BOOL' SECOND, 'FIELD' QF)'FIELD':
   #$B ADDW    B$#
   'BEGIN' 'REF'[,]'STATE' Q   = Q'OF'QF;
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           STATE Z ( QF );
           [IL-1:IU+1,JL-1:JU+1]'STATE' S;
                S[IL:IU,JL:JU]:= Q[IL:IU,JL:JU];
           'STATE' QQ;
           'FIELD' NEW = 'FIELD'((DDX'OF'QF, DDY'OF'QF,
                                  'HEAP''INT':=TYPE'OF'QF, S));
           'BOOL' PLUS = 'TRUE';
           'REAL' XX:= (IL-1.5)*DDX;
           [JL-1:JU] 'WALL'  HORI,VERI ;
           [JL-1:JU] 'POINT' COORD     ;
 
           ( 'REF'[]'STATE' QI = Q[IL,], QIP = Q[IL+1,];
             WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI);
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QQ:= QI[J]; (SECOND! QQ+:= (QQ-QIP[J])/2.0);
                  S[IL-1,J]:=NOORD(HORI[J],XXN,YY+:=DDY,QQ,'NIL')
             'OD'
           );
 
           'FOR' I 'FROM' IL 'TO' IU
           'DO' 'REF'[]'STATE'  QI = Q[I,];
                WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI);
                  QQ:= QI[JL]; (SECOND! QQ+:= (QQ-QI[JL+1])/2.0);
                  S[I,JL-1]:=WEST(VERI[JL-1],XX+:=DDX,YYW,QQ,'NIL');
                  QQ:= QI[JU]; (SECOND! QQ+:= (QQ-QI[JU-1])/2.0);
                  S[I,JU+1]:=OOST(VERI[JU  ],XX      ,YYO,QQ,'NIL')
           'OD';
 
           ( 'REF'[]'STATE' QI = Q[IU,], QIM = Q[IU-1,];
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QQ:= QI[J]; (SECOND! QQ+:= (QQ-QIM[J])/2.0);
                  S[IU+1,J]:= ZUID(HORI[J],XXZ,YY+:=DDY,QQ,'NIL')
             'OD'
           );
 
           S[IL-1,JL-1] :=  S[IL-1,JL] + S[IL,JL-1] - S[IL,JL];
           S[IL-1,JU+1] :=  S[IL-1,JU] + S[IL,JU+1] - S[IL,JU];
           S[IU+1,JL-1] :=  S[IU+1,JL] + S[IU,JL-1] - S[IU,JL];
           S[IU+1,JU+1] :=  S[IU+1,JU] + S[IU,JU+1] - S[IU,JU];
           NEW
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
 #$H#
 'PROC' ADD WALL = ('BOOL' SECOND, 'FIELD' F)'FIELD':
 #$B ADDWSH B$#
 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'F;
         'INT' L1 = 1'LWB' Q, U1 = 1'UPB'Q, L2 = 2'LWB'Q, U2 = 2'UPB'Q;
         'INT' L1S = L1-1, U1S = U1+1, L2S = L2-1, U2S = U2+1;
         'HEAP'[L1S:U1S,L2S:U2S]'STATE' S;
         S[L1:U1,L2:U2]:= Q[L1:U1,L2:U2];
         'IF'   SECOND 'AND' (U1-L1+1>=3) 'AND' (U2-L2+1>=3)
         'THEN' 'FOR' J 'FROM' L2S+1 'TO' U2S-1
                'DO'  S[L1S,J]:=
                      3.0*Q[L1,J]-3.0*Q[L1+1,J]+Q[L1+2,J];
                      S[U1S,J]:=
                      3.0*Q[U1,J]-3.0*Q[U1-1,J]+Q[U1-2,J]
                'OD'  ;
                'FOR' I 'FROM' L1S+1 'TO' U1S-1
                'DO'  S[I,L2S]:=
                      3.0*Q[I,L2]-3.0*Q[I,L2+1]+Q[I,L2+2];
                      S[I,U2S]:=
                      3.0*Q[I,U2]-3.0*Q[I,U2-1]+Q[I,U2-2]
                'OD'
         'ELSE' 'FOR' J 'FROM' L2S+1 'TO' U2S-1
                'DO'  S[L1S,J]:= 2.0*Q[L1,J]-Q[L1+1,J];
                      S[U1S,J]:= 2.0*Q[U1,J]-Q[U1-1,J]
                'OD'  ;
                'FOR' I 'FROM' L1S+1 'TO' U1S-1
                'DO'  S[I,L2S]:= 2.0*Q[I,L2]-Q[I,L2+1];
                      S[I,U2S]:= 2.0*Q[I,U2]-Q[I,U2-1]
                'OD'
         'FI'   ;
         S[L1S,L2S]:= S[L1S  ,L2S+1] + S[L1S+1,L2S  ] - S[L1S+1,L2S+1];
         S[L1S,U2S]:= S[L1S+1,U2S  ] + S[L1S  ,U2S-1] - S[L1S+1,U2S-1];
         S[U1S,L2S]:= S[U1S  ,L2S+1] + S[U1S-1,L2S  ] - S[U1S-1,L2S+1];
         S[U1S,U2S]:= S[U1S-1,U2S  ] + S[U1S  ,U2S-1] - S[U1S-1,U2S-1];
         'FIELD'(DDX'OF'F,DDY'OF'F,'HEAP''INT':=TYPE'OF'F,S)
 'END' ; #$E#
         #E$#
 
 # WALLS                          #                   'PR' EJECT 'PR'
 
 'PROC' WALLS := ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY,
                  'REF'[]'POINT' COORD,
                  'REF'[]'WALL' VERI,HORI)'VOID': (ERROR;'SKIP');
 
   #$H#
 'PROC' WALLSG = ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY,
                                  'REF'[]'POINT' COORD,
                                  'REF'[]'WALL' VERI,HORI)'VOID':
   #$B WALLSG  B$#
 'BEGIN' 'BOOL' VER = VERI 'ISNT' 'NIL';
         'INT'  JLM = 'LWB' HORI, JU = 'UPB' HORI;
         'REAL' DX,DY,DS;
 
         'POINT' PT,PB,PO;
         'REF''REAL' XT= X'OF'PT, XB= X'OF'PB, XO= X'OF'PO,
                     YT= Y'OF'PT, YB= Y'OF'PB, YO= Y'OF'PO;
         # ALLEEN AFBEELDINGEN MET POSITIEVE ORIENTATIE ZYN TOEGESTAAN #
 
         'FOR' J 'FROM' JLM 'TO' JU
         'DO' 'IF'   VER
              'THEN' PT:= COORD[J];
                     COORD[J]:= PB:= MAPPING(XX,J*DDY);
                     (FORWARD ! DX:= XT-XB; DY:= YT-YB
                              ! DX:= XB-XT; DY:= YB-YT );
                     DS:= SQRT(DX*DX+DY*DY);
                     VERI[J]:= 'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS))
              'ELSE' COORD[J]:= PB:= MAPPING(XX,J*DDY)
              'FI';
              'IF' J > JLM
              'THEN' DX:= XB-XO; DY:= YB-YO;
                     DS:= SQRT(DX*DX+DY*DY);
                     HORI[J]:= 'WALL'( DY/DS,-DX/DS,(PLUS!DS!-DS))
              'FI';  PO:= PB
         'OD'
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
 'PROC' WALLSR = ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY,
                                  'REF'[]'POINT' COORD,
                                  'REF'[]'WALL' VERT,HORI)'VOID':
   #$B WALLSR  B$#
 'IF' HORI 'IS' 'NIL'
 'THEN' 'FOR' J 'FROM' 'LWB' COORD 'TO' 'UPB' COORD
        'DO' COORD[J]:= MAPPING(XX,J*DDY) 'OD'
 'ELSE' 'BOOL'      VER = VERT 'ISNT' 'NIL';
        'INT'       JLM = 'LWB' HORI, JU = 'UPB' HORI;
        'POINT' PB;
        'REF''REAL' XB = X'OF'PB, YB = Y'OF'PB;
        'REAL' XT,YO, DX,DY,DS;
        'WALL' VERTI;
        XT:= X'OF'COORD[JLM];
        COORD[JLM]:= PB:= MAPPING(XX, JLM*DDY);
 
        'IF' VER
        'THEN' DX:= (FORWARD! XT-XB ! XB-XT);
               DS:= 'ABS'DX;
               VERT[JLM]:= VERTI:= 'WALL'(0,-'SIGN'DX,(PLUS!DS!-DS))
        'FI';
 
        'FOR' J 'FROM' JLM+1 'TO' JU
        'DO' YO:= YB;
             ( VER ! VERT[J]:= VERTI);
             COORD[J]:= PB:= MAPPING(XX,J*DDY);
             DY:= YB-YO; DS:='ABS'DY;
             HORI [J]:= 'WALL'('SIGN'DY, 0,(PLUS!DS!-DS))
        'OD'
   'FI'; #$E#
         #E$#
 
 # FLUXT                          #                   'PR' EJECT 'PR'
 
 'PROC' FLUXT:= ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX':
 (ERROR;'SKIP');
 
   #$H#
 'PROC' FLUXTR = ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX':
   #$B FLUXTR  B$#
 'IF'   'REAL' DS = DS'OF'W, ADS = 'ABS'DS'OF'W;
        C'OF'W > 0.999
 'THEN' 'FLUX' N := NFLUX(Q0,Q1,NF0,NF1);
        (NF0'ISNT''NIL'! NF0 *:= ADS );
        (NF1'ISNT''NIL'! NF1 *:= ADS );
        N*:=DS
 'ELIF' S'OF'W > 0.999
 'THEN'  'REAL' U0=C1'OF'Q0, V0=C2'OF'Q0, U1=C1'OF'Q1, V1=C2'OF'Q1;
         'FLUX' NF = NFLUX (
                (V0, -U0, C3'OF'Q0, C4'OF'Q0),
                (V1, -U1, C3'OF'Q1, C4'OF'Q1), NF0, NF1);
         'FLUX' NT;
         ( NF0 'ISNT' 'NIL'
         !  NT:= -(2'FC'NF0) ; CASS(NF0,2,1'FC'NF0) ; CASS(NF0,1,NT);
            NT:= -(3'FR'NF0) ; RASS(NF0,3,2'FR'NF0) ; RASS(NF0,2,NT);
            NF0 *:= ADS
         );
         ( NF1 'ISNT' 'NIL'
         !  NT:= -(2'FC'NF1) ; CASS(NF1,2,1'FC'NF1) ; CASS(NF1,1,NT);
            NT:= -(3'FR'NF1) ; RASS(NF1,3,2'FR'NF1) ; RASS(NF1,2,NT);
            NF1 *:= ADS
         );
         'FLUX' ( DS*C1'OF'NF,-DS*C3'OF'NF, DS*C2'OF'NF, DS*C4'OF'NF )
  'ELSE' ERROR; 'SKIP'
   'FI' ; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
 'PROC' FLUXTG = ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX':
   #$B FLUXTG  B$#
 'BEGIN' 'REAL' C  =  C'OF'W, S   =       S'OF'W,
                DS = DS'OF'W, ADS = 'ABS'DS'OF'W,
                U0=C1'OF'Q0, V0=C2'OF'Q0, U1=C1'OF'Q1, V1=C2'OF'Q1;
         'FLUX' NF = NFLUX (
                (C*U0+S*V0, -S*U0+C*V0, C3'OF'Q0, C4'OF'Q0),
                (C*U1+S*V1, -S*U1+C*V1, C3'OF'Q1, C4'OF'Q1),NF0,NF1);
         'STATE' MT,NT;
         ( NF0 'ISNT' 'NIL'
         !  MT := 1'FC'NF0; NT := 2'FC'NF0;
            CASS(NF0,1,C*MT-S*NT); CASS(NF0,2,S*MT+C*NT);
            MT := 2'FR'NF0; NT := 3'FR'NF0;
            RASS(NF0,2,C*MT-S*NT); RASS(NF0,3,S*MT+C*NT);
            NF0 *:= ADS
         );
         ( NF1 'ISNT' 'NIL'
         !  MT := 1'FC'NF1; NT := 2'FC'NF1;
            CASS(NF1,1,C*MT-S*NT);CASS(NF1,2,S*MT+C*NT);
            MT := 2'FR'NF1; NT := 3'FR'NF1;
            RASS(NF1,2,C*MT-S*NT);RASS(NF1,3,S*MT+C*NT);
            NF1 *:= ADS
         );
         'REAL' MS = DS*C2'OF'NF, NS = DS*C3'OF'NF;
         'FLUX' ( DS*C1'OF'NF, C*MS-S*NS, S*MS+C*NS, DS*C4'OF'NF )
   'END'; #$E#
          #E$#
 
 # RESIDUAL                       #                   'PR' EJECT 'PR'
 
   'BOOL' CENTRAL2:= 'FALSE', ORDER2:= 'FALSE';
 
   #$H#
   'OP' 'MINRESIDU' = ('FIELD' F,Q)'FIELD':
   #$B MINRES  B$#
   'BEGIN' 'FIELD' FF = 'COPY'F ; RESIDU('FALSE',Q,FF); FF
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'PLUSRESIDU' = ('FIELD' F,Q)'FIELD':
   #$B PLUSRES B$#
   'BEGIN' 'FIELD' FF = 'COPY'F ; RESIDU('TRUE',Q,FF); FF
   'END' ; #$E#
           #E$#
 
   #$H#
   'OP' 'RESIDU'    = ('FIELD' Q)'FIELD':
   #$B RESIDU  B$#
   'BEGIN' 'FIELD' F = 'NULRHS'Q ; RESIDU('TRUE',Q,F); F
   'END' ; #$E#
           #E$#
 
   'PROC' RESIDU  :=('BOOL'PLUS,'FIELD'QF,FF)'FIELD': 'SKIP';
   'PROC' RESIDUHO:=('BOOL'PLUS,'FIELD'QF,FF)'FIELD':RESIDU(PLUS,QF,FF);
 
   #$H#
   'PROC' RESIDU1 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD':
   #$B RESID1  B$#
   'BEGIN' 'REF'[,]'STATE' Q   = Q'OF'QF, RHS = Q'OF'RF;
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           'REAL' XX:= (IL-1.5)*DDX;
           [JL-1:JU] 'WALL' HORI,VERI;
           [JL-1:JU]'POINT' COORD;
 
           'STATE' QOUT;               STATE Z ( QF );
           'FLUX'  NF;
 
           ( 'REF'[]  'STATE'   QI =   Q[IL,],
                              RHSI = RHS[IL,];
             WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI);
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QOUT  := NOORD(HORI[J],XXN,YY+:=DDY,QI[J],'NIL');
                  NF    := FLUXT(HORI[J],QOUT,QI[J],'NIL','NIL');
                  RHSI[J]-:= NF
             'OD'
           );
 
           'FOR' I 'FROM' IL 'TO' IU
           'DO' 'REF'[]'STATE'     QI =   Q[I,],
                                 RHSI = RHS[I,];
                WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI);
 
                QOUT  := WEST(VERI[JL-1],XX+:=DDX,YYW,QI[JL],'NIL');
                NF    := FLUXT(VERI[JL-1],QOUT,QI[JL],'NIL','NIL');
                RHSI[JL]-:= NF;
 
                'FOR' J 'FROM' JL 'TO' JU-1
                'DO' # VERT. WALLS #
                     NF:= FLUXT(VERI[J],QI[J],QI[J+1],'NIL','NIL');
                     RHSI[J  ]+:= NF;
                     RHSI[J+1]-:= NF
                'OD';
 
                QOUT  := OOST(VERI[JU],  XX,YYO,QI[JU],'NIL');
                NF    := FLUXT(VERI[JU],QI[JU],QOUT,'NIL','NIL');
                RHSI [JU]+:= NF;
 
           'IF' I < IU
           'THEN' 'REF'[]'STATE'     QIP =   Q[I+1,],
                                   RHSIP = RHS[I+1,];
                  'FOR' J 'FROM' JL 'TO' JU
                  'DO' #HORIZ. WALLS #
                       NF:= FLUXT(HORI[J],QI[J],QIP[J],'NIL','NIL');
                       RHSI [J ]+:= NF;
                       RHSIP[J ]-:= NF
                  'OD'
           'FI'
           'OD';
 
           ( 'REF'[]  'STATE'   QI =   Q[IU,],
                              RHSI = RHS[IU,];
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QOUT  := ZUID(HORI[J],XXZ,YY+:=DDY,QI[J],'NIL');
                  NF    := FLUXT(HORI[J],QI[J],QOUT,'NIL','NIL');
                  RHSI[J]+:= NF
             'OD'
           );'SKIP'
   'END'; #$E#
          #E$#
 
 # SECOND ORDER RESIDUAL (I)      #                   'PR' EJECT 'PR'
 
   'PROC' RESIDU2:= ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': ( ERROR; 'SKIP');
 
   #$H#
   'PROC' CENTRALX= ('BOOL' P, 'FIELD' Q,F) 'FIELD':
   #$B CNTR.X  B$#
   'BEGIN' NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' D0,D1)'FLUX':
                   ZFLUX(0.5*(Q0+Q1), 'NIL');
 
           RESIDU1 (P,Q,F);
           NFLUX:= OSHER; 'SKIP'
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' CENTRALY= ('BOOL' P, 'FIELD' Q,F) 'FIELD':
   #$B CNTR.Y  B$#
   'BEGIN' NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' D0,D1)'FLUX':
                    0.5*(ZFLUX(Q0,'NIL') + ZFLUX(Q1,'NIL'));
 
           RESIDU1 (P,Q,F);
           NFLUX:= OSHER; 'SKIP'
   'END'; #$E#
          #E$#
 
 # RESIDUAL SMOOTHERS             #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' SMOOTHBOX = ('FIELD' A)'FIELD':
   #$B SMOBX   B$#
   'BEGIN' 'INT' TYPE = TYPE'OF'A; TYPE'OF'A:= 1;
           'FIELD' B = 'PROLON' 'RESTRICT'A;
           TYPE'OF'A:= TYPE; Q'OF'A := Q'OF'B; A
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' SMOOTH5   = ('FIELD' A)'FIELD':
   #$B SMO5    B$#
   'BEGIN' 'REF'[,]'STATE' Q=Q'OF'A;
     'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
     [JL:JU]'STATE'QIM:=Q[IL,JL:JU] ; 'STATE' QIJM,QIJOLD ;
     'FOR' I 'FROM' IL 'TO' IU-1 'DO'
       'REF'[]'STATE'QI=Q[I,],QIP=Q[I+1,] ; QIJM:=QI[JL] ;
       'FOR' J 'FROM' JL 'TO' JU-1 'DO'
         'REF''STATE' QIJ=QI[J]; QIJOLD:=QIJ ;
      (((((QIJ*:=4.0)+:=QIJM)+:=QIM[J])+:=QI[J+1])+:=QIP[J])*:=0.125;
         QIJM := QIJOLD ; QIM[J]:=QIJOLD
       'OD' ; QIJOLD:=QI[JU] ;
       ((((QI[JU]*:=5.0)+:=QIJM)+:=QIM[JU])+:=QIP[JU])*:=0.125 ;
       QIM[JU]:=QIJOLD
     'OD' ; 'REF'[]'STATE' QI=Q[IU,] ; QIJM:=QI[JL] ;
     'FOR' J 'FROM' JL 'TO' JU-1 'DO' 'REF''STATE' QIJ=QI[J] ;
       QIJOLD := QIJ ; ((((QIJ*:=5.0)+:=QIJM)+:=QIM[J])+:=
       QI[J+1])*:=0.125 ; QIJM := QIJOLD
     'OD' ; ((( QI[JU]*:=6.0)+:=QIJM)+:=QIM[JU])*:=0.125 ; A
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' JACSMOOTH  = ('REAL' DAMPING, 'FIELD' QF, RHS)'FIELD':
   #$B JACSMO  B$#
   'BEGIN' 'FIELD' JJJ = 'COPY' QF;
           'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS,
                           JJ= Q'OF'JJJ;
                               STATE Z ( QF );
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF;
 
      [JL-1:JU]   'WALL'  HORIN,HORIZ,VERI ;
      [JL-1:JU]   'POINT' COORD            ;
      [JL  :JU]   'STATE' QIM,QIO          ;
      'STATE'  UC, QOUT;  'STATAC' QAC     ;
      'FLUX'   NF,F                        ;
      'FLUXAC' NFAC0,NFAC1,FAC             ;
 
         WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIN);
         'FOR' I 'FROM' IL 'TO' IU
         'DO''REAL'XX=(I-0.5)*DDX;
             WALLS( 'TRUE','TRUE', I  *DDX,DDY,COORD, VERI,HORIZ);
             QIO:= Q[I,];
 
         'FOR' J 'FROM' JL 'TO' JU
         'DO''REAL'YY=(J-0.5)*DDY;
             'REF''STATE' QIJ=Q[I,J]; UC:=QIJ;
 
               F   :=(I=IU ! QOUT :=  ZUID(HORIZ[J],XXZ,YY,UC,QAC);
                             FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL'));
               FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC);
                             FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(HORIN[J],QIM[J],UC,'NIL',NFAC1));
               FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 );
 
               F  +:=(J=JU ! QOUT :=  OOST(VERI[J],XX,YYO,UC,QAC);
                             FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL'));
               FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(J=JL ! QOUT :=  WEST(VERI[J-1],XX,YYW,UC,QAC);
                             FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(VERI[J-1],QIO[J-1],UC,'NIL',NFAC1));
               FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 );
 
                  JJ[I,J] := (FF[I,J]-F) /FAC;
                  JJ[I,J]*:= DAMPING
 
            'OD'; QIM:= QIO; HORIN:= HORIZ
            'OD'; JJJ
   'END'; #$E#
          #E$#
 
 # SECOND ORDER RESIDUAL (II)     #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' SUPERBOX = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD':
   #$B SUPERB  B$#
   'IF'   'REF'[,]'STATE' QQ  = Q'OF'QF, RHS = Q'OF'RF;
          'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
          'INT'  IL=1'LWB'QQ,IU=1'UPB'QQ,JL=2'LWB'QQ,JU=2'UPB'QQ;
          'NOT'( 'ODD'(IU-IL) 'AND' 'ODD'(JU-JL) )
   'THEN' RESIDU1(PLUS,QF,RF)
   'ELSE' [0:2,JL-1:JU] 'WALL' HOR;
          [1:2,JL-1:JU] 'WALL' VER;
          [JL-1:JU]'POINT' COORD  ;
          [JL  :JU]'STATE' OUTW   ;
          [1:2]    'STATE' LEFT   ;
          [1:2,1:2]'STATE' NS,WE  ;
                   'STATE' QOUT   ;    STATE Z ( QF );
          'FLUX'  NF;
 
   'INT' II,JJ,JM;
   'STATE' A,B,C,D,P,Q,R,S;
 
   'REAL' XX:= (IL-1.5)*DDX,
          YY:= (JL-1.5)*DDY;
 
   WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HOR[0,]);
 
   'FOR' I 'FROM' IL 'BY' 2 'TO' IU
   'DO'
 
   WALLS('TRUE',PLUS, I   *DDX,DDY,COORD,VER[1,],HOR[1,]);
   WALLS('TRUE',PLUS,(I+1)*DDX,DDY,COORD,VER[2,],HOR[2,]);
 
   'FOR' J 'FROM' JL 'BY' 2 'TO' JU
   'DO'  JM:= J-1;
 
         # INTERPOLATIE GEDEELTE #
         P:= QQ[ I  ,J  ]  ;  R:= QQ[ I  ,J+1]  ;
         Q:= QQ[ I+1,J  ]  ;  S:= QQ[ I+1,J+1]  ;
         A:= 0.25*(P+Q+R+S);  C:= 0.25*(R+S-P-Q);
         B:= 0.25*(Q+S-P-R);  D:= 0.25*(P+S-R-Q);
         # Q = A + B.X + C.Y + D.XY #
        'FOR' K 'TO' 2    # K=1: N,W; K=2: S,E #
        'DO' 'REAL' Z = ( K !  -2.0 , 2.0 );
             P:=A+Z*B ; Q:=C+Z*D ; NS[K,1]:=P-Q ; NS[K,2]:=P+Q ;
             P:=A+Z*C ; Q:=B+Z*D ; WE[K,1]:=P-Q ; WE[K,2]:=P+Q
        'OD'; # NAGEGAAN 850315 #
 
        'FOR' L 'TO' 2
        'DO' JJ:= J+L-1;
           ( I=IL
           !OUTW[JJ] :=NOORD(HOR[0,JJ],XXN,YY+:=DDY,NS[1,L]  ,'NIL');
            NF       :=FLUXT(HOR[0,JJ],OUTW[JJ],NS[1,L],'NIL','NIL');
                       RHS[I  ,JJ] -:= NF
           !NF       :=FLUXT(HOR[0,JJ],OUTW[JJ],NS[1,L],'NIL','NIL');
                       RHS[I-1,JJ] +:= NF;
                       RHS[I  ,JJ] -:= NF
           );
           OUTW[JJ]:= NS[2,L]
        'OD';
 
        'FOR' K 'TO' 2
        'DO' II:= I+K-1;
           (J=JL
           !LEFT[K] :=WEST (VER[K,JM],XX+:=DDX,YYW,WE[1,K] ,'NIL');
            NF      :=FLUXT(VER[K,JM],LEFT[K],WE[1,K],'NIL','NIL');
                      RHS[II,J ] -:= NF
           !NF      :=FLUXT(VER[K,JM],LEFT[K],WE[1,K],'NIL','NIL');
                      RHS[II,JM] +:= NF;
                      RHS[II,J ] -:= NF
           );
           LEFT[K]:= WE[2,K]
        'OD';
 
        NF:= FLUXT( VER[1,J  ],A-B,A-B,'NIL','NIL' );
             RHS[I  ,J  ]+:= NF; RHS[I  ,J+1] -:= NF;
        NF:= FLUXT( VER[2,J  ],A+B,A+B,'NIL','NIL' );
             RHS[I+1,J  ]+:= NF; RHS[I+1,J+1] -:= NF;
        NF:= FLUXT( HOR[1,J  ],A-C,A-C,'NIL','NIL' );
             RHS[I  ,J  ]+:= NF; RHS[I+1,J  ] -:= NF;
        NF:= FLUXT( HOR[1,J+1],A+C,A+C,'NIL','NIL');
             RHS[I  ,J+1]+:= NF; RHS[I+1,J+1] -:= NF
 
   'OD'; HOR[0,]:= HOR[2,];
 
        'FOR' K 'TO' 2
        'DO' II:= I+K-1;
            WE[1,K]:= OOST (VER[K,JU],XX     ,YYO,WE[2,K],'NIL');
            NF     := FLUXT(VER[K,JU],LEFT[K],WE[1,K],'NIL','NIL');
                      RHS[II,JU]+:= NF
        'OD'
   'OD';
 
   YY:= (JL-1.5)*DDY;
   'FOR' J 'FROM' JL 'TO' JU
   'DO' QOUT  := ZUID (HOR[2,J],XXZ,YY+:=DDY,OUTW[J],'NIL');
        NF    := FLUXT(HOR[2,J],OUTW[J],QOUT,'NIL','NIL');
        RHS[IU,J]+:= NF
   'OD';'SKIP'
 
   'FI' ; #$E#
          #E$#
 
 # SECOND ORDER RESIDUAL (III)    #                   'PR' EJECT 'PR'
 
   'REAL' KAPPA:= 1/3;
 
   #$H#
   'PROC' RESIDUVL = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD':
   #$B RESIVL  B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, RHS = Q'OF'RF;
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           'REAL' XX:= (IL-1.5)*DDX, YY:= (JL-1.5)*DDY;
           [JL-1:JU]'WALL' HORIN,HORIZ,VERI;
           [JL-1:JU]'POINT' COORD;
 
           'REAL' KP = (1+KAPPA)/4, KM = (1-KAPPA)/4;
           [JL:JU]'STATE' QNO,DIM;               STATE Z ( QF );
           'STATE' DJM,DIP,DJP, QNI,QZI, QWO,QWI,QOI, QIJ;
           'FLUX'  NF;
   #
    KAPPA = -1 : ZUIVER UPWIND,
    KAPPA = +1 : CENTRALE DIFFERENTIES
                 (MET DISCONTINUE AANSLUITING OP DE RANDWAARDEN,
                  VOLGENS DE RIEMANN INVARIANTEN)
 
    BETEKENIS VAN DE VARIABELEN:
 
                           ---------> J         ^
                                                !
                                   NO           !
                 !         .-----------------   ! DIM
                 !         !       NI           !
                 !         !                    !
                 !         !                    !
                 !     WO  ! WI    QIJ   OI     V
             I   !         !                    ^
                 V         !                    !
                           !       ZI           !
                           !                    ! DIP
                   <--------------> <-----------.-------->
                          DJM             DJP   !
                                                !
                                                !
                                                V
   #
 
           WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORIN);
 
           'FOR' I 'FROM' IL 'TO' IU
           'DO' 'REF'[]'STATE' QI   =   Q[I,] ,
                               RHSI = RHS[I,] ;
                WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORIZ);
 
                'FOR' J 'FROM' JL 'TO' JU
                'DO' 'FLUX' QIJ = QI[J];
 
                     # VERT. WALLS #
                     DJP:= ( J=JU ! DJM ! QI[J+1] - QIJ);
                     ( J=JL
                     ! QWI:= QIJ -0.5*DJP;
                       QOI:= QIJ +0.5*DJP;
                       QWO:= WEST(VERI[JL-1],XX+:=DDX,YYW,QWI,'NIL');
                       NF := FLUXT(VERI[J-1],QWO,QWI,'NIL','NIL')
                     ! QWI:= QIJ - SCHEME(DJP,DJM);
                       QOI:= QIJ + SCHEME(DJM,DJP);
                       NF := FLUXT(VERI[J-1],QWO,QWI,'NIL','NIL');
                       RHSI[J-1]+:= NF
                     );RHSI[J  ]-:= NF;
                     DJM:= DJP;
                     QWO:= QOI;
 
                     # HORIZ. WALLS #
                     DIP:= ( I=IU ! DIM[J] ! Q[I+1,J] - QIJ);
                     ( I=IL
                     ! QNI:= QIJ -0.5*DIP;
                       QZI:= QIJ +0.5*DIP;
                       QNO[J]:= NOORD(HORIN[J],XXN,YY+:=DDY,QNI,'NIL');
                       NF    := FLUXT(HORIN[J],QNO[J],QNI,'NIL','NIL')
                     ! QNI:= QIJ - SCHEME(DIP   ,DIM[J]);
                       QZI:= QIJ + SCHEME(DIM[J],DIP   );
                       NF := FLUXT(HORIN[J],QNO[J],QNI,'NIL','NIL');
                       RHS[I-1,J]+:= NF
                     );RHSI   [J]-:= NF;
                     DIM[J]:= DIP;
                     QNO[J]:= QZI
                'OD';
 
                QWI := OOST (VERI[JU], XX,YYO, QWO ,'NIL');
                NF  := FLUXT(VERI[JU],QWO,QWI,'NIL','NIL');
                RHSI[JU]+:= NF;
 
                HORIN:= HORIZ
           'OD';
 
           YY:= (JL-1.5)*DDY;
           'FOR' J 'FROM' JL 'TO' JU
           'DO' QNI := ZUID (HORIZ[J],XXZ,YY+:=DDY,QNO[J],'NIL');
                NF  := FLUXT(HORIZ[J],QNO[J],QNI,'NIL','NIL');
                RHS[IU,J]+:= NF
           'OD'; 'SKIP'
   'END'; #$E#
          #E$#
 
 # DIFFERENT VAN LEER TYPE SCHEMES #                  'PR' EJECT 'PR'
 
  'PROC' SCHEME := ('STATE' DM,DP)'STATE':
                   (0.25*(1.0-KAPPA))*DM + (0.25*(1.0+KAPPA))*DP;
 
  #$H#
  'PROC' ONESIDED = ('STATE' DM,DP)'STATE':
  #$B SHM.ONE B$#
  'BEGIN' 0.5*DM
  'END' ; #$E#
          #E$#
 
  #$H#
  'PROC' FROMM    = ('STATE' DM,DP)'STATE':
  #$B SHM.FRM B$#
  'BEGIN' 0.25*(DM+DP)
  'END' ; #$E#
          #E$#
 
  #$H#
  'PROC' UPWIND BIASED = ('STATE' DM,DP)'STATE':
  #$B SHM.BIA B$#
  'BEGIN' (1/3)*DM + (2/3)*DP
  'END' ; #$E#
          #E$#
  #$H#
  'PROC' CENTERED  = ('STATE' DM,DP)'STATE':
  #$B SHM.CEN B$#
  'BEGIN' 0.5*DP
  'END' ; #$E#
          #E$#
 
  #$H#
  'PROC' ALBADA        = ('STATE' DM,DP)'STATE':
  #$B SHM.ALB B$#
  'BEGIN' [1:4]'REAL' U, 'REAL' A,B,R,F,
          []'REAL' M = (C1'OF'DM, C2'OF'DM, C3'OF'DM, C4'OF'DM),
                   P = (C1'OF'DP, C2'OF'DP, C3'OF'DP, C4'OF'DP);
          'FOR' K 'TO' 4
          'DO' A:= P[K]; B:= M[K];
               ( 'ABS'A > 'ABS'B ! R:=B/A; R:=R*R; F:= B+A*R
                                 ! R:=A/B; R:=R*R; F:= A+B*R   );
               U[K]:= F/(2.0 + 2.0*R)
          'OD'; 'STATE' (U[1],U[2],U[3],U[4])
  'END'; #$E#
         #E$#
 
 # FIRST  ORDER STAR-RESIDUAL     #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' STARRES1 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD':
   #$B STAR1   B$#
   'IF'   'REF'[,]'STATE' Q  = Q'OF'QF, R = Q'OF'RF;
          'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
          'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
          'INT'  ILH = (IL+1)'OVER'2, IUH = IU'OVER'2,
                 JLH = (JL+1)'OVER'2, JUH = JU'OVER'2;
          'NOT'(ILH*2=IL+1'AND' IUH*2=IU'AND' JLH*2=JL+1'AND' JUH*2=JU )
   'THEN' RESIDU1(PLUS,QF,RF)
   'ELSE' 'REAL' DX, DY, DS, XX,YY, DDXH:= 2*DDX, DDYH:= 2*DDY;
          'WALL'  CW;
          [JLH-1:JUH]'POINT' COORDN, COORDZ;
          [JLH-1:JUH] 'WALL' HORIN, HORIZ, VERT;
          'FLUX'  NF;
          'STATE' QOUT;                        STATE Z ( QF );
 
   WALLS('TRUE',PLUS,(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIN);
 
   'FOR' IH 'FROM' ILH 'TO' IUH
   'DO' 'INT' I = 2*IH-1; XX:= I*DDX;
 
   COORDN :=                        COORDZ;
   WALLS('TRUE',PLUS,IH*DDXH,DDYH,COORDZ,VERT,HORIZ);
 
   'FOR' JH 'FROM' JLH 'TO' JUH
   'DO' 'INT' J = 2*JH-1, JMH = JH-1;
        'REF''STATE' QA=Q[I  ,J  ], QB=Q[I  ,J+1],
                     QC=Q[I+1,J+1], QD=Q[I+1,J  ];
        'REF''FLUX'  RA=R[I  ,J  ], RB=R[I  ,J+1],
                     RC=R[I+1,J+1], RD=R[I+1,J  ];
        'POINT' MID = MAPPING(XX,YY:=J*DDY);
 
        'PROC' DIAGH = ('POINT' CP)'WALL':
        'BEGIN' DX := X'OF'CP - X'OF'MID;
                DY := Y'OF'CP - Y'OF'MID;
                DS := SQRT(DX*DX+DY*DY);
               'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS) )
        'END';
 
        ( IH=ILH
        ! QOUT :=NOORD(HORIN[JH],XXN,YY,QB,'NIL');
          RB  -:=FLUXT(HORIN[JH],QOUT    ,QB,'NIL','NIL')
        ! NF   :=FLUXT(HORIN[JH],Q[I-1,J],QB,'NIL','NIL');
                 RB       -:= NF;
                 R[I-1,J] +:= NF
        );
 
        (JH=JLH
        ! QOUT :=WEST (VERT[JMH],XX,YYW,QA,'NIL');
          RA  -:=FLUXT(VERT[JMH],QOUT,QA,'NIL','NIL')
        ! NF   :=FLUXT(VERT[JMH],Q[I+1,J-1],QA,'NIL','NIL');
                 RA         -:= NF;
                 R[I+1,J-1] +:= NF
        );
 
        CW:= DIAGH(COORDN[JMH]);
        NF:= FLUXTG(CW,QA,QB,'NIL','NIL'); RA+:= NF; RB-:= NF;
        CW:= DIAGH(COORDN[JH ]);
        NF:= FLUXTG(CW,QB,QC,'NIL','NIL'); RB+:= NF; RC-:= NF;
        CW:= DIAGH(COORDZ[JH ]);
        NF:= FLUXTG(CW,QC,QD,'NIL','NIL'); RC+:= NF; RD-:= NF;
        CW:= DIAGH(COORDZ[JMH]);
        NF:= FLUXTG(CW,QD,QA,'NIL','NIL'); RD+:= NF; RA-:= NF
 
   'OD';HORIN := HORIZ;
 
        QOUT      := OOST (VERT[JUH],XX,YYO,Q[2*IH,JU],'NIL');
        R[I+1,JU]+:= FLUXT(VERT[JUH],Q[I+1,JU],QOUT,'NIL','NIL')
   'OD';
 
   'FOR' JH 'FROM' JLH 'TO' JUH
   'DO' 'INT' J= 2*JH-1; YY:= J*DDY;
        QOUT     := ZUID (HORIZ[JH],XXZ,YY,Q[IU,J],'NIL');
        R[IU,J] +:= FLUXT(HORIZ[JH],Q[IU,J],QOUT,'NIL','NIL')
   'OD';'SKIP'
   'FI' ; #$E#
          #E$#
 
 # SECOND ORDER STAR-RESIDUAL     #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' STARRES2 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD':
   #$B STAR2   B$#
   'IF'   'REF'[,]'STATE' Q  = Q'OF'QF, R = Q'OF'RF;
          'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
          'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
          'INT'  ILH = (IL+1)'OVER'2, IUH = IU'OVER'2,
                 JLH = (JL+1)'OVER'2, JUH = JU'OVER'2;
          'NOT'(ILH*2=IL+1'AND' IUH*2=IU'AND' JLH*2=JL+1'AND' JUH*2=JU )
   'THEN' RESIDU1(PLUS,QF,RF)
   'ELSE' 'REAL' DX, DY, DS, XX,YY, DDXH:= 2*DDX, DDYH:= 2*DDY;
          'WALL'  CW;
          [JLH-1:JUH]'POINT' COORDN, COORDZ;
          [JLH-1:JUH] 'WALL' HORIN, HORIZ, VERT;
          'FLUX'  NF;
          'STATE' QIA,QIB,QIC,QID,QUO, QDX,QDY, QVER, QOUT;
          [JLH:JUH]'STATE' QHOR;
                                                     STATE Z ( QF );
 
   WALLS('TRUE',PLUS,(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIN);
 
   'FOR' IH 'FROM' ILH 'TO' IUH
   'DO' 'INT' I = 2*IH-1; XX:= I*DDX;
 
   COORDN :=                        COORDZ;
   WALLS('TRUE',PLUS,IH*DDXH,DDYH,COORDZ,VERT,HORIZ);
 
   'FOR' JH 'FROM' JLH 'TO' JUH
   'DO' 'INT' J = 2*JH-1, JMH = JH-1; YY:= J*DDY;
        'REF''STATE' QA=Q[I  ,J  ], QB=Q[I  ,J+1],
                     QC=Q[I+1,J+1], QD=Q[I+1,J  ];
        'REF''FLUX'  RA=R[I  ,J  ], RB=R[I  ,J+1],
                     RC=R[I+1,J+1], RD=R[I+1,J  ];
        'POINT' MID = MAPPING(XX,YY);
        'PROC' DIAGH = ('POINT' CP)'WALL':
        'BEGIN' DX := X'OF'CP - X'OF'MID;
                DY := Y'OF'CP - Y'OF'MID;
                DS := SQRT(DX*DX+DY*DY);
               'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS) )
        'END';
 
        QIA:= QIB:= QIC:= QID:= 0.25*(QA+QB+QC+QD);
        QDX:= 0.75*(QD-QB);
        QDY:= 0.75*(QC-QA);
        QIA-:= QDY; QIB-:= QDX;
        QIC+:= QDY; QID+:= QDX;
 
        ( IH=ILH
        ! QOUT :=NOORD(HORIN[JH],XXN,YY,QIB,'NIL');
          RB  -:=FLUXT(HORIN[JH],QOUT    ,QIB,'NIL','NIL')
        ! NF   :=FLUXT(HORIN[JH],QHOR[JH],QIB,'NIL','NIL');
                 RB       -:= NF;
                 R[I-1,J] +:= NF
        ); QHOR[JH]:= QID;
 
        (JH=JLH
        ! QOUT :=WEST (VERT[JMH],XX,YYW,QIA,'NIL');
          RA  -:=FLUXT(VERT[JMH],QOUT,QIA,'NIL','NIL')
        ! NF   :=FLUXT(VERT[JMH],QVER,QIA,'NIL','NIL');
                 RA         -:= NF;
                 R[I+1,J-1] +:= NF
        ); QVER:=QIC;
 
        NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' F0,F1)'FLUX':
                  ZFLUX(0.5*(Q0+Q1),'NIL');
        CW:= DIAGH(COORDN[JMH]);
        NF:= FLUXTG(CW,QIA,QIB,'NIL','NIL'); RA+:= NF; RB-:= NF;
        CW:= DIAGH(COORDN[JH ]);
        NF:= FLUXTG(CW,QIB,QIC,'NIL','NIL'); RB+:= NF; RC-:= NF;
        CW:= DIAGH(COORDZ[JH ]);
        NF:= FLUXTG(CW,QIC,QID,'NIL','NIL'); RC+:= NF; RD-:= NF;
        CW:= DIAGH(COORDZ[JMH]);
        NF:= FLUXTG(CW,QID,QIA,'NIL','NIL'); RD+:= NF; RA-:= NF;
        NFLUX:= OSHER
 
   'OD';HORIN := HORIZ;
 
        QOUT      := OOST (VERT[JUH],XX,YYO,QIC,'NIL');
        R[I+1,JU]+:= FLUXT(VERT[JUH],QVER,QOUT,'NIL','NIL')
   'OD';
 
   'FOR' JH 'FROM' JLH 'TO' JUH
   'DO' 'INT' J= 2*JH-1; YY:= J*DDY;
        QOUT     := ZUID (HORIZ[JH],XXZ,YY,QID,'NIL');
        R[IU,J] +:= FLUXT(HORIZ[JH],QHOR[JH],QOUT,'NIL','NIL')
   'OD';'SKIP'
   'FI' ; #$E#
          #E$#
 
 # SECOND ORDER OPERATORS         #                   'PR' EJECT 'PR'
 
   #$H#
   'OP' 'RESIDU2' = ('FIELD' Q)'FIELD':
   #$B RESD2   B$#
   'BEGIN' 'FIELD' F = 'NULRHS'Q;   RESIDU2('TRUE',Q,F); F
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'CORRECT' = ('FIELD' QL,RL)'FIELD':
   #$B CORREC  B$#
   'BEGIN' RESIDU ( 'TRUE',QL,RL);
           RESIDU2('FALSE',QL,RL);
           RL
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'TAU' = ('FIELD' Q)'FIELD':
   #$B TAU.FL B$#
   'BEGIN' 'RESTRICT''RESIDU'Q'MINRESIDU''RESTRICT'Q
   'END' ; #$E#
           #E$#
 
 # LINEAR SYSTEM CONCTRUCTION     #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' LINSYS = ('BOOL'PLUS,'FIELD' QF,RF,
                    'REF''NETMAT' JACN)'VOID':
   #$B LINSYS  B$#
   'BEGIN' 'REF'[,]'STATE' Q   = Q'OF'QF, RHS = Q'OF'RF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           JACN:= 'HEAP'[IL:IU]'LINMAT';
           'NETMAT' JAC = JACN;
           'REAL' XX:= (IL-1.5)*DDX;
           [JL-1:JU] 'WALL' HORI,VERI;
           [JL-1:JU]'POINT' COORD;
 
           'STATE'  QOUT;               STATE Z ( QF );
           'STATAC' AC;
           'FLUX'   NF;
           'FLUXAC' NFAC0,NFAC1;
 
           ( 'REF'[]  'STATE'   QI =   Q[IL,],
                              RHSI = RHS[IL,];
             'REF'[,]'FLUXAC' JACI = GEN(JAC[IL],JL,JU,-2,+2);
                                     'NUL' JACI;
 
             WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI);
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QOUT  := NOORD(HORI[J],XXN,YY+:=DDY,QI[J],AC);
                  NF    := FLUXT(HORI[J],QOUT,QI[J],NFAC0,NFAC1);
                  NFAC1+:= NFAC0*AC;
                  RHSI[  J]-:= NF;
                  JACI[J,0]-:= NFAC1
             'OD'
           );
 
           'FOR' I 'FROM' IL 'TO' IU
           'DO' 'REF'[]'STATE'     QI =   Q[I,],
                                 RHSI = RHS[I,];
                'REF'[,]'FLUXAC' JACI = GET(JAC[I]);
                WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI);
 
                QOUT  := WEST(VERI[JL-1],XX+:=DDX,YYW,QI[JL],AC);
                NF    := FLUXT(VERI[JL-1],QOUT,QI[JL],NFAC0,NFAC1);
                NFAC1+:= NFAC0*AC;
                RHSI[  JL]-:= NF;
                JACI[JL,0]-:= NFAC1;
 
                'FOR' J 'FROM' JL 'TO' JU-1
                'DO' # VERT. WALLS #
                     NF:= FLUXT(VERI[J],QI[J],QI[J+1],NFAC0,NFAC1);
                     RHSI[   J  ]+:= NF;
                     RHSI[   J+1]-:= NF;
                     JACI[J  , 0]+:= NFAC0;
                     JACI[J  , 1]+:= NFAC1;
                     JACI[J+1, 0]-:= NFAC1;
                     JACI[J+1,-1]-:= NFAC0
                'OD';
 
                QOUT  := OOST(VERI[JU],  XX,YYO,QI[JU],AC);
                NF    := FLUXT(VERI[JU],QI[JU],QOUT,NFAC0,NFAC1);
                NFAC0+:= NFAC1*AC;
                RHSI [JU  ]+:= NF;
                JACI [JU,0]+:= NFAC0;
 
           'IF' I < IU
           'THEN' 'REF'[]'STATE'     QIP =   Q[I+1,],
                                   RHSIP = RHS[I+1,];
                  'REF'[,]'FLUXAC' JACIP = GEN(JAC[I+1],JL,JU,-2,2);
                                           'NUL' JACIP;
 
                  'FOR' J 'FROM' JL 'TO' JU
                  'DO' #HORIZ. WALLS #
                       NF:= FLUXT(HORI[J],QI[J],QIP[J],NFAC0,NFAC1);
                       RHSI [J   ]+:= NF;
                       RHSIP[J   ]-:= NF;
                       JACI [J, 0]+:= NFAC0 + DIAG(QI[J]);
                       JACI [J, 2]+:= NFAC1;
                       JACIP[J, 0]-:= NFAC1;
                       JACIP[J,-2]-:= NFAC0
                  'OD';
                   PET('TRUE',JAC[I])
           'FI'
           'OD';
 
           ( 'REF'[]  'STATE'   QI =   Q[IU,],
                              RHSI = RHS[IU,];
             'REF'[,]'FLUXAC' JACI = GET(JAC[IU]);
             'REAL' YY:= (JL-1.5)*DDY;
             'FOR' J 'FROM' JL 'TO' JU
             'DO' QOUT  := ZUID(HORI[J],XXZ,YY+:=DDY,QI[J],AC);
                  NF    := FLUXT(HORI[J],QI[J],QOUT,NFAC0,NFAC1);
                  NFAC0+:= NFAC1*AC;
                  RHSI[J  ]+:= NF;
                  JACI[J,0]+:= NFAC0 + DIAG(QI[J])
             'OD';PET('TRUE',JAC[IU])
           )
   'END'; #$E#
          #E$#
 
 'PROC' DIAG := ('STATE' Q) 'FLUXAC':
 'BEGIN' PRINT((NEWLINE,">>>  DIAG STILL UNDEFINED   <<<",
                NEWLINE,">>> ADVISE : DO DIAG:=DIAG1 <<<",
                NEWLINE)); ERROR ; NULFLUXAC
 'END' ;
 
   #$H#
 'PROC' DIAG1= ('STATE' Q) 'FLUXAC':
   #$B DIAG1   B$#
 'BEGIN' 'REAL' U:=C1'OF'Q, V:=C2'OF'Q, C:=C3'OF'Q, Z:=C4'OF'Q,
                P,R,ROC,ROG,W2,CG;
         REVEAL(C,Z,P,R);     R *:= ODIAGSTEP;
         ROC:= TOGAMM1*R/C;   ROG:= -R*OGAMM1;
         W2 := 0.5*(U*U+V*V); CG := C*C*OGAMM1;
         'FLUXAC'( 0   ,0   ,ROC        ,ROG             ,
                   R   ,0   ,ROC*U      ,ROG*U           ,
                   0   ,R   ,ROC*V      ,ROG*V           ,
                   R*U ,R*V ,ROC*(W2+CG),ROG*(W2+CG*OGAM) )
 'END'; #$E#
        #E$#
 
   #$H#
 'PROC' PLUS DIAG = ('FIELD' QF,DELTAQ,RF) 'VOID':
   #$B PLUSD   B$#
 'BEGIN' 'REF'[,]'STATE' Q =Q'OF'QF, DEL = Q'OF'DELTAQ, RHS =Q'OF'RF;
     'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
 
     'FOR' I 'FROM' IL 'TO' IU
     'DO''REF'[]'STATE' QI=Q[I,], DELI=DEL[I,], RHSI=RHS[I,];
          'FOR' J 'FROM' JL 'TO' JU
          'DO' RHSI[J] +:= DIAG(QI[J])*DELI[J] 'OD'
     'OD'
 
   'END'; #$E#
          #E$#
 
 # GALERKIN OPERATOR CONSTRUCTION #                   'PR' EJECT 'PR'
 
   #$H#
 'OP' 'RAP' = ('NETMAT' F) 'NETMAT':
   #$B RAP     B$#
 'BEGIN'# BY PDZ; INSERTED 850102  #
        'INT' ILF = 'LWB'F, IUF = 'UPB'F;
        'INT' ILC = (ILF+1)'OVER'2, IUC = IUF'OVER'2;
        'NETMAT' C = 'HEAP'[ILC:IUC]'LINMAT';
        'INT' I:= ILF;
        'FOR' IC 'FROM' ILC 'TO' IUC
        'DO'  'REF'[,]'FLUXAC' FI  = GET(F[I]),
                               FIP = GET(F[I+1]);
              'INT' JLF = 1'LWB'FI,       JUF = 1'UPB'FI;
              'INT' JLC = (JLF+1)'OVER'2, JUC = JUF'OVER'2;
              'REF'[,]'FLUXAC' CIC = GEN(C[IC],JLC,JUC,-2,+2); 'NUL'CIC;
              'INT' J:= JLF;
              'FOR' JC 'FROM' JLC 'TO' JUC
              'DO'  'REF'[]'FLUXAC' FIJ   = FI [J ,], FIJP  = FI [J+1,],
                                    FIPJ  = FIP[J ,], FIPJP = FIP[J+1,],
                                    CICJC = CIC[JC,];
                    (CICJC[-2]+:=FIJ  [-2])+:=FIJP [-2];
                    (CICJC[-1]+:=FIJ  [-1])+:=FIPJ [-1];
                   ((CICJC[ 0]+:=FIJ  [ 0])+:=FIJ  [ 1])+:=FIJ  [ 2];
                   ((CICJC[ 0]+:=FIJP [ 0])+:=FIJP [-1])+:=FIJP [ 2];
                   ((CICJC[ 0]+:=FIPJ [ 0])+:=FIPJ [ 1])+:=FIPJ [-2];
                   ((CICJC[ 0]+:=FIPJP[ 0])+:=FIPJP[-1])+:=FIPJP[-2];
                    (CICJC[ 1]+:=FIJP [ 1])+:=FIPJP[ 1];
                    (CICJC[ 2]+:=FIPJ [ 2])+:=FIPJP[ 2];
                   J+:= 2
              'OD'  ;
              PET('FALSE',F[I ]); PET('FALSE',F[I+1]);
              PET('TRUE' ,C[IC]);
              I+:= 2
        'OD'; C
   'END'; #$E#
          #E$#
 
 # RELAX                          #                   'PR' EJECT 'PR'
 
   'PROC' RELAX    := ('INT' K,L, 'REF''FIELD'  Q, RHS)'VOID': ERROR;
 
   'BOOL' BACKWARD:= 'FALSE', REVERSE := 'FALSE', SYMMETRIC:= 'FALSE',
          RED     := 'FALSE', REDBLACK:= 'FALSE';
   'REAL' DAMPING :=       1, NLGSTOL :=    0.99,
          TIMESTEP:=  1.0E10, ODIAGSTEP:=   0.00;
   'REF'[]'REF'[,]'REAL' VOLUMES; # VOOR ELK LEVEL EEN 'REF'[,]'REAL'#
 
   #$H#
   'PROC' STEPPING  = ('INT' K,L, 'REF''FIELD'  Q, RHS)'VOID':
   #$B STEPP   B$#
   'TO' K 'WHILE'  WARNING
   'DO' 'FIELD' R := RHS'MINRESIDU'Q;
        'REF'[,]'STATE' QQ = Q'OF'Q, FF = Q'OF'R;
        ( VOLUMES[L] 'IS' 'NIL' ! VOLUMES[L]:= VOLUME(Q) );
        'REF'[,]'REAL' VOL = VOLUMES[L];
 
        (MONI   ! MONITOR("STEPPNG   ",R) );
        'INT'  IL=1'LWB'QQ,IU=1'UPB'QQ,JL=2'LWB'QQ,JU=2'UPB'QQ;
        'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'REAL' VOLI=VOL[I,];
                             'REF'[]'STATE'QQI=QQ[I,],FFI=FF [I,];
        'FOR' J 'FROM' JL 'TO' JU 'DO'
              QQI[J]+:= (TIMESTEP/VOLI[J]) * FFI[J]
        'OD' 'OD'
   'OD' ; #$E#
          #E$#
 
  'PROC' VOLM = ('POINT' A,B,C,D )'REAL' :
    'ABS' (0.5 * ( (X'OF'C-X'OF'A) * (Y'OF'D-Y'OF'B) -
                   (Y'OF'C-Y'OF'A) * (X'OF'D-X'OF'B) ) ) ;
 
   #$H#
   'PROC' VOLUME = ('FIELD' A) 'REF'[,]'REAL':
   #$B VOLUME  B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'A;
           'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           'REAL' DDX = DDX'OF'A, DDY = DDY'OF'A;
           'HEAP' [IL:IU,JL:JU]'REAL' VOL;
           [JL-1:JU]'POINT' COORD ; 'POINT' P1,P2 ;
           'REAL' XX:= (IL-1)*DDX,YY:=(JL-2)*DDY ;
           'FOR' J 'FROM' JL-1 'TO' JU
             'DO' COORD[J]:= MAPPING(XX,YY+:=DDY) 'OD' ;
           'FOR' I 'FROM' IL 'TO' IU
           'DO' YY:=(JL-1)*DDY ; 'REF'[]'REAL' VOLI=VOL[I,] ;
                P1:=MAPPING(XX+:=DDX,YY) ;
                'FOR' J 'FROM' JL 'TO' JU
                'DO' P2 := MAPPING(XX,YY+:=DDY) ;
                     VOLI[J]:= VOLM ( COORD[J-1],P1,P2,COORD[J] ) ;
                     COORD[J-1]:=P1 ; P1:=P2
                'OD';COORD[JU]:=P2
           'OD'; VOL
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' SOL4 = ([,]'REF''FLUXAC' A, []'REF''FLUX' R)'VOID':
 # NB: OVERSCHRIJFT A EN R !!!!!!!!!!!!!!! #
   #$B SOL4    B$#
 'BEGIN' 'FOR' I 'TO' 4
         'DO' A[I,I]:= 1/A[I,I];
              'FOR' J 'FROM' I+1 'TO' 4
              'DO'      A[I,J] := A[I,I]*A[I,J];
                   'FOR' K 'FROM' I+1 'TO' 4
                   'DO' A[K,J]-:= A[K,I]*A[I,J]
              'OD' 'OD';
                        R[I  ] := A[I,I]*R[I  ];
                   'FOR' K 'FROM' I+1 'TO' 4
                   'DO' R[K  ]-:= A[K,I]*R[I  ]
         'OD'      'OD';
 
         'FOR' I 'FROM' 3 'BY' -1 'TO' 1
         'DO''FOR' J 'FROM' I+1 'TO' 4
             'DO' R[I  ]-:= A[I,J]*R[J  ]
         'OD''OD'
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' RB4 RELAX = ('INT' K,L, 'REF''FIELD' QF,  RF)'VOID':
   #$B RB4REL  B$#
 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF;
         'INT' L1= 1'LWB'Q, U1= 1'UPB'Q,
               L2= 2'LWB'Q, U2= 2'UPB'Q;
         'BOOL'     BLACK:= 'FALSE';
         'REAL'     DAMP  = DAMPING;
         'REF''NETMAT' JAC= JACOBIAN[L];  'FIELD' RHS;
         'INT' IP,JP;
         [-2:2]'FLUXAC' JIJ,JIJP,JIPJ,JIPJP;
         'FLUXAC'       NUL1,NUL2,NUL3,NUL4;
         'FLUX'         RIJ,RIJP,RIPJ,RIPJP;
         [,]'REF''FLUXAC'JS=((JI J [ 0],JI J [ 1],JI J [ 2],  NUL1   ),
                             (JI JP[-1],JI JP[ 0],  NUL2   ,JI JP[ 2]),
                             (JIPJ [-2],  NUL3   ,JIPJ [ 0],JIPJ [ 1]),
                             (  NUL4   ,JIPJP[-2],JIPJP[-1],JIPJP[ 0]));
         []'REF''FLUX'   RS=(RIJ,RIJP,RIPJ,RIPJP);
 
         ('ODD'(1+U1-L1) 'OR' 'ODD'(1+U2-L2) ! ERROR );
 
  'FOR' KK 'TO' ( REDBLACK ! 2*K ! K )
  'DO'   # HET VOLGENDE KAN VEEL EFFICIENTER !! #
         'FIELD' QFOLD = ( MONILIN ! 'COPY'QF! 'SKIP');
 
         ( NEWLINSYS[L]
         ! LINSYS('FALSE',QF,RHS:= 'COPY'RF, JAC)
         ! RHS:= RF'MINRESIDU'QF
         );
         ( MONI ! MONITOR("RB4-RELAX ",RHS)   );
         'REF'[,]'STATE' R=Q'OF'RHS;
 
         ( REDBLACK ! BLACK:= RED:= 'NOT'RED );
 
         'FOR' I 'FROM' L1 'BY' 2 'TO' U1
         'DO' 'REF'[,]'FLUXAC' JACI = GET(JAC[I  ]),
                               JACIP= GET(JAC[I+1]);
         'FOR' J 'FROM' ( REDBLACK 'AND' BLACK ! L2+2 ! L2 )
                 'BY'   ( REDBLACK ! 4 ! 2 ) 'TO' U2
         'DO' IP:=I+1; JP:= J+1;
              JI J := JACI [J ,];  RI J := R[I ,J ];
              JI JP:= JACI [JP,];  RI JP:= R[I ,JP];
              JIPJ := JACIP[J ,];  RIPJ := R[IP,J ];
              JIPJP:= JACIP[JP,];  RIPJP:= R[IP,JP];
              NUL1:= NUL2:= NUL3:= NUL4:= NULFLUXAC;
 
              SOL4 ( JS, RS);
 
              'REAL' FACTOR:= DAMP;
              'WHILE' FACTOR*C3'OF'RI J /C3'OF'Q[I ,J ]<-0.99 'OR'
                      FACTOR*C3'OF'RIPJ /C3'OF'Q[IP,J ]<-0.99 'OR'
                      FACTOR*C3'OF'RI JP/C3'OF'Q[I ,JP]<-0.99 'OR'
                      FACTOR*C3'OF'RIPJP/C3'OF'Q[IP,JP]<-0.99
              'DO' FACTOR/:= 2 'OD';
              (FACTOR/=DAMP !PRINT((NEWLINE,"RB-WARNING",I,J,FACTOR)));
 
              Q[I ,J ]+:= FACTOR*RIJ ; Q[I ,JP]+:= FACTOR*RIJP;
              Q[IP,J ]+:= FACTOR*RIPJ; Q[IP,JP]+:= FACTOR*RIPJP
         'OD';BLACK:= 'NOT'BLACK;
              PET('FALSE',JAC[I]); PET('FALSE',JAC[I+1])
         'OD';
         ( MONILIN
         ! PRINT((NEWLINE,"LINEAR CHECK RB-RELAX",NEWLINE));
           RHS-:= JAC*(QF-QFOLD);
           'PRINT' RHS;
           MONITOR("RB-LINEAR ",RHS)
         )
  'OD'
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' RB2 RELAX = ('INT' K,L, 'REF''FIELD' QF,  RF)'VOID':
   #$B RB2REL  B$#
 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF;
         'INT' L1= 1'LWB'Q, U1= 1'UPB'Q,
               L2= 2'LWB'Q, U2= 2'UPB'Q;
         'BOOL' BLACK:= 'FALSE';
         'REAL' DAMP = DAMPING;
         'REF''NETMAT' JAC= JACOBIAN[L];  'FIELD' RHS;
         'STATE' X;
 
  'FOR' KK 'TO' ( REDBLACK ! 2*K ! K )
  'DO'   # KAN EFFICIENTER #
         'FIELD' QFOLD = ( MONILIN ! 'COPY'QF! 'SKIP');
 
         ( NEWLINSYS[L]
         ! LINSYS('FALSE',QF,RHS:= 'COPY'RF, JAC)
         ! RHS:= RF'MINRESIDU'QF
         );
         ( MONI ! MONITOR("RB2-RELAX ",RHS)   );
         'REF'[,]'STATE' R=Q'OF'RHS;
 
         ( REDBLACK ! BLACK:= RED:= 'NOT' RED );
 
         'FOR' I 'FROM' L1 'TO' U1
         'DO' 'REF'[,]'FLUXAC' JACI = GET(JAC[I  ]);
         'FOR' J 'FROM' ( REDBLACK 'AND' BLACK ! L2+1 ! L2 )
                 'BY'   ( REDBLACK ! 2 ! 1 ) 'TO' U2
         'DO' X:= R[I,J]/JACI[J,0];
 
              'REAL' FACTOR:= DAMP;
              'WHILE' FACTOR*C3'OF'X/C3'OF'Q[I,J]<-0.99
              'DO' FACTOR/:= 2 'OD';
              (FACTOR/=DAMP !PRINT((NEWLINE,"RB-WARNING",I,J,FACTOR)));
 
              Q[I,J]+:= ( FACTOR=1! X ! FACTOR*X )
         'OD';BLACK:= 'NOT'BLACK;
              PET('FALSE',JAC[I])
         'OD';
         ( MONILIN
         ! PRINT((NEWLINE,"LINEAR CHECK RB-RELAX",NEWLINE));
           RHS-:= JAC*(QF-QFOLD);
           'PRINT' RHS;
           MONITOR("RB-LINEAR ",RHS)
         )
  'OD'
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' JACOBI    = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID':
   #$B JACOBI  B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS;
                               STATE Z ( QF );
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF;
 
      [JL-1:JU]   'WALL'  HORIN,HORIZ,VERI ;
      [JL-1:JU]   'POINT' COORD            ;
      [JL  :JU]   'STATE' QIM,QIO          ;
      'STATE'  UC, QOUT;  'STATAC' QAC     ;
      'FLUX'   NF,F,DELTA                  ;
      'FLUXAC' NFAC0,NFAC1,FAC             ;
      'REAL'   ERROR,FACTOR,DAMP:=DAMPING  ;
 
      'TO' K
      'DO' WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIN);
         'FOR' I 'FROM' IL 'TO' IU
         'DO''REAL'XX=(I-0.5)*DDX;
             WALLS( 'TRUE','TRUE', I  *DDX,DDY,COORD, VERI,HORIZ);
             QIO:= Q[I,];
 
         'FOR' J 'FROM' JL 'TO' JU
         'DO''REAL'YY=(J-0.5)*DDY;
             'REF''STATE' QIJ=Q[I,J]; UC:=QIJ;
 
         'TO' 13 'WHILE'
 
               F   :=(I=IU ! QOUT :=  ZUID(HORIZ[J],XXZ,YY,UC,QAC);
                             FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL'));
               FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC);
                             FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(HORIN[J],QIM[J],UC,'NIL',NFAC1));
               FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 );
 
               F  +:=(J=JU ! QOUT :=  OOST(VERI[J],XX,YYO,UC,QAC);
                             FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL'));
               FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(J=JL ! QOUT :=  WEST(VERI[J-1],XX,YYW,UC,QAC);
                             FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(VERI[J-1],QIO[J-1],UC,'NIL',NFAC1));
               FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 );
 
               DELTA := FF[I,J]-F;
 
               ERROR := 'ABS'DELTA;
               DELTA := DELTA/FAC;
 
               FACTOR:= 1.0;
                  'WHILE' FACTOR * 'ABS'(C3'OF'DELTA/C3'OF'UC) > 0.99
                  'DO' FACTOR/:=2 'OD';
               ( FACTOR/=1.0 ! DELTA*:= FACTOR );
 
               UC+:= DELTA;
 
               (MONIGS! PRINT((NEWLINE,
                        "NEWT."+SF(ERROR)+WHOLE(I,5)+WHOLE(J,5)  ));
                        (FACTOR/=DAMP! PRINT(("  "+SF(FACTOR) )) ));
 
            ERROR > NLGSTOL
            'DO' 'SKIP' 'OD';
 
            QIJ+:= DAMP*(UC-QIJ)
            'OD';
 
            QIM:= QIO; HORIN:= HORIZ
            'OD';
 
            (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE));
                         CAVITATION:= 'FALSE')
      'OD';
           (MONI   !     MONITOR("JACOBI    ",RHS'MINRESIDU'QF))
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
   'PROC' SEIDEL    = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID':
   #$B SEIDEL  B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS;
                               STATE Z ( QF );
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF;
           'INT'  ISL,IS,ISU,JSL,JS,JSU;
           (BACKWARD! ISL:=IU;ISU:=IL;IS:=-1 ! ISL:=IL;ISU:=IU;IS:=1);
           (REVERSE ! JSL:=JU;JSU:=JL;JS:=-1 ! JSL:=JL;JSU:=JU;JS:=1);
 
      [JL-1:JU] 'WALL' HORIN,HORIZ,VERI;
      [JL-1:JU]'POINT' COORD;
      'STATE'  QOUT;  'STATAC' QAC;
      'FLUX'   NF,F,DELTA;
      'FLUXAC' NFAC0,NFAC1,FAC;
      'REAL'   ERROR,FACTOR,DAMP:=DAMPING;
 
      'FOR' KK 'TO' ( SYMMETRIC'OR'REDBLACK ! 2*K ! K)
      'DO' 'BOOL' SHIFT:= RED = 'ODD'KK;
             ( IS=1
             ! WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIZ)
             ! WALLS('FALSE','TRUE', IU   *DDX,DDY,COORD,'NIL',HORIN)
             );
         'FOR' I 'FROM' ISL 'BY' IS 'TO' ISU
         'DO''REAL'XX=(I-0.5)*DDX;
             ( IS=1
             ! HORIN:= HORIZ;
               WALLS( 'TRUE','TRUE',  I   *DDX,DDY,COORD, VERI,HORIZ)
             ! HORIZ:= HORIN;
               WALLS('FALSE','TRUE',( I-1)*DDX,DDY,COORD, VERI,HORIN)
             );
 
             'INT' JJSL:= JSL, JJS:= JS;
             (REDBLACK! (SHIFT:='NOT'SHIFT! JJSL+:=JS ); JJS+:=JS );
 
         'FOR' J 'FROM' JJSL 'BY' JJS 'TO' JSU
         'DO''REAL'YY=(J-0.5)*DDY;
 
         'TO' 13 'WHILE' 'REF''STATE' UC=Q[I,J];
 
               F   :=(I=IU ! QOUT :=  ZUID(HORIZ[J],XXZ,YY,UC,QAC);
                             FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL'));
               FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC);
                             FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(HORIN[J],Q[I-1,J],UC,'NIL',NFAC1));
               FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 );
 
               F  +:=(J=JU ! QOUT :=  OOST(VERI[J],XX,YYO,UC,QAC);
                             FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL'));
               FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(J=JL ! QOUT :=  WEST(VERI[J-1],XX,YYW,UC,QAC);
                             FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(VERI[J-1],Q[I,J-1],UC,'NIL',NFAC1));
               FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 );
 
               DELTA := FF[I,J]-F;
 
               ERROR := 'ABS'DELTA;
               DELTA := DELTA/FAC;
 
               FACTOR:= DAMP;
                  'WHILE' FACTOR * 'ABS'(C3'OF'DELTA/C3'OF'UC) > 0.99
                  'DO' FACTOR/:=2 'OD';
               ( FACTOR/=1 ! DELTA*:= FACTOR );
 
               UC+:= DELTA;
 
               'IF' MONIGS 'THEN'
                 PRINT(( NEWLINE,"NEWT." + SF( ERROR ) +
                         WHOLE(I,5) + WHOLE(J,5)      ));
                 (FACTOR/=DAMP!PRINT(("  "+SF( FACTOR ) )) )
               'FI';
            ERROR > NLGSTOL
            'DO' 'SKIP'
            'OD' 'OD' 'OD';
 
            (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE));
                         CAVITATION:= 'FALSE') ;
 
            'IF' SYMMETRIC
            'THEN' 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS;
                         IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS
            'FI'
      'OD';
           (MONI   !     MONITOR("SEIDEL    ",RHS'MINRESIDU'QF))
   'END'; #$E#
          #E$#
 
 # STAR-RELAXATIE #
                                'BOOL' BLOCK ;
   #$H#
   'PROC' STARRELAX = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID':
   #$B STRREL  B$#
   'BEGIN' 'REF'[,]'STATE' QQ = Q'OF'QF, FF = Q'OF'RHS;
                           STATE Z ( QF );
      'REAL' DDXH = 2*DDX'OF'QF, DDYH = 2*DDY'OF'QF;
      'INT'  ILH=(1'LWB'FF+1)'OVER'2,IUH=(1'UPB'FF)'OVER'2,
             JLH=(2'LWB'FF+1)'OVER'2,JUH=(2'UPB'FF)'OVER'2;
      'INT'  ISL,IS,ISU,JSL,JS,JSU;
      (BACKWARD!ISL:=IUH;ISU:=ILH;IS:=-1 !ISL:=ILH;ISU:=IUH;IS:=1);
      (REVERSE !JSL:=JUH;JSU:=JLH;JS:=-1 !JSL:=JLH;JSU:=JUH;JS:=1);
 
      [JLH-1:JUH] 'WALL' HORIN,HORIZ,VERT;
      [JLH-1:JUH]'POINT' COORDN,COORDZ;
 
      'REAL'  ERROR,FACTOR, DAMP:=DAMPING*0.5;
      [1:4]'STATE' QS,DELTA;
      [1:4]'FLUX' RS;
      [1:4,1:4]'FLUXAC' JAC ;
      'REF''STATE'  QA= QS[  1],QB= QS[  2],QC= QS[  3],QD= QS[  4];
      'REF''FLUX'   RA= RS[  1],RB= RS[  2],RC= RS[  3],RD= RS[  4];
      'REF''FLUXAC' AA=JAC[1,1],AB=JAC[1,2],AC=JAC[1,3],AD=JAC[1,4],
                    BA=JAC[2,1],BB=JAC[2,2],BC=JAC[2,3],BD=JAC[2,4],
                    CA=JAC[3,1],CB=JAC[3,2],CC=JAC[3,3],CD=JAC[3,4],
                    DA=JAC[4,1],DB=JAC[4,2],DC=JAC[4,3],DD=JAC[4,4];
 
  #
                           NOORD
                           (I-1,J)
                   O---------------------O
                   ! \     (I,J+1)     / !
                   !   \             /   !                        B
                   !     \         /     !                        !
                   !       \     / (I+1, !                        !
        (I+1,J-1)  !  (I,J)   O     J+1) !  (I,J+2)          A ---.--- C
        WEST       !       /     \       !     OOST               !
                   !     /         \     !                        !
                   !   /             \   !                        D
                   O /     (I+1,J)     \ !
                   O---------------------O
                         (I+2,J+1)
                          ZUID
  #
  #!# 'FOR' KK 'TO' ( SYMMETRIC'OR'REDBLACK ! 2*K ! K)
      'DO' 'BOOL' SHIFT:= RED = 'ODD'KK;
            ( IS=1
            !WALLS( 'TRUE','TRUE',(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIZ)
            !WALLS('FALSE','TRUE', IUH   *DDXH,DDYH,COORDN,'NIL',HORIN)
            );
  #!#   'FOR' IH 'FROM' ISL 'BY' IS 'TO' ISU
        'DO''REAL'XX=(IH-0.5)*DDXH;
            ( IS=1
            ! HORIN:= HORIZ; COORDN:= COORDZ;
             WALLS( 'TRUE','TRUE',  IH   *DDXH,DDYH,COORDZ, VERT,HORIZ)
            ! HORIZ:= HORIN; COORDZ:= COORDN;
             WALLS('FALSE','TRUE',( IH-1)*DDXH,DDYH,COORDN, VERT,HORIN)
            );
 
            'INT' JJSL:= JSL, JJS:= JS;
            (REDBLACK! (SHIFT:='NOT'SHIFT! JJSL+:=JS ); JJS+:=JS );
 
  #!#   'FOR' JH 'FROM' JJSL 'BY' JJS 'TO' JSU
        'DO''REAL'YY=(JH-0.5)*DDYH;
            'BOOL' OK:= BLOCK;
 
        'INT' I = 2*IH-1, J = 2*JH-1, JMH = JH-1;
        'POINT' MID = MAPPING(XX,YY);
        QS:= (QQ[I  ,J  ], QQ[I  ,J+1],
              QQ[I+1,J+1], QQ[I+1,J  ]);
 
        'PROC' DIAGH = ('POINT' CP)'WALL':
        'BEGIN' 'REAL' DX := X'OF'CP - X'OF'MID,
                       DY := Y'OF'CP - Y'OF'MID, DS;
                       DS := SQRT(DX*DX+DY*DY);
               'WALL'(DY/DS,-DX/DS,DS)
        'END';
 
  #$!#   'TO' 8 'WHILE'
 
 'PROC' RAJAC = 'VOID':
 'BEGIN' # BEREKENT LOKAAL:  RS:= "RHS-N(Q)"  EN  JAC:= "DN/DQ" #
        'WALL'  CW;
        'STATE' QOUT;
        'STATAC'QAC;
        'FLUX'  NF;
        'FLUXAC' NF0, NF1;
 
        RA:= FF[I  ,J  ]; RB:= FF[I  ,J+1];
        RC:= FF[I+1,J+1]; RD:= FF[I+1,J  ];
        'FOR' K 'TO' 4
        'DO' 'FOR' L 'TO' 4
             'DO' JAC[K,L]:= NULFLUXAC
        'OD' 'OD';
 
        ( IH=ILH
        ! QOUT :=NOORD(HORIN[JH],XXN,YY,QB,QAC);
          RB  +:=FLUXT(HORIN[JH],  QOUT     ,QB,NF0,NF1);
                 BB-:= (NF1+:=NF0*QAC)
        ! RB  +:=FLUXT(HORIN[JH],QQ[I-1,J],QB,'NIL',NF1);
                 BB-:=  NF1
        );
 
        ( IH=IUH
        ! QOUT :=ZUID (HORIZ[JH],XXZ,YY,QD,QAC);
          RD  -:=FLUXT(HORIZ[JH],QD,QOUT,NF0,NF1);
                 DD+:= (NF0+:=NF1*QAC)
        ! RD  -:=FLUXT(HORIZ[JH],QD,QQ[I+2,J+1],NF0,'NIL');
                 DD+:=  NF0
        );
 
        (JH=JUH
        ! QOUT :=OOST (VERT[JH],XX,YYO,QC,QAC);
          RC  -:=FLUXT(VERT[JH],QC,QOUT,NF0,NF1);
                 CC+:= (NF0+:= NF1*QAC)
        ! RC  -:=FLUXT(VERT[JH],QC,QQ[I,J+2],NF0,'NIL');
                 CC+:=  NF0
        );
 
        (JH=JLH
        ! QOUT :=WEST (VERT[JMH],XX,YYW,QA,QAC);
          RA  +:=FLUXT(VERT[JMH],QOUT,QA,NF0,NF1);
                 AA-:= (NF1+:= NF0*QAC)
        ! RA  +:=FLUXT(VERT[JMH],QQ[I+1,J-1],QA,'NIL',NF1);
                 AA-:=  NF1
        );
 
 
        CW:= DIAGH(COORDN[JMH]);
        NF:= FLUXTG(CW,QA,QB,NF0,NF1);
             RA-:=NF; RB+:=NF; AA+:=NF0; AB+:=NF1; BA-:=NF0; BB-:= NF1;
        CW:= DIAGH(COORDN[JH] );
        NF:= FLUXTG(CW,QB,QC,NF0,NF1);
             RB-:=NF; RC+:=NF; BB+:=NF0; BC+:=NF1; CB-:=NF0; CC-:= NF1;
        CW:= DIAGH(COORDZ[JH] );
        NF:= FLUXTG(CW,QC,QD,NF0,NF1);
             RC-:=NF; RD+:=NF; CC+:=NF0; CD+:=NF1; DC-:=NF0; DD-:= NF1;
        CW:= DIAGH(COORDZ[JMH]);
        NF:= FLUXTG(CW,QD,QA,NF0,NF1);
             RD-:=NF; RA+:=NF; DD+:=NF0; DA+:=NF1; AD-:=NF0; AA-:= NF1
 'END' # RAJAC# ;
 
        RAJAC;
        ERROR := 'ABS'RA + 'ABS'RB + 'ABS'RC + 'ABS'RD;
        ( MONIGS ! PRINT((NEWLINE,"GS ERROR VOOR "+WHOLE(JUH-JLH+1,5)+
                          WHOLE(IH,5)+WHOLE(JH,5), ERROR)) );
 
        'IF' ( OK   # NEWTON OR DAMPED JACOBI ? #
             ! DELTA  := RS/JAC;
               'FOR' C 'TO' 4 'WHILE' OK
               'DO' OK:= 'ABS'(C3'OF'DELTA[C]/C3'OF'QS[C]) < 0.99 'OD'
             ); OK
        'THEN' 'FOR' C 'TO' 4         # NEWTON GS #
               'DO' QS[C]+:= DELTA[C] 'OD'
        'ELSE' 'FOR' C 'TO' 4         # DAMPED JACOBI #
               'DO' FACTOR:= DAMP;
                    'WHILE' 'ABS'(C3'OF'DELTA[C]/C3'OF'QS[C])
                             > 0.9/FACTOR
                    'DO' FACTOR/:=2 'OD';
                    ( FACTOR<=0.1 'OR' ( BLOCK'AND'MONIGS)
                    ! PRINT((NEWLINE,"DAMPED JACOBI FACTOR",FACTOR)) );
                    ( FACTOR >0.1 ! QS[C]+:= ( DELTA[C]*:= FACTOR ) )
               'OD'
        'FI';
 #
 'IF'   ( 'FALSE'; MONIGS 'AND' (BLOCK'AND''NOT'OK) ); 'FALSE'
 'THEN' [1:4]'FLUX'  ROLD:= RS, RNEW;
        [1:4]'STATE' QOLD:= (QQ[I  ,J  ], QQ[I  ,J+1],
                             QQ[I+1,J+1], QQ[I+1,J  ]);
        [1:4,1:4]'FLUXAC' JACOLD:= JAC;
 
  PRINT((NEWLINE,NEWLINE,"IH,JH:"+WHOLE(IH,0)+"  "+WHOLE(JH,0),
         NEWLINE,"OPLOSSING VOOR CORRECTIE ",
         NEWLINE,QOLD,
         NEWLINE,"RESIDU VOOR CORRECTIE ",'ABS'ROLD,
         NEWLINE,ROLD));
  ('REAL' S ='ABS'(ROLD-JAC*DELTA); S > 1.0E-10
  ! PRINT((NEWLINE,"DAMPING OR BAD LIN. SOLVING",S))      );
 
  RAJAC;  RNEW:= RS;
  PRINT((NEWLINE,"RESIDU NA CORRECTIE ",'ABS'RNEW,
         NEWLINE,"CONTROLE NA DE SLAG A ",'ABS'(RNEW-ROLD),
         NEWLINE,RNEW-ROLD,
         NEWLINE,"CONTROLE NA DE SLAG B ",'ABS'(JACOLD*DELTA),
         NEWLINE,JACOLD*DELTA,
         NEWLINE,(RNEW-ROLD)/(JACOLD*DELTA),
         NEWLINE,"RELAT VERBETERING",
            'ABS'(RNEW-ROLD+JACOLD*DELTA)/'ABS'(RNEW-ROLD)     ))
 'FI';
 #
        ERROR > NLGSTOL    'AND'  OK
        'DO' 'SKIP' 'OD';
        QQ[I  ,J  ]:= QA; QQ[I  ,J+1]:= QB;
        QQ[I+1,J+1]:= QC; QQ[I+1,J  ]:= QD
  'OD'  'OD';
 
            (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE));
                         CAVITATION:= 'FALSE') ;
 
            'IF' SYMMETRIC
            'THEN' 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS;
                         IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS
            'FI'
      'OD';
           (MONI   !     MONITOR("STARRELAX ",RHS'MINRESIDU'QF))
   'END'; #$E#
          #E$#
 
 # NEWTON ITERATION               #                   'PR' EJECT 'PR'
 
   #$H#
 'PROC' NEWTON = ('INT' K,L, 'REF''FIELD' QF,RF) 'VOID':
   #$B NEWTON  B$#
 'BEGIN''REF''NETMAT' JAC= JACOBIAN[L], DEC= DECOMP[L],
        'FIELD' RHS,
        'INT' IL= 1'LWB'(Q'OF'QF), IU= 1'UPB'(Q'OF'QF);
        'REAL' ROLD:= MAXREAL, RNEW,LTOL;
 
     'FOR' KK 'TO' K
     'WHILE'( KK=1 'OR' NEW LIN SYS[L]
            ! (MONINWT !PRINT((NEWLINE,"UPDATE JACOBIANS ",NEWLINSYS)));
              LINSYS ('FALSE',QF,RHS:='COPY'RF,JAC);
              RNEW:= FLUXNORM(RHS,'LOC''FLUX','LOC''FLUX');
              NEWLINSYS[L]:= 'FALSE';
              ( ILLU ! ILLUDEC(JAC,DEC) );
 
              'FOR' LEV 'FROM' L-1 'BY'-1 'TO' 1
              'WHILE' NEWLINSYS[LEV]
              'DO'  JACOBIAN [LEV]:= 'RAP'JACOBIAN[LEV+1];
                    NEWLINSYS[LEV]:= 'FALSE';
                    ( ILLU ! ILLUDEC(JACOBIAN[LEV],DECOMP[LEV]) )
              'OD'
            );
            LTOL:= ('REAL' Z = (RNEW<1.0 ! RNEW*RNEW ! RNEW);
                               (Z>LINRESTOL!Z!LINRESTOL)     );
            ( MONINWT ! PRINT((NEWLINE,"NEWTON: "+FLOAT(RNEW,10,3,3)+
                         "DESIRED RESIDUAL = "+FLOAT(LTOL,10,3,3) )) );
     RESTOL<RNEW 'AND' RNEW<=ROLD 'AND' ACCOUNT OK
     'DO' 'FIELD' W:= 'NUL''COPY'QF, RES;      ROLD:= RNEW;
          'REAL' OLDNORM:= RNEW, ONORM:= RNEW, NORM:= RNEW;
           (MONINWT ! PRINT((NEWLINE,"LINSOLVE"+21*" "+
            FLOAT(OLDNORM,10,3,3)+" "+
            FIXED(-4*LN(OLDNORM)/LN(10.0),6,1) )));
          'INT' CM;
          'FOR' C 'TO' MLINSOLVE
          'WHILE' LTOL<NORM 'AND' NORM<=ONORM
          'DO' CM:=C; ONORM:= NORM;
               RES:= 'COPY' RHS;
               ( C>1 ! PLUS DIAG (QF,W,RES) );
               LINSOLVE(1,L,W,RES);
               NORM:= FLUXNORM(RES-:=JAC*W,'LOC''FLUX','LOC''FLUX');
               (MONINWT! PRINT((NEWLINE,"LINSOLVE "+20*" "+
                  FLOAT(NORM,10,3,3)+" "+
                  FIXED(-4*LN(NORM)/LN(10.0),6,1) )));
               ( MONI 'OR' NORM>=OLDNORM
               ! MONITOR("LINSOLV-NL",RF'MINRESIDU'(QF+W))
               )
          'OD';
 
          (MONINWT ! PRINT((NEWLINE,"LINSOLVE CONV = ",
                     FLOAT( EXP((LN(NORM/OLDNORM))/CM),10,3,3) )) );
 
          ( NORM > OLDNORM      # DIVERGENCE #
          ! ODIAGSTEP*:= 10; LINEARIZE;
            PRINT((NEWLINE, "NEWTON STEP NOT ACCEPTED! "))
          ! QF+:= W
          );
 
          ( ONORM < 1.5*NORM    # SLOW CONVERGENCE #
          ! ODIAGSTEP*:=  2; LINEARIZE;
            PRINT((NEWLINE,"STEP CHANGED:"+
                   SF( ODIAGSTEP )+50*"*",NEWLINE))
          );
 
          RHS := RF'MINRESIDU'QF;
          RNEW:= FLUXNORM(RHS,'LOC''FLUX','LOC''FLUX');
          ( # VGL. LIN-NONLIN RESIDU: # RNEW/NORM > 1.2 ! LINEARIZE )
     'OD'
   'END'; #$E#
          #E$#
 
 # APPROXIMATE LINEAR SOLVERS     #                   'PR' EJECT 'PR'
 
  'REF'[]'NETMAT' JACOBIAN,DECOMP;
  'REF'[]'BOOL'   NEW LIN SYS;
  'BOOL' ILLU:= 'TRUE';
  'INT'  MLINSOLVE:= 13;
  'REAL' LINRESTOL:= 1.0E-11, RESTOL:= 1.0E-11;
  'PROC' LINSOLVE := ('INT' K,L, 'REF''FIELD' Q,  R)'VOID':'SKIP';
 
   #$H#
  'PROC' LINEARIZE = 'VOID':
   #$B LNRIZE  B$#
    'FOR' I 'TO' 'UPB'NEWLINSYS
    'DO' NEWLINSYS[I]:='TRUE'
    'OD'; #$E#
          #E$#
 
   #$H#
  'PROC' INIZ LIN = ('INT' LEVELS)'VOID':
   #$B INIZLIN B$#
  'BEGIN' JACOBIAN := 'HEAP'[1:LEVELS]'NETMAT';
          DECOMP   := 'HEAP'[1:LEVELS]'NETMAT';
          NEWLINSYS:= 'HEAP'[1:LEVELS]'BOOL';
          LINEARIZE
  'END'; #$E#
         #E$#
 
 # LINEAR RELAXATION              #
 
 'BOOL' LINSYM:= 'TRUE', LINREVER:= 'FALSE', LINBACKW:= 'FALSE';
 'REAL' DAMP:= 1.0;
 
   #$H#
 'PROC' LU RELAX = ('INT' K, L, 'REF''FIELD' QF, RHS)  'VOID':
   #$B LUREL   B$#
 'BEGIN'# BY PDZ; INSERTED 850102  #
        'TO' K
        'DO' 'FIELD' RES:= RHS-JACOBIAN[L]*QF;
             ILLURELAX(DECOMP[L],JACOBIAN[L],QF,RES)
        'OD'
   'END'; #$E#
          #E$#
 
                                                      'PR' EJECT 'PR'
   #$H#
 'PROC' LINJACO = ('INT' K, L, 'REF''FIELD' F, RHS) 'VOID':
   #$B LINJC   B$#
 'BEGIN'# BY PDZ; INSERTED 850102  #
   'REF''NETMAT'  JAC = JACOBIAN[L],
   'REF'[,]'STATE' QR = Q'OF'RHS;
   'INT' L1 = 1'LWB'QR, U1 = 1'UPB'QR,
         L2 = 2'LWB'QR, U2 = 2'UPB'QR;
   ( L1/= 'LWB'JAC 'OR' U1/= 'UPB'JAC ! ERROR );
   'REAL' D1 = 1/DAMP, D2 = (DAMP-1)/DAMP;
   'TO'   K
   'DO'   'FIELD' G = 'COPY'RHS; TYPE'OF'G:= TYPE'OF'F;
          'REF'[,]'STATE' QG = Q'OF'G, QF = Q'OF'F;
          ( L1/=1'LWB'QF 'OR' U1/=1'UPB'QF 'OR'
            L2/=2'LWB'QF 'OR' U2/=2'UPB'QF      ! ERROR );
          'FOR' I 'FROM' L1 'TO' U1
          'DO'  'BOOL' LWB = (I=L1), UPB = (I=U1);
                'REF'[,]'FLUXAC' JACI = GET(JAC[I]),
                'REF'[]'STATE'   QRI  = QR[        I   ,],
                                 QGI  = QG[        I   ,],
                                 QFM  = QF[(LWB!L1!I-1),],
                                 QFI  = QF[        I   ,],
                                 QFP  = QF[(UPB!U1!I+1),];
                'INT' L3 = ('INT'L=2'LWB'JACI; (LWB!(L<-1!-1!L)!L)),
                      U3 = ('INT'U=2'UPB'JACI; (UPB!(U> 1! 1!U)!U));
                'FOR' J 'FROM' L2 'TO' U2
                'DO'  'BOOL' NOTL = (J>L2), NOTR = (J<U2),
                      'REF'[]'FLUXAC' JACIJ = JACI[J,],
                      'REF''STATE'    QGIJ  =  QGI[J];
                      'FOR' N 'FROM' L3 'TO' U3
                      'DO'  'CASE' N+3
                            'IN'           QGIJ-:= JACIJ[N]*QFM[J  ] ,
                                   (NOTL ! QGIJ-:= JACIJ[N]*QFI[J-1]),
                                                                'SKIP',
                                   (NOTR ! QGIJ-:= JACIJ[N]*QFI[J+1]),
                                           QGIJ-:= JACIJ[N]*QFP[J  ]
                            'OUT'  ERROR
                            'ESAC'
                      'OD'  ;
                      QGIJ*:= D1;
                      QGIJ := QGIJ/JACIJ[0]+D2*QFI[J]
                'OD'  ;
                PET('FALSE',JAC[I])
          'OD'  ;
          F:= G
   'OD'
   'END'; #$E#
          #E$#
                                                      'PR' EJECT 'PR'
   #$H#
 'PROC' LINGAUS = ('INT' K, L, 'REF''FIELD' F, RHS) 'VOID':
   #$B LINGS   B$#
 'BEGIN'# BY PDZ; INSERTED 850102  #
   'REF''NETMAT'  JAC = JACOBIAN[L],
   'REF'[,]'STATE' QF = Q'OF'F, QR = Q'OF'RHS;
   'INT' L1 = 1'LWB'QF, U1 = 1'UPB'QF,
         L2 = 2'LWB'QF, U2 = 2'UPB'QF;
   ( L1/= 'LWB'JAC 'OR' U1/= 'UPB'JAC 'OR'
     L1/=1'LWB'QR  'OR' U1/=1'UPB'QR  'OR'
     L2/=2'LWB'QR  'OR' U2/=2'UPB'QR      ! ERROR);
   'INT' ISL, IS, ISU, JSL, JS, JSU;
   (LINBACKW ! ISL:=U1; ISU:=L1; IS:=-1 ! ISL:=L1; ISU:=U1; IS:= 1);
   (LINREVER ! JSL:=U2; JSU:=L2; JS:=-1 ! JSL:=L2; JSU:=U2; JS:= 1);
   'TO'   ( LINSYM ! 2*K ! K )
   'DO'   'FOR' I 'FROM' ISL 'BY' IS 'TO' ISU
          'DO'  'BOOL' LWB = (I=L1), UPB = (I=U1);
                'REF'[,]'FLUXAC' JACI = GET(JAC[I]),
                'REF'[]'STATE'   QRI  = QR[        I   ,],
                                 QFM  = QF[(LWB!L1!I-1),],
                                 QFI  = QF[        I   ,],
                                 QFP  = QF[(UPB!U1!I+1),];
                'INT' L3 = ('INT'L=2'LWB'JACI; (LWB!(L<-1!-1!L)!L)),
                      U3 = ('INT'U=2'UPB'JACI; (UPB!(U> 1! 1!U)!U));
                'FOR' J 'FROM' JSL 'BY' JS 'TO' JSU
                'DO'  'BOOL' NOTL = (J>L2), NOTR = (J<U2),
                      'REF'[]'FLUXAC' JACIJ = JACI[J,],
                      'STATE'           RIJ:= QRI[J]  ;
                      'FOR' N 'FROM' L3 'TO' U3
                      'DO'  'CASE' N+3
                            'IN'           RIJ-:= JACIJ[N]*QFM[J  ] ,
                                   (NOTL ! RIJ-:= JACIJ[N]*QFI[J-1]),
                                                               'SKIP',
                                   (NOTR ! RIJ-:= JACIJ[N]*QFI[J+1]),
                                           RIJ-:= JACIJ[N]*QFP[J  ]
                            'OUT'  ERROR
                            'ESAC'
                      'OD'  ;
                      QFI[J]:= RIJ/JACIJ[0]
                'OD'  ;
                PET('FALSE',JAC[I])
          'OD';
          (LINSYM ! 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS;
                          IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS )
   'OD'
   'END'; #$E#
          #E$#
 
 # MDCP ROUTINES                  #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' TAU CORRECT = ('BOOL' SHOW, EXTRA, 'REF''FIELD' QL,RL)'VOID':
   #$B TAU.COR B$#
   'BEGIN' 'FIELD' TAU; 'BOOL' MORE = SHOW 'OR' EXTRA;
           Q'OF'RL:= 'NIL';  RL:= 'NULRHS' QL;
           RESIDU2('FALSE',QL,RL);
                                  ( SHOW
                                  ! PRINT((NEWLINE, "TAU CORRECT",
                                    NEWLINE, "NORM VAN N2(QH)"));
                                    'PRINTNORM' RL );
           (MORE! TAU:= 'RESTRICT' RL );
           RESIDU1( 'TRUE',QL,RL);
                                  ( SHOW ! PRINT((NEWLINE,
                                    "NORM VAN NH1(QH)-N2(QH)"));
                                    'PRINTNORM' RL );
           (MORE! RESIDU2('TRUE','RESTRICT'QL,TAU);
                                  ( SHOW ! PRINT((NEWLINE,
                                    "NORM VAN TAU[2H,H]"));
                                    'PRINTNORM' TAU);
            (EXTRA!     TAU*:= (1/3);
                        TAU := 'PROSCND' TAU;
                        RL +:= TAU
           ))
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' MDCP SMOOTH  = ('FIELD' QF )'VOID':
   #$B MDCP.SM B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF;  STATE Z ( QF );
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           [JL-1:JU]   'WALL'  HORIN,HORIZ,VERI ;
           [JL-1:JU]   'POINT' COORD            ;
           [JL  :JU]   'STATE' QIM,QIO          ;
           'STATE'  UC, QOUT;  'STATAC' QAC     ;
           'FLUX'   NF,F                        ;
           'FLUXAC' NFAC0,NFAC1,FAC             ;
 
         WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIN);
         'FOR' I 'FROM' IL 'TO' IU
         'DO''REAL'XX=(I-0.5)*DDX;
             WALLS( 'TRUE','TRUE', I  *DDX,DDY,COORD, VERI,HORIZ);
             QIO:= Q[I,];
 
         'FOR' J 'FROM' JL 'TO' JU
         'DO''REAL'YY=(J-0.5)*DDY;
             'REF''STATE' QIJ=Q[I,J]; UC:=QIJ;
 
               F   :=(I=IU ! QOUT :=  ZUID(HORIZ[J],XXZ,YY,UC,QAC);
                             FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL'));
               FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC);
                             FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(HORIN[J],QIM[J],UC,'NIL',NFAC1));
               FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 );
 
               F  +:=(J=JU ! QOUT :=  OOST(VERI[J],XX,YYO,UC,QAC);
                             FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1)
                           ! FLUXT(VERI[J],UC,QIO[J+1],NFAC0,'NIL'));
               FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 );
 
               F  -:=(J=JL ! QOUT :=  WEST(VERI[J-1],XX,YYW,UC,QAC);
                             FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1)
                           ! FLUXT(VERI[J-1],QIO[J-1],UC,'NIL',NFAC1));
               FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 );
 
               QIJ +:= -0.5 * F /FAC
 
            'OD'; QIM:= QIO; HORIN:= HORIZ
            'OD'
   'END'; #$E#
          #E$#
 
 # MULTIGRID (I)                  #                   'PR' EJECT 'PR'
 
   'PROC' LIN RELAX:= ('INT' K,L, 'REF''FIELD'  Q,   R)'VOID': 'SKIP';
   'INT' PCS:=1, SIGMACS:=1, QCS:=1;
 
   #$H#
   'PROC' MGCS      = ('INT' K,L, 'REF''FIELD'  Q, RHS)'VOID':
   #$B MGCS    B$#
   'BEGIN' 'REF''NETMAT' JAC = JACOBIAN[L];
           STATE Z(Q);
 
           LIN RELAX(PCS,L, Q, RHS);
           (MONICS ! MONITOR("MGCSPRE   ",RHS-JAC*Q));
 
           'IF' L>1 'AND' SIGMACS>0
           'THEN' 'FIELD' QC,RC;
                  QC:= 'NUL''RESTRICT' Q;
                  RC:= 'RESTRICT' (RHS-JAC*Q);
                  'TO' SIGMACS 'DO' MGCS (1,L-1,QC,RC) 'OD';
                  Q+:= 'PROLON' QC;
 
                  (MONICS ! 'FIELD' RR = RHS-JAC*Q;
                            MONITOR("MGCSCGC1  ",RR);
                            MONITOR("MGCSCGC2  ",'RESTRICT'RR) )
           'FI';
 
           LIN RELAX(QCS,L, Q, RHS);
           (MONICS ! MONITOR("MGCSPOST  ",RHS-JAC*Q))
   'END'; #$E#
          #E$#
 
 # MULTIGRID (II)                 #                   'PR' EJECT 'PR'
 
   'INT' PMG, SIGMA, QMG;
   'BOOL' FASAB:= 'FALSE', STATE2:= 'TRUE', GALERKIN:= 'FALSE';
 
   #$H#
   'PROC' FAS = ('INT' L, 'REF'[]'FIELD'Q, RHS) 'VOID':
   #$B FAS     B$#
   'BEGIN''PROC' ADDPROL = ('FIELD' QF,QC,WC)'VOID':
          'BEGIN' 'REF'[,]'STATE' QQ  = (Q'OF'QF)[@1,@1],
                                  QQC = (Q'OF'QC)[@1,@1],
                                  QWC = (Q'OF'WC)[@1,@1];
                  'INT'   N1:= 1'UPB'QQC, N2:= 2'UPB'QQC, TI,TJ,PI,PJ;
                  'REAL'  X,Y;
                  'STATE' A,B,C,D,P,Q,R,S;
 
                  'IF' GALERKIN 'OR' 'ODD'N1 'OR' 'ODD'N2
                  'THEN''FOR' I 'TO' N1 'DO' TI:= I+I;
                        'FOR' J 'TO' N2 'DO' TJ:= J+J;
                              C := QQC[I,J] - QWC[I,J];
                              QQ[TI-1,TJ-1]+:= C;
                              QQ[TI-1,TJ  ]+:= C;
                              QQ[TI  ,TJ-1]+:= C;
                              QQ[TI  ,TJ  ]+:= C
                        'OD' 'OD'
                  'ELSE''FOR' I 'FROM' 2 'BY' 2 'TO' N1 'DO' TI:= I+I;
                        'FOR' J 'FROM' 2 'BY' 2 'TO' N2 'DO' TJ:= J+J;
                              P:= QQC[I-1,J-1] - QWC[I-1,J-1];
                              Q:= QQC[I  ,J-1] - QWC[I  ,J-1];
                              R:= QQC[I-1,J  ] - QWC[I-1,J  ];
                              S:= QQC[I  ,J  ] - QWC[I  ,J  ];
                              A:= 0.25*(P+Q+R+S);
                              B:= 0.25*(Q+S-P-R);
                              C:= 0.25*(R+S-P-Q);
                              D:= 0.25*(P+S-R-Q);
                        'FOR' K 'TO' 4 'DO' X:= K-2.5; PI:= TI-4+K;
                                            P:= X*B+A; Q :=  X*D+C;
                        'FOR' L 'TO' 4 'DO' Y:= L-2.5; PJ:= TJ-4+L;
                                            QQ[PI,PJ]+:=     Y*Q+P
                        'OD' 'OD'
                        'OD' 'OD'
                  'FI'
          'END';
 
           (MONIREL! MONITOR("RELINA    ",RHS[L]'MINRESIDU'Q[L]));
           RELAX(PMG,L, Q[L], RHS[L]);
 
           'IF'    L>1  'AND' SIGMA>0
           'THEN'
           (MONIREL! MONITOR("RELOUTA   ",RHS[L]'MINRESIDU'Q[L]));
                  ( FASAB ! STATE E(Q[L]); Q[L-1]:='RESTRICT'Q[L] );
                  'FIELD' W = 'COPY' Q[L-1];
 
                  RESIDUHO('FALSE',Q[L],RHS[L-1]:= 'COPY'RHS[L]);
                  RESIDU( 'TRUE',Q[L-1],RHS[L-1]:= 'RESTRICT'RHS[L-1]);
                 #
                  RHS[L-1]:= 'RESTRICT'RHS[L];
                  'FIELD' TAU:= 'RESIDU'Q[L-1] - 'RESTRICT''RESIDU'Q[L];
                  RHS[L-1]+:= TAU;
                 #
                  'TO' SIGMA 'DO' FAS (L-1,Q,RHS) 'OD';
 
                  ADDPROL(STATE E(Q[L]), STATE E(Q[L-1]), STATE E(W));
                  # Q[L] +:= 'PROLON2' STATE E (Q[L-1]-W);  #
 
           (MONIREL! MONITOR("RELINB    ",RHS[L]'MINRESIDU'Q[L]))
           'FI';
 
           RELAX(QMG,L, Q[L], RHS[L]);
           (MONIREL! MONITOR("RELOUTB   ",RHS[L]'MINRESIDU'Q[L]))
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' FMG = ('REF'[]'FIELD' SOL,RHS)'FIELD':
   #$B FMG     B$#
   'BEGIN' 'INT' LEVELS = 'UPB'SOL;   ('LWB'SOL /=1 ! ERROR);
           ('LWB'RHS /=1 'OR'    'UPB'RHS /= LEVELS ! ERROR);
 
           'FOR' L 'TO' LEVELS
           'DO'  RHS[L]:='NULRHS'SOL[L];
                 FAS (L, SOL, RHS);
                 (L<LEVELS ! (STATE2! STATE E(SOL[L]) );
                             SOL[L+1]:= 'PROLON2' SOL[L] );
                 ( MONIFAS ! MONITOR("FAS  ",RHS[L]'MINRESIDU'SOL[L]) )
           'OD';
           SOL[LEVELS]
   'END'; #$E#
          #E$#
 
 # MULTIGRID (III)                #                   'PR' EJECT 'PR'
 
   #$H#
   'PROC' STARFAS = ('INT' L, 'REF'[]'FIELD'Q, RHS) 'VOID':
   #$B STARFAS B$#
   'BEGIN' 'PROC'('BOOL','FIELD','FIELD')'FIELD' CGRESIDU= RESIDU;
 
           RESIDU:= STARRES1;
           (MONIREL! MONITOR("RELINA*   ",RHS[L]'MINRESIDU'Q[L]));
           STARRELAX(PMG,L, Q[L], RHS[L]);
 
           'IF'    L>1  'AND' SIGMA>0
           'THEN'
           (MONIREL! MONITOR("RELOUTA*  ",RHS[L]'MINRESIDU'Q[L]));
                  'FIELD' W = 'COPY' ( Q[L-1]:= 'RESTRICT'Q[L] );
 
                  RESIDU('FALSE',Q[L],RHS[L-1]:= 'COPY'RHS[L]);
                  RESIDU:= CGRESIDU;
                  RESIDU( 'TRUE',Q[L-1],RHS[L-1]:= 'RESTRICT'RHS[L-1]);
                  'TO' SIGMA 'DO' FAS (L-1,Q,RHS) 'OD';
                  RESIDU:= STARRES1;
 
                  Q[L] +:= 'PROLON' (Q[L-1]-W);
 
           (MONIREL! MONITOR("RELINB*   ",RHS[L]'MINRESIDU'Q[L]))
           'FI';
 
           STARRELAX(QMG,L, Q[L], RHS[L]);
           (MONIREL! MONITOR("RELOUTB*  ",RHS[L]'MINRESIDU'Q[L]));
           RESIDU:= CGRESIDU
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' STARFMG = ('REF'[]'FIELD' SOL,RHS)'FIELD':
   #$B STARFMG B$#
   'BEGIN' 'PROC'('BOOL','FIELD','FIELD')'FIELD' CGRESIDU= RESIDU;
 
           'INT' LEVELS = 'UPB'SOL;   ('LWB'SOL /=1 ! ERROR);
           ('LWB'RHS /=1 'OR'    'UPB'RHS /= LEVELS ! ERROR);
 
           'FOR' L 'TO' LEVELS
           'DO'  RHS[L]:='NULRHS'SOL[L];
                 STARFAS (L, SOL, RHS);
                 (L<LEVELS ! SOL[L+1]:= 'STARPROL' SOL[L] );
                 ( MONIFAS ! RESIDU:= STARRES1;
                             MONITOR("FAS  ",RHS[L]'MINRESIDU'SOL[L]);
                             RESIDU:= CGRESIDU )
           'OD';
           SOL[LEVELS]
   'END'; #$E#
          #E$#
 
 # STANDARD PROBLEMS.             #                   'PR' EJECT 'PR'
 
 'REAL' RR,UU,VV,PP,CC,ZZ, UU1,UU2,VV1,VV2,CC1,CC2,ZZ1,ZZ2, DIST;
 
   #$H#
 'PROC' PROBLEM1 = 'FIELD':
   #$B PROB1   B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 1"+10*" "+"VSN. 850515", NEWLINE,
                "BUMP, TRANSSOON", NEWLINE,
                "EQUIDISTANT"));
   'FIELD' Q1 = COMPFIELD (1.0,1.0, -2,2, 0,2);
   MAPPING:= ('REAL' XX,YY)'POINT': (XX+0.5,YY);
   WALLS:= WALLSR; FLUXT:= FLUXTR;
 
   OOST :=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
            WALL (W, 'FALSE', 0, QIN,QAC) ;
   WEST :=('WALL'W, 'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           ('REAL' X = X'OF'MAPPING(XX,YY);
            WALL (W, 'TRUE' , ('ABS'X>0.5!0!8.0*0.042*X), QIN,QAC)
           );
   ZUID :=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
            SUB OUT PRESS(W,'FALSE',PP,QIN,QAC);
   NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
            SUB IN UVZ  (W,'TRUE',UU,VV,ZZ,QIN,QAC);
 
 
   RR := 1.0; UU:= 0.85; VV:= 0; PP:= 1/GAM;
   CC := SQRT(GAM*PP/RR);
   ZZ := LN(PP)-GAM*LN(RR);
 
   INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
   PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
   #$H#
 'PROC' PROBLEM2 = 'FIELD':
   #$B PROB2   B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 2"+10*" "+"VSN. 850515", NEWLINE,
                "BUMP, TRANSSOON", NEWLINE,
                "NON-EQUIDISTANT"));
   'FIELD' Q1 = COMPFIELD (1.0,1.0, -2,3, 0,2);
 
   OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'FALSE', 0, QIN,QAC) ;
   ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB OUT PRESS(W,'FALSE',PP,QIN,QAC);
  NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB IN UVZ  (W,'TRUE',UU,VV,ZZ,QIN,QAC);
 
   MAPPING:= ('REAL' XX,YY) 'POINT':
   'BEGIN''REAL' X,Y;
          ( XX<=-1.37 ! X:= 2.27*XX+2.54;
                        X:=-0.57-0.15*(EXP(-1.650*(X+0.57))-1.0)
          !:XX>= 2.15 ! X:= 2.88*XX-5.64;
                        X:= 0.55+0.11*(EXP( 1.286*(X-0.55))-1.0)
                      ! X:= 0.32*XX-0.14  );
          Y:= (YY>1.99 !2.0 !:YY<0.01 !0! 0.19*(EXP(1.222*YY)-1.0) );
          'POINT'(X, ('ABS'X<=0.5! Y+:= 0.084*(0.25-X*X)*(2.0-Y) ! Y))
   'END';
 
   WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'TRUE' , 0, QIN,QAC) ;
   WALLS:= WALLSG; FLUXT:= FLUXTG;
 
       RR  :=  GAM; UU:= 0.85; VV:= 0; PP:= 1.05;
       CC := SQRT(GAM*PP/RR);
       ZZ := LN(PP)-GAM*LN(RR);
     INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
   PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
   #$H#
 'PROC' PROBLEM3 = 'FIELD':
   #$B PROB3   B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 3"+10*" "+"VSN. 850515", NEWLINE,
                "CONTACT DISCONTINUITY ", NEWLINE,
                "EQUIDISTANT"));
   'FIELD' Q1 = COMPFIELD (1.0,1.0,  0,2, 0,2);
   MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY);
 
          PP:=1.0; UU1:=0.6; VV1:=-0.6; CC1:=1.0;
                   UU2:=0.3; VV2:=-0.3; CC2:=0.6;
   'REAL' RR1:=GAM*PP/(CC1*CC1),     RR2:=GAM*PP/(CC2*CC2);
          ZZ1:=LN(PP*(RR1**(-GAM))); ZZ2:=LN(PP*(RR2**(-GAM)));
          RR :=(RR1+RR2)/2.0; UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0;
 
   WALLS:= WALLSR; FLUXT:= FLUXTR;
 
   OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB IN UVZ  (W,'FALSE',UU1,VV1,ZZ1,QIN,QAC);
   WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB OUT PRESS(W,'TRUE',PP,QIN,QAC) ;
   ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB OUT PRESS(W,'FALSE',PP,QIN,QAC);
  NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB IN UVZ  (W,'TRUE',UU2,VV2,ZZ2,QIN,QAC);
 
   INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
   PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
 'REAL' PHI0;
 
   #$H#
 'PROC' PROBLEM3B = 'FIELD':
   #$B PROB3B  B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 3B"+10*" "+"VSN. 860217", NEWLINE,
                "CONTACT DISCONTINUITY,PHI=PI*"+FIXED(PHI0/PI,6,2),
                NEWLINE,"EQUIDISTANT"));
   'FIELD' Q1 = COMPFIELD (1.0,1.0,  0,2, 0,2);
   MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY);
 
          PP:=1.0; CC1:=0.6;
                   UU1:=0.3*SQRT(2.0)*SIN(PHI0);
                   VV1:=0.3*SQRT(2.0)*COS(PHI0);
          PP:=1.0; CC2:=1.0;
                   UU2:=0.6*SQRT(2.0)*SIN(PHI0);
                   VV2:=0.6*SQRT(2.0)*COS(PHI0);
 
 
         F0:= (('REAL' X,Y)'REAL': 0.0 );
         F1:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!UU1!UU2)*(1+F0(X,Y));
         F2:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!VV1!VV2)*(1+F0(X,Y));
         F3:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!CC1!CC2)*(1+F0(X,Y));
         F4: =('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!ZZ1!ZZ2)*(1+F0(X,Y));
 
         'REAL' RR1:= GAM*PP/(CC1*CC1); ZZ1:= LN(PP)-GAM*LN(RR1);
         'REAL' RR2:= GAM*PP/(CC2*CC2); ZZ2:= LN(PP)-GAM*LN(RR2);
          RR :=(RR1+RR2)/2.0; UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0;
 
         WALLS:= WALLSR; FLUXT:= FLUXTR;
 
         OOST:= WEST:= ZUID:= NOORD:=
         ('WALL' W,'REAL'X,Y,'STATE'QIN,'REF''STATAC'QAC)'STATE':
         FIXED UVCZ (W, F1(X,Y),F2(X,Y),F3(X,Y),F4(X,Y) ,QIN,QAC);
 
         INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
         PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
   #$H#
 'PROC' PROBLEM4 = 'FIELD':
   #$B PROB4   B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 4"+10*" "+"VSN. 850515", NEWLINE,
                "GLADDE SUBSONE STROMING", NEWLINE,
                "EQUIDISTANT"));
   'FIELD' Q1 = COMPFIELD (0.5,0.5, -2,2, 0,2);
   MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY);
 
   OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'FALSE', 0, QIN,QAC) ;
   ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB OUT PRESS(W,'FALSE',PP,QIN,QAC);
  NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB IN UVZ  (W,'TRUE',UU,VV,ZZ,QIN,QAC);
   WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'TRUE' , DIST*SIN(PI*XX), QIN,QAC);
   WALLS:= WALLSR; FLUXT:= FLUXTR;
 
     DIST:=0.020;
     RR  :=  1.0; UU:= 0.75; VV:= 0; PP:= 1/GAM;
     CC := SQRT(GAM*PP/RR);
     ZZ := LN(PP)-GAM*LN(RR);
 
     INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
   PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
   #$H#
 'PROC' PROBLEM5 = 'FIELD':
   #$B PROB5   B$#
 'BEGIN' PRINT(("STANDAARDPROBLEEM 5"+10*" "+"VSN. 850731", NEWLINE,
                "REFLECTERENDE SCHEVE SCHOKGOLF IN KANAAL", NEWLINE,
                "EQUIDISTANT ALS BIJ OSHER"));
   'FIELD' Q1 = COMPFIELD (4/6,1/2,  0,6, 0,2);
 
   MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY);
 
  NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           FIXED UVCZ (W,UU1,VV1,CC1,ZZ1,QIN,QAC);
   WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'TRUE' , 0, QIN,QAC) ;
   OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           FIXED UVCZ (W,UU2,VV2,CC2,ZZ2,QIN,QAC);
   ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           SUB OUT PRESS (W,'FALSE',PP,QIN,QAC);
   WALLS:= WALLSR; FLUXT:= FLUXTR;
 
  # TESTPROBLEEM: OBLIQUE SHOCK,
    EXACTE OPLOSSING:
 
   'REAL' U1:=2.9,V1:=0.0,C1:=1.0,Z1:=-0.47,P1:=1.0,R1:=1.4,S1:=0.62;
   'REAL' U2:=2.62,V2:=-0.51,C2:=1.12,Z2:=-0.45,P2:=2.14,R2:=2.38,
          S2:=0.64;
   'REAL' U3:=2.41,V3:=0.0,C3:=1.23,Z3:=-0.45,P3:=4.0,R3:=3.70,S3:=0.64;
  #
   UU1:=2.9 ; VV1:= 0.0 ; CC1:=1.0 ; ZZ1:=-0.47;
   UU2:=2.62; VV2:=-0.51; CC2:=1.12; ZZ2:=-0.45;
   PP :=4.0 ;
   UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0;
   CC:=(CC1+CC2)/2.0; ZZ:=(ZZ1+ZZ2)/2.0;
 
   INITIATE('REPORT'Q1,'EZ''STATE'(UU,VV,CC,ZZ),1);
   PRINT((NEWLINE,NEWLINE)); Q1
   'END'; #$E#
          #E$#
 
#?$H?
 'PROC' PROBLEM6 = 'FIELD':
 ?$B  PROB5 B$?
 ? NOG NIET GEBRUIKT TESTPROBLEEM, 850513 ?
 'BEGIN' PRINT(("STANDAARDPROBLEEM 6"+10*" "+"VSN. 850515", NEWLINE,
                "REFLECTERENDE SCHEVE SCHOKGOLF IN KANAAL", NEWLINE,
                "NIET-EQUIDISTANTE INLAAT, EQUIDISTANTE UITLAAT"));
   'FIELD' Q1 = COMPFIELD (1.0,1.0,  0,8, 0,1);
 
   MAPPING:= ('REAL' XX,YY) 'POINT':
   ( XX,  ( XX<1.0 ! YY*(1.0+0.525*(1.0-XX)) ! YY ) );
 
  NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           FIXED UVCZ (W,UU,VV,CC,ZZ,QIN,QAC);
   WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'TRUE' , 0, QIN,QAC) ;
   OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
           WALL (W, 'FALSE', 0, QIN,QAC) ;
   WALLS:= WALLSG; FLUXT:= FLUXTG;
 
   RR:= 1.0; UU:= 2.0; VV:= 0.0; PP:= 1.0/GAM;
   CC:= SQRT(GAM*PP/RR);
   ZZ:= LN(PP)-GAM*LN(RR);
   INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1);
   PRINT((NEWLINE,NEWLINE)); Q1
 'END';?$E?
      ?E$?
#
 
 # RESULTAAT ANALYSE              #                   'PR' EJECT 'PR'
 
 #$H#
 'PROC' PRINTORDER = ([]'INT' NOS )'VOID':
 #$B PRNT.OR B$#
 'BEGIN' ('UPB' NOS /=3 ! ERROR);
         [1:3]'FIELD' ALLF;
         PRINT((NEWLINE,"OBSERVED ORDER OF ACCURACY"));
         'FOR' L 'TO' 3 'DO' TAKE( NOS[L], ALLF[L] ) 'OD';
 
         ALLF[3]:= 'RESTRICT' ALLF[3];
         ALLF[2]-:= ALLF[3];
         ALLF[1]-:= 'RESTRICT' ALLF[3];
         ALLF[1]-:= 'RESTRICT' ALLF[2];
         [1:4 ]'REAL' DEN:= 'NORM1'ALLF[2],
                      NUM:= 'NORM1'ALLF[1];
         PRINT((NEWLINE,"ORDER="));
         'FOR' I 'TO'4
         'DO' 'REAL' FACTOR := NUM[I]/DEN[I];
              PRINT((" "+WHOLE(I,0)+" ("+FIXED(FACTOR,6,2)+
               ") "+(FACTOR>0!FIXED(LN(FACTOR)/LN(2.0),6,2)!"******") ))
         'OD'
 'END'; #$E#
        #E$#
 
    #$H#
   'PROC' TOTVAR = ('FIELD' QF )[,]'REAL':
   #$B TOTVAR  B$#
   'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF;
           'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF;
           'INT'  IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q;
           [1:4,0:1]'REAL' NRM;
 
           'FIELD' QR = (DDX,DDY,'LOC''INT':=1,Q[ IL+1:IU'AT'1,]),
                   QL = (DDX,DDY,'LOC''INT':=1,Q[ IL:IU-1'AT'1,]),
                   QU = (DDX,DDY,'LOC''INT':=1,Q[,JL+1:JU'AT'1 ]),
                   QD = (DDX,DDY,'LOC''INT':=1,Q[,JL:JU-1'AT'1 ]);
           [,]'REAL' NRMH = 'NORM' (QR-QL),
                     NRMV = 'NORM' (QU-QD);
           'FOR' C 'TO' 4
           'DO' NRM[C,1]:= (NRMH[C,1]+NRMV[C,1])*
                           0.5*SQRT((IU-IL+1)*(JU-JL+1)) ;
                NRM[C,0]:= (NRMH[C,0]<NRMV[C,0]!NRMH[C,0]!NRMV[C,0])
           'OD'; NRM
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'PRINTTV' = ('FIELD' F)'VOID':
   #$B PRINTTV B$#
   'BEGIN' [,]'REAL' N = TOTVAR( F );
           PRINT((NEWLINE,"TOT.VAR.:   ",SF(N[,1]),
                  NEWLINE,"MAX.VAR.:   ",SF(N[,0]),NEWLINE ))
   'END' ; #$E#
           #E$#
 
 # OUTPUT ROUTINES                #                   'PR' EJECT 'PR'
 
   #$H#
   'OP' 'REPORT' = ('FIELD' F)'FIELD':
   #$B REPORT  B$#
   'BEGIN' 'REAL' DDX = DDX'OF'F, DDY = DDY'OF'F;
           'REF'[,]'STATE' Q = Q'OF'F;
           'INT'  L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
           'STRING' LS2= WHOLE(L2,0), US2= WHOLE(U2,0),
                    LS3= WHOLE(L3,0), US3= WHOLE(U3,0);
           'POINT' NW= MAPPING(XXN,YYW), NO= MAPPING(XXN,YYO),
                   ZW= MAPPING(XXZ,YYW), ZO= MAPPING(XXZ,YYO);
            PRINT((NEWLINE,"VELD: [ " + LS2 + " : " + US2 +" ] * [ "+
                                        LS3 + " : " + US3 +" ]     ",
            NEWLINE,"X IN ("+ ST((L2-1)*DDX)+ST(U2*DDX) + " )  "+
                    "  IN ("+ ST(  'POINT'((XXN,XXZ))  )+ " )  "+
                    "  NW ("+ ST(         NW           )+ " )  "+
                    "  NO ("+ ST(         NO           )+ " )  ",
            NEWLINE,"Y IN ("+ ST((L3-1)*DDY)+ST(U3*DDY) + " )  "+
                    "  IN ("+ ST(  'POINT'((YYW,YYO))  )+ " )  "+
                    "  ZW ("+ ST(         ZW           )+ " )  "+
                    "  ZO ("+ ST(         ZO           )+ " )  "
                                                                )); F
   'END';#$E#
         #E$#
 
   #$H#
   'OP' 'PRINT' = ('FIELD' F)'VOID':
   #$B PRNT.1  B$#
   'FOR' K 'TO' 4
   'DO' K'PRINT'F
   'OD' ; #$E#
          #E$#
 
   #$H#
   'OP' 'PRINT' = ('INT' K,'FIELD' F)'VOID':
   #$B PRNT.2  B$#
   'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
           'STRING' S:= ( TYPE'OF'F <= 1
                ! (K!"MASS","X-MOMT","Y-MOMT","ENERGY")
                ! (K!"X-SPD","Y-SPD","SND-SPD","LN(ENTROPY)") );
           ( TYPE'OF'F = 0 ! S+:= "-FLUX");
           ( TYPE'OF'F = 1 ! S+:= "-DENSITY");
           PRINT((NEWLINE,S));
           'FOR' I 'TO' 1'UPB'Q
           'DO' PRINT((NEWLINE, ST( []'REAL' ((
           (K!C1'OF'Q[I,],C2'OF'Q[I,],C3'OF'Q[I,],C4'OF'Q[I,])
                              ))  )
                     ))
           'OD'
   'END'; #$E#
          #E$#
 
   #$H#
   'OP' 'PRINT' = ('PROC'('STATE')'REAL' P,'FIELD' F)'VOID':
   #$B OPPRINT B$#
   'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
                               STATE E (F);
           'FOR' I 'TO' 1'UPB'Q
           'DO' 'STRING' S:= "";
                'FOR' J 'TO' 2'UPB'Q 'DO' S+:=ST(P(Q[I,J])) 'OD';
                PRINT((NEWLINE,S))
           'OD'
   'END'; #$E#
          #E$#
 
   'OP' 'FIX' = ('REAL'X)'STRING' : FIXED(X,7,3  ) ;
   'OP' 'FLT' = ('REAL'X)'STRING' : FLOAT(X,7,3,1) ;
 
   'MODE' 'UNI'   = 'UNION'('REAL','POINT','WALL',[]'REAL','STATE' ) ;
 
   #$H#
   'PROC' ST = ('UNI' U)'STRING':
   #$B ST.UNI  B$#
     'CASE' U 'IN'
       ('REAL'   X): " "+'FIX' X ,
       ('POINT'  P): " "+'FIX'X'OF'P+" "+'FIX'Y'OF'P ,
       ('WALL'   W): " " + 'FIX' C'OF'W + " " + 'FIX' S'OF'W +
                     " " + 'FIX'DS'OF'W ,
       ([]'REAL' A): ('STRING'S:="";'FOR'I'FROM''LWB'A'TO''UPB'A
                        'DO' S+:=" "+'FIX' A[I] 'OD' ; S ) ,
       ('STATE'  S): " " + 'FIX'C1'OF'S + " " + 'FIX'C2'OF'S +
                     " " + 'FIX'C3'OF'S + " " + 'FIX'C4'OF'S
     'ESAC' ; #$E#
              #E$#
 
   #$H#
   'PROC' SF = ('UNI' U)'STRING':
   #$B SF.UNI  B$#
     'CASE' U 'IN'
       ('REAL'   X): " "+'FLT' X ,
       ('POINT'  P): " "+'FLT'X'OF'P+" "+'FLT'Y'OF'P ,
       ('WALL'   W): " " + 'FLT' C'OF'W + " " + 'FLT' S'OF'W +
                     " " + 'FLT'DS'OF'W ,
       ([]'REAL' A): ('STRING'S:="";'FOR'I'FROM''LWB'A'TO''UPB'A
                        'DO' S+:=" "+'FLT' A[I] 'OD' ; S ) ,
       ('STATE'  S): " " + 'FLT'C1'OF'S + " " + 'FLT'C2'OF'S +
                     " " + 'FLT'C3'OF'S + " " + 'FLT'C4'OF'S
     'ESAC' ; #$E#
              #E$#
 
   #$H#
   'PROC' FLUXNORM = ('FIELD'F,'REF''FLUX' NORM,RATIO)'REAL':
   #$B FLUXXN  B$#
   'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
       'INT' N1=1'UPB'Q, N2=2'UPB'Q;
       'FOR' K 'TO' 4
       'DO' 'REF'[,]'REAL' QK=(K!C1'OF'Q,C2'OF'Q,C3'OF'Q,C4'OF'Q);
            'REAL' T,S:= 0, S1:= 0;
            'FOR' I 'TO' N1 'DO' 'FOR'J 'TO' N2 'DO'
                  T:= 'ABS' QK[I,J]; S1+:=T; (T>S ! S:=T )
            'OD' 'OD';
            (K!C1'OF'NORM ,C2'OF'NORM ,C3'OF'NORM ,C4'OF'NORM ) := S1;
            (K!C1'OF'RATIO,C2'OF'RATIO,C3'OF'RATIO,C4'OF'RATIO) :=
               S*N1*N2/S1
       'OD';'ABS' NORM
   'END';#$E#
         #E$#
 
   #$H#
  'PROC' ACCOUNT OK = 'BOOL':
   #$B ACCOU   B$#
  'BEGIN' 'REAL' E,D,C:= CLOCK;
          D:=C-TIME;
          E:=TIMELIMIT-C;
          'BOOL' THROUGH = E > 1.2*D; TIME:= C;
          THROUGH
   'END'; #$E#
          #E$#
 
   #$H#
   'PROC' HISTO = ('PROC'('STATE')'REAL' P, 'REF'[]'STATE' S)'VOID':
   #$B HISTO   B$#
   'BEGIN' 'INT' L= 'LWB'S, U='UPB'S;
           [L:U]'REAL'X;
           'FOR'I'FROM'L'TO'U 'DO'X[I]:=P(S[I])'OD';
           'REAL' FACTOR, MAX:= 0, MIN:= 0;
           'INT' K, W:= (U-L<30!40!:U-L<60!80!120);
           'FOR'I'FROM'L'TO'U
           'DO'(X[I]>MAX!MAX:=X[I]!: X[I]<MIN! MIN:=X[I]) 'OD';
           FACTOR:= W/(MAX-MIN);
           PRINT((NEWLINE,NEWLINE,"RANGE = ",MIN," TOT ",MAX));
           'FOR'I'FROM'L'TO'U
           'DO' K:= 'ROUND'( (X[I]-MIN)*FACTOR );
                PRINT((NEWLINE,K*" ","*"))
           'OD'
   'END'; #$E#
          #E$#
 
  'PROC' MONITOR := ('STRING' T,'FIELD' F)'VOID':
  'BEGIN' PRINT((NEWLINE,">>>    MONITOR STILL UNDEFINED    <<<",
                 NEWLINE,">>> ADVISE : DO MONITOR := MONIT1 <<<",
                 NEWLINE)); ERROR ; 'SKIP'
  'END' ;
 
   #$H#
   'PROC' MONIT1  = ('STRING' TEXT, 'FIELD' R)'VOID':
   #$B MONIT1  B$#
   'BEGIN' 'FLUX' NORM,RATIO; FLUXNORM(R,NORM,RATIO);
           'REAL' S =
           ( C1'OF'NORM + C2'OF'NORM + C3'OF'NORM + C4'OF'NORM ) /4;
           PRINT((NEWLINE, TEXT+"  ",
                  ST([]'REAL'((DDX'OF'R,DDY'OF'R))),SPACE,
                  FLOAT(S,10,3,3)+" <   "+
                  SF(NORM)+ST(RATIO) ))
   'END'; #$E#
          #E$#
 
   'BOOL' MONIREL:= 'FALSE', MONICS := 'FALSE',
          MONINWT:= 'FALSE',
          MONILIN:= 'FALSE', MONIFAS:= 'FALSE',
          MONI   := 'FALSE', MONIGS := 'FALSE';
 
   'REAL' TIME:= CLOCK;
   PRINT((">> EULER A68-LIBRARY  VSN.860414 ",NEWLINE,
          ">> FLUX ETC. ARE STRUCTS (JR)    ",NEWLINE,
                                      NEWLINE,NEWLINE)) ;
   'PR' PROG 'PR' 'SKIP' ;
   (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE)) )
'END'
LIBRARY(EULIB,NEW)
ADD(*,LGO)
FINISH.
ENDRUN.
LIBRARY(EULIB,OLD)
ADD(*,LGO)
FINISH.
ENDRUN.
'BEGIN'  # EUTEST           VSN.850724    PWH/840305    #
         # REFERENTIE OPLOSSING, EERSTE ORDE, PROB. 4   #
         # GLAD TESTPROBLEEM                            #
 
   RESIDU    := RESIDU1            ;
   RESIDU2   := CENTRALX           ;
   NFLUX     := OSHER ; SOS:=-1    ;
   GALERKIN  := 'TRUE'             ;
   SYMMETRIC := 'TRUE'             ;
   FASAB     := 'FALSE'            ;
   PMG:=1;QMG:=1;SIGMA:=1          ;
   RELAX     := SEIDEL             ;
   NLGSTOL   := 0.1                ;
 
 'PROC' SHOW = ('INT' L)'VOID':
 'IF' L=LEVELS
 'THEN'  'FIELD' RR:= Q[L];
         'FLUX' NORM,RATIO;
         'REAL' N;
         PRINT((NEWLINE,NEWLINE,
                "OPLOSSING OP NIVEAU ",WHOLE(L,0) ));
         'PRINT' RR;
 
                 RR:= ADD WALL STATES('FALSE',Q[L]);
         STATE E (RR);
         PRINT((NEWLINE,NEWLINE,"PRESSURE",NEWLINE));
         PRESSURE'PRINT'RR; Q'OF'RR:= 'NIL';
 
                 RR := 'RESIDU' Q[L];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "RESIDU1 OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
                 RR := 'RESIDU2' Q[L];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "RESIDU2 OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
                 RR:= Q[L] 'CORRECT''NUL'R[L];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "CORRECTIE OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
                 RR:= 'RESTRICT' RR;
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "RESTRICTBAR VAN CORRECTIE OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
   (L>1 !        RR := 'TAU1' Q[L];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "TAU1 OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
                 RR:= 'TAU2' Q[L];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "TAU2 OP NIVEAU ",WHOLE(L,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N));
 
                 RR := R[L-1];
         FLUXNORM(RR,NORM,RATIO);
         PRINT((NEWLINE,NEWLINE,
                "RHS OP NIVEAU ",WHOLE(L-1,0),
                "   ", N:= 'ABS'NORM));
         PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0)));
         'PRINT' (RR*:= 10.0**(1-N))
   )
 'FI';
 
 'INT' LEVELS = 3;
 [1:LEVELS]'FIELD'Q,R; Q[1]:= PROBLEM4;
 
 'FOR' L 'TO'LEVELS
 'DO'  R[L]:= 'NULRHS'Q[L];
 
       'TO' 4
       'WHILE' MONITOR("FAS       ",'RESIDU'(Q[L])); ACCOUNT OK
       'DO' FAS(L, Q, R) 'OD';
 
       SHOW(L);
 
       (L<LEVELS
       !PRINT((NEWPAGE));
        Q[L+1]:= 'PROLON2' STATE E(Q[L])
       )
 'OD'
 
'END'

« March 2024 »
Su Mo Tu We Th Fr Sa
1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31
 

Powered by Plone CMS, the Open Source Content Management System

This site conforms to the following standards: