!
! Old Cartesian coordinates to centre of mass Euler conversion
! program. Updated 22/5/11 by DJW to read xyz files from CCD.
!
      PROGRAM XTOE
C
C  Convert H, H, O positions to Eulers - needs H, H, O
C  positions and centres of mass
C
      IMPLICIT REAL*8(A-H,O-Z) 
      CHARACTER(LEN=1) DUMMYS
      PARAMETER (NPQ=63,NNPQ=NPQ*(NPQ-1)/2+1,NZ=0) 
      COMMON /AR/ X(NPQ),Y(NPQ),Z(NPQ), 
     1            A(NPQ),B(NPQ),C(NPQ), 
     2            XX1(NPQ),YY1(NPQ),ZZ1(NPQ), 
     3            XX2(NPQ),YY2(NPQ),ZZ2(NPQ), 
     4            XX3(NPQ),YY3(NPQ),ZZ3(NPQ)
      DATA ANUM/6.0225D+23/ 
      DATA ANGLE,HOLEN,COLEN/104.52D0,0.9572D0,0.15D0/ 
      PI=4.0D0*ATAN(1.0D0) 
      PI2=2.0D0*PI 
      PIH=PI/2.0D0 
      RANGLE=PI*ANGLE/360.0D0 
      OHZ=HOLEN*COS(RANGLE) 
      HYL=HOLEN*SIN(RANGLE) 
      HZL=16.0D0*OHZ/18.0D0 
      OL=-OHZ+HZL 
      READ(*,*) N
      READ(*,*) DUMMYS
      N=N/3
      DO 10 I=1,N
         READ(*,*) DUMMYS,XX3(I),YY3(I),ZZ3(I)
         READ(*,*) DUMMYS,XX1(I),YY1(I),ZZ1(I)
         READ(*,*) DUMMYS,XX2(I),YY2(I),ZZ2(I)
10    CONTINUE
C
      DO 20 I=1,N 
         X(I)=(XX1(I)+XX2(I)+16.0D0*XX3(I))/18.0D0
         Y(I)=(YY1(I)+YY2(I)+16.0D0*YY3(I))/18.0D0
         Z(I)=(ZZ1(I)+ZZ2(I)+16.0D0*ZZ3(I))/18.0D0

         A(I)=DACOS((ZZ3(I)-Z(I))/OL)
         COSA=((ZZ3(I)-Z(I))/OL)
         SINA=DSIN(A(I))

         B(I)=DACOS(-(YY3(I)-Y(I))/(OL*SINA))
         COSB=-(YY3(I)-Y(I))/(OL*SINA)
         SINB=DSIN(B(I))

         C(I)=DACOS((ZZ1(I)-ZZ2(I))/(2.0D0*HYL*SINA))
         COSC=DCOS(C(I))
         SINC=DSIN(C(I))
 
         SP13T=(XX3(I)-X(I))/OL
         SP12T=(XX1(I)-XX2(I))/(2.0D0*HYL)
         SP22T=(YY1(I)-YY2(I))/(2.0D0*HYL)
         SP13=SINA*SINB
         SP12=-(SINC*COSB+COSA*SINB*COSC)
         SP22=-SINC*SINB+COSA*COSB*COSC
         IF (DABS(SP12T/SP12+1.0D0).LT.1.0D-4) THEN
            TEMP=A(I)
            A(I)=PI2-A(I)
            B(I)=PI-B(I)
            C(I)=PI-C(I)
            SINA=-SINA
            COSB=-COSB
            COSC=-COSC
            GOTO 666
         ENDIF
         IF (DABS(SP12T/SP12-1.0D0).GT.1.0D-4) THEN
            IF (DABS(SP22T/SP22-1.0D0).LT.1.0D-4) THEN
               PRINT*,'Inconsitent errors - QUIT'
               STOP
            ENDIF
            IF (DABS(SP13T/SP13+1.0D0).LT.1.0D-4) THEN
               TEMP=B(I) 
               B(I)=PI2-B(I)
               SINB=-SINB
            ELSE
               TEMP=C(I) 
               C(I)=PI2-C(I)
               SINC=-SINC
            ENDIF
         ENDIF
666      WRITE(*,30) X(I),Y(I),Z(I),A(I),B(I),C(I)
30       FORMAT(6F13.6)
20    CONTINUE
      PRINT*,'Rotating coordinates'
      DO 40 J1=1,N
         TEMP=ZZ1(J1)
         ZZ1(J1)=YY1(J1)
         YY1(J1)=XX1(J1)
         XX1(J1)=TEMP
         TEMP=ZZ2(J1)
         ZZ2(J1)=YY2(J1)
         YY2(J1)=XX2(J1)
         XX2(J1)=TEMP
         TEMP=ZZ3(J1)
         ZZ3(J1)=YY3(J1)
         YY3(J1)=XX3(J1)
         XX3(J1)=TEMP
40    CONTINUE
C
      DO 50 I=1,N 
         X(I)=(XX1(I)+XX2(I)+16.0D0*XX3(I))/18.0D0
         Y(I)=(YY1(I)+YY2(I)+16.0D0*YY3(I))/18.0D0
         Z(I)=(ZZ1(I)+ZZ2(I)+16.0D0*ZZ3(I))/18.0D0

         A(I)=DACOS((ZZ3(I)-Z(I))/OL)
         COSA=((ZZ3(I)-Z(I))/OL)
         SINA=DSIN(A(I))

         B(I)=DACOS(-(YY3(I)-Y(I))/(OL*SINA))
         COSB=-(YY3(I)-Y(I))/(OL*SINA)
         SINB=DSIN(B(I))

         C(I)=DACOS((ZZ1(I)-ZZ2(I))/(2.0D0*HYL*SINA))
         COSC=DCOS(C(I))
         SINC=DSIN(C(I))
 
         SP13T=(XX3(I)-X(I))/OL
         SP12T=(XX1(I)-XX2(I))/(2.0D0*HYL)
         SP22T=(YY1(I)-YY2(I))/(2.0D0*HYL)
         SP13=SINA*SINB
         SP12=-(SINC*COSB+COSA*SINB*COSC)
         SP22=-SINC*SINB+COSA*COSB*COSC
         IF (DABS(SP12T/SP12+1.0D0).LT.1.0D-4) THEN
            TEMP=A(I)
            A(I)=PI2-A(I)
            B(I)=PI-B(I)
            C(I)=PI-C(I)
            SINA=-SINA
            COSB=-COSB
            COSC=-COSC
            GOTO 667
         ENDIF
         IF (DABS(SP12T/SP12-1.0D0).GT.1.0D-4) THEN
            IF (DABS(SP22T/SP22-1.0D0).LT.1.0D-4) THEN
               PRINT*,'Inconsitent errors - QUIT'
               STOP
            ENDIF
            IF (DABS(SP13T/SP13+1.0D0).LT.1.0D-4) THEN
               TEMP=B(I) 
               B(I)=PI2-B(I)
               SINB=-SINB
            ELSE
               TEMP=C(I) 
               C(I)=PI2-C(I)
               SINC=-SINC
            ENDIF
         ENDIF
667      WRITE(*,30) X(I),Y(I),Z(I),A(I),B(I),C(I)
50    CONTINUE
      PRINT *,'Writing centre of mass/Euler angles in OPTIM input order'
      DO J1=1,N
         PRINT '(3G20.10)',X(J1),Y(J1),Z(J1)
      ENDDO
      DO J1=1,N
         PRINT '(3G20.10)',A(J1),B(J1),C(J1)
      ENDDO
      STOP
      END
