C
C==== CAD4(Version 4) Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 25th November 1994
C
      subroutine refrd(itype)
C**** =======================================================================
C     Routine to interpret reflection data from a CAD4.
C     All data is output by calling REFOUT
C**** =======================================================================
      integer itype,i
      real pmt
C
      include 'common.for'
C
C---- In the DIFRAC module for the CAD4 version 4.0 diffractometer, the
C---- UB, UBINV and UBDET variables contain the values reset for a negative
C---- determinant. So for the psi calculation the transformed hkls are used.
C
C---- Take care of complicated Enraf-Nonius logic for anti-reflections
C
      if (nfridl.lt.0) then
        ih = -ih
        ik = -ik
        il = -il
      endif
C
C---- Transform h k l for negative UB
C
      if (trans) then
        ih = -ih
        ik = -ik
        il = -il
      endif
C---- Calculate crystal based psi value
      call psic
C---- Set standard flag
        ptstd = -999
      if ( code6(1:1) .eq. 'N') then
        ptstd = 0
      elseif (code6(1:1) .eq. 'I') 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
      if(ptstd.eq.-999) write(*,
     1 '('' Unlisted standard '',3i5)') ih,ik,il
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---- Check up that time increases
      if ((time + bastim) .le. prvtim) bastim = prvtim - time + 1.0
      time   = time + bastim
      prvtim = time
C---- Set attenuator flag and speed
      if ( spe . ge. 0.0) then
        ptfilt = 0
      else
        ptfilt = 1
      endif
      spe = 16.48/abs(spe)
C
C---- Horizontal detector(aperture) slit opening
C
C---- Get sin theta
      sint  = sin(dar*abs(theta))
      sintl = sint/wavel
      dethz = apta + aptb * sint / sqrt(1.0 - sint*sint)
      if (dethz .lt. apmin) dethz = apmin
      if (dethz .gt. apmax) dethz = apmax
      dethz = rad * dethz / radius
C
C---- Scan width
C
      swidth = 0.66666667*abs(swidth)
C---- Calculation of transformation coefficients
C     Pmt - Peak-with-background measuring time
      pmt = swidth/spe
C     Bgt - left and right background counting time.
      bgt = 0.25 * pmt
      cnp = pmt/swidth
      cnb = pmt
      cnl = bgt
      cnr = bgt
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
      stlmin = min(sintl, stlmin)
      stlmax = max(sintl, stlmax)
C
      call refout
      return
      end
      subroutine dtflrd(jj)
C
C===== Subroutine to read a CAD4 diffractometer output file
C
      integer i,jj, mrefl, trim
      real sinsqa
C
      include 'common.for'
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
C---- Initialisations
C
      ptscl = 1
      firstr = .false.
      reoren = .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 ,readonly)
CVMS 3 ,readonly)
CUNX 3 )
C****
C     Interpretation of the .dat file
C****
 1000 continue
      line = ' '
      read(datafl, '(a64)', end = 2000) line
C---- Shift the buffer if necessary
      if (slidc4.gt.0) then
        do 1010 i = 80-slidc4,1,-1
          line(i+slidc4:i+slidc4)=line(i:i)
 1010   continue
      else if (slidc4.lt.0) then
        do 1015 i = 1-slidc4,80
          line(i+slidc4:i+slidc4) = line(i:i)
 1015   continue
      endif
C---- Find the line type
      if ( line(1:5) .eq. '     ') then
C---- Do nothing - we need to ignore entirely blank lines
C
C---- Type 21 line.
C
      elseif (line(1:3) .eq. ' 21') then
        numref = 0
        parclp = 1
        read(line,'(45x,f4.0,3f4.1)')
     1    radius, slit, apmin, apmax
C
C---- Type 22 line.
C
      elseif (line(1:3) .eq. ' 22') then
        parclp = parclp + 2
        read(line,'(25x,2f5.2,f2.0)')
     1    apta, aptb, scnrat
C
C---- Type 23 line.
C
      elseif (line(1:3) .eq. ' 23') then
        parclp = parclp + 4
C
C---- Type 24 line.
C
      elseif (line(1:3) .eq. ' 24') then
        parclp = parclp + 8
C
C---- Type 25 line.
C
      elseif (line(1:3) .eq. ' 25') then
        parclp = parclp + 16
C
C---- Type 26 line.
C
      elseif (line(1:3) .eq. ' 26') then
        read(line,'(34x,2f10.7)')
     1    sinsqa, cosalp
        sinalp = sqrt(sinsqa)
        parclp = parclp + 32
C
C---- Type 0 line. Get the check reflection indices
C
      elseif (line(1:3) .eq. '   ') then
        if (((line(15:15).eq.'H').or.(line(15:15).eq.'N')).and.
     1     ((line(19:19).eq.'I').or.(line(19:19).eq.'F'))) then
          numref = numref + 1
          read(line,'(5X,3I3)')
     1      href(numref), kref(numref), lref(numref)
        endif
C
C---- Type 31 line
C
      elseif (line(1:3) .eq. ' 31') then
        parclp = parclp + 64
        read(line,'(5x,3f9.6,2x,3f9.6)')
     1   ub(1,1),ub(1,2),ub(1,3),ub(2,1),ub(2,2),ub(2,3)
C
C---- Type 32 line
C
      elseif (line(1:3) .eq. ' 32') then
        parclp = parclp + 128
        if (parclp .eq. 255) firstr = .true.
        if (parclp .eq. 447) then
          reoren = .true.
          parclp = 255
        endif
        read(line,'(5x,3f9.6,2x,3f9.6)')
     1   ub(3,1),ub(3,2),ub(3,3),flam1,flam2,beta
C
C---- Only allow the reflection lines if control data present
C
      elseif (parclp .eq. 255) then
C
C---- Type 1 line
C
        if (line(1:3) .eq. '  1') then
          read(line,'(3x,i6,3i5,1x,a6,f7.2,f4.0,f6.0,f7.0,f6.0)')
     1      mrefl, ih, ik, il, code6, psi, spe, backl, peak, backr
C
C---- Type 2 line
C
        elseif (line(1:3) .eq. '  2') then
          read (line,'(3x,i6,f8.3,3f9.3,f7.3,f7.0,i3)')
     1      i, theta, phi, omega, kappa, swidth, time, nfridl
          if ( i .eq. mrefl) then
            if (firstr .or. reoren) call initlz
            call refrd(0)
            mrefl = 0
          endif
        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 CAD4 diffractometer parameter output file
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C
C===== On the Enraf-Nonius CAD4, with this programme installed on the
C      diffractometer itself, this subroutine could be made to read the
C      <compid>.CRY and the GONCA0.DAT files for information on the machine.
C      This is not yet implemented in this version but will come soon.
C
      subroutine prflrd
C
      return
      end
      subroutine choos
C
C----- Set up fixed characteristics for CAD4 diffractometer
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 an Enraf-Nonius diffractometer
C
      difnam = 'E'
      device = 'Enraf-Nonius CAD4 4circle diffractometer'
      softwr = 'Enraf-Nonius CAD4 version 4.0           '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = 'pck$:[cristal]address.lst                         '
CNORMAL addnam = 'address.lst                                       '
      addnam = 'address.lst                                       '
C
C---- File extension for data file name
      datxtn = '.dat'
C
C---- Set dihedral angle
      fmnhed = 90.0
C
C---- Default wavelength pointer for cell dimension determination, alpha 1
      cdefwv = 1
C
C---- Set the buffer shift for .dat files of CAD4. Enraf-Nonius has
C---- chosen to have a very funny ASCII file structure for their *.dat
C---- files. A consequence of this is that when the files are shifted from
C---- machine to machine the character get shifted one place to the right or
C---- left. The variable slidc4 sets the shift to get back to the standard
C---- position usually obtained with a DIFRAC running on the machine which
C---- produced the file.
C
      slidc4 = 0
C
C---- Options for Quesan
C
      return
      end
      subroutine initlz
C
C**** Subroutine to interpret a complete 21-26 information in memory
C     for the CAD4.
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
C---- Do not do the rest of the calculations for a reorientation only
C
      if (reoren) then
C---- Reset reoren so that INITLZ will not be called again
        reoren = .false.
      else
C
C---- Convert the scan type parameter
C
        scnrat = scnrat/3.0
        mstnn = 3
        if(abs(scnrat).lt. 0.1)       mstnn = 1
        if(abs(scnrat - 2.0).lt. 0.1) mstnn = 2
C---- Wavelength
        weigh1 = 0.6666667
        weigh2 = 0.3333333
        weigh3 = 0.0000000
        flam3  = 0.0
        wavel = weigh1*flam1 + weigh2*flam2
C---- Vertical detector opening
        detvt = rad*slit/radius
C
C---- Set Firstr to indicate that INITLZ does not need calling again.
C
        firstr = .false.
      endif
C
      return
      end
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Enraf-Nonius CAD4 machine.
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; D.Schwarzenbach and H.D.Flack (1992)
C     J. Appl. Cryst. 25, 69. 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 CAD-4 corresponds to case (a) of section 3) "Dealing with
C     Manufacturers" design of DSHDF since 2-theta, omega, phi and kappa
C     are found on the diffractometer output file as well as UB each time
C     it changes. The CAD4 is thus a DSHDF type I diffractometer and hence
C     W = 1 according to DSHDF equation (11). (W and 1 being 3x3 matrices).
C
C     The CAD4 uses the kappa geometry described in detail in the CAD4
C     User's Manual, Chapter II, "CAD4 Geometry", hereafter CAD4MII.
C     At page 6 of CAD4MII one finds the matrix expression for the setting
C     of the goniometer to angles omega, kappa and phi as:
C         Z(omega) * Y(-alpha) * Z(Kappa) * Y(+alpha) * Z(phi).
C     Thus for the kappa geometry, equation (9) of DSHDF becomes:
C
C    T        T     -1
C Psi  = s * u  * UB   * Z(-phi) * Y(-alpha) * Z(-kappa) * Y(+alpha) * Z(-omega)
C
C               -2     T     -1     -1T
C        where s   =  u  * UB   * UB    * u    and
C
C                (  1          0           0          )
C     X(alph) =  (  0          cos(alph)   sin(alph)  )
C                (  0         -sin(alph)   cos(alph)  )
C
C                (  cos(beta)  0          -sin(beta)
C     Y(beta) =  (  0          1           0          )
C                (  sin(beta)  0           cos(beta)  )
C
C                (  cos(gamm)  sin(gamm)   0          )
C     Z(gamm) =  ( -sin(gamm)  cos(gamm)   0          )
C                (  0          0           1          )
C
C
      integer  i,j
      real r(3), s(3), sinq, cosq
      include 'common.for'
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)
      sinq = sin(dar*phi)
      cosq = cos(dar*phi)
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *y(-alph)
      s(1) =  r(1)*cosalp - r(3)*sinalp
      s(2) =  r(2)
      s(3) =  r(1)*sinalp + r(3)*cosalp
C---- *z(-kappa)
      sinq = sin(dar*kappa)
      cosq = cos(dar*kappa)
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *y(+alph)
      s(1) =  r(1)*cosalp + r(3)*sinalp
      s(2) =  r(2)
      s(3) = -r(1)*sinalp + r(3)*cosalp
C---- *z(-omega)
      sinq = sin(dar*(omega-theta))
      cosq = cos(dar*(omega-theta))
      r(1) =  s(1)*cosq + s(2)*sinq
      r(2) = -s(1)*sinq + s(2)*cosq
      r(3) =  s(3)
C---- *W (DSHDF type I, W is a unit matrix, nothing to be done here
C---- Measurements at negative theta correspond to 
C          (1    0    0)
C        W (0   -1    0)
C          (0    0   -1)
C     so we must take care of that
      if (theta .lt. 0.0) then
        r(2) = -r(2)
        r(3) = -r(3)
      endif
C
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
C==== WRITE(12,'(1x,3i5,6f10.2)')
C====1 ih,ik,il, omega, 2*theta, kappa, phi, psi, psicry
C
      return
      end
