      Program HCOOR
C************************************************************************
C*  January 1998                                                        *
C*                           Program HCOOR                              *
C*                                                                      *
C*                         by Mario NARDELLI                            *
C*  Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica,  * della Universita` degli Studi di
C*  Chimica Fisica della Universita' degli Studi di Parma               *
C*  Centro di Studio per la Strutturistica Diffrattometrica del CNR,    *
C*           Viale delle Scienze 78, I-43100 Parma (Italy)              *
C*                                                                      *
C*  This routine calculates the fractional coordinates of atom X lying  *
C*  on the line joining atom A to atom B at a requested distance D from *
C*  atom A and nearest to atom B.                                       *
C*                                                                      *
C*                               INPUT                                  *
C* - Title (80A1)                                                       *
C* - a,b,c,alpha,beta,gamma (free format)                               *
C* - 'ATOM(A)',x/a,y/b,z/c     (free format; ATOM(A):A4)                *
C* - 'ATOM(B)',x/a,y/b,z/c     (free format; ATOM(B):A4)                *
C* - 'ATOM(X)',D               (free format; ATOM(X):A4)                *
C*                                                                      *
C*                               OUTPUT                                 *
C*   Title                                                              *
C*   Table of fractional coordinates                                    *
C*   Distance ATOM(A)...ATOM(X)                                         *
C*                                                                      *
C************************************************************************
      CHARACTER*1 TIT(80)
      CHARACTER*4 ATM(3)
      CHARACTER*14 FILIN,FILOUT
      DOUBLE PRECISION O(3,3),OI(3,3)
      DIMENSION X(3,3),C(3),A(3),XO(3,3),SN(3),CS(3),DDC(3),X1(3),X2(3)
      RD=57.295779513
      WRITE(*,100)
100   FORMAT(1X,'Key in the name of the input file'/)
      READ(*,'(A14)')FILIN
      WRITE(*,101)
101   FORMAT(1X,'Key in the name of the output file'/)
      READ(*,'(A14)')FILOUT
      OPEN(UNIT=13,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=13)
      OPEN(UNIT=15,FILE=FILOUT,FORM='FORMATTED')
      REWIND(UNIT=15)
      READ(13,102)TIT
102   FORMAT(80A1)
      WRITE(15,103)TIT
103   FORMAT(1X,80A1/)
      READ(13,*)C,A
      READ(13,*)(ATM(I),(X(I,K),K=1,3),I=1,2)
      READ(13,*)ATM(3),D13
      DO 1 I=1,3
      AA=A(I)/RD
      SN(I)=SIN(AA)
1     CS(I)=COS(AA)
      CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      SNAR=SQRT(1.0-CSAR**2)
      O(1,1)=C(1)
      O(1,2)=C(2)*CS(3)
      O(1,3)=C(3)*CS(2)
      O(2,1)=0.0
      O(2,2)=C(2)*SN(3)
      O(2,3)=-C(3)*SN(2)*CSAR
      O(3,1)=0.0
      O(3,2)=0.0
      O(3,3)=C(3)*SN(2)*SNAR
      DO 2 I=1,2
      DO 2 J=1,3
      XO(I,J)=0.0
      DO 2 K=1,3
2     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
      DO 3 K=1,3
      DDC(K)=XO(2,K)-XO(1,K)
3     D12=D12+DDC(K)**2
      D12=SQRT(D12)
      DO 4 K=1,3
      DDC(K)=DDC(K)/D12
      X1(K)=XO(1,K)+DDC(K)*D13
4     X2(K)=XO(1,K)-DDC(K)*D13
      DH1=0.0
      DH2=0.0
      DO 5 K=1,3
      DH1=DH1+(XO(2,K)-X1(K))**2
5     DH2=DH2+(XO(2,K)-X2(K))**2
      DH1=SQRT(DH1)
      DH2=SQRT(DH2)
      DO 7 K=1,3
      XO(3,K)=X1(K)
7     IF(DH1.GT.DH2)XO(3,K)=X2(K)
      DH=0.0
      DO 8 K=1,3
8     DH=DH+(XO(1,K)-XO(3,K))**2
      DH=SQRT(DH)
      CALL INVMAT3(O,OI)
      DO 6 J=1,3
      X(3,J)=0.0
      DO 6 K=1,3
6     X(3,J)=X(3,J)+OI(J,K)*XO(3,K)
      WRITE(15,104)(ATM(I),(X(I,K),K=1,3),I=1,3),ATM(1),ATM(3),DH
104   FORMAT(1X,'Atom',9X,'x',9X,'y',9X,'z'/3(1X,A4,3X,3F10.5/)/
     1 ' Distance ',A4,'-',A4,'  = ',F5.3)
      STOP
      END
 
      SUBROUTINE INVMAT3(A,B)
C-----Inverte una matrice 3x3 con il metodo di partizione.
      DOUBLE PRECISION A(3,3),B(3,3),D,AM1,AM2,AM3,AM4,C,X
      D=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      IF(DABS(D).LT.1D-10) GO TO 1
      AM1=A(1,2)*A(2,3)-A(1,3)*A(2,2)
      AM2=A(1,3)*A(2,1)-A(1,1)*A(2,3)
      AM3=A(2,1)*A(3,2)-A(2,2)*A(3,1)
      AM4=A(1,2)*A(3,1)-A(1,1)*A(3,2)
      X=A(3,3)+(AM1*A(3,1)+AM2*A(3,2))/D
      B(3,3)=1D0/X
      X=X*D
      B(1,3)=AM1/X
      B(1,1)=(B(1,3)*AM3+A(2,2))/D
      B(1,2)=(B(1,3)*AM4-A(1,2))/D
      B(2,3)=AM2/X
      B(2,1)=(B(2,3)*AM3-A(2,1))/D
      B(2,2)=(B(2,3)*AM4+A(1,1))/D
      B(3,1)=AM3/X
      B(3,2)=AM4/X
      GO TO 2
1     WRITE(15,100)
100   FORMAT(' The matrix cannot be inverted')
      STOP
2     RETURN
      END
