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       PKS_LIB.FOR primitives for reading and writing to peak data type
c       file formats: SQUID files.
c
c       Version 1.0 F77L3                   08May91 -- KLS
*+
c -----------------------------------------------------------------------
c PKS_WRITE - Routine to write PEAKS information to file. Recalculates the
c   relative intensities. No data, leaves the field blank.
c
*+
	integer function pks_write (filename, append)
	implicit none
	character*(*) filename  !(R) peak file to create/append
	integer append          !(R) non-zero to append
*-
	include 'pks_cmn.inc'   !(R) all peak information
	include 'raw_var.inc'   !(R) wavelengths & title information
	character fname*34, string*80, status*8
	integer i, ier, j, jer, l, lun
	real imax
c
	integer iopen, lens
c
	status = 'UNKNOWN'
	lun = iopen (filename, ' ', status, 'SEQUENTIAL',
     &    'FORMATTED', 80, 0, ier)
	if (ier.ne.0) goto 90
c
c Append mode
	if (append.ne.0) then
	  do while (.true.)
	    read (lun,'(A)',err=5,end=5) string
	  end do
05        backspace lun         ! for SVS FORTRAN !
	end if
c
c Write header information
c Siemens Analytical X-Ray Instruments. Copyright 1991. All rights reserved.
c File: \gadds\test.pks                  Created: 12May1991 12:00:00.0
	write (lun, '(2A)') 'Bruker Analytical X-Ray Systems. Copy'
     &    ,'right 1997. All rights reserved.'
	l = lens(filename)
	if (l.lt.34) then
	  fname = filename
	else
	  fname = '...'//filename(l-30:l)
	end if
	call ascdate(-1,-1,-1,string)
	call fasctime(-1,-1,-1,string(11:))
	write (lun, '(5A)') 'File: ',fname,'Created: ',string(1:24)
	write (lun, '(2A)') 'Title: ',Raw_title(:max(1,lens(Raw_title)))
	write (lun, '(A,4(2X,A2,A,F8.5))') 'Wavelengths:',
     &    Raw_target, 'Kavg:', Raw_wavelengths(1),
     &    Raw_target, 'Ka1:', Raw_wavelengths(2),
     &    Raw_target, 'Ka2:', Raw_wavelengths(3),
     &    Raw_target, 'Kbeta:', Raw_wavelengths(4)
	write (lun,'(A,F8.3,A,F8.3,A,F5.3,A,F10.1)') 'Range:',
     &    Raw_angbeg, ' to',Raw_angend, '  Step:',Raw_anginc,
     &    '  Time:',Raw_steptime
c
c Write peaks in table format. Write blank line after every 5th peak.
c------------------------------------------------------------------------------
c   d A    |    2T    | 2T corr |    Int     |  I rel |  fwhm   |   h   k   l |
c----------+----------+---------+------------+--------+---------+-------------+
c 0.12345xx| 123.567xx| 123.456 |1234567.90xx| 100.00 | 2.4567xx| -23 -23 -23x|
c------------------------------------------------------------------------------
	write (lun, '(A1)') ' '
	write (lun, '(A)') 'Table 1: Peak data'
	write (lun, '(80A1)') ('-',i=1,79)
	write (lun, '(2A)') '|   d A    |    2T    | 2T corr |    Int  '
     *    ,'   |  I rel |  fwhm   |   h   k   l |'
	write (lun, '(2A)') '|----------+----------+---------+---------'
     *    ,'---+--------+---------+-------------+'
	imax = 0.0
	do i=1,Peak_number
	  if (imax.lt.Peak_area(i)) imax = Peak_area(i)
	end do
	do i=1,Peak_number
	  j = i
	  if (mod(i,5).eq.1 .and. i.gt.1) write (lun,'(2A)')
     *      '|          |          |         |            |        |   '
     *      ,'      |             |'
	  Peak_height(j) = 100.0 * Peak_area(j) / imax
	  write (string,20,iostat=ier) Peak_dspace(j), Peak_dflag(j),
     1      Peak_2theta(j), Peak_2tflag(j), Peak_2tcorr(j),
     2      Peak_area(j), Peak_intflag(j), Peak_height(j),
     3      Peak_fwhm(i), Peak_fwhmflag(i),
     4      Peak_h(j), Peak_k(j), Peak_l(j), Peak_hklflag(j)
   20     format ('|',F8.2,A2,'|',F8.3,A2,'|',F8.3,1x,'|',F10.0,A2,'|',
     1      F7.2,1x,'|'F7.4,A2,'|',3I4,A1,'|')
c
c Write d-spacing to 5 significant digits
	  if (Peak_dspace(j).ge.100.0) then
	    write (string(2:9),'(F8.2)') Peak_dspace(j)
	  else if (Peak_dspace(j).ge.10.0) then
	    write (string(2:9),'(F8.3)') Peak_dspace(j)
	  else if (Peak_dspace(j).ge.1.0) then
	    write (string(2:9),'(F8.4)') Peak_dspace(j)
	  else
	    write (string(2:9),'(F8.5)') Peak_dspace(j)
	  end if
c
c Write intensity to variable format (F10.0 default)
	  if (imax.lt.10000.0) then
	    write (string(34:43),'(F10.2)') Peak_area(j)
	  else if (imax.lt.100.0) then
	    write (string(34:43),'(F10.4)') Peak_area(j)
	  end if
c
c Clear 2T corrected field if 0.0
	  if (Peak_2tcorr(j).eq.0.0) string(24:31) = ' '
c
c Clear hkl field if 000
	  if (Peak_h(j).eq.0 .and. Peak_k(j).eq.0 .and. Peak_l(j).eq.0)
     &      string(66:77) = ' '
c
c Output line to file
	  write (lun, '(A)') string(1:79)
	end do
	write (lun, '(80A1)') ('-',i=1,79)
c
c Write region background table to file
c----------------------     ----------------------
c|    2T   |  Lf Bkg  |     |    2T   |  Rt Bkg  |
c+---------+----------+     +---------+----------+
c| 123.567 | 1234.67V |     | 123.567 | 1234.67V |
c----------------------     ----------------------
	write (lun, '(A1)') ' '
	write (lun, '(A)') 'Table 2: Background data'
	write (lun, '(A,5x,A)') ('----------------------',i=1,2)
	write (lun, '(A,5x,A)') '|    2T   |  Lf Bkg  |',
     &    '|    2T   |  Rt Bkg  |'
	write (lun, '(A,5x,A)') ('+---------+----------+',i=1,2)
	do i=1,Region_number
	  write (lun, 70)
     &      Region_angbeg(i), Region_lfbkg(i), Region_lfbkgflag(i),
     &      Region_angend(i), Region_rtbkg(i), Region_rtbkgflag(i)
   70     format(2('|',F8.3,1x,'|',F8.2,A2,'|',:,5x))
	end do
	write (lun, '(A,5x,A)') ('----------------------',i=1,2)
c
   90   call iclose (lun, jer)
	if (ier.eq.0) ier = jer
	pks_write = ier
	return
	end
*+
c -----------------------------------------------------------------------
c PKS_READ -- reads pertinent information from *.PKS file.
c   Note: Relative intensities are recalculated.
c   Note: Missing 2T's are recalculated from 2Tcorr (or d-spacings)
c   Note: Missing d-spacings are recalculated from 2Tcorr (or 2T)
c   Note: 2Tcorr, area, FWHM, and hkl fields may be missing (zero)
c
	integer function pks_read (filename, irange,
     &     angbeg, angend, anginc)
	implicit none
	character*(*) filename  !(R) peak file to read
        integer irange          !(R) range number to read
	real angbeg             !(W) original starting 2T of raw data (or zero)
	real angend             !(W) original ending 2T of raw data (or zero)
	real anginc             !(W) original step width of raw data (or zero)
c                               ! angbeg, angend, & anginc=0.0 means unknown
*-
	include 'pks_cmn.inc'   !(W) all peak information
	include 'raw_var.inc'   !(R) Raw_valid, Raw_angbeg, Raw_anginc
c
	character line*80
	integer ier, jer, lun, n, table, nrange
        real imax
c
	integer iopen, lens
c
	lun = iopen (filename, ' ', 'OLD', 'SEQUENTIAL',
     &    'FORMATTED', 80, 0, ier)
	if (ier.ne.0) goto 90
c
c table variable: 0=header, 1=peaks, 2=background
	table = 0
	Peak_number = 0
	Region_number = 0
	angbeg = 0.0
	angend = 0.0
	anginc = 0.0
        nrange = 0      ! current range number
c
c Read header information, peak table, and background table.
c
c Siemens Analytical X-Ray Instruments. Copyright 1991. All rights reserved.
c File: \gadds\test.pks                  Created: 12May1991 12:00:00.0
c Wavelengths:
c Range:
c
cTable 1: Peak data
c------------------------------------------------------------------------------
c   d A    |    2T    | 2T corr |    Int     |  I rel |  fwhm   |   h   k   l |
c----------+----------+---------+------------+--------+---------+-------------+
c 0.12345xx| 123.567xx| 123.456 |1234567.90xx| 100.00 | 2.4567xx| -23 -23 -23x|
c------------------------------------------------------------------------------
c
cTable 2: Background data
c----------------------     ----------------------
c|    2T   |  Lf Bkg  |     |    2T   |  Rt Bkg  |
c+---------+----------+     +---------+----------+
c| 123.567 | 1234.67V |     | 123.567 | 1234.67V |
c----------------------     ----------------------
c
   10   read (lun, '(A)', end=80) line
	if (table.eq.0) then
	  if (line(1:4).eq.'Siem') then
            nrange = nrange + 1
            if (nrange.gt.irange) goto 80
          else if (nrange.ne.irange) then
            ! skip processing any unrequested range
	  else if (line(1:4).eq.'File') then
	  else if (line(1:4).eq.'Titl') then
	  else if (line(1:4).eq.'Wave') then
	  else if (line(1:4).eq.'Rang') then
	    read(line,'(6x,F8.3,3x,F8.3,7x,F5.3,7x,F10.1)',err=10)
     &        angbeg, angend, anginc
	  else if (lens(line).eq.0) then
	  else if (line(1:16).eq.'Table 1: Peak da') then
	    table = 1
	  else if (line(1:16).eq.'Table 2: Backgro') then
	    table = 2
          else
            ier = -1
            goto 80
	  end if
	else if (table.eq.1) then
	  if (line(1:4).eq.'----') then
	  else if (line(1:8).eq.'|   d A ') then
	  else if (line(1:4).eq.'|---') then
	  else if (lens(line).eq.0) then
	  else if (line(1:16).eq.'Table 2: Backgro') then
	    table = 2
	  else
	    n = Peak_number + 1
	    read (line,20,iostat=ier) Peak_dspace(n), Peak_dflag(n),
     1        Peak_2theta(n), Peak_2tflag(n), Peak_2tcorr(n),
     2        Peak_area(n), Peak_intflag(n), Peak_height(n),
     3        Peak_fwhm(n), Peak_fwhmflag(n),
     4        Peak_h(n), Peak_k(n), Peak_l(n), Peak_hklflag(n)
   20       format (1x,F8.2,A2,1x,F8.3,A2,1x,F8.3,1x,1x,F10.0,A2,1x,
     1        F7.2,1x,1x,F7.4,A2,1x,3I4,A1,1x)
            if (ier.ne.0) goto 80
            if (Peak_dspace(n).ne.0.0 .or. Peak_2theta(n).ne.0.0)
     1        Peak_number = n
	  end if
	else if (table.eq.2) then
	  if (line(1:4).eq.'----') then
	  else if (line(1:8).eq.'|    2T ') then
	  else if (line(1:4).eq.'+---') then
	  else if (lens(line).eq.0) then
	  else if (line(1:4).eq.'Siem') then
            goto 80     ! Table 2 is over, next range started
	  else
	    n = Region_number + 1
	    read (line, 70, iostat=ier)
     &        Region_angbeg(n), Region_lfbkg(n), Region_lfbkgflag(n),
     &        Region_angend(n), Region_rtbkg(n), Region_rtbkgflag(n)
   70       format(2(1x,F8.3,1x,1x,F8.2,A2,1x,5x))
            if (ier.ne.0) goto 80
            if (Raw_valid) then
	      Region_number = n
              Region_lflimit(n) =
     &          1 + nint((Region_angbeg(n) - Raw_angbeg) / Raw_anginc)
              Region_rtlimit(n) =
     &          1 + nint((Region_angend(n) - Raw_angbeg) / Raw_anginc)
            end if
	  end if
	end if
	goto 10
c
   80   call iclose (lun, jer)
	if (ier.eq.0) ier = jer

        ! requested range not found
        if (ier.eq.0 .and. nrange.lt.irange) ier = -1
c
c Fix up missing d's or 2T's
cc        do n = 1, Peak_number
cc          if (Peak_2theta(n).eq.0.0) then
cc            if (Peak_2tcorr(n).ne.0.0) then
cc              Peak_2theta(n) = Peak_2tcorr(n)
cc            else
cc              if (Peak_dspace(n).gt.0.0)
cc     &          Peak_2theta(n) = Bragg2T( wave, Peak_dspace(n) )
cc            end if
cc          end if
cc          if (Peak_dspace(n).eq.0.0) then
cc            if (Peak_2tcorr(n).ne.0.0) then
cc              Peak_dspace(n) = BraggD( wave, Peak_2tcorr(n) )
cc            else
cc              Peak_dspace(n) = BraggD( wave, Peak_2theta(n) )
cc            end if
cc          end if
cc        end do
c
c Fix up missing relative intensities
        imax = 0.0
        do n=1,Peak_number
          if (imax.lt.Peak_area(n)) imax = Peak_area(n)
        end do
        if (imax.gt.0.0) then
          do n=1,Peak_number
            Peak_height(n) = 100.0 * Peak_area(n) / imax
          end do
        end if
c
c Fix up missing start and end range of original raw data trace
cc        if (Pks_angend.eq.0.0 .and. Peak_number.gt.0) then
cc          Pks_angbeg = Peak_2theta(1) - 0.5
cc          Pks_angend = Peak_2theta(Peak_number) + 0.5
cc        end if
c
c Exit with error status
   90   pks_read = ier
	return
	end