program sod2vf
c
c  make volume fractions from an SOD file
c
c version 1,  11 v 02    adr
c
c  read orientations from TEXIN
c  symmetry elements from QUAT.SYMM.CUBIC
c  CSLs from QUAT.CSL
c
c  LOCAL VARIABLES
	integer i,j,k,indexf1,indexcf,indexf2,ii,jj,kl
	logical result,repeat
	real qq(4),qq1(4),qend(4),quint(4),PI,twopi,degrad
	integer i9,j9
c
	character*6 comp_name(12)	! names of texture components
	character*6 comp_name_xfer
	real compphi1(12,4),compcphi(12,4),compphi2(12,4)	!Euler angles
	real qcomp(12,4,4), qmis(4,2)		! quaternions for texture components
c  changed array size from 3 to 4, 24 v 02, in order to accommodate
c  the max number of sample symmetry related variants, 
c   i.e. 4 in orthorhombic sample symmetry
c
	integer comp_index(100000),comp_index2(100000)
	real thetamin,disorientmin,comp_vol(12),sumvol,criterion,sum2
c  indexes to texture components
	integer numcomps
	integer time
c
c  GLOBAL VARIABLES
	real qr1(4),qr2(4),qtwo(4,2),qresult(4),axis(3)
	integer ij,jk
c
	character*80 textitl
	character hkl*3,seclab*5,title*8,idate*8
	integer numsymm,numcsl,nsig,ngrain,numsamp
	real phi1,capphi,phi2
c
	integer intens(73,19,19),imax,iavg
	character prosa*80,prosa1*80,hh,kk,ll,lper,nomon
	character eulan*25,nomen*5,iper*6
	character*5 label,secid
	character fname*20,one,ten,fname2*20,tname*6
	real rintens(73,19,19),projs(73,19),cprojs(73,19),delph
c  NB - was 72 in the old SODCON code, which seems like a mistake
c
	real sintens(73,19,19)   ! store the intensties from the SOD here
	integer angle_type,ic,numpart(73,19,19)  ! partition numbers
	real p1,p,p2  !  Euler angles
c
	common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
	common/a2/quatcsl(4,50),numcsl,nsig(50)
	common/a3/aquat(4,100000),ngrain,phi1(100000)
     &	,capphi(100000),phi2(100000)
	common/a4/ grwt(100000),eulan
c
	common/char/ fname,hkl,seclab,title,idate,textitl
	common /t1/ rintens,projs,cprojs,intens,delph,iavg
	common /t2/ numpart
	common /rd/ sintens,angle_type,fname2   !  use for reading in data
        character*160 vtitles
c
	degrad=57.29578
	PI=3.14159265
	twopi=2.0*PI
 8	format(a)

		call textur4

	numcomps = 0
	numcomps = numcomps + 1
	comp_name(numcomps)='cube__'
	compphi1(numcomps,1)=0.
	compcphi(numcomps,1)=0.
	compphi2(numcomps,1)=0.
c		compphi1(1,2)=0.
c		compcphi(1,2)=90.
c		compphi2(1,2)=0.
c		compphi1(1,3)=90.
c		compcphi(1,3)=0.
c		compphi2(1,3)=0.
c	numcomps = numcomps + 1
c	comp_name(2)='cub15R'
c		compphi1(2,1)=0.
c		compcphi(2,1)=15.
c		compphi2(2,1)=0.
c		compphi1(2,2)=90.
c		compcphi(2,2)=90.
c		compphi2(2,2)=15.
c		compphi1(2,3)=0.
c		compcphi(2,3)=15.
c		compphi2(2,3)=90.
c	numcomps = numcomps + 1
c	comp_name(3)='cub35R'
c		compphi1(3,1)=0.
c		compcphi(3,1)=35.
c		compphi2(3,1)=0.
c		compphi1(3,2)=90.
c		compcphi(3,2)=90.
c		compphi2(3,2)=35.
c		compphi1(3,3)=35.
c		compcphi(3,3)=0.
c		compphi2(3,3)=90.
c	numcomps = numcomps + 1
c	comp_name(4)='cub20N'
c	compphi1(4,1)=20.
c	compcphi(4,1)=0.
c	compphi2(4,1)=0.
c		compphi1(4,2)=70.
c		compcphi(4,2)=90.
c		compphi2(4,2)=90.
c		compphi1(4,3)=20.
c		compcphi(4,3)=0.
c		compphi2(4,3)=90.
c	numcomps = numcomps + 1
c	comp_name(5)='cub10T'
c		compphi1(5,1)=90.
c		compcphi(5,1)=10.
c		compphi2(5,1)=0.
c		compphi1(5,2)=0.
c		compcphi(5,2)=90.
c		compphi2(5,2)=10.
c		compphi1(5,3)=0.
c		compcphi(5,3)=90.
c		compphi2(5,3)=80.
c	numcomps = numcomps + 1
c	comp_name(6)='Cube/S'
c		compphi1(6,1)=45.
c		compcphi(6,1)=15.
c		compphi2(6,1)=45.
c		compphi1(6,2)=90.
c		compcphi(6,2)=80.
c		compphi2(6,2)=10.7
c		compphi1(6,3)=0.99
c		compcphi(6,3)=79.5
c		compphi2(6,3)=79.3
c	numcomps = numcomps + 1
c	comp_name(7)='G/Brss'
c		compphi1(7,1)=15.
c		compcphi(7,1)=45.
c		compphi2(7,1)=0.
c		compphi1(7,2)=75.
c		compcphi(7,2)=90.
c		compphi2(7,2)=45.
c		compphi1(7,3)=15.
c		compcphi(7,3)=45.
c		compphi2(7,3)=90.
	numcomps = numcomps + 1
	comp_name(numcomps)='brass_'
		compphi1(numcomps,1)=35.
		compcphi(numcomps,1)=45.
		compphi2(numcomps,1)=0.

c  S {123} <63-4>  Bunge
	numcomps = numcomps + 1
	comp_name(numcomps)='123634'
	compphi1(numcomps,1) = 58.98
	compcphi(numcomps,1) = 36.70
	compphi2(numcomps,1) = 63.43
        
c	numcomps = numcomps + 1
c	comp_name(numcomps)='Goss__'
c	compphi1(numcomps,1)=0.
c	compcphi(numcomps,1)=45.
c	compphi2(numcomps,1)=0.

c	comp_name(9)='S/Brss'
c		compphi1(9,1)=42.
c		compcphi(9,1)=37.
c		compphi2(9,1)=80.
c		compphi1(9,2)=56.0157471
c		compcphi(9,2)=84.0013657
c		compphi2(9,2)=36.5793266
c		compphi1(9,3)=29.5496655
c		compcphi(9,3)=53.6531448
c		compphi2(9,3)=7.45498514
c	comp_name(10)='S___10'
c		compphi1(10,1)=50
c		compcphi(10,1)=35.
c		compphi2(10,1)=70.
c		compphi1(10,2)=56.6
c		compcphi(10,2)=78.7
c		compphi2(10,2)=33.3
c		compphi1(10,3)=26.04
c		compcphi(10,3)=57.38
c		compphi2(10,3)=18.47
c	numcomps = numcomps + 1
c	comp_name(numcomps)='C/S_11'
c	compphi1(numcomps,1)=63.
c	compcphi(numcomps,1)=28.
c	compphi2(numcomps,1)=58.

c	comp_name(12)='C___12'
c		compphi1(12,1)=81.
c		compcphi(12,1)=27.
c		compphi2(12,1)=50.
c		compphi1(12,2)=45.8
c		compcphi(12,2)=73.03
c		compphi2(12,2)=21.3
c		compphi1(12,3)=37.7
c		compcphi(12,3)=69.64
c		compphi2(12,3)=18.1
c	comp_name(4)='Dillam'
c		compphi1(4,1)=42.
c		compcphi(4,1)=71.
c		compphi2(4,1)=20.
c		compphi1(4,2)=90.
c		compcphi(4,2)=27.
c		compphi2(4,2)=45.
c		compphi1(4,3)=90.
c		compcphi(4,3)=27.
c		compphi2(4,3)=45.

	numcomps = numcomps + 1
	comp_name(numcomps)='Copper'
	compphi1(numcomps,1)=90.
	compcphi(numcomps,1)=35.
	compphi2(numcomps,1)=45.

c        comp_name(7)='110110'
c                compphi1(7,1)=85.5
c                compcphi(7,1)=50.35
c                compphi2(7,1)=4.25
c                compphi1(7,2)=87.86
c                compcphi(7,2)=39.84
c                compphi2(7,2)=5.11
c                compphi1(7,3)=1.78
c                compcphi(7,3)=86.72
c                compphi2(7,3)=39.73
c	comp_name(7)='cube2_'
c		compphi1(7,1)=90.
c		compcphi(7,1)=90.
c		compphi2(7,1)=0.
c		compphi1(7,2)=0.
c		compcphi(7,2)=90.
c		compphi2(7,2)=90.
c		compphi1(7,3)=90.
c		compcphi(7,3)=0.
c		compphi2(7,3)=90.
c	numcomps=12       !  number of texture components to check
c
c	comp_name(7)='123412'		! also known as "R" component
c	comp_name(8)='124211'
c		compphi1(8,1)=9.63
c		compcphi(8,1)=77.4
c		compphi2(8,1)=63.431
c		compphi1(8,2)=86.59
c		compcphi(8,2)=64.12
c		compphi2(8,2)=75.96
c		compphi1(8,3)=56.79
c		compcphi(8,3)=29.21
c		compphi2(8,3)=26.57
c  following Euler angles are in Kocks angles! from hkl2eul
c	comp_name(6)='123634'
c		compphi1(6,1)=80.83
c		compcphi(6,1)=74.5
c		compphi2(6,1)=56.31
c		compphi1(6,2)=0.93
c		compcphi(6,2)=57.69
c		compphi2(6,2)=18.43
c		compphi1(6,3)=31.02
c		compcphi(6,3)=36.7
c		compphi2(6,3)=26.57
c	comp_name(7)='123412'		! also known as "R" component
c		compphi1(7,1)=68.76
c		compcphi(7,1)=74.5
c		compphi2(7,1)=56.31
c		compphi1(7,2)=11.14
c		compcphi(7,2)=57.69
c		compphi2(7,2)=71.57
c		compphi1(7,3)=43.09
c		compcphi(7,3)=36.7
c		compphi2(7,3)=26.57
c	comp_name(8)='124211'
c		compphi1(8,1)=80.37
c		compcphi(8,1)=77.4
c		compphi2(8,1)=63.431
c		compphi1(8,2)=3.41
c		compcphi(8,2)=64.12
c		compphi2(8,2)=75.96
c		compphi1(8,3)=33.21
c		compcphi(8,3)=29.21
c		compphi2(8,3)=26.57
c                vtitles=''   !  start with a blank title line
c                vtitles=' Random '//comp_name(1)//' '//comp_name(2)
c     &     //' '//comp_name(3)//' '
c     $    //comp_name(4)//' '//comp_name(5)//' '//comp_name(6)
c     &    //' '//comp_name(7)//' '//comp_name(8)//' '//comp_name(9)
c     &    //' '//comp_name(10)//' '//comp_name(11)//' '//comp_name(12)

        vtitles=' Random '
	llimit=12
	do i=1,numcomps
	   ulimit=llimit+6
	   comp_name_xfer=comp_name(i)
	   do j=1,6
	      vtitles((j+llimit-1):(j+llimit-1))=comp_name_xfer(j:j)
	   enddo
	   vtitles(ulimit:ulimit)=' '
	   llimit=ulimit+1
	enddo

	write(*,*) 'RD=1 [1] or RD=2 [2] ?'
	read(*,*) idir
	if(idir.gt.2.or.idir.lt.1) stop 'invalid answer'

	do 101, i=1,numcomps
c                        vtitles=vtitles//comp_name(i)  !  and add to it
	   comp_vol(i)=0.

	   do 101, j=1,1	!  only one trip through this loop
	      if(idir.eq.1) then !  RD= 1
		 call quatB(compphi1(i,j)/degrad,compcphi(i,j)/degrad,
     &		compphi2(i,j)/degrad,qr1)
	      elseif(idir.eq.2) then !  RD= 2
		 call quatB((compphi1(i,j)+90.)/degrad,compcphi(i,j)/degrad,
     &		compphi2(i,j)/degrad,qr1)
	      endif
	      do 101, k=1,4
 101	      qcomp(i,j,k)=qr1(k)

	      do i=1,numcomps
		 do j=1,4
		    qr1(j)=qcomp(i,1,j)	!  copy quaternion into local q
		 enddo

		 if(numsamp.gt.1) then
		    do j=2,numsamp !  ignore 1st operator (identity)
		       call presamp(qr1,j,qresult) !  apply sample symmetry
		       do k=1,4
			  qcomp(i,j,k)=qresult(k) !  copy quaternion into local q
		       enddo
		    enddo
		 endif
	      enddo

	      write(*,*)
	      write(*,*) 'debug: check Euler angles:'
	      do i=1,numcomps
		 write(*,*) 'Component number ',i,', name = ',comp_name(i)
		 do j=1,numsamp
c		    write(*,*) '   Variant number ',j
		    do k=1,4
		       qr1(k)=qcomp(i,j,k) !  copy quaternion into local q
		    enddo
		    call q2eulB(p1,p,p2,qr1)
		    if(p1.lt.0.) p1 = p1 + twopi
		    if(p2.lt.0.) p2 = p2 + twopi
		    write(*,"('Variant ',i3,'  Bunge Euler angles:',3(2x,f8.1))")
     &		j,p1*degrad,p*degrad,p2*degrad
		 enddo
	      enddo

c  done with applying sample symmetry operators to each component
c
c		do 103, i=6,8
c			comp_vol(i)=0.
c			do 103, j=1,3
c				d1=compphi1(i,j)/degrad
c				d1=d1+90.
c				if(d1.gt.360.) d1=d1-360.
c				d2=compcphi(i,j)/degrad
c				d3=compphi2(i,j)/degrad
c				d3=90.0-d3
c				if(d3.lt.0.) d3=d3+360.
c				call quatB(d1,d2,d3,qr1)
c				do 103, k=1,4
c103					qcomp(i,j,k)=qr1(k)
c
c
	write(*,*) 'sod2vf: completed list of ',numcomps,' components'
	write(*,*) 'Calculates a set of misorientations etc, from a list'
	write(*,*) 'of intensities in a SOD file.'
c
c
	write(*,*) 'Enter the filename (no quotes) for the vol. frac.'
	write(*,*) ' analysis,  e.g. al.sod'
	read(*,8) fname2
        write(*,*) 'Enter the acceptance angle (degrees)'
        read(*,*) criterion 
c       disorientmin=15.
c  the value of DISORIENTMIN sets the criterion for how close
c  an orientation must be to be identified as a special component
c
	index=20
	do 5, i=20,1,-1
c		write(*,*) 'i,fname2(i:i) ',i,fname2(i:i)
		if(fname2(i:i).eq.'.') index=i
5	continue
	index=index-1
	fname=fname2(1:index)
c	fname2((index+1):(index+1))=char(0)
	write(*,*) 'fname2,index,fname'
	write(*,*) fname2,',   ',index,',   ',fname
c
	open(unit=2,file=fname(1:index)//'.comps',status='unknown')
	write(2,*) 'Acceptance angle (degrees) = ',criterion
	write(2,*) vtitles
c
c	do 9000, i9=1,9
	   one=char(48+i9)
c
	   do 10, i=1,72
	      do 10, j=1,19
		do 10, k=1,19
		  rintens(i,j,k)=0.0
		  intens(i,j,k)=0
10	continue
c
c  loop for first time counter: change if different time stepping used!
c	   do 8990, j9=1,9
c	      time=int(float(j9)*10**i9)
c	      ten=char(48+j9)
c	      fname=fname2(1:index)//'.'//ten//one
c	      inquire(file=fname(1:(index+3))//'.wts',exist=repeat)
c	      if(.not.repeat) then
c		 write(*,*) 'Did not find ',fname(1:(index+3))//'.wts'
c		 goto 8990
c	      endif
c	      write(*,*) 'Found ',fname(1:(index+3))//'.wts'
c
	do 6,i=1,100000
		phi1(i)=0.
		capphi(i)=0.
		phi2(i)=0.
		grwt(i)=0.
		aquat(1,i)=0.
		aquat(2,i)=0.
		aquat(3,i)=0.
		aquat(4,i)=0.
		comp_index(i)=0
		comp_index2(i)=0
6	continue
	do 7, i=1,numcomps
	   comp_vol(i)=0.
 7	   continue
c
		call readsod(index,numsamp)   !  get the texture data
c
		ic=0
		do ii=1,19
			do jj=1,19
				do kl=1,19
					ic=ic+1
					phi1(ic)=5.0*float(ii-1)
					capphi(ic)=5.0*float(jj-1)
					phi2(ic)=5.0*(kl-1)  ! sets angles; phi2 varies fastest
c		write(*,*) ii,jj,kl,ic,phi1(ic),capphi(ic),phi2(ic)
					call quatB(phi1(ic)/degrad
     &			,capphi(ic)/degrad,phi2(ic)/degrad,qq)
					aquat(1,ic)=qq(1)
					aquat(2,ic)=qq(2)
					aquat(3,ic)=qq(3)
					aquat(4,ic)=qq(4)
					rlower=amax1(0.,(capphi(ic)-2.5))/degrad
					rupper=amin1(90.,(capphi(ic)+2.5))/degrad
					grwt(ic)=sintens(ii,jj,kl)*(cos(rlower)-cos(rupper))
					if(jj.eq.1.or.jj.eq.19) grwt(ic)=grwt(ic)*0.5
					if(kl.eq.1.or.k.eq.19) grwt(ic)=grwt(ic)*0.5
c  take care of edge effects
				enddo
			enddo
		enddo
	ngrain=ic
c
	      write(*,*) 'number of grains = ',ngrain
c
c  get the orientations, symms etc.
c
	      sumvol=0.		! use for total volume
	      do 1000, i=1,ngrain
                disorientmin=criterion   !  set the minimum angle
		sumvol=sumvol+grwt(i)
		 qq1(1)=aquat(1,i)
		 qq1(2)=aquat(2,i)
		 qq1(3)=aquat(3,i)
		 qq1(4)=aquat(4,i)
		 qmis(1,1)=aquat(1,i)
		 qmis(2,1)=aquat(2,i)
		 qmis(3,1)=aquat(3,i)
		 qmis(4,1)=aquat(4,i)
c
c
c  now we run through all the identified texture components
c  and identify the closest one (minimum misorientation)
c
		do 300, j=1,numcomps
			do 280, k=1,numsamp	! samp. symm. variants of each component
				qmis(1,2)=qcomp(j,k,1)
				qmis(2,2)=qcomp(j,k,2)
				qmis(3,2)=qcomp(j,k,3)
				qmis(4,2)=qcomp(j,k,4)
				call misquat(qmis,thetamin)
				if(thetamin.lt.disorientmin) then
					disorientmin=thetamin
					comp_index(i)=j
					comp_index2(i)=k
				endif
280			continue
300		continue
c
			indexf1=((i-1)/19/19)+1
			indexcf=(mod((i-1)/19,19))+1
			indexf2=(mod((i-1),19))+1
			numpart(indexf1,indexcf,indexf2)=comp_index(i)
c			write(*,*) 'partition: ',indexf1,indexcf,indexf2,i,
c     &  numpart(indexf1,indexcf,indexf2)
c
		if(comp_index(i).ne.0) then
			comp_vol(comp_index(i))=comp_vol(comp_index(i))+grwt(i)
c			write(*,*) 'Component identified for grain number ',i
c			write(*,*) 'Index= ',comp_index(i),
c	1		     '   Variant Index= ',comp_index2(i)
c			write(*,*) 'Component name = ',comp_name(comp_index(i))
		endif
c
c
 1000		       continue
c
c
	write(*,*) 'Report on volume fractions..... '
	do 1005, i=1,numcomps
		sum2=sum2+comp_vol(i)
		write(*,*) 'Component ',i,', ',comp_name(i), 
	1	'  vol. frac.= ',100.*comp_vol(i)/sumvol,' %'
1005		continue
c		write(*,*) 'Time = ',time
		write(2,'(12(1x,f7.3))') 
	1	     (comp_vol(i)/sumvol,i=1,numcomps)
	write(*,*) 'Sum of volumes = ',sumvol
	write(*,*) 'Sum of components = ',sum2
	write(*,*) 'Random = ',100.*(sumvol-sum2)/sumvol,' %'
c		pause
c
			call wod(index,criterion)
c
 8990		    continue
 9000		 continue
c
	call exit
	end
c
c  ------------------
c
	subroutine readsod(index,numsamp)
c
	implicit none
	integer nval,index
	character*80 textitl
	character prosa*80,prosa1*80,hh,kk,ll,lper,nomon
	character eulan*25,nomen*5
	character*5 label,secid
	character fname*20,one,ten,fname2*20,tname*6
	character hkl*3,seclab*5,title*8,idate*8
	real delth,rm,delom,pm
	integer iw,jw,iper(3),iavg,ngr
	integer llimit,ulimit,ik,nquads,numsamp
	logical rescale
	character afile*20, file*19, out1*22
	integer nvint(20,20,20),nscale
	integer ntype,i,j,k,ii,jj,kl,angle_type,nvmx
	real scalef
	real sintens(73,19,19)   ! store the intensties from the SOD here
c
	common/char/ fname,hkl,seclab,title,idate,textitl
	common /header/ ntype,nvmx,delth,rm,delom,pm,iw,jw,iper,ngr
	common /rd/ sintens,angle_type,afile   !  use for reading in data
c
c	print *,'Name of input file (.SOD, .SMH or .COD) ?'
c	read(*,1011) afile
 1011	format(a)
	print *,' '
c
	kl=20
c	do 154, i=20,1,-1
cc  find the right-most occurrence of a "."
c	   if(afile(i:i).eq.'.') then
c	   	kl=i
c	   	goto 155
c	   endif
c154	continue
155	continue
	   file=afile(index+2:index+4)
	   if(file.eq.'sod'.or.file.eq.'SOD'.or.file.eq.'smh') then
	      ntype=1
	      goto 888
	   endif
	   if(file.eq.'cod'.or.file.eq.'COD'.or.file.eq.'cmh') then
	      ntype=2
	      goto888
	   endif
c	   if(file.eq.'son'.or.file.eq.'SON'.or.file.eq.'snh') then
c	      ntype=3
c	      goto 888
c	   endif
c	   if(file.eq.'con'.or.file.eq.'CON'.or.file.eq.'cnh') then
c	      ntype=4
c	      goto888
c	   endif
	   stop 'This program digests only SOD and COD files.'
c	   stop 'This program digests only SOD,COD,SON and CON files.'
c

 888	open(11,file=afile,status='old')
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	read in the input data
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 887	continue
c
c     read the first header lines
c
	ulimit=19
	if(ntype.eq.1) then
		ulimit=1+(4*18/numsamp)  ! in SOD files, no of sections varies
	endif
c
	write(*,*) 'Reading in ',ulimit,'  sections based on sample symm.'
	write(*,*) 'If the reading in of a data file fails, then there'
	write(*,*) ' may be a mis-match between the sample symmetry and'
	write(*,*) ' the data file.'
	do 725,k=1,ulimit
c
	read(11,"(a)") textitl
	read(11,1332) nomen,delth,rm,delom,pm,iw,jw
     &,(iper(i),i=1,3),iavg,ngr,seclab,label
 1332 format(a5,4f5.1,5i2,2i5,2a5)
c
c  on the first section, detect the type of Euler angles in use
	if(k.eq.1) then
	nquads=int((pm+0.01)/90.)  !  +.01 to be certain
c	write(*,*) 'debug: number of lines per PHI value (nquads)= ',nquads
c
	if(ntype.eq.2) then
c  now check the range of phi2 for COD files
	if(pm.eq.90.and.numsamp.ne.4) then
		write(*,*) 'Check sample symmetry!'
		write(*,*) 'Range of azimuth = ',pm
		write(*,*) 'but no. of sample symmetry operators= ',numsamp
	endif
	if(pm.eq.180.and.numsamp.ne.2) then
		write(*,*) 'Check sample symmetry!'
		write(*,*) 'Range of azimuth = ',pm
		write(*,*) 'but no. of sample symmetry operators= ',numsamp
	endif
	if(pm.eq.360.and.numsamp.ne.1) then
		write(*,*) 'Check sample symmetry!'
		write(*,*) 'Range of azimuth = ',pm
		write(*,*) 'but no. of sample symmetry operators= ',numsamp
	endif
	endif
c
		angle_type=1
		if(nomen(5:5).eq.'B'.or.nomen(4:4).eq.'B') then
			write(*,*) 'detected Bunge angles'
		endif
c  assume Bunge angles
		if(nomen(5:5).eq.'K'.or.nomen(4:4).eq.'K') then
			write(*,*) 'detected Kocks angles'
			angle_type=2
			stop 'Can only handle Bunge angles for now!'
		endif
c  Kocks angles
	endif
c
	scalef=100./float(iavg)
	rescale=(abs(scalef-1.0).gt.0.05)
	if(rescale) write(*,*) 'Rescaling intensities for k=',k,scalef
c  re-scale all the intensities in the section by the ratio
c  of the nominal average (100) divided by IAVG in the file
c  but require more than a 5% change to do it
c
	do 724, j=1,19
		llimit=1
		do ik=1,nquads
			ulimit=llimit+17
			if(ik.eq.nquads) ulimit=ulimit+1
			read(11,31,err=932)(nvint(i,j,k),i=llimit,ulimit)
 31	format(1x,19i4)
			if(rescale) then
				do i=llimit,ulimit
					nvint(i,j,k)=int(float(nvint(i,j,k))*scalef)
				enddo
			endif
		llimit=llimit+18
		enddo
c
 724	continue
	read(11,*,end=725)
c  blank line after each section
c
 725	continue
c
 	nvmx=nvint(1,1,1)
	do 722, k=1,19
	do 722, j=1,19
	do 729, i=1,19
		if(ntype.eq.1) then   !  SOD
			sintens(k,j,i)=0.01*nvint(i,j,k)  ! switch 1st/3rd indices
		elseif(ntype.eq.2) then    !  COD
			sintens(i,j,k)=0.01*nvint(i,j,k)  !  end result of subroutine
		endif
	if(nvint(i,j,k).le.0.and.nscale.eq.2) nvint(i,j,k)=1
c  reset intensity to 1 if we are using a log scale
	if(nvint(i,j,k).gt.nvmx) nvmx=nvint(i,j,k)
 729	continue
 722	continue
c
c
	write(*,492)nvmx
 492	format('Maximum intensity =',i7)
	return
c
 932	stop 'Input file is not in the right format.'
c
	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
c
	result=.true.
	if(phi1.le.0.0) result=.false.
	if(phi1.gt.90.0) result=.false.
	if(capphi.le.0.0) result=.false.
	if(capphi.gt.90.0) result=.false.
	if(phi2.le.0.0) result=.false.
	if(phi2.gt.90.0) result=.false.
	return
	end
c
c  ------------------
c
	subroutine wod(index,criterion)
c  writes out an SOD file with index numbers
c
c  LOCAL VARIABLES
	integer ns,is,iph,iran
	real cenph
c
c  GLOBAL VARIABLES
	character hkl*3,seclab*5,title*8,idate*8
	integer intens(73,19,19),imax
	character prosa*80,prosa1*80,hh,kk,ll,lper,nomon
	character eulan*25,nomen*5		! IPER is different in WOD
	character*5 label,secid
	character fname*20
	real rintens(73,19,19),projs(73,19),cprojs(73,19)
	real delph,delth,rm,delom,pm
	character*80 textitl
	integer numsymm,numcsl,nsig,ngrain,numsamp
	integer numpart(73,19,19)
	integer iw,jw,iper(3),iavg,ngr
c
	common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
	common/a2/quatcsl(4,50),numcsl,nsig(50)
	common/a3/aquat(4,100000),ngrain,phi1(100000)
     &	,capphi(100000),phi2(100000)
	common/a4/ grwt(100000),eulan
c
	common/char/ fname,hkl,seclab,title,idate,textitl
	common /header/ ntype,nvmx,delth,rm,delom,pm,iw,jw,iper,ngr
	common /t1/ rintens,projs,cprojs,intens,delph,iavg
	common /t2/ numpart
c
1	format(2a)
    2 format(i1)
    8 format(a)
c    8 format(a\)
  111 format(a1,73x,3i2)
c
	degrad=57.29578
	iper(1)=1
	iper(2)=2
	iper(3)=3
	delph=5.0
	secid=' '//'SOD'//'B'
	seclab='Phi1='
c
c	write(*,*) 'Check IAVG: ',iavg
c
c	write(*,*) ' Enter the name for the output file (xxx.sod): '
c        read(*,*) fname
	open(7,file=fname(1:(index))//'.sod.part',status='unknown')
	write(7,*) 'Acceptance angle (degrees) = ',criterion
c
	ns=1+(4*18/numsamp)
c  assume 19 sections for ortho samp symm, 37 for mono, 73 for tricl.
c
	do 1000 is=1,ns
        iph=is
        cenph=float(iph-1)*delph
        write(label,'(f5.1)') cenph
        iran=100*nph
c
	write(7,"(a78)") textitl
c	write(7,1332) secid,5.0,90.0,5.0,90.0,1,1
      write(7,1332) secid,delth,rm,delom,pm,iw,jw
     &,(iper(i),i=1,3),iavg,ngrain,seclab,label
 1332 format(a5,4f5.1,5i2,2i5,2a5)
c
	do 1312 n=1,19
		mseg=19
c
 1320				write(7,1319) (numpart(is,n,m),m=1,mseg)
 1312 continue
 1318 format(1x,18i4)
 1319 format(1x,19i4)
      write(7,*)
c****
 1000 continue
      close(7)
c
c  now for the COD part!
	secid=' '//'COD'//'B'
	seclab='Phi2='
c
	open(7,file=fname(1:(index))//'.cod.part',status='unknown')
	write(7,*) 'Acceptance angle (degrees) = ',criterion
c
	ns=19
c  assume 19 sections only
	do is=1,ns
c		write(*,*) 'section number ',is
        iph=is
        cenph=float(iph-1)*delph
        write(label,'(f5.1)') cenph
        iran=100*nph
c
	write(7,"(a78)") textitl
c	write(7,1332) secid,5.0,90.0,5.0,90.0,1,1
      write(7,1332) secid,delth,rm,delom,pm,iw,jw
     &,(iper(i),i=1,3),iavg,ngrain,seclab,label
c
	do n=1,19
		mseg=1+(4*18/numsamp)		! change this for lower sample symmetry
		nseg=(mseg/18-1)*18
c		write(*,*) 'debug COD: mseg= ',mseg,' nseg= ',nseg
c
			if(nseg.eq.0) goto 2320
				write(7,1318) (numpart(m,n,is),m=1,nseg)
 2320				write(7,1319) (numpart(m,n,is),m=nseg+1,mseg)
	enddo
      write(7,*)
c****
 	enddo
      close(7)
c
	return
	end
c
c ________________________________________
c
c
	subroutine textur4
c
c  fourth version, for use in REXGBS, WTS2POP
c  not for use in REX2D
c  SYMMETRY OPERATORS etc. ONLY!
c
	integer numsymm,numcsl,nsig,ngrain,numsamp
	real phi1,capphi,phi2
c
c  following code taken from LApp68.f for reading in TEXIN
c
	real evm,fk1(9),nstate
	character*80 textitl
	character hkl*3,seclab*5,title*8,idate*8
	character eulan*25,nomen*5,iper*6
	character prosa*80,prosa1*80,hh,kk,ll,lper,nomon
	character*5 label,secid
	character fname*20,fname2*20
	real ph,th,ps,dph,dth,dps,sph,cph,cth,sth,sps,cps,qr(4)
	real theta,r1,r2,r3
	integer lmn
	real d1,d2,d3,grwt
c
	common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
	common/a2/quatcsl(4,50),numcsl,nsig(50)
	common/a3/aquat(4,100000),ngrain,phi1(100000)
     &	,capphi(100000),phi2(100000)
	common/a4/ grwt(100000),eulan
	common/char/ fname,hkl,seclab,title,idate,textitl
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)
c
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
c		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
  	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 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
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
c
c _____________________________
c
c
	subroutine presymm(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/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 = (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  BUT, see PRESAMP elsewhere in this set of subroutines
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
	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 (degrees)
c  taken from Sutton & Baluffi
c  ASSUMES CUBIC CRYSTAL SYMMETRY
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