C
C==== STADI 4 Module of DIFRAC by H.D. Flack, University of Geneva
C==== PC version of the control programme
C==== Version 8th November 1999
      subroutine refrd(itype)
C**** =======================================================================
C     Routine to interpret reflection data from a STADI 4.
C     All data is output by calling REFOUT
C**** =======================================================================
      integer itype
      real xh, xk, xl
      real facx, aa, bb, cc, dd, x1, y1, x2, y2
      real zint3, su3
C
      include 'common.for'
C
C---- In the DIFRAC module for the Stoe Stadi 4 diffractometer, the
C---- UB, UBINV and UBDET variables contain the values used on the
C---- diffractometer even if the determinant turns out to be negative.
C---- The diffractometer hkls are transformed for output.  
C---- The "new" output file of the Stadi-4 contains the crystal-based 
C---- azimuthal angle. For the id1.eq.1 reflections the stoe file contains
C---- the inverted h k l's but the correct crystal-based azimuthal angle.
C
C---- Friedel opposite measured at negative 2 theta
      if (id1.eq.1) then
        ih = -ih
        ik = -ik
        il = -il
        psicry = -psicry
      endif
C---- Transform output h k l for negative UB
      if (trans) then
        ih = -ih
        ik = -ik
        il = -il
      endif
C
C==== Calculate the Flack, Blanc, Schwarzenbach coefficients
C
      facx = 100.0*real(nstbg)*stepwd/(real(nstep))
      cnb = (1.0 - 2.0*facx)*stptim
      cnp = stptim/(0.01*real(nstep))
      cnl = facx*stptim
      cnr = cnl
C
C===== Invent the raw data
C
      backl = (cnp*cnl*cnl*(cnp*signet*signet - intnet))
     1  /(cnb*(cnb+cnl))
      backr = backl
      peak = (cnp*(cnp*cnl*signet*signet+cnb*intnet))/(cnb+cnl)
C
C===== Round the counts and refine the value of the c's
C
      backl = anint(backl)
      if (backl.eq.0.0) then
        if (backr.gt.0.0) then
          backl = 1.0
         else
          backl = -1.0
         endif
      endif
      backr = backl
      peak  = anint(peak)
      aa = (backl*signet*signet - intnet*intnet)*cnp*cnp
      if (aa.eq.0.0)
     1  write(*,'('' Zero aa - No reset programmed'', 3i5)') ih, ik, il
      bb = 2.0*peak*intnet*cnp
      cc = -peak*(peak+backl)
      dd = bb*bb - 4.0*aa*cc
      if (dd.le.0.0) write(*,'('' Negative dd'', 3i5, e13.5)')
     1 ih, ik, il, dd
      dd = sqrt(dd)
      x2 = (-bb+dd)/(2.0*aa)
      y2 = (peak - intnet*cnp*x2)*cnl/(backl*cnb)
      if ((x2.lt.0.0).OR.(y2.lt.0.0)) then
        x1 = (-bb-dd)/(2.0*aa)
        y1 = (peak - intnet*cnp*x1)*cnl/(backl*cnb)
        if ((x1.lt.0.0).OR.(y1.lt.0.0)) then
          cnb = y2*cnb
          cnp = x2*cnp
          write (*,'('' Abnormal scales '', 4f13.5)') x2, y2, x1, y1
        else
          cnb = y1*cnb
          cnp = x1*cnp
        endif
      else
        cnb = y2*cnb
        cnp = x2*cnp
      endif
      zint3 = (peak - cnb*backl/cnl)/cnp
      su3 = sqrt(peak + cnb*cnb*backl/(cnl*cnl))/cnp
      if (
     1 (abs(intnet-zint3).gt. 0.01) .OR. (abs(signet-su3).gt. 0.01)
     8 ) write (*,'('' Reconstitution error '', 4f13.5)')
     9 intnet, signet, zint3, su3
C
C---- Mark standards and keep indices
C---- Set standard flag
      ptstd = -999
      if (id1.ne.4) then
        ptstd = 0
        stdseq = 0
      else
        stdseq = stdseq + 1
        if(stdseq.gt.numref) stdseq=1
        if ((ih.eq.href(stdseq)).and.(ik.eq.kref(stdseq)).and.
     1      (il.eq.lref(stdseq))) ptstd = stdseq
      endif
      if(ptstd.eq.-999) write(*,
     1   '('' Standard out of order'',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
C---- Update other variables
C
      time = time + stptim
      nrefl = nrefl + 1
C---- Scan width
      swidth = 0.01*real(nstep)
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---- Calculate sin(theta)/lambda
      xh=real(ih)
      xk=real(ik)
      xl=real(il)
      sintl = 0.5*sqrt(xh*xh*gr(1,1)+xk*xk*gr(2,2)+xl*xl*gr(3,3)
     1 +2.0*(xh*xk*gr(1,2)+xh*xl*gr(1,3)+xk*xl*gr(2,3)))
C---- Keep min/max values of sintl
      stlmin = min(sintl, stlmin)
      stlmax = max(sintl, stlmax)
C---- Keep on eye on the filter codes that have been used
      nfilt = max(ptfilt+1, nfilt)
C
      call refout
      return
      end
C
      subroutine dtflrd(jj)
C
C==== Subroutine to read a STADI 4 diffractometer output file
C
C
      integer i,j,k,jj
      integer trim
C
      include 'common.for'
C0000000011111111112222222222333333333344444444445555555555666666666677777777778
C2345678901234567890123456789012345678901234567890123456789012345678901234567890
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, '(a81)', end = 2000) line
C---- Look for non-reflection parameters
      if (line(1:4) .eq. 'File') then
        parclp = 1
C---- Correct for absent Attn line
        parclp = parclp + 16
      elseif (line(1:4) .eq. 'Titl') then
        title = line(7:72)
        parclp = parclp + 2
      elseif (line(1:4) .eq. 'Wave') then
        parclp = parclp + 4
        read (line, '(4x,3f10.6)') wavel, flam1, flam2
      elseif (line(1:4) .eq. 'Mono') then
        parclp = parclp + 8
        read (line, '(4x,f7.2,f4.0,2f8.4)') fmn2th, fmnhed
      elseif (line(1:4) .eq. 'Attn') then
        parclp = parclp + 16
        read (line, '(4x,f7.3,6f8.3)') (attfac(i), i=2,8)
        read(datafl, '(a80)', end = 2000) line
        read (line, '(4x,f7.3,6f8.3)') (ettfac(i), i=2,8)
        do 995 i = 2, 8
 995    nfltcd(i) = i-1 
      elseif (line(1:4) .eq. 'DTim') then
        parclp = parclp + 32
        read (line, '(4x,2f6.2)') fdeadt, efdead
      elseif (line(1:4) .eq. 'Scan') then
        parclp = parclp + 64
        read (line, '(4x,2f6.3,f7.3,6x,i1,11x,i1,11x,i1)')
     1      szwidt, scnrat, stepwd, nstbg, prfout, profit
      elseif (line(1:4) .eq. 'Slit') then
        parclp = parclp + 128
        read (line, '(4x,f6.2)') detvt
        dethz = detvt
      elseif (line(1:4) .eq. 'Stds') then
        parclp = parclp + 256
        read (line, '(5x,i2,6x,10(2x,3i2))')
     1    numref, (href(i), kref(i), lref(i), i=1,numref)
      elseif (line(1:4) .eq. 'Omat') then
        parclp = parclp + 512
        do 1010 i = 1 ,3
          if(i.ne.1) read(datafl, '(a80)', end = 2000) line
          read(line, '(4x,3f12.8)') (ub(i,k), k = 1,3)
 1010   continue
      elseif (line(1:4) .eq. 'Data') then
        parclp = parclp + 1024
        if (parclp .eq. 2047) firstr = .true.
      elseif (parclp .eq. 2047) then
C---- Check only for the characters in digit
        ok = .true.
        do 1020 i = 1, 80
          do 1015 k = 0, 11
            if (line(i:i) .eq. digit(k)) go to 1020
 1015     continue
C---- Allow . as well
          if (line(i:i) .eq. '.') go to 1020
          ok = .false.
 1020   continue
        if (ok) then
          if (firstr) call initlz
          if (profit.eq.0) then
            read (line,
     1  '(3i4,i2,1x,f8.1,f6.1,4x,i5,f5.0,i2,f7.2)')
     2        ih,ik,il,id1,intnet,signet,nstep,stptim,ptfilt,
     3         psicry
C     1  '(3i4,i2,i1,f8.1,f6.1,4x,i5,f4.0,i2,f8.0,2f6.0,f8.2)')
C     2        ih,ik,il,id1,id2,intnet,signet,nstep,stptim,ptfilt,
C     3        peak, backl, backr, psicry
          else if(profit.eq.1) then
            read (line,
     1 '(3i4,i2,i1,f8.1,f6.1,4x,i4,f4.0,i2,f8.0,2f6.0,i3,f7.4,f8.2)')
     2        ih,ik,il,id1,id2,intnet,signet,nstep,stptim,ptfilt,
     3        peak, backl, backr, nstbg, stepwd, psicry
          endif
          call refrd(0)
        else
        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 STADI 4 diffractometer parameter output file
C
      return
      end
C
      subroutine choos
C
C----- Set up fixed characteristics for STADI 4 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 Stoe Stadi 4 diffractometer
      difnam = 'T'
      device = 'Stoe Stadi-4 4-circle diffractometer  '
      softwr = 'Stoe Stadi-4                          '
      cellpg = '?                                       '
      dettyp = 'scintillation counter                   '
C
C---- Tell the programme the name of the address file
CGENEVA addnam = '[flack.cristal]address.lst                        '
CNORMAL addnam = 'address.lst                                       '
        addnam = '/export/local/bin/address.lst                     '
C
C---- Set extension for data file name
      datxtn = '.d00'
      datxto = '.d'
C
C---- Default wavelength pointer for cell dimension determination, mean lambda
      cdefwv = 0
C---- Background scanning mode is moving
      backsn = 1
C
C---- Options for Quesan
C
      ask(3) = .false.
      ask(4) = .false.
      ask(6) = .false.
      ask(7) = .false.
      return
      end
C
      subroutine initlz
C
C**** Subroutine to interpret a complete initial instructions
C     for the STADI 4.
C****
      integer i,j,k
      real ubdet
      real temp2, temp3
C
      include 'common.for'
C
C---- Set up the wavelength for output
C
      flam3 = 0.0
      weigh1 = 0.6666667
      weigh2 = 0.3333333
      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---- Convert the scan type parameter
C
        mstnn = 3
        if(abs(scnrat).lt. 0.1)       mstnn = 1
        if(abs(scnrat - 2.0).lt. 0.1) mstnn = 2
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.
      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---- Get reciprocal metric in GR from ub
      do 405 i = 1,3
      do 405 k = 1,3
        gr(i,k) = 0.0
      do 405 j = 1,3
        gr(i,k) = gr(i,k) + ub(j,i)*ub(j,k)
  405 continue
C
C---- Set Firstr to indicate that INITLZ does not need calling again.
      firstr = .false.
C
      return
 999  write(*,'('' *****Illegal values read from diffractometer'',
     1 '' file. Interpretation stopped.***** '')')
      stop
      end
