C                                                                               
C==== Rigaku Module of DIFRAC by J. Hester, National Institute for Inorganic
C==== Materials, Tsukuba, Japan               
C==== Version 25th Feb. 1997                                                 
C==== Users must modify the settings in subroutines prflrd and choos.
C                                                                               
      subroutine refrd(itype)                                                   
C**** ==========================================================================
C     Routine to interpret reflection data from a Rigaku AFC-5 FOS File                   
C     All data is output by calling REFOUT                                      
C**** ==========================================================================
      integer itype,i                                                           
      real pmt                                                                  
      logical oka                                                               
C                                                                               
      include 'common.for'                                                      
C                                                                               
C---- In the DIFRAC module for the Rigaku FOS, the           
C---- UB, UBINV and UBDET variables contain the values used on the              
C---- diffractometer even if the determinant terms out to be negative.          
C---- So for the psi calculation the diffractometer hkls are used and           
C---- are transformed for output.  It has to be done this way because           
C---- the RIGAKU sofware does not output the setting angles used for each               
C---- reflection.                                                               
C                                                                               
C---- Calculate crystal based psi value                                         
      call psic                                                                 
C---- Transform output h k l for negative UB                                    
      if (trans) then                                                           
        ih = -ih                                                                
        ik = -ik                                                                
        il = -il                                                                
      endif                                                                     
C---- Set standard flag                                                         
        ptstd = -999                                                            
        if ( itype .eq. 1 ) then                                          
          ptstd = 0                                                             
          stdseq = 0                
C
C---- Deal with standards
C                                                 
        elseif ( itype .eq. 0 ) then
          if (numref .gt. 0) then
           do 100 i = 1, numref                                                
              if ((ih.eq.href(i)).and.(ik.eq.kref(i)).and.                      
     1          (il.eq.lref(i))) ptstd = i                                      
  100        continue               
          endif
C----
C---- If first appearance of standards, put them into the array
C----
          if (ptstd .le. 0) then
            numref = numref + 1                                             
            href(numref) = ih
            kref(numref) = ik
            lref(numref) = il
            ptstd = numref                                                         
          endif                                                                 
        endif                                                                   
C---- Totalize the standards                                                    
        if(ptstd.gt.0) then                                                     
          avint(ptstd) = avint(ptstd) + peak                                    
          sdint(ptstd) = sdint(ptstd) + peak*peak                               
          count(ptstd) = count(ptstd) + 1.0                                     
        endif                                                                   
C                                                                               
C---- Horizontal detector(aperture) slit opening                                
C                                                                               
C                                                                               
C---- Scan width  
C---- Get sin theta needed for scan transition angle and swidth                 
        if (abs(theta) .le. sct) then                                           
C---- omega scans                                                               
          scnrat = 0.0                                                       
        else                                                                    
C---- omega - 2-theta scans                                                     
          scnrat = 2.0                                                       
        endif                                                                   
        swidth = 0.02*real(int(                                              
     1   50.0*(swda + swdb* sint / sqrt(1.0 - sint*sint))))                             
C                                                                               
C---- Calculation of transformation coefficients
C                                
C---- Stationary background measuring technique                                 
          cnp = timp/swidth                                                
          cnb = timp                                                       
          cnl = bgt                                                       
          cnr = bgt                                                       
          backsn = 0                                                            
C---- Count the reflections                                                     
        nrefl = nrefl + 1                                                       
C---- Keep min/max values of h, k, l                                            
        hmin = min(hmin, ih)                                                    
        kmin = min(kmin, ik)                                                    
        lmin = min(lmin, il)                                                    
        hmax = max(hmax, ih)                                                    
        kmax = max(kmax, ik)                                                    
        lmax = max(lmax, il)                                                    
C---- Keep min/max values of sintl                                              
C       
      sintl = sint/wavel                                                        
      stlmin = min(sintl, stlmin)                                               
      stlmax = max(sintl, stlmax)                                               
      
C---- Keep record of total time
C
      time=time+timp+bgt+bgt
C                                                            
      call refout                                                             
      return                                                                    
      end                                                                       
C                                                                               
      subroutine dtflrd(jj)                                                     
C                                                                               
C===== Subroutine to read a Rigaku diffractometer output file 
C     Reference: AFC-5, AFC-6 FOS Manual Number MJ201GP7 p 75                    
C                                                                               
      integer i,jj, mrefl, trim                                                 
C---- These variables for integer versions of the time
      integer itimp,ibgt,ntnet  
C                                  
      real sinsqa
C                                                                               
      include 'common.for'                                                      
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C                                                                               
C---- Initialisations                                                           
C                                                                               
      ptscl = 1                                                                 
      firstr = .false.                                                          
      mrefl = 0                                                                 
C                                                                               
C---- Get the data file name and open it                                        
C                                                                               
      open( unit = datafl,                                                      
     1 file = compid(1:trim(compid))//datxtn,                                   
     2 status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'                
     3 )                                                                        
CVAX 3 ,readonly)                                                               
CUNX 3 )                                                                        
C****                                                                           
C     Interpretation of the .dat file                                           
C****                                                                           
 1000 continue                                                                  
      line = ' '                                                                
      read(datafl, '(a64)', end = 2000) line                                    
C---- Find the line type                                                        
      if ( line(1:4) .eq. 'CELL') then
C---- Read in the unit cell                                         
      read(datafl, '(3f9.4,3f9.3,f9.2)', end = 2000)
     1  cell(1),cell(2),cell(3),cell(4),cell(5),cell(6),volum                                    
      read(datafl, '(3f9.4,3f9.3,f9.2)', end = 2000)                                           
     1  ecell(1),ecell(2),ecell(3),ecell(4),ecell(5),ecell(6),evolum                   
C                                                                               
C---- ATTF line.                              
C                                                                               
      elseif (line(1:4) .eq. 'ATTF') then                                        
      read(datafl, '(4f12.7)', end = 2000)
     1  attfac(1),attfac(2),attfac(3),attfac(4)                                             
C                                                                               
C---- UB line                                                              
C                                                                               
      elseif (line(1:2) .eq. 'UB') then                                        
        read(datafl, '(f12.7)', end = 2000) wavel                                   
        read(datafl, '(3f12.7)', end = 2000)                                   
     1   ub(1,1),ub(1,2),ub(1,3)                       
        read(datafl, '(3f12.7)', end = 2000)                                   
     1   ub(2,1),ub(2,2),ub(2,3)                        
        read(datafl, '(3f12.7)', end = 2000)                                   
     1   ub(3,1),ub(3,2),ub(3,3)
C
C---- Beginning of data lines
C     
        elseif (line(1:6) .eq. 'RDFILE') then
         firstr = .true.                       
         read(datafl,'(i5)',end=2000) mrefl
C                                                                               
C---- Type 0,1,2,3 line                                                               
C                                                                               
        elseif (line(1:2) .eq. ' 1' .or. line(1:2) .eq. ' 2' .or. 
     $      line(1:2) .eq. ' 3' .or. line(1:2) .eq. ' 0') then                           
          read(line,'(2x,3i3,20x,i2,i7,2i6,2i5,i2)')               
     1      ih, ik, il, ptfilt, ipeak, ibackl, ibackr, 
     2      itimp, ibgt, iscan
C---- Convert to scan time and background time in seconds
        timp=real(itimp)/100.
        bgt=real(ibgt)/100.
        peak=real(ipeak)
        backl=real(ibackl)
        backr=real(ibackr)
        if (firstr) call initlz
        if (line(1:2) .eq. ' 0') then 
          call refrd(0)
        else 
          call refrd(1)
        endif             
      endif                                                                     
      go to 1000                                                                
C                                                                               
C---- Finish off this data file                                                 
C                                                                               
 2000 continue                                                                  
      close (unit = datafl, status = 'KEEP')                                    
      return                                                                    
      end                                                                       
C                                                                               
C===== Subroutine to read a Rigaku diffractometer parameter output file           
C                                                                               
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C                                                                               
C                                                                               
      subroutine prflrd
      real id
      dimension id(510)             
      include 'common.for'                 
C     
C---  The following lines taken from Rigaku FOS software `file.sub'
C                                                                          
C     OPEN(50,FILE='FOR50.DAT',STATUS='OLD',FORM='UNFORMATTED',
C    1     ACCESS='DIRECT',RECL=780)
C     READ(50,REC=4,ERR=900) (ID(K),K=101,490)
C
C     Get scan speed
C
      spe=4
C
C     Get scan width parameters
C
      scanm=0.0
      swda=0.7
      swdb=0.0
      sct=90.0
C     Get apertures from user
      detvt = 2.5/315.                                                     
      dethz = 2.5/315.
      return                                                        
C
C---- Read error
C
  900 write(*,'(''file 50 read error'')')
      stop
      end
C                                                          
C----- Set up fixed characteristics for WKA diffractometer                  
C
      subroutine choos
C                                                                               
      include 'common.for'                                                      
C                                                                               
C==== Set Logical Unit Numbers for DIFRAC                                       
C                                                                               
C     The program uses extensively the standard input and output files          
C     as accessed by the FORTRAN instructions read(*,.. and write(*,.. .        
C                                                                               
C     datafl is the unit used to read the diffractometer data files             
C     scfsfl is the unit used to write the standardized output file             
C     parmfl is the unit used to read the diffractometer parameter files        
C     addrfl is the unit used to read the address list file                     
C     qaofl  is the unit used to write/read an existing saved answer file       
C     qanfl  is the unit used to read/write a new saved answer file             
C                                                                               
      datafl =  9                                                               
      scfsfl = 10                                                               
      parmfl = 11                                                               
      addrfl = 12                                                               
      qaofl  = 13                                                               
      qanfl  = 14                                                               
C                                                                               
C---- Tell them its a Huber diffractometer                             
C                                                                               
      difnam = 'H'                                                              
      device = 'Huber_Four_Circle_Diffractometer        '                       
      softwr = 'Modified_Rigaku_AFC5_Program            '                       
      cellpg = '?                                       '
      dettyp = 'scintillation counter                   '
C                                                                               
C---- Tell the programme the name of the address file                           
      addnam = 'address.lst'                               
C                                                                               
C---- Set extension for data file name                                          
      datxtn = '_87.DAT'                                                       
C                                                                               
C---- Set dihedral angle                                                        
      fmnhed = 0.0                                                              
C                                                                               
C---- Default wavelength pointer for cell dimension determination, alpha 1      
      cdefwv = 1                                                                
C                                                                               
C---- Options for Quesan                                                        
      ask(3) = .false.                                                          
C
C---  Number of filters
C
      nfilt = 4
      do 100 i=1,4 
       nfltcd(i)=i-1
  100 continue
C                                                                               
      return                                                                    
      end                                                                      
      subroutine initlz                                                         
C                                                                               
C**** Subroutine to interpret initial information             
C     for Rigaku AFC-5 FOS.                                                             
C****                                                                           
      integer i,j,k                                                             
      real ubdet                                                                
C                                                                               
      include 'common.for'                                                      
C                                                                               
C---- Get UB**-1                                                                
C                                                                               
      ubinv(1,1) = +ub(2,2)*ub(3,3) - ub(3,2)*ub(2,3)                           
      ubinv(2,1) = -ub(2,1)*ub(3,3) + ub(3,1)*ub(2,3)                           
      ubinv(3,1) = +ub(2,1)*ub(3,2) - ub(3,1)*ub(2,2)                           
      ubdet=ub(1,1)*ubinv(1,1)+ub(1,2)*ubinv(2,1)+ub(1,3)*ubinv(3,1)            
      ubinv(1,1) = ubinv(1,1)/ubdet                                             
      ubinv(2,1) = ubinv(2,1)/ubdet                                             
      ubinv(3,1) = ubinv(3,1)/ubdet                                             
      ubinv(1,2) = (-ub(1,2)*ub(3,3) + ub(3,2)*ub(1,3))/ubdet                   
      ubinv(2,2) = (+ub(1,1)*ub(3,3) - ub(3,1)*ub(1,3))/ubdet                   
      ubinv(3,2) = (-ub(1,1)*ub(3,2) + ub(3,1)*ub(1,2))/ubdet                   
      ubinv(1,3) = (+ub(1,2)*ub(2,3) - ub(2,2)*ub(1,3))/ubdet                   
      ubinv(2,3) = (-ub(1,1)*ub(2,3) + ub(2,1)*ub(1,3))/ubdet                   
      ubinv(3,3) = (+ub(1,1)*ub(2,2) - ub(2,1)*ub(1,2))/ubdet                   
C                                                                               
C**** Calculate (approximate) direct space cell dimensions from UB matrix       
C                                                                               
C---- Get real metric in GM from ubinv                                          
      do 400 i = 1,3                                                            
      do 400 k = 1,3                                                            
        gm(i,k) = 0.0                                                           
      do 400 j = 1,3                                                            
        gm(i,k) = gm(i,k) + ubinv(i,j)*ubinv(k,j)                               
  400 continue                                                                  
C---- Get cell dimensions from ub                                               
      do 401 i = 1, 3                                                           
       cellub(i) = sqrt(gm(i,i))                                                
 401  continue                                                                  
      do 402 i = 1,3                                                            
       j = 1 + i/3                                                              
       k = i -j + 2                                                             
       cellub(7-i) = rad*acos(gm(j,k)/(cellub(j)*cellub(k)))                    
 402  continue                                                                  
       volub=abs(1.0/ubdet)                                                     
C---- Make angles 90 and 120 and cell lengths equal                             
      do 404 i = 4, 6                                                           
        if(abs(cellub(i)-90.) .lt. 0.4) cellub(i) = 90.                         
 404  continue                                                                  
      if(abs(cellub(6)-120.) .lt. 0.4) cellub(6) = 120.                         
      if(abs(cellub(1)-cellub(2)).lt.0.01) then                                 
        if(abs(cellub(1)-cellub(3)).lt.0.01) then                               
          cellub(1) = (cellub(1)+cellub(2)+cellub(3))/3.0                       
          cellub(2) = cellub(1)                                                 
          cellub(3) = cellub(1)                                                 
        else                                                                    
          cellub(1) = 0.5*(cellub(1)+cellub(2))                                 
          cellub(2) = cellub(1)                                                 
        endif                                                                   
      endif                                                                     
C                                                                               
C**** Check for negative Ub mark to transform UB and hkl if necessary           
C                                                                               
      if (ubdet .gt. 0.0) then                                                  
        trans = .false.                                                         
      else                                                                      
        trans = .true.                                                          
        write(*,'('' ***** Negative ub matrix.'',                               
     1  '' Transforming UB and h k ls*****'')')                                 
        ubdet = - ubdet                                                         
        do 100 i = 1, 3                                                         
        do 100 j = 1, 3                                                         
          ub(i,j)    = - ub(i,j)                                                
          ubinv(i,j) = - ubinv(i,j)                                             
 100    continue                                                                
        do 105 i = 1, numref                                                    
          href(i) = - href(i)                                                   
          kref(i) = - kref(i)                                                   
          lref(i) = - lref(i)                                                   
 105    continue                                                                
      endif                                                                     
C---- Wavelength                                                                
      flam1 = wavel                                                             
      flam2 = 0.0                                                               
      flam3 = 0.0                                                               
      weigh1 = 1.0                                                              
      weigh2 = 0.0                                                              
      weigh3 = 0.0                                                              
      tubtar = '?'
      wl0nam = '?'
      wl1nam = '?'
      wl2nam = '?'
      do 201 i = 1, navtab                                                      
        if (abs(wavtab(1,i)-wavel).le.0.005) then                               
          wavel  = wavtab(1,i)                                                  
          flam1  = wavtab(2,i)                                                  
          flam2  = wavtab(3,i)                                                  
          weigh1 = wavtab(4,i)                                                  
          weigh2 = wavtab(5,i)                                                  
          tubtar = wavnam(1,i)
          wl0nam = wavnam(2,i)
          wl1nam = wavnam(3,i)
          wl2nam = wavnam(4,i)
        endif                                                                   
 201  continue                                                                  
        if      (cdefwv.eq.1) then                                              
          cwave = flam1                                                         
          clrnam = wl1nam
        else if (cdefwv.eq.2) then                                              
          cwave = flam2                                                         
          clrnam = wl2nam
        else if (i.ne.0) then                                                   
          cwave = wavel                                                         
          clrnam = wl0nam
        endif                                                                   
C                                                                               
C---- Set Firstr to indicate that INITLZ does not need calling again.           
C                                                                               
        firstr = .false.                                                                     
C                                                                               
      return                                                                    
      end                                                                       
C
      subroutine psic                                                           
C                                                                               
C**** Routine to calculate the crystal-based azimuth of the reflection          
C     for Rigaku FOS reading from Huber diffractometer.                    
C                                                                               
C     The equations used for the calculation of the crystal-based azimuth       
C     are given in "On the definition and Practical Use of Crystal-based        
C     Azimuthal Angles" by D.Schwarzenbach and H.D.Flack (1989)                 
C**** J. Appl. Cryst. 22, 601--605, hereafter DSHDF.                        
C     The definition of the reference vector u is given in equation (3) of      
C     DSHDF and equation (9) gives the psi-crystal.                             
C                                                                               
C     The Rigaku/Huber pair corresponds to case (c) of section 3) "Dealing with            
C     Manufacturers" design of DSHDF since 2-theta, omega, phi and kappa        
C     are not found on the diffractometer output file but the current UB is     
C     available.  The Rigaku software assumes a         
C     DSHDF type I diffractometer and hence according to DSHDF equation (11):  
C                                                                               
C                   (  1   0   0)                                               
C             W =   (  0   1   0)                                               
C                   (  0   0   1).                                              
C                         
C     The Rigaku software implements clockwise = positive motion when viewed from
C     the direction of the source (chi) or from above (phi, omega, theta).  The Rigaku
C     software model has z vertical and y parallel to the incident beam (a 90 rotation 
C     compared to DSHDF).  
C
C     Thus, the rotation matrices required for DSHDF equation (9) are, using DSHDF 
C     definitions for x,y,z laboratory axes:
C                                                      
C                (  1          0           0          )                         
C     X(chi)  =  (  0          cos(chi)   -sin(chi)   )                         
C                (  0          sin(chi)    cos(chi)   )                         
C                                                                               
C                (  cos(phi)  -sin(phi)    0          )                         
C     Z(phi)  =  (  sin(phi)   cos(phi)    0          )                         
C                (  0          0           1          )                         
C                                                                               
C                (  cos(omega)  -sin(omega)    0      )                         
C     Z(omega) = (  sin(omega)   cos(omega)    0      )                         
C                (  0            0             1      )                         
C                                                                               
      integer  i,j                                                              
      real x(3), r(3), s(3), d, d2, p, p2, temp                                 
      real sinchb, coschb, sinphb, cosphb                                       
      real sinom, cosom, sinchi, coschi, sinphi, cosphi                         
      real sindph, cosdph, sinpsi, cospsi, cost                                 
      include 'common.for'                                                      
C---- Calculate the angles for bisecting position                               
      do 205 i = 1 ,3                                                           
        x(i) = (ub(i,1)*real(ih)+ub(i,2)*real(ik)+ub(i,3)*real(il))
     1          *wavel               
 205  continue                                                                  
C---- Get p and d                                                               
      p2 = x(1)*x(1)+x(2)*x(2)                                                  
      d2 = p2 + x(3)*x(3)                                                       
      p  = sqrt(p2)                                                             
      d  = sqrt(d2)                                                             
C---- Get sint and cost, theta and two theta                                    
      sint = 0.5*d                                                              
      if( sint.gt.1.0 ) then                                                    
        write(*,'('' Sint > 1.0 for reflection '',3i5)') ih, ik, il             
        return                                                                  
      endif                                                                     
      cost = sqrt(1.0 - 0.25*d2)                                                
      theta = rad*atan2( sint, cost)                                            
      twoth = 2.0*theta                                                         
C==== Bisecting Position Angle Calculation
C---- This calculation based on routine CLANG from the RD5 control
C---- source code
C
      chb=atan2(x(3),p)
      chb=chb/dar
      if(chb .gt. 89.995) then
        phb=0.0
      else
        phb = atan2(x(2),x(1))
      endif
      sinchb = sin(chb*dar)                                                          
      coschb = cos(chb*dar)                                                         
      sinphb = sin(phb)                                                         
      cosphb = cos(phb)                                                         
C---- Is the reflection set out of the bisecting position                       
      if (psi.ne.0.0) then                                                      
        sinpsi = sin(dar*psi)                                                   
        cospsi = cos(dar*psi)                                                   
C     Make omega. Temp will contain tan(omega). Force cos(omega) positive       
C     so abs(omega)=< 90 degrees.                                               
        temp = sinpsi*coschb/sinchb                                             
        cosom = 1.0/sqrt(1.0+temp*temp)                                         
        sinom = cosom*temp                                                      
C     Make chi. Temp will contain scale factor to ensure that                   
C     sinchb and sinchi have the same sign                                      
        sinchi = sinchb                                                         
        coschi = cosom*coschb*cospsi                                            
        temp = 1.0 / sqrt (sinchi*sinchi + coschi*coschi)                       
        sinchi = sinchi * temp                                                  
        coschi = coschi * temp                                                  
C     Make delta phi (dph). Temp will contain scale to ensure that              
C     sindph and sinom have opposite signs                                      
        sindph = -sinom                                                         
        cosdph =  cosom*coschi                                                  
        temp = 1.0 / sqrt (sindph*sindph + cosdph*cosdph)                       
        sindph = sindph * temp                                                  
        cosdph = cosdph * temp                                                  
C     Make phi from phi(bisecting) and delta phi                                
        sinphi = sinphb*cosdph + cosphb*sindph                                  
        cosphi = cosphb*cosdph - sinphb*sindph                                  
      else                                                                      
C---- Bisecting position                                                        
        sinom = 0.0                                                             
        cosom = 1.0                                                             
        sinchi = sinchb                                                         
        coschi = coschb                                                         
        sinphi = sinphb                                                         
        cosphi = cosphb                                                         
      endif                                                                     
C                                                                               
C---- Form the reference vector u: equation (3) of DSHDF                        
C                                                                               
      if ((ih.eq.ik).and.(ik.eq.il)) then                                       
        r(1) = real(ih)                                                         
        r(2) = real(-ih)                                                        
        r(3) = 0.0                                                              
      else                                                                      
        r(1) = real(ik-il)                                                      
        r(2) = real(il-ih)                                                      
        r(3) = real(ih-ik)                                                      
      endif                                                                     
C---- u(transpose)*ubinv                                                        
      do 110 i = 1 ,3                                                           
        s(i) = 0.0                                                              
        do 105 j = 1, 3                                                         
          s(i) = s(i) + r(j)*ubinv(j,i)                                         
 105    continue                                                                
 110  continue                                                                  
C---- *z(-phi)                                                                  
      r(1) =  s(1)*cosphi - s(2)*sinphi                                         
      r(2) =  s(1)*sinphi + s(2)*cosphi                                         
      r(3) =  s(3)                                                              
C---- *x(-chi)                                                                  
      s(1) =  r(1)                                                              
      s(2) =  r(2)*coschi - r(3)*sinchi                                         
      s(3) =  r(2)*sinchi + r(3)*coschi                                         
C---- *z(-omega)                                                                
      r(1) =  s(1)*cosom - s(2)*sinom                                           
      r(2) =  s(1)*sinom + s(2)*cosom                                           
      r(3) =  s(3)                                                              
C---- *W (DSHDF Type I diffractometer)                                         
      r(1) = r(1)                                                              
      r(3) = r(3)                                                              
C---- Crystal based azimuth                                                     
      psicry = rad*atan2(r(1),r(3))                                             
C                                                                               
C==== r(2) = r(2)/ SQRT(r(1)*r(1) + r(2)*r(2) + r(3)*r(3))                      
C==== WRITE(12,'(/'' psi(2) '', e15.6)') r(2)                                   
C==== Look at them to see if they agree                                         
      omega = rad*ATAN2(sinom,cosom) + theta                                    
      chi   = rad*ATAN2(sinchi,coschi)                                          
      phi   = rad*ATAN2(sinphi,cosphi)                                          
      WRITE(12,'(1x,3i5,6f10.2)')                                               
     1 ih,ik,il, omega, twoth, chi, phi, psi, psicry
      write(12,'(1x,3f10.5)') x(1),x(2),x(3)                            
C                                                                               
      return                                                                    
      end                                                                       



