program eul2eul
c  program to accept Euler angles and output equivalent Euler angles
c  to Bunge Euler angles
c
c   ADR  v 02
c
      real hkl(3),uvw(3),t1(3)
      real a(3,3),d1,d2,d3
      integer i,j,k,iopt,ior,kerr
      character nomen
      logical result
      real qq1(4),quint(4),qend(4)
      integer numsymm,numsamp
c
      common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
      degrad=57.29578
      PI=3.14159265
      twopi=2.0*PI
      nomen='B'
      call textur4              !reads in symmetry operators
      write(*,*) 'Found ',numsymm,' xtal symms',' and '
      write(*,*) numsamp,' sample symms'
      open(unit=2,file='Bunge.angles.txt',status='unknown')
c
      write(*,*) 'To accept Euler (Bunge) angles and',
     &' output Bunge Euler angles'
      write(*,*)
      write(*,*) 'Enter three Bunge angles in degrees'
      read(*,*) d1,d2,d3
c
         call quatB(d1/degrad,d2/degrad,d3/degrad,qq1)
		 do 900, j=1,numsymm ! crystal symmetry
		   call postsymm(qq1,j,quint)
c
			do 880, k=1,numsamp	! sample symmetry
			   call presamp(quint,k,qend)
			   call q2eulB(d1,d2,d3,qend)
			   if(d1.lt.0.0) d1=d1+twopi
			   if(d2.lt.0.0) d2=d2+twopi
			   if(d3.lt.0.0) d3=d3+twopi
			   d1=d1*degrad
			   d2=d2*degrad
			   d3=d3*degrad
c  convert to Bunge Euler angles
c				write(*,*) d1, d2, d3
                           call testangs(d1,d2,d3,result)
                           if(result) then
                              write(*,*) d1, d2, d3
                              write(2,'(3(1x,f8.2)," all<90")') 
     &                          d1, d2, d3
                              else
                              write(2,'(3(1x,f8.2))') 
     &                          d1, d2, d3
                              endif
c                           write(*,*) 'textangs result: ',result
c                           if(result) then
c  this should yield 3 points per orientation
c  thereby enforcing the cubic crystal symmetry within the 90x90x90 space
c
c  DEBUG LINES ->
c	write(*,*)
c	write(*,*) 'Angles for grain number ',i
c	write(*,*) phi1(i),capphi(i),phi2(i)
c                              write(*,*) 'Angles returned from sort '
c                           endif
c
 880                    continue
 900                 continue
c
      call exit
      end
c
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
	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 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 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
	real pi,rad
      character nomen
      save kor
	PI=3.14159265
	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
        if(dph.lt.0.) dph=dph+360.
      elseif(nomen.eq.'B') then
        dps=d1-90.
        dph=90.-d3
        if(dps.lt.0.) dps=dps+360.
        if(dph.lt.0.) dph=dph+360.
      else
        dth=d1
        dph=d2-90.
        if(dph.lt.0.) dph=dph+360.
        dps=90.-d3
        if(dps.lt.0.) dps=dps+360.
      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
	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 testangs(phi1, capphi, phi2,result)
c  tests to see if the angles fall within the specified range
c  which for now is just 0<angle<=90
c
c  LOCAL VARIABLES
	integer i
	real phi1,capphi,phi2
	logical result
        real eps,low_limit,high_limit
        data eps/1.e-4/
        data low_limit /0./
        data high_limit /90./
c
        low_limit=low_limit-eps
        high_limit=high_limit+eps
	result=.true.
	if(phi1.le.low_limit) result=.false.
	if(phi1.ge.high_limit) result=.false.
	if(capphi.le.low_limit) result=.false.
	if(capphi.ge.high_limit) result=.false.
	if(phi2.le.low_limit) result=.false.
	if(phi2.ge.high_limit) result=.false.
	return
	end
c
c ________________________________________
c
c
	subroutine textur4
c
c  second version, for use in REXGBS, WTS2POP
c  not for use in REX2D
c
	integer numsymm,numsamp
c
	common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
	rad=57.29578
c	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(numsymm.gt.24) then
		numsymm=24
		write(*,*) 'too many symmetry operators, cutting back to 24!'
		else
c		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
c	write(*,*) 'reading sample symmetry operators from "quat.symm.ort"'
	open(unit=3,name='quat.symm.ort',status='old')
	read(3,*)	! skip over title line
	read(3,*) numsamp	! read number of sample symmetry operators
	if(numsamp.gt.4) then
		numsamp=4
		write(*,*) 'too many symmetry operators, cutting back to 4!'
		else
c		write(*,*) 'reading ',numsamp,' symmetry operators'
	endif
c
	do 1810, i=1,numsamp
		read(3,*) (quatsamp(ijk,i),ijk=1,4)
c		write(*,*) 'samp. symm. no. ',i,(quatsamp(ijk,i),ijk=1,4)
1810	continue
	close(unit=3)
        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
	double precision dp1,dp,dp2,dq1,dq2,dq3,dq4,sum,diff,pi,dtmp
c
       PI=3.14159265d0
       dq1=qq(1)
       dq2=qq(2)
       dq3=qq(3)
       dq4=qq(4)
c
	if((dabs(dq2).lt.1d-35).and.(dabs(dq1).lt.1d-35)) then
		diff=pi/4.d0
	else
		diff=datan2(dq2,dq1)
	endif
c  ATAN2 is unhappy is both arguments are exactly zero!
c
	if((dabs(dq3).lt.1d-35).and.(dabs(dq4).lt.1d-35)) then
		sum=pi/4.d0
	else
		sum=datan2(dq3,dq4)
	endif
c
	dp1=(diff+sum)
	dp2=(sum-diff)
	dtmp=dsqrt(dq3**2+dq4**2)
	if(dtmp.gt.1.0d0) dtmp=1.0d0
	dp=2.d0*dacos(dtmp)
        p1=dp1
        p=dp
        p2=dp2
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
c
	subroutine presamp(qq,lindex,qresult)
	real qq(4),qresult(4)
	integer lindex
c
	common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
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, and the input quaternion, QQ, is an active
c  rotation, then Qresult = (O x QQ) 
c  the "PRE" in the name refers to writing the operator before the
c  rotation/orientation in vector/tensor notation:  Q' = O x Q
c  or in conventional quaternion notation, Qresult = QQ Ä Qsymm
c
c  Thus, for active rotations (standard definition of orientation in
c  mechanics) this is suitable for applying SAMPLE symmetry
c
	if(lindex.gt.numsamp) stop 'error in PRESAMP, lindex>numsamp'
	if(lindex.lt.1) stop 'error in presamp, lindex<1'
	qresult(1)=qq(1)*quatsamp(4,lindex)+qq(4)*quatsamp(1,lindex)
     &	-qq(2)*quatsamp(3,lindex)+qq(3)*quatsamp(2,lindex)
	qresult(2)=qq(2)*quatsamp(4,lindex)+qq(4)*quatsamp(2,lindex)
     &	-qq(3)*quatsamp(1,lindex)+qq(1)*quatsamp(3,lindex)
	qresult(3)=qq(3)*quatsamp(4,lindex)+qq(4)*quatsamp(3,lindex)
     &	-qq(1)*quatsamp(2,lindex)+qq(2)*quatsamp(1,lindex)
	qresult(4)=qq(4)*quatsamp(4,lindex)-qq(1)*quatsamp(1,lindex)
     &	-qq(2)*quatsamp(2,lindex)-qq(3)*quatsamp(3,lindex)
c
c	write(*,*) 'qresult ',qresult
c
	return
	end
c
c _____________________________
c
c
        subroutine postsymm(qq,lindex,qresult)
        real qq(4),qresult(4)
        integer lindex
c
c        include 'common.f'
        common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion/rotation
c  from applying a symmetry operator, QUATSYMM(n,lindex)
c  to the first quaternion/rotation, QQ
c
c  If the symm oper == O, and the input quaternion, QQ, is an active
c  rotation, then Qresult = (QQ x O) 
c  the "POST" in the name refers to writing the operator after the
c  rotation/orientation in vector/tensor notation:  Q' = Q x O
c  or in conventional quaternion notation, Qresult = Qsymm Ä QQ
c
c  Thus, for active rotations (standard definition of orientation in
c  mechanics) this is suitable for applying CRYSTAL symmetry
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