program wts2pop ! compile with ! g77 -O3 -o wts2pop wts2pop-[date].f ! if any problems with the code are suspected, use ! g77 -ffortran-bounds-check -g -O0 -o wts2pop wts2pop-[date].f ! if your PATH is not correctly set up, try /usr/local/bin/g77 ! bin orientations from a ".WTS" file into a discrete ! S[ample] O[rientation] D[istribution], i.e. SOD in popLA format ! which means slices in phi1 ! also COD, slices in phi2 ! also pole figures: ! edit the list of generating poles in order to add different ! pole figures ! version 1, 23 xii 01 adr ! v2, iv 03 adr ! v3, vi 03, adr - cleaned up sample symmetry so that only one ! specification per component required ! v4, xii 04, adr - added binning into triclinic sample symmetry ! v5 ix 05, adr - added binning into inverse pole figures ! v6 xi 05, adr - added SOD input ! v7 ix 06, adr - longer file names ! v8 xii 06, adr - command line arguments ! v9 i 07, adr - consolidate near-identical points on input (efficiency) ! vii 07 - partial preparation for orthorhombic xtals ! xii 07 - cleaned up some details in normalization of ODs (NORMOD) ! ? 08 - inserted HCP capability ! vi 08 - cleaned up ATAN2 problem for PF and IPF calc. ! vii 08 - added CTF read capability (need to do ANG also!) ! vii 08 - checked out HCP functionality ! vii 08 - added output of discrete PF files (for GMT plotting) ! vii 08 minor cleanup - no lines >72 (for John Butler) ! viii 08 - added ability to read quaternion sets (Joel Bernier) ! i 09 - cleaned up titles line for output of vol fracs ! iii 09 - allowed reading WTS file with arbitrary extension (e.g. TXT) ! iv 09 - starts sequence with 01, sorted arb. extension ! vii 09 - changed default output in q2eulB to give 0 instead of 45 degrees ! ix 09 - corrected a problem in NORMOD with SODs out to 360 in phi1 ! DON"T FORGET TO UPDATE THE "LAST REVISED" print!! ! read orientations from a ".wts" file with Bunge or Kocks angles ! symmetry elements from QUAT.SYMM.CUBIC ! CSLs from QUAT.CSL ! LOCAL VARIABLES implicit none real scalar integer i,j,k,l,ipoint,pc,dot_index integer indexf1,indexcf,indexf2,ijk,irange real theta,phi,d1,d2,d3,discomps logical result,repeat,only_one real qq(4),qq1(4),qend(4),quint(4),PI,twopi,degrad integer i9 , j9 , i9lower , i9upper , j9lower , j9upper character*6 comp_name(12) ! names of texture components real compphi1(12,3),compcphi(12,3),compphi2(12,3) !Euler angles real qcomp(12,3,4), qmis(4,2) ! quaternions for texture components integer comp_index(2000000),comp_index2(2000000) real thetamin,disorientmin,comp_vol(12),sumvol c indexes to texture components integer numcomps integer time 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 c integer numsymm,numcsl,nsig,ngrain,numsamp c real phi1,capphi,phi2 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character eulan*25,nomen,iper*6 character*5 label,secid character fname*50,one,ten,fname2*50,tname*6 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph real triclinic(73,19,19) ! for tricilinic sample symmetry SODs integer numpolfigs,numpoles(10) real hklgen(3,10) ! allow for 10 pole figures real hklarray(3,10,48) ! all possible 48 general poles real tmp(3),tmp1(3),tmp2(3),rtmp(3) logical already real qqsym(4) c data numpolfigs /3/ c data hklgen(1,1),hklgen(2,1),hklgen(3,1) /1.,0.,0./ c data hklgen(1,2),hklgen(2,2),hklgen(3,2) /1.,1.,1./ c data hklgen(1,3),hklgen(2,3),hklgen(3,3) /1.,0.,1./ real polearray(73,19,10) ! accumulator for pole fig integer ipolearray(73,19,10),iavgpf(10) ! integer form of pole fig real invpfgen(3,10) ! allow for 10 pole figures real ivparray(3,10,16) ! store poles for IVPs real invpolearray(73,19,10) ! accumulator for pole fig integer iinvpolearray(73,19,10),iavginvpf(10) ! integer form of pole fig character type*3 ! pass to WPF c }}}} everything the same as for pole figures, just invert the rotations! integer iq,iq2 logical read_another real rhkl(3),ruvw(3),dismin(12),wtcomp(12),tolerns,totwt,bbtmp c RHKL is for reading in HKL, RUVW for UVW for specifying components c DISMIN is used to track min. mis. angle for each component c WTCOMP is for the weighting assigned to each component c TOTWT is the total weight integer ncount ! use for component count integer llimit,ulimit ! loop counters for names logical Kocks ! if TRUE, then output in Kocks angles (for use with WEIGHTS) integer iq_sod ! see below logical read_sod ! if TRUE then input is from .SOD file(s) integer index_sod ! for file name indexing with SOD logical l_NoWeight ! for using uniform weight with each line (WTS only) integer NoWeight ! see above logical read_ctf ! for reading CTF files (from Channel/HKL) logical read_wquat ! for reading WQUAT (from Joel Bernier) logical read_any ! for reading arb. extension WTS file integer numcsl,nsig(50),numsymm,numsamp,ngrain real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50) real aquat(4,2000000),phi1(2000000),capphi(2000000) $ ,phi2(2000000) real grwt(2000000) common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg ! common /t1/ rintens,projs,cprojs,intens,delph,iavg character vtitles*130 integer vtitles_upper data vtitles_upper / 130 / ! edited 9 Jan 09 character*6 comp_name_xfer common /poles/ polearray,ipolearray,hklgen,iavgpf,numpolfigs character pfnum*1 ! use for writing discrete PFs ! character secid*5,hh*1,kk*1,ll*1 ! also for discrete PFs integer ns integer num_args character inline*80 integer date_time(8) character stamp*99 character tim*8,dat*9 integer iq_symm ! crystal symmetry integer equivorth , equivtri ! for counting how many copies in the space integer equivorthLim , equivtriLim data equivtriLim / 1 / ! set limit for no. of copies data equivorthLim / 1 / integer iqsampoff ! CODE:: call newdate(dat,tim) stamp = 'wts2pop: '//DAT//' '//TIM num_args = iargc() ! arguments: ! 1 IQ ! 2 IQ_SOD reading from WTS or SOD ! 3 filename base ! 4 weighting of points ! 5 symmetry if(num_args.le.0) then print*,'WTS2POP Usage: Example:' print*,'/code/texture/wts2pop 0 0 generic.wts 0 cubic' print* print*,'In more detail ...' print*,'./wts2pop Built-in_comps? [0=yes]' print*,' Read .WTS [ 0 ] or .SOD [ = 1 ] ' $ ,'or .CTF [ = 2 ] or .WQUAT [ = 3 ]' $ ,'or .any-extension [ = 4 ] ' print*,' filename ' print*,' read weights (for .WTS files) [0=yes, 1=no] ' print*,' crystal symmetry [ cubic | hex | orth ] ' print* endif ! are we using the components in the code? iq = -1 ! use built-in components, or read in from STDIN if(num_args.gt.0) then CALL GetArg(1,inline) c print*,'debug: afile = ',afile if(inline.ne.'') then read(inline,"(i4)") iq if(iq.lt.0 .and. iq.gt.1) stop 'bad value for IQ' endif endif ! normal is to read WTS files, but 1= read .SOD iq_sod = -1 if(num_args.gt.1) then CALL GetArg(2,inline) read(inline,"(i4)") iq_sod if( iq_sod .eq. 0 ) print*,'reading WTS' if( iq_sod .eq. 1 ) print*,'reading SOD' if( iq_sod .eq. 2 ) print*,'reading CTF' if( iq_sod .eq. 3 ) print*,'reading WQUAT (Bernier)' if( iq_sod .eq. 4 ) print*,'reading WTS, arb extension (TXT)' if( iq_sod.lt.0 .and. iq_sod.gt.4 ) then stop 'bad value for IQ_SOD' endif endif ! now get the filename for the data file(s) fname2 = '' if(num_args.gt.2) then CALL GetArg(3,inline) ipoint = len_trim(inline) if( ipoint.gt.1 ) fname2 = inline(1:ipoint) endif ! weights in col 4? NoWeight = -1 if(num_args.gt.3) then CALL GetArg(4,inline) c print*,'debug: inline = ',inline l_NoWeight = .false. if(inline.eq.'1') then l_NoWeight = .true. NoWeight = 1 elseif(inline.eq.'0') then c l_NoWeight = .true. NoWeight = 0 endif endif iq_symm = 0 if(num_args.gt.4) then CALL GetArg(5,inline) ! print*,'5th arg: ',inline ipoint = 0 ipoint = index(inline,'cubic') if(ipoint.gt.0) then iq_symm = 1 print*,'cubic xtal symmetry requested' endif if (ipoint .eq. 0 ) then ipoint = index(inline,'hex') if(ipoint.gt.0) then iq_symm = 2 print*,'hexagonal xtal symmetry requested' endif endif if (ipoint .eq. 0 ) then ipoint = index(inline,'orth') if(ipoint.gt.0) then iq_symm = 3 print*,'orthorhombic xtal symmetry requested' endif endif if(ipoint.eq.0) stop 'bad symmetry spec.?' else print*,'Which crystal symmetry ?' print*,'1=cubic 2=hex 3=orth ...' read*,iq_symm if( iq_symm .lt. 1 .or. iq_symm .gt. 3 ) then iq_symm = 1 print*,'cubic xtal symmetry assumed' endif endif ! print*,'iq,NoWeight,iq_sod,fname2,iq_symm' ! print*,iq,NoWeight,iq_sod,fname2,iq_symm Kocks = .false. tolerns = 15. ! just in case if(iq_symm.eq.1) then numpolfigs=3 hklgen(1,1)=1. hklgen(2,1)=0. hklgen(3,1)=0. hklgen(1,2)=1. hklgen(2,2)=1. hklgen(3,2)=1. hklgen(1,3)=1. hklgen(2,3)=0. hklgen(3,3)=1. elseif(iq_symm.eq.2) then numpolfigs=3 hklgen(1,1) = 1. hklgen(2,1) = 0. hklgen(3,1) = 0. c gives (100) or (2 -1 -1 0) PF hklgen(1,2) = 0. hklgen(2,2) = 0. hklgen(3,2) = 1. c gives (001) or (0001) PF hklgen(1,3) = 0. hklgen(2,3) = 1. hklgen(3,3) = 0. c gives (1 0 -1 0) PF endif degrad=57.29578 PI=3.14159265 twopi=2.0*PI 8 format(a) call textur4(iq_symm) ! read in the symmetry operators ncount=0 print*,'List of built-in components:' if(iq_symm.eq.1) then ncount=ncount+1 comp_name(ncount)='cube ' compphi1(ncount,1)=0. compcphi(ncount,1)=0. compphi2(ncount,1)=0. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='NDcube' c$$$ compphi1(ncount,1) = 45. c$$$ compcphi(ncount,1) = 0. c$$$ compphi2(ncount,1) = 0. c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1) c$$$ $ ,compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='10RDcb' c compphi1(ncount,1)=0. c compcphi(ncount,1)=10. c compphi2(ncount,1)=0. c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='20RDcb' c compphi1(ncount,1)=0. c compcphi(ncount,1)=20. c compphi2(ncount,1)=0. c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='30RDcb' c compphi1(ncount,1)=0. c compcphi(ncount,1)=30. c compphi2(ncount,1)=0. c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='rot_Gs' c compphi1(ncount,1)=90.0 c compcphi(ncount,1)=45.0 c compphi2(ncount,1)=0.0 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='20r321' c$$$ compphi1(ncount,1)=216.388153 c$$$ compcphi(ncount,1)=19.26544 c$$$ compphi2(ncount,1)=329.008026 c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='30r321' c compphi1(ncount,1)=217.786163 c compcphi(ncount,1)= 28.8845215 c compphi2(ncount,1)= 330.406036 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='20r123' c compphi1(ncount,1)=71.48 c compcphi(ncount,1)=11.91 c compphi2(ncount,1)=-55.39 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='30r123' c compphi1(ncount,1)=75.56 c compcphi(ncount,1)=17.8 c compphi2(ncount,1)=-51.31 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='10Copp' c compphi1(ncount,1)=85.86 c compcphi(ncount,1)=29.27 c compphi2(ncount,1)=130.80 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c comp_name(ncount)='20Copp' c compphi1(ncount,1)=81.81 c compcphi(ncount,1)=23.14 c compphi2(ncount,1)=126.74 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='Goss ' c$$$ compphi1(ncount,1)=0. c$$$ compcphi(ncount,1)=45. c$$$ compphi2(ncount,1)=0. c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1) c$$$ $ ,compphi2(ncount,1) c$$$ c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='brass ' c$$$ compphi1(ncount,1)=35. c$$$ compcphi(ncount,1)=45. c$$$ compphi2(ncount,1)=0. c$$$ print*,'Comp ',ncount,': ' c$$$ & ,comp_name(ncount),compphi1(ncount,1) c$$$ & ,compcphi(ncount,1),compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='Dillam' ! 4 4 11 11 11 -8 compphi1(ncount,1) = 45. compcphi(ncount,1) = 0. compphi2(ncount,1) = 30. print*,'Comp ',ncount,': ',comp_name(ncount) $ ,compphi1(ncount,1) & ,compcphi(ncount,1),compphi2(ncount,1) c$$$! Copper in Bunge angles c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='Copper' c$$$ compphi1(ncount,1) = 90. c$$$ compcphi(ncount,1) = 35. c$$$ compphi2(ncount,1) = 45. c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1) c$$$ $ ,compphi2(ncount,1) c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='echCop' c$$$ compphi1(ncount,1) = 90. c$$$ compcphi(ncount,1) = 10. c$$$ compphi2(ncount,1) = 45. c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1) c$$$ $ ,compphi2(ncount,1) c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='231124' ! 1 2 3 6 3 -4 c$$$ compphi1(ncount,1)=64.934 c$$$ compcphi(ncount,1)=74.499 c$$$ compphi2(ncount,1)=33.69 c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1) c$$$ $ ,compphi2(ncount,1) c$$$ c$$$ ncount=ncount+1 c$$$ comp_name(ncount)='111110' ! 1 1 1 1 -1 0 c$$$ compphi1(ncount,1)=0.0 c$$$ compcphi(ncount,1)=54.7 c$$$ compphi2(ncount,1)=45.0 c$$$ print*,'Comp ',ncount,': ',comp_name(ncount) c$$$ & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='110110' ! 1 1 0 1 -1 0 c compphi1(ncount,1)=85.5 c compcphi(ncount,1)=50.35 c compphi2(ncount,1)=4.25 c print*,'Comp ',ncount,': ',comp_name(ncount) c & ,compphi1(ncount,1),compcphi(ncount,1),compphi2(ncount,1) c ncount=ncount+1 c comp_name(ncount)='001110' ! 0 0 1 1 1 0 c compphi1(ncount,1)=0. c compcphi(ncount,1)=0. c compphi2(ncount,1)=45. 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 elseif(iq_symm.eq.2) then ncount=ncount+1 comp_name(1)='TiCmp1' compphi1(ncount,1) = 0. compcphi(ncount,1) = 40. compphi2(ncount,1) = 0. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='TiCube' compphi1(ncount,1) = 0. compcphi(ncount,1) = 0. compphi2(ncount,1) = 0. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='Ti__90' compphi1(ncount,1) = 0. compcphi(ncount,1) = 90. compphi2(ncount,1) = 0. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='TiP_15' compphi1(ncount,1) = 15. compcphi(ncount,1) = 40. compphi2(ncount,1) = 60. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='TiP_30' compphi1(ncount,1) = 30. compcphi(ncount,1) = 40. compphi2(ncount,1) = 55. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='TiP_45' compphi1(ncount,1) = 45. compcphi(ncount,1) = 40. compphi2(ncount,1) = 45. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) ncount=ncount+1 comp_name(ncount)='Ti_Rex' ! Rex component compphi1(ncount,1) = 0. compcphi(ncount,1) = 30. compphi2(ncount,1) = 30. print*,'Comp ',ncount,': ',comp_name(ncount) & ,compphi1(ncount,1),compcphi(ncount,1) $ ,compphi2(ncount,1) endif numcomps=ncount ! number of texture components to check vtitles=' Time(MCS) ' 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 comp_vol(i)=0. do j=1,3 call quatB(compphi1(i,j)/degrad,compcphi(i,j)/degrad, 1 compphi2(i,j)/degrad,qr1) do k=1,4 qcomp(i,j,k)=qr1(k) enddo enddo enddo do i = ulimit , vtitles_upper vtitles(ulimit:ulimit) = ' ' end do print*,'Title line = ' print* , vtitles(1:vtitles_upper) do i=1,numpolfigs c write(*,*) c write(*,*) 'pole figure set number ',i ipoint=1 tmp(1)=hklgen(1,i) tmp(2)=hklgen(2,i) tmp(3)=hklgen(3,i) call donorm(tmp) hklarray(1,i,1)=tmp(1) hklarray(2,i,1)=tmp(2) hklarray(3,i,1)=tmp(3) c write(*,"(i5,3(2x,f5.2))") ipoint,(tmp(j),j=1,3) do j=1,numsymm do k=1,4 qqsym(k)=quatsymm(k,j) enddo call qrvec(qqsym,tmp,tmp1) c write(*,*) 'symm oper no, tmp1: ',j,tmp1 c apply a symmetry operator already=.false. do k=1,ipoint tmp2(1)=hklarray(1,i,k) tmp2(2)=hklarray(2,i,k) tmp2(3)=hklarray(3,i,k) if(abs(1.-scalar(tmp1,tmp2)).lt.1e-2) already=.true. c write(*,*) 'k,dot product,tmp1,tmp2 ',k,scalar(tmp1,tmp2),tmp1,tmp2 c we've found a match for this pole enddo c write(*,*) 'tmp1 ',tmp1 if(.not.already) then ipoint=ipoint+1 do k=1,3 c write(*,*) 'k,tmp1(k) ',k,tmp1(k) hklarray(k,i,ipoint)=tmp1(k) enddo c hklarray(2,i,ipoint)=tmp1(2) c hklarray(3,i,ipoint)=tmp1(3) c write(*,"('ipoint= ',i3,3(2x,f5.2))") c & ipoint,(hklarray(l,i,ipoint),l=1,3) c write(*,*) 'CHECK/symm oper no, tmp1: ',j,tmp1 endif enddo numpoles(i)=ipoint enddo c }}}}} do i = 1,10 invpfgen(1,i) = 0 invpfgen(2,i) = 0 invpfgen(3,i) = 0 do j = 1,16 ! max number of eq poles for orthorhombic symmetry ivparray(1,i,j) = 0. ivparray(2,i,j) = 0. ivparray(3,i,j) = 0. enddo enddo invpfgen(1,1) = 1 invpfgen(2,2) = 1 invpfgen(3,3) = 1 ivparray(1,1,1) = 1. ivparray(1,1,2) = -1. c gives sample 100 ivparray(2,2,1) = 1. ivparray(2,2,2) = -1. c gives sample 010 ivparray(3,3,1) = 1. ivparray(3,3,2) = -1. c gives sample 001 c note for inv pole figures, for now: since we are going to calculate c only the RD, TD, ND (sample 1, 2 & 3) directions, there's no need c for this apparatus, other than the set of individual points. write(*,*) write(*,*) 'WTS2POP is primarily for binning lists of discrete ' write(*,*) 'orientations into discretized orientation ' write(*,*) 'distributions in popLA format for subsequent' write(*,*) 'analysis or plotting.' print* write(*,*) 'It also performs a volume fraction analysis' print*,' Look for the volume fraction output in a *.comps file' print*,' also outputs standard popLA .COD, .SOD' print*,' also pole figures, and inverse PFs as .ppf, .ipf' print*,' also the first 10000 data points as .dpf and dip' print* print*,'Sample symmetry is treated as follows:' print*,'Each orientation is binned assuming triclinic (none)' print*,' symmetry, and also assuming orthorhombic (mmm)' print*,'The output for triclinic has _360 added to the name' print*,'whereas the orthorhombic just has the rootname' write(*,*) print*,' Written by A.D. Rollett, Carnegie Mellon University' print*,' Please acknowledge use of the program in whatever' print*,' fashion seems appropriate and inform me of errors' print*,' at rollett@andrew.cmu.edu ' print* print*,' last revised : September 2009' print* print*,'Note that the capture/tolerance angle = ' $ , tolerns , ' degrees' print*,'To compare against random, use random200k.wts as input' print* if(iq_symm.eq.2) then print* print*,' IMPORTANT: HKL in the pole figures is w.r.t.' print*,' the Cartesian reference axes: ' print*,' There are two possible settings:' print*,'A: in Metz, e.g., with Bunge Euler angles, ' print*,' the X-axis is parallel to 1 0 -1 0 , NOT 2 -1 -1 0' print*,' Thus, (100) represents (1 0 -1 0)' print*,' Thus, (010) represents (2 -1 -1 0)' print* print*,'B: in CMU, e.g., with Bunge Euler angles, ' print*,' the X-axis is parallel to 2 -1 -1 0, NOT 1 0 -1 0' print*,' Thus, (100) represents (2 -1 -1 0)' print*,' Thus, (010) represents (1 0 -1 0)' print* endif if(iq.eq.-1) then print*,'Use the built-in set of components' print*,' (see the source code for this) [=0]' print*,'or enter the individual components by hand [=1]?' read(*,*) iq endif if(iq.eq.0) then continue ! nothing new to do elseif(iq.eq.1) then write(*,*) 'Do you want to enter Miller Indices [=0]' write(*,*) 'or, Euler angles (Bunge!) [=1]' read(*,*) iq2 pc=1 read_another=.true. if(iq2.eq.0) then write(*,*) 'Keep entering h k l u v w (max 12)' write(*,*) 'When you are done, enter 0 0 0 0 0 0' do jk=1,12 numcomps=pc ! just to be sure if(read_another) then write(*,*) 'Enter indices for ' $ ,jk,'th comp. (avoid h=k=0)' read(*,*) rhkl,ruvw if(rhkl(1).eq.999.) then read_another=.false. numcomps=pc-1 else write(*,*) 'and 6-char name for this component' read(*,"(a)") comp_name(pc) call donorm(rhkl) call donorm(ruvw) bbtmp = rhkl(3) / (sqrt(rhkl(1)**2 $ + rhkl(2)**2+rhkl(3)**2)) if(bbtmp.gt.1.0) bbtmp=1.0 if(bbtmp.lt.-1.0) bbtmp=-1.0 compcphi(pc,1)=degrad*acos(bbtmp) bbtmp=rhkl(2)/(sqrt(rhkl(1)**2+rhkl(2)**2)) if(bbtmp.gt.1.0) bbtmp=1.0 if(bbtmp.lt.-1.0) bbtmp=-1.0 compphi2(pc,1)=degrad*acos(bbtmp) bbtmp = ruvw(3) * sqrt((rhkl(1)**2 $ +rhkl(2)**2+rhkl(3)**2) & / (rhkl(1)**2+rhkl(2)**2) $ / (ruvw(1)**2+ruvw(2)**2+ruvw(3)**2)) if(bbtmp.gt.1.0) bbtmp=1.0 if(bbtmp.lt.-1.0) bbtmp=-1.0 compphi1(pc,1)=degrad*acos(bbtmp) write(*,*) 'check the Euler Angles: ', & compphi1(pc,1),compcphi(pc,1),compphi2(pc,1) pc=pc+1 endif c number of texture components to check endif enddo elseif(iq2.eq.1) then write(*,*) 'Keep entering phi1 PHI phi2 (max 12)' write(*,*) 'When you are done, enter 999. 999. 999.' do jk=1,12 numcomps=pc ! just to be sure if(read_another) then write(*,*) 'Enter angles for ',jk,'th comp.' read(*,*) compphi1(pc,1),compcphi(pc,1) $ ,compphi2(pc,1) write(*,*) 'and 6-char name for this component' read(*,"(a)") comp_name(pc) if(compphi1(pc,1).eq.999.) then read_another=.false. numcomps=pc-1 endif c number of texture components to check pc=pc+1 endif enddo else stop 'not an option!' endif else stop 'not an option!' endif iq2 = 0 c$$$ write(*,*) c$$$ print*,' Output in Bunge [default] or Kocks angles [=1]?' c$$$ read(*,*) iq2 c$$$ if(iq2.eq.1) then c$$$ Kocks = .true. c$$$ print*,'Kocks angles will be used for output!' c$$$ endif ! so rarely need the Kocks angle output that I'm discarding it, xi 05 if(iq_sod.lt.0) then print* print*,'Read .wts [ = 0 ], or .sod [ = 1 ]' $ ,', or .CTF [ = 2 ] , or .WQUAT [ = 3] ' $ ,'or .any-extension [ = 4 ] ?' read(*,*) iq_sod endif if( iq_sod .lt. 0 .or. iq_sod .gt. 4 ) stop 'illegal entry' read_sod = .false. if( iq_sod .eq. 1 ) then read_sod = .true. print*,'The output files will have _Z at the end of the root' print*,' to differentiate them from the input set' endif read_ctf = .false. if( iq_sod .eq. 2 ) then read_ctf = .true. endif read_wquat = .false. if( iq_sod .eq. 3 ) then read_wquat = .true. endif read_any = .false. if( iq_sod .eq. 4 ) then read_any = .true. endif print*,'__________' if(fname2.eq.'') then print*,'Enter the KEYWORD for this run, e.g. cube55' print*, 'Operates on a series of .wts files from a gr gr/rex' print*, 'simulation with names like cube55.11.wts (or .sod),' print*, ' cube55.21.wts, etc.' print*, 'If there is only one file with the name cube55.wts,' print*, ' wts2pop will try to analyze that one first.' only_one=.false. read(*,"(a)") fname2 endif c$$$ pointx=20 c$$$ do 5, i=20,1,-1 c$$$c print*, 'i,fname2(i:i) ',i,fname2(i:i) c$$$ if(fname2(i:i).eq.' ') pointx=i c$$$ 5 enddo dot_index = 0 dot_index = index(fname2,'.') if(dot_index.eq.0) dot_index = index(fname2,' ') dot_index = dot_index - 1 c fname = fname2(1:dot_index) fname = fname2(1:dot_index)//'_Z' c the extra 2 characters will be used for output when using SOD as input index_sod = dot_index if(read_sod) index_sod = index_sod+2 c see above for filename modification for writing out results from SOD input ! fname2((dot_index+1):(dot_index+1)) = char(0) ! print*, 'fname2,dot_index,fname' ! print*, fname2,', ',dot_index,', ',fname open(unit=2,file=fname(1:index_sod)//'.comps',status='unknown') write(2,*) vtitles(1:vtitles_upper) ! do we need to enquire about weighting? if(NoWeight.lt.0) then l_NoWeight = .false. if(.not.read_sod) then print*,'weights from .WTS input [=0], or unweighted [=1]?' read(*,*) NoWeight if(NoWeight.eq.1) then l_NoWeight = .true. else NoWeight = 0 endif endif endif ! now deal with crystal symmetry if not specified in arguments if ( iq_symm .eq. 0 ) then print*,'Enter xtal symmetry [cubic/hex/orth] :' read"(a)",inline ipoint = 0 ipoint = index(inline,'cubic') if(ipoint.gt.0) then iq_symm = 1 print*,'cubic xtal symmetry requested' endif if (ipoint .eq. 0 ) then ipoint = index(inline,'hex') if(ipoint.gt.0) then iq_symm = 2 print*,'hexagonal xtal symmetry requested' endif endif if (ipoint .eq. 0 ) then ipoint = index(inline,'orth') if(ipoint.gt.0) then iq_symm = 3 print*,'orthorhombic xtal symmetry requested' endif endif if(ipoint.eq.0) stop 'bad symmetry spec.?' endif ! check to see what's out there ... if( read_sod ) then inquire(file=fname2(1:(dot_index))//'.sod',exist=only_one) endif if ( iq_sod .eq. 0 ) then inquire(file=fname2(1:(dot_index))//'.wts',exist=only_one) endif if ( iq_sod .eq. 4 ) then inquire( file=fname2 , exist=only_one ) endif if( read_ctf ) then inquire(file=fname2(1:(dot_index))//'.ctf',exist=only_one) endif if( read_wquat ) then inquire(file=fname2(1:(dot_index))//'.wquat',exist=only_one) endif if(only_one) then time = 0 if( iq_sod .eq. 4 ) then print*, 'found ' , fname2 endif if( iq_sod .eq. 3 ) then print*, 'found ',fname2(1:(dot_index))//'.wquat' endif if( iq_sod .eq. 2 ) then print*, 'found ',fname2(1:(dot_index))//'.ctf' endif if( iq_sod .eq. 1 ) then print*, 'found ',fname2(1:(dot_index))//'.sod' endif if( iq_sod .eq. 0 ) then print*, 'found ',fname2(1:(dot_index))//'.wts' endif else print*, 'did not find: ' $ , fname2(1:(dot_index))//'.wts/.sod/.ctf/arb' print*, 'Will look for a series of WTS/SOD files' endif i9upper = 9 i9lower = 0 if(only_one) then i9upper = 1 i9lower = 1 end if do 9000, i9 = i9lower , i9upper one = char(48+i9) ! loop for first time counter: change if different time stepping used! j9upper = 9 j9lower = 0 if(only_one) then j9upper=1 j9lower=1 end if do 8990, j9 = j9lower , j9upper do i=1,73 do j=1,19 do k=1,10 polearray(i,j,k)=0. ipolearray(i,j,k)=0 invpolearray(i,j,k)=0. iinvpolearray(i,j,k)=0 end do end do end do ! zero out the pole & inv. pole figure accumulators do 10, i=1,72 do j=1,19 do k=1,19 rintens(i,j,k)=0.0 intens(i,j,k)=0 triclinic(i,j,k) = 0. end do end do 10 end do ! zero out the OD accumulators if(.not.only_one) then time=int(float(j9)*10**i9) ten = char(48+j9) fname = fname2(1:dot_index)//'.'//ten//one print*,'debug: fname= ',fname if( read_wquat ) then inquire(file=fname(1:(dot_index+3))//'.wquat' $ ,exist=repeat) endif if(read_ctf) then inquire(file=fname(1:(dot_index+3))//'.ctf' $ ,exist=repeat) endif if(read_sod) then inquire(file=fname(1:(dot_index+3))//'.sod' $ ,exist=repeat) endif if ( iq_sod .eq. 0 ) then inquire(file=fname(1:(dot_index+3))//'.wts' $ ,exist=repeat) endif if(.not.repeat) then print*, 'Did not find ' $ ,fname(1:(dot_index+3))//'.wts/sod' goto 8990 endif print*, 'Found ',fname(1:(dot_index+3))//'.wts/sod/ctf' endif c figure out file names 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 if ( read_any ) then call texture_any( fname2 , only_one , l_NoWeight ) endif if ( read_wquat ) then call texture_wquat( dot_index $ , only_one , l_NoWeight ) endif if ( read_ctf ) then call texture_ctf(dot_index,only_one,l_NoWeight) endif if ( read_sod ) then call readsod(dot_index,only_one) ! get the texture data endif if ( iq_sod .eq. 0 ) then call textur5(dot_index,only_one,l_NoWeight) endif print*, 'number of grains read in = ',ngrain c get the orientations, symms etc. ! open files for discrete pole figures do j = 1 , numpolfigs write(pfnum,"(i1.1)") j open(unit = j+10 $ , file = fname(1:index_sod)//"_"//pfnum//'.dpf' $ , status = 'unknown' ) ns=hklgen(1,j) hh=char(48+ns) ns=hklgen(2,j) kk=char(48+ns) ns=hklgen(3,j) ll=char(48+ns) secid='('//hh//kk//ll//')' write( j+10 , "(5a)" ) secid , ' ' $ , stamp(1:20) , ' / ' , textitl(1:66) enddo ! num. PFs do j = 1 , 3 write(pfnum,"(i1.1)") j open(unit = j+20 $ , file = fname(1:index_sod)//"_"//pfnum//'.dip' $ , status = 'unknown' ) ns = invpfgen(1,j) hh = char(48+ns) ns = invpfgen(2,j) kk = char(48+ns) ns = invpfgen(3,j) ll = char(48+ns) secid = '('//hh//kk//ll//')' write( j+20 , "(5a)" ) secid , ' ' $ , stamp(1:20) , ' / ' , textitl(1:66) enddo ! num. IPFs ! start calcs. sumvol=0. ! use for total volume ! print*,'turn off sample symmetry? yes = 1' ! read*,iqsampoff ! if( iqsampoff .eq. 1 ) numsamp = 1 ! leaving these lines in for now, but, part of the attempt to ! be able to read in weak textures correctly do 1000, i=1,ngrain equivorth = 0 equivtri = 0 ! counting no of copies if(mod(i,10000).eq.0) then print("('Working on grain no. ',i9)"),i endif 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 now to calculate each pole figure do j=1,numpolfigs c print*, 'pole figure no. ',j do k=1,numpoles(j) do l=1,3 tmp(l)=hklarray(l,j,k) enddo ! print*, 'pole no. ',k,'; dirn= ',tmp call qrvec ( qq1 , tmp , tmp1 ) ! call invqrvec( qq1 , tmp , tmp1 ) ! print*,'k: ',k,' rotated vec: ',tmp1 c rotate each pole by the grain's orientation call donorm(tmp1) ! just to be sure if( tmp1(3) .ge. 0. ) then ! don't index if on lower hemisph. if( tmp1(3) .gt. 1.0 ) tmp1(3) = 1.0 ! just in case theta = acos(tmp1(3))*degrad ! declination if((abs(tmp1(2)).lt.1e-35) $ .and. (abs(tmp1(1)).lt.1e-35)) then phi = pi/4. else phi = atan2(tmp1(2),tmp1(1))*degrad ! azimuth ! changed the order of 1 and 2, because in fortran it's atan2(y-component, x-component) ! see L3_Components_pt2_16Jan08_v6.ppt for notes; ADR, jun 08 endif if(phi.lt.0.) phi = phi+360. ! write out to discrete PFs - no more than 10000 points if ( i .lt. 10000 ) then write ( j+10 , "(3(1x,f8.3))") $ phi , 90.-theta , grwt(i) endif ! now we calc indexes indexf1 = int(theta/5.0+0.5)+1 if(indexf1.lt.1.or.indexf1.gt.19) then print*, 'problem: indexf1, theta = ' $ ,indexf1,theta stop endif indexcf = int(phi/5.0+0.5)+1 if(indexcf.lt.1.or.indexcf.gt.73) then print*, 'problem: indexcf, phi = ' $ ,indexcf,phi stop endif c$$$ print*, 'th,ph,i1,i2 ',theta,phi,indexf1,indexcf polearray(indexcf,indexf1,j) $ = polearray(indexcf,indexf1,j)+grwt(i) endif enddo ! loop over k, the set of poles enddo ! loop over j, each PF c END OF POLE FIGURE CALC c }}}} BEGIN INVERSE POLE FIGURE - RD, TD and ND do j=1,3 c print*, 'inverse pole figure no. ',j do k=1,2 c always 2 poles for a sample direction, + and - do l=1,3 tmp(l) = ivparray(l,j,k) enddo c print*, 'pole no. ',k,'; dirn= ',tmp do ijk = 1,numsymm ! crystal symmetry call postsymm(qq1,ijk,quint) call invqrvec(quint,tmp,tmp1) c rotate each inverse pole by the grain's orientation back to crystal axes call donorm(tmp1) ! just to be sure if(tmp1(3).ge.0.) then ! don't index if on lower hemisph. if(tmp1(3).gt.1.0) tmp1(3)=1.0 ! just in case theta=acos(tmp1(3))*degrad ! declination if((abs(tmp1(2)).lt.1e-35).and. 1 (abs(tmp1(1)).lt.1e-35)) then phi=pi/4. else phi=atan2(tmp1(2),tmp1(1))*degrad ! azimuth ! changed the order of args 1 and 2 - see above endif if(phi.lt.0.) phi=phi+360. ! write out to discrete PFs - no more than 500 points if ( i .lt. 500 ) then write ( j+20 , "(3(1x,f8.3))" ) $ phi , 90.-theta , grwt(i) endif ! now we calc indexes indexf1=int(theta/5.0+0.5)+1 if(indexf1.lt.1.or.indexf1.gt.19) then print*, 'problem: indexf1, theta = ' $ ,indexf1,theta stop endif indexcf=int(phi/5.0+0.5)+1 if(indexcf.lt.1.or.indexcf.gt.73) then print*, 'problem: indexcf, phi = ' $ ,indexcf,phi stop endif c print*, 'th,ph,i1,i2 ',theta,phi,indexf1,indexcf invpolearray(indexcf,indexf1,j) = 1 invpolearray(indexcf,indexf1,j) $ + grwt(i) endif enddo enddo ! loop over the set of IPFs enddo ! loop over each IPF ! END OF INVERSE POLE FIGURE CALCULATION c now we run through all the identified texture components c and identify the closest one (minimum misorientation) disorientmin=15. do j=1,numcomps dismin(j)=999999. wtcomp(j)=0. enddo c ### the value of DISORIENTMIN sets the criterion for how close c an orientation must be to be identified as a special component c do 300, j=1,numcomps c do 280, k=1,3 ! variants of each component c qmis(1,2)=qcomp(j,k,1) c qmis(2,2)=qcomp(j,k,2) c qmis(3,2)=qcomp(j,k,3) c qmis(4,2)=qcomp(j,k,4) c call misquat(qmis,thetamin) c if(thetamin.lt.disorientmin) then c disorientmin=thetamin c comp_index(i)=j c comp_index2(i)=k c endif c 280 continue c 300 continue c if(comp_index(i).ne.0) then c comp_vol(comp_index(i))=comp_vol(comp_index(i))+grwt(i) c print*, 'Component identified for grain number ',i c print*, 'Index= ',comp_index(i), c 1 ' Variant Index= ',comp_index2(i) c print*, 'Component name = ',comp_name(comp_index(i)) c endif do 900, j=1,numsymm ! crystal symmetry call postsymm(qq1,j,quint) c**** this next section bins the result with ONLY crystal symm applied call q2eulB(d1,d2,d3,quint) 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 print*, d1, d2, d3 if(Kocks) then d1 = d1 - 90. if(d1.lt.0.) d1 = d1 + 360. d3 = 90. - d3 if(d3.lt.0.) d3 = d3 + 360. endif ! convert to Kocks if(iq_symm.eq.1) then call testangs_tri(d1,d2,d3,result) elseif(iq_symm.eq.2) then call testangs_tri_hex(d1,d2,d3,result) elseif(iq_symm.eq.3) then call testangs_tri_ort(d1,d2,d3,result) endif if( result .and. equivtri.lt.equivtriLim ) then ! this should yield 1 point per orientation ! thereby enforcing the cubic crystal symmetry within the 360x90x90 space equivtri = equivtri + 1 indexf1 = int(d1/5.0+0.5)+1 indexcf = int(d2/5.0+0.5)+1 indexf2 = int(d3/5.0+0.5)+1 ! DEBUG LINES -> c$$$ print* c$$$ print*, 'At triclinic section' c$$$ print*, 'Angles for grain number ',i c$$$ print*, phi1(i),capphi(i),phi2(i) c$$$ print*, 'Angles returned from sort ' c$$$ print*, d1, d2, d3 c$$$ print*, 'Loop Indices i , j = ',i,j c$$$ print*, 'Binning for this grain:' c$$$ print*, indexf1,indexcf,indexf2 triclinic(indexf1,indexcf,indexf2)= 1 triclinic(indexf1,indexcf,indexf2)+grwt(i) endif c**** END triclinic section do 880, k=1,numsamp ! sample symmetry call presamp(quint,k,qend) c now we try each component and see if we find a smaller misorientation do ijk=1,numcomps do jk=1,4 qmis(jk,1)=qend(jk) ! gets the current orient qmis(jk,2)=qcomp(ijk,1,jk) ! gets the component enddo ! if(iq_symm.eq.1) then call misquat(qmis,thetamin) elseif(iq_symm.eq.2) then call comquat(qmis,qr1) call disquat_ang_hex(qr1,thetamin) ! this should speed up the calculation (uses one only copy of xtal symmetry) elseif(iq_symm.eq.3) then call comquat(qmis,qr1) call disquat_orth(qr1,thetamin) ! this should speed up the calculation (uses one only copy of xtal symmetry) endif ! new for v10 with hex if(thetamin.lt.dismin(ijk)) then dismin(ijk)=thetamin endif enddo ! ijk loop on components 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 print*, d1, d2, d3 if(Kocks) then d1 = d1 - 90. if(d1.lt.0.) d1 = d1 + 360. d3 = 90. - d3 if(d3.lt.0.) d3 = d3 + 360. endif ! convert to Kocks if(iq_symm.eq.1) then call testangs(d1,d2,d3,result) elseif(iq_symm.eq.2) then call testangs_hex(d1,d2,d3,result) elseif(iq_symm.eq.3) then call testangs_ort(d1,d2,d3,result) endif ! new in v10 to combine hex with cubic if( result .and. equivorth.lt.equivorthLim ) then ! this should yield 3 points per orientation ! thereby enforcing the cubic crystal symmetry within the 90x90x90 space equivorth = equivorth + 1 indexf1 = int(d1/5.0+0.5)+1 indexcf = int(d2/5.0+0.5)+1 indexf2 = int(d3/5.0+0.5)+1 c DEBUG LINES -> c$$$ print*, c$$$ print*, 'Angles for grain number ',i c$$$ print*, phi1(i),capphi(i),phi2(i) c$$$ print*, 'Angles returned from sort ' c$$$ print*, d1, d2, d3 c$$$ print*, 'Loop Indices = ',i,j,k c$$$ print*, 'Binning for this grain:' c$$$ print*, indexf1,indexcf,indexf2 rintens(indexf1,indexcf,indexf2)= 1 rintens(indexf1,indexcf,indexf2)+grwt(i) endif 880 continue 900 continue c so we've gone thru all the components; now evaluate what to do c with this grain c c print*, c print*, 'For grain number ',i comp_index(i)=0 ! make sure discomps=997. totwt=0. do ijk=1,numcomps c write(*,"('comp no, min ang',i4,f8.3)") ijk,dismin(ijk) if(dismin(ijk).lt.discomps) then discomps=dismin(ijk) if(dismin(ijk).lt.tolerns) comp_index(i)=ijk endif if(dismin(ijk).lt.tolerns) then if(dismin(ijk).gt.1.e-3) then wtcomp(ijk)=1./(dismin(ijk)**4) else wtcomp(ijk)=1.e+12 endif totwt=totwt+wtcomp(ijk) c write(*,"('comp no, wt.',i4,f8.3)") ijk,wtcomp(ijk) endif enddo if(comp_index(i).ne.0) then c print*, 'Component ',comp_index(i), c & ' identified for grain number ',i endif if(totwt.gt.0.) then do ijk=1,numcomps if(dismin(ijk).lt.tolerns) then wtcomp(ijk)=wtcomp(ijk)/totwt ! normalize c print*, 'comp no , wt ',ijk,wtcomp(ijk) comp_vol(ijk)=comp_vol(ijk) $ + (grwt(i)*wtcomp(ijk)) c print*, 'Component name = ',comp_name(comp_index(i)) endif enddo endif ! done with this grain (at last!) ! print*,'gr, equivs TRI, ORTHO: ',i,equivtri,equivorth 1000 continue c c print*, 'Report on volume fractions..... ' do 1005, i=1,numcomps write(*, $ "('For component ',a $ ,' vol. frac.= ',f8.3,' %')") 1 comp_name(i),100.*comp_vol(i)/sumvol 1005 continue print*, 'Time = ',time write(2,'(2x,i12,10(1x,f7.3))') time 1 ,(comp_vol(i)/sumvol,i=1,numcomps) c pause call normpf(polearray,ipolearray,hklgen $ , iavgpf,numpolfigs) call normpf(invpolearray,iinvpolearray $ ,invpfgen,iavginvpf,3) c normalize the pole and inverse pole figures ! print*, 'done with normpf' type = 'ppf' call wpf(index_sod,only_one,type,polearray,ipolearray & ,hklgen,iavgpf,numpolfigs) c write out the pole figures ! call wpf_GMT(index,only_one) c write out GMT formatted PFs: longitude, latitude, intensity ! write(*,*) 'done with wpf for use with GMT plotting' ! implement this later on! type = 'ipf' print*,'IPF output disabled for now' call wpf(index_sod,only_one,type $ , invpolearray,iinvpolearray & ,invpfgen,iavginvpf,3) c write out the inverse pole figures ! print*, 'done with wpf, ipf' ! call wpf_GMT(index,only_one) c write out GMT formatted PFs: longitude, latitude, intensity ! write(*,*) 'done with wpf for use with GMT plotting' ! implement this later on! c triclinic sample symm irange = 4 call normod(irange) !normalizes call wod(index_sod,only_one,irange,Kocks) ! writes SOD file call wodcod(index_sod,only_one,irange,Kocks) ! writes COD file ! orthorhombic sample symm irange = 1 call normod(irange) !normalizes call wod(index_sod,only_one,irange,Kocks) ! writes SOD file call wodcod(index_sod,only_one,irange,Kocks) ! writes COD file 8990 continue 9000 continue do j = 1 , numpolfigs close( unit = j+10 ) enddo ! num. PFs do j = 1 , 3 open( unit = j+20 ) enddo ! num. IPFs print*,'_________' print*,'3 degr smoothing recommended for pole figures:' print*,'/code/popLA/popLA_ADR/smooth ' & ,fname(1:(index_sod))//'.ppf 3. n' print*,' ... similar for .ipf' print* print*,'Use this command to convert PF files to GMT format:' print*, $ '/code/popLA/Batchfiles.SVEN/pf2GMT ' & , fname(1:(index_sod))//'.mpf' print*,' and then this (mpf example) to generate plots:' print*, 1 '/code/popLA/Batchfiles.SVEN/Draw_stereograms 3 ' & ,fname(1:(index_sod))//'_mpf',' PF' print* print*,'_________' print*,'For the .COD or .SOD files, try smoothing with:' print*,'/code/paul.lee.codes/smoothsod' print* print*,'then, plot it with ' print*,'[/code/paul.lee.codes/]sodcon ' & , fname(1:(index_sod))//'.cmh 1 1 1' call exit end ! ______________________________ C*****NEWDATE / borrowed from psplot.txt subroutine newdate(dat,tim) integer date_time(8) character*10 dummyarg character*(*)dat,tim character*36 monnam data monnam/'JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC'/ call date_and_time(dummyarg,dummyarg,dummyarg,date_time) c Convert new date/time to the old way write(dat(1:2),'(i2.2)')date_time(3) mind=(date_time(2)-1)*3+1 dat(4:6)=monnam(mind:mind+2) iyronly=date_time(1)-(date_time(1)/100)*100 write(dat(8:9),'(i2.2)')iyronly dat(3:3)='-' dat(7:7)='-' write(tim(1:2),'(i2.2)')date_time(5) write(tim(4:5),'(i2.2)')date_time(6) write(tim(7:8),'(i2.2)')date_time(7) tim(3:3)=':' tim(6:6)=':' return end c c ------------------ c subroutine testangs(phi1, capphi, phi2,result) c tests to see if the angles fall within the specified range ! which for now is just 0<=angle<=90 ! changed it to accept 0. exactly, vii 09 ! LOCAL VARIABLES integer i real phi1,capphi,phi2 logical result ! CODE:: result=.true. if(phi1.lt.0.0) result=.false. if(phi1.gt.90.0) result=.false. if(capphi.lt.0.0) result=.false. if(capphi.gt.90.0) result=.false. if(phi2.lt.0.0) result=.false. if(phi2.gt.90.0) result=.false. return end c c ------------------ c subroutine testangs_tri(phi1, capphi, phi2,result) c tests to see if the angles fall within the specified range c appropriate for the specified xtal symmetry and no sample symmetry c LOCAL VARIABLES implicit none integer i real phi1,capphi,phi2 logical result result=.true. if(phi1.le.0.0) result = .false. if(phi1.gt.360.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 testangs_hex(phi1, capphi, phi2,result) c tests to see if the angles fall within the specified range c which for now is just 0 c$$$ do 4000, i = 1 , nuphi1 c$$$ print*, 'NORMOD: Block for phi1 = ',(i-1)*5. c$$$ do 3900, j=1,19 c$$$ write(*,3800) (intens(i,j,k),k=1,19) c$$$ 3800 format(1x,19i5) c$$$ 3900 end do c$$$ 4000 end do ! pause c print*, 'Real Numbers:' c do 5000, i=1,19 c print*, 'Block for phi1 = ',(i-1)*5. c do 4900, j=1,19 c write(*,4800) (xintens(i,j,k),k=1,19) c4800 format(1x,19f6.1) c4900 continue c5000 continue c pause c now to deal with the projections c first the SOP do j=1,19 do k=1,19 projs(j,k)=0.0 do is = 1,nuphi1 projs(j,k)=projs(j,k)+xintens(is,j,k) enddo projs(j,k)=projs(j,k)/float(nuphi1) enddo enddo c and now the COP do j=1,nuphi1 do k=1,19 cprojs(j,k)=0.0 do is = 1,19 cprojs(j,k)=cprojs(j,k)+xintens(j,k,is) enddo cprojs(j,k)=cprojs(j,k)/19. enddo enddo c DEBUG LINES -> c print*, 'Check SOP' c do 2000, i=1,19 c write(*,1800) (projs(i,j),j=1,19) c1800 format(1x,19f4.1) c2000 continue c print*, 'Check COP' c do 3000, i=1,19 c write(*,1800) (cprojs(i,j),j=1,19) c3000 continue return end ! ------------------ subroutine normpf(polearray,ipolearray,hklgen,iavgpf,numpolfigs) c normalizes an Pole Figure array c cells centered on the value are assumed (jw=kw=1 in popLA) c which means that the cells at 0 and 90 have to be treated c differently from the rest c LOCAL VARIABLES integer ns,is,iph,iran,nuphi1,nuphi2 real cenph,phi1,capphi,phi2,oldcapphi,newcapphi,darea(19) real sum,pfmax(10),scale,dvol(73,19,19),voltot real sumpf(10) c DVOL contains the volume of each cell c c GLOBAL VARIABLES integer intens(73,19,19),imax,iavg real rintens(73,19,19),projs(73,19),cprojs(73,19),delph real hklgen(3,10),check(10) ! allow for 10 pole figures integer numpolfigs real polearray(73,19,10) ! accumulator for pole fig integer ipolearray(73,19,10),iavgpf(10) ! integer form of pole fig real triclinic(73,19,19) common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg c common /poles/ polearray,ipolearray,hklgen,iavgpf,numpolfigs c changed from common to arguments so as to re-use for IVPs c CODE:: degrad=57.29578 oldcapphi=0. do i=1,19 capphi = float(i-1)*5.0 newcapphi = capphi+2.5 newcapphi = amin1(newcapphi,90.0) darea(i) = cos(oldcapphi/degrad)-cos(newcapphi/degrad) oldcapphi=newcapphi 100 enddo c print*, 'completed 100 loop; DAREA()= ',darea sum = 0.0 nuphi1 = 73 !changed for full pole figures nuphi2 = 19 voltot = 360. do 200, i = 1,nuphi1 do 190, j = 1,nuphi2 k=1 dvol(i,j,k) = darea(j)*5./voltot if(i.eq.1.or.i.eq.nuphi1) dvol(i,j,k)=dvol(i,j,k)*0.5 c if(j.eq.1.or.j.eq.nuphi2) dvol(i,j,k)=dvol(i,j,k)*0.5 c this was dealt with above sum=sum+dvol(i,j,k) c deal with the edge and corner cells 190 enddo 200 enddo c$$$ print*, 'completed 200 loop; check VOL=1? vol.sum=',sum c print*, 'no. of pole figures =',numpolfigs do i=1,numpolfigs !make this 72 for no sample symm sumpf(i)=0.0 do j=1,72 do k=1,19 sumpf(i) = sumpf(i)+polearray(j,k,i) enddo enddo enddo c$$$ print*, 'completed 1000 loop; SUM= ',(sumpf(i),i=1,numpolfigs) do i=1,numpolfigs !make this 72 for no sample symm if(sumpf(i).eq.0.0) stop 'failed - no intensity?!' scale = 1.0/sumpf(i) pfmax(i)=0.0 do j=1,72 do k=1,19 polearray(j,k,i)=scale*polearray(j,k,i)/dvol(j,k,1) if(polearray(j,k,i).gt.pfmax(i)) then pfmax(i)=polearray(j,k,i) endif enddo enddo c$$$ print*, 'For PF #',i,' MAX value = ',pfmax(i) 2000 enddo do i = 1,numpolfigs check(i) = 0. do j = 1,72 do k = 1,19 check(i) = check(i)+polearray(j,k,i)*dvol(j,k,1) enddo enddo c$$$ print*,'check sum= for PF #',i,' = ',check(i),' =1?' enddo do i=1,numpolfigs scale=1.0 iavgpf(i)=100 if(pfmax(i).gt.99.98) then scale=scale*99.98/pfmax(i) iavgpf(i)=max0(1,int(9998./pfmax(i))) endif c this brings everything back down into i4 format! pfmax(i)=0.0 do j=1,72 do k=1,19 polearray(j,k,i)=scale*polearray(j,k,i) ipolearray(j,k,i)=int(100.0*polearray(j,k,i)) if(polearray(j,k,i).gt.pfmax(i)) then pfmax(i)=polearray(j,k,i) endif c print*, 'j,k,polearray,ipolearray ' c & ,j,k,polearray(j,k,i),ipolearray(j,k,i) 2800 enddo 2900 enddo 3000 enddo do i=1,numpolfigs if(pfmax(i).gt.9999.) stop 'error in scaling!' enddo c pause 'hit return to continue' c DEBUG LINES -> c do 4000, i=1,19 c print*, 'Block for phi1 = ',(i-1)*5. c do 3900, j=1,19 c write(*,3800) (intens(i,j,k),k=1,19) c3800 format(1x,19i4) c3900 continue c4000 continue c pause 'hit return to continue' c c print*, 'Real Numbers:' c do 5000, i=1,19 c print*, 'Block for phi1 = ',(i-1)*5. c do 4900, j=1,19 c write(*,4800) (rintens(i,j,k),k=1,19) c4800 format(1x,19f4.1) c4900 continue c5000 continue return end c c ------------------ c subroutine sop_calc c calculates the SOP projection c from RINTENS into PROJS c c LOCAL VARIABLES integer ns,is,iph,iran real cenph integer iper(3) c c GLOBAL VARIABLES character hkl*3,seclab*5,title*8,idate*8 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh*1,kk*1,ll*1,lper,nomon character eulan*25,nomen ! IPER is different in WOD character*5 label,secid character fname*50 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph character*80 textitl integer numsymm,numcsl,nsig,ngrain,numsamp real triclinic(73,19,19) common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,50),numcsl,nsig(50) common/a3/aquat(4,2000000),ngrain,phi1(2000000) & ,capphi(2000000),phi2(2000000) common/a4/ grwt(2000000),eulan common/strings/ fname,hkl,seclab,title,idate,textitl common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg 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=' '//'SOP'//'B' seclab='PROJ=' ns=19 c assume 19 sections plus a projection for now c so, this version no longer usable for general sample symmetry do 1312 j=1,19 do 1200, k=1,19 projs(j,k)=0.0 cprojs(j,k)=0.0 do 1000 is=1,ns projs(j,k)=projs(j,k)+rintens(is,j,k) cprojs(j,k)=cprojs(j,k)+rintens(k,j,is) 1000 continue projs(j,k)=projs(j,k)/float(ns) cprojs(j,k)=cprojs(j,k)/float(ns) 1200 continue 1312 continue c c DEBUG LINES -> c print*, 'Check SOP' c do 2000, i=1,19 c write(*,1800) (projs(i,j),j=1,19) 1800 format(1x,19f4.1) c2000 continue c c print*, 'Check COP' c do 3000, i=1,19 c write(*,1800) (cprojs(i,j),j=1,19) c3000 continue c return end c c ------------------ c subroutine wpf(pointx,only_one,type,polearray $ ,ipolearray,hklgen,iavgpf,numpolfigs) c writes out an pole figure file, or inverse pole figure file c LOCAL VARIABLES integer n,ns,is,iph,iran,pointx real cenph,degrad integer iper(3) logical only_one c GLOBAL VARIABLES character hkl*3,seclab*5,title*8,idate*8 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh*1,kk*1,ll*1,lper,nomon character eulan*25,nomen ! IPER is different in WOD character*5 label,secid character fname*50 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph character*80 textitl integer numsymm,numcsl,nsig,ngrain,numsamp character type*3 real hklgen(3,10) ! allow for 10 pole figures real polearray(73,19,10) ! accumulator for pole fig integer ipolearray(73,19,10),iavgpf(10) ! integer form of pole fig integer numpolfigs real triclinic(73,19,19) common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,50),numcsl,nsig(50) common/a3/aquat(4,2000000),ngrain,phi1(2000000) & ,capphi(2000000),phi2(2000000) common/a4/ grwt(2000000),eulan common/strings/ fname,hkl,seclab,title,idate,textitl common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg c common /poles/ polearray,ipolearray,hklgen,iavgpf,numpolfigs c CODE:: 1 format(2a) 2 format(i1) 8 format(a) c 8 format(a\) 111 format(a1,73x,3i2) ! print*,'Entering WPF with type = ',type degrad=57.29578 iper(1)=1 iper(2)=2 iper(3)=3 delph=5.0 seclab='PROJ ' label=' ' c print*, 'Check IAVG: ',iavg if(only_one) then if(type.eq.'ppf') then open(7,file=fname(1:(pointx))//'.ppf',status='unknown') elseif(type.eq.'ipf') then open(7,file=fname(1:(pointx))//'.ipf',status='unknown') endif else if(type.eq.'ppf') then open(7,file=fname(1:(pointx+3))//'.ppf',status='unknown') elseif(type.eq.'ipf') then open(7,file=fname(1:(pointx+3))//'.ipf',status='unknown') endif endif nseg=72 ! 4 rows, 18 in each row do n=1,numpolfigs ! 1000 ! print*, 'writing pole figure no. ',n c sort out the labels ns=hklgen(1,n) hh=char(48+ns) ns=hklgen(2,n) kk=char(48+ns) ns=hklgen(3,n) ll=char(48+ns) secid='('//hh//kk//ll//')' write(7,8) textitl if(ngrain.lt.99999) then write(7,1332) secid,5.0,90.0,5.0,360.0,1,1 c write(7,1332) secid,delth,rm,delom,pm,iw,jw & ,(iper(i),i=1,3),iavgpf(n),ngrain,seclab,label else write(7,1332) secid,5.0,90.0,5.0,360.0,1,1 & ,(iper(i),i=1,3),iavgpf(n),1,seclab,label c have to fudge if dealing with a very large set of grains, e.g. from OIM 1332 format(a5,4f5.1,5i2,2i5,2a5) endif mseg=19 do m=1,mseg write(7,1318) (ipolearray(is,m,n),is=1,nseg) enddo 1318 format(1x,18i4) 1319 format(1x,19i4) write(7,*) c* ********** (commented out:) ring intensities: c* if(is.lt.ns) goto 1350 c* write(7,1) ' ' c* write(7,1332) secid,delth,rm,delom,pm,iw,jw c* &,(iper(i),i=1,3),iavg,ngr,seclab,label c* write(7,1335) c* 1335 format('ring centered at; avg.grain fraction; std.dev.:') c* do 1359 n=1,nrad c* centh=delth*(n-.5-iw*.5) c* 1359 write(7,1333) centh,frad(n),fdev(n) c* write(7,*) c* 1350 continue c* 1333 format(3f8.5) c**** 1000 enddo close(7) return end c c ------------------ c subroutine wpf_GMT(pointx,only_one,type,polearray $ ,ipolearray,hklgen,iavgpf,numpolfigs) c writes out an pole figure file c in GMT format, longitude, latitude, intensity c file extension = GPF c c LOCAL VARIABLES integer n,ns,is,iph,iran real cenph,degrad integer iper(3) logical only_one c c GLOBAL VARIABLES character hkl*3,seclab*5,title*8,idate*8 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh*1,kk*1,ll*1,lper,nomon character eulan*25,nomen ! IPER is different in WOD character*5 label,secid character fname*50 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph character*80 textitl integer numsymm,numcsl,nsig,ngrain,numsamp c real hklgen(3,10) ! allow for 10 pole figures real polearray(73,19,10) ! accumulator for pole fig integer ipolearray(73,19,10),iavgpf(10) ! integer form of pole fig integer numpolfigs real triclinic(73,19,19) character one*1 ! for file name index common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,50),numcsl,nsig(50) common/a3/aquat(4,2000000),ngrain,phi1(2000000) & ,capphi(2000000),phi2(2000000) common/a4/ grwt(2000000),eulan common/strings/ fname,hkl,seclab,title,idate,textitl ! common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg ! common /poles/ polearray,ipolearray,hklgen,iavgpf,numpolfigs c CODE:: 1 format(2a) 2 format(i1) 8 format(a) c 8 format(a\) 111 format(a1,73x,3i2) degrad=57.29578 iper(1)=1 iper(2)=2 iper(3)=3 delph=5.0 seclab='PROJ ' label=' ' c write(*,*) 'Check IAVG: ',iavg c write(*,*) ' Enter the name for the output file (xxx.sod): ' c read(*,*) fname nseg=72 ! 4 rows, 18 in each row do n=1,numpolfigs ! 1000 one=char(48+n) write(*,*) 'writing pole figure no. ',n if(only_one) then open(7,file=fname(1:(pointx))//'_'//one//'.gpf' & ,status='unknown') else open(7,file=fname(1:(pointx+3))//'_'//one//'.gpf' & ,status='unknown') endif c sort out the labels ns=hklgen(1,n) hh=char(48+ns) ns=hklgen(2,n) kk=char(48+ns) ns=hklgen(3,n) ll=char(48+ns) secid='('//hh//kk//ll//')' c write(7,8) textitl if(ngrain.lt.99999) then write(7,1332) secid,5.0,90.0,5.0,360.0,1,1 c write(7,1332) secid,delth,rm,delom,pm,iw,jw & ,(iper(i),i=1,3),iavgpf(n),ngrain,seclab,label else write(7,1332) secid,5.0,90.0,5.0,360.0,1,1 & ,(iper(i),i=1,3),iavgpf(n),1,seclab,label c have to fudge if dealing with a very large set of grains, e.g. from OIM 1332 format(a5,4f5.1,5i2,2i5,2a5) endif c mseg = 19 do m = 1,mseg do is = 1,nseg write(7,*) 5.*float(is-1),' ' $ ,5.*float(mseg-m),' ' ! 19-m gives geog. latitude $ ,float(ipolearray(is,m,n))/100. enddo enddo 1318 format(1x,18i4) 1319 format(1x,19i4) close(7) 1000 enddo return end c c ------------------ c subroutine wod(pointx,only_one,irange,Kocks) c writes out an SOD file c LOCAL VARIABLES integer ns,is,iph,iran , pointx real cenph integer iper(3) logical only_one integer irange ! signals range of phi1 c GLOBAL VARIABLES character hkl*3,seclab*5,title*8,idate*8 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character eulan*25,nomen ! IPER is different in WOD character*5 label,secid character fname*50 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph character*80 textitl integer numsymm,numcsl,nsig,ngrain,numsamp real triclinic(73,19,19) logical Kocks c common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,50),numcsl,nsig(50) common/a3/aquat(4,2000000),ngrain,phi1(2000000) & ,capphi(2000000),phi2(2000000) common/a4/ grwt(2000000),eulan c common/strings/ fname,hkl,seclab,title,idate,textitl common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg c 1 format(2a) 2 format(i1) 8 format(a) c 8 format(a\) 111 format(a1,73x,3i2) print*,'Entering WOD to print out SOD' degrad=57.29578 iper(1)=1 iper(2)=2 iper(3)=3 delph=5.0 secid=' '//'SOD'//'B' seclab='Phi1=' if(Kocks) then secid=' '//'SOD'//'K' seclab='PSI =' endif c print*, 'Check IAVG: ',iavg c print*, ' Enter the name for the output file (xxx.sod): ' c read(*,*) fname c DEBUG LINES -> c do i=1,19 c do i=1,1 c print*, 'Block for phi1 = ',(i-1)*5. c do j=1,19 c write(*,"(1x,19i5)") (intens(i,j,k),k=1,19) c enddo c enddo c pause 'debug print out of INTENS' if(irange.eq.1) then if(only_one) then open(7,file=fname(1:(pointx))//'.sod',status='unknown') else open(7,file=fname(1:(pointx+3))//'.sod',status='unknown') endif c next, shear textures elseif(irange.eq.2) then if(only_one) then open(7,file=fname(1:(pointx))//'_180.sod' $ ,status='unknown') else open(7,file=fname(1:(pointx+3))//'_180.sod' $ ,status='unknown') endif c next triclinic sample symmetry elseif(irange.eq.4) then if(only_one) then open(7,file=fname(1:(pointx))//'_360.sod' $ ,status='unknown') else open(7,file=fname(1:(pointx+3))//'_360.sod' $ ,status='unknown') endif else stop 'invalid value of irange' endif c ns=20 ns = 18*irange + 2 c number of sections determined by range of phi1 c plus a projection do 1000 is=1,ns if(is.eq.ns) then iph=1 label='PROJ ' iran=100 else iph=is cenph=float(iph-1)*delph write(label,'(f5.1)') cenph iran=100*nph endif c write(7,8) textitl if(ngrain.le.99999) then write(7,1332) secid,5.0,90.0,5.0,90.0,1,1 c write(7,1332) secid,delth,rm,delom,pm,iw,jw &,(iper(i),i=1,3),iavg,ngrain,seclab,label else write(7,1332) secid,5.0,90.0,5.0,90.0,1,1 &,(iper(i),i=1,3),iavg,1,seclab,label 1332 format(a5,4f5.1,5i2,2i5,2a5) c have to fudge for v large data sets, e.g. OIM endif do 1312 n=1,19 mseg=19 ! change this for lower sample symmetry c do 1312 n=1,nrad c do 1313 m=1,mseg c 1313 if(intens(is,m,n).lt.1) intens(is,m,n)=1 ccccc to avoid 0 (taking logs in plotting routine) nseg=(mseg/18-1)*18 if(is.eq.ns) then ! test for PROJ if(nseg.eq.0) goto 1420 write(7,1318) (int(100.*projs(n,m)),m=1,nseg) 1420 write(7,1319) (int(100.*projs(n,m)),m=nseg+1,mseg) else if(nseg.eq.0) goto 1320 write(7,1318) (intens(is,n,m),m=1,nseg) 1320 write(7,1319) (intens(is,n,m),m=nseg+1,mseg) endif 1312 continue 1318 format(1x,18i4) 1319 format(1x,19i4) write(7,*) c* ********** (commented out:) ring intensities: c* if(is.lt.ns) goto 1350 c* write(7,1) ' ' c* write(7,1332) secid,delth,rm,delom,pm,iw,jw c* &,(iper(i),i=1,3),iavg,ngr,seclab,label c* write(7,1335) c* 1335 format('ring centered at; avg.grain fraction; std.dev.:') c* do 1359 n=1,nrad c* centh=delth*(n-.5-iw*.5) c* 1359 write(7,1333) centh,frad(n),fdev(n) c* write(7,*) c* 1350 continue c* 1333 format(3f8.5) c**** 1000 continue close(7) return end c c ------------------ c subroutine wodcod(pointx,only_one,irange,Kocks) c writes out an COD file implicit none c LOCAL VARIABLES integer i,j,k,pointx integer ns,is,iph,iran real cenph integer iper(3) logical only_one integer irange real degrad integer nph,mseg,nseg,n,m real rrange c GLOBAL VARIABLES character hkl*3,seclab*5,title*8,idate*8 integer intens(73,19,19),imax,iavg character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character eulan*25,nomen ! IPER is different in WOD character*5 label,secid character fname*50 real rintens(73,19,19),projs(73,19),cprojs(73,19),delph character*80 textitl integer numsymm,ngrain,numsamp real triclinic(73,19,19) real quatsymm(4,48),quatsamp(4,4),quatcsl(4,50) integer numcsl,nsig(50) real aquat(4,2000000),phi1(2000000) & ,capphi(2000000),phi2(2000000) real grwt(2000000) logical Kocks common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl common /t1/ rintens,triclinic,projs,cprojs,intens,delph,iavg c common /t1/ rintens,projs,cprojs,intens,delph,iavg c 1 format(2a) 2 format(i1) 8 format(a) c 8 format(a\) 111 format(a1,73x,3i2) print*,'Entering WODCOD to print out COD' degrad=57.29578 iper(1)=1 iper(2)=2 iper(3)=3 delph=5.0 secid=' '//'COD'//'B' seclab='Phi2=' if(Kocks) then secid=' '//'COD'//'K' seclab='PHI =' endif c DEBUG LINES -> goto 3172 c do i=1,19 do i=1,1 print*, 'Block for phi1 = ',(i-1)*5. do j=1,19 write(*,"(1x,19i5)") (intens(i,j,k),k=1,19) enddo enddo 3172 continue c pause 'debug print out of INTENS' c print*, 'Check IAVG: ',iavg c print*, ' Enter the name for the output file (xxx.cod): ' c read(*,*) fname if(irange.eq.1) then if(only_one) then open(7,file=fname(1:(pointx))//'.cod',status='unknown') else open(7,file=fname(1:(pointx+3))//'.cod',status='unknown') endif c next, shear textures elseif(irange.eq.2) then if(only_one) then open(7,file=fname(1:(pointx))//'_180.cod',status='unknown') else open(7,file=fname(1:(pointx+3))//'_180.cod' $ ,status='unknown') endif c next triclinic sample symmetry elseif(irange.eq.4) then if(only_one) then open(7,file=fname(1:(pointx))//'_360.cod',status='unknown') else open(7,file=fname(1:(pointx+3))//'_360.cod' $ ,status='unknown') endif else stop 'invalid value of irange' endif ns=20 c assume 19 sections plus a projection for now do 1000 is=1,ns if(is.eq.ns) then iph=1 label='PROJ ' iran=100 else iph=is cenph=float(iph-1)*delph write(label,'(f5.1)') cenph iran=100*nph endif rrange = float(irange)*90. write(7,8) textitl if(ngrain.le.99999) then write(7,1332) secid,5.0,90.0,5.0,rrange,1,1 c write(7,1332) secid,delth,rm,delom,pm,iw,jw & ,(iper(i),i=1,3),iavg,ngrain,seclab,label else write(7,1332) secid,5.0,90.0,5.0,rrange,1,1 & ,(iper(i),i=1,3),iavg,1,seclab,label 1332 format(a5,4f5.1,5i2,2i5,2a5) c have to fudge for v large data sets, e.g. OIM endif c mseg=19 ! change this for lower sample symmetry mseg = 18*irange + 1 nseg = (mseg/18-1)*18 do 1312 n=1,19 c do 1312 n=1,nrad c do 1313 m=1,mseg c 1313 if(intens(is,m,n).lt.1) intens(is,m,n)=1 ccccc to avoid 0 (taking logs in plotting routine) if(is.eq.ns) then ! test for PROJ if(nseg.eq.0) goto 1420 c$$$ write(7,1318) (int(100.*cprojs(n,m)),m=1,nseg) c$$$ 1420 write(7,1319) (int(100.*cprojs(n,m)),m=nseg+1,mseg) ! jul 08 - discovered that the order (n,m) should be (m,n) write(7,1318) (int(100.*cprojs(m,n)),m=1,nseg) 1420 write(7,1319) (int(100.*cprojs(m,n)),m=nseg+1,mseg) else if(nseg.eq.0) goto 1320 write(7,1318) (intens(m,n,is),m=1,nseg) 1320 write(7,1319) (intens(m,n,is),m=nseg+1,mseg) endif 1312 continue 1318 format(1x,18i4) 1319 format(1x,19i4) write(7,*) c* ********** (commented out:) ring intensities: c* if(is.lt.ns) goto 1350 c* write(7,1) ' ' c* write(7,1332) secid,delth,rm,delom,pm,iw,jw c* &,(iper(i),i=1,3),iavg,ngr,seclab,label c* write(7,1335) c* 1335 format('ring centered at; avg.grain fraction; std.dev.:') c* do 1359 n=1,nrad c* centh=delth*(n-.5-iw*.5) c* 1359 write(7,1333) centh,frad(n),fdev(n) c* write(7,*) c* 1350 continue c* 1333 format(3f8.5) c**** 1000 continue close(7) return end c c ________________________________________ c c subroutine textur4(iq_symm) c second version, for use in REXGBS, WTS2POP c not for use in REX2D integer numsymm,numsamp integer iq_symm,n_max common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp rad = 57.29578 if(iq_symm.lt.1 .or. iq_symm.gt.3) stop 'bad symmetry class?' if(iq_symm.eq.1) then n_max = 24 elseif(iq_symm.eq.2) then n_max = 12 ! Hex elseif(iq_symm.eq.3) then n_max = 4 ! Orth endif c print*, 'reading symmetry operators from "quat.symm.cubic"' if(iq_symm.eq.1) then open(unit=3,name='quat.symm.cubic',status='old') elseif(iq_symm.eq.2) then open(unit=3,name='quat.symm.hex',status='old') elseif(iq_symm.eq.3) then open(unit=3,name='quat.symm.orth',status='old') endif read(3,*) ! skip over title line read(3,*) numsymm ! read number of symmetry operators if(numsymm.gt.n_max) then numsymm = n_max print*, 'too many symmetry operators, cutting back to ',n_max else c print*, 'reading ',numsymm,' symmetry operators' endif do 1800, i=1,numsymm read(3,*) (quatsymm(ijk,i),ijk=1,4) c print*, 'symm. no. ',i,(quatsymm(ijk,i),ijk=1,4) 1800 continue close(unit=3) c print*, 'reading sample symmetry operators from "quat.symm.ort"' c which contains: c quat.symm.ort file with 4 proper rotations in quaternion form c 4 c 0 0 0 1 c 1 0 0 0 c 0 1 0 0 c 0 0 1 0 c open(unit=3,name='quat.symm.ort',status='old') c read(3,*) ! skip over title line c read(3,*) numsamp ! read number of sample symmetry operators c if(numsamp.gt.4) then c numsamp=4 c print*, 'too many symmetry operators, cutting back to 4!' c else c print*, 'reading ',numsamp,' symmetry operators' c endif c do 1810, i=1,numsamp c read(3,*) (quatsamp(ijk,i),ijk=1,4) c print*, 'samp. symm. no. ',i,(quatsamp(ijk,i),ijk=1,4) c 1810 continue c close(unit=3) c now we hard-wire the four operators into the code c so that choices by user are easier to deal with numsamp = 4 do i = 1,numsamp do ijk = 1,4 quatsamp(ijk,i) = 0. enddo enddo quatsamp(4,1) = 1. quatsamp(1,2) = 1. quatsamp(2,3) = 1. quatsamp(3,4) = 1. return end c c ________________________________________ c subroutine textur5(pointx,only_one,l_NoWeight) implicit none c fifth version, for use in REXGBS, WTS2POP c only reads in orientations; symm opers in textur4 above c not for use in REX2D c following code taken from LApp68.f for reading in TEXIN integer pointx,ng,i,j real evm,fk1(9),nstate character*80 textitl character hkl*3,seclab*5,title*8,idate*8 character eulan*25,nomen,iper*6 character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character*5 label,secid character fname*50,fname2*50 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 logical only_one logical l_NoWeight ! for unweighted input integer numcsl,nsig(50),numsymm,numsamp,ngrain real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50) real aquat(4,2000000),phi1(2000000),capphi(2000000) $ ,phi2(2000000) real grwt(2000000) common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl real tmp_wt,last_phi1,last_capphi,last_phi2 integer count ! used to consolidate orientations real rad c CODE:: rad=57.29578 c print*, 'Enter the KEYWORD for this run, e.g. "cube55"' c read(*,'(a)') fname2 c pointx=20 c do 5, i=20,1,-1 c print*, 'i,fname2(i:i) ',i,fname2(i:i) c if(fname2(i:i).eq.' ') pointx=i c5 continue c pointx=pointx-1 c fname=fname2(1:pointx) c fname2((pointx+1):(pointx+1))=char(0) c print*, 'fname2,pointx,fname' c print*, fname2,', ',pointx,', ',fname c c print*, 'reading orientations from fname....',fname(1:(pointx+3)) if(only_one) then open(unit=3,name=fname(1:(pointx))//'.wts',status='old') else open(unit=3,name=fname(1:(pointx+3))//'.wts',status='old') endif read(3,7) textitl 7 format(a) read(3,*) read(3,*) ! don't need the strain etc. c read(3,70) evm,fk1,nstate 70 format(10f7.3,i3) read(3,71) eulan,iper 71 format(a25,49x,a6) nomen=eulan if(nomen.eq.'B') print*, 'Detected Bunge angles' if(nomen.eq.'K') print*, 'Detected Kocks angles' c print 76,textitl 76 format(/ #' TEXIN file= '/1x,a/) ng = 0 count = 0 do 529 i = 1,2000000 c do 520 i=1,18 c 520 state(i,ng)=0. c read(3,50,err=530,end=530) d1,d2,d3,grwt(ng),(state(i,ng),i=1,nstate) if(l_NoWeight) then read(3,*,end=530,err=530) d1,d2,d3 tmp_wt = 1. ! grwt(ng) = 1. else read(3,*,end=530,err=530) d1,d2,d3,tmp_wt endif ng = ng + 1 c$$$ print*,'i,d1,d2,d3,tmp_wt,ng:',i,d1,d2,d3,tmp_wt,ng if ( d1.ge.0. .and. d1.lt.1.e-3 ) d1 = 1.0e-3 if ( d1.gt.89.999 .and. d1.le.90. ) d1 = 89.999 if ( d1.gt.179.999 .and. d1.le.180. ) d1 = 179.999 if ( d1.gt.359.999 .and. d1.le.360. ) d1 = 359.999 if ( d2.ge.0. .and. d2.lt.1.e-3 ) d2 = 1.0e-3 if ( d2.gt.89.999 .and. d2.le.90. ) d2 = 89.999 if ( d2.gt.179.999 .and. d2.le.180. ) d2 = 179.999 if ( d2.gt.359.999 .and. d2.le.360. ) d2 = 359.999 if ( d3.ge.0. .and. d3.lt.1.e-3 ) d3 = 1.0e-3 if ( d3.gt.89.999 .and. d3.le.90. ) d3 = 89.999 if ( d3.gt.179.999 .and. d3.le.180. ) d3 = 179.999 if ( d3.gt.359.999 .and. d3.le.360. ) d3 = 359.999 ! avoid values on the edges of the space if(i.gt.1) then if(abs(d1 - last_phi1).lt.0.1 .and. & abs(d2 - last_capphi).lt.0.1 .and. & abs(d3 - last_phi2).lt.0.1 ) then ng = ng - 1 count = count + 1 grwt(ng) = grwt(ng) + tmp_wt c$$$ print*,'SAME: i,ng,count,d1,d2,d3' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) ! angles are sufficiently similar to consoliate the point ! with the previous one else grwt(ng) = tmp_wt count = 1 c$$$ print*,'DIFFNT: i,ng,count,d1,d2,d3,wt' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(abs(... else grwt(1) = tmp_wt c$$$ print*,'NG=1: i,ng,count,d1,d2,d3,wt' c$$$ & ,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(ng.gt.1 last_phi1 = d1 last_capphi = d2 last_phi2 = d3 ! save the values ! print*, 'ng, angles.. ',ng,d1,d2,d3 c if((d1.eq.999.).and.(d2.eq.999.).and.(d3.eq.999.)) goto 530 50 format(6f8.2,24f6.2) c read(3,*,end=530) d1,d2,d3,grwt(ng) c .: for old unformatted TEXIN files if(nomen.eq.'B') then phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 elseif(nomen.eq.'K') then phi1(ng)=d1+90. if(phi1(ng).gt.360.) phi1(ng)=phi1(ng)-360. ! corrected 24 Jul 09 capphi(ng)=d2 phi2(ng)=90.-d3 if(phi2(ng).lt.0.) phi2(ng)=phi2(ng)+360. else stop 'cannot deal with that type! ' endif ! 20 dth=d2 ! dps=d1 ! if(nomen.eq.'B') then call quatB(d1/rad,d2/rad,d3/rad,qr) ! elseif(nomen.eq.'K') then ! d1=d1+90. ! if(d1.gt.360.) d1=d1-360. ! d3=90.0-d3 ! if(d3.lt.0.) d3=d3+360. ! call quatB((d1/rad),d2/rad,(d3/rad),qr) ! else ! stop 'cannot deal with that type! ' ! endif c use the new conversion routine based on Altmann if Bunge c else, use longer routine ! print*,'textur5: ng, quat: ',ng,qr do j=1,4 aquat(j,ng)=qr(j) enddo 529 continue 530 continue ! print*,'end ... ng: ',ng ! ngrain = ng-1 ngrain = ng return end c c ________________________________________ c subroutine texture_any(filename,only_one,l_NoWeight) implicit none c fifth version, for use in REXGBS, WTS2POP c only reads in orientations; symm opers in textur4 above c not for use in REX2D c following code taken from LApp68.f for reading in TEXIN ! same as regular WTS version, but arb. extension integer pointx,ng,i,j real evm,fk1(9),nstate character*80 textitl character hkl*3,seclab*5,title*8,idate*8 character eulan*25,nomen,iper*6 character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character*5 label,secid character fname*50,fname2*50 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 logical only_one logical l_NoWeight ! for unweighted input integer numcsl,nsig(50),numsymm,numsamp,ngrain real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50) real aquat(4,2000000),phi1(2000000),capphi(2000000) $ ,phi2(2000000) real grwt(2000000) common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl real tmp_wt,last_phi1,last_capphi,last_phi2 integer count ! used to consolidate orientations real rad character*50 filename ! must match fname2 ! CODE:: rad=57.29578 if(only_one) then open( unit=3 , name= filename , status='old' ) else print*,'Cannot use arbitrary extension with multiple inputs' stop ! open(unit=3,name=fname(1:(pointx+3))//'.wts',status='old') endif read(3,7) textitl 7 format(a) read(3,*) read(3,*) ! don't need the strain etc. c read(3,70) evm,fk1,nstate 70 format(10f7.3,i3) read(3,71) eulan,iper 71 format(a25,49x,a6) nomen=eulan if(nomen.eq.'B') print*, 'Detected Bunge angles' if(nomen.eq.'K') print*, 'Detected Kocks angles' c print 76,textitl 76 format(/ #' TEXIN file= '/1x,a/) ng = 0 count = 0 do 529 i = 1,2000000 if ( mod(i,10000).eq.0 ) print*,'reading line ',i if(l_NoWeight) then read(3,*,end=530,err=530) d1,d2,d3 tmp_wt = 1. ! grwt(ng) = 1. else read(3,*,end=530,err=530) d1,d2,d3,tmp_wt endif ng = ng + 1 c$$$ print*,'i,d1,d2,d3,tmp_wt,ng:',i,d1,d2,d3,tmp_wt,ng if ( d1.ge.0. .and. d1.lt.1.e-3 ) d1 = 1.0e-3 if ( d1.gt.89.999 .and. d1.le.90. ) d1 = 89.999 if ( d1.gt.179.999 .and. d1.le.180. ) d1 = 179.999 if ( d1.gt.359.999 .and. d1.le.360. ) d1 = 359.999 if ( d2.ge.0. .and. d2.lt.1.e-3 ) d2 = 1.0e-3 if ( d2.gt.89.999 .and. d2.le.90. ) d2 = 89.999 if ( d2.gt.179.999 .and. d2.le.180. ) d2 = 179.999 if ( d2.gt.359.999 .and. d2.le.360. ) d2 = 359.999 if ( d3.ge.0. .and. d3.lt.1.e-3 ) d3 = 1.0e-3 if ( d3.gt.89.999 .and. d3.le.90. ) d3 = 89.999 if ( d3.gt.179.999 .and. d3.le.180. ) d3 = 179.999 if ( d3.gt.359.999 .and. d3.le.360. ) d3 = 359.999 ! avoid values on the edges of the space if(i.gt.1) then if(abs(d1 - last_phi1).lt.0.1 .and. & abs(d2 - last_capphi).lt.0.1 .and. & abs(d3 - last_phi2).lt.0.1 ) then ng = ng - 1 count = count + 1 grwt(ng) = grwt(ng) + tmp_wt c$$$ print*,'SAME: i,ng,count,d1,d2,d3' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) ! angles are sufficiently similar to consoliate the point ! with the previous one else grwt(ng) = tmp_wt count = 1 c$$$ print*,'DIFFNT: i,ng,count,d1,d2,d3,wt' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(abs(... else grwt(1) = tmp_wt c$$$ print*,'NG=1: i,ng,count,d1,d2,d3,wt' c$$$ & ,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(ng.gt.1 last_phi1 = d1 last_capphi = d2 last_phi2 = d3 ! save the values ! print*, 'ng, angles.. ',ng,d1,d2,d3 c if((d1.eq.999.).and.(d2.eq.999.).and.(d3.eq.999.)) goto 530 50 format(6f8.2,24f6.2) c read(3,*,end=530) d1,d2,d3,grwt(ng) c .: for old unformatted TEXIN files if(nomen.eq.'B') then phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 elseif(nomen.eq.'K') then phi1(ng)=d1+90. if(phi1(ng).lt.0.) phi1(ng)=phi1(ng)-360. capphi(ng)=d2 phi2(ng)=90.-d3 if(phi2(ng).lt.0.) phi2(ng)=phi2(ng)+360. endif 20 dth=d2 dps=d1 if(nomen.eq.'B') then call quatB(d1/rad,d2/rad,d3/rad,qr) elseif(nomen.eq.'K') then d1=d1+90. if(d1.gt.360.) d1=d1-360. d3=90.0-d3 if(d3.lt.0.) d3=d3+360. call quatB((d1/rad),d2/rad,(d3/rad),qr) else stop 'cannot deal with that type! ' endif c use the new conversion routine based on Altmann if Bunge c else, use longer routine ! print*,'textur5: ng, quat: ',ng,qr do j=1,4 aquat(j,ng)=qr(j) enddo 529 continue 530 continue ! print*,'end ... ng: ',ng ! ngrain = ng-1 ngrain = ng return end c c ________________________________________ c subroutine texture_wquat( pointx ,only_one , l_NoWeight ) implicit none c fifth version, for use in REXGBS, WTS2POP c only reads in orientations; symm opers in textur4 above c not for use in REX2D c following code taken from LApp68.f for reading in TEXIN ! based on textur5 integer pointx,ng,i,j real evm,fk1(9),nstate character*80 textitl character hkl*3,seclab*5,title*8,idate*8 character eulan*25,nomen,iper*6 character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character*5 label,secid character fname*50,fname2*50 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 logical only_one logical l_NoWeight ! for unweighted input integer numcsl,nsig(50),numsymm,numsamp,ngrain real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50) real aquat(4,2000000),phi1(2000000),capphi(2000000),phi2(2000000) real grwt(2000000) common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl real tmp_wt,last_phi1,last_capphi,last_phi2 integer count ! used to consolidate orientations real rad logical last_header character inline*132 integer ind real dum1,dum2,dum3,dum4,dum5 c CODE:: rad = 57.29578 if(only_one) then open(unit=3,name=fname(1:(pointx))//'.wquat',status='old') else open( unit = 3 , name = fname(1:(pointx+3))//'.wquat' $ , status = 'old' ) endif count = 0 last_header = .false. do while( .not.last_header ) read(3,"(a)") inline ! print"('inline:',132a)",inline ind = 0 ind = index( inline , "#" ) ! allows for comment lines (although none in my examples) ! print*,'ind = ',ind, ' last_header ',last_header if ( ind .eq. 1 ) then last_header = .false. count = count + 1 else last_header = .true. ! make it continue unless # found endif enddo print*,'read_wquat: read ',count,' header lines' rewind ( 3 ) if ( count .gt. 0 ) then do i = 1 , count read(3,*) enddo endif ! print*,'read_wquat: reading data lines' nomen = 'B' ! print*, 'Bunge Euler angles assumed' ng = 0 count = 0 do 529 i = 1,9000000 if(l_NoWeight) then ! read(3,*,end=530,err=530) dum1,dum2,dum3,dum4,dum5,d1,d2,d3 read( 3 , * , end=530 ) dum1,dum2,dum3,dum4,dum5,d1,d2,d3 313 format(10f8.3) tmp_wt = 1. else ! read(3,*,end=530,err=530) read( 3 , * , end=530 ) $ qr(4) , qr(1) , qr(2) , qr(3) , tmp_wt endif c$$$ do j = 1 , 3 c$$$ qr(j) = -1. * qr(j) ! active rotation??? c$$$ enddo ng = ng + 1 ! print*,'i,d1,d2,d3,tmp_wt,ng:',i,d1,d2,d3,tmp_wt,ng grwt(ng) = tmp_wt call q2eulB( d1 , d2 , d3 , qr ) print*,'textur_wquat: ng, angles: ', ng , d1 , d2 , d3 ! debug d1 = d1 * rad d2 = d2 * rad d3 = d3 * rad if ( d1.ge.0. .and. d1.lt.1.e-3 ) d1 = 1.0e-3 if ( d1.gt.89.999 .and. d1.le.90. ) d1 = 89.999 if ( d1.gt.179.999 .and. d1.le.180. ) d1 = 179.999 if ( d1.gt.359.999 .and. d1.le.360. ) d1 = 359.999 if ( d2.ge.0. .and. d2.lt.1.e-3 ) d2 = 1.0e-3 if ( d2.gt.89.999 .and. d2.le.90. ) d2 = 89.999 if ( d2.gt.179.999 .and. d2.le.180. ) d2 = 179.999 if ( d2.gt.359.999 .and. d2.le.360. ) d2 = 359.999 if ( d3.ge.0. .and. d3.lt.1.e-3 ) d3 = 1.0e-3 if ( d3.gt.89.999 .and. d3.le.90. ) d3 = 89.999 if ( d3.gt.179.999 .and. d3.le.180. ) d3 = 179.999 if ( d3.gt.359.999 .and. d3.le.360. ) d3 = 359.999 ! avoid values on the edges of the space phi1(ng) = d1 capphi(ng) = d2 phi2(ng) = d3 do j = 1 , 4 aquat(j,ng) = qr(j) enddo 529 continue 530 continue ! print*,'TEXTURE_WQUAT: ' , ng , ' grains read in' ! ngrain = ng-1 ngrain = ng return end c c ________________________________________ c subroutine texture_ctf(pointx,only_one,l_NoWeight) implicit none c fifth version, for use in REXGBS, WTS2POP c only reads in orientations; symm opers in textur4 above c not for use in REX2D c following code taken from LApp68.f for reading in TEXIN ! based on textur5 integer pointx,ng,i,j real evm,fk1(9),nstate character*80 textitl character hkl*3,seclab*5,title*8,idate*8 character eulan*25,nomen,iper*6 character prosa*80,prosa1*80,hh,kk,ll,lper,nomon character*5 label,secid character fname*50,fname2*50 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 logical only_one logical l_NoWeight ! for unweighted input integer numcsl,nsig(50),numsymm,numsamp,ngrain real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50) real aquat(4,2000000),phi1(2000000),capphi(2000000),phi2(2000000) real grwt(2000000) common/a1/ quatsymm,numsymm,quatsamp,numsamp common/a2/quatcsl,numcsl,nsig common/a3/aquat,ngrain,phi1,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl real tmp_wt,last_phi1,last_capphi,last_phi2 integer count ! used to consolidate orientations real rad logical last_header character inline*132 integer ind real dum1,dum2,dum3,dum4,dum5 c CODE:: rad=57.29578 if(only_one) then open(unit=3,name=fname(1:(pointx))//'.ctf',status='old') else open(unit=3,name=fname(1:(pointx+3))//'.ctf',status='old') endif count = 0 last_header = .true. do while(last_header) read(3,"(a)") inline ! print"('inline:',132a)",inline count = count + 1 ind = 0 ind = index(inline,"Euler1") ! print*,'ind = ',ind if ( ind .ne. 0 ) last_header = .false. enddo print*,'read_ctf: read ',count,' header lines' print*,'read_ctf: reading data lines' nomen = 'B' print*, 'Bunge Euler angles in degrees assumed' ng = 0 count = 0 do 529 i = 1,2000000 if(l_NoWeight) then ! read(3,*,end=530,err=530) dum1,dum2,dum3,dum4,dum5,d1,d2,d3 read( 3 , * , end=530 ) dum1,dum2,dum3,dum4,dum5,d1,d2,d3 313 format(10f8.3) tmp_wt = 1. else ! read(3,*,end=530,err=530) read( 3 , * , end=530 ) $ dum1,dum2,dum3,dum4,dum5,d1,d2,d3,tmp_wt endif ng = ng + 1 ! print*,'i,d1,d2,d3,tmp_wt,ng:',i,d1,d2,d3,tmp_wt,ng if ( d1.ge.0. .and. d1.lt.1.e-3 ) d1 = 1.0e-3 if ( d1.gt.89.999 .and. d1.le.90. ) d1 = 89.999 if ( d1.gt.179.999 .and. d1.le.180. ) d1 = 179.999 if ( d1.gt.359.999 .and. d1.le.360. ) d1 = 359.999 if ( d2.ge.0. .and. d2.lt.1.e-3 ) d2 = 1.0e-3 if ( d2.gt.89.999 .and. d2.le.90. ) d2 = 89.999 if ( d2.gt.179.999 .and. d2.le.180. ) d2 = 179.999 if ( d2.gt.359.999 .and. d2.le.360. ) d2 = 359.999 if ( d3.ge.0. .and. d3.lt.1.e-3 ) d3 = 1.0e-3 if ( d3.gt.89.999 .and. d3.le.90. ) d3 = 89.999 if ( d3.gt.179.999 .and. d3.le.180. ) d3 = 179.999 if ( d3.gt.359.999 .and. d3.le.360. ) d3 = 359.999 ! avoid values on the edges of the space if(i.gt.1) then if(abs(d1 - last_phi1).lt.0.1 .and. & abs(d2 - last_capphi).lt.0.1 .and. & abs(d3 - last_phi2).lt.0.1 ) then ng = ng - 1 count = count + 1 grwt(ng) = grwt(ng) + tmp_wt c$$$ print*,'SAME: i,ng,count,d1,d2,d3' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) ! angles are sufficiently similar to consoliate the point ! with the previous one else grwt(ng) = tmp_wt count = 1 c$$$ print*,'DIFFNT: i,ng,count,d1,d2,d3,wt' c$$$ print*,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(abs(... else grwt(1) = tmp_wt c$$$ print*,'NG=1: i,ng,count,d1,d2,d3,wt' c$$$ & ,i,ng,count,d1,d2,d3,grwt(ng) endif ! if(ng.gt.1 last_phi1 = d1 last_capphi = d2 last_phi2 = d3 ! save the values ! print*, 'ng, angles.. ',ng,d1,d2,d3 c if((d1.eq.999.).and.(d2.eq.999.).and.(d3.eq.999.)) goto 530 50 format(6f8.2,24f6.2) c read(3,*,end=530) d1,d2,d3,grwt(ng) c .: for old unformatted TEXIN files phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 20 dth=d2 dps=d1 call quatB(d1/rad,d2/rad,d3/rad,qr) c use the new conversion routine based on Altmann if Bunge c else, use longer routine ! print*,'textur5: ng, quat: ',ng,qr do j=1,4 aquat(j,ng)=qr(j) enddo 529 continue 530 continue ! print*,'end ... ng: ',ng ! ngrain = ng-1 ngrain = ng return end c c ------------------ c subroutine readsod(pointx,only_one) c taken from SOD2VF_17JUN05.F implicit none integer nval,pointx character*80 textitl character prosa*80,prosa1*80,hh,kk,ll,lper character eulan*25,nomen*5 character*5 label,secid character fname*50,one,ten,fname2*50,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(73,20,73),nscale integer ntype,i,j,k,ii,jj,kl,angle_type,nvmx real scalef c real sintens(73,19,19) ! store the intensities from the SOD here integer num_sections logical only_one real rad real oldcapphi,newcapphi,darea(20) integer ngrain,ng real d1,d2,d3,dth,dps,qr(4) real aquat(4,2000000),phi1(2000000) & ,capphi(2000000),phi2(2000000) real grwt(2000000) common/a3/aquat,ngrain,phi1 & ,capphi,phi2 common/a4/ grwt,eulan common/strings/ fname,hkl,seclab,title,idate,textitl c CODE:: 1011 format(a) rad = 57.29578 print*,' ' print*,'In Read_SOD' oldcapphi=0. do 100, i=1,19 newcapphi = (float(i-1)*5.0)+2.5 newcapphi = amin1(newcapphi,90.0) darea(i) = cos(oldcapphi/rad)-cos(newcapphi/rad) oldcapphi = newcapphi 100 continue c print*, 'completed 100 loop; DAREA()= ',darea if(only_one) then open(unit=11,name=fname(1:(pointx))//'.sod',status='old') else open(unit=11,name=fname(1:(pointx+3))//'.sod',status='old') endif ntype = 1 goto 888 c$$$ file=afile(pointx+2:pointx+4) c$$$ if(file.eq.'sod'.or.file.eq.'SOD'.or.file.eq.'smh') then c$$$ ntype=1 c$$$ goto 888 c$$$ endif c$$$ if(file.eq.'cod'.or.file.eq.'COD'.or.file.eq.'cmh') then c$$$ ntype=2 c$$$ goto 888 c$$$ 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.' cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c read in the input data c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c read the first header lines 888 continue c$$$ if(ntype.eq.1) then c$$$ ulimit=1+(4*18/numsamp) ! in SOD files, no of sections varies c$$$ endif c$$$ print*, 'Reading in ',ulimit,' sections based on sample symm.' c$$$ print*, 'If the reading in of a data file fails, then there' c$$$ print*, ' may be a mis-match between the sample symmetry and' c$$$ print*, ' the data file.' c in this version, we will read until we run out of sections! do 725, k = 1,73 read(11,"(a)",err=932,end=932) textitl print*,k,' title=',textitl read(11,1332) nomen,delth,rm,delom,pm,iw,jw 1 ,(iper(i),i=1,3),iavg,ngr,seclab,label 1332 format(a5,4f5.1,5i2,2i5,2a5) if(k.eq.1) then nquads=int((pm+0.01)/90.) ! +.01 to be certain c print*, 'debug: number of lines per PHI value (nquads)= ',nquads if(ntype.eq.2) then c now check the range of phi2 for COD files if(pm.eq.90.and.numsamp.ne.4) then print*, 'Check sample symmetry!' print*, 'Range of azimuth = ',pm print*, 'but no. of sample symmetry operators= ' $ ,numsamp endif if(pm.eq.180.and.numsamp.ne.2) then print*, 'Check sample symmetry!' print*, 'Range of azimuth = ',pm print*, 'but no. of sample symmetry operators= ' $ ,numsamp endif if(pm.eq.360.and.numsamp.ne.1) then print*, 'Check sample symmetry!' print*, 'Range of azimuth = ',pm print*, 'but no. of sample symmetry operators= ' $ ,numsamp endif endif c on the first section, detect the type of Euler angles in use angle_type=1 if(nomen(5:5).eq.'B'.or.nomen(4:4).eq.'B') then print*, 'detected Bunge angles' endif c assume Bunge angles if(nomen(5:5).eq.'K'.or.nomen(4:4).eq.'K') then print*, 'detected Kocks angles' angle_type=2 c$$$ stop 'Can only handle Bunge angles for now!' endif c Kocks angles endif scalef=100./float(iavg) rescale=(abs(scalef-1.0).gt.0.05) if(rescale) print*, '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 do 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 ! ik= 724 enddo ! j= read(11,*,err=932,end=932) c blank line after each section 725 continue num_sections = 73 933 ng = 0 nvmx = nvint(1,1,1) do k = 1,num_sections do 722, j=1,19 do 729, i=1,19 ng = ng+1 if(ntype.eq.1) then ! SOD c sintens(k,j,i)=0.01*nvint(i,j,k) ! switch 1st/3rd indices d1 = 5. * float(k-1) d2 = 5. * float(j-1) d3 = 5. * float(i-1) elseif(ntype.eq.2) then ! COD c sintens(i,j,k)=0.01*nvint(i,j,k) ! end result of subroutine d1 = 5. * float(i-1) d2 = 5. * float(j-1) d3 = 5. * float(k-1) endif grwt(ng) = float(nvint(i,j,k))*darea(j) c correction for size of area is critical! call adjust(d1) call adjust(d2) call adjust(d3) ! ********** this will need to be different for HEX input!!!! c 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) c if(nomen.eq.'B') then if(angle_type.eq.1) then phi1(ng) = d1 capphi(ng) = d2 phi2(ng) = d3 call quatB(d1/rad,d2/rad,d3/rad,qr) c elseif(nomen.eq.'K') then elseif(angle_type.eq.2) then d1=d1+90. if(d1.gt.360.) d1=d1-360. d3=90.0-d3 if(d3.lt.0.) d3=d3+360. phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 call quatB((d1/rad),d2/rad,(d3/rad),qr) else stop 'cannot deal with that type! ' endif c use the new conversion routine based on Altmann if Bunge c else, use longer routine do ik = 1,4 aquat(ik,ng)=qr(ik) enddo print"('angles,grwt= ',4(2x,f10.2))" 1 ,phi1(ng),capphi(ng),phi2(ng),grwt(ng) 729 enddo ! i= 722 enddo ! j= c pause 'check angle values ...' enddo ! k= ngrain = ng write(*,492) nvmx 492 format('READ_SOD: Maximum intensity from SOD =',i7) return c 932 stop 'Input file is not in the right format.' 932 continue if(j.eq.20) num_sections = j-1 if(j.eq.21) num_sections = j-2 ! allow for projection if(j.eq.38) num_sections = j-1 if(j.eq.39) num_sections = j-2 ! allow for projection print*,'j, num_sections,label =',j,num_sections,label if(num_sections.eq.19) then print*,'SOD has 19 sections ' goto 933 elseif(num_sections.eq.37) then print*,'SOD has 37 sections ' goto 933 elseif(num_sections.eq.73) then print*,'SOD has 73 sections ' goto 933 else print*,'SOD has unusable size ' stop endif end c ________________________________ subroutine adjust(angle) real angle c edited 8 x 06, adr c CODE:: if(abs(angle).le.1.e-3) then angle = 1.e-3 endif if(abs(angle-90.).lt.1.e-3) then c angle = 90.001 angle = 89.999 endif if(abs(angle-180.).lt.1.e-3) then c angle = 180.001 angle = 179.999 endif if(abs(angle-270.).lt.1.e-3) then c angle = 270.001 angle = 269.999 endif if(abs(angle-360.).lt.1.e-3) then c angle = 360.001 angle = 359.999 endif return end 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) real p1,p,p2,q(4) double precision c1,c,c2,s1,s,s2 PI=3.14159265 ! 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 print*, 's1,c1,s,c,s2,c2' c print*, s1,c1 c print*, s,c c print*, 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) ! converts quaternion to Bunge Euler angles in radians ! based on Altmann's solution for Euler->quat real qq(4),p1,p,p2 real sum,diff ! CODE:: PI=3.14159265 if((abs(qq(2)).lt.1e-35).and.(abs(qq(1)).lt.1e-35)) then ! diff=pi/4. diff = 0. ! different choice!!! else diff = atan2(qq(2),qq(1)) endif ! ATAN2 is unhappy is both arguments are exactly zero! if((abs(qq(3)).lt.1e-35).and.(abs(qq(4)).lt.1e-35)) then ! sum=pi/4. sum = 0. ! different choice!!! else sum = atan2(qq(3),qq(4)) endif 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$$$ print*, 'q2EulB: quaternion input= ',qq c$$$ print*, 'q2EulB: Bunge angles output= ',p1,p,p2 c$$$ print*, 'q2EulB: angles [degrees]= ' c$$$ & ,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 print*, '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 print*, '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 print*, 'postsymm: qresult ',qresult c return end c 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 donorm(a) real a(3),b(3),rsum,rnorm rsum=0. do 100, i=1,3 rsum=rsum+a(i)**2 100 continue rnorm=sqrt(rsum) do 200, i=1,3 b(i)=a(i)/rnorm a(i)=b(i) 200 continue return end c c ______________________ subroutine qrvec(qq,din,dout) c uses quaternion to rotate a vector from DIN to DOUT ! note the minus sign in the line with the permutation tensor ! we use this version to generate pole figures, which means that ! we are effectively performing the reverse transformation ! from crystal to sample axes (frames) ! cleaned up code vi 08, ADR c note that -q gives the same result as it should (represents same rotation) real qq(4),din(3),dout(3) real t1 , tmp c CODE:: ! print*,'qrvec: qq: ',qq t1 = qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2 do i = 1,3 dout(i) = t1*din(i) tmp = 0. do j = 1,3 tmp = tmp + (qq(j)*din(j)) enddo dout(i) = dout(i) + 2. * qq(i) * tmp do j = 1,3 if(i.ne.j) then do 80, ijk = 1,3 if(i.ne.ijk.and.j.ne.ijk) k=ijk 80 continue kl = j-i if( kl .eq. 2 ) kl = -1 if( kl .eq. -2 ) kl = 1 c poor man's permutation tensor! ! print*,'i,j,k,kl: ',i,j,k,kl dout(i) = dout(i) - (2.*qq(k)*qq(4)*float(kl)*din(j)) endif ! if(i.ne.j) then enddo enddo return end c c _____________________________ c }}}} subroutine invqrvec( qq , din , dout ) c uses quaternion to rotate a vector from DIN to DOUT; c but in the inverse sense so as to give sample to crystal axes c derived from qrvec, above routine c note that -q gives the same result as it should (represents same rotation) real qq(4),din(3),dout(3) real t1 c CODE:: t1 = qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2 do i = 1,3 dout(i) = t1*din(i) tmp = 0. do j = 1,3 tmp = tmp + (qq(j)*din(j)) enddo dout(i) = dout(i) + 2. * qq(i) * tmp do j = 1,3 if(i.ne.j) then do 80, ijk = 1,3 if(i.ne.ijk.and.j.ne.ijk) k=ijk 80 continue kl = j-i if( kl .eq. 2 ) kl = -1 if( kl .eq. -2 ) kl = 1 c poor man's permutation tensor! ! print*,'i,j,k,kl: ',i,j,k,kl dout(i) = dout(i) + (2.*qq(k)*qq(4)*float(kl)*din(j)) ! the ONLY difference between this and QRVEC is the + instead of - above endif ! if(i.ne.j) then enddo enddo return end 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 print*, '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 print*,'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 q1R2>R3 c return end c c _____________________________ c subroutine comquat(qq,qresult) real qq(4,2),qresult(4) c PI=3.14159265 c algorithm for forming resultant quaternion as Q1 x Q2^-1 c where Q2 follows Q1, i.e. is the second rotation; c note change of signs to get inverse of second orientation qresult(1)=qq(1,1)*qq(4,2)-qq(4,1)*qq(1,2) & +qq(2,1)*qq(3,2)-qq(3,1)*qq(2,2) qresult(2)=qq(2,1)*qq(4,2)-qq(4,1)*qq(2,2) & +qq(3,1)*qq(1,2)-qq(1,1)*qq(3,2) qresult(3)=qq(3,1)*qq(4,2)-qq(4,1)*qq(3,2) & +qq(1,1)*qq(2,2)-qq(2,1)*qq(1,2) qresult(4)=qq(4,1)*qq(4,2)+qq(1,1)*qq(1,2) & +qq(2,1)*qq(2,2)+qq(3,1)*qq(3,2) c c write(*,*) 'qresult ',qresult c return end c _____________________________ c c subroutine disquat_hex(qq,index1,index2,qresult,index1a,index2a) c ***************************************** HCP !! c this version for HCP, ADR 25 vi 05 c the shape of the fundamental zone (FZ) for hcp is c a prismatic box, with opt and bottom faces at +/- pi/2*6 = 30 degrees c such that the limits on R3 are independent of R1 and R2 c and square prism faces at a distance 1 (+/- pi/2*4 = 45 degrees) from the center c algorithm for finding disorientation corresponding to c input quaternion, QQ (from COMQUAT): resultant quat. is in QRESULT c indices for symmetry operators are in INDEX1, INDEX2 implicit none real qq(4),qqn(4),qresult(4),quint(4),quintn(4),qmax integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp integer ijk,kjl real switch real pi,rtmp,eps real tan60,sin60,cos60,tan30,sin30,cos30 real rodr(3) logical tests(20) c CODE:: eps = 1.0e-5 PI=3.14159265 tan60 = tan (pi*60./180.) sin60 = sin (pi*60./180.) cos60 = cos (pi*60./180.) tan30 = tan (pi*30./180.) sin30 = sin (pi*30./180.) cos30 = cos (pi*30./180.) index1=0 index2=0 index1a=0 index2a=0 qmax=0. ip=1 jp=1 c do 200, i=1,24 do 200, i=1,12 ! HCP call presymm(qq,i,quint) c do 190, j=1,24 do 190, j=1,12 ! HCP c loop over the two sides of the grain call postsymm(quint,j,quintn) do 185, kjl=1,2 switch=1.0 if(kjl.eq.2) then switch=-1.0 endif do 180, k=1,4 qqn(k)=quintn(k)*switch 180 continue c can take negative of Q because identity = +/-(0,0,0,1) c do ijk=1,2 qresult(1)=qqn(1) qresult(2)=qqn(2) qresult(3)=qqn(3) qresult(4)=qqn(4) iq4neg=ijk if(ijk.eq.2) qresult(4)=-1.*qqn(4) c take the inverse rotation, i.e. apply the switching symmetry call q2rod(rodr,qresult) c$$$ print"('I/J/KJL/Rodr V: ',3i4,3(1x,g10.3))",i,j,kjl,rodr do k = 1,20 tests(k) = .false. enddo if(rodr(2).ge.0.0) then tests(1) = .true. if(rodr(1).le.1.0) then tests(2) = .true. c if((qresult(1)-qresult(2)).lt.1e-5) then c if((qresult(2)-qresult(3)).lt.1e-5) then c if((qresult(3)-qresult(4)).lt.1e-5) then c if(qresult(2).le.qresult(1)) then if(rodr(3).ge.0.) then tests(3) = .true. if(rodr(3).le.tan30) then c these two IFs put the point inside the FZ for R3 tests(4) = .true. rtmp = rodr(1)*tan30 - rodr(2) if(rtmp.ge.0.) then tests(5) = .true. c line from origin at 30 degr rtmp = rodr(1)*cos30 + rodr(2)*sin30 if(rtmp.le.1.) then tests(6) = .true. c 1st sector rtmp = rodr(1)*cos60 $ + rodr(2)*sin60 if(rtmp.le.1.) then tests(7) = .true. c 2nd sector if((qresult(4)-qmax) $ .gt.eps) then tests(8) = .true. qmax=qresult(4) c so now we test to see if we have a smaller angle index1=i index2=j c index3=0 index1a=0 if(ijk.eq.2) index1a=1 index2a=0 if(kjl.eq.2) index2a=1 c record the result if we have found a quat in the fund. zone, and a smaller angle endif endif endif endif endif endif endif endif c$$$ print*,'TESTS: ',(tests(k),k=1,8) c$$$ if(tests(8)) print*,'ALL = TRUE' 210 enddo 185 continue 190 continue 200 continue if(index1.eq.0.or.index2.eq.0) stop 'DISQUAT failed!' call presymm(qq,index1,quint) call postsymm(quint,index2,qqn) c write(*,"('DISQUAT: qqn = ',4f8.3)") qqn do 280, k=1,4 if(index2a.eq.1) then c if we had to take the negative of the Q, then do so c to the answer qresult(k)=qqn(k)*-1.0 c can take negative of Q because identity = +/-(0,0,0,1) else qresult(k)=qqn(k) endif 280 continue if(index1a.eq.1) then qresult(4)=-1.*qresult(4) c restore the 4th component if disorient found for ijk=1 endif c write(*,*) 'DISQUAT: qmax= ',qmax return end c c _____________________________ c c subroutine disquat_ang_hex(qq,angle) c ***************************************** HCP !! c this version for HCP, ADR 25 vi 05 c returning only an angle [degrees], not a quaternion c the shape of the fundamental zone (FZ) for hcp is c a prismatic box, with opt and bottom faces at +/- pi/2*6 = 30 degrees c such that the limits on R3 are independent of R1 and R2 c and square prism faces at a distance 1 (+/- pi/2*4 = 45 degrees) from the center c algorithm for finding disorientation corresponding to c input quaternion, QQ (from COMQUAT): resultant quat. is in QRESULT c indices for symmetry operators are in INDEX1, INDEX2 implicit none real qq(4),qqn(4),qresult(4),quint(4),quintn(4),qmax integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp integer ijk,kjl real switch real pi,rtmp,eps real tan60,sin60,cos60,tan30,sin30,cos30 real rodr(3) logical tests(20) real angle c CODE:: eps = -1.0e-5 PI=3.14159265 tan60 = tan (pi*60./180.) sin60 = sin (pi*60./180.) cos60 = cos (pi*60./180.) tan30 = tan (pi*30./180.) sin30 = sin (pi*30./180.) cos30 = cos (pi*30./180.) index1=0 index2=0 index1a=0 index2a=0 qmax = 0. ip = 1 jp = 1 c do 200, i=1,24 do 200, i=1,12 ! HCP call presymm(qq,i,quint) c do 190, j=1,24 c$$$ do 190, j=1,12 ! HCP c$$$c loop over the two sides of the grain c$$$ call postsymm(quint,j,quintn) do 185, kjl=1,2 switch=1.0 if(kjl.eq.2) then switch=-1.0 endif do 180, k=1,4 c$$$ qqn(k)=quintn(k)*switch qqn(k)=quint(k)*switch 180 continue c can take negative of Q because identity = +/-(0,0,0,1) c do ijk=1,2 qresult(1)=qqn(1) qresult(2)=qqn(2) qresult(3)=qqn(3) qresult(4)=qqn(4) iq4neg=ijk if(ijk.eq.2) qresult(4)=-1.*qqn(4) c take the inverse rotation, i.e. apply the switching symmetry c$$$ call q2rod(rodr,qresult) c$$$ print"('I/J/KJL/Rodr V: ',3i4,3(1x,g10.3))",i,j,kjl,rodr c$$$ c$$$ do k = 1,20 c$$$ tests(k) = .false. c$$$ enddo c$$$ if(rodr(2).ge.0.0) then c$$$ tests(1) = .true. c$$$ if(rodr(1).le.1.0) then c$$$ tests(2) = .true. c if((qresult(1)-qresult(2)).lt.1e-5) then c if((qresult(2)-qresult(3)).lt.1e-5) then c if((qresult(3)-qresult(4)).lt.1e-5) then c if(qresult(2).le.qresult(1)) then c$$$ if(rodr(3).ge.0.) then c$$$ tests(3) = .true. c$$$ if(rodr(3).le.tan30) then c$$$c these two IFs put the point inside the FZ for R3 c$$$ tests(4) = .true. c$$$ c$$$ rtmp = rodr(1)*tan30 - rodr(2) c$$$ if(rtmp.ge.0.) then c$$$ tests(5) = .true. c$$$c line from origin at 30 degr c$$$ rtmp = rodr(1)*cos30 + rodr(2)*sin30 c$$$ if(rtmp.le.1.) then c$$$ tests(6) = .true. c$$$c 1st sector c$$$ rtmp = rodr(1)*cos60 + rodr(2)*sin60 c$$$ if(rtmp.le.1.) then c$$$ tests(7) = .true. c 2nd sector if((qresult(4)-qmax).gt.eps) then c$$$ tests(8) = .true. qmax=qresult(4) c so now we test to see if we have a smaller angle c$$$ index1=i c$$$ index2=j c index3=0 c$$$ index1a=0 c$$$ if(ijk.eq.2) index1a=1 c$$$ index2a=0 c$$$ if(kjl.eq.2) index2a=1 c record the result if we have found a quat in the fund. zone, and a smaller angle endif c$$$ endif c$$$ endif c$$$ endif c$$$ endif c$$$ endif c$$$ endif c$$$ endif c$$$ print*,'TESTS: ',(tests(k),k=1,8) c$$$ if(tests(8)) print*,'ALL = TRUE' enddo 185 continue c$$$ 190 continue 200 continue c$$$ if(index1.eq.0.or.index2.eq.0) stop 'DISQUAT failed!' c$$$ c$$$ call presymm(qq,index1,quint) c$$$ c$$$ call postsymm(quint,index2,qqn) c$$$c write(*,"('DISQUAT: qqn = ',4f8.3)") qqn c$$$ do 280, k=1,4 c$$$ if(index2a.eq.1) then c$$$c if we had to take the negative of the Q, then do so c$$$c to the answer c$$$ c$$$ qresult(k)=qqn(k)*-1.0 c$$$c can take negative of Q because identity = +/-(0,0,0,1) c$$$ else c$$$ qresult(k)=qqn(k) c$$$ endif c$$$ 280 continue c$$$ c$$$ if(index1a.eq.1) then c$$$ qresult(4)=-1.*qresult(4) c$$$c restore the 4th component if disorient found for ijk=1 c$$$ endif c write(*,*) 'DISQUAT_ANG: qmax= ',qmax if(qmax.gt.1.) qmax = 1. angle = 360. / pi * acos(qmax) return end c c _____________________________ c subroutine disquat_orth(qq,thetamin) real qq(4),qqn(4),qresult(4),quint(4),quintn(4),qmax integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp real thetamin c algorithm for finding disorientation corresponding to c input quaternion, QQ (from COMQUAT): resultant quat. is in QRESULT c indices for symmetry operators are in INDEX1, INDEX2 ! edited for orthorhombic crystal symmetry, Jul 07 ! CODE:: index1=0 index2=0 index1a=0 index2a=0 qmax=0. ip=1 jp=1 do 200, i=1,4 ! ortho call presymm(qq,i,quint) do 190, j=1,4 ! ortho c loop over the two sides of the grain call postsymm(quint,j,quintn) do 185, kjl=1,2 switch=1.0 if(kjl.eq.2) then switch=-1.0 endif do 180, k=1,4 qqn(k)=quintn(k)*switch 180 continue c can take negative of Q because identity = +/-(0,0,0,1) do 210, ijk=1,2 qresult(1)=qqn(1) qresult(2)=qqn(2) qresult(3)=qqn(3) qresult(4)=qqn(4) iq4neg=ijk if(ijk.eq.2) qresult(4)=-1.*qqn(4) c take the inverse rotation, i.e. apply the switching symmetry if(qresult(1).ge.0.0) then if(qresult(2).ge.0.0) then if(qresult(3).ge.0.0) then c these IFs ensure that 0<=q1, 0<=q2, 0<=q3 to place it in the fundamental zone ! which is a cube of size=1., between 0 and 1 in all 3 axes if(qresult(4).gt.qmax) then c so now we test to see if we have a smaller angle, which closes the box qmax=qresult(4) index1=i index2=j c index3=0 index1a=0 if(ijk.eq.2) index1a=1 index2a=0 if(kjl.eq.2) index2a=1 c record the result if we have found a quat in the fund. zone, and a smaller angle endif endif endif endif 210 continue 185 continue 190 continue 200 continue if(index1.eq.0.or.index2.eq.0) stop 'DISQUAT failed!' call presymm(qq,index1,quint) call postsymm(quint,index2,qqn) c write(*,"('DISQUAT: qqn = ',4f8.3)") qqn do 280, k=1,4 if(index2a.eq.1) then c if we had to take the negative of the Q, then do so c to the answer qresult(k)=qqn(k)*-1.0 c can take negative of Q because identity = +/-(0,0,0,1) else qresult(k)=qqn(k) endif 280 continue if(index1a.eq.1) then qresult(4)=-1.*qresult(4) c restore the 4th component if disorient found for ijk=1 endif c write(*,*) 'DISQUAT: qmax= ',qmax if (qmax .gt. 1. ) qmax = 1. thetamin = 360. / pi * acos(qmax) return end c c ____________________________ c c subroutine q2rod(d,quat) real d(3),quat(4),rtmp,rnorm,rnorm2 c converts a quaternion to a Rodrigues vector rnorm=sqrt(quat(1)**2+quat(2)**2+quat(3)**2+quat(4)**2) if(rnorm.lt.1e-5) stop 'Q2ROD: too small a quaternion!' if(abs(rnorm-1.0).gt.1.e-3) then write(*,*) 'Q2ROD: quat required normalization' do i=1,4 quat(i)=quat(i)/rnorm enddo endif if(quat(4).gt.1.e-2) then d(1)=quat(1)/quat(4) d(2)=quat(2)/quat(4) d(3)=quat(3)/quat(4) else c we calculate sine(theta/2) from quat1-3 c and scale by tangent/sine in order to make c the system deal with large numbers! c rtmp=sqrt(quat(1)**2+quat(2)**2+quat(3)**2) c this will be close to unity stheta=asin(rtmp) ttheta=tan(stheta) ! and this will be large d(1)=quat(1)*ttheta/stheta d(2)=quat(2)*ttheta/stheta d(3)=quat(3)*ttheta/stheta endif return end c c ------------------- c subroutine euler(a,iopt,nomen,d1,d2,d3,ior,kerr) c Last revision 20nov90 UFK c common a(3,3),grvol(1152),epsga(5),ist1,ist2,sqr3,sqrh,ident(3,3) c SPECIAL VERSION WITHOUT COMMON BLOCK real a(3,3),d1,d2,d3,th,sth,cth,sph,cph,sps,cps,ps,ph,dth,dph,dps character nomen save kor rad=57.29578 c c *** this subroutine calculates the euler angles associated with the c transformation matrix a(i,j) if iopt=1 and viceversa if iopt=2 c *** Note that a is sample (rows) in terms of crystal (columns); c -- opposite of standard g (e.g.Bunge) - this is Canova's c *** Note that in this version, the Euler angles are defined symmetrically: c so that interchanging phi and psi means transposing a. c ("Kocks" nomen: defined going from +X to +Y in both COD and SOD) c *** However, other angle conventions are translated, according to c nomen="K" - Kocks (as internally) -- also sometimes "N"... c "R" - Roe (Psi=psi, Phi=180-phi) c "B" - Bunge (phi1=90+psi, Phi=Theta, phi2=90-phi) c any other - Canova (Theta first, phiC=90+phi, omega=90-psi) c *** Note: only in symm. notation does a point with all Euler angles c between 0 and 90 deg appear in the +x,+y quadrant! c If you want to see an individual point in this quadrant and are: c using Roe, the third angle must be between 90 and 180; c Bunge, first ; c Canova, second . c *** Input and output Euler angles d1,d2,d3 in degrees c goto(5,20),iopt 5 if(abs(a(3,3)) .ge. 0.999) goto 10 th=acos(a(3,3)) sth=sin(th) c if((abs(a(2,3)/sth).lt.1e-35).and. & (abs(a(1,3)/sth).lt.1e-35)) then ps=pi/4. else ps=atan2(a(2,3)/sth,a(1,3)/sth) endif c if((abs(a(3,2)/sth).lt.1e-35).and. & (abs(a(3,1)/sth).lt.1e-35)) then ph=pi/4. else ph=atan2(a(3,2)/sth,a(3,1)/sth) endif ccccc if it bombs out here, both a-components in one arg. are zero c (this should not be possible, but has happened, probably fixed) c ADR: i 01 - above is the fix!! c go to 15 c 10 if((abs(a(1,2)/sth).lt.1e-35).and. & (abs(a(1,1)/sth).lt.1e-35)) then ps=pi/4. else ps=0.5*atan2(a(1,2),-a(1,1)) endif c ph=-ps c The above still have the problem that they give too many c equivalents of the same grain in DIOROUT and density file. Therefore: if(kerr.eq.1.and.kor.ne.ior) then print*,'NOTE: grain',ior,' has Theta < 1 deg.: sometimes problems' print* kor=ior endif c 15 dth=th*rad dph=ph*rad dps=ps*rad d1=dps d2=dth if(nomen.eq.'K'.or.nomen.eq.'N') then d3=dph elseif(nomen.eq.'R') then d3=180.-dph elseif(nomen.eq.'B') then d1=dps+90. d3=90.-dph else d1=dth d2=dph+90. d3=90.-dps endif if(d1.ge.360.) d1=d1-360. if(d3.ge.360.) d3=d3-360. if(d1.lt.0.) d1=d1+360. if(d3.lt.0.) d3=d3+360. return c ************************************* 20 dth=d2 dps=d1 if(nomen.eq.'K') then dph=d3 elseif(nomen.eq.'R') then dph=180.-d3 elseif(nomen.eq.'B') then dps=d1-90. dph=90.-d3 else dth=d1 dph=d2-90. dps=90.-d3 endif ph=dph/rad th=dth/rad ps=dps/rad sph=sin(ph) cph=cos(ph) sth=sin(th) cth=cos(th) sps=sin(ps) cps=cos(ps) a(1,1)=-sps*sph-cph*cps*cth a(2,1)=cps*sph-cph*sps*cth a(3,1)=cph*sth a(1,2)=sps*cph-sph*cps*cth a(2,2)=-cph*cps-sph*sps*cth a(3,2)=sth*sph a(1,3)=sth*cps a(2,3)=sps*sth a(3,3)=cth return end c