[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Reply to: [list | sender only]
Re: thermal expansion
- Subject: Re: thermal expansion
- From: Kristin <[email protected]>
- Date: Thu, 21 Mar 2002 21:28:01 GMT
Thanks Paul - I can handle that attempt:
I am sending it through the list in case someone else is interested.
Kristin
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
"Paul D. Boyle" wrote:
> Jolanta Holband wrote:
> > Dear Kristin,
> >
> > Unfortunately, your second try to attach the program was
> > also unsuccessful. I have not found "ALPHA" in the net. If
> > you put this file on your server and I will download it from
> > there myself (from your www page or by ftp).
>
> FORTRAN source files are plain text. There is no need to attach them.
> Just use your text editor to read the "alpha.f" file directly into the body
> of your message and send it. Alternatively, if you want to save space
> you can compress(1) or gzip(1) the file, and then attach it (compressed
> files are binary and need to be MIME encoded when emailed.
>
> Paul
>
> --
> Paul D. Boyle | [email protected]
> Director, X-ray Structural Facility | phone: (919) 515-7362
> Department of Chemistry - Box 8204 | FAX: (919) 515-5079
> North Carolina State University |
> Raleigh, NC, 27695-8204
> http://laue.chem.ncsu.edu/web/xray.welcome.html
Reply to: [list | sender only]
- Prev by Date: Re: thermal expansion
- Next by Date: Re: thermal expansion
- Prev by thread: Re: thermal expansion
- Next by thread: Re: thermal expansion
- Index(es):