C
C==== Niv411 Module of DIFRAC by E. Blanc, University of Lausanne
C==== Version 25th November 1994
C
      subroutine refrd(itype)
C**** ==========================================================================
C     Routine to interpret reflection data from a NICOLET P21
C     All data is output by calling REFOUT
C**** ==========================================================================
      integer itype,i
      real pmt
      logical oka
C
      include 'common.for'
C---- Get theta
      theta = twoth/2.0
C
C---- In the DIFRAC module for the Siemens P21 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---- 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 (nrefl .gt. 0) then
        ptstd = 0
        stdseq = 0
      else
        do 10 i = 1, numref
          if ((ih .eq. href(i)) .and. (ik .eq. kref(i)) .and.
     $        (il .eq. lref(i))) ptstd = i
   10   continue
      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---- Check up that time increases
      time = time*60.0
      if ((time+bastim) .lt. prvtim) bastim = prvtim-time+1.0
      time = time+bastim
      prvtim = time
C---- Set attenuator flag and speed
      if (swidth .ge. 0) then
        ptfilt = 0
      else
        ptfilt = 1
      endif
C
C---- Horizontal and vertical detector(aperture) slit opening
C
C
C---- Scan width
C
      swidth = abs(swidth)
C---- Calculation of transformation coefficients
C---- Scanned background measuring technique
C     Pmt - Total measuring time
      pmt = swidth/spe
C     Bgt - left and right background counting time.
      bgt = 0.5 * ratio * pmt
      cnp = 1.0/spe
      cnb = pmt
      cnl = bgt
      cnr = bgt
      backsn = 1
      nrefl = iabs(nrefl)
C---- Keep one eye on the filter codes that have been used
      oka = .false.
      do 2010 i = 1, nfilt
        if (ptfilt .eq. nfltcd(i)) oka = .true.
 2010 continue
      if (.not. oka) then
        nfilt = nfilt + 1
        nfltcd(nfilt) = ptfilt
      endif
C
      call refout
      return
      end
C
      subroutine dtflrd(jj)
C
C===== Subroutine to read a NICOLET diffractometer output file
C
      integer jj, i, trim
C
      include 'common.for'
C
      ptscl=1
C
C---- Open input file _REF.DCF)
C
      open( unit=datafl,  
     1  file = compid(1:trim(compid))//datxtn, status = 'OLD', 
     2  access = 'SEQUENTIAL', form = 'UNFORMATTED'
     3  )
CVMS 3  ,readonly)
CUNX 3  )
      if (mstnn .eq. 6) then
        nstep = nstepw
      else
        nstep = 96
      endif
      call initlz
C
C Reads & write all the data file
C
   20 continue
        read(datafl, end = 9999, err = 30) nrefl, ih, ik, il, twoth,
     $  omega, phi, chi, swidth, spe, backl, peak, backr, firl, sir,
     $  time, (profil(i), i = 17, nstep+16)
        call refrd
      go to 20
C
C Error encountered
C
   30 write(*,'(''### Error during reading datafile ###'')')
C
C Finish !!!
C
 9999 close (datafl, status = 'KEEP')
      return
      end
C
C===== Subroutine to read a SYNTEX diffractometer parameter output file
C
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
C===============================================================================
C
C Le diffractometre SYNTEX cree, pour chaque acquisition de donnees, deux
C fichiers, l'un contenant les parametres de la prise de donnee (ouverture
C des detecteurs, ...), et l'autre contenant les intensites mesurees et
C les profils.
C
C Cette subroutine lit les parametres
C
C===============================================================================
      subroutine prflrd

C     Print anything
      character cdum*15
C     File name
      character dcfnme*10
C     Open file type
      character dcopen*10
      character fsd*2
C     User hkl file name
      character hklnme*10
      character lscons*1
C     Anode metal
      character target*2
      character xrysrc*2

C     Loop indices
      integer i, j, trim
      integer ianum
C     Check flag for reference reflections
      integer chknum(12)
      integer dcmode
C     Direction to collect data
      integer idir
C     Exposure time in min
      integer iexp
      integer ifit
      integer iflip
      integer ifrdl
      integer imode
      integer indx
      integer iord(4)
      integer iordr(3)
C     Nb of page in log file
      integer ipage
      integer iprfl
      integer iprint
      integer iprmp
      integer iprof
      integer ipstp
      integer iqus
      integer ircnt
      integer irstp
      integer iscnt
      integer scnhgb
      integer ismpl
      integer ispn
      integer itotl
      integer itype
C     Indices from file flag
      integer iuser
      integer ivect(3)
C     Standard every iwhen refl
      integer iwhen
      integer iwl
      integer iwntme
      integer izm
      integer jtemp
C     Indices of next reflection
      integer ks(3)
      integer ksp(6)
C     Lattice nb
      integer latt
      integer maxfnd
      integer maxp
      integer mlts
      integer msoln
      integer mtape
C     Number of standard
      integer ncheck
      integer ncstp
      integer ngrp
C     Number of reflections left
      integer norefl
      integer nrcnt
C     Integer part of centering reflections
      integer crefl(5, 50)
C     Number of centering reflections
      integer nrf
      integer nsoln
      integer ntrc
C
C     Filter  WARNING- Printed as character
      integer*2 filter
      integer*2 he, hp
      integer*2 iword
      integer*2 lstype
      integer*2 lt
C
C     Used to print systematic absences
      logical print
C
C     Attenuator
      real atnbem
      real axlth(6)
      real baschi
      real basphi
      real baseom
      real basett
      real cfor
C     Chi collision protection limits (NOT READ BY NICOLET.SUB)
      real chibg, chism
C     Chi zero (NOT READ BY NICOLET.SUB)
      real chizro
C     Check value for standard reflections
      real chk(12)
      real chlf
      real cscan
      real ctim
C     Date of data collection
      real date
      real del(4)
      real dev(6)
C     Dispersion (see formula)
      real disp
      real dm
      real dpsi
      real edgmon
C     Low bound used in Wyckoff scans
      real fintm
C     Low & hich correction bounds
      real flcut, fhcut
      real flim1(3), flim2(3)
      real fm3
      real fnl
      real fnono(3)
C     Speed on each circle (NOT READ BY NICOLET.SUB)
      real frat(4)
      real fsig
C     Fastest scan rate
      real fsrt
      real fstt
      real ftpe
      real ofsmic
      real omg1, omg2
C     Omega collision protection limits (NOT READ BY NICOLET.SUB)
      real omgbg, omgsm
C     Omega zero (NOT READ BY NICOLET.SUB)
      real omgzro
      real omspd
      real p3d(10)
      real pcent
C     Phi collision protection limits (NOT READ BY NICOLET.SUB)
      real phibg, phism
      real phis
      real pltim
      real prat
      real pscan
      real psi1, psi2
      real ptim
      real qm
      real ratem
C     Fp part of centering reflections
      real refl(5, 50)
      real refps
      real resmx
C     Range around peak (see formula)
      real rng1, rng2
C     Systematic absences
      real sahkl(12)
      real sahhl(12)
      real sahk(8), sahl(8), sakl(8)
      real sah(6), sak(6), sal(12)
      real sahhlh(2), sahhlr(2), sahblh(2)
      real sclplt
      real slow
C     Slowest scan rate
      real slrt
C     Scan rate
      real srat
C     Estimated time left
      real timest
C     Temporary storage
      real tmp1, tmp2
C     Tolerances (NOT READ BY NICOLET.SUB)
      real tol(4)
C     2 theta collision protection limits (NOT READ BY NICOLET.SUB)
      real tthbg, tthsm
C     2 theta zero (NOT READ BY NICOLET.SUB)
      real tthzro
C     P3 version number
      real versnb
      real xfrat(3)
      real xm
      real ym
      real zm
C
      include 'common.for'
C
C****
C     Interpretation of the parameters' file
C****
C
C---- Open input files (_PAR.DCF and _REF.DCF)
C
      open( unit = parmfl, 
     1  file = compid(1:trim(compid))//'_PAR.DCF', status = 'OLD',
     2  form = 'UNFORMATTED')
C
C---- Read title
C
      do 1 i = 1, 60
        title(i:i) = ' '
   1  continue
      read(parmfl) title(1:60), date, versnb
C
C---- Read cell parameters
C
      read(parmfl) (cell(i), i = 1, 6),
     $  (ecell(i), i = 1, 6), ((ub(i, j), j = 1, 3), i = 1, 3),
     $  ((ubinv(i, j), j = 1, 3), i = 1, 3), tthzro, omgzro, chizro,
     $  target, wavel, flam1, flam2, beta, volum, evolum,
     $  edgmon, filter
      if (versnb .ge. 4.12) then
        read(parmfl, err = 90) iprmp, iqus, mstnn, xrysrc, fsd, hp,
     1    he, ispn
      else
        read(parmfl, err = 90) iprmp, iqus, mstnn, fnl, ftpe, ispn
      endif
      go to 91
   90 versnb = 3.43
      mstnn = 0
      ispn = 0
   91 continue
      if (versnb .ge. 3.99) then
        read(parmfl) jtemp, lscons, lstype
      else if (versnb .ne. 3.43) then
        read(parmfl) jtemp
      else
        read(parmfl) iword
        jtemp = iword
      endif
      read(parmfl) fnono, frat, tthsm, tthbg, omgsm,
     $  omgbg, phism, phibg, chism, chibg, imode, lt, basphi,
     $  baschi, basett, baseom, ofsmic
      if ((tthsm .eq. 0.0) .and. (tthbg .eq. 0.0) .and.
     $  (omgsm .eq. 0.0) .and. (omgbg .eq. 0.0)) then
        omgbg = 58.0
        omgsm = 312.0
        tthbg = omgbg + fnono(3)
        tthsm = amod(omgsm - fnono(3), 360.0)
      endif
      read(parmfl) itype, pscan, pltim, iprint,
     $  sclplt, ratem
      read(parmfl) iwl, iord, tol, del, iscnt, cscan,
     $  phis, iprof
      read(parmfl) zm, qm, dm, ptim, xm, ym, prat
      read(parmfl) ivect, ifit, indx, axlth, dev, fm3,
     $  izm, msoln, nsoln
      read(parmfl) omg1, omg2, omspd, ianum
C
C---- Read data collection parameters
C
      read(parmfl) stlmin, stlmax, hmin, kmin, lmin, hmax, kmax, lmax,
     $  ntrc, latt, iuser, nrefl, ks, ipage, iexp, disp, rng1, rng2,
     $  ratio, srat, slrt, fsrt, chlf, cfor, iordr, norefl, timest,
     $  idir
      norefl = iabs(norefl)
      if (versnb .ge. 4.01) then
        read(parmfl, err = 92) ncheck, iwhen, itotl, mtape, ismpl,
     $    iprfl, maxp, dcfnme, dcopen, iwntme, dcmode, hklnme
      else
        read(parmfl, err = 92) ncheck, iwhen, itotl, mtape, ismpl,
     $    iprfl, maxp, dcfnme, dcopen, iwntme, dcmode
      endif
      go to 93
   92 dcfnme = 'SCRATCH'
      dcopen = 'NEW'
      iwntme = 0
   93 continue
      read(parmfl) psi1, psi2, dpsi
      read(parmfl) pcent, ircnt, nrcnt
      read(parmfl) ksp
      read(parmfl) p3d
      read(parmfl) resmx, refps, mlts
      read(parmfl) ncstp, ctim, ngrp, slow
C
C---- Read dead time parameters
C
      read(parmfl) fdeadt, flcut, fhcut, atnbem
      fdeadt = fdeadt * 10**6
      read(parmfl) scnhgb
      read(parmfl) nstepw, ipstp, irstp, fintm,
     $  ifrdl, iflip
C
C---- Read systematic absences
C
      read(parmfl) (sahkl(i), i = 1, 12)
      read(parmfl) (sahk(i), i = 1, 8),
     $  (sahl(i), i = 1, 8), (sakl(i), i = 1, 8)
      read(parmfl) (sah(i), i = 1, 6),
     $  (sak(i), i = 1, 6), (sal(i), i = 1, 12)
      read(parmfl) (sahhl(i), i = 1, 12)
      read(parmfl) (sahhlh(i), i = 1, 2),
     $  (sahblh(i), i = 1, 2)
      read(parmfl) (sahhlr(i), i = 1, 2)
      tmp1 = flim1(1)
      tmp2 = flim1(2)
      read(parmfl, err = 94) maxfnd, fstt, fsig, flim1, flim2, xfrat
      go to 95
   94 if (maxfnd .le. 0) maxfnd = 25
      xfrat(3) = flim1(1)
      flim2(3) = flim1(2)
      flim1(3) = flim1(1)
      flim1(2) = tmp2
      flim1(1) = tmp1
   95 continue
C
C---- Standard reflexions
C
      do 96 i = 1, 10
        href(i) = 0
        kref(i) = 0
        lref(i) = 0
   96 continue
      numref = 0
      do 97 i = 1, 12
        read(parmfl, err = 97) ih, ik, il, chk(i), chknum(i)
        if ((numref .lt. 10) .and. (chknum(i) .gt. 0)) then
          numref = numref+1
          href(numref) = ih
          kref(numref) = ik
          lref(numref) = il
        endif
   97 continue
C
C---- Centering reflexions
C
      nrf = 50
      do 98 i = 1, nrf
        read(parmfl, end = 9999) (crefl(j, i), j = 1, 5),
     $    (refl(j, i), j = 1, 5)
   98 continue
      if (iordr(1) .ne. 0) then
        do 99 i = 1, 3
          iordr(i) = i
   99   continue
      endif
      if (iabs(idir) .ne. 1) idir = 1
      if (edgmon .eq. 0.0) edgmon = 3.1
 9999 close(parmfl, status = 'KEEP')
C
C---- Print parameter file for user
C
      write(*,'('' Do you wish to display the contents of '',
     1 ''the parameter file (Y,N)? '')')
      call fread(000080,'XNNNNNNNNN')
      inchar = charin(1)(1:1)
      if (inchar .eq. 'Y') then
        write(*, '('' Title '', a60, 3X, ''Software version '', F5.2)')
     1    title, versnb
C---- Write cell parameters
        write(*, '(/'' Cell parameters & orientation matrix'')')
        write(*, '(5X, F11.4, 2X, ''a in Angstroems'')') cell(1)
        write(*, '(5X, F11.4, 2X, ''Esd of a'')') ecell(1)
        write(*, '(5X, F11.4, 2X, ''b in Angstroems'')') cell(2)
        write(*, '(5X, F11.4, 2X, ''Esd of b'')') ecell(2)
        write(*, '(5X, F11.4, 2X, ''c in Angstroems'')') cell(3)
        write(*, '(5X, F11.4, 2X, ''Esd of c'')') ecell(3)
        write(*, '(5X, F11.4, 2X, ''alpha in degrees'')') cell(4)
        write(*, '(5X, F11.4, 2X, ''Esd of alpha'')') ecell(4)
        write(*, '(5X, F11.4, 2X, ''beta  in degrees'')') cell(5)
        write(*, '(5X, F11.4, 2X, ''Esd of beta'')') ecell(5)
        write(*, '(5X, F11.4, 2X, ''gamma in degrees'')') cell(6)
        write(*, '(5X, F11.4, 2X, ''Esd of gamma'')') ecell(6)
        write(*, '(5X, F11.4, 2X, ''Volume in Angstroems**3'')')
     1    volum
        write(*, '(5X, F11.4, 2X, ''Esd of volume'')') evolum
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 1)'')')
     1    ub(1, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 2)'')')
     1    ub(1, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (1, 3)'')')
     1    ub(1, 3)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 1)'')')
     1    ub(2, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 2)'')')
     1    ub(2, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (2, 3)'')')
     1    ub(2, 3)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 1)'')')
     1    ub(3, 1)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 2)'')')
     1    ub(3, 2)
        write(*, '(5X, F11.6, 2X, ''Orientation matrix, term (3, 3)'')')
     1    ub(3, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 1)'')') ubinv(1, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 2)'')') ubinv(1, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (1, 3)'')') ubinv(1, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 1)'')') ubinv(2, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 2)'')') ubinv(2, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (2, 3)'')') ubinv(2, 3)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 1)'')') ubinv(3, 1)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 2)'')') ubinv(3, 2)
        write(*, '(5x, f11.4, 2x, ''Inverse orientation matrix, '',
     1    ''term (3, 3)'')') ubinv(3, 3)
        cdum = '    Constrained'
        if ((lscons .ne. 'C') .and. (lscons .ne. 'c'))
     1    cdum = '  Unconstrained'
        write(*, '(1X, A15, 2X, ''Orientation matrix'')') cdum
C---- Write radiation & filter stuff
        write(*, '(/'' Radiation and filter'')')
        write(*, '(14X, A2, 2X, ''Anode element'')') target
        write(*, '(5X, F11.6, 2X, ''Mean wavelength'')') wavel
        write(*, '(5X, F11.6, 2X, ''Kalpha1 wavelength'')') flam1
        write(*, '(5X, F11.6, 2X, ''Kalpha2 wavelength'')') flam2
        write(*, '(5X, F11.6, 2X, ''Kbeta'')') beta
        if (filter .ne. 0) then
          write(*, '(14X, A2, 2X, ''Absorption filter element'')')
     1      filter
        endif
        write(*, '(14X, I2, 2X, ''Wavelength pointer'')') iwl
C---- Write lattice type & systematic absences
        write(*, '(/'' Lattice type & systematic absences'')')
        cdum = '              ?'
        go to (111, 112, 113, 114, 115, 116, 117, 118) latt+1
  111   cdum = '              P'
        go to 119
  112   cdum = '              A'
        go to 119
  113   cdum = '              B'
        go to 119
  114   cdum = '              C'
        go to 119
  115   cdum = '              F'
        go to 119
  116   cdum = '              I'
        go to 119
  117   cdum = '           Robv'
        go to 119
  118   cdum = '           Rrev'
  119   write(*, '(1X, A15, 2X, ''Lattice type'')') cdum
        print = .false.
        do 130 i = 1, 12
          if (sahkl(i) .ne. 0.0) print = .true.
  130   continue
        if (print) then
          do 135 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hkl'')')
     1        sahkl(i)
  135     continue
        end if
        print = .false.
        do 140 i = 1, 8
          if (sahk(i) .ne. 0.0) print = .true.
  140   continue
        if (print) then
          do 145 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hk'')')
     1        sahk(i)
  145     continue
        end if
        print = .false.
        do 150 i = 1, 8
          if (sahl(i) .ne. 0.0) print = .true.
  150   continue
        if (print) then
          do 155 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hl'')')
     1        sahl(i)
  155     continue
        end if
        print = .false.
        do 160 i = 1, 8
          if (sakl(i) .ne. 0.0) print = .true.
  160   continue
        if (print) then
          do 165 i = 1, 8
            write(*, '(5X, F11.2, 2X, ''Systematic absence type kl'')')
     1        sakl(i)
  165     continue
        end if
        print = .false.
        do 170 i = 1, 6
          if (sah(i) .ne. 0.0) print = .true.
  170   continue
        if (print) then
          do 175 i = 1, 6
            write(*, '(5X, F11.2, 2X, ''Systematic absence type h'')')
     1        sah(i)
  175     continue
        end if
        print = .false.
        do 180 i = 1, 6
          if (sak(i) .ne. 0.0) print = .true.
  180   continue
        if (print) then
          do 185 i = 1, 6
            write(*, '(5X, F11.2, 2X, ''Systematic absence type k'')')
     1        sak(i)
  185     continue
        end if
        print = .false.
        do 190 i = 1, 12
          if (sal(i) .ne. 0.0) print = .true.
  190   continue
        if (print) then
          do 195 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type l'')')
     1        sal(i)
  195     continue
        end if
        print = .false.
        do 200 i = 1, 12
          if (sahhl(i) .ne. 0.0) print = .true.
  200   continue
        if (print) then
          do 205 i = 1, 12
            write(*, '(5X, F11.2, 2X, ''Systematic absence type hhl'')')
     1        sahhl(i)
  205     continue
        end if
        print = .false.
        do 210 i = 1, 2
          if (sahhlh(i) .ne. 0.0) print = .true.
  210   continue
        if (print) then
          do 215 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhl_H'')') sahhlh(i)
  215     continue
        end if
        print = .false.
        do 220 i = 1, 2
          if (sahblh(i) .ne. 0.0) print = .true.
  220   continue
        if (print) then
          do 225 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhbl_H'')') sahblh(i)
  225     continue
        end if
        print = .false.
        do 230 i = 1, 2
          if (sahhlr(i) .ne. 0.0) print = .true.
  230   continue
        if (print) then
          do 235 i = 1, 2
            write(*, '(5x, f11.2, 2x, ''Systematic absence type '',
     1        ''hhl_R'')') sahhlr(i)
  235     continue
        end if
C---- Write data collection parameters
        write(*, '(/'' Data collection parameters'')')
        cdum = '        Unknown'
        if (mstnn .eq. 1) cdum = '          Omega'
        if (mstnn .eq. 2) cdum = '   Theta:2Theta'
        if (mstnn. eq. 6) cdum = '          Wycoff'
        write(*, '(1X, A15, 2X, ''Data collection type'')') cdum
        cdum = '              ?'
        go to (251, 252, 250, 253, 254, 255) imode+3
  251     cdum = ' Parallel(user)'
          go to 250
  252     cdum = 'Bisecting(user)'
          go to 250
  253     cdum = '      Bisecting'
          go to 250
  254     cdum = '      Bisecting'
          go to 250
  255     cdum = '       Parallel'
  250   continue
        write(*, '(1X, A15, 2X, ''Data collectiom mode '')') cdum
        write(*, '(14x, i2, 2x, ''Low temperature mode '',
     1    ''(0=no low temp)'')') lt
        write(*, '(6X, A10, 2X, ''Data collection file name'')') dcfnme
        write(*, '(10X, A6, 2X, ''Open mode'')') dcopen
C---- Data collection bounds
        write(*, '(/'' Data collection bounds'')')
        write(*, '(5X, F11.2, 2X, ''2 theta min'')') stlmin
        write(*, '(5X, F11.2, 2X, ''2 theta max'')') stlmax
        write(*, '(12X, I4, 2X, ''H min'')') hmin
        write(*, '(12X, I4, 2X, ''K min'')') kmin
        write(*, '(12X, I4, 2X, ''L min'')') lmin
        write(*, '(12X, I4, 2X, ''H max'')') hmax
        write(*, '(12X, I4, 2X, ''K max'')') kmax
        write(*, '(12X, I4, 2X, ''L max'')') lmax
        write(*, '(12x, i4, 2x, ''Number of last reflection '',
     1    ''proceeded'')') nrefl
        write(*, '(12X, I4, 2X, ''H of last reflection proceeded'')')
     1    ks(1)
        write(*, '(12X, I4, 2X, ''K of last reflection proceeded'')')
     1    ks(2)
        write(*, '(12X, I4, 2X, ''L of last reflection proceeded'')')
     1    ks(3)
        write(*, '(8X, I8, 2X, ''Number of reflections proceeded'')')
     1    norefl
        write(*, '(8X, I8, 2X, ''Total exposure time in min'')') iexp
        write(*, '(5X, F11.4, 2X, ''Remaining time in hours'')') timest
        write(*, '(12X, I4, 2X, ''H of first reflection'')') ksp(1)
        write(*, '(12X, I4, 2X, ''K of first reflection'')') ksp(2)
        write(*, '(12X, I4, 2X, ''L of first reflection'')') ksp(3)
        write(*, '(12X, I4, 2X, ''H of last  reflection'')') ksp(4)
        write(*, '(12X, I4, 2X, ''K of last  reflection'')') ksp(5)
        write(*, '(12X, I4, 2X, ''L of last  reflection'')') ksp(6)
        write(*, '(14X, I2, 2X, ''H variation order'')') iordr(1)
        write(*, '(14X, I2, 2X, ''K variation order'')') iordr(2)
        write(*, '(14X, I2, 2X, ''L variation order'')') iordr(3)
        write(*, '(14X, I2, 2X, ''Direction of index variation'')')
     1    idir
        write(*, '(14X, I2, 2X, ''Index restriction code'')') ntrc
        if (iuser .eq. 1) then
          write(*, '(6x, a10, 2x, ''Selected hkl data collection '',
     1      ''file name'')') hklnme
        endif
        if ((dpsi .ne. 0.0) .or. (psi1 .ne. 0.0)) then
          write(*, '(5X, F11.2, 2X, ''Low  psi bound'')') psi1
          write(*, '(5X, F11.2, 2X, ''High psi bound'')') psi2
          write(*, '(5X, F11.4, 2X, ''Delta psi angle'')') dpsi
        endif
C---- Data collection parameters
        write(*, '(/'' Data collection parameters (part 2)'')')
        write(*, '(5X, F11.4, 2X, ''Multiplicative factor for range'')')
     1    disp
        write(*, '(5X, F11.4, 2X, ''Below Kalpha1'')') rng1
        write(*, '(5X, F11.4, 2X, ''Above Kalpha2'')') rng2
        if (srat .eq. 0.0) then
          write(*, '(5X, F11.4, 2X, ''Slowest scan speed'')') slrt
          write(*, '(5X, F11.4, 2X, ''Fastest scan speed'')') fsrt
        else
          write(*, '(5X, F11.4, 2X, ''Fixed scan speed'')') srat
        endif
        write(*, '(5X, F11.6, 2X, ''Background/scan measuring time'')')
     1    ratio
        write(*, '(5X, F11.4, 2X, ''Low  I or |F|'')') chlf
        write(*, '(5X, F11.4, 2X, ''High I or |F|'')') cfor
        if (mstnn .eq. 6) then
          write(*, '(12X, I4, 2X, ''Total steps in Wycoff scan'')')
     1      nstepw
          write(*, '(12X, I4, 2X, ''Calculated steps'')') ipstp
          write(*, '(12X, I4, 2X, ''Extra steps'')') irstp
          write(*, '(5X, F11.2, 2X, ''Intensity'')') fintm
          write(*, '(14X, I2, 2X, ''Collect Friedel pair'')') ifrdl
          write(*, '(14X, I2, 2X, ''Start with hkl or -h-k-l'')') iflip
        endif
        write(*, '(14x, i2, 2x, ''1=profiled prescan, 2=intensity, '',
     1    ''<0=Profiled postscan'')') ispn
C---- Diffractometer parameters
        write(*, '(/'' Diffractometer parameters'')')
        write(*, '(12X, I4, 2X, ''Number of steps'')') ncstp
        write(*, '(12X, I4, 2X, ''Number of group'')') ngrp
        write(*, '(5X, F11.4, 2X, ''Second/step'')') ctim
        write(*, '(5X, F11.4, 2X, ''Min drive speed'')') slow
        write(*, '(5X, F11.4, 2X, ''Dead time in microsec'')') fdeadt
        write(*, '(5X, F11.2, 2X, ''Correction low  bound'')') flcut
        write(*, '(5X, F11.2, 2X, ''Correction high bound'')') fhcut
        if (atnbem .ne. 0.0) then
          write(*, '(5X, F11.6, 2X, ''Attenuation'')') atnbem
        end if
        write(*, '(5x, f11.6, 2x, ''2 theta zero '',
     1    ''(calibration value) '')') tthzro
        write(*, '(5x, f11.6, 2x, ''omega   zero '',
     1    ''(calibration value) '')') omgzro
        write(*, '(5x, f11.6, 2x, ''chi     zero '',
     1    ''(calibration value) '')') chizro
        write(*, '(5X, F11.2, 2X, ''2 Theta collision low  bound'')')
     1    tthsm
        write(*, '(5X, F11.2, 2X, ''2 Theta collision high bound'')')
     1    tthbg
        write(*, '(5X, F11.2, 2X, ''Omega   collision low  bound'')')
     1    omgsm
        write(*, '(5X, F11.2, 2X, ''Omega   collision high bound'')')
     1    omgbg
        write(*, '(5X, F11.2, 2X, ''Phi     collision low  bound'')')
     1    phism
        write(*, '(5X, F11.2, 2X, ''Phi     collision high bound'')')
     1    phibg
        write(*, '(5X, F11.2, 2X, ''Chi     collision low  bound'')')
     1    chism
        write(*, '(5X, F11.2, 2X, ''Chi     collision high bound'')')
     1    chibg
        write(*, '(5X, F11.4, 2X, ''2 Theta fastest speed'')') frat(1)
        write(*, '(5X, F11.4, 2X, ''Omega   fastest speed'')') frat(2)
        write(*, '(5X, F11.4, 2X, ''Phi     fastest speed'')') frat(3)
        write(*, '(5X, F11.4, 2X, ''Chi     fastest speed'')') frat(4)
        write(*, '(5X, F11.4, 2X, ''Microscope offset'')') ofsmic
        write(*, '(5x, f11.2, 2x, ''2 theta home value for '',
     1    ''centering crystal in beam'')') basett
        write(*, '(5x, f11.2, 2x, ''Omega   home value for '',
     1    ''centering crystal in beam'')') baseom
        write(*, '(5x, f11.2, 2x, ''Phi     home value for '',
     1    ''centering crystal in beam'')') basphi
        write(*, '(5x, f11.2, 2x, ''Chi     home value for '',
     1    ''centering crystal in beam'')') baschi
C---- Write standard reflections
        write(*, '(/'' Standard reflections'')')
        write(*, '(14X, I2, 2X, ''Number of check reflections'')')
     1    numref
        write(*, '(12X, I4, 2X, ''Interval between checks'')') iwhen
        write(*, '(12X, I4, 2X, ''Interval between checks in min'')')
     1    iwntme
        write(*, '(5X, F11.4, 2X, ''Recenter if intensity below'')')
     1    pcent
        write(*, '(/''   h   k   l'', 5x, ''low bound'', 3x, ''check''
     1    )')
        do 300 i = 1, numref
          write(*, '(3(1X, I3), 4X, F10.2, 5X, I3)')
     1      href(i), kref(i), lref(i), chk(i), chknum(i)
  300   continue
C---- Write centering reflections
        write(*, '(/'' Centering reflections'')')
        write(*, '(14x, i2, 2x, ''Centering: 1st pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(1)
        write(*, '(14x, i2, 2x, ''Centering: 2nd pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(2)
        write(*, '(14x, i2, 2x, ''Centering: 3rd pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(3)
        write(*, '(14x, i2, 2x, ''Centering: 4th pass '',
     1    ''(1=2T, 2=O, 3=P, 4=C)'')') iord(4)
        write(*, '(5X, F11.4, 2X, ''Tolerance on 2 Theta'')') tol(1)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Omega'')') tol(2)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Phi'')') tol(3)
        write(*, '(5X, F11.4, 2X, ''Tolerance on Chi'')') tol(4)
        write(*, '(5X, F11.4, 2X, ''2 Theta scan range'')') del(1)
        write(*, '(5X, F11.4, 2X, ''Omega   scan range'')') del(2)
        write(*, '(5X, F11.4, 2X, ''Phi     scan range'')') del(3)
        write(*, '(5X, F11.4, 2X, ''Chi     scan range'')') del(4)
        write(*, '(5X, F11.4, 2X, ''Centering speed'')') cscan
        write(*, '(5X, F11.4, 2X, ''Phi search speed'')') phis
        write(*, '(10X, I6, 2X, ''Min intensity for phi search'')')
     1    iscnt
        write(*, '(10X, I6, 2X, ''Background intensity'')') scnhgb
        go to (311, 312, 313, 314) ircnt+1
  311     cdum = '           None'
        go to 310
  312     cdum = '       Recenter'
        go to 310
  313     cdum = '  Least squares'
        go to 310
  314     cdum = 'Data Collection'
  310   write(*, '(1X, A15, 2X, ''First operation'')') cdum
        write(*, '(10X, I6, 2X, ''Interval between recenter'')') nrcnt
        write(*, '(14X, I2, 2X, ''Max nb of reflections'')') maxfnd
        write(*, '(5X, F11.6, 2X, ''Search time [sec/step]'')') fstt
        write(*, '(5X, F11.4, 2X, ''Refl significance fact'')') fsig
        write(*, '(5X, F11.2, 2X, ''2 Theta min'')') flim1(1)
        write(*, '(5X, F11.2, 2X, ''Omega   min'')') flim1(2)
        write(*, '(5X, F11.2, 2X, ''Chi     min'')') flim1(3)
        write(*, '(5X, F11.2, 2X, ''2 Theta max'')') flim2(1)
        write(*, '(5X, F11.2, 2X, ''Omega   max'')') flim2(2)
        write(*, '(5X, F11.2, 2X, ''Chi     max'')') flim2(3)
        write(*, '(5X, F11.4, 2X, ''2 Theta max speed'')') xfrat(1)
        write(*, '(5X, F11.4, 2X, ''Omega   max speed'')') xfrat(2)
        write(*, '(5X, F11.4, 2X, ''Chi     max speed'')') xfrat(3)
C---- Orientation matrix calculation
        write(*, '(/'' Orientation matrix calculation'')')
        write(*, '(12X, I4, 2X, ''Number to reflections to be used'')')
     1    ifit
        write(*, '(12X, I4, 2X, ''Max value allowed for indexing'')')
     1    indx
        write(*, '(5X, F11.4, 2X, ''Min axis a value'')') axlth(1)
        write(*, '(5X, F11.4, 2X, ''Max axis a value'')') axlth(2)
        write(*, '(5X, F11.4, 2X, ''Min axis b value'')') axlth(3)
        write(*, '(5X, F11.4, 2X, ''Max axis b value'')') axlth(4)
        write(*, '(5X, F11.4, 2X, ''Min axis c value'')') axlth(5)
        write(*, '(5X, F11.4, 2X, ''Max axis c value'')') axlth(6)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis a value'')')
     1    dev(1)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis a value'')')
     1    dev(2)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis b value'')')
     1    dev(3)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis b value'')')
     1    dev(4)
        write(*, '(5X, F11.4, 2X, ''Max deviation min axis c value'')')
     1    dev(5)
        write(*, '(5X, F11.4, 2X, ''Max deviation max axis c value'')')
     1    dev(6)
        write(*, '(5X, F11.4, 2X, ''Max index deviation for integer'')')
     1    fm3
        write(*, '(12X, I4, 2X, ''Number of angles'')') izm
        write(*, '(12X, I4, 2X, ''Number solutions'')') msoln
        write(*, '(12X, I4, 2X, ''Number of input solutions'')') nsoln
        write(*, '(14X, I2, 2X, ''Offset for centering reflections'')')
     1    jtemp
        write(*, '(/''  c    h    k    l   wl     2 theta'', 7x,
     1    ''Omega         Phi         Chi           I'')')
        do 350 i = 1, nrf
          write(*, '(5(I3, 2X), 5(F10.4, 2X))')
     1      (crefl(j, i), j = 1, 5), (refl(j, i), j = 1, 5)
  350   continue
C---- Control parameters
        write(*, '(/'' Control parameters'')')
        go to (321, 322, 323) ismpl+1
  321     cdum = '            Off'
        go to 320
  322     cdum = '            Log'
        go to 320
  323     cdum = '          Print'
  320   write(*, '(1X, A15, 2X, ''Log file'')') cdum
        write(*, '(14X, I2, 2X, ''Number of steps to be printed'')')
     1    iprfl
        write(*, '(5X, F11.6, 2X, ''Plot scale'')') sclplt
        go to (331, 332, 333) iprof+1
  331     cdum = '           None'
        go to 330
  332     cdum = '           Line'
        go to 330
  333     cdum = '         Screen'
  330   write(*, '(1X, A15, 2X, ''Profile flag '')') cdum
        if (ismpl .eq. 1) then
          write(*, '(10X, I6, 2X, ''Number of pages in log file '')')
     1      ipage
        endif
C---- Check profiles
        write(*, '(/'' Plot one reflection parameters'')')
        go to (361, 362, 363, 364), itype+1
  361     cdum = '           None'
        go to 360
  362     cdum = '   Omega:2Theta'
        go to 360
  363     cdum = '         2Theta'
        go to 360
  364     cdum = '          Omega'
        go to 360
C 365     cdum = '            Phi'
C       go to 360
C 366     cdum = '            Chi'
  360   write(*, '(1X, A15, 2X, ''Scan type'')') cdum
        go to (341, 342, 343, 344) iprint+1
  341     cdum = '      Peak Plot'
        go to 340
  342     cdum = '         Counts'
        go to 340
  343     cdum = '      Ratemeter'
        go to 340
  344     cdum = '    Powers of 2'
  340   write(*, '(1X, A15, 2X, ''Screen output'')') cdum
        write(*, '(5X, F11.4, 2X, ''Scan range'')') pscan
        write(*, '(5X, F11.6, 2X, ''Sec/step'')') pltim
        write(*, '(5X, F11.4, 2X, ''Drive speed'')') ratem
C---- Photo parameters
        write(*, '(/'' Rotation photo parameters'')')
        write(*, '(5X, F11.4, 2X, ''Exposure time in min'')') ptim
        write(*, '(5X, F11.4, 2X, ''Distance form crystal to film'')')
     1    dm/2.0
        write(*, '(5X, F11.2, 2X, ''2 Theta position'')') zm
        write(*, '(5X, F11.2, 2X, ''Chi     position'')') qm
        write(*, '(5X, F11.2, 2X, ''Phi start'')') xm
        write(*, '(5X, F11.2, 2X, ''Phi end'')') ym
        write(*, '(5X, F11.4, 2X, ''Phi rate'')') prat
        write(*, '(/'' Axial photo parameters'')')
        write(*, '(5X, F11.2, 2X, ''Omega start'')') omg1
        write(*, '(5X, F11.2, 2X, ''Omega end'')') omg2
        write(*, '(5X, F11.4, 2X, ''Omega rate'')') omspd
C        WRITE(formfl, '('' Iprmp      I4  1 '', I6, 74X,
C     $    ''?'')') Iprmp
C        WRITE(formfl, '('' Iqus       I4  1 '', I6, 74X,
C     $    ''?'')') Iqus
C        WRITE(formfl, '('' XrySrc     C   2 '', A2, 78X,
C     $    ''?'')') XrySrc
C        WRITE(formfl, '('' FSD        C   2 '', A2, 78X,
C     $    ''?'')') FSD
C        WRITE(formfl, '('' HP         I2  1 '', I6, 74X,
C     $    ''?'')') HP
C        WRITE(formfl, '('' HE         I2  1 '', I6, 74X,
C     $    ''?'')') HE
C        WRITE(formfl, '('' Fnl        R   1 '', F10.4 ,70X,
C     $    ''?'')') Fnl
C        WRITE(formfl, '('' Ftpe       R   1 '', F10.4, 70X,
C     $    ''?'')') Ftpe
C        WRITE(formfl, '('' LSType     I   2 '', I4, 76X,
C     $    ''?'')') LSType
C      WRITE(formfl, '('' Fnono      R   3 '', 3(F6.2, 2X), 56X,
C     $  ''?'')') Fnono
C      WRITE(formfl, '('' Ivect      I4  3 '', 3(I6, 2X), 56X,
C     $  ''?'')') ivect
C      WRITE(formfl, '('' Ianum      I4  1 '', I6, 74X,
C     $  ''?'')') Ianum
C      WRITE(formfl, '('' Itotl      I4  1 '', I6, 74X,
C     $  ''?'')') Itotl
C      WRITE(formfl, '('' Mtape      I4  1 '', I6, 74X,
C     $  ''?'')') Mtape
C      WRITE(formfl, '('' Maxp       I4  1 '', I6, 74X,
C     $  ''?'')') Maxp
C      WRITE(formfl, '('' P3d        R  10 '', 10(F6.2, 2x),
C     $  ''?'')') p3d
C      WRITE(formfl, '('' Resmx      R   1 '', F10.4, 70X,
C     $  ''?'')') Resmx
C      WRITE(formfl, '('' Refps      R   1 '', F10.4, 70X,
C     $  ''?'')') Refps
C      WRITE(formfl, '('' Mlts       I4  1 '', I6, 74X,
C     $  ''?'')') Mlts
C
C---- End of normal job
C
      endif
      do 10 i = 61, 66
        title(i:i) = ' '
   10 continue
      if (iwl .eq. 2) wavel = flam1
      if (iwl .eq. 3) wavel = flam2
      stlmin = sin(dar*stlmin/2.0)/wavel
      stlmax = sin(dar*stlmax/2.0)/wavel
      return
      end
      subroutine choos
C
C----- Set up fixed characteristics for NICOLET 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 a Nicolet/Syntex diffractometer
C
      difnam = 'S'
      device = 'Siemens Nicolet P21 4 circle diffracto  '
      softwr = 'Siemens Nicolet P3R3 V4.11              '
C
C---- Tell the programme the name of the address file
      addnam = 'address.lst                                       '
C
C---- Set extension for data file name
      datxtn = '_REF.DCF'
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(2) = .false.
      ask(3) = .false.
      ask(4) = .false.
C
      return
      end
      subroutine initlz
C
C**** Subroutine for the NICOLET P21 to convert diffractometer parameter
C     information into the standardized form used by these programmes
C     The UB will be converted to the International Tables Volume 4
C     definition.
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(i) = 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---- Set up the wavelength for output
C
        do 201 i = 1, navtab
          if (abs(wavtab(1,i)-wavel).le.0.005) then
            flam1  = wavtab(2,i)
            flam2  = wavtab(3,i)
            weigh1 = wavtab(4,i)
            weigh2 = wavtab(5,i)
          endif
 201    continue
        cwave = flam1
C
C---- Horizontal and vertical apertures
C     ouverture : 2.5 mm environ, distance 135 mm environ
C
      detvt = 2.5/135.0
      dethz = 2.5/135.0
C
      return
      end
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Siemens-Nicolet P21 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 Nicolet P21 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 Nicolet P21 is thus a DSHDF type II diffractometer and
C     hence W = 2y according to DSHDF equation (11). (W and 1 being 3x3
C     matrices).
C
C     The Nicolet P21 uses the Eulerian cradle geometry described in the 
C     following
C
C    T        T     -1
C Psi  = s * u  * UB   * Z(-phi) * X(-chi) * 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     The matrices represent active rotation in the positive sense, i.e.
C     X:   e1 -> e1,   e2 -> e3,   e3 -> -e2
C     Y:   e1 -> -e3,  e2 -> e2,   e3 -> e1
C     Z:   e1 -> e2,   e2 -> -e3,  e3 -> e3
C
      integer  i,j
      real r(3), s(3), sinq, cosq
      real w(9)
      data w/-1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0/
      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---- *x(-chi)
      sinq = sin(dar*chi)
      cosq = cos(dar*chi)
      s(1) =  r(1)
      s(2) =  r(2)*cosq + r(3)*sinq
      s(3) = -r(2)*sinq + r(3)*cosq
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 II)
      do 115 i = 1, 3
        s(i) = 0.0
        do 115 j = 1, 3
          s(i) = s(i) + w(3*(i-1)+j)*r(j)
  115 continue
C
C---- Crystal based azimuth
      psicry = rad*atan2(s(1),s(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
