C Program to calculate crystal orientation resulting from rotation about
C    lab axes, expressed as total rotation from std position in standard
C    rotx, roty, rotz.
	real*4 RMX(3,3),RMY(3,3),RMZ(3,3),PROD(3,3)
	real*4 TEMP(3,3) ,A,B,X
	character*5 astring
5 	pi=4*ATAN(1.)
c	write(6,*)'pi=',pi
50 	write(6,*) 'standard rotation angles for current position (rx,ry,rz)?'
	read(5,*) rotx, roty,rotz
56 	rotx=rotx*pi/180
	roty=roty*pi/180
	rotz=rotz*pi/180
c Construct rotation matrix for rotation about x (rmx)
58 	rmx(1,1)=1
	rmx(1,2)=0
	rmx(1,3)=0
	rmx(2,1)=0
	rmx(3,1)=0 
	rmx(2,2)=COS(rotx)
	rmx(3,2)=SIN(rotx)
	rmx(2,3)=-rmx(3,2)
	rmx(3,3)=rmx(2,2)
c Construct rotation matrix for rotation about y (rmy)
60 	rmy(2,2)=1
	rmy(1,2)=0
	rmy(3,2)=0
	rmy(2,1)=0
	rmy(2,3)=0
	rmy(1,1)=COS(roty)
	rmy(3,1)=SIN(roty)
	rmy(1,3)=-rmy(3,1)
	rmy(3,3)=rmy(1,1)
c Construct rotation matrix for rotation about z (rmz)
62 	rmz(3,3)=1
	rmz(1,3)=0
	rmz(2,3)=0
	rmz(3,1)=0
	rmz(3,2)=0
	rmz(1,1)=COS(rotz)
	rmz(2,1)=SIN(rotz)
	rmz(1,2)=-rmz(2,1)
	rmz(2,2)=rmz(1,1)
C multiply the matrices together, put result in prod().
c multiply rmz() by rmy() and put in temp()
 	DO 63 i=1,3
	Do 63 K=1,3
c zero the temp matrix:
	temp(i,K)=0
c accumulate product of 
	Do 63 j=1,3
63	temp(i,K)=temp(i,K)+rmy(i,j)*rmz(j,K)

c multiply temp() by rmx() and put in prod()
	do 64 i=1,3
	do 64 K=1,3
	prod(i,K)=0
	do 64 j=1,3
64	prod(i,K)=prod(i,K)+rmx(i,j)*temp(j,K)

C print the matrix:
65 	write(6,*) 'Initial rot matrix '
 	do 66 i=1,3
66	write(6,'(3F12.6)')(prod(i,j),j=1,3)


c	postrot:
67 	write(6,*) 'further rotate about x,y, or z?'
	read(5,*)astring
	if ((astring(1:1).eq.'q').or.(astring(1:1).eq.'Q')) stop
	write(6,*) 'how many degrees?'
	read(5,*)x
	x=x*pi/180
C. make the new rotation matrix in rmx, even if roty or rotz
	if ((astring(1:1).eq.'x').or.(astring(1:1).eq.'X')) then
	 rmx(1,2)=0
	 rmx(1,3)=0
	 rmx(2,1)=0
	 rmx(3,1)=0 
	 rmx(1,1)=1
	 rmx(2,2)=COS(x)
	 rmx(3,2)=SIN(x)
	 rmx(2,3)=-rmx(3,2)
	 rmx(3,3)= rmx(2,2)
	elseif ((astring(1:1).eq.'y').or.(astring(1:1).eq.'Y')) then
	 rmx(1,2)=0
	 rmx(3,2)=0
	 rmx(2,1)=0
	 rmx(2,3)=0
	 rmx(2,2)=1
	 rmx(1,1)=COS(x)
	 rmx(3,1)=SIN(x)
	 rmx(1,3)=-rmx(3,1)
	 rmx(3,3)= rmx(1,1)
	elseif ((astring(1:1).eq.'z').or.(astring(1:1).eq.'Z')) then
	 rmx(1,3)=0
	 rmx(2,3)=0
	 rmx(3,1)=0
	 rmx(3,2)=0
	 rmx(3,3)=1
	 rmx(1,1)=COS(x)
	 rmx(2,1)=SIN(x)
	 rmx(1,2)=-rmx(2,1)
	 rmx(2,2)= rmx(1,1)
	 ELSE 
	 write(6,*) 'I did not understand which angle to rotate about'
	 goto 67
	 endif

C. multiply the matrix describing current position
c  by the new rotation matrix:
	DO 70 i=1,3
	DO 70 K=1,3
	temp(i,K)=0
	DO 70 j=1,3
70	temp(i,K)=temp(i,K)+rmx(i,j)*prod(j,K)

 	DO 80 i=1,3
	DO 80 j=1,3
80	prod(i,j)=temp(i,j)

	write(6,*) 'New rotation matrix:'
 	do 86 i=1,3
86	write(6,'(3F12.6)')(prod(i,j),j=1,3)
c Factor prod() into 3 rotation matrices
 	  b=ASIN( -prod(1,3)) 
c	(b is angle of rotation about y)
  	  IF (b.GT.pi) b=b-2*pi
  	a=ATAN(-prod(2,3)/prod(3,3))
  	  IF (prod(3,3)/COS(b).LT.0.0) a=a+pi
  	  IF (a.gt.pi) a=a-2*pi
  	c=ATAN(-prod(1,2)/prod(1,1))
  	IF (prod(1,1)/COS(b).lt.0.0) c=c+pi
  	IF (c.gt.pi)  c=c-2*pi
91	FORMAT ('rotation about x,y,z',3F9.3)
	write(6,91) 180*a/pi,180*b/pi,180*c/pi

 	GOTO 67
	end