      PROGRAM MAGCAL1
C                     24 JANUARY 1993
C     FOR FURTHER INFORMATION, CONTACT:
C         A. L. Vampola
C         P. O. Box 10225
C         Torrance, CA 90505
C         (310) 375-1012
C         vampola@envnet.gsfc.nasa.gov
C
C     A MAGNETIC FIELD FILE OF THE .FPM TYPE IS GENERATED AND SAVED.
      INTEGER*2 IPRC,IA1,IA2,IA3,I1,I2,IFLAG
      INTEGER*4 KTIM1,KTIM2
      CHARACTER*50 MAGFIL,FPMFIL
      CHARACTER*1 DIGITS(10),JNC(3),ANS
      CHARACTER*3 MAG,FPM,FIX
      DATA FPM,MAG,FIX/'FPM','MAG','FIX'/
      DATA DIGITS/'0','1','2','3','4','5','6','7','8','9'/
C     MAG FILE IS UNIT 4
C     FPM FILE IS UNIT 7
      WRITE(*,1000)
1000  FORMAT(1X,'CRRES FPM Stripout Program VER. 1.1   1-24-93')
6790  PRINT*,' INPUT FILE NAME?  '
      READ(*,1004)MAGFIL
1004  FORMAT(A)
      FPMFIL=MAGFIL
      PRINT*,'DOES THE FILE NEED TO BE INVERTED? (Y,N)   '
      READ(*,1004)ANS
      CALL SUFIX(MAGFIL,MAG)
      CALL SUFIX(FPMFIL,FPM)
      PRINT*,' NUMBER OF FIRST ORBIT?    '
      READ(*,4900)I1
      PRINT*,' NUMBER OF ORBITS?  (I3)  '
      READ(*,4900)I2
4900  FORMAT(I4)
      CALL TIMER(KTIM1)                         !FOR TIMING PURPOSES
      I2=I1+I2-1
      DO 5432 IPRC=I1,I2
      IA1=IPRC/100
      IA2=(IPRC-IA1*100)/10
      IA3=I2MOD(IPRC,10)
      IF(IA1.EQ.10)IA1=0
      JNC(1)=DIGITS(IA1+1)
      JNC(2)=DIGITS(IA2+1)
      JNC(3)=DIGITS(IA3+1)
      CALL PRFIX(MAGFIL,JNC)
      CALL PRFIX(FPMFIL,JNC)
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')THEN
        CALL SUFIX(MAGFIL,MAG)
        CALL FIXMAGFL(MAGFIL)  !THIS CALL PRODUCES AN INVERTED .FIX FILE
        CALL SUFIX(MAGFIL,FIX)
      ENDIF
      CALL MAGCAL(MAGFIL,FPMFIL,IPRC,IFLAG)
      IF(IFLAG.EQ.1)GO TO 9338
5432  CONTINUE
      CALL TIMER(KTIM2)
      TME=.01*FLOAT(KTIM2-KTIM1)
      WRITE(*,4828)TME
4828  FORMAT(' ELAPSED TIME = ',F8.2)
      STOP
9338  WRITE(*,9339)
9339  FORMAT(' ERROR IN CALIBRATING MAG DATA, PROGRAM HALTED')
      STOP
      END

      SUBROUTINE MAGCAL(MAGFIL,FPMFIL,IPROC,IFLAG)
C     THIS PROGRAM READS THE RAW CRRES MAGNETOMETER DATA FILE ON THE AGENCY
C     TAPE, CONVERTS IT TO FIELD COMPONENTS, AND WRITES IT OUT.  THE OUTPUT
C     FORMAT IS UT, BX, BY, BZ IN INTEL FLOATING POINT.  UT IS IN SECONDS,
C     BX, BY, AND BZ ARE IN NANOTESLA
C     THIS SOFTWARE WAS ORIGINALLY OBTAINED FROM H. SINGER.  IT HAS BEEN
C     MODIFIED BY A. L. VAMPOLA TO RUN ON A PC AND PRODUCE AN OUTPUT FILE
C     FOR USE WITH THE REST OF HIS CRRES SOFTWARE
      INTEGER*4 ITIME(160),LEN,IORB
      INTEGER*2 NDX,NRY(148),IPROC,MRY(148,160)
      CHARACTER*2 CALFIL(10),DIGITS(21)
      CHARACTER*50 MAGFIL,FPMFIL
      CHARACTER*42 DUMY1
      CHARACTER*20 DUMY2
      EQUIVALENCE(DIGITS,DUMY1),(CALFIL,DUMY2)
      DATA DUMY2/'D:\MAGCAL\MAG000.CAL'/
      DATA DUMY1/'010203040506070809101112131415161718192021'/
C   The following is a driver routine that reads the CMF magnetic field
C     file and produces calibrated signals.  The file MAG0nn.CAL is the
C     calibration file.  Several print statements in this code are included
C     only for diagnostic purposes.  They are indicated by *****
      CALL INITMAG
      NDX=160
      IUNIT=10
      IORB=IPROC
      INDX=IORB/50
      IF(INDX.EQ.0)INDX=1
      CALFIL(8)=DIGITS(INDX)
      WRITE(*,8729)DUMY2
8729  FORMAT(' CAL FILE WILL BE: ',A)
      OPEN(UNIT=10,FILE=DUMY2,FORM='FORMATTED',STATUS='OLD')
      WRITE(*,8728)CHARNB(DUMY2)
8728  FORMAT(' OPENED ',A)
      OPEN(UNIT=4,FILE=MAGFIL,FORM='UNFORMATTED',STATUS='OLD',
     1 ACCESS='TRANSPARENT')
      INQUIRE(FILE=MAGFIL,FLEN=LEN)
      LEN=LEN/300
      WRITE(*,8728)CHARNB(MAGFIL)
      OPEN(UNIT=7,FILE=FPMFIL,STATUS='NEW',FORM='UNFORMATTED',
     1 ACCESS='TRANSPARENT')
      WRITE(*,8728)CHARNB(FPMFIL)
C      IOPT      - THE OPTION FOR CONVERSION OF THE MAGNETIC FIELD SIGNALS
C                    0=CALIBRATE BUT LEAVE IN SENSOR COORDINATES
C                    1=DUE TO THE DEPLOYMENT OF THE MAGNETOMETER BOOM TO A
C                      POSITION SIGNIFICANTLY DIFFERENT FROM GROUND MEASURED
C                      ALIGNMENT, THIS OPTION SHOULD NO LONGER BE USED
C                    2=CALIBRATE AND CONVERT TO S/C FRAME USING ON-ORBIT
C                      CALCULATED ALIGNMENT FACTORS
      IOPT=0
      IAMP=0
      IREC=0
1     IF(LEN.LT.160)NDX=LEN
      LEN=LEN-NDX
      READ(4,ERR=4000)(ITIME(L),(MRY(M,L),M=1,148),L=1,NDX)
      IREC=IREC+1
      DELT=1.024/8.
      DO 30 IMF=1,NDX
      UTSEC=ITIME(IMF)/1000.
      ISHIFT=0
      DO 31 LREC=1,4
      DO 32 L=1,30
   32 NRY(L)=MRY(L+ISHIFT,IMF)
      ISHIFT=ISHIFT+32
      K=4
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+6
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      K=K+3
      UTSEC=UTSEC+DELT
      CALL CALINT(IUNIT,IORB,UTSEC,NRY(K),NRY(K+1),NRY(K+2),IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
      UTSEC=UTSEC+DELT
   31 CONTINUE
   30 CONTINUE
C     OUTPUT TO SCREEN TO SHOW ACTIVITY AND AS A QUALITY CHECK
      SUMB2=BX*BX+BY*BY+BZ*BZ
      BTOTLY=SQRT(SUMB2)
      PRINT*,' UT AND B ARE: ',UTSEC,BTOTLY
C     IF THE .MAG FILE INCLUDES EOR's, COMMENT THIS NEXT CARD
C     THE FILES THAT COME OFF OF THE 8mm TAPES DO NOT HAVE EMBED EOR's
      IF(LEN.GT.0)GO TO 1
4000  CALL ENDCALS
      CLOSE(UNIT=10)
      CLOSE(UNIT=4)
      RETURN
      END

      SUBROUTINE CALINT(IUNIT,IORB,UTSEC,ISX,ISY,ISZ,IAMP,IOPT,BX,BY,BZ,
     1                  IFLAG)
      INTEGER*2 ISX,ISY,ISZ,IMASK,IONE,ITMP,IRY
      COMMON/SCRAT/ARRAY(15000)
      DATA IMASK/4095/,IONE/1/,IRY/1/
C******************************************************************************
C   THIS IS A FRONT-END ROUTINE THAT UNPACKS THE 12-BIT TWOS-COMPLIMENT
C     MAGNETOMETER WORDS THEN PERFORMS CALIBRATION WITH THE 
C     CALMAG.  THE USER SHOULD TAKE CARE THAT THE MAGNETOMETER WORDS
C     ISX,ISY,ISZ ARE DECLARED AS INTEGER*2 IN THE CALLING PROGRAM
C     INPUT PARAMETERS:
C       ISX,ISY,ISZ - INTEGER*2 MAGNETOMETER WORDS IN THE FOLLOWING FORMAT
C    VALUE:   0   0   (6X)   (gain)  (-- telemetered 12-Bit count --)
C    BIT:    16  15    14      13    12 11 10  9  8  7  6  5  4  3  2  1
C     SEE THE ROUTINE CALMAG FOR DEFINITIONS OF OTHER PARAMETERS
C*******************************************************************************
C   EXTRACT MODE, GAIN AND CONVERT 2-S COMPLIMENT
      NX=IAND(ISX,IMASK)
      IF(NX.GT.2047)NX=NX-4096
      NY=IAND(ISY,IMASK)
      IF(NY.GT.2047)NY=NY-4096
      NZ=IAND(ISZ,IMASK)
      IF(NZ.GT.2047)NZ=NZ-4096
      ITMP=ISX/4096
      IMX=IAND(ITMP,IONE)
      ITMP=ISY/4096
      IMY=IAND(ITMP,IONE)
      ITMP=ISZ/4096
      IMZ=IAND(ITMP,IONE)
      ITMP=ISX/8192
      IX6=IAND(ITMP,IONE)
      CALL CALMAG(IUNIT,IORB,UTSEC,NX,IMX,NY,IMY,IX6,NZ,IMZ,IAMP,IOPT,
     1            BX,BY,BZ,IFLAG)
C***************************************************************************
      ARRAY(IRY)=UTSEC
      ARRAY(IRY+1)=BX
      ARRAY(IRY+2)=BY
      ARRAY(IRY+3)=BZ
      IRY=IRY+4
      IF(IRY.LT.15000)RETURN
      WRITE(7)ARRAY
      IRY=1
      RETURN
      ENTRY ENDCALS
      IF(IRY.EQ.1)GO TO 767
      IRY=IRY-1
      WRITE(UNIT=7)(ARRAY(I),I=1,IRY)
      IRY=1
767   CLOSE(UNIT=7)
      RETURN
      END

      SUBROUTINE CALMAG(IUNIT,IORB,UTSEC,NX,IMX,NY,IMY,IX6,NZ,IMZ,
     1                  IAMP,IOPT,BX,BY,BZ,IFLAG)
C*******************************************************************************
C   SUBROUTINE TO CALIBRATE CRRES FLUXGATE MAGNETOMETER DATA
C      VERSION 0 CREATED 16 JUNE 1990
C      VERSION 1 CREATED 19 NOVEMBER 1990 FOR ORIGINAL RELEASE
C   INPUT PARAMTERS:
C      IUNIT     - THE UNIT NUMBER ON WHICH THE MAGNETOMETER CALIBRATION FILE
C                  HAS BEEN OPENED
C      IORB      - THE NUMBER OF THE ORBIT
C      UTSEC     - THE UT OF THE DATA POINT IN SECONDS
C      NX,NY,NZ  - THE SENSOR SIGNALS IN COUNTS
C      IMX,IMY,
C       IMZ      - THE MODE INDICATORS FOR THE THREE SENSORS (0=LOW 1=HIGH)
C      IX6       - THE X6 MULTIPLICATION MODE FOR THE Y-SENSOR (0=X1 1=X6)
C      IAMP      - THE STATE OF THE POWER AMPLIFIERS (B2, HIGHEST 2 BITS)
C                    0=NEITHER PA1 NOR PA2 IS ON 1=EITHER PA1 OR PA2 IS ON
C      IOPT      - THE OPTION FOR CONVERSION OF THE MAGNETIC FIELD SIGNALS
C                    0=CALIBRATE BUT LEAVE IN SENSOR COORDINATES
C                    1=DUE TO THE DEPLOYMENT OF THE MAGNETOMETER BOOM TO A
C                      POSITION SIGNIFICANTLY DIFFERENT FROM GROUND MEASURED
C                      ALIGNMENT, THIS OPTION SHOULD NO LONGER BE USED
C                    2=CALIBRATE AND CONVERT TO S/C FRAME USING ON-ORBIT
C                      CALCULATED ALIGNMENT FACTORS
C   OUTPUT PARAMETERS:
C      BX,BY,BZ   - THE MAGNETIC FIELD COMPONENTS IN NANO-TESLA
C      IFLAG      - THE STATUS FLAG FOR THE CONVERSION
C                    0=CALIBRATION FACTORS FOUND FOR SPECIFIED TIME PERIOD
C                   -1=CALIBRATION FACTORS NOT FOUND, CONVERSION DONE WITH
C                      FACTORS FROM ADJACENT ORBIT
C                   -2=REQUESTED TIME AFTER LAST ENTRY IN CALIBRATION FILE,
C                      CONVERSION DONE WITH LAST AVAILABLE SET OF CONSTANTS
C                      BUT CHECK TO SEE IF LATEST VERSION OF CALIBRATION FILE
C                      IS BEING USED
C*******************************************************************************
      DIMENSION A(3,3),GHIGH(4),GLOW(4),DHIGH(4),DLOW(4),PAHIGH(4),
     1          PALOW(4)
      CHARACTER*80 LINE
C   THESE VARIABLES ARE DECLARED AS DATA IN ORDER TO ENSURE THAT THEY
C     RETAIN THEIR VALUES BETWEEN CALLS
      DATA A/9*0./,GHIGH/4*0./,GLOW/4*0./DHIGH/4*0./,DLOW/4*0./
     1     JORB/-1/,JFLAG/-1/,JOPT/-1/
C   THE NOMINAL VALUES OF THE ALIGNMENTS ARE HARD WIRED HERE
      DATA THX,THY,THZ,PHX,PHY,PHZ/-1.,-1.,-1.,-1.,-1.,-1./
C     DON'T GO THROUGH CAL ROUTINE MORE THAN ONCE PER ORBIT
      DATA IFRST/0/
      IF(IFRST.EQ.1)GO TO 20
      IFRST=1
C   CHECK AND ABORT IF IOPT=1
      IF(IOPT.EQ.1)THEN
        PRINT*,'PLEASE USE IOPT=0 FOR SENSOR SIGNALS OR'
        PRINT*,'           IOPT=2 FOR SPACECRAFT FRAME SIGNALS'
        STOP
      ENDIF
C   SEE IF THE CORRECT ORBIT IS ALREADY LOADED
      IF(JORB.NE.IORB)THEN
        REWIND(IUNIT)
        READ(IUNIT,100)LINE
        READ(IUNIT,100)LINE
        READ(IUNIT,100)LINE
        READ(IUNIT,100)LINE
  100   FORMAT(A80)
        NEAR=99999
        IFOUND=0
   10   READ(IUNIT,*,END=11)KORB,KYR,KDAY,KUT1,KUT2,GXL,GYL,GZL,GXH,
     1      GYH,GZH,DXL,DYL,DZL,DXH,DYH,DZH,GYL6,GYH6,DYL6,DYH6,PAXL,
     2      PAYL,PAZL,PAXH,PAYH,PAZH,PA6L,PA6H,TTX,TTY,TTZ,PPX,PPY,PPZ,
     3      ITIME,IDATE,IQUAL

        IF(KORB.EQ.IORB.OR.ABS(KORB-IORB).LT.NEAR)THEN
          NEAR=ABS(KORB-IORB)
          IF(KORB.EQ.IORB)THEN
            JFLAG=0
            JORB=KORB
          ELSE
            JFLAG=-1
            JORB=IORB
          ENDIF
          GHIGH(1)=GXH
          GHIGH(2)=GYH
          GHIGH(3)=GYH6
          GHIGH(4)=GZH
          DHIGH(1)=DXH
          DHIGH(2)=DYH
          DHIGH(3)=DYH6
          DHIGH(4)=DZH
          GLOW(1)=GXL
          GLOW(2)=GYL
          GLOW(3)=GYL6
          GLOW(4)=GZL
          DLOW(1)=DXL
          DLOW(2)=DYL
          DLOW(3)=DYL6
          DLOW(4)=DZL
          PAHIGH(1)=PAXH
          PAHIGH(2)=PAYH
          PAHIGH(3)=PA6H
          PAHIGH(4)=PAZH
          PALOW(1)=PAXL
          PALOW(2)=PAYL
          PALOW(3)=PA6L
          PALOW(4)=PAZL
          THX=TTX
          THY=TTY
          THZ=TTZ
          PHX=PPX
          PHY=PPY
          PHZ=PPZ
        ENDIF
        GOTO10
C   IF WE'VE FOUND ONE ALREADY, O.K. ELSE LABEL JFLAG AS -2 AND
C      TAKE THE LAST SET IN THE FILE
   11   IF(IORB.LE.KORB)GOTO13
        JORB=IORB
        JFLAG=-2
        GHIGH(1)=GXH
        GHIGH(2)=GYH
        GHIGH(3)=GYH6
        GHIGH(4)=GZH
        DHIGH(1)=DXH
        DHIGH(2)=DYH
        DHIGH(3)=DYH6
        DHIGH(4)=DZH
        GLOW(1)=GXL
        GLOW(2)=GYL
        GLOW(3)=GYL6
        GLOW(4)=GZL
        DLOW(1)=DXL
        DLOW(2)=DYL
        DLOW(3)=DYL6
        DLOW(4)=DZL
        PAHIGH(1)=PAXH
        PAHIGH(2)=PAYH
        PAHIGH(3)=PA6H
        PAHIGH(4)=PAZH
        PALOW(1)=PAXL
        PALOW(2)=PAYL
        PALOW(3)=PA6L
        PALOW(4)=PAZL
        THX=TTX
        THY=TTY
        THZ=TTZ
        PHX=PPX
        PHY=PPY
        PHZ=PPZ
   13   CONTINUE
C     RESET JFLAG TO RECALCULATE ALIGNMENT MATRIX
C     (SEE BELOW)
        JOPT=-1
        IFLAG=JFLAG
C*****************************************************************************
        PRINT*,'IOPT IS ',IOPT
        PRINT*,'HIGH GAIN GAINS  ',GHIGH
        PRINT*,'LOW GAIN GAINS   ',GLOW
        PRINT*,'HIGH GAIN OFFSTS ',DHIGH
        PRINT*,'LOW GAIN OFFSTS  ',DLOW
        PRINT*,'HIGH GAIN P.A.   ',PAHIGH
        PRINT*,'LOW GAIN P.A.    ',PALOW
        PRINT*,'ALIGNMENTS       ',THX,THY,THZ,PHX,PHY,PHZ
        PRINT*,'        FLAG.....',IFLAG
C*****************************************************************************
      ENDIF
C   SEE IF IT'S NECESSARY TO RECALCULATE THE ALIGNMENTS MATRIX
C     (A) ON THE FIRST CALL
C     (B) IF JOPT HAS CHANGED
C     (C) IF A NEW SET OF COEFFICIENTS HAS BEEN LOADED
C   BUT ONLY IF IOPT .NE. 0
      IF(IOPT.EQ.0)GOTO20
      IF(IOPT.NE.JOPT)THEN
        JOPT=IOPT
         A(1,1)=COS(THX)*COS(PHX)
         A(1,2)=COS(THX)*SIN(PHX)
         A(1,3)=SIN(THX)
         A(2,1)=COS(THY)*COS(PHY)
         A(2,2)=COS(THY)*SIN(PHY)
         A(2,3)=SIN(THY)
         A(3,1)=COS(THZ)*COS(PHZ)
         A(3,2)=COS(THZ)*SIN(PHZ)
         A(3,3)=SIN(THZ)
C   INVERT THE MATRIX
        CALL INVDET(A,3)
C*****************************************************************************
        PRINT300,((A(JJ,KK),KK=1,3),JJ=1,3)
  300   FORMAT(1X,'S/C CONVERSION MATRIX:',//,10X,3F12.5,/,10X,3F12.5,/,
     1      10X,3F12.5)
C*****************************************************************************
      ENDIF
   20 CONTINUE
C   CONVERT COUNTS TO nT
      IF(IMX.EQ.0)THEN
        BBX=NX-DLOW(1)
        IF(IAMP.EQ.1)BBX=BBX-PALOW(1)
        BBX=BBX*GLOW(1)
      ELSE
        BBX=NX-DHIGH(1)
        IF(IAMP.EQ.1)BBX=BBX-PAHIGH(1)
        BBX=BBX*GHIGH(1)
      ENDIF
      IF(IX6.EQ.0)THEN
        IF(IMY.EQ.0)THEN
          BBY=NY-DLOW(2)
          IF(IAMP.EQ.1)BBY=BBY-PALOW(2)
          BBY=BBY*GLOW(2)
        ELSE
          BBY=NY-DHIGH(2)
          IF(IAMP.EQ.1)BBY=BBY-PAHIGH(2)
          BBY=BBY*GHIGH(2)
        ENDIF
      ELSE
        IF(IMY.EQ.0)THEN
          BBY=NY-DLOW(3)
          IF(IAMP.EQ.1)BBY=BBY-PALOW(3)
          BBY=BBY*GLOW(3)
        ELSE
          BBY=NY-DHIGH(3)
          IF(IAMP.EQ.1)BBY=BBY-PAHIGH(3)
          BBY=BBY*GHIGH(3)
        ENDIF
      ENDIF
      IF(IMZ.EQ.0)THEN
        BBZ=NZ-DLOW(4)
        IF(IAMP.EQ.1)BBZ=BBZ-PALOW(4)
        BBZ=BBZ*GLOW(4)
      ELSE
        BBZ=NZ-DHIGH(4)
        IF(IAMP.EQ.1)BBZ=BBZ-PAHIGH(4)
        BBZ=BBZ*GHIGH(4)
      ENDIF
      IF(IOPT.EQ.0)THEN
        BX=BBX
        BY=BBY
        BZ=BBZ
      ELSE
        BX=A(1,1)*BBX+A(1,2)*BBY+A(1,3)*BBZ
        BY=A(2,1)*BBX+A(2,2)*BBY+A(2,3)*BBZ
        BZ=A(3,1)*BBX+A(3,2)*BBY+A(3,3)*BBZ
      ENDIF
      RETURN
      ENTRY INITMAG
      IFRST=0
      RETURN
      END


      SUBROUTINE INVDET(C,N)
C   GENERAL MATRIX INVERSION ROUTINE
C   DIMENSIONS OF C MUST MATCH THAT IN CALLING
C     ROUTINE
C   INVERSE IS RETURNED IN C
      DIMENSION C(3,3),J(200)
      PD=1.
      DO 124 L=1,N
      DD=0.
      DO 123 K=1,N
  123 DD=DD+C(L,K)*C(L,K)
      DD=SQRT(DD)
  124 PD=PD*DD
      DETM=1.
      DO 125 L=1,N
  125 J(L+20)=L
      DO 144 L=1,N
      CC=0.
      M=L
      DO 135 K=L,N
      IF((ABS(CC)-ABS(C(L,K))).GE.0.)GOTO135
  126 M=K
      CC=C(L,K)
  135 CONTINUE
  127 IF(L.EQ.M)GOTO138
  128 K=J(M+20)
      J(M+20)=J(L+20)
      J(L+20)=K
      DO 137 K=1,N
      S=C(K,L)
      C(K,L)=C(K,M)
  137 C(K,M)=S
  138 C(L,L)=1.
      DETM=DETM*CC
      DO 139 M=1,N
  139 C(L,M)=C(L,M)/CC
      DO 142 M=1,N
      IF(L.EQ.M)GOTO142
  129 CC=C(M,L)
      IF(CC.EQ.0.)GOTO142
  130 C(M,L)=0.
      DO 141 K=1,N
  141 C(M,K)=C(M,K)-CC*C(L,K)
  142 CONTINUE
  144 CONTINUE
      DO 143 L=1,N
      IF(J(L+20).EQ.L)GOTO143
  131 M=L
  132 M=M+1
      IF(J(M+20).EQ.L)GOTO133
  136 IF(N.GT.M)GOTO132
  133 J(M+20)=J(L+20)
      DO 163 K=1,N
      CC=C(L,K)
      C(L,K)=C(M,K)
  163 C(M,K)=CC
      J(L+20)=L
  143 CONTINUE
      DETM=ABS(DETM)
C     THE FOLLOWING STATEMENT IS IN THE ORIGINAL SINGER CODE BUT DOESN'T
C     DO ANYTHING.  IT HAS BEEN COMMENTED UNTIL WE FIGURE OUT WHAT IT IS FOR
C      DTNRM=DETM/PD
      RETURN
      END

      SUBROUTINE PRFIX(NAME,SFX)
      CHARACTER*1 NAME(51),SFX(3),SFX1
      DATA SFX1/'.'/
      DO 1 I=1,51
      IF(NAME(I).EQ.'.')GOTO4
1     IF(NAME(I).EQ.' ')GOTO3
      WRITE(*,2)NAME
2     FORMAT(' ***** ERROR ***** ERROR ***** ERROR *****',//,
     1 A,//,' exceeds the limit on file names (50 characters)')
      STOP
3     NAME(I)=SFX1
4     NAME(I-3)=SFX(1)
      NAME(I-2)=SFX(2)
      NAME(I-1)=SFX(3)
      RETURN
      END

      SUBROUTINE SUFIX(NAME,SUF)
      CHARACTER*1 NAME(51),SUF(3),SUF1
      DATA SUF1/'.'/
      DO 1 I=1,51
      IF(NAME(I).EQ.'.')GOTO4
1     IF(NAME(I).EQ.' ')GOTO3
      WRITE(*,2)NAME
2     FORMAT(' ***** ERROR ***** ERROR ***** ERROR *****',//,
     1 A,//,' exceeds the limit on file names (50 characters)')
      STOP
3     NAME(I)=SUF1
4     NAME(I+1)=SUF(1)
      NAME(I+2)=SUF(2)
      NAME(I+3)=SUF(3)
      RETURN
      END

      SUBROUTINE FIXMAGFL(FILNAM)
      INTEGER*4 NLEN,FLE$N,ISIZE,JSIZE
      INTEGER*2 IOPEN2,IOPEN3
      LOGICAL*1 IRCD(48020),DATOUT(48020),EXIS$T
C     IRCD HOLDS ONE RECORD OF CRRES MAGNETOMETER FILE DATA
      CHARACTER*30 OUTFILE,FILNAM
      CHARACTER*1 SUF(3)
      DATA SUF/'F','I','X'/
      VER=1.1
      WRITE(*,1000)VER
1000  FORMAT(1X,'CRRES MAGNETOMETER FILE INVERSION PROGRAM VER. ',F3.1)
C     OPEN THE FILES
1     IOPEN3=0
1004  FORMAT(A)
      INQUIRE(FILE=FILNAM,EXIST=EXIS$T,FLEN=FLE$N)
      IF(EXIS$T)GOTO222
      WRITE(*,221)CHARNB(FILNAM)
221   FORMAT(' FILE ',A,' DOES NOT EXIST')
      STOP
222   OPEN (UNIT=53,FILE=FILNAM,STATUS='OLD',ACCESS='TRANSPARENT',
     1 FORM='UNFORMATTED')
223   IOPEN3=1
      NLEN=FLE$N/4800
      WRITE(*,2252)CHARNB(FILNAM),NLEN
2252  FORMAT(' FILE ',A,' IS ',I8,' FRAMES LONG',/)
      OUTFILE=FILNAM
      CALL SUFIX(OUTFILE,SUF)
      WRITE(*,1010)CHARNB(OUTFILE)
1010  FORMAT(X,'THE OUTPUT FILE WILL BE NAMED ',A,/)
      INQUIRE(FILE=OUTFILE,EXIST=EXIS$T)
      IF(EXIS$T)WRITE(*,2003)CHARNB(OUTFILE)
2003  FORMAT(X,A,' ALREADY EXISTS.',/)
      IF(EXIS$T)STOP
6     OPEN (UNIT=52,FILE=OUTFILE,STATUS='NEW',ACCESS='TRANSPARENT',
     1 FORM='UNFORMATTED')
      IOPEN2=1
      ISIZE=4800
C     NOW WE PROCESS THE DATA
      JSIZE=10*ISIZE
100   IF(NLEN.LT.10)JSIZE=NLEN*ISIZE
      READ(UNIT=53,END=4000)(IRCD(JJ),JJ=1,JSIZE)
      DO 510 IKL=1,JSIZE,ISIZE
      L=IKL
      NLEN=NLEN-1
      DO 500 I=1,16                    !MASTER FRAME CYCLING
      DO 320 J=1,4                     !MINOR FRAME CYCLING
      DATOUT(L)=IRCD(L+3)
      DATOUT(L+1)=IRCD(L+2)
      DATOUT(L+2)=IRCD(L+1)
      DATOUT(L+3)=IRCD(L)
      L=L+4
      DO 300 M=1,30
      DATOUT(L)=IRCD(L+1)
      DATOUT(L+1)=IRCD(L)
300   L=L+2
320   CONTINUE
      DO 350 M=1,44
      DATOUT(L)=IRCD(L)
350   L=L+1
400   CONTINUE
500   CONTINUE
510   CONTINUE
      WRITE(*,519)NLEN
519   FORMAT('+FRAMES TO GO: ',I3)
      WRITE(UNIT=52,ERR=410)(DATOUT(JJ),JJ=1,JSIZE)
      IF(NLEN.GT.0)GO TO 100
      GO TO 70
410   WRITE(*,1014)
1014  FORMAT(' ERROR DURING WRITING DATA',/,' JOB ABORTED IN MAGFIX')
      GO TO 70
4000  PRINT*,' END OF DATA, PROCESSING STOPPED'
70    IF(IOPEN3.EQ.1)CLOSE(UNIT=53)
      IF(IOPEN2.EQ.1)CLOSE(UNIT=52)
      RETURN
      END
