//JOB(ERCC07,STEPHENS,T=200,L=10000)
//OPTIONS(QUOTES,OPT)
//ALGOL
'COMMENT'
                         ***************
                         *             *
                         *   ALGG047   *
                         *             *
                         ***************
;
 'BEGIN'
'BEGIN'
'COMMENT' A PROGRAM TO FIND THE STARTING VALUE OF HMIN;
'COMMENT' A30;
      'REAL' A,B,BE,BO,C,CR,D,DU,EL,EM,EN,EX,ESS,         FN,FN2,FN3,HO,
           H1,H2,OMEGA,P1,P2,RT,RU1,RU2,S,VE,VEL,WY,X,DIF,BIF,INT;
      'INTEGER' I,J,N;
      'ARRAY' HMIN[1:20];
      'REAL' 'PROCEDURE' SIMP(P,Q,PX,N);
           'VALUE' P,Q,N;
           'REAL' P,Q;
           'INTEGER' N;
           'REAL' 'PROCEDURE' PX;
'COMMENT' (X):'VALUE' X: 'REAL' X;
           'BEGIN'
                 'INTEGER' R; 'REAL' Z,DX;
                 DX:=(P-Q)/N; Z:=PX(P)+PX(Q);
                 'FOR' R:=1'STEP' 2 'UNTIL' N-1 'DO' Z:=Z+4*PX(P+DX*R);
                 'FOR' R:=2'STEP' 2 'UNTIL' N-2 'DO' Z:=Z+2*PX(P+DX*R);
                 SIMP:=Z*DX/3;
                PRINT(Z*DX/3,0,3);  NEWLINES(1);
           'END';
      'REAL' 'PROCEDURE' BX(X);
           'VALUE' X; 'REAL' X;
           'BEGIN'
                  FN:=A*X**2+B*X+C; FN2:=FN**2; FN3:=FN**3;
                 D:=2*A*X+B;
                    'IF' S<0 'THEN' BX:=(2/RT)*ARCTAN(D/RT)
                 'ELSE' 'IF' S=0 'THEN' BX:=-2/D
                 'ELSE' BX:=(1/RT)*LN(ABS((D-RT)/(D+RT)));
           'END';
      'REAL' 'PROCEDURE' GX(X);
           'VALUE' X; 'REAL' X;
           'BEGIN'
                 GX:=D/2*(-S)*FN2-3*A*D/S**2*FN+(6*A**2/S**2)*BX(X);
           'END';
      'REAL' 'PROCEDURE' HX(X);
           'VALUE' X; 'REAL' X;
           'BEGIN'
'REAL' Y;
Y:=BX(X);
                 HX:=D/(-S)*FN +(2*A/(-S))*Y;
           'END';
      'REAL' 'PROCEDURE' FX(X);
           'VALUE' X; 'REAL' X;
           'BEGIN'
'REAL' Y;
Y:=HX(X);
                 FX:=(B*X+2*C)/3*(-S)*FN2-(3*B/2*(-S))*Y;
           'END';


      HO:=READ; H1:=READ; H2:=READ; BE:=READ; P1:=READ;
      P2:=READ; OMEGA:=READ; VEL:=READ; CR:=READ;
      N:=READ;
      'FOR' I:=1 'STEP' 1 'UNTIL' 6 'DO' HMIN[I]:=READ;

      'FOR' J:=1 'STEP' 1 'UNTIL' 6 'DO'
           'BEGIN'
                 A:=(H1+H2-2*HO)/2*BE**2;
                 B:=(H1-H2)/2*BE;
                 C:=HO+HMIN[J];
                 'IF' ABS(-B/2*A)>BE 'THEN' BO:=-BE 'ELSE' BO:=-B/2*A;
                 S:=B**2-4*A*C;
                 RT:=SQRT(ABS(S));
                 EL:=FX(BE)-FX(BO); EM:=GX(BE)-GX(BO); EN:=HX(BE)-HX(BO)
;
                 DIF:=(BE-BO);
                 INT:=SIMP(BE,BO,GX,N);
                 BIF:=DIF*GX(BE)-INT;
                 EX:=SIMP(BE,BO,FX,N)-DIF*FX(BE)+(EL/EM)*BIF;
                 WY:=SIMP(BE,BO,HX,N)-DIF*HX(BE)+(EN/EM)*BIF;
                 VE:=(INT-DIF*GX(BO))/EM;
                 DU:=(DIF*GX(BE)-INT)/EM;
                 ESS:=VEL/OMEGA*CR;
                 RU1:=ESS*WY-P1*VE-P2*DU-P1*(BE+BO)+2*P1*BE;
       LOOP:     N:=2*N;
                 INT:=SIMP(BE,BO,GX,N);
                 BIF:=DIF*GX(BE)-INT;
                 EX:=SIMP(BE,BO,FX,N)-DIF*FX(BE)+(EL/EM)*BIF;
                 WY:=SIMP(BE,BO,HX,N)-DIF*HX(BE)+(EN/EM)*BIF;
                 VE:=(INT-DIF*GX(BO))/EM;
                 DU:=(DIF*GX(BE)-INT)/EM;
                 RU2:=ESS*WY-P1*VE-P2*DU-P1*(BE+BO)+2*P1*BE;
                 'IF' ABS((RU2-RU1)/RU2) > 1 'THEN'
                 'BEGIN'
          PRINT(RU1,0,3);
          PRINT(RU2,0,3);
                        RU1:=RU2;
                        'GOTO' LOOP
                 'END';
      WRITETEXT('(''('3S')'HMIN')'); WRITETEXT('(''('19S')'RESIDUAL'('C'
)'')');
                 PRINT(HMIN[J],0,3);
                 SPACES(15);   PRINT(RU2,0,3);   NEWLINES(1);
           'END' OF FOR STATEMENT;
      WRITETEXT('(''('S')'NO.OF STEPS=')');
PRINT(N,0,3);  NEWLINES(1);

'END'
 'END';
//RUN
0.4@-3;1.3@-3;0.00;138.3@-3;200.00;116.00;52.5;460.00;7.25;
10.00;1@-6;1@-5;1@-4;1@-3;1@-2;1@-1;