"PROCEDURE" FIND GRADIENTS AT NODES(X,Y,GRAD,LBX,UBX);
  "VALUE" LBX,UBX;   "INTEGER" LBX,UBX;
  "ARRAY" X,Y,GRAD;
"COMMENT" X AND Y ARE ARRAYS WITH SUBSCRIPTS FROM LBX
   TO UBX GIVING THE ABSCISSAE AND ORDINATES RESPECT-
   IVELY OF THE DATA POINTS.  THESE SHOULD BE IN
   ASCENDING ORDER OF ABSCISSAE, WITH NO TWO ABSCISSAE
   EQUAL.  THE ARRAYS MUST CONTAIN AT LEAST FOUR
   POINTS.
   GRAD IS AN ARRAY HAVING THE SAME SUBSCRIPTS AS X
   AND Y INTO WHICH WILL BE PLACED THE CALCULATED
   GRADIENTS OF THE REQUIRED CUBICS AT EACH OF THE
   DATA POINTS.
   LBX AND UBX ARE INTEGERS GIVING THE MINIMUM AND
   MAXIMUM VALUES RESPECTIVELY OF THE SUBSCRIPTS OF
   THE ARRAYS X, Y, AND GRAD;
"BEGIN"
  "INTEGER" I,ILESS2,ILESS1,IPLUS1,IPLUS2;
  "REAL" X0,X1,X2,X3,X4,Y2,PROD1,PROD2,NUM,DENOM,G,
     COEFF2,XDIFF,XPROD,WEIGHT;
  "FOR" I:= LBX "STEP" 1 "UNTIL" UBX "DO"
  "BEGIN"
"COMMENT" SPECIAL TREATMENT IS NEEDED AT END POINT;
     ILESS1:= "IF" I>LBX "THEN" I-1 "ELSE" I+3;
     IPLUS1:= "IF" I<UBX "THEN" I+1 "ELSE" I-3;
     X2:= X[I];    Y2:= Y[I];
     X1:= X[ILESS1]-X2;    X3:= X[IPLUS1]-X2;
"COMMENT" FIRST FIT A QUADRATIC THROUGH X1,X2,X3;
     PROD1:= (Y[ILESS1]-Y2)*X3;
     PROD2:= (Y[IPLUS1]-Y2)*X1;
     DENOM:= X1*X3*(X[ILESS1]-X[IPLUS1]);
     G:= (X1*PROD2-X3*PROD1)/DENOM;
     COEFF2:= (PROD1-PROD2)/DENOM;
"COMMENT" IF X0 EXISTS, FIND ITS CONTRIBUTION TO THE
   CUBIC ADJUSTMENT;
    "IF" I"LE"LBX+1 "THEN" NUM:= DENOM:= 0.0 "ELSE"
    "BEGIN"
       ILESS2:= I-2;    X0:= X[ILESS2]-X2;
       XDIFF:= X[ILESS2]-X[ILESS1];
       XPROD:= X0*XDIFF*(X[ILESS2]-X[IPLUS1]);
       WEIGHT:= XPROD/(XDIFF*XDIFF);
       NUM:= WEIGHT*(Y[ILESS2]-Y2-X0*(G+X0*COEFF2));
       DENOM:= WEIGHT*XPROD
    "END";
"COMMENT" IF X4 EXISTS, FIND ITS CONTRIBUTION TO THE
   CUBIC ADJUSTMENT;
    "IF" I<UBX-1 "THEN"
    "BEGIN"
       IPLUS2:= I+2;    X4:= X[IPLUS2]-X2;
       XDIFF:= X[IPLUS2]-X[IPLUS1];
       XPROD:= X4*XDIFF*(X[IPLUS2]-X[ILESS1]);
       WEIGHT:= XPROD/(XDIFF*XDIFF);
       NUM:= NUM+WEIGHT*(Y[IPLUS2]-Y2-X4*(G+X4*COEFF2));
       DENOM:= DENOM+WEIGHT*XPROD
    "END";
     GRAD[I]:= G+NUM*X1*X3/DENOM
  "END"
"END" OF PROCEDURE FIND GRADIENTS AT NODES;
"PROCEDURE" COEFFS(X,Y,GRAD,I,C0,C1,C2,C3);
  "VALUE" I;   "INTEGER" I;
  "REAL" C0,C1,C2,C3;
  "ARRAY" X,Y,GRAD;
"COMMENT" THIS PROCEDURE CALCULATES THE COEFFICIENTS
   OF THE CUBIC (IN THE RANGE X[I] TO X[I+1]) WHICH
   HAS VALUES Y[I], Y[I 1] AND GRADIENTS GRAD[I],
   GRAD[I+1] AT X[I], X[I+1] RESPECTIVELY.
   THE CUBIC TAKES THE FORM
     C0 + C1*(X-X[I]) + C2*(X-X[I])'2 + C3*(X-X[I])'3;
"BEGIN"
  "REAL" H,DY;
  C0:= Y[I];
  H:= X[I+1]-X[I];    DY:= Y[I+1]-C0;
  C1:= GRAD[I];
  C2:= (3.0*DY-H*(2.0*C1+GRAD[I+1]))/(H*H);
  C3:= (H*(C1+GRAD[I+1])-2.0*DY)/H'3
"END" OF PROCEDURE COEFFS;

.