c
c  Catalogue of Subroutines in TEXTUR.F
c
c	TEXTUR - sets up required arrays
c	CHARLOW - creates properties for low angle bdies
c	QUAT  convert euler angles to a quaternion
c	MISQUAT  combines two quaternions and finds the minimum angle
c	CSLQUAT  gets distance (angle) from a CSL boundary: misorient must be in 
c		the fundamental zone,  q1<q2<q3<q4
c	Q2MAT  converts quaternion to a 3x3 matrix
c	QRVEC  uses a quaternion to rotate a vector
c	COMQUAT  combines two quaternions (no disorientation)
c	PRESYMM  applies symmetry element from QUATSYMM as 1st q
c	POSTSYMM applies symm elem as 2nd q
c	DISQUAT  finds the disorientation the hard way but also determines
c		which symmetry operators put you in the fundamental zone
c	DISQUAT2  finds the disorientation, using two quaterions as input
c       the hard way but also determines
c		which symmetry operators put you in the fundamental zone
c
	subroutine textur
	include 'common.f'
	common/a1/ quatsymm(4,48),numsymm
	common/a2/quatcsl(4,50),numcsl,nsig(50)
c  following code taken from LApp68.f for reading in TEXIN
c
	real evm,fk1(9),nstate
	character*78 textitl
	character eulan*25,nomen,iper*6
	real ph,th,ps,dph,dth,dps,sph,cph,cth,sth,sps,cps,qr(4)
	real theta,r1,r2,r3
	integer lmn
c
	rad=57.29578
	write(*,*) 'reading symmetry operators from "quat.symm.cubic"'
	write(*,*) 'reading texture for unrecrystallized grains from "texin1"'
	write(*,*) 'reading texture for recrystallized grains from "texin2"'
	open(unit=3,name='quat.symm.cubic',status='old')
	read(3,*)	! skip over title line
	read(3,*) numsymm	! read number of symmetry operators
	if(numsymm.gt.24) then
		numsymm=24
		write(*,*) 'too many symmetry operators, cutting back to 24!'
		else
		write(*,*) 'reading ',numsymm,' symmetry operators'
		endif
c
	do 1800, i=1,numsymm
		read(3,*) (quatsymm(ijk,i),ijk=1,4)
c		write(*,"('symm. no. ',i3,4f8.3)") i,(quatsymm(ijk,i),ijk=1,4)
1800	continue
	close(unit=3)
c
	write(*,*) 'reading CSL information from "quat.csl.data"'
	open(unit=3,name='quat.csl.data',status='old')
	read(3,*)	! skip over title line
	read(3,*) numcsl	! read number of CSL types
	if(numcsl.gt.50) then
		numcsl=50
		write(*,*) 'too many CSL types, cutting back to 50!'
		else
		write(*,*) 'reading ',numcsl,' boundary types'
		endif
	read(3,*)	! skip over comment
	do 1900, i=1,numcsl
		read(3,*) nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4)
c		write(*,*) i,nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4)
1900	continue
	close(unit=3)
c
c	write(*,*) 'reading grain orientations from TEXIN'
	do 529 ng=1,q
	if(ng.eq.1.or.ng.eq.(q2+1)) then
c  adr, 22 vii 00; actually need independent orientations for all Q
		if(ng.eq.1) open(unit=3,name='texin1',status='old')
		if(ng.eq.(q2+1)) open(unit=3,name='texin2',status='old')
		read(3,7) textitl
    7 	format(a)
		read(3,*)
		read(3,70) evm,fk1,nstate
   70 	format(10f7.3,i3)
 		read(3,71) eulan,iper
   71 	format(a25,49x,a6)
 		nomen=eulan
 		print 76,textitl
   76 	format(/
     #	' TEXIN file= '/1x,a/)
	endif
c
c	do 529 ng=1,q2
c      do 520 i=1,18
c 520  state(i,ng)=0.
c      read(3,50,end=530) d1,d2,d3,grwt(ng),(state(i,ng),i=1,nstate)
c      read(3,50,end=530) d1,d2,d3
      read(3,*,end=530) d1,d2,d3,dummy,tayf(ng)
      if(tayf(ng).eq.0.) tayf(ng)=1.
c  now let's read in the taylor factor for each orientation (skip wt.)
c
   50 format(6f8.2,24f6.2)
c      read(3,*,end=530) d1,d2,d3,grwt(ng)
c            .: for old unformatted TEXIN files
c
   20 dth=d2
      dps=d1
c      if(nomen.eq.'K') then
c        dph=d3
c      elseif(nomen.eq.'R') then
c        dph=180.-d3
c      elseif(nomen.eq.'B') then
c        dps=d1-90.
c        dph=90.-d3
c      else
c        dth=d1
c        dph=d2-90.
c        dps=90.-d3
c	endif
c	ph=dph/rad
c	th=dth/rad
c	ps=dps/rad
	if(nomen.eq.'B') then 
		call quatB(d1/rad,d2/rad,d3/rad,qr)
	elseif(nomen.eq.'K') then
		call  quatB((d1+90.)/rad,d2/rad,(90.-d3)/rad,qr)
	else
		write(*,*) 'cannot deal with type ',nomen,' !!!'
		stop
	endif
c  use the new conversion routine based on Altmann if Bunge
c  else, use longer routine
c
c      if(grwt(ng).eq.0.) grwt(ng)=1.
c      volume=volume+grwt(ng)
c      do 524 i=1,18
c  524 if(state(i,ng).eq.0.) state(i,ng)=1.
c      if(evm.ge.0.01) goto 526
c        state(1,ng)=harpar(1)
c        state(2,ng)=harpar(3)
c
c         *** define initial crss for each mode in each grain:
c  526 if(namodes.gt.1) then
c         do 7770 im=1,namodes
c         state(2*im+1,ng)=taux(im,1)
c 7770    state(2*im+2,ng)=taux(im,2)
c      endif
c     ///////////
c      call euler(2,nomen,d1,d2,d3,ng,0)
c     ///////////
c
	do 527 j=1,4
c	aquat(j,ng+q2)=qr(j)
  527 aquat(j,ng)=qr(j)
c  NO!  Only once through!
c  copy result into the quaternion array for both recrystallized
c  and unrecrystallized grains
c
	write(*,"('orient., quat ',i4,5(x,f8.3))") ng,qr,tayf(ng)
  529 continue
530  	if(ng.lt.q) stop 'INSUFFICIENT GRAINS IN TEXIN!'
c530  	if(ng.lt.q2) stop 'INSUFFICIENT GRAINS IN TEXIN!'
	return
	end
c
c ________________________________________
c
c
	subroutine textur2
c
c  second version, not for use in REX2D
c
	common/a1/ quatsymm(4,48),numsymm
	common/a2/quatcsl(4,50),numcsl,nsig(50)
	common/a3/aquat(4,100000),ngrain
     &	,phi1(100000),capphi(100000),phi2(100000)
	integer numsymm,numcsl,nsig,ngrain
	real phi1,capphi,phi2
c
c	include 'common.f'
c  following code taken from LApp68.f for reading in TEXIN
c
	real evm,fk1(9),nstate
	character*78 textitl
	character eulan*25,nomen,iper*6
	real ph,th,ps,dph,dth,dps,sph,cph,cth,sth,sps,cps,qr(4)
	real theta,r1,r2,r3
	integer lmn
c
	rad=57.29578
	write(*,*) 'reading symmetry operators from "quat.symm.cubic"'
	open(unit=3,name='quat.symm.cubic',status='old')
	read(3,*)	! skip over title line
	read(3,*) numsymm	! read number of symmetry operators
	if(numcsl.gt.24) then
		numcsl=24
		write(*,*) 'too many symmetry operators, cutting back to 24!'
		else
		write(*,*) 'reading ',numsymm,' symmetry operators'
		endif
c
	do 1800, i=1,numsymm
		read(3,*) (quatsymm(ijk,i),ijk=1,4)
c		write(*,*) 'symm. no. ',i,(quatsymm(ijk,i),ijk=1,4)
1800	continue
	close(unit=3)
c
	write(*,*) 'reading CSL information from "quat.csl.data"'
	open(unit=3,name='quat.csl.data',status='old')
	read(3,*)	! skip over title line
	read(3,*) numcsl	! read number of CSL types
	if(numcsl.gt.50) then
		numcsl=50
		write(*,*) 'too many CSL types, cutting back to 50!'
		else
		write(*,*) 'reading ',numcsl,' boundary types'
		endif
	read(3,*)	! skip over comment
	do 1900, i=1,numcsl
		read(3,*) nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4)
c		write(*,*) i,nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4)
1900	continue
	close(unit=3)
c
	write(*,*) 'reading grain orientations from TEXIN'
	open(unit=3,name='texin',status='old')
	read(3,7) textitl
    7 format(a)
	read(3,*)
	read(3,70) evm,fk1,nstate
   70 format(10f7.3,i3)
 	read(3,71) eulan,iper
   71 format(a25,49x,a6)
 	nomen=eulan
 	print 76,textitl
   76 format(/
     #' TEXIN file= '/1x,a/)
c
      do 529 ng=1,100000
c      do 520 i=1,18
c 520  state(i,ng)=0.
c      read(3,50,err=530,end=530) d1,d2,d3,grwt(ng),(state(i,ng),i=1,nstate)
      read(3,50,end=530,err=530) d1,d2,d3
c	write(*,*) 'angles..  ',d1,d2,d3
c	if((d1.eq.999.).and.(d2.eq.999.).and.(d3.eq.999.)) goto 530
   50 format(6f8.2,24f6.2)
c      read(3,*,end=530) d1,d2,d3,grwt(ng)
c            .: for old unformatted TEXIN files
c
	if(nomen.eq.'B') then 
	   phi1(ng)=d1
	   capphi(ng)=d2
	   phi2(ng)=d3
	elseif(nomen.eq.'K') then
	   phi1(ng)=d1+90.
	   capphi(ng)=d2
	   phi2(ng)=90.-d3
	endif
   20 dth=d2
      dps=d1
	if(nomen.eq.'B') then 
		call quatB(d1/rad,d2/rad,d3/rad,qr)
	elseif(nomen.eq.'K') then
		call  quatB((d1+90.)/rad,d2/rad,(90.-d3)/rad,qr)
	else
		stop 'cannot deal with that type! '
	endif
c  use the new conversion routine based on Altmann if Bunge
c  else, use longer routine
c
c
c      if(grwt(ng).eq.0.) grwt(ng)=1.
c      volume=volume+grwt(ng)
c      do 524 i=1,18
c  524 if(state(i,ng).eq.0.) state(i,ng)=1.
c      if(evm.ge.0.01) goto 526
c        state(1,ng)=harpar(1)
c        state(2,ng)=harpar(3)
c
c         *** define initial crss for each mode in each grain:
c  526 if(namodes.gt.1) then
c         do 7770 im=1,namodes
c         state(2*im+1,ng)=taux(im,1)
c 7770    state(2*im+2,ng)=taux(im,2)
c      endif
c     ///////////
c      call euler(2,nomen,d1,d2,d3,ng,0)
c     ///////////
c
	do 527 j=1,4
  527 aquat(j,ng)=qr(j)
c
  529 continue
  530 	continue
	ngrain=ng-1
c
  	return
	end
c
c ________________________________________
c
	subroutine tex_out(one,ten,sum)
c
c  writes out a "*.wts" file
c
c	ADR,  viii 00
c
	include 'common.f'
	character header*80
	character*1 one,ten
c  following code taken from LApp68.f for reading in TEXIN
c
	real evm,fk1(9),nstate
	character*78 textitl
	character eulan*25,nomen,iper*6
	real ph,th,ps,dph,dth,dps
	real theta,r1,r2,r3,d1,d2,d3
	integer lmn
	real sum,tmp,angle
	real qr(4)
c  QR  contains local quaternion for conversion to Euler angles
	real pquat(4,2)
c
	pquat(1,2)=0.
	pquat(2,2)=0.
	pquat(3,2)=0.
	pquat(4,2)=1.0
c  set up a quaternion for the cube orientation
c
	open(17,file=keyword//'.'//ten//one//'.wts',recl=132,status='unknown')
c
	sin60=0.8660254
	pi=3.141592654
c
	sum=0.0
c  SUM is used to estimate nearness to an orientation
c  1 iv 01: for now, just use the CUBE, so can use the stored quat.
c
c	write(*,*) 'reading grain orientations from TEXIN'
	do 529 ng=1,q
	if(ng.eq.1.or.ng.eq.(q2+1)) then
c  adr, 22 vii 00; actually need independent orientations for all Q
		if(ng.eq.1) open(unit=3,name='texin1',status='old')
		if(ng.eq.(q2+1)) then
			close(3)
			open(unit=3,name='texin2',status='old')
			endif
		read(3,7) textitl
    7 	format(a)
		if(ng.eq.1) write(17,7) textitl
		read(3,*)
	header='Evm    F11    F12    F13    F21    '
     &//'F22    F23    F31    F32    F33 nstate'
		if(ng.eq.1) write(17,"(a)") header
		read(3,70) evm,fk1,nstate
   70 	format(10f7.3,i3)
		if(ng.eq.1) write(17,70) evm,fk1,nstate
		read(3,71) eulan,iper
   71 	format(a25,49x,a6)
		if(ng.eq.1) write(17,71) eulan,iper
 		nomen=eulan
c 		print 76,textitl
   76 	format(/
     #	' TEXIN file= '/1x,a/)
	endif
c
      read(3,*,err=530) d1,d2,d3
c
c  as of 9 iv 01, no longer use orientations in the file
c  but convert the quaternions back to Euler angles for
c  storage with q2eulB;  adr
c
		do 527 j=1,4
			pquat(j,1)=aquat(j,ng)
  527		qr(j)=aquat(j,ng)
c  copy quat into QR and into PQUAT (see below)
		call q2eulB(d1,d2,d3,qr)
c  convert to Bunge Euler angles
		d1=d1*180./pi
		d2=d2*180./pi
		d3=d3*180./pi
		if(areatex(ng).ne.0) then
		write(17,50) d1,d2,d3,10.*float(areatex(ng))/float(mn),tayf(ng)
		endif
c  and write the orientation out to file if there is a non-zero area
c
   50 format(3f8.2,f8.4,2f8.2,24f6.2)
c
c
c      if(grwt(ng).eq.0.) grwt(ng)=1.
c      volume=volume+grwt(ng)
c      endif
c     ///////////
c
c
		call misquat(pquat,angle)
c  use MISQUAT to find out how near the orientation is to cube
c  or whatever we use for the orientation of interest
c  CAUTION: MISQUAT is for misorientations; better to use the
c  standard symmetry operations for general orientations!!!!!!
c
		tmp=exp(-1.*((0.05*angle)**2))*float(areatex(ng))/float(mn)
		sum=sum+tmp
c
  529 continue
c
	close(17)
	write(*,*) 'TEX_OUT: weighted sum, grains near cube = ',sum
	return
c
530	stop 'error in reading from texin*'
	end
c
c ________________________________________
c
	subroutine vecnorm(vec)
	real vec(3),rnorm
	rnorm=0.
	do 10, i=1,3
	rnorm=rnorm+vec(i)**2
10	continue
	if(rnorm.le.0.0) return
	rnorm=sqrt(rnorm)
	do 20, i=1,3
	vec(i)=vec(i)/rnorm
20	continue
	return
	end
c
c ________________________________________
c
	subroutine charlow(axis,aenerg,amobility,numtypes)
c
	integer numtypes
	real axis(20,3),aenerg(20),amobility(20),rnorm,remax,rmmax
c  AXIS describes the rotation axes, AENERG the energy
c  AMOBILITY the mobility
c
	include 'common.f'
c
	numtypes=13
c
	do 10, i=1,20
		aenerg(i)=0.
		amobility(i)=0.
		do 10, j=1,3
			axis(i,j)=0.
10			continue
c
	axis(1,1)=0.
	axis(1,2)=0.
	axis(1,3)=1.
	aenerg(1)=0.33
	amobility(1)=0.01
c
	axis(2,1)=0.
	axis(2,2)=1.
	axis(2,3)=5.
	aenerg(2)=0.3
	amobility(2)=0.02
c
	axis(3,1)=0.
	axis(3,2)=2.
	axis(3,3)=5.
	aenerg(3)=0.3
	amobility(3)=0.1
c
	axis(4,1)=0.
	axis(4,2)=2.
	axis(4,3)=3.
	aenerg(4)=0.26
	amobility(4)=0.02
c
	axis(5,1)=0.
	axis(5,2)=1.
	axis(5,3)=1.
	aenerg(5)=0.3
	amobility(5)=0.17
c
	axis(6,1)=2.
	axis(6,2)=7.
	axis(6,3)=7.
	aenerg(6)=0.26
	amobility(6)=0.1
c
	axis(7,1)=2.
	axis(7,2)=3.
	axis(7,3)=3.
	aenerg(7)=0.23
	amobility(7)=0.05
c
	axis(8,1)=1.
	axis(8,2)=1.
	axis(8,3)=1.
	aenerg(8)=0.23
	amobility(8)=0.95
c
	axis(9,1)=3.
	axis(9,2)=3.
	axis(9,3)=5.
	aenerg(9)=0.23
	amobility(9)=0.02
c
	axis(10,1)=1.
	axis(10,2)=1.
	axis(10,3)=3.
	aenerg(10)=0.23
	amobility(10)=0.02
c
	axis(11,1)=1.
	axis(11,2)=1.
	axis(11,3)=7.
	aenerg(11)=0.33
	amobility(11)=0.05
c
	axis(12,1)=1.
	axis(12,2)=2.
	axis(12,3)=5.
	aenerg(12)=0.3
	amobility(12)=0.1
c
	axis(13,1)=4.
	axis(13,2)=8.
	axis(13,3)=11.
	aenerg(13)=0.26
	amobility(13)=0.17
c
	do 100, i=1,numtypes
		rnorm=0.
		do 90, j=1,3
			rnorm=rnorm+axis(i,j)**2
90		continue
		rnorm=sqrt(rnorm)
		do 95, j=1,3
			axis(i,j)=axis(i,j)/rnorm
95		continue
100	continue
c  make unit vectors
c
		remax=aenerg(1)
		rmmax=amobility(1)
	do 200, i=2,numtypes
		if(aenerg(i).gt.remax) remax=aenerg(i)
		if(amobility(i).gt.rmmax) rmmax=amobility(i)
200		continue
	do 210, i=1,13
		aenerg(i)=aenerg(i)/remax
		amobility(i)=amobility(i)/rmmax
210		continue
c  normalize the energies and mobilities
c
	return
	end
c
c
c ________________________________________
c
	subroutine quatB(p1,p,p2,q)
c  convert Bunge Euler angles (radians) to quaternion; 
c  direct conversion from angles, see Altmann's book
c  the rotation is a vector transformation (active rotation)
c
	real p1,p,p2,q(4)
	double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm
c
c
       PI=3.14159265
c  form cosine, sine of Phi, and sum & diff of phi1, phi2
       S=DSIN(0.5d0*P)
       C=DCOS(0.5d0*P)
       S1=DSIN(0.5d0*(P1-P2))
       C1=DCOS(0.5d0*(P1-P2))
       S2=DSIN(0.5d0*(P1+P2))
       C2=Dcos(0.5d0*(P1+P2))
c      write(*,*) 's1,c1,s,c,s2,c2'
c      write(*,*) s1,c1
c      write(*,*) s,c
c      write(*,*) s2,c2
	q(1)=s*c1
	q(2)=s*s1
	q(3)=c*s2
	q(4)=c*c2
	return
	end
c
c ________________________________________
c
	subroutine quat(p1,p,p2,q)
c  convert Euler  angles to quaternion; note that this is derived from LApp
c  so the rotation is a vector transformation (active rotation)
c
	real p1,p,p2,q(4)
	double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm
c
c
       PI=3.14159265
c
c first we get the matrix
       C1=DCOS(1d0*P1)
       C=DCOS(1d0*P)
       C2=DCOS(1d0*P2)
       S1=DSIN(1d0*P1)
       S=DSIN(1d0*P)
       S2=DSIN(1d0*P2)
c      write(*,*) 's1,c1,s,c,s2,c2'
c      write(*,*) s1,c1
c      write(*,*) s,c
c      write(*,*) s2,c2
      d(1,1)=-s1*s2-c2*c1*c
      d(1,2)=c1*s2-c2*s1*c
      d(1,3)=c2*s
      d(2,1)=s1*c2-s2*c1*c
      d(2,2)=-c2*c1-s2*s1*c
      d(2,3)=s*s2
      d(3,1)=s*c1
      d(3,2)=s1*s
      d(3,3)=c
c  try the transpose!
c
c      d(1,1)=-s1*s2-c2*c1*c
c      d(2,1)=c1*s2-c2*s1*c
c      d(3,1)=c2*s
c      d(1,2)=s1*c2-s2*c1*c
c      d(2,2)=-c2*c1-s2*s1*c
c      d(3,2)=s*s2
c      d(1,3)=s*c1
c      d(2,3)=s1*s
c      d(3,3)=c
c	do 10, i=1,3
c	write(*,*) 'matrix= ',(d(i,j),j=1,3)
c10	continue
c  taken from LApp68:euler
c
c       D(1,1)=C1*C2-S1*S2*C
c       D(1,2)=S1*C2+C1*S2*C
c       D(1,3)=S2*S
c       D(2,1)=-C1*S2-S1*C2*C
c       D(2,2)=-S1*S2+C1*C2*C
c       D(2,3)=C2*S
c       D(3,1)=S1*S
c       D(3,2)=-C1*S
c       D(3,3)=C
c
	rtmp=((d(1,1)+d(2,2)+d(3,3))-1.d0)/2.d0
	if(rtmp.lt.-1.d0) rtmp=-1.d0
	if(rtmp.gt.1.d0) rtmp=1.d0
	tmp(4)=dacos(rtmp)
c	write(*,*) 'angle, theta = ',tmp(4)*180./pi
	sin2=dsin(tmp(4)/2.d0)
	cos2=dcos(tmp(4)/2.d0)
c  magnitude of rotation
	tmp(1)=(d(2,3)-d(3,2))
	tmp(2)=(d(3,1)-d(1,3))
	tmp(3)=(d(1,2)-d(2,1))
	rnorm=dsqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
	if(rnorm.lt.1.d-5) then
		tmp(1)=dsqrt(1.d0+d(1,1))
		tmp(2)=dsqrt(1.d0+d(2,2))
		tmp(3)=dsqrt(1.d0+d(3,3))
		rnorm=dsqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
	endif
c	write(*,*) rnorm, tmp(1),tmp(2),tmp(3),tan(tmp(4)/2.)
c
c  now to form the quaternion
	q(1)=sin2*tmp(1)/rnorm
	q(2)=sin2*tmp(2)/rnorm
	q(3)=sin2*tmp(3)/rnorm
	q(4)=cos2
c
c	write(*,*) 'quaternion= ',q
	return
	end
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
c
	subroutine misquatinv(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  "inv" indicates reverse order of quaternions
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)=-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)=-1.*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)=-1.*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
	subroutine cslquat(qq,numsig,theta,qresult)
	real qq(4),theta,qresult(4),tmp(2)
	real pi
	integer numsig	! index of CSL boundary - must be supplied
	include 'common.f'
	common/a2/quatcsl(4,50),numcsl,nsig(50)
c
       PI=3.14159265
c
c  algorithm for forming resultant quaternion
c  and determining angle based on fourth component
c
c  note change of signs to get inverse of second orientation
c
	qresult(1)=qq(1)*quatcsl(4,numsig)-qq(4)*quatcsl(1,numsig)
     &	+qq(2)*quatcsl(3,numsig)-qq(3)*quatcsl(2,numsig)
	qresult(2)=qq(2)*quatcsl(4,numsig)-qq(4)*quatcsl(2,numsig)
     &	+qq(3)*quatcsl(1,numsig)-qq(1)*quatcsl(3,numsig)
	qresult(3)=qq(3)*quatcsl(4,numsig)-qq(4)*quatcsl(3,numsig)
     &	+qq(1)*quatcsl(2,numsig)-qq(2)*quatcsl(1,numsig)
	qresult(4)=qq(4)*quatcsl(4,numsig)+qq(1)*quatcsl(1,numsig)
     &	+qq(2)*quatcsl(2,numsig)+qq(3)*quatcsl(3,numsig)
c
c	write(*,*) 'qresult ',qresult
	if(qresult(4).gt.1.0) qresult(4)=1.0
	if(qresult(4).lt.-1.0) qresult(4)=-1.0
	theta=acos(qresult(4))*360./pi
c	write(*,*) 'theta= ',theta
c
	return
	end
c
c  ______________________
c
	subroutine q2eulB(p1,p,p2,qq)
c
c  converts quaternion to Bunge Euler angles
c  based on Altmann's solution for Euler->quat
c
	real qq(4),p1,p,p2
	real sum,diff
c
       PI=3.14159265
c
	if(abs(qq(2).lt.1e-35).and.abs(qq(1).lt.1e-35)) then
		diff=pi/4.
	else
		diff=atan2(qq(2),qq(1))
	endif
c  ATAN2 is unhappy is both arguments are exactly zero!
c
	if(abs(qq(3).lt.1e-35).and.abs(qq(4).lt.1e-35)) then
		sum=pi/4.
	else
		sum=atan2(qq(3),qq(4))
	endif
c
	p1=(diff+sum)
	p2=(sum-diff)
	tmp=sqrt(qq(3)**2+qq(4)**2)
	if(tmp.gt.1.0) tmp=1.0
	p=2.*acos(tmp)
c	write(*,*) ' quaternion input= ',qq
c	write(*,*) 'Bunge angles output= ',p1,p,p2
c	write(*,*) ' angles [degrees]= ',180.*p1/pi,180.*p/pi,180.*p2/pi
	return
	end
c
c  ______________________
c
	subroutine q2mat(qq,dd)
c
c  converts quaternion to matrix
c
	real qq(4),dd(3,3)
	real t1
c
	t1=qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
	do 100, i=1,3
	do 90,  j=1,3
	dd(i,j)=0.
	if(i.eq.j) then
		dd(i,j)=dd(i,j)+t1
		endif
	dd(i,j)=dd(i,j)+(2.*qq(i)*qq(j))
	if(i.ne.j) then
		do 80, ijk=1,3
			if(i.ne.ijk.and.j.ne.ijk) k=ijk
80			continue
		kl=j-i
		if(kl.eq.2) kl=-1
		if(kl.eq.-2) kl=1
c  poor man"s permutation tensor!
		dd(i,j)=dd(i,j)-(2.*qq(k)*qq(4)*float(kl))
		endif
90		continue
100	continue
	return
	end
c
c +++++++++++++++++++
c
	subroutine vsort(a)
	real a(3),rmax,rmin
	integer i,imax,imin
c	write(*,*) 're-sorting components of vector'
	a(1)=abs(a(1))
	a(2)=abs(a(2))
	a(3)=abs(a(3))
	rmax=0.0
	rmin=1.0
	do 10, i=1,3
		if(a(i).ge.rmax) then
			rmax=a(i)
			imax=i
			endif
		if(a(i).le.rmin) then
			rmin=a(i)
			imin=i
			endif
10	continue
c
	do 15, i=1,3
		if(i.ne.imax.and.i.ne.imin)	a(2)=a(i)
15	continue
	a(1)=rmax
	a(3)=rmin
c	write(*,*) 'sorted a= ',a
c  above code sorts a into descending order
c
	return
	end
c
c
c  ______________________
c
	subroutine rfrvec(rf,din,dout)
c
c  uses RF-vector to rotate a vector from DIN to DOUT; 
c  Eq. 1.36 in Sutton & Balluffi
c
	real rf(3),din(3),dout(3)
	real tmp,tmp2,cross(3)
c
	tmp=2.*scalar(rf,din)
	tmp2=rf(1)**2+rf(2)**2+rf(3)**2
	call vecpro2(din,rf,cross)
	dout(1)=((tmp*rf(1))+(din(1)*(1.-tmp2)))-(2.*cross(1))/(1.+tmp2)
	dout(2)=((tmp*rf(2))+(din(2)*(1.-tmp2)))-(2.*cross(2))/(1.+tmp2)
	dout(3)=((tmp*rf(3))+(din(3)*(1.-tmp2)))-(2.*cross(3))/(1.+tmp2)
	return
	end
c
c  ______________________
c
	subroutine qrvec(qq,din,dout)
c
c  uses quaternion to rotate a vector from DIN to DOUT; 
c  note similarity to above routine
c  note that -q gives the same result as it should (represents same rotation)
c
	real qq(4),din(3),dout(3)
	real t1
c
	t1=qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
	do 100, i=1,3
	dout(i)=t1*din(i)
c
		do 90, j=1,3
			dout(i)=dout(i)+(2.*qq(i)*qq(j)*din(j))
c
			if(i.ne.j) then
				do 80, ijk=1,3
					if(i.ne.ijk.and.j.ne.ijk) k=ijk
80				continue
				kl=j-i
				if(kl.eq.2) kl=-1
				if(kl.eq.-2) kl=1
c  poor man"s permutation tensor!
				dout(i)=dout(i)-(2.*qq(k)*qq(4)*float(kl)*din(j))
		endif
90		continue
100	continue
	return
	end
c
c
c _____________________________
c
c
	subroutine comquat2(qq1,qq2,qresult)
	real qq1(4),qq2(4),qresult(4)
c
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion as Q1 x Q2^-1
c  where Q2 follows Q1, i.e. is the second rotation;
c  note change of signs to get inverse of second orientation
c
	qresult(1)=qq1(1)*qq2(4)-qq1(4)*qq2(1)
     &	+qq1(2)*qq2(3)-qq1(3)*qq2(2)
	qresult(2)=qq1(2)*qq2(4)-qq1(4)*qq2(2)
     &	+qq1(3)*qq2(1)-qq1(1)*qq2(3)
	qresult(3)=qq1(3)*qq2(4)-qq1(4)*qq2(3)
     &	+qq1(1)*qq2(2)-qq1(2)*qq2(1)
	qresult(4)=qq1(4)*qq2(4)+qq1(1)*qq2(1)
     &	+qq1(2)*qq2(2)+qq1(3)*qq2(3)
c
c	write(*,*) 'qresult ',qresult
c
	return
	end
c
c
c _____________________________
c
c
	subroutine comquat(qq,qresult)
	real qq(4,2),qresult(4)
c
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion as Q1 x Q2^-1
c  where Q2 follows Q1, i.e. is the second rotation;
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
c
	return
	end
c
c _____________________________
c
c
	subroutine comquatinv(qq,qresult)
	real qq(4,2),qresult(4)
c
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion as Q1^1 x Q2
c  where Q2 follows Q1, i.e. is the second rotation;
c  note change of signs to get inverse of second orientation
c
	qresult(1)=(-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)=(-1.*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)=(-1.*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
c
	return
	end
c
c _____________________________
c
c
	subroutine presymm(qq,lindex,qresult)
	real qq(4),qresult(4)
	integer lindex
c
	include 'common.f'
	common/a1/ quatsymm(4,48),numsymm
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion
c  from applying a symmetry operator, QUATSYMM(n,lindex)
c  to the first quaternion/rotation
c
c  If the symm oper == O, then
c  Qresult = (O x QQ) where QQ= Q1 x Q2^-1
c  i.e. Qresult = (O x Q1) x Q2^-1
c  note change of signs to get inverse of first orientation
c
	if(lindex.gt.numsymm) stop 'error in presymm, lindex>numsymm'
	if(lindex.lt.1) stop 'error in presymm, lindex<1'
	qresult(1)=quatsymm(4,lindex)*qq(1)+quatsymm(1,lindex)*qq(4)
     &	+quatsymm(3,lindex)*qq(2)-quatsymm(2,lindex)*qq(3)
	qresult(2)=quatsymm(4,lindex)*qq(2)+quatsymm(2,lindex)*qq(4)
     &	+quatsymm(1,lindex)*qq(3)-quatsymm(3,lindex)*qq(1)
	qresult(3)=quatsymm(4,lindex)*qq(3)+quatsymm(3,lindex)*qq(4)
     &	+quatsymm(2,lindex)*qq(1)-quatsymm(1,lindex)*qq(2)
	qresult(4)=quatsymm(4,lindex)*qq(4)-quatsymm(1,lindex)*qq(1)
     &	-quatsymm(2,lindex)*qq(2)-quatsymm(3,lindex)*qq(3)
c
c	write(*,*) 'qresult ',qresult
c
	return
	end
c
c
c _____________________________
c
c
	subroutine postsymm(qq,lindex,qresult)
	real qq(4),qresult(4)
	integer lindex
c
	include 'common.f'
	common/a1/ quatsymm(4,48),numsymm
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion
c  from applying a symmetry operator, QUATSYMM(n,lindex)
c  to the second quaternion
c
c  If the symm oper == O, then
c  Qresult = (QQ x O^-1) where QQ= Q1 x Q2^-1
c  i.e. Qresult = Q1 x (O x Q2)^-1
c  note change of signs to get inverse of first orientation
c
	if(lindex.gt.numsymm) stop 'error in postsymm, lindex>numsymm'
	if(lindex.lt.1) stop 'error in postsymm, lindex<1'
	qresult(1)=qq(1)*quatsymm(4,lindex)-qq(4)*quatsymm(1,lindex)
     &	+qq(2)*quatsymm(3,lindex)-qq(3)*quatsymm(2,lindex)
	qresult(2)=qq(2)*quatsymm(4,lindex)-qq(4)*quatsymm(2,lindex)
     &	+qq(3)*quatsymm(1,lindex)-qq(1)*quatsymm(3,lindex)
	qresult(3)=qq(3)*quatsymm(4,lindex)-qq(4)*quatsymm(3,lindex)
     &	+qq(1)*quatsymm(2,lindex)-qq(2)*quatsymm(1,lindex)
	qresult(4)=qq(4)*quatsymm(4,lindex)+qq(1)*quatsymm(1,lindex)
     &	+qq(2)*quatsymm(2,lindex)+qq(3)*quatsymm(3,lindex)
c
c	write(*,*) 'postsymm: qresult ',qresult
c
	return
	end
c
c
c _____________________________
c
c
	subroutine disquat2(qq1,qq2,index1,index2,qresult,index1a,index2a)
	real qq1(4),qq2(4),qqn(4),qresult(4),quint(4),quintn(4),qmax
	real q4neg
	integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp
c
c       PI=3.14159265
c
c  algorithm for finding disorientation corresponding to 
c  input quaternion, QQ (from COMQUAT):  resultant quat. is in QRESULT
c  indices for symmetry operators are in INDEX1, INDEX2
c
	index1=0
	index2=0
	index1a=0
	index2a=0
	qmax=0.
	ip=1
	jp=1
c
	do 200, i=1,24
		 call postsymm(qq1,i,quint)
c
		do 195, jjm=1,2
			if(jjm.eq.2) then
				do 170, k=1,4
					quint(k)=quint(k)*-1.0
170					continue
c  can take negative of Q because identity = +/-(0,0,0,1)
				endif
c
		do 190, j=1,24
c  loop over the two sides of the grain boundary
			call postsymm(qq2,j,quintn)
c
			call comquat2(quint,quintn,qqn)
		do 185, kjl=1,2
			if(kjl.eq.2) then
				do 180, k=1,4
					qresult(k)=qqn(k)*-1.0
180					continue
c  can take negative of Q because identity = +/-(0,0,0,1)
				endif
c
	do 210, ijk=1,2
			q4neg=1.0
			if(ijk.gt.1) q4neg=-1.0
					qresult(4)=qresult(4)*(q4neg)
c  take the inverse rotation, i.e. apply the switching symmetry
c
c	write(*,"('i,j,kjl,q4neg,qresult ',3i3,5f8.3)") i,j,kjl,q4neg,qresult
			if(qresult(1).ge.0.0) then
			  if((qresult(1)-qresult(2)).lt.1e-5) then
			    if((qresult(2)-qresult(3)).lt.1e-5) then
			      if((qresult(3)-qresult(4)).lt.1e-5) then
c   these IFs ensure that 0<=q1<q2<q3<q4 to place it in the fundamental zone
			        if(qresult(4).gt.qmax) then
c  so now we test to see if we have a smaller angle
			        qmax=qresult(4)
	write(*,"('i,j,kjl,q4neg,qresult ',3i3,5f8.3)") i,j,kjl,q4neg,qresult
			        index1=i
			        index2=j
c			        index3=0
				index1a=0
			        if(ijk.eq.2) index1a=1
			        index2a=0
			        if(kjl.eq.2) index2a=1
c  record the result if we have found a quat in the fund. zone, and a smaller angle
				  endif
				endif
			       endif
			      endif
			     endif
210		continue
185				continue
190				continue
195			continue
200			continue
c
	if(index1.eq.0.or.index2.eq.0) stop 'DISQUAT2 failed!'
c
		 call postsymm(qq1,index1,quint)
c
		 call postsymm(qq2,index2,quintn)
c
			call comquat2(quint,quintn,qresult)
c
			if(index1a.eq.1) qresult(4)=-1.*qresult(4)
c
			if(index2a.eq.1) then
c  if we had to take the negative of the Q, then do so
c  to the answer
c
				do 280, k=1,3
					qresult(k)=qresult(k)*-1.0
280					continue
				if(index1a.eq.1) qresult(4)=-1.*qresult(4)
c
c  can take negative of Q because identity = +/-(0,0,0,1)
			endif
	return
	end
c
c
c _____________________________
c
c
	subroutine disquat(qq,index1,index2,qresult,index1a,index2a)
	real qq(4),qqn(4),qresult(4),quint(4),quintn(4),qmax
	integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp
c
c       PI=3.14159265
c
c  algorithm for finding disorientation corresponding to 
c  input quaternion, QQ (from COMQUAT):  resultant quat. is in QRESULT
c  indices for symmetry operators are in INDEX1, INDEX2
c
	index1=0
	index2=0
	index1a=0
	index2a=0
	qmax=0.
	ip=1
	jp=1
c
	do 200, i=1,24
		 call presymm(qq,i,quint)
c
		do 190, j=1,24
c  loop over the two sides of the grain
			call postsymm(quint,j,quintn)
c
		do 185, kjl=1,2
		    switch=1.0
			if(kjl.eq.2) then
			switch=-1.0
			endif
				do 180, k=1,4
					qqn(k)=quintn(k)*switch
180					continue
c  can take negative of Q because identity = +/-(0,0,0,1)
c
	do 210, ijk=1,2
		qresult(1)=qqn(1)
		qresult(2)=qqn(2)
		qresult(3)=qqn(3)
		qresult(4)=qqn(4)
			iq4neg=ijk
	if(ijk.eq.2) qresult(4)=-1.*qqn(4)
c  take the inverse rotation, i.e. apply the switching symmetry
c
			if(qresult(1).ge.0.0) then
c			  if((qresult(1)-qresult(2)).lt.1e-5) then
c			    if((qresult(2)-qresult(3)).lt.1e-5) then
c			      if((qresult(3)-qresult(4)).lt.1e-5) then
			  if(qresult(1).le.qresult(2)) then
			    if(qresult(2).le.qresult(3)) then
			      if(qresult(3).le.qresult(4)) then
c   these IFs ensure that 0<=q1<q2<q3<q4 to place it in the fundamental zone
			        if(qresult(4).gt.qmax) then
c  so now we test to see if we have a smaller angle
			        qmax=qresult(4)
			        index1=i
			        index2=j
c			        index3=0
				index1a=0
			        if(ijk.eq.2) index1a=1
			        index2a=0
			        if(kjl.eq.2) index2a=1
c  record the result if we have found a quat in the fund. zone, and a smaller angle
				  endif
				endif
			       endif
			      endif
			     endif
210		continue
185				continue
190				continue
200			continue
c
	if(index1.eq.0.or.index2.eq.0) stop 'DISQUAT failed!'
c
		 call presymm(qq,index1,quint)
c
		 call postsymm(quint,index2,qqn)
c		 write(*,"('DISQUAT: qqn = ',4f8.3)") qqn
		do 280, k=1,4
			if(index2a.eq.1) then
c  if we had to take the negative of the Q, then do so
c  to the answer
c
				qresult(k)=qqn(k)*-1.0
c  can take negative of Q because identity = +/-(0,0,0,1)
			else
				qresult(k)=qqn(k)
			endif
280		continue
c
	if(index1a.eq.1) then
		qresult(4)=-1.*qresult(4)
c  restore the 4th component if disorient found for ijk=1
	endif
c
c	write(*,*) 'DISQUAT: qmax= ',qmax
	return
	end
c
c
c  ----------------
c
	subroutine axisang2rot(axis,ang,rot)
c
c  AXIS  is the normalized axis vector
c  ANG is the rotation angle in radians
c  ROT is the 3x3 matrix
	real cc,ss,rnorm,rsum,axis(3),axis2(3),ang,rot(3,3)
c
	rsum=0.
	do 100, i=1,3
	axis2(i)=axis(i)**2
	rsum=rsum+axis2(i)
100	continue
	rnorm=sqrt(rsum)
	do 200, i=1,3
	axis(i)=axis(i)/rnorm
	axis2(i)=axis(i)**2
200	continue
c  normalize anyway!
	cc=cos(ang)
	ss=sin(ang)
	rot(1,1)=axis2(1)+cc*(1-axis2(1))
	rot(2,2)=axis2(2)+cc*(1-axis2(2))
	rot(3,3)=axis2(3)+cc*(1-axis2(3))
	rot(1,2)=(1.-cc)*axis(1)*axis(2)+ss*axis(3)
	rot(2,1)=(1.-cc)*axis(1)*axis(2)-ss*axis(3)
	rot(2,3)=(1.-cc)*axis(2)*axis(3)+ss*axis(1)
	rot(3,2)=(1.-cc)*axis(2)*axis(3)-ss*axis(1)
	rot(3,1)=(1.-cc)*axis(1)*axis(3)+ss*axis(2)
	rot(1,3)=(1.-cc)*axis(1)*axis(3)-ss*axis(2)
c
	return
	end
c
c  --------------------
c
	function scalar(a,b)
c  return SCALAR PRODUCT of A and B
	real scalar,a(3),b(3)
	scalar=0.
	do 100, i=1,3
	scalar=scalar+a(i)*b(i)
100	continue
	return
	end
c
c  -------------------
c
	subroutine axvecb2vecc(a,b,c)
c  transforms vector B by A to get vector C
	real a(3,3),b(3),c(3)
c
	do 100, i=1,3
	c(i)=0.
	do 100, j=1,3
	c(i)=c(i)+a(i,j)*b(j)
100	continue
	return
	end
c
c  -------------------
c
	subroutine donorm(a)
	real a(3),b(3),rsum,rnorm
	rsum=0.
	do 100, i=1,3
	rsum=rsum+a(i)**2
100	continue
	rnorm=sqrt(rsum)
	do 200, i=1,3
	b(i)=a(i)/rnorm
	a(i)=b(i)
200	continue
	return
	end
c
c  -------------------
c
	subroutine vecneg(a)
	real a(3)
	do 100, i=1,3
	a(i)=-1.*a(i)
100	continue
	return
	end
c
c  -------------------
c
	subroutine raneuler(iisave,pi2,phi1,phi,phi2)
c  assumes that the RAN1 function is properly set up
c  in the calling program
c
	phi1=4*pi2*ran1(iisave)
	phi2=4*pi2*ran1(iisave)
	phi=acos(ran1(iisave))
	return
	end
c
c  -------------------
c
      subroutine euler(a,iopt,nomen,d1,d2,d3,ior,kerr)
c                          Last revision 20nov90 UFK
c      common a(3,3),grvol(1152),epsga(5),ist1,ist2,sqr3,sqrh,ident(3,3)
c  SPECIAL VERSION WITHOUT COMMON BLOCK
	real a(3,3),d1,d2,d3,th,sth,cth,sph,cph,sps,cps,ps,ph,dth,dph,dps
      character nomen
      save kor
      rad=57.29578
c
c *** this subroutine calculates the euler angles associated with the
c     transformation matrix a(i,j) if iopt=1 and viceversa if iopt=2
c ***  Note that a is sample (rows) in terms of crystal (columns);
c      -- opposite of standard g (e.g.Bunge) - this is Canova's
c ***  Note that in this version, the Euler angles are defined symmetrically:
c      so that interchanging phi and psi means transposing a.
c      ("Kocks" nomen: defined going from +X to +Y in both COD and SOD)
c ***  However, other angle conventions are translated, according to
c  nomen="K" - Kocks (as internally) -- also sometimes "N"...
c        "R" - Roe          (Psi=psi, Phi=180-phi)
c        "B" - Bunge        (phi1=90+psi, Phi=Theta, phi2=90-phi)
c        any other - Canova (Theta first, phiC=90+phi, omega=90-psi)
c ***   Note: only in symm. notation does a point with all Euler angles
c             between 0 and 90 deg appear in the +x,+y quadrant!
c       If you want to see an individual point in this quadrant and are:
c       using Roe,    the third angle must be between 90 and 180;
c             Bunge,      first                                 ;
c             Canova,     second                                .
c ***  Input and output Euler angles d1,d2,d3 in degrees
c
      goto(5,20),iopt
    5 if(abs(a(3,3)) .ge. 0.999) goto 10
      th=acos(a(3,3))
      sth=sin(th)
c
	if(abs(a(2,3)/sth.lt.1e-35).and.abs(a(1,3)/sth.lt.1e-35)) then
		ps=pi/4.
	else
      ps=atan2(a(2,3)/sth,a(1,3)/sth)
	endif
c
	if(abs(a(3,2)/sth.lt.1e-35).and.abs(a(3,1)/sth.lt.1e-35)) then
		ph=pi/4.
	else
      ph=atan2(a(3,2)/sth,a(3,1)/sth)
	endif
ccccc   if it bombs out here, both a-components in one arg. are zero
c       (this should not be possible, but has happened, probably fixed)
c  ADR: i 01 - above is the fix!!
c
      go to 15
c
   10 	if(abs(a(1,2)/sth.lt.1e-35).and.abs(a(1,1)/sth.lt.1e-35)) then
		ps=pi/4.
	else
		ps=0.5*atan2(a(1,2),-a(1,1))
	endif
c
      ph=-ps
c       The above still have the problem that they give too many
c       equivalents of the same grain in DIOROUT and density file. Therefore:
      if(kerr.eq.1.and.kor.ne.ior) then
      print*,'NOTE: grain',ior,' has Theta < 1 deg.: sometimes problems'
        print*
        kor=ior
      endif
c
   15 dth=th*rad
      dph=ph*rad
      dps=ps*rad
      d1=dps
      d2=dth
      if(nomen.eq.'K'.or.nomen.eq.'N') then
        d3=dph
      elseif(nomen.eq.'R') then
        d3=180.-dph
      elseif(nomen.eq.'B') then
        d1=dps+90.
        d3=90.-dph
      else
        d1=dth
        d2=dph+90.
        d3=90.-dps
      endif
      if(d1.ge.360.) d1=d1-360.
      if(d3.ge.360.) d3=d3-360.
      if(d1.lt.0.) d1=d1+360.
      if(d3.lt.0.) d3=d3+360.
      return
c  *************************************
   20 dth=d2
      dps=d1
      if(nomen.eq.'K') then
        dph=d3
      elseif(nomen.eq.'R') then
        dph=180.-d3
      elseif(nomen.eq.'B') then
        dps=d1-90.
        dph=90.-d3
      else
        dth=d1
        dph=d2-90.
        dps=90.-d3
      endif
      ph=dph/rad
      th=dth/rad
      ps=dps/rad
      sph=sin(ph)
      cph=cos(ph)
      sth=sin(th)
      cth=cos(th)
      sps=sin(ps)
      cps=cos(ps)
      a(1,1)=-sps*sph-cph*cps*cth
      a(2,1)=cps*sph-cph*sps*cth
      a(3,1)=cph*sth
      a(1,2)=sps*cph-sph*cps*cth
      a(2,2)=-cph*cps-sph*sps*cth
      a(3,2)=sth*sph
      a(1,3)=sth*cps
      a(2,3)=sps*sth
      a(3,3)=cth
      return
      end
c
c  &&&&&&&&&&&&&&
c
c  these subroutines copied from map5d
c ____________________________
c
c
	subroutine m2r(e,rodr,axis)
	real e(3,3),rodr(4),tmp(3),rnorm
	real det,detm,pi,axis(3)
c
c  converts a MISOR matrix to Rodrigues vector
c
      pi=4.*atan(1.)
c  calc. determinant
c
c	det=e(1,1)*e(2,2)*e(3,3)
c     &	+e(2,1)*e(3,2)*e(1,3)
c     &	+e(3,1)*e(1,2)*e(2,3)
c     &	-e(1,1)*e(3,2)*e(2,3)
c     &	-e(1,2)*e(2,1)*e(3,3)
c     &	-e(1,3)*e(2,2)*e(3,1)
c	if(det.lt.-0.95) det=-1.0
c	if(det.gt.0.95) det=1.0
c	detm=det
c
c  we are only interested in proper vs. improper rotations here
c
c	write(*,*) 'rotation matrixes:  ',(e(ii,1),ii=1,3)
c	write(*,*) 'rotation matrixes:  ',(e(ii,2),ii=1,3)
c	write(*,*) 'rotation matrixes:  ',(e(ii,3),ii=1,3)
c	write(*,*) 'index,  determinant= ',i,det
c	rtmp=(((e(1,1)+e(2,2)+e(3,3))/det)-1.)/2.
	rtmp=(((e(1,1)+e(2,2)+e(3,3)))-1.)/2.
	if(rtmp.lt.-1.) rtmp=-1.
	if(rtmp.gt.1.) rtmp=1.
	rodr(4)=acos(rtmp)
c  magnitude of rotation
c	tmp(1)=(e(2,3)-e(3,2))/det
c	tmp(2)=(e(1,3)-e(3,1))/det
c	tmp(3)=(e(1,2)-e(2,1))/det
	axis(1)=(e(2,3)-e(3,2))
	axis(2)=(e(1,3)-e(3,1))
	axis(3)=(e(1,2)-e(2,1))
c	rnorm=sqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
	rnorm=sqrt(axis(1)**2+axis(2)**2+axis(3)**2)
	if(rnorm.lt.1e-5) then
c		tmp(1)=sqrt(1.+e(1,1))
c		tmp(2)=sqrt(1.+e(2,2))
c		tmp(3)=sqrt(1.+e(3,3))
		axis(1)=sqrt(1.+e(1,1))
		axis(2)=sqrt(1.+e(2,2))
		axis(3)=sqrt(1.+e(3,3))
c		rnorm=sqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
		rnorm=sqrt(axis(1)**2+axis(2)**2+axis(3)**2)
	endif
c	write(*,*) rnorm, tmp(1),tmp(2),tmp(3),tan(rodr(4)/2)
c	rodr(1)=tmp(1)/rnorm*tan(rodr(4)/2.)
c	rodr(2)=tmp(2)/rnorm*tan(rodr(4)/2.)
c	rodr(3)=tmp(3)/rnorm*tan(rodr(4)/2.)
	rodr(1)=axis(1)/rnorm*tan(rodr(4)/2.)
	rodr(2)=axis(2)/rnorm*tan(rodr(4)/2.)
	rodr(3)=axis(3)/rnorm*tan(rodr(4)/2.)
c  components of the Rodrigues vector
c	write(*,190)  rodr(4)*180./pi
c     &	,rodr(1),rodr(2),rodr(3)
190	format(' angle, vector....',4f9.3)
c
100	continue
	return
	end
c
c
c  ____________________________
c
c
	subroutine eul2rod(d1,d2,d3,rodr)
	real d1,d2,d3,rodr(3),sum,diff,csum,cdiff,sdiff
c  converts Bunge EUler angles to Rodrigues vector
c  d1== phi1,  d2== Phi,  d3== phi2
c
	sum=(d1+d3)/2.
	diff=(d1-d3)/2.
	csum=cos(sum)
	cdiff=cos(sum)
	sdiff=sin(diff)
	t2=tan(d2/2.)
c
	rodr(1)=t2*cdiff/csum
	rodr(2)=t2*sdiff/csum
	rodr(3)=tan(sum)
c
	return
	end
c
c
c  ____________________________
c
c
	subroutine rod2eul(d1,d2,d3,rodr)
	real d1,d2,d3,rodr(3),sum,diff
c  converts Rodrigues vector to Bunge Euler angles 
c  d1== phi1,  d2== Phi,  d3== phi2
c
	sum=atan(rodr(3))
	diff=atan(rodr(2)/rodr(1))
	d1=sum+diff
	d2=2.*atan(rodr(1)*cos(sum)/cos(diff))
	d3=sum-diff
c
	return
	end
c
c
c  ____________________________
c
	subroutine symm_apply(eee,f,index)
	real r(3,3,48)
	real eee(3,3),f(3,3)
	integer index,ia,ib,ic
c  matrix multiply, for symmetry operations
c
	do 1940, ia=1,3
	do 1940, ib=1,3
	f(ia,ib)=0.
	do 1940, ic=1,3
	f(ia,ib)=f(ia,ib)+eee(ia,ic)*r(ib,ic,index)
c  postmultiply by transpose of the symmetry operator (c1)
c  
1940	continue
c
	return
	end
c
c  ____________________________
c
	subroutine symm_pre_apply(eee,f,index)
	real r(3,3,48)
	real eee(3,3),f(3,3)
	integer index,ia,ib,ic
c  matrix multiply, for symmetry operations
c
	do 1940, ia=1,3
	do 1940, ib=1,3
	f(ia,ib)=0.
	do 1940, ic=1,3
	f(ia,ib)=f(ia,ib)+r(ia,ic,index)*eee(ic,ib)
c  premultiply by the symmetry operator (c1)
c  
1940	continue
c
	return
	end
c
c
c  -----------------------------
c
	subroutine dis2mat(ein1,ein2,eout,rout,angle,axis,index1,index2,invert)
c
c  analyzes EIN as a rotation matrix (assumed!) to find
c  the disorientation and return the corresponding matrix, EOUT,
c  the rotation ANGLE (degrees) and the AXIS of rotation (in the SST)
c  INDEX1 is the index to the first symm oper, INDEX2, the second
c  INVERT tells whether the rotations were taken in the opposite order
c  ROUT is the Rodrigues vector of the result
c
c  rotation matrices supplied by EULER and are therefore active
c  rotations from ref. to crystal, or axis transformations from
c  crystal to ref
c
c  a misorientation in the local frame must be calculated as 
c  E1^T.E2, so with symmetry operators applied,
c  (E1.O1)^T(E2.O2)
c
	real ein1(3,3),ein2(3,3),eout(3,3),angle,axis(3)
	real ea(3,3),eb(3,3),easym(3,3),ebsym(3,3),ec(3,3)
	real tracemax,anglemin,rout(3),rodr(4)
	integer index1,index2,invert
c
	pi=4.*atan(1.)
	degrad=180./pi
	anglemin=pi
	index1=0
	index2=0
	invert=0
c
	do 3010, ii=1,2
c
		do 100, i=1,3
		do 100, j=1,3
100		ea(i,j)=ein1(i,j)
c
		do 110, i=1,3
		do 110, j=1,3
110		eb(i,j)=ein2(i,j)
c
	if(ii.eq.2) then
			do 200, i=1,3
			do 200, j=1,3
200			ea(i,j)=ein2(i,j)
c
			do 210, i=1,3
			do 210, j=1,3
210			eb(i,j)=ein1(i,j)
c
	endif
c
		do 3000, i=1,24
			call symm_apply(ea,easym,i)
c
			do 2990, j=1,24
				call symm_apply(eb,ebsym,i)
c
				call atxb2c(easym,ebsym,ec)
				call m2r(ec,rodr,axis)
				if(rodr(4).lt.anglemin) then
		if(((rodr(1).ge.0..and.rodr(2).ge.0..and.rodr(3).ge.0.).and.
     &	(rodr(1).ge.rodr(2).and.rodr(2).ge.rodr(3)))
     &	) then
					anglemin=rodr(4)
					index1=i
					index2=j
					invert=ii
					angle=rodr(4)*degrad
					rout(1)=rodr(1)
					rout(2)=rodr(2)
					rout(3)=rodr(3)
				endif
			endif
c
2990	continue
c  end of inner loop
3000	continue
c  end of outer loop
3010	continue
c  loop for inverse of rotation
c
	if(index1.eq.0) stop 'dis2mat failed!'
c
	write(*,*) 'DIS2MAT result: '
	write(*,"('DIS2MAT angle, rodrigues vec: ',4f8.3)")  angle,rout
	write(*,"('DIS2MAT  indices: ',3i4)") index1,index2,invert
c
	return
	end
c
c  -----------------------------
c
	subroutine readsigmas
      real qq(3,3,50),oc(50)
      integer numsigs,i
      common /m2/ qq,oc,numsigs
c
      open(4,file='misdat',status='old')
	do 50, i=1,8
	read(4,*) 
	read(4,*) 
	read(4,*) 
c  skip over the symmetry matrices!
c
50	continue
c
	do 100, i=1,50
	write(*,*) ' sigma no.  ', i
	read(4,*,end=200) ((qq(ii,jj,i),ii=1,3),jj=1,3),oc(i)
c
c   NOTE the transposition of the input matrices
c  this is done so that when we form a Rodrigues vector from them
c  , if needed, the result is still in the fundamental zone!
c   but it means that one must (post-)multiply with the transpose
c  of the CSL matrix!
c   ADR  vi 99
c
100	continue
	write(*,*) 'There may have been more than 50 Sigma matrices....'
200	numsigs=i-1
	write(*,*) 'number of sigma values entered = ',numsigs
      return
      end
c
c ------------------------------------------------
c
	subroutine rodcom(r1,r2,rout)
c
c  combines two rotations as Rodrigues vectors
c  ROUT= R1+R2- R1 cross R2 / 1 - (R1dotR2)
	real r1(3),r2(3),rout(3)
	real dot,cross(3)
	dot=1.-scalar(r1,r2)
	call vecpro2(r1,r2,cross)
	if(abs(dot).gt.1e-5) then
		rout(1)=(r1(1)+r2(1)-cross(1))/dot
		rout(2)=(r1(2)+r2(2)-cross(2))/dot
		rout(3)=(r1(3)+r2(3)-cross(3))/dot
	else
		rout(1)=1e38
		rout(2)=1e38
		rout(3)=1e38
		write(*,*) 'warning: RODCOM ->infinite result'
	endif
	return
	end
c
c
c ------------------------------------------------
c
	subroutine rodcomABT(r1,r2,rout)
c
c  combines two rotations as Rodrigues vectors
c  ROUT= R1+R2- R1 cross R2 / 1 - (R1dotR2)
c  except use the transpose of second rotation
c
	real r1(3),r2(3),rout(3)
	real dot,cross(3)
	call vecneg(r2)
	dot=1.-scalar(r1,r2)
	call vecpro2(r1,r2,cross)
	if(abs(dot).gt.1e-5) then
		rout(1)=(r1(1)+r2(1)-cross(1))/dot
		rout(2)=(r1(2)+r2(2)-cross(2))/dot
		rout(3)=(r1(3)+r2(3)-cross(3))/dot
	else
		rout(1)=1e38
		rout(2)=1e38
		rout(3)=1e38
		write(*,*) 'warning: RODCOM ->infinite result'
	endif
	return
	end
c
c
c ------------------------------------------------
c
	subroutine rod2ang(r1,angle)
	real r1(3),angle
c  extracts a rotation angle from a Rodrigues vector
c
	angle=2.*atan((sqrt(r1(1)**2+r1(2)**2+r1(3)**2)))
	return
	end
c
c
c ------------------------------------------------
c
	subroutine rodcomATB(r1,r2,rout)
c
c  combines two rotations as Rodrigues vectors
c  ROUT= R1+R2- R1 cross R2 / 1 - (R1dotR2)
c  except use the transpose of first rotation
c
	real r1(3),r2(3),rout(3)
	real dot,cross(3)
	call vecneg(r1)
	dot=1.-scalar(r1,r2)
	call vecpro2(r1,r2,cross)
	if(abs(dot).gt.1e-5) then
		rout(1)=(r1(1)+r2(1)-cross(1))/dot
		rout(2)=(r1(2)+r2(2)-cross(2))/dot
		rout(3)=(r1(3)+r2(3)-cross(3))/dot
	else
		rout(1)=1e38
		rout(2)=1e38
		rout(3)=1e38
		write(*,*) 'warning: RODCOM ->infinite result'
	endif
	return
	end
c
c ------------------------------------------------
c
	subroutine CSL(rodrvec,toler,sigmin,sigma)
c
c	ADR  19 vi 99
c
c  NOW FOR SIGMA analysis; form a misorientation matrix from the results
c
c  NEEDS READSIGMAS to be executed before CSL is called
c
      real qq(3,3,50),oc(50)
      integer numsigs,i
      common /m2/ qq,oc,numsigs
c
	real rodrvec(4),s1,s2,s3,u(3,3)
	real rj1,rj2,rj3,rnorm,rb,rn1,rn2,rn3,rn4,rn5
	real axis(3),qsig(3,3),cresult(3,3)
	real sigmin,toler,sigma,pi,pi2,degrad
	real rodr(4)
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
	sigmin=180.0
	sigma=0.0
c
	do 10, i=1,3
	axis(i)=rodrvec(i)
10	continue
c
	call  axisang2rot(axis,rodrvec(4),u)
c
	do 3000, ijk=1,numsigs
c	write(*,*) 'working on sigma = ',int(oc(ijk))
c
	do 20, i=1,3
	do 20,j=1,3
20	qsig(i,j)=qq(i,j,ijk)
c
c	write(*,320) (u(1,j1),j1=1,3),(qsig(1,j1),j1=1,3)
c	write(*,320) (u(2,j1),j1=1,3),(qsig(2,j1),j1=1,3)
c	write(*,320) (u(3,j1),j1=1,3),(qsig(3,j1),j1=1,3)
320	format(x,3(x,f8.3),6x,3(x,f8.3))
		call axbt2c(u,qsig,cresult)
	rnorm=(cresult(1,1)+cresult(2,2)+cresult(3,3)-1.)/2.
c  rnorm has the trace of the resulting rotation matrix
	if(rnorm.lt.-1.) rnorm=-1.
	if(rnorm.gt.1.) rnorm=1.
	rn4=acos(rnorm)
c  RN4 has the magnitude of rotation in radians
c
c	call m2r(cresult,rodr,det)
c	write(*,190)  rodr(4)*180./pi
c     &	,rodr(1),rodr(2),rodr(3)
190	format('  after combining:  angle, vector....',4f9.3)
c
      rb=(15.0*1./sqrt(oc(ijk)))
c  RB  is for Brandon's criterion in degrees
      rn5=rn4*degrad/rb
c  outputs CSL and deviation
c      if (rn5.lt.1.0) then 
      if (rn5.lt.toler) then 
c      write(2,*) 'Sigma is ',int(oc(ijk)),'; V/Vm =  ',rn5
c        write(*,*) 'Sigma is ',int(oc(ijk)),'; V/Vm =  ',rn5
      endif
c
	if((rn5.lt.toler).and.(rn5.lt.sigmin)) then
		sigma=oc(ijk)
		sigmin=rn5
	endif
3000	continue
	return
	end
c
c  -------------------
c
	subroutine axb2c(a,b,c)
c  multiply a by b to get c
	real a(3,3),b(3,3),c(3,3)
c
	do 100, i=1,3
	do 100, j=1,3
	c(i,j)=0.0
	do 100,ijk=1,3
	c(i,j)=c(i,j)+a(i,ijk)*b(ijk,j)
100	continue
	return
	end
c
c  -------------------
c
	subroutine atxb2c(a,b,c)
c  multiply a-TRANSPOSE by b to get c
	real a(3,3),b(3,3),c(3,3)
c
	do 100, i=1,3
	do 100, j=1,3
	c(i,j)=0.0
	do 100,ijk=1,3
	c(i,j)=c(i,j)+a(ijk,i)*b(ijk,j)
100	continue
	return
	end
c
c  -------------------
c
	subroutine axbt2c(a,b,c)
c  multiply a by b-TRANSPOSE to get c
	real a(3,3),b(3,3),c(3,3)
c
	do 100, i=1,3
	do 100, j=1,3
	c(i,j)=0.0
	do 100,ijk=1,3
	c(i,j)=c(i,j)+a(i,ijk)*b(j,ijk)
100	continue
	return
	end
c
c ____________________________
c
      subroutine vecpro2(v1,v2,vout)
c  vector product
	real v1(3),v2(3),vout(3)
	vout(3)=v1(1)*v2(2)-v1(2)*v2(1)
	vout(1)=v1(2)*v2(3)-v1(3)*v2(2)
	vout(2)=v1(3)*v2(1)-v1(1)*v2(3)
	return
	end
c
c ____________________________
c
      subroutine vecpro(k,t1)
c  based on LApp routine
c     *****************
      real t1(3,3)
c
c *** calculates axis x(k) as vector product of x(k1),x(k2) taken cyclically
c
      k1=k+1
      if(k1.gt.3) k1=k1-3
      k2=k+2
      if(k2.gt.3) k2=k2-3
      t1(k,1)=t1(k1,2)*t1(k2,3)-t1(k1,3)*t1(k2,2)
      t1(k,2)=t1(k1,3)*t1(k2,1)-t1(k1,1)*t1(k2,3)
      t1(k,3)=t1(k1,1)*t1(k2,2)-t1(k1,2)*t1(k2,1)
      return
c  look for vector product in the kth column of t1
      end
c
c -----------------------------------------
c
C subroutine ZUR BERECHNUNG VON ORIENTIERUNGSBEZIEHUNGEN ZWISCHEN
C ZWEI VORGEGEBENEN ORIENTIERUNGEN
c
c  subroutine provided by Dierk Raabe (to Paul Lee, 1996)
c
      subroutine misori(p1,p,p2,al_min,nsymm,ac_min)
c  
c  p1,p,p2 are (2 sets of) Bunge Euler angles
c  nsymm controls the sample symmetry
c  al_min is the disorientation angle
c
       character meu,isy
       integer ac_min(3)
       DIMENSION D(3,3,2),SYM(3,3,24),P1(2),P(2),P2(2)
     $,DM(3,3),DRR(3,3,24),A(3),IA(3)
     $,DR(3,3)
       PI=3.14159265

	al_min=360.
	ac_min(1)=0
	ac_min(2)=0
	ac_min(3)=0
C ERSTELLEN DER SYMMETRIEMATRIZEN 
C **************************************************************
       DO 2 I=1,3
       DO 2 J=1,3
       DO 2 K=1,24
2      SYM(I,J,K)=0.
C  1
       SYM(1,1,1)=1.
       		SYM(2,2,1)=1.
       				SYM(3,3,1)=1.
C  5
       SYM(1,1,2)=1.
       				SYM(2,3,2)=-1.
       		SYM(3,2,2)=1.
C  2
       SYM(1,1,3)=1.
       		SYM(2,2,3)=-1.
       				SYM(3,3,3)=-1.
C  11
       SYM(1,1,4)=1.
       				SYM(2,3,4)=1.
       		SYM(3,2,4)=-1.
C  7
       				SYM(1,3,5)=-1.
       		SYM(2,2,5)=1.
       SYM(3,1,5)=1.
C  12
       				SYM(1,3,6)=1.
       		SYM(2,2,6)=1.
       SYM(3,1,6)=-1.
C  3
       SYM(1,1,7)=-1.
       		SYM(2,2,7)=1.
       				SYM(3,3,7)=-1.
C  4
       SYM(1,1,8)=-1.
       		SYM(2,2,8)=-1.
       				SYM(3,3,8)=1.
C  13
       		SYM(1,2,9)=1.
       SYM(2,1,9)=-1.
       				SYM(3,3,9)=1.
C  6
       		SYM(1,2,10)=-1.
       SYM(2,1,10)=1.
       				SYM(3,3,10)=1.
C  20
       		SYM(1,2,11)=-1.
       				SYM(2,3,11)=1.
       SYM(3,1,11)=-1.
C  23
       				SYM(1,3,12)=1.
       SYM(2,1,12)=-1.
       		SYM(3,2,12)=-1.
C  19
       		SYM(1,2,13)=-1.
       				SYM(2,3,13)=-1.
       SYM(3,1,13)=1.
C  21
       				SYM(1,3,14)=-1.
       SYM(2,1,14)=1.
       		SYM(3,2,14)=-1.
C  18
       		SYM(1,2,15)=1.
       				SYM(2,3,15)=-1.
       SYM(3,1,15)=-1.
C  22
       				SYM(1,3,16)=-1.
       SYM(2,1,16)=-1.
       		SYM(3,2,16)=1.
C  9
       		SYM(1,2,17)=1.
       				SYM(2,3,17)=1.
       SYM(3,1,17)=1.
C  8
       				SYM(1,3,18)=1.
       SYM(2,1,18)=1.
       		SYM(3,2,18)=1.
C  17
       		SYM(1,2,19)=1.
       SYM(2,1,19)=1.
       				SYM(3,3,19)=-1.
C  24
       SYM(1,1,20)=-1.
       				SYM(2,3,20)=1.
       		SYM(3,2,20)=1.
C  10
       				SYM(1,3,21)=1.
       		SYM(2,2,21)=-1.
       SYM(3,1,21)=1.
C  15
       SYM(1,1,22)=-1.
       				SYM(2,3,22)=-1.
       		SYM(3,2,22)=-1.
C  16
       				SYM(1,3,23)=-1.
       		SYM(2,2,23)=-1.
       SYM(3,1,23)=-1.
C  14
       		SYM(1,2,24)=-1.
       SYM(2,1,24)=-1.
       				SYM(3,3,24)=-1.
C **********************************
C ERSTELLEN DER BEIDEN DMATRIZEN
c 205    FORMAT(A1)

      meu='E'
      isy='E'
	I31MAX=nsymm

       DO 31 I31=1,I31MAX
       DO 1 I1=1,2
       FAKT=PI/180.
       PP1=P1(I1)*FAKT
       PP=P(I1)*FAKT
       PP2=P2(I1)*FAKT
       C1=COS(PP1)
       C=COS(PP)
       C2=COS(PP2)
       S1=SIN(PP1)
       S=SIN(PP)
       S2=SIN(PP2)
       D(1,1,I1)=C1*C2-S1*S2*C
       D(1,2,I1)=S1*C2+C1*S2*C
       D(1,3,I1)=S2*S
       D(2,1,I1)=-C1*S2-S1*C2*C
       D(2,2,I1)=-S1*S2+C1*C2*C
       D(2,3,I1)=C2*S
       D(3,1,I1)=S1*S
       D(3,2,I1)=-C1*S
       D(3,3,I1)=C
       GOTO 1
1      CONTINUE
C **********************************

       DO 24 I24=1,3
       DO 25 I25=1,3
       DM(I24,I25)=D(I24,I25,1)
       DR(I24,I25)=D(I24,I25,2)
25     CONTINUE
24     CONTINUE
       F=PI/180.
C
C ****************************************************
C BESTIMMUNG DER INVERSEN MATRIX ZUR ORIENTIERUNG 1:DM
C ****************************************************
       DO 3 I=1,3
       DO 3 J=1,3
3      DM(I,J)=D(J,I,1)
C *************************
C FELDER GLEICH NULL SETZEN
C *************************
       DO 4 I=1,3
       DO 4 J=1,3
       DO 4 K=1,24
       DR(I,J)=0.
4      DRR(I,J,K)=0.
C ***********************************************************
C MATRIZENMULTIPLIKATION DER MATRIZEN D(I,J,2) UND DM=DR(I,J)
C ***********************************************************
       DO 5 I=1,3
       DO 5 J=1,3
       DO 6 K=1,3
6      DR(I,J)=D(I,K,2)*DM(K,J)+DR(I,J)
5      CONTINUE
C ***********************************************
C MATRIZENMULITIPLIKATION DER MATRIZEN SYM UND DR
C ***********************************************
       DO 7 IZ=1,24
	DO 8 I=1,3
       DO 8 J=1,3
       DO 9 K=1,3
9      DRR(I,J,IZ)=SYM(I,K,IZ)*DR(K,J)+DRR(I,J,IZ)
8      CONTINUE
C *******************************
C BESTIMMUNG DES ROTATIONSWINKELS
C *******************************
       SPUR=0.
       DO 10 I=1,3
10     SPUR=SPUR+DRR(I,I,IZ)
       SP=(SPUR-1.)/2.
       IF(SP.GE.1.)SP=1.-1.E-5
       IF(SP.LE.-1.)SP=-1.+1.E-5
       OMEGA=PI/2.-ASIN(SP)
       IOMEGA=INT(OMEGA*180./PI+0.5)
       ALPHA=OMEGA*180./PI
       IF(IOMEGA.EQ.180.or.iomega.eq.0)GOTO 12
       X=2.*SIN(OMEGA)

C ********************************************************
C BESTIMMUNG DER ROTATIONSACHSE, INCL SONDERFALL OMEGA = 180
C ********************************************************
       A(1)=100.*(DRR(2,3,IZ)-DRR(3,2,IZ))/X
       A(2)=100.*(DRR(3,1,IZ)-DRR(1,3,IZ))/X
       A(3)=100.*(DRR(1,2,IZ)-DRR(2,1,IZ))/X
       GOTO 13
12     A(1)=SQRT((DRR(1,1,IZ)+1.)/2.)*100.
       A(2)=SQRT((DRR(2,2,IZ)+1.)/2.)*100.
       A(3)=SQRT((DRR(3,3,IZ)+1.)/2.)*100.
13     DO 11 I=1,3
       IF(A(I).EQ.0.)GOTO 19
       IA(I)=INT(A(I)+0.5*A(I)/ABS(A(I)))
       GOTO 11
19     IA(I)=0
11     CONTINUE

      if(abs(alpha).lt.al_min) then
	al_min=abs(alpha)
	ac_min(1)=ia(1)
	ac_min(2)=ia(2)
	ac_min(3)=ia(3)
	end if

7      CONTINUE
       IF(I31.EQ.1)P1(2)=180.+P1(2)
       IF(I31.EQ.2)P1(2)=360.-P1(2)
       IF(I31.EQ.2)P2(2)=90.-P2(2)
       IF(I31.EQ.3)P1(2)=180.+P1(2)


31     CONTINUE
      i_min=100
	if(abs(ac_min(1)).lt.i_min.and.abs(ac_min(1)).gt.0) 
     $i_min=abs(ac_min(1))
	if(abs(ac_min(2)).lt.i_min.and.abs(ac_min(2)).gt.0) 
     $i_min=abs(ac_min(2))
	if(abs(ac_min(3)).lt.i_min.and.abs(ac_min(3)).gt.0) 
     $i_min=abs(ac_min(3))

	do 5501 ii=1,3
5501  ac_min(ii)=nint(float(ac_min(ii))/float(i_min))


	return
	end
c
c
c -----------------------------------------
c
C subroutine ZUR BERECHNUNG VON ORIENTIERUNGSBEZIEHUNGEN ZWISCHEN
C ZWEI VORGEGEBENEN ORIENTIERUNGEN
c
c  subroutine provided by Dierk Raabe (to Paul Lee, 1996)
c
      subroutine misorinv(p1,p,p2,al_min,nsymm,ac_min)
c  
c  p1,p,p2 are (2 sets of) Bunge Euler angles
c  nsymm controls the sample symmetry
c  al_min is the disorientation angle
c
       character meu,isy
       integer ac_min(3)
       DIMENSION D(3,3,2),SYM(3,3,24),P1(2),P(2),P2(2)
     $,DM(3,3),DRR(3,3,24),A(3),IA(3)
     $,DR(3,3)
       PI=3.14159265

	al_min=360.
	ac_min(1)=0
	ac_min(2)=0
	ac_min(3)=0
C ERSTELLEN DER SYMMETRIEMATRIZEN 
C **************************************************************
       DO 2 I=1,3
       DO 2 J=1,3
       DO 2 K=1,24
2      SYM(I,J,K)=0.
C  1
       SYM(1,1,1)=1.
       		SYM(2,2,1)=1.
       				SYM(3,3,1)=1.
C  5
       SYM(1,1,2)=1.
       				SYM(2,3,2)=-1.
       		SYM(3,2,2)=1.
C  2
       SYM(1,1,3)=1.
       		SYM(2,2,3)=-1.
       				SYM(3,3,3)=-1.
C  11
       SYM(1,1,4)=1.
       				SYM(2,3,4)=1.
       		SYM(3,2,4)=-1.
C  7
       				SYM(1,3,5)=-1.
       		SYM(2,2,5)=1.
       SYM(3,1,5)=1.
C  12
       				SYM(1,3,6)=1.
       		SYM(2,2,6)=1.
       SYM(3,1,6)=-1.
C  3
       SYM(1,1,7)=-1.
       		SYM(2,2,7)=1.
       				SYM(3,3,7)=-1.
C  4
       SYM(1,1,8)=-1.
       		SYM(2,2,8)=-1.
       				SYM(3,3,8)=1.
C  13
       		SYM(1,2,9)=1.
       SYM(2,1,9)=-1.
       				SYM(3,3,9)=1.
C  6
       		SYM(1,2,10)=-1.
       SYM(2,1,10)=1.
       				SYM(3,3,10)=1.
C  20
       		SYM(1,2,11)=-1.
       				SYM(2,3,11)=1.
       SYM(3,1,11)=-1.
C  23
       				SYM(1,3,12)=1.
       SYM(2,1,12)=-1.
       		SYM(3,2,12)=-1.
C  19
       		SYM(1,2,13)=-1.
       				SYM(2,3,13)=-1.
       SYM(3,1,13)=1.
C  21
       				SYM(1,3,14)=-1.
       SYM(2,1,14)=1.
       		SYM(3,2,14)=-1.
C  18
       		SYM(1,2,15)=1.
       				SYM(2,3,15)=-1.
       SYM(3,1,15)=-1.
C  22
       				SYM(1,3,16)=-1.
       SYM(2,1,16)=-1.
       		SYM(3,2,16)=1.
C  9
       		SYM(1,2,17)=1.
       				SYM(2,3,17)=1.
       SYM(3,1,17)=1.
C  8
       				SYM(1,3,18)=1.
       SYM(2,1,18)=1.
       		SYM(3,2,18)=1.
C  17
       		SYM(1,2,19)=1.
       SYM(2,1,19)=1.
       				SYM(3,3,19)=-1.
C  24
       SYM(1,1,20)=-1.
       				SYM(2,3,20)=1.
       		SYM(3,2,20)=1.
C  10
       				SYM(1,3,21)=1.
       		SYM(2,2,21)=-1.
       SYM(3,1,21)=1.
C  15
       SYM(1,1,22)=-1.
       				SYM(2,3,22)=-1.
       		SYM(3,2,22)=-1.
C  16
       				SYM(1,3,23)=-1.
       		SYM(2,2,23)=-1.
       SYM(3,1,23)=-1.
C  14
       		SYM(1,2,24)=-1.
       SYM(2,1,24)=-1.
       				SYM(3,3,24)=-1.
C **********************************
C ERSTELLEN DER BEIDEN DMATRIZEN
c 205    FORMAT(A1)

      meu='E'
      isy='E'
	I31MAX=nsymm

       DO 31 I31=1,I31MAX
       DO 1 I1=1,2
       FAKT=PI/180.
       PP1=P1(I1)*FAKT
       PP=P(I1)*FAKT
       PP2=P2(I1)*FAKT
       C1=COS(PP1)
       C=COS(PP)
       C2=COS(PP2)
       S1=SIN(PP1)
       S=SIN(PP)
       S2=SIN(PP2)
       D(1,1,I1)=C1*C2-S1*S2*C
       D(1,2,I1)=S1*C2+C1*S2*C
       D(1,3,I1)=S2*S
       D(2,1,I1)=-C1*S2-S1*C2*C
       D(2,2,I1)=-S1*S2+C1*C2*C
       D(2,3,I1)=C2*S
       D(3,1,I1)=S1*S
       D(3,2,I1)=-C1*S
       D(3,3,I1)=C
       GOTO 1
1      CONTINUE
C **********************************

       DO 24 I24=1,3
       DO 25 I25=1,3
       DM(I24,I25)=D(I24,I25,1)
       DR(I24,I25)=D(I24,I25,2)
25     CONTINUE
24     CONTINUE
       F=PI/180.
C
C ****************************************************
C BESTIMMUNG DER INVERSEN MATRIX ZUR ORIENTIERUNG 1:DM
C ****************************************************
       DO 3 I=1,3
       DO 3 J=1,3
3      DM(I,J)=D(J,I,1)
C *************************
C FELDER GLEICH NULL SETZEN
C *************************
       DO 4 I=1,3
       DO 4 J=1,3
       DO 4 K=1,24
       DR(I,J)=0.
4      DRR(I,J,K)=0.
C ***********************************************************
C MATRIZENMULTIPLIKATION DER MATRIZEN D(I,J,2) UND DM=DR(I,J)
C ***********************************************************
       DO 5 I=1,3
       DO 5 J=1,3
       DO 6 K=1,3
6      DR(I,J)=D(k,i,2)*DM(j,k)+DR(I,J)
c  !!!!!  inverted the transposition here !!!!
5      CONTINUE
C ***********************************************
C MATRIZENMULITIPLIKATION DER MATRIZEN SYM UND DR
C ***********************************************
       DO 7 IZ=1,24
	DO 8 I=1,3
       DO 8 J=1,3
       DO 9 K=1,3
9      DRR(I,J,IZ)=SYM(I,K,IZ)*DR(K,J)+DRR(I,J,IZ)
8      CONTINUE
C *******************************
C BESTIMMUNG DES ROTATIONSWINKELS
C *******************************
       SPUR=0.
       DO 10 I=1,3
10     SPUR=SPUR+DRR(I,I,IZ)
       SP=(SPUR-1.)/2.
       IF(SP.GE.1.)SP=1.-1.E-5
       IF(SP.LE.-1.)SP=-1.+1.E-5
       OMEGA=PI/2.-ASIN(SP)
       IOMEGA=INT(OMEGA*180./PI+0.5)
       ALPHA=OMEGA*180./PI
       IF(IOMEGA.EQ.180.or.iomega.eq.0)GOTO 12
       X=2.*SIN(OMEGA)

C ********************************************************
C BESTIMMUNG DER ROTATIONSACHSE, INCL SONDERFALL OMEGA = 180
C ********************************************************
       A(1)=100.*(DRR(2,3,IZ)-DRR(3,2,IZ))/X
       A(2)=100.*(DRR(3,1,IZ)-DRR(1,3,IZ))/X
       A(3)=100.*(DRR(1,2,IZ)-DRR(2,1,IZ))/X
       GOTO 13
12     A(1)=SQRT((DRR(1,1,IZ)+1.)/2.)*100.
       A(2)=SQRT((DRR(2,2,IZ)+1.)/2.)*100.
       A(3)=SQRT((DRR(3,3,IZ)+1.)/2.)*100.
13     DO 11 I=1,3
       IF(A(I).EQ.0.)GOTO 19
       IA(I)=INT(A(I)+0.5*A(I)/ABS(A(I)))
       GOTO 11
19     IA(I)=0
11     CONTINUE

      if(abs(alpha).lt.al_min) then
	al_min=abs(alpha)
	ac_min(1)=ia(1)
	ac_min(2)=ia(2)
	ac_min(3)=ia(3)
	end if

7      CONTINUE
       IF(I31.EQ.1)P1(2)=180.+P1(2)
       IF(I31.EQ.2)P1(2)=360.-P1(2)
       IF(I31.EQ.2)P2(2)=90.-P2(2)
       IF(I31.EQ.3)P1(2)=180.+P1(2)


31     CONTINUE
      i_min=100
	if(abs(ac_min(1)).lt.i_min.and.abs(ac_min(1)).gt.0) 
     $i_min=abs(ac_min(1))
	if(abs(ac_min(2)).lt.i_min.and.abs(ac_min(2)).gt.0) 
     $i_min=abs(ac_min(2))
	if(abs(ac_min(3)).lt.i_min.and.abs(ac_min(3)).gt.0) 
     $i_min=abs(ac_min(3))

	do 5501 ii=1,3
5501  ac_min(ii)=nint(float(ac_min(ii))/float(i_min))


	return
	end
c