!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! This program evaluates the performance and regression rate !
! ballistic coefficients of a hybrid rocket motor based on user- !
! defined specs and raw pressure data taken within the motor and !
! upstream from the oxidizer injector. !
! !
! It implements a Runge-Kutta fourth order scheme using the !
! averaged Characteristic Velocity method to evaluate the !
! time-dependent port radius, regression rate, and oxidizer mass !
! flux. These values are used in a least squares algorith for !
! determining the ballistic coefficients that best model the fuel !
! regression behavior. !
! !
! NOTES: !
! (1) The program reads tab-separated pressure data from a !
! file named FiringData.txt located within the same !
! folder as this program. The document is headed with !
! the total number of entries. The columns are !
! time, upstream P, and Motor P respectively. !
! (2) Operation of the program requires the user to have !
! determined: !
! * Nozzle Throat Area in m^2 !
! * Total Fuel Mass cosumed during firing in Kg !
! * Total Oxidizer Mass consumed during firing in Kg !
! * Fuel Density in Kg/m^3 !
! * Fuel Port Length in m !
! * Initial Fuel Port Diameter in m !
! **These values are prompted from the command line** !
! (3) Time, Radius, Regression Rate, Oxidizer Mass Flux, !
! Chamber Pressure, and Oxidizer to Fuel Ratio !
! respectively are recorded in tab-separated columns !
! of the outputed data file FiringDataOutput.dat !
! (4) The determined Characteristic Velocity, Orifice !
! Pressure Drop Multiplier and ballistic coefficients !
! in the equation !
! dr/dt = a * Gox^n * Pc^m !
! are written to the command line. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM RocketDataAnalyzer
IMPLICIT NONE
DOUBLE PRECISION Cstar, Astar, Mox1, Mox2, Mfuel, rho, L,
. dt, C, Tmax, dnot, temp, ac, mc, nc
DOUBLE PRECISION K(5), Add(5)
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: T, P1, P2,
. mdotox, mdotf, D, Ddot, r, rdot, OF, Gox, Pdrop
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, AT,
. ATA, ATAinv, ATAinvAT, b, X
INTEGER i, n, Nmax, start, done, size
INTEGER, DIMENSION(:), ALLOCATABLE :: INDX
LOGICAL test
!
! Prompt user for motor information
!
WRITE(*, *) "Throat Area(m^2), ", "Total Fuel Mass(Kg), ",
. "Total Oxidizer Mass(Kg), ", "Fuel Density(Kg/m^3), ",
. "Fuel Length(m), ", "and Initial Diameter(m)"
READ(*, *) Astar, Mfuel, Mox1, rho, L, dnot
!
! Read pressure and time data into allocated arrays
!
OPEN (unit=99, file='FiringData.txt', status='old',
. action='read')
READ(99, *) Nmax
ALLOCATE(T(Nmax), P1(Nmax), P2(Nmax), Mdotox(Nmax), D(Nmax),
. Ddot(Nmax), r(Nmax), rdot(Nmax), OF(Nmax), Gox(Nmax),
. Pdrop(Nmax), Mdotf(Nmax))
DO n=1, Nmax
READ(99, *) T(n), P1(n), P2(n)
T(n) = T(n) - T(1)
P1(n) = P1(n) * 6894.757
P2(n) = P2(n) * 6894.757
END DO
!
! Calculate time step
!
Tmax = T(Nmax)
dt = Tmax/(Nmax-1)
start = 1
done = Nmax
!
! Interatively determine firing start and ending time indices
! based on mass conservation within the motor. In the process,
! critical operating parameters are also iteratively determined.
!
10 test = .FALSE.
Pdrop = SQRT(P1 - P2)
CALL Integrate(Pdrop, dt, start, Nmax, temp)
C = Mox1 / temp
Mdotox(:) = C * Pdrop(:)
CALL Integrate(Mdotox, dt, start, done, Mox2)
CALL Integrate(P2, dt, start, done, temp)
Cstar = Astar * temp / (Mfuel + Mox2)
DO n=start, done
Mdotf(n) = Astar * P2(n) / Cstar - Mdotox(n)
END DO
DO n=start, done
!
! Check mass conservation
!
IF (Mdotf(n) .LT. 0) THEN
Mdotf(n) = 0
Mdotox(n) = 0
IF (n .EQ. start .OR. (n .GT. 1 .AND. Mdotf(n-1)
. .LT. 0)) THEN
start = n + 1
ELSE IF (n .GT. 1 .AND. Mdotf(n-1) .GT. 0) THEN
done = n - 1
ENDIF
test = .TRUE. !Mass conservation NOT satisfied
ENDIF
END DO
if (test) THEN
GOTO 10
ENDIF
OF = Mdotox / Mdotf ! Oxidizer to Fuel Ratio array
!
! Initial Conditions
!
D(start) = dnot
r(start) = dnot / 2
Gox(start) = 4 * Mdotox(start) / (3.14159 * D(start)**2)
Add(1) = 0
Add(2) = 0
Add(3) = .5
Add(4) = .5
Add(5) = 1
!
! RK-4 algorithm for determining instantaneous port diameter
!
DO n=start, done-1
K(1) = 0
DO i=2, 5
K(i) = 2 * (Mdotf(n)) / (3.14159 * rho * L * (D(n) +
. Add(i)*K(i-1)))
END DO
D(n+1) = D(n) + dt * (K(2) + 2*K(3) + 2*K(4) + K(5)) / 6
r(n+1) = D(n+1) / 2
rdot(n) = K(2) / 2
Gox(n+1) = 4 * Mdotox(n+1) / (3.14159 * D(n+1)**2)
END DO
rdot(done) = (Mdotf(done)) / (3.14159 * rho * L * D(done))
!
! Write instantaneous operating conditions to file
!
OPEN(1, file = 'FiringDataOutput.dat')
WRITE(1, *) 'time, radius, rdot, Gox, Pc, MdotF, MdotOx OF'
DO n=start, done
WRITE(1, *) T(n), r(n), rdot(n),
. Gox(n), P2(n), MdotF(n), MdotOx(n), OF(n)
END DO
!
! Allocate matrices to be used in Least Squares Algorithm
!
size = done-start + 1
ALLOCATE(A(size,3), AT(3,size), ATA(3,3), ATAinv(3,3),
. ATAinvAT(3,size), b(size,1), X(3,1))
ALLOCATE(INDX(3))
DO n=1, size
b(n,1) = LOG(rdot(n + start -1))
A(n,1) = LOG(P2(n + start -1))
AT(1,n)= A(n,1)
A(n,2) = LOG(Gox(n + start -1))
AT(2,n)= A(n,2)
A(n,3) = 1
AT(3,n)= A(n,3)
END DO
!
! Find best x (Xhat) for approximating Ax = b using Least Squares
! Xhat = inv(A'A)A'b
!
CALL Multiply(AT, A, 3, size, 3, ATA)
CALL MIGS(ATA,3,ATAinv,INDX)
CALL Multiply(ATAinv, AT, 3, 3, size, ATAinvAT)
CALL Multiply(ATAinvAT, b, 3, size, 1, X)
mc = X(1,1)
nc = X(2,1)
ac = EXP(X(3,1))
!
! Write results to command line
!
WRITE(*, *) "Combusted Oxidizer Mass"
WRITE(*, *) Mox2
WRITE(*, *) "Average Characteristic Velocity"
WRITE(*, *) Cstar
WRITE(*, *) "Orifice Discharge Multiplier A*Cd*sqrt(rho)"
WRITE(*, *) C
WRITE(*, *) " "
WRITE(*, *) "dr/dt = (a) * Gox^(n) * Pc^(m)"
WRITE(*, *) "a"
WRITE(*, *) ac
WRITE(*, *) "m"
WRITE(*, *) mc
WRITE(*, *) "n"
WRITE(*, *) nc
END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Subroutine implementing a trapezoidal algorithm for numerically !
! integrating a discretized function f over the discrete set of !
! values differing by dx between indices Istart and Imax. !
! !
! Input Variables: !
! DOUBLE PRECISION Array f: Discrete function !
! DOUBLE PRECISION dx: The constant increment between values !
! at adjacent indices in the discrete domain of f !
! INTEGER Istart and Imax: starting and ending indices !
! Istart < Imax !
! Output Variables: !
! DOUBLE PRECISION S: The area under f A.K.A. the integral !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Integrate(f, dx, Istart, Imax, S)
IMPLICIT NONE
INTEGER Istart, Imax, i
DOUBLE PRECISION f(Imax), dx, S
S=0
DO i=Istart, Imax-1
S = S + .5 * dx * (f(i) + f(i+1))
END DO
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Subroutine used to evaluate the product of m by n matrix A and !
! n by p matix B. The result is the m by p matrix C. Values !
! contained in matrices A, B, and C are of type DOUBLE PRECISION. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Multiply(A, B, m, n, p, C)
IMPLICIT NONE
INTEGER m, n, p, i, j, k
DOUBLE PRECISION sum
DOUBLE PRECISION A(m, n), B(n, p), C(m, p), tmp(n)
DO i=1, m
DO k=1, p
sum = 0
DO j=1, n
sum = sum + A(i, j) * B(j, k)
END DO
C(i,k) = sum
END DO
END DO
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Please Note: !
! !
! (1) The following subroutines are written by Tao Pang in conjunction !
! with his book, "An Introduction to Computational Physics," !
! published by Cambridge University Press in 1997. ! !
! !
! (2) No warranties, express or implied, are made for this program. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
SUBROUTINE MIGS (A,N,X,INDX)
!
! Subroutine to invert matrix A(N,N) with the inverse stored
! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J,K
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL, INTENT (IN), DIMENSION (N,N):: A
REAL, INTENT (OUT), DIMENSION (N,N):: X
REAL, DIMENSION (N,N) :: B
!
DO I = 1, N
DO J = 1, N
B(I,J) = 0.0
END DO
END DO
DO I = 1, N
B(I,I) = 1.0
END DO
!
CALL ELGS (A,N,INDX)
!
DO I = 1, N-1
DO J = I+1, N
DO K = 1, N
B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)
. *B(INDX(I),K)
END DO
END DO
END DO
!
DO I = 1, N
X(N,I) = B(INDX(N),I)/A(INDX(N),N)
DO J = N-1, 1, -1
X(J,I) = B(INDX(J),I)
DO K = J+1, N
X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
END DO
X(J,I) = X(J,I)/A(INDX(J),J)
END DO
END DO
END SUBROUTINE MIGS
!
SUBROUTINE ELGS (A,N,INDX)
!
! Subroutine to perform the partial-pivoting Gaussian elimination.
! A(N,N) is the original matrix in the input and transformed matrix
! plus the pivoting element ratios below the diagonal in the output.
! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J,K,ITMP
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL :: C1,PI,PI1,PJ
REAL, INTENT (INOUT), DIMENSION (N,N) :: A
REAL, DIMENSION (N) :: C
!
! Initialize the index
!
DO I = 1, N
INDX(I) = I
END DO
!
! Find the rescaling factors, one from each row
!
DO I = 1, N
C1= 0.0
DO J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
END DO
C(I) = C1
END DO
!
! Search the pivoting (largest) element from each column
!
DO J = 1, N-1
PI1 = 0.0
DO I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ENDIF
END DO
!
! Interchange the rows via INDX(N) to record pivoting order
!
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
!
! Record pivoting ratios below the diagonal
!
A(INDX(I),J) = PJ
!
! Modify other elements accordingly
!
DO K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
END DO
END DO
END DO
!
END SUBROUTINE ELGS
! !
! This program evaluates the performance and regression rate !
! ballistic coefficients of a hybrid rocket motor based on user- !
! defined specs and raw pressure data taken within the motor and !
! upstream from the oxidizer injector. !
! !
! It implements a Runge-Kutta fourth order scheme using the !
! averaged Characteristic Velocity method to evaluate the !
! time-dependent port radius, regression rate, and oxidizer mass !
! flux. These values are used in a least squares algorith for !
! determining the ballistic coefficients that best model the fuel !
! regression behavior. !
! !
! NOTES: !
! (1) The program reads tab-separated pressure data from a !
! file named FiringData.txt located within the same !
! folder as this program. The document is headed with !
! the total number of entries. The columns are !
! time, upstream P, and Motor P respectively. !
! (2) Operation of the program requires the user to have !
! determined: !
! * Nozzle Throat Area in m^2 !
! * Total Fuel Mass cosumed during firing in Kg !
! * Total Oxidizer Mass consumed during firing in Kg !
! * Fuel Density in Kg/m^3 !
! * Fuel Port Length in m !
! * Initial Fuel Port Diameter in m !
! **These values are prompted from the command line** !
! (3) Time, Radius, Regression Rate, Oxidizer Mass Flux, !
! Chamber Pressure, and Oxidizer to Fuel Ratio !
! respectively are recorded in tab-separated columns !
! of the outputed data file FiringDataOutput.dat !
! (4) The determined Characteristic Velocity, Orifice !
! Pressure Drop Multiplier and ballistic coefficients !
! in the equation !
! dr/dt = a * Gox^n * Pc^m !
! are written to the command line. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM RocketDataAnalyzer
IMPLICIT NONE
DOUBLE PRECISION Cstar, Astar, Mox1, Mox2, Mfuel, rho, L,
. dt, C, Tmax, dnot, temp, ac, mc, nc
DOUBLE PRECISION K(5), Add(5)
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: T, P1, P2,
. mdotox, mdotf, D, Ddot, r, rdot, OF, Gox, Pdrop
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, AT,
. ATA, ATAinv, ATAinvAT, b, X
INTEGER i, n, Nmax, start, done, size
INTEGER, DIMENSION(:), ALLOCATABLE :: INDX
LOGICAL test
!
! Prompt user for motor information
!
WRITE(*, *) "Throat Area(m^2), ", "Total Fuel Mass(Kg), ",
. "Total Oxidizer Mass(Kg), ", "Fuel Density(Kg/m^3), ",
. "Fuel Length(m), ", "and Initial Diameter(m)"
READ(*, *) Astar, Mfuel, Mox1, rho, L, dnot
!
! Read pressure and time data into allocated arrays
!
OPEN (unit=99, file='FiringData.txt', status='old',
. action='read')
READ(99, *) Nmax
ALLOCATE(T(Nmax), P1(Nmax), P2(Nmax), Mdotox(Nmax), D(Nmax),
. Ddot(Nmax), r(Nmax), rdot(Nmax), OF(Nmax), Gox(Nmax),
. Pdrop(Nmax), Mdotf(Nmax))
DO n=1, Nmax
READ(99, *) T(n), P1(n), P2(n)
T(n) = T(n) - T(1)
P1(n) = P1(n) * 6894.757
P2(n) = P2(n) * 6894.757
END DO
!
! Calculate time step
!
Tmax = T(Nmax)
dt = Tmax/(Nmax-1)
start = 1
done = Nmax
!
! Interatively determine firing start and ending time indices
! based on mass conservation within the motor. In the process,
! critical operating parameters are also iteratively determined.
!
10 test = .FALSE.
Pdrop = SQRT(P1 - P2)
CALL Integrate(Pdrop, dt, start, Nmax, temp)
C = Mox1 / temp
Mdotox(:) = C * Pdrop(:)
CALL Integrate(Mdotox, dt, start, done, Mox2)
CALL Integrate(P2, dt, start, done, temp)
Cstar = Astar * temp / (Mfuel + Mox2)
DO n=start, done
Mdotf(n) = Astar * P2(n) / Cstar - Mdotox(n)
END DO
DO n=start, done
!
! Check mass conservation
!
IF (Mdotf(n) .LT. 0) THEN
Mdotf(n) = 0
Mdotox(n) = 0
IF (n .EQ. start .OR. (n .GT. 1 .AND. Mdotf(n-1)
. .LT. 0)) THEN
start = n + 1
ELSE IF (n .GT. 1 .AND. Mdotf(n-1) .GT. 0) THEN
done = n - 1
ENDIF
test = .TRUE. !Mass conservation NOT satisfied
ENDIF
END DO
if (test) THEN
GOTO 10
ENDIF
OF = Mdotox / Mdotf ! Oxidizer to Fuel Ratio array
!
! Initial Conditions
!
D(start) = dnot
r(start) = dnot / 2
Gox(start) = 4 * Mdotox(start) / (3.14159 * D(start)**2)
Add(1) = 0
Add(2) = 0
Add(3) = .5
Add(4) = .5
Add(5) = 1
!
! RK-4 algorithm for determining instantaneous port diameter
!
DO n=start, done-1
K(1) = 0
DO i=2, 5
K(i) = 2 * (Mdotf(n)) / (3.14159 * rho * L * (D(n) +
. Add(i)*K(i-1)))
END DO
D(n+1) = D(n) + dt * (K(2) + 2*K(3) + 2*K(4) + K(5)) / 6
r(n+1) = D(n+1) / 2
rdot(n) = K(2) / 2
Gox(n+1) = 4 * Mdotox(n+1) / (3.14159 * D(n+1)**2)
END DO
rdot(done) = (Mdotf(done)) / (3.14159 * rho * L * D(done))
!
! Write instantaneous operating conditions to file
!
OPEN(1, file = 'FiringDataOutput.dat')
WRITE(1, *) 'time, radius, rdot, Gox, Pc, MdotF, MdotOx OF'
DO n=start, done
WRITE(1, *) T(n), r(n), rdot(n),
. Gox(n), P2(n), MdotF(n), MdotOx(n), OF(n)
END DO
!
! Allocate matrices to be used in Least Squares Algorithm
!
size = done-start + 1
ALLOCATE(A(size,3), AT(3,size), ATA(3,3), ATAinv(3,3),
. ATAinvAT(3,size), b(size,1), X(3,1))
ALLOCATE(INDX(3))
DO n=1, size
b(n,1) = LOG(rdot(n + start -1))
A(n,1) = LOG(P2(n + start -1))
AT(1,n)= A(n,1)
A(n,2) = LOG(Gox(n + start -1))
AT(2,n)= A(n,2)
A(n,3) = 1
AT(3,n)= A(n,3)
END DO
!
! Find best x (Xhat) for approximating Ax = b using Least Squares
! Xhat = inv(A'A)A'b
!
CALL Multiply(AT, A, 3, size, 3, ATA)
CALL MIGS(ATA,3,ATAinv,INDX)
CALL Multiply(ATAinv, AT, 3, 3, size, ATAinvAT)
CALL Multiply(ATAinvAT, b, 3, size, 1, X)
mc = X(1,1)
nc = X(2,1)
ac = EXP(X(3,1))
!
! Write results to command line
!
WRITE(*, *) "Combusted Oxidizer Mass"
WRITE(*, *) Mox2
WRITE(*, *) "Average Characteristic Velocity"
WRITE(*, *) Cstar
WRITE(*, *) "Orifice Discharge Multiplier A*Cd*sqrt(rho)"
WRITE(*, *) C
WRITE(*, *) " "
WRITE(*, *) "dr/dt = (a) * Gox^(n) * Pc^(m)"
WRITE(*, *) "a"
WRITE(*, *) ac
WRITE(*, *) "m"
WRITE(*, *) mc
WRITE(*, *) "n"
WRITE(*, *) nc
END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Subroutine implementing a trapezoidal algorithm for numerically !
! integrating a discretized function f over the discrete set of !
! values differing by dx between indices Istart and Imax. !
! !
! Input Variables: !
! DOUBLE PRECISION Array f: Discrete function !
! DOUBLE PRECISION dx: The constant increment between values !
! at adjacent indices in the discrete domain of f !
! INTEGER Istart and Imax: starting and ending indices !
! Istart < Imax !
! Output Variables: !
! DOUBLE PRECISION S: The area under f A.K.A. the integral !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Integrate(f, dx, Istart, Imax, S)
IMPLICIT NONE
INTEGER Istart, Imax, i
DOUBLE PRECISION f(Imax), dx, S
S=0
DO i=Istart, Imax-1
S = S + .5 * dx * (f(i) + f(i+1))
END DO
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Subroutine used to evaluate the product of m by n matrix A and !
! n by p matix B. The result is the m by p matrix C. Values !
! contained in matrices A, B, and C are of type DOUBLE PRECISION. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Multiply(A, B, m, n, p, C)
IMPLICIT NONE
INTEGER m, n, p, i, j, k
DOUBLE PRECISION sum
DOUBLE PRECISION A(m, n), B(n, p), C(m, p), tmp(n)
DO i=1, m
DO k=1, p
sum = 0
DO j=1, n
sum = sum + A(i, j) * B(j, k)
END DO
C(i,k) = sum
END DO
END DO
END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Please Note: !
! !
! (1) The following subroutines are written by Tao Pang in conjunction !
! with his book, "An Introduction to Computational Physics," !
! published by Cambridge University Press in 1997. ! !
! !
! (2) No warranties, express or implied, are made for this program. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
SUBROUTINE MIGS (A,N,X,INDX)
!
! Subroutine to invert matrix A(N,N) with the inverse stored
! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J,K
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL, INTENT (IN), DIMENSION (N,N):: A
REAL, INTENT (OUT), DIMENSION (N,N):: X
REAL, DIMENSION (N,N) :: B
!
DO I = 1, N
DO J = 1, N
B(I,J) = 0.0
END DO
END DO
DO I = 1, N
B(I,I) = 1.0
END DO
!
CALL ELGS (A,N,INDX)
!
DO I = 1, N-1
DO J = I+1, N
DO K = 1, N
B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)
. *B(INDX(I),K)
END DO
END DO
END DO
!
DO I = 1, N
X(N,I) = B(INDX(N),I)/A(INDX(N),N)
DO J = N-1, 1, -1
X(J,I) = B(INDX(J),I)
DO K = J+1, N
X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
END DO
X(J,I) = X(J,I)/A(INDX(J),J)
END DO
END DO
END SUBROUTINE MIGS
!
SUBROUTINE ELGS (A,N,INDX)
!
! Subroutine to perform the partial-pivoting Gaussian elimination.
! A(N,N) is the original matrix in the input and transformed matrix
! plus the pivoting element ratios below the diagonal in the output.
! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J,K,ITMP
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL :: C1,PI,PI1,PJ
REAL, INTENT (INOUT), DIMENSION (N,N) :: A
REAL, DIMENSION (N) :: C
!
! Initialize the index
!
DO I = 1, N
INDX(I) = I
END DO
!
! Find the rescaling factors, one from each row
!
DO I = 1, N
C1= 0.0
DO J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
END DO
C(I) = C1
END DO
!
! Search the pivoting (largest) element from each column
!
DO J = 1, N-1
PI1 = 0.0
DO I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ENDIF
END DO
!
! Interchange the rows via INDX(N) to record pivoting order
!
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
!
! Record pivoting ratios below the diagonal
!
A(INDX(I),J) = PJ
!
! Modify other elements accordingly
!
DO K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
END DO
END DO
END DO
!
END SUBROUTINE ELGS
No comments:
Post a Comment