! For use with single processor Monte Carlo output files

! Author: Chris Roberts, Anthony Rollett

!To compile: {F90 compiler} MC_mdf_check.f -o mdf.exe
!To run: ./mdf.exe 'keyword'
! where keyword is 6letter simulation name

! Input Files (previously generated by myCA.f and annealfinal.c)
!       cellIdealization.xml  EAorts.txt
! Output files
!       MDF and WTS files

! nbors variable will determine the search radius
! nbors=6 is 1st NN  nbors=26 is out to 3rd NN
! The differences between the distributions is minor
!
! MODS
! 9/28/06   added subroutine to calculate average mobility and energy
! 11/16/06  fixed issue with MDF calculations. problem arose from inert
!           spins and the attempts to obtain a misorientation
! 1/31/07   edited to calculate MDF for digital microstructure.
!outputs an WTS file to observe PF,IPF,ODF in POPLA
! 7/24/07   V2 is applicable to large structures >=500^3 where
!           memory demands become very critical

!____________________________________________________
! Define global variables in this module
      module global

!Enter MDF resolution and max disorientation angle possible
!Default assumes a cubic symmetry
      real,parameter::binwidth=5.0
      real,parameter::maxdisorient=62.8

      integer,allocatable,dimension(:)::freq
      integer::numbins

      integer::Q

! Define for Parallel MC (=0 or =-1) or serial MC (= -1)
      integer,parameter::inertspin=-1

! Define boundary condions
      logical,parameter::periodic=.true.

! Storage array for Texture info
      type aniso
      real, dimension(4)::quat
      integer::volume
      end type
      type(aniso),allocatable,dimension(:)::grain

      type aniso2
      real::phi1,phi,phi2
      end type
      type(aniso2),allocatable,dimension(:)::bunge

      integer::gbsegments,badtally

      end module
!_______________________________________

      program XML_TEXTURE_FIT

      use global
      implicit none

      character::keyword*40,fname*40,ones*1,tens*1
      integer::i,j
      real::num1,num2,num3
      integer::x,y,z,iargc,tally
      integer::tvol=0 ! total volume
      character::junk*20,junk2*20

      if(iargc()/=1)stop 'no filename on commandline'
      call getarg(1,fname)

      open(1,file=fname)
      read(1,*)junk,Q,junk2
      close(1)
      Q=Q-1  !!! There are Q grains but GID=0 is one of them 
             !!! so the largest spin # is actually Q-1

! Allocate memory for storage of spin values and orientations
      allocate(grain(0:Q),bunge(0:Q),stat=i)

      do x=0,Q
        grain(x)%volume=0
      enddo

! Allocate memory for MDF frequency
      numbins=int(maxdisorient/binwidth)
      numbins=numbins+1
      allocate(freq(numbins),stat=i)
      if(i.ne.0)stop 'error allocating memory'
      call texture() ! reads in texture info for each grain
  
! Reset all counters
      freq=0
      badtally=0
      gbsegments=0
! end reset counters

      call xml_mdf_calculation(fname,tvol)

      keyword='ellipsoid.yourCA'
      call mdf_output(keyword)

      !write(*,*)'Creating WTS and ANG files'
      !write(*,*)keyword
      call rod2eul(tvol,keyword)

 998  continue


      deallocate(grain,bunge,freq,stat=i)

      goto 996
 997  stop 'error opening file'
 996  continue

      end program
!___________________________________________________
      subroutine mdf_output(keyword)

      use global
      implicit none
      integer::i,j,sum,sum2
      double precision::aveE,totalE,aveM,totalM
      character::fname*40,keyword*(*)

      write(*,*) 'Calling MDF_OUTPUT'
      sum=0

      do i=1,numbins
         sum=freq(i)+sum
      enddo

      write(*,*)'gbsegments=',gbsegments
      write(*,*)'sum+badtally=',(sum+badtally)
      write(*,*)'A total of',sum,'grain boundary segments were analyzed'
      fname=trim(keyword)//'.mdf'
      open(1,file=fname)
      write(1,'(a)')'MisO    Frequency'

      do i=1,numbins
         write(1,*)(i-1)*binwidth,(1.0*freq(i))/(1.0*sum)
      enddo

      close(1)

      return
      end
!___________________________________________________
      subroutine texture()
! Reads in texin1 and texin2 orientation info,assigns it to spin number
! and converts to quaternion format

      use global
      implicit none
      integer:: i,j,id
      real:: phi1,phi,phi2,qresult(4)

      open(1,file='EAorts.txt',err=999)
!      open(2,file='texin2',err=999)      

      do i=0,Q
!         read(1,*,end=899)id,phi1,phi,phi2
         read(1,*,end=899)id,bunge(id)%phi1,bunge(id)%phi,bunge(id)%phi2
!         call b2quat(phi1,phi,phi2,qresult)
         call b2quat(bunge(id)%phi1,bunge(id)%phi,
     &        bunge(id)%phi2,qresult)
         do j=1,4
            grain(i)%quat(j)=qresult(j)
         enddo
      enddo  

 899  continue

      close(1)
      goto 897
 999  stop 'Error opening texin1 or texin2 files'
 897  return
      end
!_____________________________________________

      subroutine b2quat(phi1,phi,phi2,qresult)
! Converts Bunge Angles to Quaternions
! Expects phi1,phi,phi2 in RADIANS
      implicit none
      real::phi1,phi,phi2,qresult(4)

      qresult(1)=sin(phi/2)*cos((phi1-phi2)/2)
      qresult(2)=sin(phi/2)*sin((phi1-phi2)/2)
      qresult(3)=cos(phi/2)*sin((phi1+phi2)/2)
      qresult(4)=cos(phi/2)*cos((phi1+phi2)/2)

      return
      end
!_______________________________________________________
      subroutine misquat(qq,thetamin)

      implicit none
      real:: qq(4,2),thetamin,qresult(4),rquat(4)
      real:: disor,pi
      real:: qmax,q1max,q2max
      integer::iqindex,iq1index,iq2index,i
c     
      PI=3.14159265
c ADR
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 find the largest q value
	qmax=0.
	iqindex=0
	do  i=1,4
           qresult(i)=abs(qresult(i)) !Convert to + values
           if(qresult(i).gt.qmax) then
              qmax=qresult(i)
              iqindex=i
           endif
        enddo
c  find the next highest q component
	q1max=0.
	iq1index=0
	do  i=1,4
           if(i.eq.iqindex)then
              continue
           elseif(qresult(i).gt.q1max) then
              q1max=qresult(i)
              iq1index=i
           endif
	enddo
c  find the next highest q component
	q2max=0.
	iq2index=0
	do  i=1,4
           if(i.eq.iqindex.or.i.eq.iq1index)then
              continue
           elseif(qresult(i).gt.q2max) then
              q2max=qresult(i)
              iq2index=i
           endif
	enddo

        rquat(1)=minval(qresult)
	rquat(4)=qmax
	rquat(3)=q1max
	rquat(2)=q2max

	disor=amax1(qmax,(qmax+q1max)/sqrt(2.0),
     &	(qresult(1)+qresult(2)+qresult(3)+qresult(4))/2.0)
	if(disor.gt.1.0) disor=1.0
	if(disor.lt.-1.0) disor=-1.0
	thetamin=acos(disor)*360./pi

c  CAUTION:  note that q1<q2<q3<q4
c  whereas typical Rodrigues sorting is R1>R2>R3
	return !returns disorintation angle in DEGREES
	end
c
c _____________________________
c
	subroutine norm(q ,qnorm)

        implicit none
	real :: q(4),qnorm
        integer::i

	qnorm=0
	do i=1,4
		qnorm=qnorm+q(i)*q(i)
	end do

	qnorm=sqrt(qnorm)

        return
	end
!__________________________________________________________

	subroutine rod2eul(tvol,fname)

! Does not CONVERT Rod --> Euler 	
c       converts Rodrigues vector to Bunge Euler angles 
c       d1== phi1,  d2== Phi,  d3== phi2
c
c       input file : orts.txt (obtained from D. Saylor's annealfinal.c)
c       output file : fridyTexin (Texin1/2 format for rex3d)

        use global
        implicit none
	integer:: i,j,index,tvol,io=10
	real:: d1,d2,d3,rodr(3),r2d
        real:: sum,diff,phi1,Phi,phi2	
        real,parameter:: pi=3.1415926535
        character::fname*40,fname2*40
        
        fname2=trim(fname)//'.odf'
	open(22,file=fname2, status='unknown')

	
	write(22,*) 'EAorts.txt'  
	write(22,*) 'Evm    F11    F12    F13    F21    F22    F23    F31
     &    F32    F33    nstate'
	
	write(22,*) '0.000  1.000  0.000  0.000  0.000  1.000  0.000  0.00
     &0  0.000  1.000  1'
      
        write(22,601) 
601	format('Bunge:phi1   PHI   phi2  ,gr.wt.,  tau,    taus;taumo
     &des/tau;   XYZ= 1 2 3')

        r2d=180.0D0/pi
	do i=0,Q
		!read(21,*,iostat=io) index, rodr(1),rodr(2),rodr(3)
		!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
		!phi1=(180/pi)*d1
		!Phi=(180/pi)*d2
		!phi2=(180/pi)*d3
           write(22,1000) bunge(i)%phi1*r2d,bunge(i)%phi*r2d,
     &          bunge(i)%phi2*r2d,grain(i)%volume/(1.0*tvol),1.0
	enddo

        close(22)

 1000   format(2x,3(g10.4,1x),2(g8.3,1x))
        write(*,*) "fridyTexin should be moved to texin1 and texin2."
        
        return
        end

!__________________________________________________________________

      subroutine xml_mdf_calculation(fname,tvol)
      ! Calculates MDF based on xml measurements

      use global
      implicit none
      integer::i,j,io=10,tvol
      integer::id,neigh_id,area,neighbors,tempval
      real::qq(4,2),thetamin,rvol,rarea
      character::junk*20,junk2*2,fname*(*)
      
      write(*,*)'Entering XML_mdf_calculation'
      open(1,file=fname,iostat=io)
      if(io/=0)stop 'Error opening cellIdealization file'
      
      read(1,*)!<nregions>
      read(1,*)!

      tvol=0

      !write(*,*)'Q=',Q
      do i=0,Q
         
         read(1,*,end=765) junk! <region>

         read(1,*)junk,id

         read(1,*)junk,rvol

         grain(id)%volume=int(rvol)
         tvol=tvol+int(rvol)

         do j=1,9 
            read(1,*)           ! Dummy line reads
         enddo
         
         read(1,*)junk,neighbors,junk2

         ! Orientation of Grain
         qq(1,1)=grain(id)%quat(1)
         qq(2,1)=grain(id)%quat(2)
         qq(3,1)=grain(id)%quat(3)
         qq(4,1)=grain(id)%quat(4)
         
         do j=1,neighbors
            read(1,*)        ! <patch>
            read(1,*)junk,neigh_id
            read(1,*) ! <average_outward_normal>
            read(1,*)junk,rarea
            area=int(rarea)
            read(1,*) ! <nfacets>
            read(1,*)  ! /<patch>

            qq(1,2)=grain(neigh_id)%quat(1)
            qq(2,2)=grain(neigh_id)%quat(2)
            qq(3,2)=grain(neigh_id)%quat(3)
            qq(4,2)=grain(neigh_id)%quat(4)

            call misquat(qq,thetamin)
            
            tempval=int(thetamin/binwidth)+1
            
            if(tempval.lt.1.or.tempval.gt.numbins)then
               badtally=badtally+1
            endif

            freq(tempval)=freq(tempval)+1

         enddo

         read(1,*) junk ! /<region>

      enddo

 765  close(1)
      
      return
      end subroutine xml_mdf_calculation
      
