CalcPolynomCoef Subroutine

public subroutine CalcPolynomCoef(OrderedPair, PolynomCoef)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=r64),			DIMENSION(2,MaxOrderedPairs)	::	OrderedPair
real(kind=r64),			DIMENSION(MaxPolynomOrder+1)	::	PolynomCoef
Called By

Help
Source Code

CalcPolynomCoef
Source Code

SUBROUTINE CalcPolynomCoef(OrderedPair,PolynomCoef)
          ! SUBROUTINE INFORMATION:
          !       AUTHOR   Unknown
          !       DATE WRITTEN   Unknown
          !       DATE REWRITTEN  April 1997 by Russell D. Taylor, Ph.D.
          !       MODIFIED
          !       RE-ENGINEERED

          ! PURPOSE OF THIS SUBROUTINE:
          ! Fits polynomial of order from 1 to MaxPolynomOrder to the
          ! ordered pairs of data points X,Y

          ! USE STATEMENTS:
          ! na

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
    REAL(r64), DIMENSION(MaxPolynomOrder+1) :: PolynomCoef
    REAL(r64), DIMENSION(2,MaxOrderedPairs) :: OrderedPair

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    LOGICAL :: Converged
    REAL(r64), DIMENSION(2,10) :: OrdPairSum
    REAL(r64), DIMENSION(10,10) :: OrdPairSumMatrix
    REAL(r64) B
    INTEGER I,II,J
    INTEGER :: PolynomOrder
    INTEGER :: CurrentOrder
    INTEGER :: CurrentOrdPair
    REAL(r64) S1,S2

    OrdPairSum(1,1) = MaxOrderedPairs
    OrdPairSum(1,2:3) = 0.0d0
    OrdPairSum(2,:) = 0.0d0
    PolynomCoef = 0.0d0
    DO CurrentOrdPair = 1, MaxOrderedPairs
      OrdPairSum(1,2) = OrdPairSum(1,2) + OrderedPair(1, CurrentOrdPair)
      OrdPairSum(1,3) = OrdPairSum(1,3) + OrderedPair(1, CurrentOrdPair) * &
                        OrderedPair(1, CurrentOrdPair)
      OrdPairSum(2,1) = OrdPairSum(2,1) + OrderedPair(2, CurrentOrdPair)
      OrdPairSum(2,2) = OrdPairSum(2,2) + OrderedPair(1, CurrentOrdPair) * &
                        OrderedPair(2, CurrentOrdPair)
    END DO
    PolynomOrder = 1
    Converged = .FALSE.
    DO WHILE (.NOT. Converged)
      DO CurrentOrder = 1, PolynomOrder + 1
        DO J = 1, PolynomOrder + 1
          OrdPairSumMatrix(CurrentOrder, J) = OrdPairSum(1,J - 1 + CurrentOrder)
        END DO !End of J loop
        OrdPairSumMatrix(CurrentOrder, (PolynomOrder + 2)) = OrdPairSum(2,CurrentOrder)
      END DO !End of CurrentOrder loop

      DO CurrentOrder = 1, PolynomOrder + 1
        OrdPairSumMatrix((PolynomOrder + 2), CurrentOrder)= -1.d0
        DO J = CurrentOrder + 1,(PolynomOrder + 2)
          OrdPairSumMatrix((PolynomOrder + 2), J) = 0.0d0
        END DO !End of J loop

        DO II = 2, (PolynomOrder + 2)
          DO J = CurrentOrder + 1, (PolynomOrder + 2)
            OrdPairSumMatrix(II, J) = OrdPairSumMatrix(II, J) - OrdPairSumMatrix(1, J) &
                                      * OrdPairSumMatrix(II, CurrentOrder) / &
                                      OrdPairSumMatrix(1, CurrentOrder)
          END DO !End of J loop
        END DO !End of II loop
        DO II = 1, PolynomOrder + 1
          DO J = CurrentOrder + 1, (PolynomOrder + 2)
            OrdPairSumMatrix(II, J)=OrdPairSumMatrix(II + 1, J)
          END DO !End of J loop
        END DO !End of II loop
      END DO !End of CurrentOrder loop

      S2 = 0.0d0
      DO CurrentOrdPair = 1, MaxOrderedPairs
        S1 = OrdPairSumMatrix(1, (PolynomOrder + 2))
        DO CurrentOrder = 1, PolynomOrder
          S1 = S1 + OrdPairSumMatrix(CurrentOrder + 1,(PolynomOrder + 2)) &
               * OrderedPair(1, CurrentOrdPair)**CurrentOrder
        END DO !End of CurrentOrder loop
        S2 = S2 + (S1 - OrderedPair(2, CurrentOrdPair)) * &
             (S1 - OrderedPair(2, CurrentOrdPair))
      END DO !End of CurrentOrdPair loop
      B = MaxOrderedPairs - (PolynomOrder + 1)
      IF(S2 .GT. 0.0001d0) S2 = SQRT(S2 / B)
      DO CurrentOrder = 1, PolynomOrder + 1
        PolynomCoef(CurrentOrder) = OrdPairSumMatrix(CurrentOrder, (PolynomOrder + 2))
      END DO !End of CurrentOrder loop

      IF (((PolynomOrder - MaxPolynomOrder) .LT. 0.0d0) .AND. &
          ((S2 - PolyConvgTol) .GT. 0.0d0)) THEN
        PolynomOrder = PolynomOrder + 1
        J = 2 * PolynomOrder
        OrdPairSum(1, J:J+1) = 0.0d0
        OrdPairSum(2, PolynomOrder + 1)=0.0d0
        DO I = 1, MaxOrderedPairs
          OrdPairSum(1, J) = OrdPairSum(1, J) + OrderedPair(1, I)**(J-1)
          OrdPairSum(1, J+1)=OrdPairSum(1, J+1) + OrderedPair(1, I)**J
          OrdPairSum(2, PolynomOrder+1) = OrdPairSum(2, PolynomOrder+1) + &
                                          OrderedPair(2, I) * OrderedPair(1, I)**PolynomOrder
        END DO
      ELSE
        Converged = .TRUE.
      END IF

    END DO

    RETURN
END SUBROUTINE CalcPolynomCoef
All Procedures
