C LAPOD
C -----
C Least-Squares Refinement of Cell Dimensions from Powder Data by
C Cohen's Method
C J. Ian Langford, School of Physics & Space Research, University of
C Birmingham, Birmingham B15 2TT, UK.
C Tel. (44) 21 414 4662/4666; FAX (44) 21 414 4719
C E-Mail:  J.I.LANGFORD@BHAM.AC.UK
C Langford, J.I., J. Appl. Cryst. 4, 259-260 (1971) & 6, 190-196 (1973)
C
C 13.ii.70: Original version
C 17.i.84:  FORT77 Version
C 12.v.94:  PC version (Ron Pflaumer)
C 12.i.98:  Minor upgrade
C 17 May 1999 - Put on CCP14 site: http://www.ccp14.ac.uk/ccp/web-mirrors/lapod-langford/
C
C Input:
C Input data via unit IR (fn.LAP)
C Line 1:     Title
C Line 2:     Wavelength (A), Diffractometer radius (mm)
C Line 3:     Crystal System Code, 1 Triclinic, 2 Mononclinic,
C             3 Orthorhombic, 4 Hexagonal, 5 Rhombohedral, 6 Tetragonal,
C             7 Cubic
C Line 4:     Trial cell dimensions, ALL SIX MUST BE INPUT, values not
C             critical  
C Line 5:     h k l      2Theta     weight
C et seq.
C              
C Output:
C Output is to Unit IW (fn.OUT)
C NB Vc = Q(Obs) - Q(Cal), where Q = 1/d*d
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*8 FLNAME
      CHARACTER*16 INAME,ONAME
      CHARACTER*80 TITLE
      DIMENSION CD(6),IH(100),IK(100),IL(100),T(100),DE(100),
     1P(10),QE(100),QP(100),PP(300,10),Z(3),X(9,9),RP(9),RD(10),RO(9),
     2RC(6),CO(6),RX(6),DR(100),VO(100),ED(100),Y(9),COF(9,9),
     3SD(10),TT(100),ET(100),WT(100),PWTD(10),TTCOR(100),
     4TTOBS(100)
C
C Input & Output Devices
      IR = 1
      IT = 5
      IW = 7
C      
      PI = 4.0*DATAN(1.0D0)
C For B'ham Picker diffractometer, R = 457.2/pi mm
C     RADIUS = 457.2/pi
C Set Instrumental Corrections & Correction Code to Zero
      ZEROC = 0.0
      SSDC = 0.0
      EMUC = 0.0
      L = 0
C Some program information & references
      PRINT *,' Least-Squares Refinement of Cell Dimensions from Powder
     +Data'
      PRINT *,' by the Cohen Method'
      PRINT *,' J. Ian Langford, School of Physics & Space Research, Uni
     +versity of Birmingham'
      PRINT *,' Birmingham B15 2TT, UK.'
      PRINT *,' Tel. (44) 21 414 4662/4666; FAX (44) 21 414 4719'
C
C Open files for input and output
      PRINT *,'Input Data Filename (Extension must be .LAP)'
      PRINT *,'Default: LAPOD.LAP'
      READ (IT,'(A)') FLNAME
      IF (FLNAME.EQ.' ') FLNAME = 'LAPOD'
      INAME = FLNAME//'.LAP'
      OPEN (UNIT=IR,FILE=INAME,STATUS='OLD',ERR=10000,FORM='FORMATTED',
     +ACCESS='SEQUENTIAL')
      READ (IR,'(A)',ERR=10000) TITLE
      REWIND IR
 1900 PRINT *,'Output Data Filename (Extension will be .OUT)'
      PRINT *,'Existing file will be overwritten!'
      PRINT *,'Default: LAPOD.OUT'
      READ (IT,'(A)') FLNAME
      IF (FLNAME.EQ.' ') FLNAME = 'LAPOD'
      ONAME = FLNAME//'.OUT'
      OPEN (UNIT=IW,FILE=ONAME,STATUS='UNKNOWN',ERR=10000,
     +FORM='FORMATTED',ACCESS='SEQUENTIAL')
      PRINT *,' '
      PRINT *,' '
 2010 PRINT *,'Opening file ',ONAME,' for output'
C
C Read Correction Code
   90 PRINT *,' Input Instrumental Correction Codes'
      PRINT *,'   1: Correct TwoTheta Values for known Zero Error, etc'
      PRINT *,'   2: No Instrumental Correction'
      PRINT *,'   3: Include Zero Error'
      PRINT *,'   4: Include SSD'
      PRINT *,'   5: Include Zero Error & SSD'
      PRINT *,' '
      READ  (IT,*) ICODE
      IF (ICODE.GT.5) GO TO 90
      L = ICODE - 2
      IF (L.EQ.2) L = 4
      IF (L.EQ.3) L = 2
      IF (ICODE.EQ.1) THEN
         PRINT *,' Correction data for obs. TwoTheta values:'
         PRINT *,' Input Zero Error, SSD (mm) & mu (1/mm)'
         READ (IT,*) ZEROC,SSDC,EMUC
      END IF
C
    1 READ  (IR,'(A)') TITLE
      READ  (IR,*) ALG,RADIUS
      READ  (IR,*) M
      READ  (IR,*) (CD(I), I = 1,6)
      WRITE (IW,50) TITLE
      WRITE (IW,51)
      WRITE (IW,52)
      WRITE (IW,53)
      WRITE (IW,'(A)') ' Refinement Code:'
      IF (L.EQ.-1) WRITE (IW,'(2X,I2,A)') ICODE,': 2Theta corrected for
     &pre-determined zero error, etc'
      IF (L.EQ.0) WRITE (IW,'(2X,I2,A)') ICODE,': No correction to 2Thet
     &a'
      IF (L.EQ.1) WRITE (IW,'(2X,I2,A)') ICODE,': Zero error refined'
      IF (L.EQ.2) WRITE (IW,'(2X,I2,A)') ICODE,': Zero + SSD refined'
      IF (L.EQ.4) WRITE (IW,'(2X,I2,A)') ICODE,': SSD refined'
      IF (M.EQ.1) WRITE (IW,'(/A)') ' Triclinic (1)'
      IF (M.EQ.2) WRITE (IW,'(/A)') ' Monoclinic (2)'
      IF (M.EQ.3) WRITE (IW,'(/A)') ' Orthorhombic (3)'
      IF (M.EQ.4) WRITE (IW,'(/A)') ' Hexagonal (4)'
      IF (M.EQ.5) WRITE (IW,'(/A)') ' Rhombohedral (5)'
      IF (M.EQ.6) WRITE (IW,'(/A)') ' Tetragonal (6)'
      IF (M.EQ.7) WRITE (IW,'(/A)') ' Cubic (7)'
      WRITE (IW,54) ALG
      IF (M.LT.8) THEN
         WRITE (IW,62)
         WRITE (IW,57)
      ELSE
         WRITE (IW,60)
         WRITE (IW,55)
      END IF
      WRITE (IW,63) (CD(I), I = 1,3)
      WRITE (IW,64) (CD(I), I = 4,6)
      I = 1
    6 READ (IR,*,END=120) IH(I),IK(I),IL(I),TTOBS(I),WT(I)
      TTCOR(I) = TTOBS(I)
      IF (L.EQ.-1) THEN
         TOBSR = PI*TTOBS(I)/360.0
         A = ZEROC
         B = -360.0*SSDC/(PI*RADIUS)
         C = 0.0
         IF (EMUC.GT.0.0) C = 180.0/(2.0*PI*EMUC*RADIUS)
         TTCOR(I) = TTCOR(I) + A + B*DCOS(TOBSR) + C*DSIN(2.0*TOBSR)
      END IF
      T(I) = TTCOR(I)/2.0
      I = I + 1
      GO TO 6
  120 REWIND IR
      N = I - 1
    8 IF (M.EQ.1) KK = 6
      IF (M.EQ.2) KK = 4
      IF (M.EQ.3) KK = 3
      IF (M.GE.4.AND.M.NE.7) KK = 2
      IF (M.EQ.7) KK = 1
      IF (M.EQ.8) KK = 0
      IF (L.LT.4) THEN
         MM = 6 + L
      ELSE IF (L.EQ.4) THEN
         MM = 7
      END IF
      NM = N - KK - L
C
C FORM RECIPROCAL CELL PARAMETERS
C
      DO 10 I = 4,6
      CD(I) = CD(I)*PI/180.0
   10 CONTINUE
      CALL RECIP(CD,RC,V,RV,W)
      J = 2
      K = 3
      DO 11 I = 1,3
      RP(I) = RC(I)**2
      RP(I+3) = 2.0*RC(J)*RC(K)*RC(I+3)
      CALL TEST3(J)
      CALL TEST3(K)
   11 CONTINUE
      RP(7) = 0.0
      RP(8) = 0.0
      RP(9) = 0.0
C
C FORM COEFFICIENTS OF NORMAL EQUATIONS
C
      DO 21 I = 1,MM
      P(I) = 0.0
      Y(I) = 0.0
      DO 20 J = 1,MM
      X(I,J) = 0.0
   20 CONTINUE
   21 CONTINUE
      DO 28 I = 1,N
      THETA = T(I)*PI/180.0
      P(1) = FLOAT(IH(I)*IH(I))
      P(2) = FLOAT(IK(I)*IK(I))
      P(3) = FLOAT(IL(I)*IL(I))
      IF (M.EQ.1.OR.M.EQ.5.OR.M.EQ.8) P(4) = FLOAT(IK(I)*IL(I))
      IF (M.LE.2.OR.M.EQ.5) P(5) = FLOAT(IL(I)*IH(I))
      IF (M.EQ.1.OR.M.EQ.4.OR.M.EQ.5.OR.M.EQ.8) P(6)= FLOAT(IH(I)*IK(I))
      IF (M.LE.3) GO TO 23
      P(1) = P(1) + P(2)
      P(2) = 0.0
      IF (M.EQ.6) GO TO 23
      IF (M.GE.5) GO TO 22
      P(1) = P(1) + P(6)
      P(6) = 0.0
      GO TO 23
   22 P(1) = P(1) + P(3)
      P(3) = 0.0
      IF (M.EQ.7) GO TO 23
      P(4) = P(4) + P(5) + P(6)
      P(5) = 0.0
      P(6) = 0.0
   23 IF (L.GT.0.AND.L.LT.4) THEN
         DO 36 LI = 1,L
         P(6+LI) = DEL(THETA,LI)
   36    CONTINUE
      ELSE IF (L.EQ.4) THEN
         P(7) = DEL(THETA,4)
      END IF
      QE(I)=(2.0*DSIN(THETA)/ALG)**2
      DE(I)=1.0/DSQRT(QE(I))
      QP(I) = 0.0
      DO 24 J = 1,6
      QP(I) = QP(I) + RP(J)*P(J)
   24 CONTINUE
      P(10) = (QE(I) - QP(I))*100000.0
      DO 37 J = 1,10
      PWTD(J) = P(J)*WT(I)
   37 CONTINUE
      DO 26 JJ = 1,MM
      DO 25 KK = JJ,MM
      X(JJ,KK) = X(JJ,KK) + PWTD(JJ)*PWTD(KK)
   25 CONTINUE
      Y(JJ) = Y(JJ) + PWTD(JJ)*PWTD(10)
   26 CONTINUE
      DO 27 J = 1,10
      PP(I,J) = P(J)
   27 CONTINUE
   28 CONTINUE
      DO 30 I = 1,MM
      DO 29 J = 1,MM
      X(J,I) = X(I,J)
   29 CONTINUE
   30 CONTINUE
      F = 1.0
      CALL TEST1(X,F,M,L)
      CALL COFAC(X,COF,9)
      DELTA = 0.0
      DO 31 J = 1,MM
      DELTA = DELTA + COF(1,J)*X(1,J)
   31 CONTINUE
      DO 33 I = 1,MM
      DO 32 J = 1,MM
      COF(J,I) = COF(I,J)
   32 CONTINUE
   33 CONTINUE
      F = 0.0
      CALL TEST1(COF,F,M,L)
      DO 35 I = 1,MM
      RD(I) = 0.0
      DO 34 J = 1,MM
      RD(I) = RD(I) + COF(I,J)*Y(J)
   34 CONTINUE
      RD(I) = RD(I)/(DELTA*100000.0)
   35 CONTINUE
      IF (M.EQ.8) THEN
         DO 39 I=1,6
         RD(I) = 0.0
   39    CONTINUE
      END IF
C
C CALCULATE REFINED CELL DIMENSIONS
C
      CALL TEST2(M,RD)
      DO 40 I = 1,MM
      RO(I) = RP(I) + RD(I)
   40 CONTINUE
      DO 41 I = 1,3
      RX(I) =DSQRT(RO(I))
   41 CONTINUE
      Z(1) = RO(4)/(2.0*RX(2)*RX(3))
      Z(2) = RO(5)/(2.0*RX(3)*RX(1))
      Z(3) = RO(6)/(2.0*RX(1)*RX(2))
      IF(M.EQ.4)Z(3)=-1.0*DCOS(CD(6))
      DO 42 I = 1,3
      RX(I+3) = DACOS(Z(I))
      IF (RX(I+3).LE.0.0) RX(I+3) = RX(I+3) + PI
   42 CONTINUE
      CALL RECIP(RX,CO,RVOL,VOL,RW)
      IF (M.LT.8) THEN
         WRITE (IW,65)
         WRITE (IW,57)
         WRITE (IW,63) (CO(I), I = 1,3)
      END IF
      EE = 0.0
      DO 44 I = 1,N
      QO = 0.0
      DO 43 J = 1,6
      QO = QO + PP(I,J)*RO(J)
   43 CONTINUE
      IF (M.LT.8) THEN
         DR(I) = 1.0/DSQRT(QO)
      ELSE
         DR(I) = 1.0/DSQRT(QP(I))
      END IF
      TT(I) = DASIN(ALG/(2.0*DR(I)))
      QCOR = QE(I) - (PP(I,7)*RD(7) + PP(I,8)*RD(8) + PP(I,9)*RD(9))
      E = QO - QCOR
      IF (L.GE.0) TTCOR(I) = 2.0*DASIN(SQRT(QCOR)*ALG/2.0)
      EE = EE + WT(I)*E*E
      VO(I) = -E*100000.0
      ED(I) = DR(I) - DE(I)
   44 CONTINUE
      EE = EE/REAL(NM)
      F = EE/DELTA
      CALL ERR(CO,F,COF,SD,VOL,M,MM)
      IF (M.LT.8) WRITE (IW,66) (SD(I), I = 1,3)
      DO 45 I = 4,6
      CO(I) = DACOS(CO(I))
      IF (CO(I).LE.0.0) CO(I) = CO(I) + PI
      CD(I) = CD(I)*180.0/PI
      CO(I) = CO(I)*180.0/PI
      SD(I) = SD(I)*180.0/PI
   45 CONTINUE
      CONST1 = ALG*ALG/4.0
      DO 38 I=7,9
      RD(I) = CONST1*RD(I)
      SD(I) = CONST1*SD(I)
   38 CONTINUE
      IF (M.EQ.4) CO(4) = 90.0
      IF (M.EQ.4) CO(5) = 90.0
C
C OUTPUT DATA
C
      IF (M.EQ.8) GO TO 72
      WRITE (IW,64) (CO(I), I = 4,6)
      WRITE (IW,67) (SD(I), I = 4,6)
      WRITE (IW,68) VOL,SD(10)
      IF (M.GE.4) GO TO 70
      AB = CO(1)/CO(2)
      CB = CO(3)/CO(2)
      SDAB = AB*SQRT((SD(1)/CO(1))**2 + (SD(2)/CO(2))**2)
      SDCB = CB*SQRT((SD(3)/CO(3))**2 + (SD(2)/CO(2))**2)
      WRITE (IW,69) AB,CB
      WRITE (IW,85) SDAB,SDCB
      GO TO 72
   70 IF (M.EQ.7.OR.M.EQ.5) GO TO 72
      CA = CO(3)/CO(1)
      SDCA = CA*SQRT((SD(1)/CO(1))**2 + (SD(3)/CO(3))**2)
      WRITE (IW,71) CA
      WRITE (IW,86) SDCA
   72 IF (L.EQ.4) THEN
         RD(8) = RD(7)
         RD(7) = 0.0
         SD(8) = SD(7)
         SD(7) = 0.0
      END IF
      IF (L.GT.0) THEN
         A = -360.0*RD(7)/PI
         B = -360.0*RD(8)/PI
         C = -360.0*RD(9)/PI
         SDA = 360.0*SD(7)/PI
         SDB = 360.0*SD(8)/PI
         SDC = 360.0*SD(9)/PI
      END IF
      WRITE (IW,'(/A)') ' Coefficients of Correction Function A + BCosTh
     1eta + CSin2Theta are:'
      WRITE (IW,'(3(A,F7.4))') ' A      = ',A,' B      = ',B,
     &' C      = ',C
      WRITE (IW,'(3(A,F7.4)/)') ' S.D. A = ',SDA,' S.D. B = ',SDB,
     &' S.D. C = ',SDC
      IF (L.GT.0) THEN
         CONST2 = 360.0/PI
         ZERO = -CONST2*RD(7)
         SDZERO = CONST2*SD(7)
         SSD = RADIUS*RD(8)
         SDSSD = RADIUS*SD(8)
         IF (L.EQ.3) THEN
            CONST3 = 0.25/RADIUS
            EMU = CONST3/RD(9)
            SDEMU = CONST3*SD(9)/(RD(9)*RD(9))
         END IF
      ELSE
         ZERO = ZEROC
         SSD = SSDC
         EMU = EMUC
      END IF
      WRITE (IW,'(A,F6.3,A,A,F6.3,A)') ' Zero = ',ZERO,' Deg 2Th, ',
     1'S.D. Zero = ',SDZERO,' Deg 2Th'
      WRITE (IW,'(A,F6.3,A,A,F6.3,A)') ' SSD  = ',SSD,' mm,',
     1'      S.D. SSD  = ',SDSSD,' mm'
      WRITE (IW,'(A,F6.1,A,A,F6.1,A/)') ' mu   = ',EMU,' 1/mm',
     1'     S.D. mu   = ',SDEMU,' 1/mm'
      WRITE (IW,73)
      WRITE (IW,74)
      WRITE (IW,75) N
      WRITE (IW,76)
      WRITE (IW,77)
      SET = 0.0
      SDET = 0.0
      VOSUM = 0.0
      DO 79 I = 1,N
      IPP = -INT(PP(I,10))
      IVO = -INT(VO(I))
      VOSUM = VOSUM - VO(I)
      IF (L.GE.0) TTCOR(I) = 180.0*TTCOR(I)/PI
      TTCAL = TT(I)*360.0/PI
      ET(I) = TTCAL - TTCOR(I)
      SET = SET + ET(I)
      SDET = SDET + ET(I)**2
      WRITE (IW,78)IH(I),IK(I),IL(I),TTOBS(I),TTCOR(I),TTCAL,ET(I),
     1IVO,DE(I),DR(I),ED(I)
   79 CONTINUE
      SDET = DSQRT(SDET/(N-1))
      WRITE (IW,80)
      WRITE (IW,83) SET
      WRITE (IW,84) SDET
      WRITE (IW,82) VOSUM
      STOP
   50 FORMAT (1X//(A))
   51 FORMAT (33X,',')
   52 FORMAT (1X,'Refined Cell Dimensions by Cohen s Method of Least',
     1' Squares')
   53 FORMAT (1X,58('-')/)
   54 FORMAT (1X,'Wavelength = ',F8.6,' A'/)
   57 FORMAT (1X,23('-'))
   55 FORMAT (1X,15('-'))
   60 FORMAT (1X,'Cell Dimensions')
   62 FORMAT (1X,'Trial  Cell  Dimensions')
   63 FORMAT(1X,'a     = ',F8.5,' A,    b    = ',F8.5,' A,    c     =
     1 ',F8.5,' A,')
   64 FORMAT (1X,'Alpha = ',F7.3,' Deg.,  Beta = ',F7.3,' Deg.,  Gamma =
     1 ',F7.3,' Deg.')
   65 FORMAT (/1X,'Refined Cell Dimensions')
   66 FORMAT (1X,'Std. Errors in a, b and c are',F8.5,',',F8.5,' and',
     1F8.5,' A')
   67 FORMAT (1X,'Std. Errors in Alpha, Beta and Gamma are',F6.3,',',
     1F6.3,' and',F6.3,' Deg.'/)
   68 FORMAT (1X,'Volume = ',F8.3,' Cu A,',1X,'Std. Error = ',F7.3,
     1 ' Cu A')
   69 FORMAT (1X,'a/b = ',F8.5,'  c/b = ',F8.5)
   71 FORMAT (1X,'c/a = ',F8.5)
   73 FORMAT (/1X,'d-Spacings, etc.')
   74 FORMAT (1X,16('-'))
   75 FORMAT (1X,'n = ',I3)
   76 FORMAT (50X,'5')
   77 FORMAT (3X,'h',2X,'k',2X,'l',2X,'2THobs',2X,'2THcor',2X,'2THcal',
     12X,'Del2TH',3X,'VcX10',3X,'Do A',6X,'Dc A',6X,'ED A'/)
   78 FORMAT (1X,3I3,3F8.3,F7.3,I7,3F10.5,3X,F5.2)
   80 FORMAT (1X/)
   82 FORMAT (1X,'Sum Vc X 100,000  =',F6.1/)
   83 FORMAT (1X,'Sum Delta 2Theta  = ',F6.3)
   84 FORMAT (1X,'S.D. Delta 2Theta = ',F6.3)
   85 FORMAT (1X,'S.D. a/b = ',F7.5,', S.D. c/b = ',F7.5)
   86 FORMAT (1X,'S.D. c/a = ',F7.5)
******
******
******
10000 print *,' File I/O Error'
22000 END
******
******
******
C/1
C CORRECTION FUNCTION (THETA,L)
C
      FUNCTION DEL(THETA,L)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      STT = SIN(2.0*THETA)
      IF (L.EQ.1) DEL = STT
      IF (L.EQ.2.OR.L.EQ.4) DEL = STT*COS(THETA)
      IF (L.EQ.3) DEL = STT*STT
      RETURN
      END
C/2
C SUBROUTINE FOR EVALUATING COFACTORS,  A(I,J) = A(J,I)
C
      SUBROUTINE COFAC(X,COF,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(9,9),XX(9,9),COF(9,9)
      DO 202 I = 1,N
      DO 201 J = I,N
      CALL COPY(XX,X,N)
      DO 200 L = 1,N
      XX(I,L) = 0.0
      XX(L,J) = 0.0
      XX(I,J) = 1.0
  200 CONTINUE
      COF(I,J) = DET(XX,N)
  201 CONTINUE
  202 CONTINUE
      RETURN
      END
C/3
C SUBROUTINE FOR CALCULATING RECIPROCAL CELL DIMENSIONS
C
      SUBROUTINE RECIP(CD,RC,V,RV,W)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CD(6),RC(6)
      PI=4.0D0*DATAN(1.0D0)
      A = CD(1)
      B = CD(2)
      C = CD(3)
      SINAL =DSIN(CD(4))
      SINBE =DSIN(CD(5))
      SINGA =DSIN(CD(6))
      COSAL =DCOS(CD(4))
      COSBE =DCOS(CD(5))
      COSGA =DCOS(CD(6))
      W =DSQRT(1.0 - COSAL**2 - COSBE**2 - COSGA**2 + 2.0*COSAL*COSBE
     1 *COSGA)
      V = A*B*C*W
      RV = 1.0/V
      RC(1) = SINAL/(A*W)
      RC(2) = SINBE/(B*W)
      RC(3) = SINGA/(C*W)
      RC(4) = (COSBE*COSGA - COSAL)/(SINBE*SINGA)
      RC(5) = (COSGA*COSAL - COSBE)/(SINGA*SINAL)
      RC(6) = (COSAL*COSBE - COSGA)/(SINAL*SINBE)
      RETURN
      END
C/4
C SUBROUTINE FOR COPYING X(I,J)
C
      SUBROUTINE COPY(Y,X,MM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(9,9),X(9,9)
      DO 401 I = 1,MM
      DO 400 J = 1,MM
      Y(I,J) = X(I,J)
  400 CONTINUE
  401 CONTINUE
      RETURN
      END
C/5
C EVALUATION OF DETERMINANT (X,N)
C
      FUNCTION DET(ARRAY,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARRAY(1)
      DO 1 I=1,N
      DO 2 J=1,N
      K=N*(J-1)+I
      IF(ARRAY(K))1,2,1
    2 CONTINUE
      DET=0
      RETURN
    1 CONTINUE
      DO 3 J=1,N
      DO 4 I=1,N
      K=N*(J-1)+I
      IF(ARRAY(K))3,4,3
    4 CONTINUE
      DET=0
      RETURN
    3 CONTINUE
      NM=N-1
      DO 5 J=1,NM
      K=N*(J-1)+J
      IF(ARRAY(K))10,6,10
    6 K=J+1
      DO 7 I=K,N
      K=N*(J-1)+I
      IF(ARRAY(K))8,7,8
    7 CONTINUE
    8 DO 9 JJ=1,N
      K=N*(JJ-1)+J
      L=N*(JJ-1)+I
    9 ARRAY(K)=ARRAY(K)+ARRAY(L)
   10 JP=J+1
      DO 5 I=JP,N
      K=N*(I-1)+J
      L=N*(J-1)+J
      F=-ARRAY(K)/ARRAY(L)
      DO 13 K=JP,N
      L=N*(I-1)+K
      M=N*(J-1)+K
   13 ARRAY(L)=F*ARRAY(M)+ARRAY(L)
    5 CONTINUE
      DETS = 1.0
      DO 14 I=1,N
      K=N*(I-1)+I
   14 DETS = ARRAY(K)*DETS
      DET = DETS
      RETURN
      END
C/6
C SUBROUTINE FOR CALCULATING STANDARD ERRORS
C
      SUBROUTINE ERR(CO,F,COF,SD,VOL,M,MM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION COF(9,9),SD(10),CX(9,9),COV(9,9),KK(3),S(3),Z(3),ZZ(3),
     1A(6),CO(6)
C
      IF (M.EQ.8) GO TO 617
      DO 600 I = 1,3
      A(I) = CO(I)
      Z(I) = CO(I+3)
      ZZ(I) =DSQRT(1.0 - Z(I)**2)
  600 CONTINUE
C
C FORM DERIVATIVES
      KK(1) = 4
      KK(2) = 6
      KK(3) = 5
      DO 603 I = 1,3
      DO 601 J = 1,3
      L = KK(J)
      S(J) = A(J)
      IF (J.NE.I) S(J) = S(J)*Z(L-3)
  601 CONTINUE
      K = 2
      L = 3
      DO 602 J = 1,3
      CX(I,J) = -0.5*(S(J)**2)*A(I)
      CX(I,J+3) = -0.5*S(K)*S(L)*A(I)
      CALL TEST3(K)
      CALL TEST3(L)
  602 CONTINUE
      NN    = KK(1)
      KK(1) = KK(2)
      KK(2) = KK(3)
      KK(3) = NN
  603 CONTINUE
      DD = 1.0 - Z(1)**2 - Z(2)**2 - Z(3)**2 + 2.0*Z(1)*Z(2)*Z(3)
      IJ = 2
      IK = 3
      IL = 0
      DO 606 I = 1,3
      A(I+3) = DD/(ZZ(I)*(Z(I) - Z(IJ)*Z(IK)))
      DO 605 J = 1,3
      CALL TEST3(IL)
      IF (I.NE.J) GO TO 604
      CX(I+3,J) = -0.5*(A(I)**2)*(Z(IK)*(Z(IJ) - Z(IK)*Z(I))
     1 + Z(IJ)*(Z(IK) - Z(IJ)*Z(I)))/ZZ(I)
      CX(I+3,J+3) = -0.5*A(IJ)*A(IK)*ZZ(I)
      GO TO 605
  604 CX(I+3,J) = -0.5*(A(IL)**2)*ZZ(I)*Z(I)
      IM = IJ
      IF (I.EQ.1.AND.J.EQ.2) IM = IK
      IF (I.EQ.2.AND.J.EQ.3) IM = IK
      IF (I.EQ.3.AND.J.EQ.1) IM = IK
      CX(I+3,J+3) = -0.5*A(I)*ZZ(I)*A(IM)*Z(IM)
  605 CONTINUE
      CALL TEST3(IJ)
      CALL TEST3(IK)
  606 CONTINUE
  617 CONTINUE
      CX(7,7) = 1.0
      CX(8,8) = 1.0
      CX(9,9) = 1.0
C
C FORM COVARIANCES
      DO 610 IL = 1,MM
      DO 609 IM = IL,MM
      C = 0.0
      DO 608 II = 1,MM
      DO 607 IJ = II,MM
      B = CX(IL,II)*CX(IM,IJ)*COF(II,IJ)
      IF (II.NE.IJ) B = 2.0*B
      C = C + B
  607 CONTINUE
  608 CONTINUE
      COV(IL,IM) = C
  609 CONTINUE
  610 CONTINUE
      IF (M.LE.3) GO TO 613
      DO 612 I = 1,3
      DO 611 J = I,3
      IF (CO(J).EQ.CO(1)) COV(I,J) = COV(1,1)
      IF (M.EQ.5) COV(I+3,J+3) = COV(4,4)
  611 CONTINUE
  612 CONTINUE
      IF (M.EQ.4.OR.M.EQ.6) COV(2,3) = COV(1,3)
C
C FORM STANDARD DEVIATIONS
  613 DO 614 I = 1,MM
      SD(I) =DSQRT(F*COV(I,I))
  614 CONTINUE
      IF (M.EQ.4) SD(6) = 0.0
C SD in Volume
      C = 0.0
      DO 616 I = 1,6
      DO 615 J = I,6
      IF (DABS(A(I)).GT.1.0D10.OR.DABS(A(J)).GT.1.0D10) GO TO 615
      IF (DABS(A(I)).LT.1.0D-10.OR.DABS(A(J)).LT.1.0D-10) GO TO 615
      B = COV(I,J)/(A(I)*A(J))
      IF (I.NE.J) B = 2.0*B
      C = C + B
  615 CONTINUE
  616 CONTINUE
      SD(10)=VOL*DSQRT(DABS(F*C))
      RETURN
      END
C/7
C SYSTEM TEST SUBROUTINE NO.1
C
      SUBROUTINE TEST1(X,T,M,L)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(9,9)
      IF (L.EQ.4) THEN
         LL = 1
      ELSE
         LL = L
      END IF
      IF (M.EQ.1) GO TO 700
      IF (M.EQ.8) X(1,1) = T
      IF (M.GE.4) X(2,2) = T
      IF (M.EQ.5.OR.M.EQ.7.OR.M.EQ.8) X(3,3) = T
      IF (M.GE.2.AND.M.NE.5) X(4,4) = T
      IF (M.GE.3) X(5,5) = T
      IF (M.GE.2) X(6,6) = T
  700 IF (LL.LT.3) X(9,9) = T
      IF (LL.LT.2) X(8,8) = T
      IF (LL.LT.1) X(7,7) = T
      RETURN
      END
C/8
C SYSTEM TEST SUBROUTINE NO.2
C
      SUBROUTINE TEST2(M,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(9)
      IF (M.LE.3) GO TO 800
      IF (M.EQ.5.OR.M.EQ.7) A(3) = A(1)
      IF (M.GE.4) A(2) = A(1)
      IF (M.EQ.5) A(5) = A(4)
      IF (M.EQ.5) A(6) = A(4)
  800 RETURN
      END
C/9
C TEST SUBROUTINE NO.3
C
      SUBROUTINE TEST3(I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      I = I + 1
      IF (I.GT.3) I = 1
      RETURN
      END


