C00000000111111111122222222223333333333444444444455555555556666666666777
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C==== PW1100 Module of DIFRAC by H.D. Flack, University of Geneva
C==== Version 5th February 1997
C
      subroutine refrd(itype)
C**** ==================================================================
C     Routine to interpret reflection data from a PW1100. Itype 
C     indicates 0 for an ordinary reflection and 1 for a reference 
C     reflection. All data is output by calling REFOUT
C**** ==================================================================
      integer itype,i
      real pmt
      logical oka
C
      include 'common.for'
C
      if (itype .eq. 0) go to 2005
      go to 3005
C---- Ordinary reflection line
 2005 continue
C
C---- In the DIFRAC module for the Philips PW1100 diffractometer, 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 PW1100 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
      peak  = values(ptcsum)
      backl = values(ptcb1)
      backr = values(ptcb2)
      nrefl = nrefl + 1
      ptstd = 0
C---- Calculation of transformation coefficients
C     Pmt - Peak-with-background measuring time
      pmt = values(ptnsc)*values(ptswd)/spe
C     Bgt - left and right background counting time.
C     For bmo = 0 the time is that in values(ptbgt).
C     For bmo = 1 its half the peak scan time
C     For bmo = 2 the time is half that in values(ptbgt).
      if (bmo. eq. 0) then
        bgt = values(ptbgt)
      else if (bmo .eq. 1) then
        bgt = 0.5 * pmt
      else if (bmo .eq. 2) then
        bgt = 0.5 * values(ptbgt)
      endif
C     time = time + 1.0
      time = time + pmt + bgt + bgt
      cnp = pmt/values(ptswd)
      cnb = pmt
      cnl = bgt
      cnr = bgt
C---- Scan width
      swidth = values(ptswd)
C---- Scan ratio
      if (abs(theta) .le. sct) then
C---- omega scans
        scnrat = 0.0
        mstnn  = 1
      else
C---- omega - 2-theta scans
        scnrat = 2.0
        mstnn  = 2
      endif
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
      sintl = sint/wavel
      stlmin = min(sintl, stlmin)
      stlmax = max(sintl, stlmax)
C---- Keep on 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
C---- Reference reflection line
 3005 continue
      do 3095 i = 1 , numref
        ih = href(i)
        ik = kref(i)
        il = lref(i)
        peak = refint(i)
        backl = 0.0
        backr = 0.0
        psicry = psirf(i)
C       time = time + 1.0
        time = time + cnbrf(i)
        nrefl = nrefl + 1
        ptfilt = fltref(i)
        ptstd = i
        cnp = cnprf(i)
        cnb = cnbrf(i)
        cnl = 0.0
        cnr = 0.0
C---- Scan width and scan ratio
        swidth = swdrf(i)
        scnrat = scnrf(i)
        if (abs(scnrat) .le. 0.01) mstnn = 1
        if (abs(scnrat-2.0) .le. 0.01) mstnn = 2
        avint(i) = avint(i) + peak
        sdint(i) = sdint(i) + peak*peak
        count(i) = count(i) + 1.0
C
        call refout
 3095 continue
      return
      end
C
      subroutine dtflrd(jj)
C
C==== Subroutine to read a PW11OO diffractometer output file
C
C
      integer i,j,k,jj, trim
C
      include 'common.for'
C
C---- Initialisations
C
      firstr = .false.
C
C---- Get the data file name and open it
C
      if (jj .eq. 0) then
        ptscl = 1
        open( unit = datafl, file = compid(1:trim(compid))//datxtn,
     1  status = 'OLD', access = 'SEQUENTIAL', form = 'FORMATTED'
     2  )
CVMS 2  ,readonly)
CUNX 2  )
        write(*,'(/'' Reading file: '',A)')
     1   compid(1:trim(compid))//datxtn
      else
        ptscl = jj
        i = jj/10
        j = mod(jj, 10)
        open( unit = datafl, 
     1  file = compid(1:trim(compid))//datxto(1:trim(datxto))
     2    //digit(i)//digit(j), status = 'OLD', access = 'SEQUENTIAL', 
     3    form = 'FORMATTED'
     4  )
CVMS 4  ,readonly)
CUNX 4  )
        write(*,'(/'' Reading file: '',A)') compid(1:trim(compid))
     1    //datxto(1:trim(datxto))//digit(i)//digit(j)
      endif
C****
C     Interpretation of the .DAT file
C****
 1000 continue
      line = ' '
      read(datafl, '(a80)', end = 2000) line
C---- Look for PAR parameters
      if (line(36:38) .eq. 'SCT') then
        parclp = 1
        read (line, '(7x,f7.4,7x,f5.2,1x,f5.2,7x,f7.2)')
     1   spe, swda, swdb, sct
      elseif (line(32:34) .eq. 'MXN') then
        parclp = parclp + 2
        read (line, '(35x,i3)') mxn
      elseif (line(14:16) .eq. 'BMO') then
        parclp = parclp + 4
        read (line, '(9x,i1,9x,i1)') smo, bmo
        if (bmo. eq. 0) read (line, '(20x,f6.0)') bgt
      elseif (line(4:6) .eq. 'REF') then
        parclp = parclp + 8
        read (line, '(9x,i1,1x,f5.0)') numref, sinter
        read(datafl, '(a80)', end = 2000) line
        do 1005 i = 1 ,numref
          read(datafl, '(a80)', end = 2000) line
          read(line, '(6x,3i4,f8.4,i4)')
     1    href(i), kref(i), lref(i), speref(i), fltref(i)
 1005   continue
      elseif (line(4:6) .eq. 'SUB') then
        parclp = parclp + 16
        do 1010 i = 1 ,3
          read(datafl, '(a80)', end = 2000) line
          read(datafl, '(a80)', end = 2000) line
          read(line, '(3f12.6)') (ub(i,k), k = 1,3)
 1010   continue
      elseif (line(18:20) .eq. 'DET') then
        parclp = parclp + 32
        read(line, '(6x,f8.4,6x,2f6.2)') wavel, dethz, detvt
        if (parclp .eq. 63) then
          firstr = .true.
        else
        endif
      elseif (parclp .eq. 63) then
        if (line(8:10) .eq. 'RFA') then
          read(line, '(18x,f8.0,12x,f8.0,12x,f8.0)')
     1     (refint(i), i = 1 , numref)
          if (firstr) call initlz
          call refrd(1)
        else
C---- Check only for the characters in digit
          ok = .true.
          jj = 0
          do 1020 i = 1, 80
          do 1015 k = 0, 11
            if (line(i:i) .eq. digit(k)) go to 1019
 1015     continue
C---- Make a concession for column 22 when there is a variable scan 
C     width
          if ((line(22:22) .eq. '.') .and. (swdb .ne. 0.0)) go to 1020
          ok = .false.
          go to 1020
 1019     if (k .eq. 11) jj = jj+1
 1020     continue
C---- Skip blank lines
          if (jj.eq.80) ok = .false.
          if (ok) then
            if (firstr) call initlz
            read (line, fmt) ih, ik, il, psi, ptfilt,
     1                       (values(i),i=1,lfnorm)
            call refrd(0)
          else
          endif
        endif
      else
      endif
      go to 1000
C
C---- Finish off this data file
C
 2000 continue
      close (unit = datafl, status = 'KEEP')
      return
      end
C
      subroutine prflrd
C
C===== Subroutine to read a PW1100 diffractometer parameter output file
C
C===== On the Philips PW1100 there is no separate parameter output file
C      produced directly by the diffractometer. However we edit the 
C      .LOG file to produce a file containing input lines for standalone
C      a LATCON.
C
C      The file has two lines to skip at the top and then lines:
C      h, k, l, omega positive, omega negative
C
      include 'common.for'
      integer nrec, i, trim
      real omegp, omegn
      character*80  zharin
      logical zresnt
C
C==== Read the <Compid>.tth file of a PW1100
C
C---- Inquire <Compid>.tth and open it if present
C
      inquire(file=compid(1:trim(compid))//'.tth',
     1 exist=zresnt)
      if(.not.zresnt) then
        write(*,'(''   File: '',A,''.tth not available'')')
     1   compid(1:trim(compid))
      else
        write(*,'(/'' Reading file: '',A,''.tth'')')
     1   compid(1:trim(compid))
        open( unit=parmfl,file=compid(1:trim(compid))//'.tth',
     1  status='OLD', access='SEQUENTIAL', form='FORMATTED',err = 70
     2  ) 
CVMS 2  ,readonly) 
CUNX 2  )
C
C---- Header record
C
        read(unit=parmfl,fmt=*)
C
C---- Default wavelength index
C     0 = already defined default, 1 = alpha 1, 2 = alpha 2,
C     anything else = mean alpha
C
        zharin = ' '
        read(unit=parmfl,fmt='(a80)',end=70) zharin
        zharin(78:80) = ' 0 '
        read(unit=zharin,fmt=*) i
        if      (i.eq.1) then
          cdefwv = 1
        else if (i.eq.2) then
          cdefwv = 2
        else if (i.ne.0) then
          cdefwv = 0
        endif
C
C---- Treat  Reflection records
C
        do 2000 nrec = 1, 100
          zharin = ' '
          read(unit=parmfl,fmt='(a80)',end=70) zharin
          zharin(78:80) = ' 0 '
          read(unit=zharin,fmt=*)
     1     hang(nrec), kang(nrec), lang(nrec), omegp, omegn, wang(nrec)
C---- Pick up the two theta angles for "Good" angle reflections
          numang = numang + 1
          twoang(nrec) = abs(omegp) + abs(omegn)
          if (wang(nrec).eq.0) then
            wang(nrec) = cdefwv
          else if (wang(nrec).gt.2) then
            wang(nrec) = 0
          endif
          ctthn = amin1(twoang(nrec),ctthn)
          ctthx = amax1(twoang(nrec),ctthx)
 2000   continue
 70     close(unit=parmfl,status='KEEP')
        nmcell = numang
      endif
C
      return
      end
C
      subroutine choos
C
C----- Set up fixed characteristics for PW1100 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 Philips PW1100 diffractometer
      difnam = 'W'
      device = 'Philips PW1100 4-circle diffractometer  '
      softwr = 'Philips PW1100 April 1974               '
      cellpg = '?                                       '
      dettyp = 'scintillation counter                   '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = 'pck$:[cristal]address.lst                         '
CGENEVA addnam = '/home1/bin/address.lst                            '
CGENEVA addnam = '/users/crita/flack/bin/address.lst                '
CNORMAL addnam = 'address.lst                                       '
        addnam = '/home1/bin/address.lst                           '
C
C---- Set extension for data file name
      datxtn = '.dat'
      datxto = '.d'
C
C---- Set dihedral angle
      fmnhed = 90.0
C
C---- Default wavelength pointer for cell dimension determination, 
C     mean lambda
      cdefwv = 0
C
C---- Options for Quesan
C
      return
      end
C
      subroutine initlz
C
C**** Subroutine to interpret a complete PAR instruction in memory
C     for the PW1100.
C****
      integer i,j,k,iz
      real ubdet, temp2, temp3
      logical oka
C
      include 'common.for'
C
C---- Check PAR values
C
      if (spe .le. 0.0)  goto 999
      if (spe .gt. 0.32) goto 999
      if (swda .le. 0.0) goto 999
      if ((smo .lt. 1) .or. (smo .gt. 2)) goto 999
      if ((bmo. lt. 0) .or. (bmo .gt. 2)) goto 999
      if (bgt .lt. 0.0) goto 999
C
C---- Set up the format of the data collection lines and count number of
C---- items not including h,k,l,psi and F
C
      fmt = ' '
      fmt = '(3i4,F4.0,i2'
      lfnorm = 0
C---- Variable scan width?
      if (swdb .ne. 0.0) then
        fmt = fmt(1:index(fmt,' ')-1)//',f6.2'
        lfnorm = lfnorm + 1
        ptswd  = lfnorm
      else
        ptswd  = 7
        values(ptswd) = swda
      endif
C---- Variable number of scans?
      if (mxn .ne. 1) then
        fmt = fmt(1:index(fmt,' ')-1)//',f4.0'
        lfnorm = lfnorm + 1
        ptnsc  = lfnorm
      else
        ptnsc  = 8
        values(ptnsc) = 1.0
      endif
C---- Background counting time on file?
      if (bmo .eq. 2) then
        fmt = fmt(1:index(fmt,' ')-1)//',f6.0'
        lfnorm = lfnorm + 1
        ptbgt  = lfnorm
      elseif (bmo .eq. 1) then
        ptbgt = 9
C---- You don't know the length of the background counting time for 
C     bmo 1 at this point in time.
      elseif (bmo .eq. 0) then
        ptbgt = 9
        values(ptbgt) = bgt
      endif
C---- Background 1
      lfnorm = lfnorm + 1
      ptcb1  = lfnorm
C---- Peak with background
      lfnorm = lfnorm + 1
      ptcsum = lfnorm
C---- Background 2
      lfnorm = lfnorm + 1
      ptcb2  = lfnorm
C---- Finish off the format statement
      fmt = fmt(1:index(fmt,' ')-1)//',3f8.0)'
C
C---- Set up the wavelength for output
C
      flam1 = wavel
      flam2 = 0.0
      flam3 = 0.0
      weigh1 = 1.0
      weigh2 = 0.0
      weigh3 = 0.0
      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 up scan type for output: mstnn points into meas
C
      mstnn = 4
      if (sct .lt.  1.0) mstnn = 2
      if (sct. gt. 89.0) mstnn = 1
C---- Time interval in seconds between reference reflections
      sinter = 60.*sinter
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 
C     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
C---- Scale cell dimensions for wavelength
      do 403 i = 1, 3
       cellub(i) = wavel*cellub(i)
 403  continue
      volub=abs(wavel*wavel*wavel/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.
      temp3 = (cellub(1)+cellub(2)+cellub(3))/3.0
      if((abs(cellub(1)-temp3).lt.0.002*temp3).and.
     1   (abs(cellub(2)-temp3).lt.0.002*temp3).and.
     2   (abs(cellub(3)-temp3).lt.0.002*temp3)) then
        cellub(1) = temp3
        cellub(2) = cellub(1)
        cellub(3) = cellub(1)
      else
        temp2 = 0.5*(cellub(1)+cellub(2))
        if(abs(cellub(1)-cellub(2)).lt.0.002*temp2) then
          cellub(1) = temp2
          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 h k ls*****'')')
      endif
C
C---- Set up fixed values of reference reflections
C
      do 1005 iz = 1 , numref
C---- Get the crystal-based azimuthal angle
        ih = href(iz)
        ik = kref(iz)
        il = lref(iz)
        psi = 0.0
        call psic
        psirf(iz)=psicry
        if(trans) then
C---- Transform the stored reference hkl's
          href(iz) = -href(iz)
          kref(iz) = -kref(iz)
          lref(iz) = -lref(iz)
        endif
C---- Get sin theta needed for scan transition angle and swidth
        if (abs(theta) .le. sct) then
C---- omega scans
          scnrf(iz) = 0.0
        else
C---- omega - 2-theta scans
          scnrf(iz) = 2.0
        endif
        swdrf(iz) = 0.02*real(int(
     1   50.0*(swda + swdb* sint / sqrt(1.0 - sint*sint))))
        cnprf(iz) = 1.0/speref(iz)
        cnbrf(iz) = swdrf(iz)/speref(iz)
C---- Keep on eye on the filter codes that have been used
        oka = .false.
        do 2010 i = 1, nfilt
          if (fltref(iz) .eq. nfltcd(i)) oka = .true.
 2010   continue
        if (.not. oka) then
          nfilt = nfilt + 1
          nfltcd(nfilt) = fltref(iz)
        endif
 1005 continue
C
C---- Set Firstr to indicate that INITLZ does not need calling again.
      firstr = .false.
C
      return
 999  write(*,'('' *****Illegal par values read from diffractometer'',
     1 '' file. Interpretation stopped.***** '')')
      stop
      end
C
      subroutine psic
C
C**** Routine to calculate the crystal-based azimuth of the reflection
C     for the Philips PW1100 machine.
C
C     The equations used for the calculation of the crystal-based 
C     azimuth are given in "On the definition and Practical Use of 
C     Crystal-based Azimuthal Angles" by D.Schwarzenbach and 
C     H.D.Flack (1989) J. Appl. Cryst. 22, 601-605; D.Schwarzenbach and
C     H.D.Flack (1992) J. Appl. Cryst. 25, 69. hereafter DSHDF.
C     The definition of the reference vector u is given in equation (3)
C     of  DSHDF and equation (9) gives the psi-crystal.
C
C     The PW1100 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
C     is available (at least in the set up in Geneva where the output to
C     the teletype is captured on the hard disk of an IBM-PC look alike
C     instead of the Philips paper tape or magnetic tape output. The 
C     PW1100 is a DSHDF type II diffractometer and hence according to 
C     DSHDF equation (11):
C
C                   ( -1   0   0)
C             W =   (  0   1   0)
C                   (  0   0  -1).
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)
 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
      sinchb = -x(3)/d
      coschb =     p/d
      if ( p .eq. 0.0)  then
C---- Set phi to 45 degrees if both X and Y are ZERO
        sinphb = 0.707107
        cosphb = 0.707107
      else
        sinphb = x(1)/p
        cosphb = x(2)/p
      endif
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)
C     positive 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 II 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
C==== omega = rad*ATAN2(sinom,cosom) + theta
C==== chi   = rad*ATAN2(sinchi,coschi)
C==== phi   = rad*ATAN2(sinphi,cosphi)
C==== WRITE(12,'(1x,3i5,6f10.2)')
C====1 ih,ik,il, omega, twoth, chi, phi, psi, psicry
C
      return
      end
