c ########################################################################
c  Copyright (c) 1991  Siemens Analytical X-ray Instruments (SAXI)
c                      6300 Enterprise Lane
c                      Madison, WI 53719 USA          All rights reserved.
c ########################################################################
c
c RAW_LIB.FOR - Primitives to read or write powder raw data traces in PLOTSO
c   file format as used by various SAXI programs (GADDS, SQUID, etc.)
c
c Version 1.0 F77L3                   18Oct91 -- KLS
c
*+
c raw_plotso_read -  reads a single range of a PLOTSO style file
c
	integer function raw_plotso_read ( rawfile, nrange, raw, angbeg,
     &     angend, anginc, wave, title)
	implicit none
	character*(*) rawfile   !(R) Output file name
	integer nrange          !(R) Range number
	real raw(*)             !(W) Intensity array
	real angbeg             !(W) Starting angle
	real angend             !(W) Ending angle
	real anginc             !(W) Step width (angle increment)
	real wave(3)            !(W) Wavelengths (avg, a1, a2)
	character*(*) title     !(W) Title
*-
	character line*80
	integer ier, jer, iopen, npts, lens, ncurve, lun
	real x, y
c
c Input file must exist
	lun = iopen (rawfile, ' ', 'OLD', 'SEQUENTIAL',
     &     'FORMATTED', 0, 0, ier)
	if (ier.ne.0) goto 90
c
c Read PLOTSO file
      npts = 0
      ncurve = 1
      do while (.true.)
	read (lun, '(A)', end=80) line
	if (lens(line).le.1) then
c         ! null
	else if (line(1:3).eq.'!@!') then
	  if (line(1:10).eq.'!@!!Title') then
	    title = line(11:)
	  else if (line(1:16).eq.'!@!!Wavelengths') then
	    read (line(17:),'(3F10.5)') wave
	  else if (line(1:16).eq.'!@!A') then
	    title = line(5:)
	  else if (line(1:4).eq.'!@!N' .or. line(1:4).eq.'!@!P') then
	    if (npts.gt.0) then
	      ncurve = ncurve + 1
	      if (ncurve.gt.nrange) goto 80
	      npts = 0
	    endif
	  end if
	else
	  npts = npts + 1
	  if (nrange.eq.ncurve) then
	     read (line(1:lens(line)),*) x, y
	     raw(npts) = y
	     if (npts.eq.1) angbeg = x
	     angend = x
	  end if
	end if
      end do
c
   80   if (ncurve.lt.nrange) then
	  call prg_set_error_string('Requested range not found')
	  ier = -2
	  goto 90
	end if
c
	anginc = 0.02
	if (npts.gt.1) anginc = (angend - angbeg) / (npts - 1)
c
   90   call iclose (lun, jer)
	if (ier.eq.0) ier = jer
	raw_plotso_read = ier
	return
	end
*+
c RAW_PLOTSO_WRITE - writes a single range to a PLOTSO style file
c
	integer function raw_plotso_write ( rawfile, append, raw,
     &     angbeg, angend, anginc, wave, title, mode)
	implicit none
	character*(*) rawfile   !(R) Output file name
	integer append          !(R)non-zero to append data to existing file
	real raw(*)             !(R) Intensity array
	real angbeg             !(R) Starting angle
	real angend             !(R) Ending angle
	real anginc             !(R) Step width (angle increment)
	real wave(3)            !(R) Wavelengths (avg, a1, a2)
	character*(*) title     !(R) Title
	character*(*) mode      !(R) Mode: 'Chi ', '2th ', 'Slic', or '    '
*-
	character string*24, systyp$p
	integer ier, lun, iopen, i, npts, lens, jer, iprecision
	real angle, quantize
c
	ier = 0
	if (rawfile(1:1).eq.' ') goto 9999
c
	lun = iopen (rawfile, ' ', 'UNKNOWN', 'SEQUENTIAL',
     &     'FORMATTED', 0, 0, ier)
	if (ier.ne.0) goto 90
c
c Append file: read to end of file
	if (append.ne.0) then
    3      read (lun,'(A)', end=5) string
	   goto 3
    5      if (systyp$p().eq.'N') backspace lun    ! for SVS-FORTRAN !
	end if
c
	if (mode.ne.'    ') write (lun,10) mode(1:lens(mode))
   10   format('!@!!GADDS-I-Created by ',A,' integration of area ',
     &     'detector frame')
	write (lun,11) title(1:max(1,lens(title)))
   11   format('!@!!Title ',A)
	write (lun,12) wave
   12   format('!@!!Wavelengths ',3F10.5)
	write (lun,14)
   14   format('!@!N')
	write (lun,15)
   15   format('!@!SS')
	write (lun,16)
   16   format('!@!M')
	write (lun,17)
   17   format('!@!L 0.0 0.0 0.0 0.0')
	write (lun,18)
   18   format('!@!XDegrees')
	write (lun,20)
   20   format('!@!YCounts')
c
c Determine precision (number of places to right of decimal point) for angles
        iprecision = 4
        write (string, '(F10.4)', iostat=ier) anginc
        if (string(10:10).eq.'0') iprecision = 3
        if (string(9:10).eq.'00') iprecision = 2
        if (string(8:10).eq.'000') iprecision = 1
c
	npts = 1 + nint( (angend - angbeg) / anginc )
	angle = angbeg
	do 80  i = 1, npts
           if (iprecision.eq.4) then
    	      write (string,75,iostat=ier) angle, raw(i)
           else if (iprecision.eq.3) then
    	      write (string,76,iostat=ier) angle, raw(i)
           else if (iprecision.eq.2) then
    	      write (string,77,iostat=ier) angle, raw(i)
           else
    	      write (string,78,iostat=ier) angle, raw(i)
           end if
   75      format(F10.4,1x,F10.2)
   76      format(F10.3,1x,F10.2)
   77      format(F10.2,1x,F10.2)
   78      format(F10.1,1x,F10.2)
c          ! remove trailing zeros to right of decimal point from intensity
           if (string(20:21).eq.'00')  string(20:21) = '  '
           if (string(21:21).eq.'0')   string(21:21) = ' '
           if (string(18:20).eq.' . ') string(20:20) = '0'
	   call strcompress (string)
	   write (lun,'(A)') string(:lens(string))
	   angle = angle + anginc
	   angle = quantize( angle, anginc )
   80   continue
c
c Truncate and close or just close data file
	if (append.eq.0) endfile (unit=lun, iostat=ier)
   90   call iclose (lun, jer)
	if (ier.eq.0) ier = jer
 9999   raw_plotso_write = ier
	return
	end