FTN,L,B 
C 
C HP92407A INTEGRATION PACKAGE
C 
C SOURCE TAPE 92407-80001 REV. A
C RELOC. TAPE 92407-60001 REV. A
C 
C AUTHOR - T.A. SAPONAS 
C 
C VERSION 2 OCTOBER 1973
C 
C 
      SUBROUTINE STRTA(A,XDLTA,Y1)
C "STRTA" INITIALIZES ARRAY "A" FOR THE START OF RUNNING INTEGRATION. 
C  THE CONTENTS OF "A" ARE AS FOLLOWS:
C   A(1) - PRESENT VALUE OF NUMERIC INTEGRAL
C   A(2) - 1ST ORDER INTEGRATION - STEP SIZE/2. 
C          2ND ORDER INTEGRATION - MINUS STEP SIZE/12.
C   A(3) - PREVIOUS VALUE OF INTEGRAND
C   A(4) - NOT NEEDED FOR 1ST ORDER INTEGRATION 
C          2ND ORDER INTEGRATION - PREVIOUS VALUE BEFORE A(3) 
      DIMENSION A(4)
C INITIALIZE INTEGRAL 
      A = 0.
C INITIALIZE PREVIOUS INTEGRAND 
      A(3) = Y1 
C 
C IF XDLTA IS NEGATIVE THEN 2ND ORDER INTEGRATION IS ASSUMED
      IF(XDLTA)10,20
C 
C INITIALIZATION OF A(2) AND A(4) FOR 2ND ORDER INTEGRATION 
   10 A(2) = XDLTA/12.
      A(4) = Y1 
      RETURN
C 
C INITIALIZATION OF A(2) FOR 1ST ORDER INTEGRATION
   20 A(2) = XDLTA/2. 
      RETURN
      END 
      FUNCTION AREA(DATA,A) 
C "AREA" COMPUTES THE RUNNING NUMERIC INTEGRAL OF "DATA"
C EACH TIME AREA IS CALLED IT UPDATES ARRAY "A" AND COMPUTES
C THE NEW VALUE OF THE INTEGRAL.  TO AVOID ROUNDOFF ERROR, A(1) 
C SHOULD BE SET TO ZERO EVERY 1000 CALLS WHICH RESETS THE CURRENT 
C VALUE OF THE INTEGRAL.
      DIMENSION A(4)
C 
C IF A(2) < 0 THEN PERFORM 2ND ORDER INTEGRATION
      IF(A(2))10,20 
C 
C EVALUATION OF 2ND ORDER ALGORITHM 
   10 A(1) = A(1)-(5.*DATA+8.*A(3)-A(4))*A(2) 
C 
C  UPDATE PREVIOUS VALUE OF DATA
      A(4) = A(3) 
      GO TO 30
C 
C EVALUATION OF 1ST ORDER ALGORITHM 
   20 A(1) = A(1)+(A(3)+DATA)*A(2)
C 
C  UPDATE CURRENT VALUE OF DATA 
30    A(3) = DATA 
C  RETURN CURRENT VALUE OF INTEGRAL 
      AREA = A(1) 
      RETURN
      END 
      FUNCTION FAREA(Y,N,DX1) 
C  "FAREA" COMPUTES THE NUMERICAL INTEGRAL OF THE DATA POINTS 
C   IN ARRAY "Y".  THE ABSOLUTE VALUE OF "DX" IS THE INTEGRATION
C   STEP SIZE.  IF DX IS NEGATIVE THEN 2ND ORDER INTEGRATION IS 
C   PERFORMED INSTEAD OF 1ST ORDER. 
C   THE 1ST ORDER ALGORITHM IS TRAPEZOIDAL RULE INTEGRATION.
C   THE 2ND ORDER ALGORITHM IS SIMPSON'S RULE INTEGRATION WITH
C   A CORRECTION IN THE LAST INTERVAL IF AN EVEN NUMBER OF POINTS 
C   ARE GIVEN.
      DIMENSION Y(1)
      FAREA = 0.
      NM1 = N-1 
C 
C   DETERMINE IF 1ST OR 2ND ORDER INTEGRATION 
      IF(DX1)10,20
C 
C  EVALUATION OF 2ND ORDER ALGORITHM
   10 DX = DX1/(-3.)
      IFODD = 1 
      DO 12 I = 2,NM1 
      IFODD = -IFODD
      IF(IFODD)11,12
   11 FAREA = FAREA+Y(I)
   12 FAREA = FAREA+Y(I)
C 
C  IF AN EVEN NUMBER OF POINTS COMPUTE LAST INTERVAL CORRECTION 
      IF(IFODD)22,14
   14 FAREA = FAREA+Y(NM1)/2.+(Y(N)-Y(N-2))/8.
      GO TO 22
C 
C  EVALUATION OF 1ST ORDER ALGORITHM
   20 DO 21 I = 2,NM1 
   21 FAREA = FAREA+Y(I)
      DX = DX1/2. 
   22 FAREA = (Y+Y(N)+2.*FAREA)*DX
      RETURN
      END 
      END$
                                                              