c
c _____________________________
c
c
	subroutine misquat(qq,thetamin)
	real qq(4,2),thetamin,qresult(4),tmp(2),rquat(4)
	real disor,pi
	real qmax,q1max,q2max
c
       PI=3.14159265
c
c  algorithm for forming resultant quaternion
c  and determining minimum angle taken from Sutton & Baluffi
c
c  note that the resultant quaternion is not returned
c  because it is not in the fundamental zone
c
c  note change of signs to get inverse of second orientation
c
	qresult(1)=qq(1,1)*qq(4,2)-qq(4,1)*qq(1,2)
     &	+qq(2,1)*qq(3,2)-qq(3,1)*qq(2,2)
	qresult(2)=qq(2,1)*qq(4,2)-qq(4,1)*qq(2,2)
     &	+qq(3,1)*qq(1,2)-qq(1,1)*qq(3,2)
	qresult(3)=qq(3,1)*qq(4,2)-qq(4,1)*qq(3,2)
     &	+qq(1,1)*qq(2,2)-qq(2,1)*qq(1,2)
	qresult(4)=qq(4,1)*qq(4,2)+qq(1,1)*qq(1,2)
     &	+qq(2,1)*qq(2,2)+qq(3,1)*qq(3,2)
c
c	write(*,*) 'qresult ',qresult
	qmax=0.
	iqindex=0
	do 10, i=1,4
	qresult(i)=abs(qresult(i))
	if(qresult(i).gt.qmax) then
		qmax=qresult(i)
		iqindex=i
	endif
10	continue
c
	q1max=0.
	iq1index=0
c  find the next highest q component
	do 20, i=1,4
	if(i.eq.iqindex) goto 20
	if(qresult(i).gt.q1max) then
		q1max=qresult(i)
		iq1index=i
	endif
20	continue
c
	disor=amax1(qmax,(qmax+q1max)/sqrt(2.),
     &	(qresult(1)+qresult(2)+qresult(3)+qresult(4))/2.)
	if(disor.gt.1.0) disor=1.0
	if(disor.lt.-1.0) disor=-1.0
	thetamin=acos(disor)*360./pi
c	write(*,*) 'thetamin ',thetamin
c
	q2max=0.
	iq2index=0
c  find the next highest q component
	do 30, i=1,4
	if(i.eq.iqindex.or.i.eq.iq1index) goto 30
	if(qresult(i).gt.q2max) then
		q2max=qresult(i)
		iq2index=i
	endif
30	continue
c
	rquat(4)=qmax
	rquat(3)=q1max
	rquat(2)=q2max
	do 40, i=1,4
		if(i.ne.iqindex.and.
     &	i.ne.iq1index.and.
     &	i.ne.iq2index)  rquat(1)=qresult(i)
40	continue
c  supply the sorted quaternion 
c  CAUTION:  note that q1<q2<q3<q4
c  whereas typical Rodrigues sorting is R1>R2>R3
c
	return
	end
c
c _____________________________
c