      PROGRAM VOLCAL
C VOLCAL.FOR
C
C PROGRAM TO CALCULATE POLYHEDRAL VOLUMES AND DISTORTION PARAMETERS
C        WITH ERRORS (IF DESIRED).
C
C ORIGINAL PROGRAM BY L. W. FINGER,  9/21/71
C MODIFIED 10/24/73 BY Y. OHASHI
C MODIFIED 6/20/79  BY L. W. FINGER - MADE INTERACTIVE AND ERRORS ADDED
C
C Modified 10/18/96 by L.W.F. to fix problem with more than 3 points on face.
C
C     *****************************************************************
C     INPUT DATA FORMAT IF READ FROM FILE
C     *****************************************************************
C
C     1. TITLE CARD (A80)
C     2. CELL CARD (6F8.0,2I4)
C        COLS.1-8,9-16,17-24 = A,B AND C
C        COLS.25-32,33-40,41-48 = ALPHA,BETA AND GAMMA IN DEG.
C        COLS.49-52  MP  =  NUMBER OF POLYHEDRA TO BE PROCESSED
C        COLS.53-56  IER =  NON-ZERO IF ERRORS TO BE CALCULATED
C     3. CELL VARIANCE-COVARIANCE MATRIX (6F10.0) (NEEDED ONLY IF IER NON-ZERO)
C        THIS INPUT IS SIX LINES WITH ONE ROW OF MATRIX/LINE.
C     4. CENTRAL ATOM CARD (A4,A2,3F10.5,I4)
C        COLS.1-6 SITE NAME
C        COLS.7-16,17-26,27-36 =  X,Y AND Z
C        COLS.37-40 NP  =  NUMBER OF COORDINATING ATOMS
C     5. CENTRAL ATOM SIGMA CARDS (3F9.5) (NEEDED ONLY IF IER NON-ZERO)
C        COLS.1-9,10-18,19-27 = SIGMA X, Y, AND Z
C     6. COORDINATING ATOM CARDS (A4,A2,3F10.5)
C        COLS. 1-6 SITE NAME
C        COLS.7-16,17-26,27-36 =  X,Y AND Z
C     7. COORDINATING ATOM SIGMA CARDS (3F9.5) NEEDED ONLY IF IER NON-ZERO)
C        COLS.1-9,10-18,19-27 = SIGMA X, Y, AND Z
C
C    CARDS 6 AND 7 SHOULD BE REPEATED FOR EACH COORDINATING ATOM
C
C     IF DESIRED, ANOTHER SET OF DATA CAN FOLLOW
C
C     ******************************************************************
C
        DIMENSION A(6),VARA(6,6),SITE(2),XC(3),SXC(3)
        DIMENSION XP(3,20),SXP(3,20),DERDA(6,3),AX(3),SIGMA(3)
     1,X(3),PRM(3)
        CHARACTER TITLE*80,IYES*1,INO*1,ILYES*1,ILNO*1,IPAGE*1,IRESP*1
        DATA IYES,INO,ILYES,ILNO,IPAGE,IN/'Y','N','y','n',' ',4/
C
90      WRITE(*,1)
1       FORMAT('0POLYHEDRAL VOLUME AND DISTORTION PARAMETER PROGRAM'
     1/' ARE DATA TO BE READ FROM FILE(Y,N)?',$)
        READ(*,'(a)')IRESP
        IFILE=1
        IF(IRESP.EQ.IYES.OR.IRESP.EQ.ILYES)GO TO 1000
        IF(IRESP.NE.INO.AND.IRESP.NE.ILNO)GO TO 90
        IFILE=0
        WRITE(*,2)
2       FORMAT(' TYPE TITLE?',$)
        READ(*,'(A)')TITLE
        WRITE(*,6)
6       FORMAT(' ARE ERRORS TO BE CALCULATED(Y,N)?',$)
100     READ(*,'(A)')IRESP
8       FORMAT(A1)
        IER=1
        IF(IRESP.EQ.IYES.OR.IRESP.EQ.ILYES)GO TO 110
        IF(IRESP.NE.INO.AND.IRESP.NE.ILNO)GO TO 100
        IER=0
110     WRITE(*, 10)
10      FORMAT(' TYPE CELL CONSTANTS?',$)
        READ(*, 12)A
12      FORMAT(6F10.0)
        IF(IER.EQ.0)GO TO 120
        WRITE(*, 14)
14      FORMAT(' TYPE VARIANCE MATRIX FOR CELL DATA (1 ROW PER LINE)')
        READ(*, 12)VARA
120     WRITE(*, 16)
16      FORMAT(' TYPE NO. OF COORDINATING IONS?',$)
        READ(*, 18)NP
18      FORMAT(I6)
        IF(NP.EQ.0)STOP
        WRITE(6,5)IPAGE,TITLE
5       FORMAT(A1,a)
        IPAGE='1'
        WRITE(6,13)A
13      FORMAT('0CELL CONSTANTS',6F10.5)
        IF(IER.NE.0)WRITE(6,15)VARA
15      FORMAT('0CELL COVARIANCE MATRIX'/(6E13.6))
        WRITE(*, 20)
20      FORMAT(' TYPE CATION COORDS. AND SITE NAME?',$)
        READ(*, 22)XC,SITE
22      FORMAT(3F9.0,2A4)
        WRITE(6,23)SITE,XC
23      FORMAT('0CATION: ',2A4,3F10.5)
        IF(IER.EQ.0)GO TO 130
        WRITE(*, 24)
24      FORMAT(' TYPE SIGMAS FOR SITE?',$)
        READ(*, 22)SXC
        WRITE(6,25)SXC
25      FORMAT(5X,'SIGMAS',6X,3F10.5)
130     WRITE(6,27)
27      FORMAT('0ANIONS')
        DO 140 I=1,NP
        WRITE(*, 26)I
26      FORMAT(' TYPE ANION COORDS. AND LABEL FOR NO.',I3,'?',$)
        READ(*, 22)(XP(J,I),J=1,3),SITE
        WRITE(6,29)I,SITE,(XP(J,I),J=1,3)
29      FORMAT(3X,I5,1X,2A4,3F10.5)
        IF(IER.EQ.0)GO TO 140
        WRITE(*, 28)
28      FORMAT(' SIGMAS?',$)
        READ(*, 22)(SXP(J,I),J=1,3)
        WRITE(6,25)(SXP(J,I),J=1,3)
140     CONTINUE
        GO TO 1410
C
C READ DATA FROM FILE
C
1000    WRITE(*, 40)
40      FORMAT(' TYPE INPUT FILE NAME?',$)
        READ(*, 42)TITLE
42      FORMAT(A)
        OPEN(UNIT=IN,FILE=TITLE,STATUS='OLD')
1010    READ(IN,'(A)',END=260)TITLE
        READ(IN,44)A,MP,IER
        NCOUNT=0
44      FORMAT(6F8.0,2I4)
        WRITE(6,5)IPAGE,TITLE
        IPAGE='1'
        WRITE(6,13)A
        IF(IER.EQ.0)GO TO 1020
        READ(IN,12)VARA
        WRITE(6,15)VARA
1020    READ(IN,46)SITE,XC,NP
46      FORMAT(A4,A2,3F10.0,I4)
        WRITE(6,23)SITE,XC
        IF(IER.EQ.0)GO TO 1030
        READ(IN,22)SXC
        WRITE(6,25)SXC
1030    DO 1040 I=1,NP
        READ(IN,46)SITE,(XP(J,I),J=1,3)
        WRITE(6,29)I,SITE,(XP(J,I),J=1,3)
        IF(IER.EQ.0)GO TO 1040
        READ(IN,22)(SXP(J,I),J=1,3)
        WRITE(6,25)(SXP(J,I),J=1,3)
1040    CONTINUE
C
1410    CALL PRMCAL(0,A,A(4),XC,NP,XP,PRM(1),PRM(2))
        IF(IER.EQ.0)GO TO 250
C FUNCTIONS CALULATED - NOW GET ERRORS
        DO 142 I=1,3
142     SIGMA(I)=0.0
        DO 155 I=1,6
        DX=SQRT(VARA(I,I))*.1
        IF(DX.EQ.0.0)GO TO 155
        A(I)=A(I)+DX
        CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2))
        A(I)=A(I)-DX
        DO 150 J=1,3
150     DERDA(I,J)=(X(J)-PRM(J))/DX
155     CONTINUE
        DO 170 I=1,6
        DO 170 J=1,3
        AX(I)=0.0
        DO 160 K=1,6
160     AX(I)=AX(I)+DERDA(I,J)*VARA(K,I)
170     SIGMA(J)=SIGMA(J)+AX(I)*DERDA(I,J)
        DO 225 I=1,3
        DX=XC(I)+.1*SXC(I)
        IF(DX.EQ.0.0)GO TO 225
        XC(I)=XC(I)+DX
        CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2))
        XC(I)=XC(I)-DX
        DO 220 J=1,3
220     SIGMA(J)=SIGMA(J)+((X(J)-PRM(J))/DX*SXC(I))**2
225     CONTINUE
        DO 231 I=1,NP
        DO 231 J=1,3
        DX=XP(J,I)+0.1*SXP(J,I)
        IF(DX.EQ.0.0)GO TO 231
        XP(J,I)=XP(J,I)+DX
        CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2))
        XP(J,I)=XP(J,I)-DX
        DO 230 K=1,3
230     SIGMA(K)=SIGMA(K)+((X(K)-PRM(K))/DX*SXP(J,I))**2
231     CONTINUE
        DO 232 I=1,3
232     SIGMA(I)=SQRT(SIGMA(I))
        IF(NP.NE.4.AND.NP.NE.6)GO TO 240
        WRITE(6,30)SIGMA
30      FORMAT(' SIGMA(S)',3F12.5)
        GO TO 250
240     WRITE(6,30)SIGMA(1)
250     IF(IFILE.EQ.0)GO TO 120
        NCOUNT=NCOUNT+1
        IF(NCOUNT-MP)1020,1010,1010
260     STOP
        END
        SUBROUTINE PRMCAL(IORG,A,ANG1,XC,NP,YP,VOL,QE)
        DIMENSION XC(3),A(3),ANG1(3),ANG(3),XP(3,20),VXYZ(6),D(3)
     1,DM(20),IH(3)
        DIMENSION YP(3,20)
C
      DATA IOUT/6/
      DATA PI/3.141592/
C
    6 FORMAT(1H0,'POLYHEDRAL VOLUME',F13.5/5X,'QUAD. ELONG.',F14.5/5X,
     1'ANGLE VARIANCE',F12.5)
    8 FORMAT (1H0,'POLYHEDRAL VOLUME',F12.5)
   15 FORMAT(1H0,10X,'ATOMS',13X,'DISTANCES',20X,'ANGLE'/12X,'I   J',
     1 7X,'0-I',9X,'0-J',9X,'I-J',7X,'I-0-J'/)
   20 FORMAT(1H ,8X,2I4,3F12.5,F12.3)
   25 FORMAT(1H0,12X,'FACE DEFINED',10X,'PLANE NORMAL',5X,'AREA'/13X,
     1 'BY POINTS',7X,'OX',6X,'OY',6X,'OZ'/)
   26 FORMAT (1H ,10X,3I3,4X,3F8.5,4X,F10.5)
C
      DO 100 I=1,3
  100 ANG(I)=COS(ANG1(I)*PI/180.0)
      SINA=SQRT(1.0-ANG(1)**2)
      ASTAR=SINA/(A(1)*SQRT(1.0-ANG(1)**2-ANG(2)**2-ANG(3)**2+2.0*ANG(1)
     1 *ANG(2)*ANG(3)))
C        CONVERT POINTS TO ORTHOGONAL COORDS.
      DO 130 I=1,NP
C        REFER TO ORIGIN AT CENTRAL ATOM
      DO 120 J=1,3
  120 XP(J,I)=YP(J,I)-XC(J)
      XP(3,I)=A(3)*XP(3,I)+XP(2,I)*A(2)*ANG(1)+XP(1,I)*A(1)*ANG(2)
      XP(2,I)=XP(2,I)*A(2)*SINA+XP(1,I)*A(1)*(ANG(3)-ANG(2)*ANG(1))/SINA
      XP(1,I)=XP(1,I)/ASTAR
      DM(I)=SQRT(XP(1,I)**2+XP(2,I)**2+XP(3,I)**2)
  130 CONTINUE
C     DISTANCES AND ANGLES
      IF(IORG.EQ.0)WRITE (IOUT,15)
      DO 145 I=1,NP-1
      DO 145 J=I+1,NP
      TEMPA=0.
      TEMPB=0.
      DO 135 K=1,3
      TEMPA=TEMPA+(XP(K,I)-XP(K,J))**2
  135 TEMPB=TEMPB+XP(K,I)*XP(K,J)
      DIS=SQRT(TEMPA)
      ANGL=TEMPB/(DM(I)*DM(J))
      SUMTH=1.0-ANGL*ANGL
      IF(SUMTH.LT.0.0)SUMTH=0.0
      ANGL=57.2958*ATAN2(SQRT(SUMTH),ANGL)
      IF(IORG.EQ.0)WRITE (IOUT,20) I,J,DM(I),DM(J),DIS,ANGL
  145 CONTINUE
      IF(IORG.EQ.0)WRITE(IOUT,25)
C        FIND FACES   - CHOOSE ANY THREE POINTS
      VOL=0.0
      NP1=NP-1
      NP2=NP-2
      SUMTH =0.0
      SUMTH2=0.0
      NA=0
      DO 180 I=1,NP2
      IH(1)=I
      I1=I+1
      DO 180 J=I1,NP1
      J1=J+1
      IH(2)=J
      VXYZ(1)=XP(1,J)-XP(1,I)
      VXYZ(2)=XP(2,J)-XP(2,I)
      VXYZ(3)=XP(3,J)-XP(3,I)
      DO 180 K=J1,NP
      IH(3)=K
      NS=0
      VXYZ(4)=XP(1,K)-XP(1,I)
      VXYZ(5)=XP(2,K)-XP(2,I)
      VXYZ(6)=XP(3,K)-XP(3,I)
      D(1)=VXYZ(2)*VXYZ(6)-VXYZ(5)*VXYZ(3)
      D(2)=VXYZ(4)*VXYZ(3)-VXYZ(1)*VXYZ(6)
      D(3)=VXYZ(1)*VXYZ(5)-VXYZ(4)*VXYZ(2)
      AREA=0.5*SQRT(D(1)**2+D(2)**2+D(3)**2)
      Z0=0.5*(XP(1,I)*D(1)+XP(2,I)*D(2)+XP(3,I)*D(3))/AREA
C        CHECK FOR AND AVOID PLANE THROUGH ORIGIN
      IF(ABS(Z0).LT.1.0E-5) GO TO 180
      Factor = 3.0
      DO 170 L=1,NP
      IF(L.EQ.I.OR.L.EQ.J.OR.L.EQ.K) GO TO 170
C        CALCULATE DISTANCE OF POINT L FROM PLANE OF IJK
      Z=0.5*((XP(1,I)-XP(1,L))*D(1)+(XP(2,I)-XP(2,L))*D(2)+(XP(3,I)-XP(3
     1,L))*D(3))/AREA
c
c        Code added 18-Oct-1996
C
C Z and Z0 must have the same sign
c
      if (z * z0 .lt. -0.001)go to 180
      if (abs(z * z0) .lt. 0.001)then
c If more than 3 corners on this face, the area will be counted twice.
c  Change factor to handle this case.
      	factor = 6.0
	endif
c
  170 CONTINUE
c
C        ALL POINTS ON SAME SIDE,  THUS IJK ARE FACE
C
C     DIRECTION COSINES OF PLANE NORMAL
      D(1)=D(1)/(AREA*2.)
      D(2)=D(2)/(AREA*2.)
      D(3)=D(3)/(AREA*2.)
      IF(IORG.EQ.0)WRITE (IOUT,26) I,J,K,D,AREA
      VOL=VOL+AREA*ABS(Z0)/FACTOR
      DO 172 L=1,2
      L1=L+1
      NM=IH(L)
      DO 172 M=L1,3
      MN=IH(M)
      TEMP=0.0
      DO 171 N=1,3
  171 TEMP=TEMP+XP(N,NM)*XP(N,MN)
      TEMP=TEMP/(DM(NM)*DM(MN))
      TEMP=57.2958*ATAN2(SQRT(1.0-TEMP*TEMP),TEMP)
      SUMTH=SUMTH+TEMP
  172 SUMTH2=SUMTH2+TEMP**2
      NA=NA+3
  180 CONTINUE
      IF((NP-4)*(NP-6))230,190,230
C        ALL FACES FOUND, CALCULATE QUADRATIC ELONGATION FOR TETR AND OC
  190 SUM=0.0
      DO 200 I=1,NP
  200 SUM=SUM+DM(I)**2
      IF(NP.EQ.4) GO TO 210
C        SET UP FOR OCTAHEDRON
      CONS=0.75
      CONS2=90.0
      GO TO 220
C        SET UP FOR TETRAHEDRON
  210 CONS=9.0*SQRT(3.0)/8.0
      CONS2=109.45
  220 VLO=EXP(2.0*ALOG(CONS*VOL)/3.0)
      QE=SUM/(NP*VLO)
      NA=(NA+1)/2
      SUMTH2=0.5*SUMTH2
      SUMTH =0.5*SUMTH
      SIGA=(SUMTH2-2.0*CONS2*SUMTH+NA*CONS2**2)/(NA-1)
      IF(IORG.EQ.0)WRITE(IOUT,6) VOL,QE,SIGA
      GO TO 240
  230 IF(IORG.EQ.0)WRITE(IOUT,8) VOL
  240 CONTINUE
      RETURN
      END
