qiuni.for

计算卫星位置的fortran程序

	SUBROUTINE inv(A,N)
	implicit none
!     矩阵求逆
        integer :: N, k, i, j
        real*8  A(N,N)   !存放矩阵A
        INTEGER :: IP(N)      !记录主列号
        REAL*8 :: P, S         !工作单元，放主元
        INTEGER :: I0, R      !工作单元，放主列号
        real*8 :: EPS = 0.00000000000000001

!        write(*,*)'gjcp  ok0000000'
        DO K = 1, N
          P = 0
          I0 = K
          IP(K) = K

          DO I = K, N
            IF(ABS(A(I,K)) .GT. ABS(P))  THEN
              P = A(I,K)
              I0 = I
              IP(K) = I
            END IF
          END DO

          IF(ABS(P)+1.0.eq.1.0)  THEN
             WRITE(*,*) 'DET=0'
             stop !后来加的
             GOTO 10
          END IF

          IF(I0 .NE. K)  THEN
            DO J = 1, N
              S = A(K,J)
              A(K,J) =A (I0,J)
              A(I0,J) = S
            END DO
          END IF

          A(K, K) = 1./P

          DO I = 1, N
            IF(I .NE. K)  THEN
              A(I,K) = -A(I,K) * A(K,K)
              DO J = 1, N
                IF(J .NE. K)THEN
                  A(I,J) = A(I,J) + A(I,K) * A(K,J)
                END IF
              END DO
            END IF
          END DO

          DO J = 1, N
            IF(J .NE. K)THEN
                A(K,J) = A(K,K) * A(K,J)
            END IF
          END DO
        END DO

!        write(*,*) 'gjcp  ok01'

        DO K = N - 1, 1, -1
          R = IP(K)
          IF(R.NE.K)  THEN
            DO I = 1, N
               S = A(I,R)
               A(I,R) = A(I,K)
               A(I,K) = S
            END DO
          END IF
        END DO

10      continue
!  write(*,*)'the end'

      END