c***********************************************************************
c                                                                      *
c    Program       S Y S T E R        version  0.9       Aug 24, 2000  *
c                                                                      *
c***********************************************************************
c                                                                      *
c    J.M.M. Smits  and  R. de Gelder                                   *
c    Department of Inorganic Chemistry                                 *
c    University of Nijmegen                                            *
c    Toernooiveld 1                                                    *
c    6525 ED  Nijmegen                                                 *
c    The Netherlands                                                   *
c    smits@sci.kun.nl   and   rdg@sci.kun.nl                           *
c                                                                      *
c    This program is freeware, which means that the user is invited    *
c    to adapt this program to his/her situation. However, some         *
c    subroutines are independent of any such situation, and should     *
c    be used unaltered. Proper credits would be appreciated, though.   *
c                                                                      *
c***********************************************************************
c
c
      program syster
      include "syster.inc"
      character*64    argument
      integer*4       argm_length
c
      isize = 20000
      call callcode ( argument, argm_length )
      call files ( argument, argm_length )
      call get_scale
      call get_isyst
      call get_fcf
      write (11,10)
   10 format ('    13')
   20 read (1,30,end=80) ihs, iks, ils, fos, sfos
   30 format (3i4,2f8.2)
      if ( ihs .eq. 0 .and. iks .eq. 0 .and. ils .eq. 0 ) goto 80
      read (2,40,end=80) nrefl, ihs2, iks2, ils2, ixrt, ibgl, ibgr,
     *                   theta, phik, omk, fkap
   40 format (i5,3i3,17x,i7,8x,2i6,f8.3,3f9.3)
      read (4,50,end=80) ihs2, iks2, ils2, drift
      read (5,50,end=80) ihs2, iks2, ils2, abscor
   50 format (3i4,f10.5)
      bckr = float( ibgl ) / float( ibgr )
      if ( bckr .lt. 1.0 ) bckr = 1.0 / bckr
      if ( ihs .ne. ihs2 .or. iks .ne. iks2 .or. ils .ne. ils2 ) then
         write (6,60) i, ihs, iks, ils, ihs2, iks2, ils2
   60    format ('fatal error, input data out of sync',/,
     *           'nrefl: ',i6,4x,3i5,'  --  ',3i5)
         stop
      end if
      ihn = ihs
      ikn = iks
      iln = ils
      call norm_hkl ( ihn, ikn, iln )
      jhkl = (ihn+500)*1000000 + (ikn+500)*1000 + (iln+500)
      call find_fc ( ihkl, nfcf, jhkl, index )
c
c     if index = 0: jhkl not found in ihkl, i.e. systematically absent
c
      if ( index .eq. 0 ) goto 20
      pfac = ( max( fos, 0.0 ) + 2.0 * fc( index ) * fc( index ) ) / 3.0
      weight = sfos + wfac1 * wfac1 + pfac * pfac + wfac2 * pfac
      weight = 1.0 / weight
      if ( fos .ge. 0.0 ) then
         fos  = sqrt( fos )
         sfos = 0.5 * sfos / fos
      else
         fos  = -sqrt( -fos )
         sfos = -0.5 * sfos / fos
      end if
c
c     reset smallest values of Fo and Fc to 0.05 and sig(Fo) = 0.15
c     i.e. if -0.05 <= Fo**2 => 0.05. Keep the sign of Fo.
c
      if ( fos .ge. 0.0 .and. fos .le. 0.05 )  then
         fos  = 0.0500
         sfos = 0.1500
      end if
      if ( fos .lt. 0.0 .and. fos .ge. -0.05 )  then
         fos  = -0.0500
         sfos =  0.1500
      end if
      fcs = max( fc( index ), 0.05000 )
      write (11,70) float(nrefl), fos, fcs, sfos, float(ixrt),
     *              theta, phik, omk, fkap, bckr,
     *              abscor, drift, weight
   70 format (f10.0,3f10.4,f10.0,4f10.3,4f10.5)
      goto 20
c
c
   80 write (6,90)
   90 format ('syster finished')
      stop
      end
c
c
      subroutine get_scale
      include "syster.inc"
      logical*4  found_scale, found_weight
c
c     subroutine get_scale gets the scale factor to scale fc to fo
c     from the SHELXL *.res file
c
      scale = 1.0
      wfac1 = 1.0
      wfac2 = 0.0
      found_scale  = .false.
      found_weight = .false.
   10 read (8,20,end=50) line
   20 format (a)
      if ( line(1:6) .eq. 'FVAR  ' ) then
         read (line(7:18),30) scale
   30    format (f12.5)
         found_scale = .true.
      else if ( line(1:6) .eq. 'WGHT  ' ) then
         read (line(7:28),40) wfac1, wfac2
   40    format (f10.6,f12.6)
         found_weight = .true.
      end if
      if ( found_scale .and. found_weight ) return
      goto 10
c
   50 return
      end
c
c
      subroutine get_isyst
      include "syster.inc"
c
c     subroutine get_isyst gets several data from the *.crysda file
c
      icent = 0
      isyst = 8
   10 read (9,20,end=40) line
   20 format (a)
      if      ( line(1:6) .eq. 'ICENT ' ) then
         read (line(19:20),30) icent
   30    format (i2)
      else if ( line(1:6) .eq. 'ISYST ' ) then
         read (line(19:20),30) isyst
      end if
      goto 10
c
   40 return
      end
c
c
      subroutine get_fcf
      include "syster.inc"
c
c     subroutine get_fcf reads the calculated data from the *.fcf
c     file, packs the hkl-values into one word, and sorts the data
c     on these packed hkl-values and the calculated data accordingly.
c
      nfcf = 0
   10 read (3,20,end=60) line
   20 format (a)
      if ( line(1:23) .ne. ' _refln_observed_status' ) goto 10
   30 read (3,40,end=60) ih, ik, il, fc2
   40 format (3i4,f12.2)
      if ( ih .eq. 0 .and. ik .eq. 0 .and. il .eq. 0 ) goto 30
      nfcf = nfcf + 1
      if ( nfcf .gt. isize ) then
         write (6,50) isize
   50    format ('error stop: more than ',i5,' calculated data',/,
     *           'resize arrays, reset parameter isize')
         stop
      end if
      itmp(nfcf) = (ih+500)*1000000 + (ik+500)*1000 + (il+500)
      ftmp(nfcf) = scale * sqrt( fc2 )
      goto 30
c
c Sort ihkl in ascending order and fc accordingly
c
   60 call sorti ( itmp, nfcf, ip )
      do 70 i = 1, nfcf
         ihkl(nfcf+1-i) = itmp( ip(i) )
         fc(nfcf+1-i)   = ftmp( ip(i) )
   70 continue
      return
      end
c
c
      subroutine norm_hkl ( ihs, iks, ils )
      include "syster.inc"
c
c     subroutine norm_hkl converts the hkl-values ihs, iks and ils
c     to the 'normalized' values as used by SHELXL
c
      if ( isyst .eq. 1 ) then
         if ( icent .eq. 1 ) return
         if ( ils .lt. 0 ) then
            ihs = -ihs
            iks = -iks
            ils = -ils
         end if
         if ( ils .eq. 0 .and. iks .lt. 0 ) then
            ihs = -ihs
            iks = -iks
         end if
         if ( ils .eq. 0 .and. iks .eq. 0 .and. ihs .lt. 0 ) then
            ihs = -ihs
         end if
      else if ( isyst .eq. 2 ) then
         if ( ils .lt. 0 ) then
            ihs = -ihs
            iks = -iks
            ils = -ils
         end if
         if ( iks .lt. 0 ) then
            iks = -iks
         end if
         if ( ihs .lt. 0 .and. ils .eq. 0 ) then
            ihs = -ihs
         end if
      else if ( isyst .eq. 3 ) then
         ihs = iabs( ihs )
         iks = iabs( iks )
         ils = iabs( ils )
      end if
c
      return
      end
c
c
      subroutine find_fc ( ihkl, nfcf, jhkl, index )
      dimension ihkl(nfcf)
c
c     subroutine find_fc finds the position of jhkl in an ordered list
c     of integer values ihkl(1..nfcf), using a binary search strategy.
c     on return, index is set to the position of jhkl in array ihkl;
c     index = 0 if jkhl is not found in array ihkl
c
      nlow = 1
      nhigh = nfcf
   10 if ( nhigh - nlow .le. 1 ) then
         if ( ihkl(nlow) .eq. jhkl ) then
            index = nlow
            goto 20
         else if ( ihkl(nhigh) .eq. jhkl ) then
            index = nhigh
            goto 20
         else
            index = 0
            goto 20
         end if
      end if
      index = ( nhigh + nlow ) / 2
      if ( ihkl(index) .gt. jhkl ) then
         nhigh = index
         goto 10
      else if ( ihkl(index) .lt. jhkl ) then
         nlow = index
         goto 10
      end if
c
   20 return
      end
c
c
      subroutine sorti ( ix, nx, ip )
      dimension ix(nx), ip(nx)
c
c  sort ix(nx) creating pointers array ip
c  on completion ix(ip(1)) will be the largest value of ix()
c
      do 10 i = 1, nx
         ip(i) = i
   10 continue
      int = 2
   20 int = int + int
      if ( int .lt. nx ) goto 20
      int = min0( nx, (3*int)/4-1 )
   30 int = int/2
      ifin = nx-int
      do 70 ii = 1, ifin
         i = ii
         j = i+int
         if( ix(ip(i)) .ge. ix(ip(j)) ) goto 70
         k = ip(j)
   40    ip(j) = ip(i)
         j = i
         i = i-int
         if ( i ) 60, 60, 50
   50    if( ix(ip(i) ) .lt. ix(k) ) goto 40
   60    ip(j) = k
   70 continue
      if ( int .gt. 1 ) goto 30
      return
      end
c
c
      subroutine callcode ( argument, argm_length )
      character*(*)   argument
      integer*4       argm_length
c
c     subroutine callcode gets the first argument given in the program
c     calling sequence; if no argument is given it is asked for.
c     it is assumed that this argument is the compound code, which is
c     subsequently used to construct several file names.
c
c     this subroutine uses two special statements which are / may be
c     platform dependent:
c     iargc()      integer function, returns the number of arguments
c                  given in the calling sequence, excluding the program
c                  name itself; returns 0 if no arguments are given.
c     getarg(i,s)  this subroutine gets the i-th argument from the
c                  calling sequence, excluding the program name, and
c                  stores the argument value in character string s.
c
      numarg = iargc()
      if ( numarg .eq. 0 ) then
         write (6,10)
   10    format ('Enter the compound code: ',$)
         read (5,20) argument
   20    format (a)
      else
         call getarg (1, argument)
      end if
      argm_length = leng( argument )
c
      return
      end
c
c
      integer function leng ( string )
      character*(*) string
c
c     integer function leng returns the length in bytes of the
c     meaningful part of a string, i.e. excluding trailing blanks;
c     leng = 0 if an empty string is supplied.
c     to be used if no such intrinsic is supplied by the compiler.
c     n.b. the standard intrinsic integer function len returns the
c     defined length of a string in bytes.
c
      ls = len( string )
      do 10 i = ls, 1, -1
         if ( string(i:i) .ne. ' ' ) goto 20
   10 continue
      leng = 0
      return
c
   20 leng = i
      return
      end
c
c
      subroutine oldfile ( iunit, filename )
      character*(*) filename
      logical*1 exists
c
c     subroutine oldfile tests the presence of file 'filename', assumed
c     to be present, and subsequently opens it on unit iunit.
c
      lenf = leng( filename )
      inquire ( file = filename(1:lenf), exist = exists,
     *          err = 20, iostat = ios )
      if ( exists ) then
         open ( file = filename(1:lenf), unit = iunit, status = 'old',
     *          err = 40, iostat = ios )
      return
      else
         write (6,10) filename(1:lenf)
   10    format ('Fatal error, file ',a,' does not exist')
      stop
      end if
c
   20 write (6,30) filename(1:lenf), ios
   30 format ('Fatal error, processor error during inquire on file ',
     * a,'; ios = ',i3)
      stop
c
   40 write (6,50) filename(1:lenf), ios
   50 format ('Fatal error, processor error during opening of file ',
     * a,'; ios = ',i3)
      stop
      end
c
c
      subroutine newfile ( iunit, filename )
      character*(*) filename
c
c     subroutine newfile opens a new file 'filename' on unit iunit;
c     status 'unknown' is used, i.e. an axisting file is overwritten.
c
      lenf = leng( filename )
      open ( file = filename(1:lenf), unit = iunit, status = 'unknown',
     *       err = 10, iostat = ios )
      return
c
   10 write (6,20) filename(1:lenf), ios
   20 format ('Fatal error, processor error during opening of file ',
     * a,'; ios = ',i3)
      stop
      end
c
c
      subroutine files ( argument, argm_length )
      character*(*)   argument
      integer*4       argm_length
      character*64    filename
c
c     subroutine files opens all the necessary input and output
c     files, using subroutines oldfile and newfile as appropriate.
c
c     unit   extension
c                        input files containing reflection data
c                        -----------------------------------------------
c       1    *.hkl       data reduction output, source of Fo and sig(Fo)
c       2    *.rf1       measurement output, source of angular info
c       3    *.fcf       SHELXL output, source of Fc
c       4    *.drift     data reduction output, source of drift info
c       5    *.empabs    data reduction output, source of abscor info
c                        -----------------------------------------------
c                        input files containing crystal data
c                        -----------------------------------------------
c       8    *.res       SHELXL output, source of scale and weight
c       9    *.crysda    local program output, source of symmetry info
c                        (see manual for specifics and how to get it)
c                        -----------------------------------------------
c                        output file
c                        -----------------------------------------------
c      11    *.syster    current program output, input to SYSTERPLOT
c
      filename = argument(1:argm_length) // '.hkl'
      call oldfile ( 1, filename )
      filename = '../data/' // argument(1:argm_length) // '.rf1'
      call oldfile ( 2, filename )
      filename = argument(1:argm_length) // '.fcf'
      call oldfile ( 3, filename )
      filename = '../data/' // argument(1:argm_length) // '.drift'
      call oldfile ( 4, filename )
      filename = '../data/' // argument(1:argm_length) // '.empabs'
      call oldfile ( 5, filename )
      filename = argument(1:argm_length) // '.res'
      call oldfile ( 8, filename )
      filename = argument(1:argm_length) // '.crysda'
      call oldfile ( 9, filename )
      filename = argument(1:argm_length) // '.syster'
      call newfile ( 11, filename )
      return
      end
