linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
ef5a0354f3edfc92ac90fc462f77df3478bad641
sics/difrac/matrix.f
cvs cbc7fdf334 Compiled under Redhat Linux<Right
2000-02-18 15:54:23 +00:00

235 lines
8.6 KiB
Fortran
Raw Blame History

C-----------------------------------------------------------------------
C Library of matrix operations for crystal geometry
C-----------------------------------------------------------------------
SUBROUTINE MATRIX(A,B,C,D,IWHAT)
COMMON /IOUASS/ IOUNIT(12)
CHARACTER COUT*132
COMMON /IOUASC/ COUT(20)
CHARACTER*6 IWHAT
DIMENSION A(3,3),B(3,3),C(3,3),D(3,3),E(3,3),V(3)
DATA RA/57.29578/
IF (IWHAT .EQ. 'INVERT') GO TO 100
IF (IWHAT .EQ. 'MATMUL') GO TO 120
IF (IWHAT .EQ. 'MATVEC') GO TO 150
IF (IWHAT .EQ. 'VECMAT') GO TO 180
IF (IWHAT .EQ. 'SCALPR') GO TO 210
IF (IWHAT .EQ. 'LENGTH') GO TO 230
IF (IWHAT .EQ. 'ORTHOG') GO TO 250
IF (IWHAT .EQ. 'DETERM') GO TO 270
IF (IWHAT .EQ. 'MVMULT') GO TO 290
IF (IWHAT .EQ. 'VMMULT') GO TO 320
IF (IWHAT .EQ. 'TRNSPS') GO TO 340
IF (IWHAT .EQ. 'SYMOPR') GO TO 370
IF (IWHAT .EQ. 'VECPRD') GO TO 400
IF (IWHAT .EQ. 'COPRIM') GO TO 410
IF (IWHAT .EQ. 'INTRCH') GO TO 440
IF (IWHAT .EQ. 'SUMVEC') GO TO 460
IF (IWHAT .EQ. 'DIFVEC') GO TO 480
IF (IWHAT .EQ. 'COPVEC') GO TO 500
ITP = IOUNIT(6)
WRITE (COUT,10000) IWHAT
CALL GWRITE (ITP,' ')
STOP
C-----------------------------------------------------------------------
C Invert 3x3 matrix A, put the result in B
C-----------------------------------------------------------------------
100 E(1,1) = A(2,2)*A(3,3) - A(2,3)*A(3,2)
E(2,1) = -(A(2,1)*A(3,3) - A(2,3)*A(3,1))
E(3,1) = A(2,1)*A(3,2) - A(2,2)*A(3,1)
E(1,2) = -(A(1,2)*A(3,3) - A(1,3)*A(3,2))
E(2,2) = A(1,1)*A(3,3) - A(1,3)*A(3,1)
E(3,2) = -(A(1,1)*A(3,2) - A(1,2)*A(3,1))
E(1,3) = A(1,2)*A(2,3) - A(1,3)*A(2,2)
E(2,3) = -(A(1,1)*A(2,3) - A(1,3)*A(2,1))
E(3,3) = A(1,1)*A(2,2) - A(1,2)*A(2,1)
DMAT = A(1,1)*E(1,1) + A(1,2)*E(2,1) + A(1,3)*E(3,1)
DO 115 I=1,3
DO 110 J = 1,3
110 B(I,J) = E(I,J)/DMAT
115 CONTINUE
GO TO 520
C-----------------------------------------------------------------------
C Multiply 3x3 matrices A and B, store result in C
C-----------------------------------------------------------------------
120 DO 135 I = 1,3
DO 132 J = 1,3
E(I,J) = 0.0
DO 130 K = 1,3
130 E(I,J) = E(I,J) + A(I,K)*B(K,J)
132 CONTINUE
135 CONTINUE
DO 145 I = 1,3
DO 140 J = 1,3
140 C(I,J) = E(I,J)
145 CONTINUE
GO TO 520
C-----------------------------------------------------------------------
C Multiply matrix A by vector B, store dir. cosines of result in C
C-----------------------------------------------------------------------
150 DO 165 I = 1,3
V(I) = 0.
DO 160 J = 1,3
160 V(I) = V(I) + A(I,J)*B(J,1)
165 CONTINUE
IF(V(1)**2 + V(2)**2 +V(3)**2 .GT. 0) THEN
VMOD = SQRT(V(1)**2 + V(2)**2 + V(3)**2)
ELSE
VMOD = 1
ENDIF
DO 170 I = 1,3
170 C(I,1) = V(I)/VMOD
GO TO 520
C-----------------------------------------------------------------------
C Multiply vector A by matrix B, store dir. cosines of result in C
C-----------------------------------------------------------------------
180 DO 195 I = 1,3
V(I) = 0.
DO 190 J = 1,3
190 V(I) = V(I) + B(J,I)*A(J,1)
195 CONTINUE
VMOD = SQRT(V(1)**2 + V(2)**2 + V(3)**2)
DO 200 I = 1,3
200 C(I,1) = V(I)/VMOD
GO TO 520
C-----------------------------------------------------------------------
C Scalar product of vectors A and B
C-----------------------------------------------------------------------
210 S = 0
DO 220 I = 1,3
220 S = S + A(I,1)*B(I,1)
C(1,1) = S
GO TO 520
C-----------------------------------------------------------------------
C length of vector B when A is the metric matrix
C-----------------------------------------------------------------------
230 DO 245 I = 1,3
V(I) = 0.
DO 240 J = 1,3
240 V(I) = V(I) + A(I,J)*B(J,1)
245 CONTINUE
C(1,1) = SQRT(V(1)**2 + V(2)**2 + V(3)**2)
GO TO 520
C-----------------------------------------------------------------------
C Get the metric matrix C corresponding to cell edges A & angles B
C-----------------------------------------------------------------------
250 COSGAS = (COS(B(1,1)/RA)*COS(B(2,1)/RA) - COS(B(3,1)/RA))
COSGAS = COSGAS/(SIN(B(1,1)/RA)*SIN(B(2,1)/RA))
SINGAS = SQRT(1.0 - COSGAS**2)
E(1,1) = A(1,1)*SIN(B(2,1)/RA)*SINGAS
E(1,2) = 0
E(1,3) = 0
E(2,1) = -A(1,1)*SIN(B(2,1)/RA)*COSGAS
E(2,2) = A(2,1)*SIN(B(1,1)/RA)
E(2,3) = 0
E(3,1) = A(1,1)*COS(B(2,1)/RA)
E(3,2) = A(2,1)*COS(B(1,1)/RA)
E(3,3) = A(3,1)
DO 265 I = 1,3
DO 260 J = 1,3
260 C(I,J) = E(I,J)
265 CONTINUE
GO TO 520
C-----------------------------------------------------------------------
C Calculate the determinant D of the vectors A,B,C
C-----------------------------------------------------------------------
270 DET = 0.
DO 280 I = 1,3
J = I + 1
IF (J .EQ. 4) J = 1
K = 6 - I - J
280 DET = DET + A(I,1)*(B(J,1)*C(K,1) - B(K,1)*C(J,1))
D(1,1) = DET
GO TO 520
C-----------------------------------------------------------------------
C Multiply matrix A by vector B, store result in C
C-----------------------------------------------------------------------
290 DO 305 I = 1,3
E(I,1) = 0
DO 300 J = 1,3
300 E(I,1) = E(I,1) + A(I,J)*B(J,1)
305 CONTINUE
DO 310 I = 1,3
310 C(I,1) = E(I,1)
GO TO 520
C-----------------------------------------------------------------------
C Multiply vector A by matrix B, store result in C
C-----------------------------------------------------------------------
320 DO 335 I = 1,3
C(I,1) = 0.
DO 330 J = 1,3
330 C(I,1) = C(I,1) + A(J,1)*B(J,I)
335 CONTINUE
GO TO 520
C-----------------------------------------------------------------------
Ctranspose matrix A and put it in B
C-----------------------------------------------------------------------
340 DO 355 I = 1,3
DO 350 J = 1,3
350 E(I,J) = A(J,I)
355 CONTINUE
DO 365 I = 1,3
DO 360 J = 1,3
360 B(I,J) = E(I,J)
365 CONTINUE
GO TO 520
C-----------------------------------------------------------------------
C Get the symmetry-equivalent of an atom
C-----------------------------------------------------------------------
370 DO 390 I = 1,3
C(I,1) = 0.
DO 380 J = 1,3
380 C(I,1) = C(I,1) + A(I,J)*B(J,1)
J = 4
390 C(I,1) = C(I,1) + A(I,J)/12.
GO TO 520
C-----------------------------------------------------------------------
C Vector product C = A x B
C-----------------------------------------------------------------------
400 C(1,1) = A(2,1)*B(3,1) - A(3,1)*B(2,1)
C(2,1) = A(3,1)*B(1,1) - A(1,1)*B(3,1)
C(3,1) = A(1,1)*B(2,1) - A(2,1)*B(1,1)
GO TO 520
C-----------------------------------------------------------------------
C Make coprime integers (the smallest non-zero integer will be 1)
C-----------------------------------------------------------------------
410 SMALL = 2.
DO 420 I = 1,3
IF (ABS(A(I,1)) .LE. 0.1 .OR. ABS(A(I,1)) .GE. SMALL) GO TO 420
SMALL = ABS(A(I,1))
420 CONTINUE
DO 430 I = 1,3
INDEX = A(I,1)/SMALL + 0.5
IF (A(I,1) .LT. 0.) INDEX = A(I,1)/SMALL - 0.5
430 B(I,1) = INDEX
GO TO 520
C-----------------------------------------------------------------------
C Interchange two vectors A and B
C-----------------------------------------------------------------------
440 DO 450 I = 1,3
SAVE = A(I,1)
A(I,1) = B(I,1)
450 B(I,1) = SAVE
GO TO 520
C-----------------------------------------------------------------------
C Sum of vectors C = A + B
C-----------------------------------------------------------------------
460 DO 470 I = 1,3
470 C(I,1) = A(I,1) + B(I,1)
GO TO 520
C-----------------------------------------------------------------------
C Vector difference C = A - B
C-----------------------------------------------------------------------
480 DO 490 I = 1,3
490 C(I,1) = A(I,1) - B(I,1)
GO TO 520
C-----------------------------------------------------------------------
C Vector copy B = A
C-----------------------------------------------------------------------
500 DO 510 I = 1,3
510 B(I,1) = A(I,1)
GO TO 520
520 RETURN
10000 FORMAT(' Matrix operation ',A6,' is not programmed')
END
Reference in New Issue View Git Blame Copy Permalink