C
C     CALCULATION OF DISTANCE AND ANGLE
C                  STANDARD DEVIATIONS OF BOND LENGTH AND ANGLES
C         FOR THE IBM 7090     R. SNYDER   MARCH 1969
C  CONVERTED TO 360/75  12/9/69
C        PROGRAM WITH SYMMETRY AND TRANSLATION
C
C  ITEM 1 CONTROL CARD   (3I5)
C        COL 5 = ITAP = 0  INPUT DATA FROM TAPE 9 PREPARED BY ORFLS
C                     = 1 INPUT DATA FROM CARDS
C  NOTE THE NEXT TWO QUANTITIES (ISYMT AND ISYMC) MAY BOTH BE ZERO OR
C  MAY BOTH BE NON ZERO.  IF BOTH ARE ZERO THEN THE SYMMETRY SELECTION
C  NUMBER ON THE SPECIFICATION CARD MAY BE LEFT BLANK AND AN IDENTITY
C  SYMMETRY ELEMENT WILL BE USED BY PROGRAM.
C  IF BOTH ARE NON ZERO THEN ISYMT ELEMENTS WILL BE READ FROM TAPE AND
C  ISYMC ELEMENTS WILL BE READ FROM CARDS.  OF COURSE ALL OR NONE MAY
C  BE READ FROM EITHER TAPE OR CARDS
C        COL  6-10 = ISYMT = THE NUMBER OF SYMMETRY TO BE INPUT FROM TAP
C                          = -1  ALL SYMMETRY ELEMENTS USED IN L.S. WILL
C                          =  0 SYMMETRY ELEMENTS ON TAPE IGNORED
C                          =  N  THE FIRST N ELEMENTS FROM L.S. WILL BE
C        COL 11-15 = ISYMC = THE NUMBER OF SYMMETRY TO BE INPUT FROM CAR
C  REMBER THAT THE SYM. ELEMENTS USED IN ORFLS DO NOT INCLUDE CENTROSYMM
C        OR LATTICE TRANSLATION RELATIONS.  IF THE COMPLETE SET IS
C        DESIRED THE ADDITIONAL (OR ALL) SHOULD BE READ FROM CARDS
C        BY MAKING ISYMC  NON ZERO
C
C
C*****  IF ITAP = 0 AND ISYMC = 0 THEN ITEM 6 DIRECTLY FOLLOWS ITEM 1 **
C   IF ITAP= 0 AND ISYMC = N THEN ORDER IS ITEM1,ITEM3,ITEM6
C
C
C
C  ITEM 2  CELL CONSTANTS  (3F7.3,3F7.4) ONLY USED WHEN ITAP = 1
C        A, B, C, COS(ALPHA), COS(BETA), COS(GAMMA)
C
C  ITEM 3   ISYM SYMMETRY CARDS IN STANDARD SHIONO FORMAT
C        COL 1-11  TRANSLATION FOR X                       WHERE
C        COL 13  SYMMETRY FOR X                            / = BLANK
C        COL 15  SYMMETRY FOR X (IF TWO TERMS ARE NEEDED)  X= 1
C        COL 16-26 TRANSLATION FOR Y                       -X = -1
C        COL 28  SYMMETRY FOR Y                            Y= 2
C        COL 30  SECOND SYMMETRY FOR Y(IF NEEDED)          -Y= -2
C        COL 31-41  TRANSLATION FOR Z                      Z = 3
C        COL 43 SYMMETRY FOR Z                             -Z = -3
C        COL 45 SECOND SYMMETRY FOR Z (IF NEEDED)
C
C  ITEM 4  COORDINATES AND ERRORS PUNCHED BY L.S.  (A6,1X1H ,2F13.7)
C        COL 1-6   ATOM NAME
C        COL 8-13  LETTER X,Y, OR Z AS PUNCHED BY L.S.
C        COL 14-27 COORDINATE
C        COL 28-41 STANDARD ERROR
C    COORDINATES MUST BE INPUT IN THE ORDER X, Y, Z
C
C  ITEM 5  TERMINATE COORDINATE CARDS
C        COL 14-27 999.0
C
C  ITEM 6  SPECIFICATION CARDS   (3(I4,1X3I1,I2),I5)
C  THESE SPECIFY THE ATOMS TO BE USED IN THE BOND LENGTH AND ANGLES CALC
C        COL 1-3   I  THE SEQUENTIAL NUMBER OF ATOM DESIRED
C        COL 4-6   CELL DESIRED  EG 555
C                       IF BLANK 555 CELL WILL BE ASSUMED
C        COL 7-8   SYMMETRY ELEMENT NUMBER TO BE APPLIED
C                       IF BLANK A DUMMY WILL BE USED
C        COL 11-13 J  ANOTHER SEQUENTIAL ATOM NUMBER
C        COL 14-16 CELL
C        COL 17-18 SYM. ELEMENT NO.
C        COL 21-23 K  A THIRD SEQUENTIAL ATOM NO.
C                       IF BLANK ONLY DIST I-J WILL BE COMPUTED
C        COL 24-26 CELL
C        COL 27-28 SYM. ELEMENT NO.
C        COL 31-35 =0   NORMAL
C                  =1  BONDS I-J AND J-K ARE CHEMICALLY EQUIVILENT
C
C  ITEM 7  BLANK CARD TO TERMINATE RUN
C
C  THERE IS NO LIMIT TO THE NUMBER OF CONTROL CARDS WHICH MAY BE USED
C  BUT ONLY 200 COORDINATES (ATOMS) MAY BE READ AND ONLY 99 SYMMETRY CAR
C
C
	REAL*8 ALP(200)
	DIMENSION X(200),Y(200),Z(200),SD(200),C1(99),C2(99),C3(99),A1(99)
	1,A2(99),A3(99),B1(99),B2(99),B3(99)
	DIMENSION I5(4),J5(4),K5(4)
	COMMON X,Y,Z,C1,C2,C3,A1,A2,A3,B1,B2,B3
	DATA XL,YL,ZL/'X','Y','Z'/
	DISTF(P,Q,R,S,T,U) = SQRT(P*P + Q*Q + R*R + 2.0*(P*Q*U + Q*R*S +
	1R*P*T))
	NTP=9
	I = 0
C
C  READ CONTROL CARD
C
	READ 1,  ITAP,ISYMT,ISYMC
1	FORMAT(3I5)
	IF(ITAP)1000,1000,1010
C
C  READ CELL CONSTANT RECORD FROM ORFLS TAPE
C
1000	REWIND NTP
	READ(NTP)  CELA,CELB,CELC,COSA,COSB,COSC
	GO TO 1020
C
C  READ CELL CONSTANT CARD
C
1010	READ 5,        CELA,CELB,CELC,COSA,COSB,COSC
5	FORMAT ( 3F7.3,3F7.4 )
1020	PRINT  6,         CELA,CELB,CELC,COSA,COSB,COSC
6	FORMAT ( 3H1A=F7.3, 4H, B=,F7.3,4H, C=F7.3,8H, COS A=F7.4,
	1  8H, COS B= F7.4, 8H, COS C= F7.4  // )
C
C  ESTABLISH SYMMETRY ELEMENT 5 AS IDENTITY DEFAULT VALUE
C
	C1(5) = 0.0
	C2(5) = 0.0
	C3(5) = 0.0
	A1(5) = +1.0
	A2(5) = +1.0
	A3(5) = +1.0
	B1(5) = +1.0
	B2(5) = +1.0
	B3(5) = +1.0
	IF(ITAP)1030,1030,1050
C
C  READ SYMMETRY ELEMENTS FROM ORFLS TAPE
C
1030	READ(NTP) ISYM
	PRINT 1031, ISYM
1031	FORMAT(' THERE ARE ',I5,' SYMMETRY ELEMENTS ON ORFLS TAPE'/)
	DO 1035  K=1,ISYM
1035	READ(NTP) C1(K),A1(K),B1(K),C2(K),A2(K),B2(K),C3(K),A3(K),B3(K)
	IF( ISYMT )1040,1050,1050
1040	ISYMT = ISYM
1050	ISTRT = ISYMT + 1
	ISYM = ISYMT
C
C  READ SYMMETRY CARDS
C
	IF( ISYMC ) 1060,1060,80
80	ISYM = ISYMT + ISYMC
	READ  7, (C1(K),A1(K),B1(K),C2(K),A2(K),B2(K),C3(K),A3(K),B3(K),K=
	1ISTRT,ISYM)
7	FORMAT(3(F11.6,2F2.0))
1060	IF ( ISYM ) 93,93,1070
1070	PRINT 9
9	FORMAT(4X6H TRANS,4X1HX,5X5HTRANS,4X1HY,5X5HTRANS,4X1HZ /)
C
C  SET BLANKS (INTERPRETED AS -0.0) TO +0.0
C
	DO 92 K = 1, ISYM
	IF(A1(K))82,81,82
81	A1(K) = +0.0
82	IF(A2(K))84,83,84
83	A2(K) = +0.0
84	IF(A3(K))86,85,86
85	A3(K) = +0.0
86	IF(B1(K))88,87,88
87	B1(K) = +0.0
88	IF(B2(K))90,89,90
89	B2(K) = +0.0
90	IF(B3(K))92,91,92
91	B3(K) = +0.0
92	CONTINUE
	PRINT17, (C1(K),A1(K),B1(K),C2(K),A2(K),B2(K),C3(K),A3(K),B3(K),K=
	11,ISYM)
17	FORMAT( 3( F8.3,2XF3.0,1XF3.0))
	GO TO 95
93	PRINT 8
8	FORMAT(40H0SYMMETRY ELEMENT 5 ON OUTPUT IS A DUMMY  //)
C
C  READ AND PRINT COORDINATES AND ERRORS
C
95	PRINT 10
10	FORMAT(40H0NO. ATOM     X         Y         Z       ,
	140H  SIG X     SIG Y     SIG Z    SIGMA XYZ  )
100	I = I + 1
	IF(ITAP) 101,101,102
C
C  READ COORDINATES AND ERRORS FROM TAPE
C
101	READ(NTP,END=250) ALP(I) , X(I) , SX , Y(I) , SY , Z(I) , SZ
	GO TO 130
C
C  READ COORDINATE AND ERROR CARDS
C
102	READ 15,ALP(I),XP,X(I),SX
15	FORMAT(A6,1XA1 ,5X2F13.7)
	IF(X(I) - 999.0)110,250,250
110	IF(XP - XL)900,120,900
120	READ 20,YP,Y(I),SY
20	FORMAT(7XA1 ,5X2F13.7)
	IF(YP - YL)900,125,900
125	READ 20 , ZP,Z(I),SZ
	IF(ZP -ZL)900,130,900
130	SD(I) = CELA*SX + CELB*SY + CELC*SZ
	IF(SX)150,150,140
140	M = 1
	GO TO 160
150	M = 0
160	IF(SY)180,180,170
170	M = M + 1
180	IF(SZ)200,200,190
190	M = M + 1
200	IF(M)220,220,210
210	FM = M
	GO TO 230
220	FM = 3.0
230	SD(I) = SD(I) / FM
	PRINT 25,I,ALP(I),X(I),Y(I),Z(I),SX,SY,SZ,SD(I)
25	FORMAT(I4,2XA6,3F9.6,4F10.7)
	SD(I) = SD(I) * SD(I)
	GO TO 100
250	N = I - 1
C
C  PRINT HEADER
C
255	IPAGE = 0
	PRINT 30
30	FORMAT(  96H1   I    ATOM       J    ATOM       K    ATOM   DIST I
	1-J   ESD  DIST J-K   ESD     ANGLE     ESD  /)
C
C  READ CONTROL (DIRECTION) CARDS AND CALCULATE DISTS AND ANGLES
C
260	READ 35,I,(I5(II),II=1,4),J,(J5(II),II=1,4),K,(K5(II),II=1,4),KSM
35	FORMAT(3(I3,3I1,I2,2X),I5)
	IPAGE = IPAGE + 1
	IF(I) 999,999,265
C
C  ADJUST SYMMETRY AND TRANSLATIONS IF ZERO
C
265	DO 270 II = 1,4
	IF(I5(II))271,271,272
271	I5(II) = 5
272	IF(J5(II))273,273,274
273	J5(II) = 5
274	IF(K5(II))275,275,270
275	K5(II) = 5
270	CONTINUE
	CALL TRANS(XI,YI,ZI,I5,I)
	CALL TRANS(XJ,YJ,ZJ,J5,J)
	A = CELA * (XI - XJ)
	B = CELB * (YI - YJ)
	C = CELC * (ZI - ZJ)
	DIJ = DISTF(A,B,C,COSA,COSB,COSC)
	IF(K)290,280,290
C
C     DISTANCE CALCULATION ONLY
C
280	V1 = SQRT( SD(I) + SD(J) )
	PRINT 40 ,I,(I5(II),II=1,4),ALP(I),J,(J5(II),II=1,4),ALP(J),DIJ,V1
40	FORMAT(2(I3,3I1,I2,2XA6),16X2F8.4/)
	GO TO 260
C
C     CALCULATE THE DISTANCE AND ANGLE
C
290	CALL TRANS(XK,YK,ZK,K5,K)
	A = CELA * (XK - XJ)
	B = CELB * (YK - YJ)
	C = CELC * (ZK - ZJ)
	DJK = DISTF(A,B,C,COSA,COSB,COSC)
	A = CELA * (XI - XK)
	B = CELB * (YI - YK)
	C = CELC * (ZI - ZK)
	DIK = DISTF(A,B,C,COSA,COSB,COSC)
	CANG=0.5*(DIJ*DIJ+DJK*DJK-DIK*DIK)/(DIJ*DJK)
	QA=1.0/(DIJ*DIJ)
	QC=1.0/(DJK*DJK)
	QB=QA+QC-2.0*CANG/(DIJ*DJK)
	IF(CANG) 310,300,320
300	ANG = 90.0
	GO TO 350
310	KANG = 1
	GO TO 330
320	KANG = 2
330	ANG = ATAN (SQRT (1.0 - CANG * CANG) / CANG) * 57.29578
	GO TO(340,350),KANG
340	ANG = ANG + 180.0
C
C     CALCULATE STANDARD DEVIATIONS
C
350	V2 = SQRT  ( SD(J) + SD(K) )
	V1 = SQRT  ( SD(I) + SD(J) )
	VA = SQRT  ( SD(I) * QA + SD(J) * QB + SD(K) * QC ) * 57.29578
	IF( KSM ) 370,370,360
360	V1 = SQRT  ( (SD(I) + 2.0 * SD(J) * (1.0 + CANG) + SD(K) ) * .25 )
	V2 = 0.0
370	PRINT 45,I,(I5(II),II=1,4),ALP(I),J,(J5(II),II=1,4),ALP(J),K,
	1(K5(II),II=1,4),ALP(K),DIJ,V1,DJK,V2,ANG,VA
45	FORMAT(3(I3,3I1,I2,2XA6),4F8.4,F10.3,F8.4/)
	IF(IPAGE - 27)260,255,260
900	PRINT 50
50	FORMAT(33H0COORDINATE CARDS OUT OF SEQUENCE   )
999	STOP
	END

	SUBROUTINE TRANS(U,V,W,L5,L)
	DIMENSION X(200),Y(200),Z(200),SD(200),C1(99),C2(99),C3(99),A1(99)
	1,A2(99),A3(99),B1(99),B2(99),B3(99)
	DIMENSION L5(4)
	COMMON X,Y,Z,C1,C2,C3,A1,A2,A3,B1,B2,B3
	II = L5(4)
	XT = L5(1) - 5
	YT = L5(2) - 5
	ZT = L5(3) - 5
	U = C1(II) + X(L)*SIGN (1.,A1(II))*SIGN (1.,B1(II)) + XT
	V = C2(II) + Y(L)*SIGN (1.,A2(II))*SIGN (1.,B2(II)) + YT
	W = C3(II) + Z(L)*SIGN (1.,A3(II))*SIGN (1.,B3(II)) + ZT
	RETURN
	END
