C
C
C     ALPHA.FOR   11.90
C
C
C     MINERALOGICAL INSTITUTE            DEPARTMENT OF CRYSTALLOGRAPHY
C     UNIVERSITY OF KIEL                 JESSEN,KUEPPERS
C     GERMANY                            (1989)
C
C #################################################################
C
C
C     CALCULATES THERMAL EXPANSION TENSOR FROM BRAGG-ANGLE SHIFTS
C
C     1.) OPTIONAL : POLYNOM-FIT OF D-VALUES
C
C     2.) CALCULATION OF LATTICE CONSTANTS (A,B,C,ALPHA,BETA,GAMMA)
C
C                  *  *  *      *     *      *
C         LSQ --> A, B, C, ALPHA, BETA, GAMMA
C         (1/D(HKL))**2=(H*AS+K*BS+L*CS)**2
C
C     3.) CALCULATION OF TENSORCOMPONENTS  ALPHA(I,J)
C
C                      3   3
C         ALPHA'(1,1)=SUM(SUM(A(1,I)*A(1,J)*ALPHA(I,J)))
C                      I   J
C                    =(DELTA D)/(D*DELTA(T))
C
C     4.) DETERMINATION OF PRINCIPAL AXES
C
C     5.) OPTIONAL : POLYNOM-FIT OF PRINCIPAL VALUES
C
C ################################################################
C
C     INPUT:
C
C     1. CARD (A80)  : TITLE
C
C     2. CARD (6I5)  : CONTROL PARAMETERS
C
C              NT    : NUMBER OF TEMPERATURES(MAX 20)
C              I2T   : 2, IF 2-THETA-VALUES ARE USED, ELSE 0
C              NR    : NUMBER OF REFLEXES (MAX 8000)
C         200->8000 TMARTIN 1997)
C
C              IPOL  : IF =/= 0 --> D-VALUES  FITTED BY
C                      POLYNOM OF N TH ORDER (MAX 3)
C              ICRYS : CRYSTAL SYSTEM
C                      1 : TRICLINIC
C                      2 : MONOCLINIC
C                      3 : ORTHORHOMBIC
C                      4 : TETRAGONAL
C                      5 : TRIGONAL/HEXAGONAL
C              IFIT  : IF =/= 0 --> PRINCIPAL VALUES  FITTED BY
C                      POLYNOM OF N TH ORDER (MAX 3)
C
C     3. CARD (1F)   : LAMBDA
C
C     4. CARD (F)    : TEMPERATURES
C
C     5. CARD (3I4)  : H K L
C
C     6. CARD (F)    : THETA-VALUES FOR DIFF. TEMPERATURES
C
C     REPEAT CARDS 5,6 FOR EACH H K L
C
C     7. CARD (A4,3I4) OPTIONAL : OMIT H K L  (REPEAT, IF NECESSARY)
C
C     LAST CARD      : END
C
C
C     MEANING OF VARIABLES :
C     SOL0 : RECIPROCAL LATTICE CONSTANTS
C     SOL1 : LATTICE CONSTANTS
C     SOL2 : TENSOR COMPONENTS
C     SOL3 : EIGENVALUES
C     SOL4 : EIGENVECTORS
C     SOL5 : STANDARD DEVIATIONS OF TENSOR COMPONENTS
C     SOL6 : RADII OF STEREOGRAPHIC PROJECTION OF EIGENVECTORS
C     SOL7 : ANGLES OF STEREOGRAPHIC PROJECTION OF EIGENVECTORS
C
      DIMENSION IREF(8000,3),THETA(20,8000),D(20,8000),DEL(20,8000)
     1          ,IFKJ(20),SOL0(20,8),SOL1(20,8),SOL2(20,8),AT(8000,4)
     2          ,W(3),Z(3,3),A(8000,6),B(8000),TP(20),SOL5(20,7)
     3          ,LSG(6),SOL3(20,3,2),SOL4(20,3,3,2)
     4          ,V(6,6),FEHLE(6),SOL6(20,3,2),SOL7(20,3,2)
     5          ,VA(3),VN(3),DW(3),DZ(3,3),DVA(3)
     6          ,DVN(3),SW(3),SZ(3,3),SVA(3),SVN(3),TITL(20)
      REAL LSG
      READ(51,87)TITL
      WRITE(53,86)TITL
      READ(51,89)NT   ,I2T  ,NR   ,IPOL ,ICRYS,IFIT
      WRITE(53,89)NT   ,I2T  ,NR   ,IPOL ,ICRYS,IFIT
      IF(I2T  .EQ.0)I2T  =1
      READ(51,*)RLAM
      WRITE(53,*)RLAM
      READ(51,*)(TP(K),K=1,NT)
      WRITE(53,*)(TP(K),K=1,NT)
      DO 20 I=1,NR
      READ(51,101)(IREF(I,J),J=1,3)
      READ(51,*)(THETA(J,I),J=1,NT)
   20 CONTINUE
      DO 1 IFT=1,NT
      DO 2 N=1,NR
      D(IFT,N)=RLAM/(2.*SIN(THETA(IFT,N)/I2T  *ASIN(1.0)/90.))
    2 CONTINUE
    1 CONTINUE
      IPAR3=1
C
C   OMIT-ROUTINE:
C
  200 READ(51,88)TEXT,IRH,IRK,IRL
      IF (TEXT.EQ.'END ') GOTO 201
      IOR=1
  202 IFK=((IREF(IOR,1)-IRH)**2+(IREF(IOR,2)-IRK)**2
     1+(IREF(IOR,3)-IRL)**2)
      IF(IFK.EQ.0) GOTO 207
      IOR=IOR+1
      IF((NR   -IOR)+1.EQ.0) GOTO 200
      GOTO 202
  207 DO 205 IX=IOR,NR
      DO 206 IY=1,3
      IREF(IX,IY)=IREF(IX+1,IY)
  206 CONTINUE
      DO 204 IZV=1,NT
      D(IZV,IX)=D(IZV,IX+1)
      THETA(IZV,IX)=THETA(IZV,IX+1)
  204 CONTINUE
  205 CONTINUE
      NR   =NR   -1
      GOTO 200
C
C ******** END OF OMIT ; **********************
C
C     D-VALUE FITTING
  201 IF(IPOL .EQ.0)GOTO 71
      WRITE(53,4715)IPOL
 4715 FORMAT(//'  D-VALUES FITTED BY POLYNOM OF ORDER',I3)
      CALL POLYNOM(TP,NT   ,IPOL ,D,NR   ,DEL,AT)
  71  IF(IPAR3.EQ.1)GOTO 4020
      IF(IPAR3.EQ.3)GOTO 4020
      GOTO 4068
 4020 CALL OUT2  (IREF,NR   ,THETA,IPAR3,0,1,IPOL ,DEL)
      IF(NT   -10) 4068,4068,4014
 4014 CALL OUT2  (IREF,NR   ,THETA,IPAR3,10,1,IPOL ,DEL)
 4068 CALL OUT2  (IREF,NR   ,D,IPAR3,0,0,IPOL ,DEL)
      IF (NT   -10) 3999,3999,4012
 4012 CALL OUT2  (IREF,NR   ,D,IPAR3,10,0,IPOL ,DEL)
 3999 DO 7 IFL=1,NT
      IF(IFL.EQ.1)GOTO 6
C
C     B: ALPHA'(I,J)
C
      IF(IPOL .EQ.0)GOTO 300
      DO 310 M=1,NR
      T=(TP(IFL)+TP(IFL-1))/2.
      B(M)=(AT(M,2)+2*AT(M,3)*T+3*AT(M,4)*T*T)/(D(IFL,M)+D(IFL-1,M))*2.
  310 CONTINUE
      GOTO 350
  300 DO 8 MB=1,NR
      B(MB)=(D(IFL,MB)-D(IFL-1,MB))/(((D(IFL-
     11,MB)+D(IFL,MB))*(TP(IFL)-TP(IFL-1)))/2.)
    8 CONTINUE
C
C     COEFFICIENT  MATRIX A
C     AEE,AEZ,AED : DIRECTION COSINES U(1,1),U(1,2),U(1,3)
C
  350 DO 14 M=1,NR
      AEE=D(IFL-1,M)*VOLS/BS*(IREF(M,1)*CR-IREF(M,3)*AR*CBER)
      AEZ=D(IFL-1,M)*(IREF(M,1)*AS*CGAS+IREF(M,2)*BS+
     &                                  IREF(M,3)*CS*CALS)
      AED=D(IFL-1,M)*IREF(M,3)/CR
      A(M,1)=AEE**2
      A(M,2)=AEZ**2
      A(M,3)=AED**2
      A(M,4)=2.*AEZ*AED
      A(M,5)=2.*AEE*AED
      A(M,6)=2.*AEE*AEZ
   14 CONTINUE
      CALL LSQ (NR,ICRYS,A,B,LSG,FEHLE,0)
C
C     LSG(1) = ALPHA(1,1)
C         2          2,2
C         3          3,3
C         4          2,3
C         5          1,3
C         6          1,2
C
C     STANDARDDEVIATION  OF  TENSORCOMPONENTS  -> SOL5
C
      DO 738 K=1,NR
      D(IFL-1,K)=B(K)*1000000.
  738 CONTINUE
      DO 739 IL=1,6
      SOL5(IFL,IL)=FEHLE(IL)
      SOL2(IFL,IL)=LSG(IL)
  739 CONTINUE
C
    6 DO 3 MA=1,NR
      DO 4 NA=1,3
      A(MA,NA)=FLOAT(IREF(MA,NA))**2
    4 CONTINUE
      A(MA,4)=2.*FLOAT(IREF(MA,1)*IREF(MA,2))
      A(MA,5)=2.*FLOAT(IREF(MA,1)*IREF(MA,3))
      A(MA,6)=2.*FLOAT(IREF(MA,2)*IREF(MA,3))
    3 CONTINUE
      DO 5 MB=1,NR
      B(MB)=1./D(IFL,MB)**2
    5 CONTINUE
      CALL LSQ (NR,ICRYS,A,B,LSG,FEHLE,1)
C
C     S : STAR  = RECIPROCAL LATTICE
C
C     LSG(1) = AS**2
C         2    BS**2
C         3    CS**2
C         4    AS*BS*COS(GAMMA S)
C         5    AS*CS*COS(BETA S)
C         6    BS*CS*COS(ALPHA S)
C
      AS=SQRT(ABS(LSG(1)))
      BS=SQRT(ABS(LSG(2)))
      CS=SQRT(ABS(LSG(3)))
      CALS=LSG(6)/(BS*CS)
      CBES=LSG(5)/(AS*CS)
      CGAS=LSG(4)/(AS*BS)
      IF(ICRYS.EQ.5) CGAS=0.5
C
      SOL0(IFL,1)=AS
      SOL0(IFL,2)=BS
      SOL0(IFL,3)=CS
      SOL0(IFL,4)=ACOS(CALS)*90./ASIN(1.0)
      SOL0(IFL,5)=ACOS(CBES)*90./ASIN(1.0)
      SOL0(IFL,6)=ACOS(CGAS)*90./ASIN(1.0)
C
C      R : REAL
C
      VOLS=AS*BS*CS*SQRT(1.+2.*CALS*CBES*CGAS-CALS**2-CBES**2-CGAS**2)
      SOL0(IFL,7)=VOLS
      SOL1(IFL,7)=1./VOLS
      AR=BS*CS*SQRT(1.-CALS**2)/VOLS
      BR=CS*AS*SQRT(1.-CBES**2)/VOLS
      CR=AS*BS*SQRT(1.-CGAS**2)/VOLS
      CALR=(CBES*CGAS-CALS)/(SQRT(1.-CBES**2)*SQRT(1.-CGAS**2))
      CBER=(CGAS*CALS-CBES)/(SQRT(1.-CGAS**2)*SQRT(1.-CALS**2))
      CGAR=(CALS*CBES-CGAS)/(SQRT(1.-CALS**2)*SQRT(1.-CBES**2))
C
      SOL1(IFL,1)=AR
      SOL1(IFL,2)=BR
      SOL1(IFL,3)=CR
      SOL1(IFL,4)=ACOS(CALR)*90./ASIN(1.0)
      SOL1(IFL,5)=ACOS(CBER)*90./ASIN(1.0)
      SOL1(IFL,6)=ACOS(CGAR)*90./ASIN(1.0)
C
    7 CONTINUE
C
      DO 5750 I57=1,20
      DO 5750 J57=1,7
      SOL2(I57,J57)=SOL2(I57,J57)*1000000.
      SOL5(I57,J57)=SOL5(I57,J57)*1000000.
 5750 CONTINUE
      IF(ICRYS.GT.2)GOTO 4000
C
C     EIGENVALUES AND EIGENVECTORS
C     IN:         SOL2(20,6)
C                 SOL5(20,6)
C                 NT
C     OUT:        SOL3(20,3,2)
C                 SOL4(20,3,3,2)
C                 SOL6(20,3,2)
C                 SOL7(20,3,2)
C
      DO 5999 IZZL=2,NT
      DO 5002 IXG=1,3
      V(IXG,IXG)=SOL2(IZZL,IXG)
 5002 CONTINUE
      V(1,2)=SOL2(IZZL,6)
      V(1,3)=SOL2(IZZL,5)
      V(2,3)=SOL2(IZZL,4)
      V(2,1)=V(1,2)
      V(3,1)=V(1,3)
      V(3,2)=V(2,3)
C
C     V(I,J)=ALPHA(I,J)
C
      CALL EIGEN (V,W,Z,VA,VN)
      DO 5005 I=1,3
      SW(I)=0.
      SVA(I)=0.
      SVN(I)=0.
      DO 5005 J=1,3
      SZ(I,J)=0.
 5005 CONTINUE
      DO 5010 IP=1,6
      SOL2(IZZL,IP)=SOL2(IZZL,IP)+SOL5(IZZL,IP)
      DO 5007 IXG=1,3
      V(IXG,IXG)=SOL2(IZZL,IXG)
 5007 CONTINUE
      V(1,2)=SOL2(IZZL,6)
      V(1,3)=SOL2(IZZL,5)
      V(2,3)=SOL2(IZZL,4)
      V(2,1)=V(1,2)
      V(3,1)=V(1,3)
      V(3,2)=V(2,3)
      CALL EIGEN (V,DW,DZ,DVA,DVN)
      DO 5030 I=1,3
      DO 5020 J=1,3
      SZ(I,J)=SZ(I,J)+(ABS(Z(I,J))-ABS(DZ(I,J)))**2
 5020 CONTINUE
      SW(I)=SW(I)+(W(I)-DW(I))**2
      SVA(I)=SVA(I)+(VA(I)-DVA(I))**2
      IF((DVN(I)-VN(I)).GT. 100.)DVN(I)=DVN(I)-180.
      IF((DVN(I)-VN(I)).LT.-100.)DVN(I)=DVN(I)+180.
      SVN(I)=SVN(I)+(VN(I)-DVN(I))**2
 5030 CONTINUE
      SOL2(IZZL,IP)=SOL2(IZZL,IP)-SOL5(IZZL,IP)
 5010 CONTINUE
C
C     W : EIGENVALUES, Z : MATRIX OF EIGENVECTORS
C     VA : RADII, VN : ANGLES OF STEREOGRAPHIC PROJECTION
C
      DO 5003 ICG=1,3
      SOL3(IZZL,ICG,1)=W(ICG)
      SOL6(IZZL,ICG,1)=VA(ICG)
      SOL7(IZZL,ICG,1)=VN(ICG)
      SOL3(IZZL,ICG,2)=SQRT(SW(ICG))
      SOL6(IZZL,ICG,2)=SQRT(SVA(ICG))
      SOL7(IZZL,ICG,2)=SQRT(SVN(ICG))
      DO 5004 ICV=1,3
      SOL4(IZZL,ICG,ICV,1)=Z(ICG,ICV)
      SOL4(IZZL,ICG,ICV,2)=SQRT(SZ(ICG,ICV))
 5004 CONTINUE
 5003 CONTINUE
 5999 CONTINUE
C
C
C     OUTPUT-ROUTINE:
C
C     SOL0 = RECIPROCAL A,B,C,...
C     SOL1 =            A,B,C,...
C
 4000 CALL OUT1  (TP,SOL1,1,IPAR3,NT   ,SOL5)
      IF(IPAR3.EQ.0)GOTO 4022
      IF(IPAR3.EQ.2)GOTO 4022
      CALL OUT1  (TP,SOL0,0,IPAR3,NT   ,SOL5)
 4022 CALL OUT1  (TP,SOL2,2,IPAR3,NT   ,SOL5)
      IF(ICRYS.LT.3)CALL OUT3  (SOL3,SOL4,SOL6,SOL7,NT   ,IPAR3)
      IF(IFIT.EQ.0)GOTO 4030
      DO 4025 I=1,19
      TP(I)=(TP(I)+TP(I+1))/2.
      DO 4025 J=1,3
      THETA(I,J)=SOL3(I+1,J,1)
 4025 CONTINUE
      CALL POLYNOM(TP,NT-1,IFIT,THETA,3,DEL,AT)
      WRITE(53,4027)IFIT
 4027 FORMAT(//'  PRINCIPAL VALUES FITTED BY POLYNOM OF ORDER',I3)
      DO 4028 I=1,3
      WRITE(53,4029)I,I,AT(I,1),AT(I,2),AT(I,3),AT(I,4)
 4028 CONTINUE
 4029 FORMAT(/'  ALPHA(',I1,',',I1,') =',F10.3,' + ',F10.5,'*T + ',
     &F10.7,'*(T**2) + ',F12.9,'*(T**3)')
 4030 DO 4035 J=1,20
      IFKJ(J)=J
 4035 CONTINUE
      WRITE(53,513)
      WRITE(53,506)(IFKJ(J),J=1,NT   -1)
      DO 4049 KFL=1,NR
      WRITE(53,514)KFL,(D(J,KFL),J=1,NT   -1)
 4049 CONTINUE
      WRITE(53,99)
   86 FORMAT(1H ,20A4)
   87 FORMAT(20A4)
   88 FORMAT(1A4,3I4)
   99 FORMAT(/' END OF ALPHA-EXECUTION')
   89 FORMAT(6I5)
  101 FORMAT(3I4)
  500 FORMAT(A5)
  504 FORMAT(1H )
  505 FORMAT(1H*)
  506 FORMAT(4H*NR.,20(1I3,3H.  ))
  507 FORMAT(1H*,131(1H-))
  508 FORMAT(1H*,1I3,20F6.0,F8.0)
  513 FORMAT(/' RESIDUALS OF ALPHA''(1,1) '/)
  514 FORMAT(1H*,1I5,20F6.1)
      STOP
      END
C******************************************************************
      SUBROUTINE OUT1  (TPM,SOLU,NI1,NI2,NI3,SOL5)
      DIMENSION SOLU(20,8),TPM(20),SOL5(20,7)
      IF(NI2-1)2000,2000,2001
 2000 WRITE(53,2101)
      IF(NI1.EQ.2)GOTO 2306
      GOTO 2305
 2306 WRITE(53,2108)
      GOTO 2314
 2305 IF(NI1.EQ.0)WRITE(53,2102)
      WRITE(53,2103)
      DO 2200 LAUF=1,NI3
      IPH=INT(TPM(LAUF))
      WRITE(53,2106)IPH,(SOLU(LAUF,IIJ),IIJ=1,7)
 2200 CONTINUE
      WRITE(53,2101)
      GOTO 2999
 2001 WRITE(53,2100)
      IF(NI1-2)2303,2302,2303
 2302 WRITE(53,2108)
      GOTO 2314
 2303 IF(NI1.EQ.0)WRITE(53,2102)
      WRITE(53,2103)
 2314 WRITE(53,2104)
      IF(NI1.EQ.2)GOTO 2500
      DO 2300 LAUF=1,NI3
      IPH=INT(TPM(LAUF))
      WRITE(53,2105)IPH,(SOLU(LAUF,IJJ),IJJ=1,7)
      WRITE(53,2107)
 2300 CONTINUE
      GOTO 2499
 2500 DO 2502 LAUF=2,NI3
      IPH=INT(TPM(LAUF-1))
      WRITE(53,2109)IPH
      LAU25=LAUF-1
      WRITE(53,2110)LAU25,(SOLU(LAUF,IJ),IJ=1,6)
      WRITE(53,2111)(SOL5(LAUF,IJ),IJ=1,6)
      IPH=INT(TPM(LAUF))
      WRITE(53,2109)IPH
 2502 CONTINUE
 2499 WRITE(53,2104)
 2999 RETURN
 2100 FORMAT(1H1///)
 2101 FORMAT(1H*)
 2102 FORMAT(11H*          ,3(14H       *      ),4(14H         *    ))
 2103 FORMAT(9H* TEMP   ,16H        A       ,14H     B        ,
     114H      C       ,14H    ALPHA     ,14H    BETA      ,
     214H    GAMMA     ,14H   VOLUME     )
 2104 FORMAT(1H ,126(1H-))
 2105 FORMAT(1I6,3H  !,7(3H   E11.5))
 2106 FORMAT(1H*,1I5,3H   ,7(3H   F11.5))
 2107 FORMAT(9H        !)
 2108 FORMAT(9H* TEMP   ,16H    ALPHA(1,1)  ,14H  ALPHA(2,2)  ,14
     1H  ALPHA(3,3)  ,14H  ALPHA(2,3)  ,14H  ALPHA(1,3)  ,
     214H  ALPHA(1,2)  )
 2109 FORMAT(1I6,3H  !)
 2110 FORMAT(1I3,5(1H ),1H!,6(3H   F11.2))
 2111 FORMAT(8H        ,1H!,6(3H   F11.2),(10H          F10.2))
      END
C********************************************************************
      SUBROUTINE OUT2  (IRD,IP4,AW,IP1,IP2,IP3,IP6,DEL)
      DIMENSION IRD(8000,3),AW(20,8000),IP2A(10),DEL(20,8000)
      DO 1004 LJ=1,10
      IP2A(LJ)=LJ+IP2
 1004 CONTINUE
      IF(IP1-2)1000,1000,1001
 1000 WRITE(53,1101)
      IF(IP3)1006,1005,1006
 1005 WRITE(53,1109)
      GOTO 1007
 1006 WRITE(53,1108)
 1007 WRITE(53,1102)(IP2A(J),J=1,10)
      DO 1002 LI=1,IP4
      WRITE(53,1103)(IRD(LI,JI),JI=1,3),(AW(JK,LI),JK=1+IP2,10+IP2)
      IF((IP3.EQ.0).AND.(IP6.NE.0))WRITE(53,1113)
     &(DEL(JK,LI),JK=1+IP2,10+IP2)
 1002 CONTINUE
      GOTO 1999
 1001 WRITE(53,1100)
      IF(IP3)1009,1008,1009
 1008 WRITE(53,1111)
      GOTO 1010
 1009 WRITE(53,1110)
 1010 WRITE(53,1107)
      WRITE(53,1102)(IP2A(JJ),JJ=1,10)
      WRITE(53,1104)
      DO 1003 LI=1,IP4
      WRITE(53,1106)(IRD(LI,JI),JI=1,3),(AW(JK,LI),JK=IP2+1,IP2+10)
      IF((IP3.EQ.0).AND.(IP6.NE.0))WRITE(53,1116)
     &(DEL(JK,LI),JK=1+IP2,10+IP2)
      WRITE(53,1105)
 1003 CONTINUE
      WRITE(53,1104)
 1999 RETURN
 1100 FORMAT(1H1///)
 1101 FORMAT(1H*)
 1102 FORMAT(18H*   H    K    L   ,10(6H     T,I2,3H   ))
 1103 FORMAT(1H*,1I4,2I5,3H   ,10F11.5)
 1104 FORMAT(128(1H-))
 1105 FORMAT(16(1H ),1H!)
 1106 FORMAT(3I5,3H ! ,10F11.5)
 1107 FORMAT(///)
 1108 FORMAT(16H* THETA-LISTING:)
 1109 FORMAT(16H* D-LISTING    :)
 1110 FORMAT(16H  THETA-LISTING:)
 1111 FORMAT(16H  D-LISTING    :)
 1113 FORMAT(18X,10F11.5)
 1116 FORMAT(15X,3H ! ,10F11.5)
      END
C*********************************************************************
      SUBROUTINE OUT3  (SOLT,SO,RAD,WIN,NQ,NP)
      DIMENSION SOLT(20,3,2),SO(20,3,3,2),RAD(20,3,2),WIN(20,3,2)
      WRITE(53,6105)
      WRITE(53,6106)
      DO 6010 IYF=2,NQ
      WRITE(53,6101)
      WRITE(53,6102)SOLT(IYF,1,1),SOLT(IYF,1,2),(SO(IYF,1,IUY,1),
     &SO(IYF,1,IUY,2),IUY=1,3),RAD(IYF,1,1),RAD(IYF,1,2),WIN(IYF,1,1),
     &WIN(IYF,1,2)
      WRITE(53,6104)IYF-1,SOLT(IYF,2,1),SOLT(IYF,2,2),(SO(IYF,2,IUY,1),
     &SO(IYF,2,IUY,2),IUY=1,3),RAD(IYF,2,1),RAD(IYF,2,2),WIN(IYF,2,1),
     &WIN(IYF,2,2)
      WRITE(53,6102)SOLT(IYF,3,1),SOLT(IYF,3,2),(SO(IYF,3,IUY,1),
     &SO(IYF,3,IUY,2),IUY=1,3),RAD(IYF,3,1),RAD(IYF,3,2),WIN(IYF,3,1),
     &WIN(IYF,3,2)
 6010 CONTINUE
      RETURN
 6100 FORMAT(120(1H-))
 6101 FORMAT(1H*,19(1H ),1H!)
 6102 FORMAT(1H*,5X,F6.1,1H(,F4.1,4H)  !,
     13(F9.4,1H(,F5.4,1H)),2(F7.1,1H(,F4.1,1H)))
 6103 FORMAT(1H1//)
 6104 FORMAT(1H*,1I3,2X,F6.1,1H(,F4.1,4H)  !,
     13(F9.4,1H(,F5.4,1H)),2(F7.1,1H(,F4.1,1H)))
 6105 FORMAT(1H*)
 6106 FORMAT(4H*NR.,3X,11HEIGENVALUES,3H  !,3(1H ),
     113HEIGENVECTORS ,36X,2HR',9X,3HETA,5X,'FOR  STEREOGR. ',
     2'PROJECTION (R=10CM)')
      END
      SUBROUTINE INVERS(A,N)
C     GAUSS-JORDAN-RUTISHAUSER MATRIX INVERSION WITH DOUBLE PIVOTING.
      DIMENSION A(6,6),B(20),C(20),P(20),Q(20)
      INTEGER P,Q
      XPS=0.1E-10
      MSING=1
      DO 10 K=1,N
C     DETERMINATION OF THE PIVOT ELEMENT
      PIVOT=0.
      DO 20 I=K,N
      DO 20 J=K,N
      IF(ABS(A(I,J))-ABS(PIVOT))20,20,30
   30 PIVOT=A(I,J)
      P(K)=I
      Q(K)=J
   20 CONTINUE
      IF(ABS(PIVOT)-XPS) 40,40,50
C     EXCHANGE OF THE PIVOTAL ROW WITH THE KTH ROW
   50 IF(P(K)-K)60,80,60
   60 DO 70 J=1,N
      L=P(K)
      Z=A(L,J)
      A(L,J)=A(K,J)
   70 A(K,J)=Z
C     EXCHANGE OF THE PIVOTAL COLUMN WITH THE KTH COLUMN
   80 IF(Q(K)-K)85,90,85
   85 DO 100 I=1,N
      L=Q(K)
      Z=A(I,L)
      A(I,L)=A(I,K)
  100 A(I,K)=Z
   90 CONTINUE
C     JORDAN STEP
      DO 110 J=1,N
      IF(J-K)130,120,130
  120 B(J)=1./PIVOT
      C(J)=1.
      GO TO 140
  130 B(J)=-A(K,J)/PIVOT
      C(J)=A(J,K)
  140 A(K,J)=0.
  110 A(J,K)=0.
      DO 10 I=1,N
      DO 10 J=1,N
   10 A(I,J)=A(I,J)+C(I)*B(J)
C     REORDERING THE MATRIX
      DO 155 M=1,N
      K=N-M+1
      IF(P(K)-K)160,170,160
  160 DO 180 I=1,N
      L=P(K)
      Z=A(I,L)
      A(I,L)=A(I,K)
  180 A(I,K)=Z
  170 IF(Q(K)-K) 190,155,190
  190 DO 150 J=1,N
      L=Q(K)
      Z=A(L,J)
      A(L,J)=A(K,J)
  150 A(K,J)=Z
  155 CONTINUE
  151 RETURN
   40 WRITE(53,45)P(K),Q(K),PIVOT
   45 FORMAT(16H0SINGULAR MATRIX3H I=I3,3H J=I3,7H PIVOT=E16.8/)
      MSING=2
      GO TO 151
      END
C************************************************************************
      SUBROUTINE POLYNOM(TP,IPAR,IK,SOL,IVA,DEL,AT)
      DIMENSION A(8000,6),TP(20),Q(6,6),PH(6,8000),LSG(6),
     &SOL(20,8000),DEL(20,8000),AT(8000,4)
      REAL LSG
      IG=IK+1
      DO 10999 IP=1,IPAR
      DO 10100 IQ=1,IG
      A(IP,IQ)=1.0
      DO 10101 I=2,IQ
      A(IP,IQ)=A(IP,IQ)*TP(IP)
10101 CONTINUE
10100 CONTINUE
10999 CONTINUE
      DO 10422 L=1,IG
      DO 10421 M=1,IG
      Q(M,L)=0.0
      DO 10420 N=1,IPAR
      Q(M,L)=Q(M,L)+A(N,M)*A(N,L)
10420 CONTINUE
10421 CONTINUE
10422 CONTINUE
      CALL INVERS(Q,IG)
      DO 10932 L=1,IPAR
      DO 10831 K=1,IG
      PH(K,L)=0.0
      DO 10830 M=1,IG
      PH(K,L)=PH(K,L)+Q(K,M)*A(L,M)
10830 CONTINUE
10831 CONTINUE
10932 CONTINUE
      DO 10888 IVAR=1,IVA
      DO 10041 J=1,IG
      LSG(J)=0.0
      DO 10040 L=1,IPAR
      LSG(J)=LSG(J)+PH(J,L)*SOL(L,IVAR)
10040 CONTINUE
10041 CONTINUE
      DO 10500 IAT=1,4
      AT(IVAR,IAT)=LSG(IAT)
10500 CONTINUE
      DO 10887 IVG=1,IPAR
      DEL(IVG,IVAR)=LSG(1)+LSG(2)*TP(IVG)+LSG(3)*(TP(IVG))**2
     &+LSG(4)*(TP(IVG))**3-SOL(IVG,IVAR)
      SOL(IVG,IVAR)=SOL(IVG,IVAR)+DEL(IVG,IVAR)
10887 CONTINUE
10888 CONTINUE
      RETURN
      END
C************************************************************************
C
C     CALCULATION OF STEREOGRAPHIC PROJECTION FROM CARTESIAN COORDINATES
C
      SUBROUTINE STEREOGR (K,R,W)
      INTEGER I
      REAL K,R,W
      DIMENSION K(3,3),R(3),W(3)
      DO 10 I=1,3
      IF(K(I,3).GT.0.)GOTO 5
      DO 3 J=1,3
      K(I,J)=-K(I,J)
   3  CONTINUE
   5    W(I) = (ATAN(K(I,2)/K(I,1)))*180./3.1415926
        IF (K(I,1).LT.0.0) W(I) = W(I)-180.0
        IF (W(I).LT.0.0) W(I) = W(I)+360.0
        R(I) = 10.0*TAN(ACOS(K(I,3))/2)
   10 CONTINUE
      RETURN
      END
C*********************************************************************
C     ROUTINE ACCORDING TO NYE J. F. : PHYSICAL PROPERTIES OF CRYSTALS
C*********************************************************************
      SUBROUTINE EIGEN (V,W,Z,VA,VN)
      DIMENSION V(6,6),W(3),Z(3,3),VA(3),VN(3),AL(3,3),V1(3),V2(3)
      DO 5009 I50=1,3
      DO 5009 J50=1,3
      AL(I50,J50) = V(I50,J50)
 5009 CONTINUE
      NZ50=-1
 5010 VA(1) = 1.0
      VA(2) = 0.0
      VA(3) = 0.0
 5011 DO 5012 I50=1,3
      VN(I50) = 0.0
      DO 5012 J50=1,3
      VN(I50) = VN(I50) + V(I50,J50)*VA(J50)
 5012 CONTINUE
      BETR = SQRT(VN(1)**2 + VN(2)**2 + VN(3)**2)
      SU = 0.0
      DO 5013 I50=1,3
      VN(I50) = VN(I50) / BETR
      SU = SU + ABS(VN(I50)) - ABS(VA(I50))
      SU = ABS(SU)
 5013 CONTINUE
      SUN = SU - 0.0002
      DO 5014 I50=1,3
      VA(I50) = VN(I50)
 5014 CONTINUE
      IF (SUN) 5015,5015,5011
 5015 DO 5016 I50=1,3
      VN(I50) = 0.0
 5016 CONTINUE
      DO 5017 I50=1,3
      DO 5017 J50=1,3
      VN(I50) = VN(I50) + AL(I50,J50)*VA(J50)
 5017 CONTINUE
      LAU50=NZ50 + 2
      HA1 = SQRT(VN(1)**2 + VN(2)**2 + VN(3)**2)
      DO 5118 I50=1,3
      VN(I50) = VN(I50) / HA1
 5118 CONTINUE
      IF(VA(1)*VN(1)+VA(2)*VN(2)+VA(3)*VN(3))   5219,5220,5220
 5219 HA1=-HA1
 5220 W(LAU50)=HA1
      DO 5119 I50=1,3
      Z(LAU50,I50) = VN(I50)
 5119 CONTINUE
      IF (NZ50) 5201,5203,5203
 5201 DO 5202 I50=1,3
      V1(I50) = VA(I50)
 5202 CONTINUE
      GO TO 5205
 5203 DO 5204 I50=1,3
 5204 V2(I50) = VA(I50)
 5205 CALL INVERS(V,3)
      NZ50 = NZ50 + 1
      IF (NZ50) 5010,5010,5019
 5019 VA(1) = V1(2)*V2(3) - V1(3)*V2(2)
      VA(2) = V1(3)*V2(1) - V1(1)*V2(3)
      VA(3) = V1(1)*V2(2) - V1(2)*V2(1)
      DO 5020 I50=1,3
      VN(I50) = 0.0
 5020 CONTINUE
      DO 5021 I50=1,3
      DO 5021 J50=1,3
      VN(I50) = VN(I50) + AL(I50,J50)*VA(J50)
 5021 CONTINUE
      HA1 = SQRT(VN(1)**2 + VN(2)**2 + VN(3)**2)
      IF(VA(1)*VN(1)+VA(2)*VN(2)+VA(3)*VN(3)) 5022,5023,5023
 5022 HA1=-HA1
 5023 LAU50 = 3
      DO 5120 I50=1,3
      VN(I50) = VN(I50) / HA1
 5120 CONTINUE
      W(LAU50) = HA1
      DO 5121 I50=1,3
      Z(LAU50,I50)=VN(I50)
 5121 CONTINUE
      IF (W(1).LT.W(2)) GOTO 6010
      X=W(1)
      W(1)=W(2)
      W(2)=X
      DO 6000 I=1,3
      X=Z(1,I)
      Z(1,I)=Z(2,I)
      Z(2,I)=X
 6000 CONTINUE
 6010 IF (W(1).LT.W(3)) GOTO 6030
      X=W(1)
      W(1)=W(3)
      W(3)=X
      DO 6020 I=1,3
      X=Z(1,I)
      Z(1,I)=Z(3,I)
      Z(3,I)=X
 6020 CONTINUE
 6030 IF (W(2).LT.W(3)) GOTO 6050
      X=W(2)
      W(2)=W(3)
      W(3)=X
      DO 6040 I=1,3
      X=Z(2,I)
      Z(2,I)=Z(3,I)
      Z(3,I)=X
 6040 CONTINUE
 6050 CALL STEREOGR(Z,VA,VN)
      RETURN
      END
C***********************************************************************
      SUBROUTINE LSQ (NY,NSYM,A,Y,XC,SIGMA,IE)
      DIMENSION A(8000,6),PH(6,8000),Y(8000),DIV(8000)
      DIMENSION XC(6),SIGMA(6),V(6,6)
      NX=2
      IF(NSYM.EQ.1)NX=6
      IF(NSYM.EQ.2)NX=4
      IF(NSYM.EQ.3)NX=3
      DO 10 M=1,NY
      IF(NSYM.EQ.1)GOTO 10
      AK=A(M,4)
      A(M,4)=A(M,5)
      IF(NSYM.LT.4)GOTO 10
      A(M,1)=A(M,1)+A(M,2)
      A(M,2)=A(M,3)
      IF((NSYM.EQ.4).OR.(IE.EQ.0))GOTO 10
      A(M,1)=A(M,1)+AK/2.
 10   CONTINUE
      DO 20 L=1,NX
      DO 20 M=1,NX
      V(M,L)=0.0
      DO 20 N=1,NY
      V(M,L)=V(M,L)+A(N,M)*A(N,L)
 20   CONTINUE
      CALL INVERS(V,NX)
      DO 30 L=1,NY
      DO 30 K=1,NX
      PH(K,L)=0.0
      DO 30 M=1,NX
      PH(K,L)=PH(K,L)+V(K,M)*A(L,M)
 30   CONTINUE
      DO 41 J=1,NX
      XC(J)=0.0
      SIGMA(J)=0.0
      DO 40 L=1,NY
      XC(J)=XC(J)+PH(J,L)*Y(L)
 40   SIGMA(J)=SIGMA(J)+PH(J,L)*PH(J,L)
 41   CONTINUE
      F=0.0
      DO 60 I=1,NY
      DIV(I)=0.0
      DO 50 J=1,NX
 50   DIV(I)=DIV(I)+A(I,J)*XC(J)
      DIV(I)=DIV(I)-Y(I)
 60   F=F+DIV(I)*DIV(I)
      F=SQRT(F/(NY-NX))
      DO 70 K=1,NX
 70   SIGMA(K)=F*SQRT(SIGMA(K))
      DO 80 K=1,NY
 80   Y(K)=DIV(K)
      IF(NSYM.EQ.1)GOTO 100
      XC(5)=XC(4)
      SIGMA(5)=SIGMA(4)
      XC(4)=0.0
      SIGMA(4)=0.0
      IF(NSYM.LT.4)GOTO 100
      XC(3)=XC(2)
      SIGMA(3)=SIGMA(2)
      XC(2)=XC(1)
      SIGMA(2)=SIGMA(1)
 100  RETURN
      END