program mapenergy
c
c  uses energy/angle ratios to iterate towards a trial distribution of energies
c  (excess free energy) of grain boundaries, derived from measurements of
c  triple junction dihedral angles in local equilibrium
c
c  a.d.rollett, copyright 1998, carnegie mellon university
c  version for reading DIHANG.DAT from Adam Morawiec
c
c  this version based on disorientation only for now
c
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
	common /four/ r,nsym
	common /eula/ e
	real r(3,3,48),e(3,3)
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
c  DELTAR is used to index CORRECT, DELTAE is used to index ENERGY
c  E is the 3x3 matrix with misorientation, calculated in EULER
c
	parameter (m=4)
	parameter (n=10000)
	real d1,d2,d3
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2),dom,dth,dph,sig(3)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata,jndex,index,i,j,k,ii
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
	real deltar,error,tend
	real rangle,rtmp,tmp(3)
c  variables for reading in disorientation information
	real dummy,aplane(3,3),bplane(3,3)
c
c
	write(*,*) 'Arrays dimensioned to n= ',n
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
	mresolve=2
c	deltae=tend/float(m-1)
c  cells centered on zero
	deltae=tend/float(m)
c  cells start at zero
	do 5, i=1,n
	do 5, j=1,3
	do 5, k=1,2
5	ratio(i,j,k)=1.0
c
	do 10, i=1,m
	do 10, j=1,m
	do 10, k=1,m
	energy(i,j,k)=1.0
	correct(i,j,k)=1.0
10	continue
c  this is inefficient for cubic systems, but does allow for less symmetric crystal symmetries
c
	call insym
c  gets the crystal symmetry elements for use in INDATA
	call indata
c  fetch the data to be plotted
c
	write(*,*) 'angles   ',((angle(i,j),j=1,3),i=1,5)
c
	write(*,*) 'Enter lower and upper resolution [0,4]'
	write(*,*) 'The resolution is 2^n, so 2^4=16, e.g.  '
	read(*,*) ilow,ihigh
	write(*,*) 'Enter number of times to iterate at each level  [5] '
	read(*,*) iit
	write(*,*) 'Do you want the RATIO calculation? [0=NO]'
	read(*,*) iqq
	do 1000, ijk=ilow,ihigh
c	mresolve=ijk
	mresolve=2**ijk
c  for now, set the resolution to be stepwise increase (cf. powers of 2, possible in 1-parameter system)
c
c	deltar=tend/float(mresolve-1)
c  cells centered on zero
	deltar=tend/float(mresolve)
c  cells start at zero
c  DELTAR sets the increment in Rodrigues units
	npoints=mresolve**3
	write(*,*) 'number of points =',npoints
	write(*,*) 'compare with no. of data points, =', ndata
	write(*,*) 'R vector delta= ',deltar
c	if(ndata.lt.npoints) stop 'not enough data!'
c
	do 1000, kk=1,iit
c  now iterate at this resolution
	if(iqq.ne.0) then
		call ratiocalc
c  calculates RATIOs
		call errorcalc
c  calculates the error
		write(*,*) 'ERROR at iteration ',ijk,':',kk,'  resolution= ',mresolve
		write(*,*) ' =  ',error
		call corrcalc(kk)
c  calculates the correction array (same resolution as energy distribution)
		call energyfix
c   fixes up the energy distribution (applies the corrections to the distribution)
	endif
	call rodplot(0)
c  plot the misorientations with their ratios
	if(iqq.ne.0) then
		call rodplot(2)
c  plot the energy distribution
		call rodplot(3)
c  plot the correction array
c  plot various arrays in postscript files for viewing
		call enerout
c  writes out the current version of the energy distribution
	endif
	call fibout
c  writes out particular zone axes in the energy distribution
1000	continue
c
	stop 'done!'
	end
c
c  __________________________
c
	subroutine insym
c
	common /four/ r,nsym
	common /eula/ e
	real r(3,3,48),e(3,3)
	integer type
	real  pi,pi2,sqr2
c  AX used for points along the arc at the RH edge of the unit triangle
      logical repeat
c
c**** ******************************************************************
c
c**** read input files:
c
c**** ******************************************************************
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
    8 format(a)
c    8 format(a\)
      write(*,*) 'opening files'
c
      inquire(file='crysym',exist=repeat)
      if(repeat) then
        open(1,file='crysym',status='old')
      else
  131   write(*,8)' Enter the filename for CRYstal SYMmetries: '
        read 1,fname
    1 format(2a)
    2 format(i1)
        open(1,file=fname,status='old',err=131)
      endif
      read(1,*) nsym
      print *,'no.of crystal symmetry operators = ',nsym
      do 10 i=1,nsym
      do 10 j=1,3
      read(1,*) (r(j,k,i),k=1,3)
   10 continue
c *** above: crystal symmetries in R()
      close(1)
	do 11, i=nsym+1,nsym*2
	do 11, j=1,3
	do 11, k=1,3
	r(j,k,i)=-1.*r(j,k,(i-nsym))
11	continue
c	nsym=nsym*2
	write(*,*) ' only interested in proper rotations, so use 1st 24 '
      print *,'no.of crystal symmetry operators w/inversion = ',nsym
c  add the inversion center!
c
	return
	end
c
c  __________________________
c
	subroutine indata
c
	parameter (m=4)
	parameter (n=10000)
c
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
	common /four/ r,nsym
	common /eula/ e
	real energy(m,m,m),correct(m,m,m)
	real r(3,3,48),e(3,3)
	integer type,kindex(3)
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2),dom,dth,dph,sig(3)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata,jndex,index,i,j,k,ii
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
	real deltar,error,tend
	real rodr(4,3),thetmin(3),tmp(3)
	real d1,d2,d3
	integer iopt,ior,kerr
	character nomen
c
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
c
c
	real v100(3),v101(3),v110(3),tmpmis(3,3),detm(3)
	real misor(3,3,3),fmat(3,3,3)
	real aplane(3),bplane(3)
	real aeul(3,3),aaeul(3,3,3)
c  AEUL is for reading in 3 sets of euler angles from saylor files
c  AAEUL holds the 3x3 matrices from the euler angles
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	sqr2=sqrt(2.)
 	sqr3=sqrt(3.)
	nomen='B'
c  bunge angles assumed for both Morawiec and Saylor inputs
	iopt=2
c  calcs. 3x3 from euler angles
	ior=0
	kerr=0
	v100(1)=0.
	v100(2)=1.
	v100(3)=0.
	v101(1)=-1.*sqr2
	v101(2)=0.
	v101(3)=sqr2
	v110(1)=sqr2
	v110(2)=-1.*sqr2
	v110(3)=0.
c
	write(*,*) 'Is this a Morawiec file [=0],'
     &	,'a Yang file [=1], or a Saylor file [=2]?'
	read(*,*) iqans
	if(iqans.eq.0) open(unit=2,file='dihang.dat',status='old')
	if(iqans.eq.1) then
	    open(unit=2,file='energy.in',status='old')
	    read(2,*)
	elseif(iqans.eq.2) then
	    open(unit=2,file='COMBINE.TXT',status='old')
	    read(2,*)
	    endif
	if(iqans.lt.0.or.iqans.gt.2) stop 'invalid repsonse!'
	open(unit=1,file='check.out',status='unknown')
c	read(2,*) 
c  skip title line
	write(*,*) 'How many triple junctions to read in [max=10,000]?'
	read(*,*) iqtop
c
	do 20, i=1,iqtop
c	write(*,*) 'TJ set no. ',i
	   if(iqans.eq.0) read(2,*,end=30) inum
	   if(iqans.eq.2) read(2,*,end=30) inum
c  skip over line with number in it for morawiec or saylor files
	   if(i.lt.10) write(1,*) ' triple junction no.  ',inum
c
	if(iqans.eq.2) then
		do 40, ijk=1,3
		 read(2,*,end=30) angle(i,ijk),(aeul(ijk,ik),ik=1,3)
	if(i.lt.10) write(*,*) 'input= ',angle(i,ijk),(aeul(ijk,ik),ik=1,3)
		call euler(iopt,nomen,aeul(ijk,1),aeul(ijk,2),aeul(ijk,3),ior,kerr)
		do 38, ik=1,3
		do 38, jk=1,3
38		aaeul(ijk,ik,jk)=e(jk,ik)
c  take the transpose of E because of Canova convention to go from
c  crystal to sample, as passive transformation of axes
c
40		continue
		endif
c
	do 580, k=1,3
	if(iqans.eq.0) then
c  first let's do the Adam M style files
	   read(2,*,end=30) angle(i,k),d1,d2,d3
575	format(20f12.4)
	if(i.lt.10) write(1,*) 'input line:  ',angle(i,k),d1,d2,d3
c  read in the disorientation angle (degrees), true dihedral angle, disorientation matrix (3x3)
c
	  call euler(iopt,nomen,d1,d2,d3,ior,kerr)
c  now for the CCY style files
	  elseif(iqans.eq.1) then
	read(2,*,end=30) disorient(i,k,5),angle(i,k),disorient(i,k,1)
     &		,(bplane(jj),jj=1,3),(aplane(jj),jj=1,3)
     &		,((e(ii,jj),jj=1,3),ii=1,3)
c	      read(2,575,end=30) disorient(i,k,5),angle(i,k),disorient(i,k,1)
c     &		,(aplane(jj),jj=1,3),(bplane(jj),jj=1,3)
c     &		,((e(ii,jj),jj=1,3),ii=1,3),dummy
c  DISORIENT(i,k,5) contains the minimum angle, DISORIENT(i,k,1) contains the angle associated with 
c the 3x3 matrix
c  Aplane and Bplane are already in crystal axes
c
	elseif(iqans.eq.2) then
c  now we deal with the saylor files (we already read in the angles)
		ij=k+1
		if(ij.gt.3) ij=ij-3
		jj=k+2
		if(jj.gt.3) jj=jj-3
		do 65, ik=1,3
		do 65, jk=1,3
		e(ik,jk)=0.
		do 65, kk=1,3
		e(ik,jk)=e(ik,jk)+aaeul(ij,ik,kk)*aaeul(jj,jk,kk)
65		continue
c  obtain misorientation matrix (not disorientation) for opposite boundary
c
	endif
c
	do 577, ii=1,3
	do 577, jj=1,3
577	misor(ii,jj,k)=e(ii,jj)
c
	if(i.lt.10) write(1,*) 'dihedral, misorientation',angle(i,k)
     &		,((misor(ii,jj,k),ii=1,3),jj=1,3)
c
c  have read a set of 3 boundaries+dihedral angles
c
580	continue
c  now we have three sets of boundary info
c
	do 620, ijk=1,3
	thetmin(ijk)=360.
620	continue
c
c	write(*,*) 'about to start 1800 loop'
	do 1800, isort=1,nsym
c
c  apply symm element to MISOR to get FMAT
c
	do 630, ijk=1,3
	call symm_apply(misor,fmat,ijk,isort)
630	continue
c
c	if(i.lt.10) write(1,*) 'FMAT= ',(((fmat(ii,jj,kk),ii=1,3),jj=1,3),kk=1,3)
c
c  now to convert misor to rodrigues vectors for all 3 boundaries at the TJ
	call mat2rod(fmat,rodr,detm)
c
	do 640, ijk=1,3
	if((rodr(4,ijk)*degrad).lt.thetmin(ijk)) then
		thetmin(ijk)=rodr(4,ijk)*degrad
		kindex(ijk)=isort
		endif
c  identify minimum angle and index of symmetry element
640	continue
c
1800	continue
c  end of loop to identify min. misorient.
c
	if(i.lt.10) then 
	write(1,*) ' min angles, indexes'
	write(1,*) (thetmin(ijk),ijk=1,3)
	write(1,*) (kindex(ijk),ijk=1,3)
	endif
c
	do 680, ijk=1,3
c
	call symm_apply(misor,fmat,ijk,kindex(ijk))
c
680	continue
c
c  now to convert misor to rodrigues vectors for all 3 boundaries at the TJ
	call mat2rod(fmat,rodr,detm)
c
	do 3010, k=1,3
	do 3008, ii=1,3
	disorient(i,k,ii+1)=rodr(ii,k)
3008	continue
	disorient(i,k,1)=rodr(4,k)*degrad
c  this is the misorientation angle, now in degrees in DISORIENT
c
	tmp(1)=0.
	imax=0
	tmp(3)=1.
	imin=0
	do 194, ijk=1,3
		if(abs(disorient(i,k,1+ijk)).gt.tmp(1)) then
			tmp(1)=abs(disorient(i,k,1+ijk))
			imax=ijk
		endif
		if(abs(disorient(i,k,1+ijk)).lt.tmp(3)) then
			tmp(3)=abs(disorient(i,k,1+ijk))
			imin=ijk
		endif
194	continue
	do 196, ijk=1,3
		if(ijk.ne.imax.and.ijk.ne.imin) tmp(2)=abs(disorient(i,k,1+ijk))
196	continue
	do 198, ijk=1,3
		disorient(i,k,1+ijk)=tmp(ijk)
198	continue
	if(i.lt.10) write(1,*) 'Rodrigues vector=',(disorient(i,k,ijk),ijk=1,4)
c  all above sorts the components of the Rodrigues vector in descending order
c  equivalent to using the xtal symms on the other xtal!  see Moraviec article
c
3010	continue
c
20	continue
	write(*,*) ' used up all ',iqtop,' points!'
30	ndata=i-1
	write(*,*) ' number of actual triple junctions =',ndata
c
c
	write(*,*) ' done with input!'
	close(unit=1)
	close(unit=2)
c
	open(unit=3,file='mape.angles.out',status='unknown')
c  we will use this file to write out the angles from the input analysis
c
	write(3,*)  '  TJno  GBno  dihedral  disorient   R1  R2  R3'
	do 710, i=1,ndata
	do 710, k=1,3
		write(3,700) i,k,angle(i,k),(disorient(i,k,ijk),ijk=1,4)
700	format(2i6,6f10.4)
710		continue
c
	close(unit=3)
	return
	end
c
c  ____________________________
c
	subroutine symm_apply(eee,f,ijk,index)
	common /four/ r,nsym
	real r(3,3,48)
	real eee(3,3,3),f(3,3,3)
	integer ijk,index,ia,ib,ic
c  matrix multiply, for symmetry operations
c
	do 1940, ia=1,3
	do 1940, ib=1,3
	f(ia,ib,ijk)=0.
	do 1940, ic=1,3
	f(ia,ib,ijk)=f(ia,ib,ijk)+eee(ia,ic,ijk)*r(ib,ic,index)
c  postmultiply by transpose of the symmetry operator (c1)
c  for A crystal since E goes from 1st (A) to 2nd (B) crystal
1940	continue
c
	return
	end
c
c  ____________________________
c
	subroutine mat2rod(e,rodr,detm)
	real e(3,3,3),rodr(4,3),tmp(3),rnorm
	real det,detm(3)
c
c  converts a MISOR matrix to Rodrigues vector, 3 at a time
c
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
	do 100, i=1,3
c  calc. determinant
c
	det=e(1,1,i)*e(2,2,i)*e(3,3,i)
     &	+e(2,1,i)*e(3,2,i)*e(1,3,i)
     &	+e(3,1,i)*e(1,2,i)*e(2,3,i)
     &	-e(1,1,i)*e(3,2,i)*e(2,3,i)
     &	-e(1,2,i)*e(2,1,i)*e(3,3,i)
     &	-e(1,3,i)*e(2,2,i)*e(3,1,i)
	if(det.lt.-0.95) det=-1.0
	if(det.gt.0.95) det=1.0
	detm(i)=det
c
c  we are only interested in proper vs. improper rotations here
c
c	write(1,*) 'rotation matrixes:  ',(e(ii,1,i),ii=1,3)
c	write(1,*) 'rotation matrixes:  ',(e(ii,2,i),ii=1,3)
c	write(1,*) 'rotation matrixes:  ',(e(ii,3,i),ii=1,3)
c	write(1,*) 'index,  determinant= ',i,det
	rtmp=(((e(1,1,i)+e(2,2,i)+e(3,3,i))/det)-1.)/2.
	if(rtmp.lt.-1.) rtmp=-1.
	if(rtmp.gt.1.) rtmp=1.
	rodr(4,i)=acos(rtmp)
c  magnitude of rotation
	tmp(1)=(e(2,3,i)-e(3,2,i))/det
	tmp(2)=(e(1,3,i)-e(3,1,i))/det
	tmp(3)=(e(1,2,i)-e(2,1,i))/det
	rnorm=sqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
	if(rnorm.lt.1e-5) then
		tmp(1)=sqrt(1.+e(1,1,i))
		tmp(2)=sqrt(1.+e(2,2,i))
		tmp(3)=sqrt(1.+e(3,3,i))
		rnorm=sqrt(tmp(1)**2+tmp(2)**2+tmp(3)**2)
	endif
c	write(1,*) rnorm, tmp(1),tmp(2),tmp(3),tan(rodr(4,i)/2)
	rodr(1,i)=tmp(1)/rnorm*tan(rodr(4,i)/2.)
	rodr(2,i)=tmp(2)/rnorm*tan(rodr(4,i)/2.)
	rodr(3,i)=tmp(3)/rnorm*tan(rodr(4,i)/2.)
c  components of the Rodrigues vector
c	write(1,190)  rodr(4,i)*180./pi
c     &	,rodr(1,i),rodr(2,i),rodr(3,i)
c190	format(' angle, vector....',4f9.3)
c
100	continue
	return
	end
c
c  _____________________________
c
	subroutine ratiocalc
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
c
	integer mresolve,ndata
	integer jcntnum(0:10),num10
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
	integer index(3),jcount
	real etmp1,etmp2
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	num10=0
c  record number of ratios>10
c
	do 100, i=1,ndata
	do 50, j=1,3
	   jcount=0
	do 50, k=1,3
	if(j.eq.k) goto 50
	jcount=jcount+1
c
c  note the use of DELTAE here, so that all the array is used but low resolution means that many
c  cells have the same entry
c
	do 45,l=1,3
c	index(l)=1+int(0.5+disorient(i,j,(1+l))/deltae)
	index(l)=1+int(disorient(i,j,(1+l))/deltae)
	if(index(l).gt.m) index(l)=m
	if(index(l).lt.1) index(l)=1
45	continue
	etmp1=energy(index(1),index(2),index(3))
	if(etmp1.le.0.) etmp1=1.0
	do 46,l=1,3
c	index(l)=1+int(0.5+disorient(i,k,(1+l))/deltae)
	index(l)=1+int(disorient(i,k,(1+l))/deltae)
	if(index(l).gt.m) index(l)=m
	if(index(l).lt.1) index(l)=1
46	continue
	etmp2=energy(index(1),index(2),index(3))
	if(etmp2.le.0.) etmp2=1.0
	ratio(i,j,jcount)=(etmp1/etmp2)*(sin(angle(i,k)/degrad)
	1    /sin(angle(i,j)/degrad))
	if(ratio(i,j,jcount).gt.10.) then
c	    write(*,*) 'RATIOCALC:index... ',index(1),index(2),index(3)
c	    write(*,*) 'angles ',angle(i,k),angle(i,j)
c	    write(*,*) 'i,j,k,ratiojl... ',i,j,jcount,ratio(i,j,jcount)
	    ratio(i,j,jcount)=10.
	    num10=num10+1
	endif
c
50	continue
100	continue
	write(*,*) ' (Warning) Number of ratios > 10.0 = ',num10
	return
	end
c
c  _________________________
c
	subroutine errorcalc
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
c
	integer mresolve,ndata
	real error
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
	do 100, i=1,ndata
	do 100, j=1,3
	do 100, k=1,2
	error=error+abs(alog10(ratio(i,j,k)))
100	continue
	error=error/ndata/3./2.
	return
	end
c
c  _____________________________
c
	subroutine corrcalc(kiter)
c  calculates the averages of the ratios
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
c
	integer mresolve,ndata
	integer ix(3),jcount(m,m,m)
c  JCOUNT is used to count how many ratios were included in each entry in CORRECT
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
	real dummy
	integer kiter
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	do 10, i=1,m
	do 10, j=1,m
	do 10, k=1,m
	correct(i,j,k)=0.0
	jcount(i,j,k)=0
10	continue
c
	do 11, i=0,10
	jcntnum(i)=0
 11	continue
c
c  now we bin the RATIOs and calculate a product (Pi(Ratio[i]))
	do 100, i=1,ndata
	do 50, j=1,3
	do 45,l=1,3
c	ix(l)=1+int(0.5+disorient(i,j,(1+l))/deltar)
	ix(l)=1+int(disorient(i,j,(1+l))/deltar)
c  use of DELTAR means that we are using the coarse grid to bin ratios
	if(ix(l).gt.m) ix(l)=m
	if(ix(l).lt.1) ix(l)=1
45	continue
c	write(*,*) 'i,ix1,ix2,ix3  ',i,ix(1),ix(2),ix(3)
c
c	correct(ix(1),ix(2),ix(3))=
c	1    correct(ix(1),ix(2),ix(3))*ratio(i,j,1)*ratio(i,j,2)
c  geometric mean
	correct(ix(1),ix(2),ix(3))=
	1    correct(ix(1),ix(2),ix(3))+log(ratio(i,j,1))+log(ratio(i,j,2))
c	correct(ix(1),ix(2),ix(3))=
c	1    correct(ix(1),ix(2),ix(3))+ratio(i,j,1)
c	correct(ix(1),ix(2),ix(3))=
c	1    correct(ix(1),ix(2),ix(3))+ratio(i,j,2)
c  arithmetic mean
c
c  accumulate the ratios in the correct array
c  and count how many are there
	jcount(ix(1),ix(2),ix(3))=jcount(ix(1),ix(2),ix(3))+2
c
50	continue
100	continue
c
c  now we take the products of ratios and take a 1/N power to get the average
c  note the inclusion of a factor of 2 to reduce the correction in each cycle
c
	do 200, k=1,mresolve
	do 200, i=k,mresolve
	do 200, j=k,i
c	write(*,*) 'i,j,k,correct(i,j,k)  ',i,j,k,correct(i,j,k)
c	write(*,*) 'jcount  ',jcount(i,j,k)
	if(jcount(i,j,k).gt.0) then
c	correct(i,j,k)=(exp(correct(i,j,k))/
c	1	float(jcount(i,j,k)))
c  geometric mean
c
	    correct(i,j,k)=exp((correct(i,j,k))/
	1	float(jcount(i,j,k)))
c(	1	float(jcount(i,j,k))/float(kiter))
c	correct(i,j,k)=(correct(i,j,k)/
c	1	float(jcount(i,j,k)))
	    else
	        correct(i,j,k)=1.0
	endif
200	continue
c
c	do 300, k=1,mresolve
c	do 300, i=k,mresolve
c	do 300, j=k,mresolve
c	write(*,290) i,k,(jcount(i,j,k),j=k,i)
c 290	   format( 'jcount(i,',2(i4,x),')= ',32i6)
c300	continue
c
	do 400, k=1,mresolve
	r3=tend*((float(k)-0.5)/float(mresolve))
	do 400, i=k,mresolve
	r1=tend*((float(i)-0.5)/float(mresolve))
	do 400, j=k,i
	r2=tend*((float(j)-0.5)/float(mresolve))
	if(r1+r2+r3.gt.1.0) goto 400
c  skip if we are beyond the triangular cutoff in Rodrigues space
	index=jcount(i,j,k)/2
c  pairs of ratios, therefore divide by 2
	if(index.lt.0) goto 400
c  back-stop: there should not be any negative counts!
	if(index.gt.10) index=10
	    jcntnum(index)=jcntnum(index)+1
 400	continue
c
	return
	end
c
c  ________________________
c
	subroutine energyfix
c  calculates the averages of the ratios
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata,ix,jx,kx,jcount
	real enorm
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
	enorm=0.0
	jcount=0
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	do 10, k=1,m
	do 10, i=k,m
	do 10, j=k,i
	r3=tend*((float(k)-0.5)/float(m))
	r2=tend*((float(j)-0.5)/float(m))
	r1=tend*((float(i)-0.5)/float(m))
	if(r1+r2+r3.gt.1.0) goto 10
c  skip if we are beyond the triangular cutoff in Rodrigues space
c  this is not important for the correction of the energy
c  but IS important to the normalization which should only be performed
c inside the fundamental zone!
c
c	ix=1+int(0.5+float((mresolve-1)*(i-1))/float(m-1))
	ix=1+int((float(mresolve)*(float(i)-0.5))/float(m))
		if(ix.gt.m) ix=mresolve
		if(ix.lt.1) ix=1
c	jx=1+int(0.5+float((mresolve-1)*(j-1))/float(m-1))
	jx=1+int((float(mresolve)*(float(j)-0.5))/float(m))
		if(jx.gt.m) jx=mresolve
		if(jx.lt.1) jx=1
c	kx=1+int(0.5+float((mresolve-1)*(k-1))/float(m-1))
	kx=1+int((float(mresolve)*(float(k)-0.5))/float(m))
c  cells start at zero
		if(kx.gt.mresolve) kx=mresolve
		if(kx.lt.1) kx=1
	if(correct(ix,jx,kx).ne.0.) then
	    energy(i,j,k)=energy(i,j,k)/correct(ix,jx,kx)
	    else
	    write(*,*) 'i,j,k,ix,jx,kx..  ',i,j,k,ix,jx,kx
	    write(*,*) energy(i,j,k),correct(ix,jx,kx)
	endif
c  if the average ratio is > 1, reduce the energy at that location
	enorm=enorm+energy(i,j,k)
	jcount=jcount+1
c	correct(i,j,k)=1.0
c  restore the correction array: not now that we use different resolutions
c	write(*,*) 'energyfix: energy(',i,j,k,')=',energy(i,j,k) 
10	continue
c
	enorm=enorm/float(jcount)
c
	do 100, k=1,m
	do 100, i=k,m
	do 100, j=k,i
	energy(i,j,k)=energy(i,j,k)/enorm
100	continue
c
	return
	end
c
c  ________________________
c
	subroutine enerout
c  writes out the energy array to a file
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
	character*1 ten,one
	character blankout*80,blank*80
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	blank='                                                            '//
     &	'                                               '
c
	ten=char(48+(mresolve/10))
	one=char(48+mod(mresolve,10))
	open(unit=3,file='energy.out.'//ten//one,status='unknown')
	write(3,*) m,'            energy.out'//ten//one
	do 10, k=1,m
	do 10, j=1,m
	do 10, i=k,m
c	write(*,*) 'energy(',i,j,k,')=',energy(i,j,k) 
	write(3,20) i,j,k,  energy(i,j,k)
10	continue
20	format(3i8,2f12.4)
	close(unit=3)
c
	open(unit=4,file='energy.out',status='unknown',recl=132)
	write(4,*) ' energy.out     M=',m,'    error=',error
	do 25, k=m,1,-1
	write(4,*) 'energy(i,j[1..M],',k,')='
     &	     , '  R3= ',(float(k-1)+0.5)*tend/float(m)
	do 25, i=m,k,-1
c	itmp=max0(1,8*(k-1))
c	blankout=blank(1:itmp)
c
	write(4,24) (energy(j,i,k) ,j=i,m)
24	format(12f8.3)
25	continue
	close(unit=3)
c
	open(unit=3,file='correct.out.'//ten//one,status='unknown')
	do 30, k=1,mresolve
	do 30, i=k,mresolve
	write(3,*) 'correct(',i,',',j,'k,)=  '
	write(3,40) (correct(i,j,k) ,j=k,i)
30	continue
40	format(10f8.3/)
	close(unit=3)
c
	open(unit=3,file='jcntnum.out.'//ten//one,status='unknown')
	write(3,*) '  number  count(ratios)  '
	write(*,*) '  number  count(ratios)  '
	do 50, k=0,10
	write(3,45) k,jcntnum(k)
	write(*,45) k,jcntnum(k)
 45	format(x,2i9)
 50	continue
	close(unit=3) 
	write(*,*)
c
	return
	end
c
c  ________________________
c
	subroutine fibout
c  writes out certain zone axes of the energy array to a file
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
	character*1 ten,one
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	ten=char(48+(mresolve/10))
	one=char(48+mod(mresolve,10))
	open(unit=3,file='fib.100.out.'//ten//one,status='unknown')
	write(3,*) '  num  R1    R2  R3  gamma(100)'
	do 10, k=1,m
	write(3,20) k,(float(k)-0.5)*tend/float(m),0.,0.,energy(k,1,1) 
10	continue
20	format(2x,i5,4(2x,f8.3))
	close(unit=3)
c
	open(unit=3,file='fib.110.out.'//ten//one,status='unknown')
	write(3,*) '  num  R1    R2  R3      gamma(110)'
	do 30, k=1,m
	rtmp=(float(k)-0.5)*tend/float(m)
	write(3,20) k,rtmp,rtmp,0.,energy(k,k,1) 
30	continue
	close(unit=3)
c
	open(unit=3,file='fib.111.out.'//ten//one,status='unknown')
	write(3,*) '  num  R1    R2  R3      gamma(111)'
	do 50, k=1,m
	rtmp=(float(k)-0.5)*tend/float(m)
	write(3,20) k,rtmp,rtmp,rtmp,energy(k,k,k) 
50	continue
	close(unit=3)
c
	return
	end
c
c  __________________________________
c
	subroutine rodplot(type)
c
c  plots sets of Rodrigues triples in 7 triangles
c  output= postscript file
c
	integer type
	common /one/ energy,tj,angle,disorient,ratio,mresolve,deltar,correct
	common /two/ ndata,pi,pi2,degrad,deltae,error
	common/three/jcntnum,jcntfreq
c  MRESOLVE sets the resolution, NDATA is the number of exptl TJs
	parameter (m=4)
	parameter (n=10000)
	real energy(m,m,m),correct(m,m,m)
c  the gb energy is a function of R1,R2,R3: norm(energy)=1
c  CORRECT contains the binned & averaged RATIOS for correcting the energy distribution
c
	real tj(n)
	real angle(n,3)
c  3 dihedral angles per TJ
	real disorient(n,3,5)
c  describe the disorientation as Theta ,R1,R2,R3
	real ratio(n,3,2)
c  each angle has 2 (reciprocally related) energy/angle ratios associated
	integer mresolve,ndata
	integer jcntnum(0:10)
	real jcntfreq(20)
c  JCNTNUM contains number of cells with "index" counts; zero counts is significant!
c  JCNTFREQ contains frequency data normalized by mean number
c
	character*1 ten,one,ctype
c
	real   tx(10),ty(10)
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      degrad=180./pi
c
	tx(1)=1
	ty(1)=1
	tx(2)=1
	ty(2)=3.1
	tx(3)=1
	ty(3)=5.2
	tx(4)=1
	ty(4)=7.3
	tx(5)=3.5
	ty(5)=3.1
	tx(6)=3.5
	ty(6)=5.2
	tx(7)=3.5
	ty(7)=7.3
c  7 triangles, ll corners
c
      pi=4.*atan(1.)
c
	ten=char(48+(mresolve/10))
	one=char(48+mod(mresolve,10))
	ctype=char(48+type)
	open(unit=9,file='m2.'//ctype//'.'//ten//one//'.ps',status='unknown')
c
	write(9,4010) 
4010	format('%!PS-Adobe-3.0')
	write(9,4020) 
4020	format('%%BoundingBox: 36 36 599 792')
c  box for .5" margins on a 8.5x11 page
	write(9,4030) 
4030	format('%%Creator:mape.f')
	write(9,4040) 
4040	format('%%EndComments')
c
	write(9,5010) 
5010	format('/rod2pix {347.646753 mul} def')
c  transforms RF units (0.41 = 45 degr about 100) to 2" or 144 pixels in PS space
	write(9,5020) 
5020	format('/inch {72 mul} def')
c  transforms inches to 72 pixels in PS space
	tri_size=2.0
c  sets size of triangle
	write(9,5030) tri_size,tri_size,-1*tri_size,-1*tri_size
5030	format('/triangle    %stack: x y == ll corner'/
     &	'{ newpath moveto'/
     &	' ',f8.3,' inch 0 rlineto'/
     &	' 0 ',f8.3,' inch rlineto'/
     &	' ',f8.3,' inch ',f8.3,' inch rlineto'/
     &	' closepath } def')
c  defines a triangle of size 2 inches, for plotting within
	box_size=0.2
	write(9,5040) box_size
5040	format('/bsz ',f8.3,' def')
	write(9,5050) 
5050	format('/box    %stack: x y == corner'/
     &	'{ newpath moveto'/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	write(9,5060) 
5060	format('/rbox    %draws a box wherever we are!'/
     &	'{ '/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	WRITE(9,4050) 
4050	format('/Times-Roman findfont 18 scalefont setfont')
	write(9,4090) 
4090	format('%%EndProlog')
	rfincr=1/3./7.
c  increment of R3 between sections
	if(type.eq.0) call pltpts(disorient,ratio,ndata,rfincr,tx,ty,tri_size)
c	if(type.eq.1) call pltpts2()
c  not yet implemented
	if(type.eq.2) call pltdist(energy,m,deltae,tx,ty,tri_size,m)
	if(type.eq.3) call pltdist(correct,m,deltar,tx,ty,tri_size,mresolve)
	do 100, i=1,7
	rflabel=float(i-1)*rfincr
c  current R3 value
	rtmp=tri_size*rfincr
	write(9,7010) tx(i),ty(i)
7010	format(f8.3,' inch ',f8.3,' inch triangle '/
     &	'2 setlinewidth stroke ')
	write(9,7020) tx(i)+0.2,ty(i)+1.6,rflabel
7020	format(f8.3,' inch ',f8.3,' inch moveto (R3 =',f6.3,') show')
c  should put the label above and to the right of the ll corner
c
	write(9,7030) tx(i),ty(i)
7030	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7035) rflabel,rflabel
7035	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7040) tend-rflabel
7040	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [] 0 setdash 2 setlinewidth stroke')
c  draw dashed line to the edge
	write(9,7050) tx(i),ty(i)
7050	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7060) rflabel+rfincr,rflabel+rfincr
7060	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7070) tend-rflabel-rfincr
7070	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [2] 1 setdash 2 setlinewidth stroke [] 0 setdash')
c  draw solid line to the edge
c
	if(rflabel.gt.0.1715) then
c  we need the cut at the top right corner
	write(9,7110) tx(i),ty(i)
7110	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-rflabel-tend
c   based on the cutoff plane being defined by X+Y+Z=1
	write(9,7115) tend,rod2
7115	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-rflabel)/2.
	write(9,7120) (rod2p-tend),(rod2p-rod2)
7120	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	     ' [] 0 setdash 2 setlinewidth stroke')
c  draw dashed line to the diagonal
	write(9,7130) tx(i),ty(i)
7130	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-(rflabel+rfincr)-tend
	write(9,7135) tend,rod2
7135	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-(rflabel+rfincr))/2.
	write(9,7140) rod2p-tend,rod2p-rod2
7140	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	    ' [2] 1 setdash 2 setlinewidth stroke')
c  draw solid line to the diagonal
	endif
c    draws lines in upper right corner of triangle
c	  
100	continue
c
	write(9,8010) tx(1)+2.7,ty(1)+.5,type
8010	format(f8.3,' inch ',f8.3,' inch moveto (plot type = ',i3,') show')
	write(9,8020) tx(1)+2.7,ty(1)+.9,mresolve
8020	format(f8.3,' inch ',f8.3,' inch moveto (resol. = ',i4,') show')
	write(9,8030) tx(1)+2.7,ty(1)+1.3,error
8030	format(f8.3,' inch ',f8.3,' inch moveto (error = ',f10.6,') show')
	write(9,9010) 
9010	format('showpage  '/'%%PageEnd'/'%%EOF')
	close(9)
	return
	end
c
c ______________________________
c
	subroutine pltpts(angle,ratio,ndata,rfincr,tx,ty,tri_size)
	parameter (n=10000)
	real tx(10),ty(10)
	real angle(n,3,5),ratio(n,3,2),rfincr,tri_size
	integer ndata
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
c
	write(9,*) '/bsz 0.05 def'
c  sets boxsize to 0.1 to be certain
c
	scale=0.1
	do 100, i=1,ndata
	do 100, j=1,3
	ntriang=1+int(0.5+angle(i,j,4)/rfincr)
	if(ntriang.lt.1) ntriang=1
	if(ntriang.gt.7) ntriang=7
	write(9,6010) tx(ntriang),ty(ntriang)
6010	format('newpath ',f8.3,' inch ',f8.3,' inch moveto ')
c  puts us at the ll corner of the triangle
	write(9,6020) angle(i,j,2),angle(i,j,3)
6020	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto rbox')
	tmp=alog10(sqrt(ratio(i,j,1)*ratio(i,j,2)))
c  for now, assume that we are interested in a scale of 0.5 to 2 in ratios
	if(tmp.gt.scale) tmp=scale
	if(tmp.lt.-1*scale) tmp=-1*scale
	thue=scale-tmp
c  colour (hue) ranges from about red to blue
	tbright=4.*abs(tmp)/scale
	if(tbright.gt.1.0) tbright=1.0
	write(9,6030) thue,1.0,tbright
c  arrange for bright colors unless close to ratio=1
6030	format(' gsave ',3f8.3,' sethsbcolor fill grestore ')
c  outline a square symbol, set color and grey level, fill
100	continue
c
	write(9,*) '/bsz 0.5 def'
c  sets boxsize to 0.5" for legend
c
	do 200, i=-4,4
	rtmp=(i+4)*0.5+0.25
	tmp=scale*float(i)/4.
	if(mod(i,2).eq.0) write(9,8010) 
     &	    tx(4)+rtmp-0.25,ty(4)+2.15,exp(tmp*2.3026)
8010	format(f8.3,' inch ',f8.3,' inch moveto (',f6.2,') show')
	write(9,8020) tx(4)+rtmp,ty(4)+2.6
8020	format(f8.3,' inch ',f8.3,' inch box')
	thue=scale-tmp
c  colour (hue) ranges from about red to blue
	tbright=4.*abs(tmp)/scale
	if(tbright.gt.1.0) tbright=1.0
	write(9,6030) thue,1.0,tbright
c  arrange for bright colors unless close to ratio=1
c  outline a square symbol, set color and grey level, fill

200	continue
	return
	end
c
c  _____________________________
c
	subroutine pltdist(array,m,delta,tx,ty,tri_size,mres)
	real array(m,m,m),delta,rod1,rod2,thue,tbright
	real tx(10),ty(10)
	integer m,mres
c
c
	tend=sqrt(2.)-1.
c  defines edge of the R-F space
c
      pi=4.*atan(1.)
      pi2=0.5*pi
      degrad=180./pi
	cell=tend/float(mres-1)
	deltax=tri_size/float(mres)
	do 150,  k=1,7
c  k indexes R3 by taking a section through RF space
	kindex=1+int(float(k-1)*float(mres-1)/3./tend/7.)
c  using 32-3=29 as the upper limit gives a few more points in each section
c  especially at the upper end
c
	write(9,*) ' gsave'
	write(9,5010) tx(k),ty(k)
5010	format(f8.3,' inch ',f8.3,' inch triangle clip')
c  clip the triangle to keep the plot clean since we are using filled boxes
c
	write(9,*) '/bsz ',deltax,'  def'
c	jindex=max0(kindex-int(float(mres)/14.),1)
	jindex=kindex
	do 100, i=jindex,mres
c  i indexes R1
	rod1=float(i-1)*deltax
c  calculate the x coordinate based on a triangle with 2" side
	do 100, j=jindex,i
c  j indexes R2
	rod2=float(j-1)*deltax
	write(9,6020) tx(k)+rod1,ty(k)+rod2
6020	format(f8.3,' inch ',f8.3,' inch box')
c  boxes are centered on each point, but clipping takes care of spill-over
c
	rztmp=-1.*alog10(array(i,j,kindex))
	scale=0.1
c  for now, assume that we are interested in a scale of 0.5 to 2 in energies
c
	rtmp=rztmp/scale
	tmp=amin1(1.0,rtmp)
	tmp=amax1(-1.0,rtmp)
	thue=((tmp+1.)/2.9)
c  colour (hue) ranges from about red to blue
c	tbright=(tmp+1.)/2.
	tsat=amax1(0.1,(2.*abs(tmp)/scale))
c	tbright=4.*abs(tmp)/scale
	if(tsat.gt.1.0) tsat=1.0
c
	write(9,6030) thue,tsat,1.0
6030	format(' gsave ',3f8.3,' sethsbcolor fill grestore ')
c  outline a square symbol, set color and grey level, fill
	if(rztmp.gt.0.0) write(9,6040) tx(k)+rod1,ty(k)+rod2
6040	format(' newpath',f8.3,' inch ',f8.3,' inch ',
     &	'0.08 0 360 arc gsave 0 setgray fill grestore')
c   put a small dot in the box if the value is < 0
100	continue
	write(9,*) ' grestore'
150	continue
c
	write(9,*) '/bsz 0.5 def'
c  sets boxsize to 0.5" for legend
c
	do 200, i=-4,4
	rtmp=(4-i)*0.5
	tmp=float(i)/4.
	if(mod(i,2).eq.0) write(9,8010) 
     &	     tx(4)+rtmp,ty(4)+2.15,exp(-1.*scale*tmp*2.3026)
8010	format(f8.3,' inch ',f8.3,' inch moveto (',f6.2,') show')
	write(9,6020) tx(4)+rtmp,ty(4)+2.6
	thue=((tmp+1.)/2.9)
c	thue=scale-tmp
c  colour (hue) ranges from about red to blue
	tsat=amax1(0.1,(2.*abs(tmp)/scale))
c	tbright=4.*abs(tmp)/scale
	if(tsat.gt.1.0) tsat=1.0
	write(9,6030) thue,tsat,1.0
c  arrange for bright colors unless close to ratio=1
c  outline a square symbol, set color and grey level, fill

200	continue
	return
	end
ccccc
      subroutine euler(iopt,nomen,d1,d2,d3,ior,kerr)
c                          Last revision 20nov90 UFK
      common /eula/ a
      real a(3,3)
      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.9999) goto 10
      th=acos(a(3,3))
      sth=sin(th)
      ps=atan2(a(2,3)/sth,a(1,3)/sth)
      ph=atan2(a(3,2)/sth,a(3,1)/sth)
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)
      go to 15
c
   10 ps=0.5*atan2(a(1,2),-a(1,1))
      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 *****************************************