      program trans                                                
c  calculate orientation angles of reduced cell from old cell and phis,
c  and new cell and transformation matrix
      double precision car(3,3), ar(3,3), c(3,3)                  
      double precision at1(3,3),att1(3,3),at1inv(3,3),t(3,3),det    
      character*80 fname                              
      pi=3.14159                                                        
      cpi=pi/180.0                                                      
      det=0.0
      write(*,100)
100   format(1x,'name of input data file ?',$)
      read(*,110) fname
110   format(a)
      open(unit=8,file=fname,status='old')
c
c read in old cell and orientation
      read(8,*)as,bs,cs,alps,bets,gams
      write(6,*)
      write(6,*)'gn3 cell a*,b*,c*,alpha*,beta*,gamma*'
      write(6,200)as,bs,cs,alps,bets,gams
      alps=alps*cpi
      bets=bets*cpi
      gams=gams*cpi
      call setm(as,bs,cs,alps,bets,gams,ar)
      read(8,*)px,py,pz
      write(6,*)'gn3 orientation PX,PY,PZ'
      write(6,150)px,py,pz
150   format(3f10.3)
      call cwrite(px,py,pz,c)
      call matmult(c,ar,car)
c
c read in cell from tracer 
c
      read(8,*)as,bs,cs,alps,bets,gams
      write(6,*)
      write(6,*)
      write(6,*)'reduced (tracer) cell a*,b*,c*,alpha*,beta*,gamma*'
      write(6,200)as,bs,cs,alps,bets,gams
200   format(1x,6f10.5)
      alps=alps*cpi
      bets=bets*cpi
      gams=gams*cpi
      call setm(as,bs,cs,alps,bets,gams,at1)
CD    write(6,*)' matrix at1..standard orientation of tracer cell'
CD    call matprint(at1)
c
c read in transformation matrix and calc orientation angles
c
      read(8,*)(t(1,i),i=1,3)
      read(8,*)(t(2,i),i=1,3)
      read(8,*)(t(3,i),i=1,3)
      write(6,*)'tracer transformation matrix (cell(1) -> cell(4))'
      call matprint(t)
      call matmult(at1,t,att1)
      call invert3(att1,at1inv,det)
CD    write(6,*)' matrix at1inv....inverse of [at1].[t]'
CD    call matprint(at1inv)
      call matmult(car,at1inv,c)
CD    write(6,*)' matrix c'
CD    call matprint(c)
      call solvps(c,px,py,pz)
      write(6,*)
      write(6,*)'missorientation of reduced cell '
      write(6,900)px,py,pz
900   format(' PX, PY, PZ =',3f12.3) 
      write(6,*)' '
      write(6,*)'with respect to c* along spindle and a*'
      write(6,*)'nearest to x-ray beam in plane defined by'
      write(6,*)'the beam axis and rotation axis'
c
      stop
      end                                 
C                                                                       
C  *************************************************************        
      SUBROUTINE SOLVPS(AM,PX,PY,PZ)                                    
C GIVEN A ROTATION MATRIX AM, EXTRACT THE CORRESPONDING PX,PY,PZ        
      DOUBLE PRECISION AM(3,3),det,aminv(3,3)             
      PI=3.14159                                                        
      CPI=PI/180
      det=0.0 
c
c check determinant if .lt. 0 assume =-1 and multiply all
c  matrix elements by -1.0 to get correct orientation
c
c if not will get 180+px,-py,180+pz
c   
      call invert3(am,aminv,det)
       if (det.lt.0.0) then
        do 10 i=1,3
         do 10 j=1,3
          am(i,j)=-am(i,j)
10     continue
       endif
c
c          
      SY=-AM(3,1)                                                        
      IF (SY.GE.1.000) SY=0.999999                                      
      IF (SY.LE.-1.00) SY=-0.999999                                     
      CY=SQRT(1-SY*SY)                                                  
      IF (CY.LT.0.000001) CY=0.000001                                   
      CZ=AM(1,1)/CY                                                     
      SZ=AM(2,1)/CY                                                     
      SX=AM(3,2)/CY                                                     
      CX=AM(3,3)/CY                                                     
      PX=ATAN2(SX,CX)/CPI                                               
      PY=ATAN2(SY,CY)/CPI                                               
      PZ=ATAN2(SZ,CZ)/CPI                      
C             
      RETURN                                                            
      END                                                               
C                                                                       
C *******************************************************               
      SUBROUTINE SETM(AS,BS,CS,ALPS,BETS,GAMS,AM)                       
C SETS UP ORTHOGONALISATION MATRIX FROM ASTAR TO GAMMASTAR              
C  CSTAR PARALLEL TO Z, ASTAR NEAR x,angles in radians                  
      DOUBLE PRECISION AM(3,3)                                          
      CA=COS(ALPS)                                                      
      SA=SIN(ALPS)                                                      
      CB=COS(BETS)                                                      
      SB=SIN(BETS)                                                      
      CG=COS(GAMS)                                                      
      SG=SIN(GAMS)                                                      
      COSALP=(CB*CG-CA)/(SB*SG)                                         
      SINALP=SQRT(1-COSALP**2)                                          
      COSBET=(CG*CA-CB)/(SG*SA)                                         
      SINBET=SQRT(1-COSBET**2)                                          
      COSGAM=(CA*CB-CG)/(SA*SB)                                         
      SINGAM=SQRT(1-COSGAM**2)                                          
      AM(1,1)=AS*SB                                                     
      AM(2,1)=0                                                         
      AM(3,1)=AS*CB                                                     
      AM(1,2)=-BS*SA*COSGAM                                             
      AM(2,2)=BS*SA*SINGAM                                              
      AM(3,2)=BS*CA                                                     
      AM(1,3)=0                                                         
      AM(2,3)=0                                                         
      AM(3,3)=CS                                                        
      RETURN                                                            
      END                                                               
C                                                                       
C  ***************************************************                  
      SUBROUTINE MATPRINT(AM)                                           
      DOUBLE PRECISION AM(3,3)                                          
      DO 10 I=1,3                                                       
      WRITE(6,900)(AM(I,J),J=1,3)                                       
  900 FORMAT(3F12.6)                                                    
   10 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
C **********************************************************            
      SUBROUTINE MATMULT(A,B,C)                                         
C 3 X 3 MATRICES, AXB INTO C                                            
      DOUBLE PRECISION A(3,3), B(3,3), C(3,3)                           
      DO 20 I=1,3                                                       
      DO 10 J=1,3                                                       
      C(I,J)=A(I,1)*B(1,J) + A(I,2)*B(2,J) + A(I,3)*B(3,J)              
   10 CONTINUE                                                          
   20 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      FUNCTION   ANGLE(     X,Y,Z,I,J,K)                                
C                                                                       
C     USES                 X,Y,Z ARRAYS, WORKS AND PRINTS ANGLE         
C     ATOM I-ATOM J-ATOM K                                              
C      (TAKEN FROM PROGRAM COCALC AND MODIFIED, JUNE 89)                
C                                                                       
      DOUBLE PRECISION X(0:10),Y(0:10),Z(0:10)                          
      PI=3.14159                                                        
C     WRITE(6,921)      I, X(I),Y(I),Z(I)                               
C     WRITE(6,921) J,      X(J),Y(J),Z(J)                               
C     WRITE(6,921) K,      X(K),Y(K),Z(K)                               
c  921 FORMAT(2X,I6,3F10.4)                                              
      DISTIJ = (X(I)-X(J))**2+(Y(I)-Y(J))**2+(Z(I)-Z(J))**2             
      DISTIJ = SQRT(DISTIJ)                                             
      DISTJK = (X(J)-X(K))**2+(Y(J)-Y(K))**2+(Z(J)-Z(K))**2             
      DISTJK = SQRT(DISTJK)                                             
      DISTKI = (X(K)-X(I))**2+(Y(K)-Y(I))**2+(Z(K)-Z(I))**2             
      DISTKI = SQRT(DISTKI)                                             
C     WRITE(6,930) DISTIJ,DISTJK,DISTKI                                 
C 930 FORMAT(3F10.3)                                                    
      COSANG = DISTIJ**2+DISTJK**2-DISTKI**2                            
      COSANG = COSANG/(2*DISTIJ*DISTJK)                                 
      IF (COSANG.LE.-1.0) COSANG=-0.999999                              
      IF (COSANG.GE.1.0)  COSANG=0.9999999                              
      ANG = ACOS(COSANG)                                                
      WRITE(6,900)I,J,K                   ,ANG*180/PI                   
  900 FORMAT(10X,' ANGLE',2X,I2,'-',I2,'-',I2,'IS  ',F6.1)              
      WRITE(6,910) DISTIJ,DISTJK,DISTKI                                 
  910 FORMAT(' DISTANCES',3F10.3)                                       
C                                                                       
      ANGLE=ANG                                                         
      RESULT=ANG                                                        
      RETURN                                                            
      END                                                               
C   ******************************************************              
C                                                                       
      subroutine invert3(a,ainv,det)                                       
c invert a 3x3 matrix ..return inverted matrix and determinant
c 
      double precision a(3,3),ainv(3,3),acof(3,3),det
c
c  calculate cofactors:
c
      acof(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))
      acof(1,2)=-1.0*(a(2,1)*a(3,3)-a(2,3)*a(3,1))
      acof(1,3)=(a(2,1)*a(3,2)-a(2,2)*a(3,1))
      acof(2,1)=-1.0*(a(1,2)*a(3,3)-a(1,3)*a(3,2))
      acof(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))
      acof(2,3)=-1.0*(a(1,1)*a(3,2)-a(1,2)*a(3,1))
      acof(3,1)=(a(1,2)*a(2,3)-a(1,3)*a(2,2))
      acof(3,2)=-1.0*(a(1,1)*a(2,3)-a(1,3)*a(2,1))
      acof(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))
c
c  calculate determinant
c
      det=a(1,1)*acof(1,1)+a(1,2)*acof(1,2)+a(1,3)*acof(1,3)
c
c  transpose and divide by determinant
c
      do 10 i=1,3
       do 10 j=1,3
        ainv(i,j)=acof(j,i)/det
10    continue
c
      return 
      end        
C                                                                                                                                      
c***************************************************************
      subroutine cwrite(px,py,pz,c)
c     writes the missetting matrix [c] from a knowledge of px,py,pz
c
      real px,py,pz,cpx,spx,cpy,spy,cpz,spz
      double precision cx(3,3),cy(3,3),cz(3,3),c(3,3),czy(3,3)
      piby180=3.1415926536/180.0
c
      cpx=cos(px*piby180)
      spx=sin(px*piby180)
      cpy=cos(py*piby180)
      spy=sin(py*piby180)
      cpz=cos(pz*piby180)
      spz=sin(pz*piby180)
c
      cx(1,1)=1.0
      cx(1,2)=0.0
      cx(1,3)=0.0
      cx(2,1)=0.0
      cx(2,2)=cpx
      cx(2,3)=-spx
      cx(3,1)=0.0
      cx(3,2)=spx
      cx(3,3)=cpx
c
      cy(1,1)=cpy
      cy(1,2)=0.0
      cy(1,3)=spy
      cy(2,1)=0.0
      cy(2,2)=1.0
      cy(2,3)=0.0
      cy(3,1)=-spy
      cy(3,2)=0.0
      cy(3,3)=cpy
c
      cz(1,1)=cpz
      cz(1,2)=-spz
      cz(1,3)=0.0
      cz(2,1)=spz
      cz(2,2)=cpz
      cz(2,3)=0.0
      cz(3,1)=0.0
      cz(3,2)=0.0
      cz(3,3)=1.0
c
      call matmult(cz,cy,czy)
      call matmult(czy,cx,c)
c
      return
      end
