From:	BIGVAX::edu%"toby%xtal.dnet@ttown.msc.edu" 10-MAR-1992 17:32:38.98
To:	snyder@ttown.apci.com
CC:	
Subj:	type sys$input ! FILE(S) [.GSAS]RRW*.* ENCLOSED

Received: From UMNACVX(SMTPUSER) by CERAMICS with Jnet id 7900
          for SNYDER@CERAMICS; Tue, 10 Mar 1992 17:31 EST
Received: from noc.msc.edu by vx.acs.umn.edu; Tue, 10 Mar 92 16:34 CST
Received: from ttown.apci.com by noc.msc.edu (5.65/MSC/v3.0(901107)) id
 AA03099; Tue, 10 Mar 92 16:30:06 -0600
Received: from XTAL.DECnet MAIL11D_V3 by ttown.apci.com (5.57/Ultrix3.0-C) id
 AA24344; Tue, 10 Mar 92 17:29:30 -0500
Date: Tue, 10 Mar 92 17:29:30 -0500
From: toby%xtal.dnet@ttown.msc.edu
Subject: type sys$input ! FILE(S) [.GSAS]RRW*.* ENCLOSED
To: snyder@ttown.apci.com
Message-id: <9203102229.AA24344@ttown.apci.com>
X-Envelope-to: snyder@CERAMICS.BITNET
 
$ copy/log sys$input RRWSAS.FOR;8
$deck/dollars=$EODBHTCOPY**
        SUBROUTINE RRWSAS
C
C  PROGRAMMER:  R. L. Snyder/ ALFRED UNIVERSITY   Jul-1988
C               B. H. Toby / Air Products 10-91, 3-92
C
C       THIS ROUTINE READS FROM EITHER A GSAS HISTOGRAM FILE OR
C       AN ASCII GSAS RAW DATA FILE
C       AND STORES IN THE CORRESPONDING COMMON BLOCKS
C
        INCLUDE 'FILECONV.CBS'
 
        INTEGER*4 IUEXP,IUPRF,IUREF
        CHARACTER*60 TEMP
        INTEGER*4 IHIST
        CHARACTER*4 hhead(99),HTYPE
        CHARACTER*4 label*20
        INTEGER*4 READPRF,REDREFP
        INTEGER*4 EFLAG
        CHARACTER*6 KEY,HSTKEY
        LOGICAL NORM,IANS
        REAL*4 TTOF,YO,YC,YI,YB,YW,CHWDT
        REAL*4 DIFC,DIFA,ZERO
        REAL TMAX, BCOF(12), DMIN
        INTEGER IBAK, NCOF
        INTEGER I, J, K, j1, j3
        CHARACTER*1 digit(0:9)/'0','1','2','3','4','5','6','7','8','9'/
 
        CLOSE( UNIT=IINP )
1       CALL CLEARS( 0)                 !CLEAR THE SCREEN
        REFRSH = .true.
        CALL COLORS(80)                 !INTENSE PRINT
        WRITE(6,'(A)')
        1  ' This routine can read from a GSAS format ASCII data file',
        1  ' or from a GSAS Histogram (Rietveld output) file'
        CALL COLORS(1)
        WRITE(6,'(A)') ' ',
        1       ' Choose the code for the desired operation: ',
        1       '       0) read ASCII file [default]',
        2       '       1) read Observed histogram',
        3       '       2) read Computed histogram',
        3       '       3) read Histogram Background ',
        3       '       4) read Observed-Background Histogram',
        3       '       5) read Computed Histogram w/o Background',
        3       '       6) read Observed-Computed Histogram'
        CALL COLORS(-1)                 !REVERSE VIDEO
        WRITE(6,'(A)') '$Enter a number (0-6): '
        CALL COLORS(1)
        READ (5,'(I)',err=1,end=883) iread
        IF (iread .lt. 0 .or. iread .gt. 6) goto 883
        IF (iread .gt. 0) goto 500
C----------------------------------------------------------------------
C read an ASCII file
        OPEN( UNIT=IINP, FILE=FILINP, STATUS='OLD', ACCESS='SEQUENTIAL',
     1 FORM='FORMATTED', ERR=888 )
C
C  READ THE FILE HEADER...
C
        READ(IINP,1000,ERR=777,END=777) SLID
1000    FORMAT(20A4)
C
C  READ THE SCAN PARAMETERS...
C
        READ(IINP,1001,ERR=777,END=777) NPTS, NREC, BANG(1),
     1  BANG(3), BANG(10)
        BANG(1) = BANG(1) / 100.
        BANG(3) = BANG(3) / 100.
        BANG(2) = (NPTS - 1) * BANG(3) + BANG(1)
1001    FORMAT(6X,2I8,6X,3F10.0)
C
C  START READING THE INTENSITY DATA...
C
        do i=1, npts, 10
          READ(IINP,1200,ERR=777,END=776) (YOBS(j),j=i,i+9)
        enddo
*
        QMAX=YOBS(1)
 776    DO I=1,NPTS
          IF(YOBS(I).GT.qMAX)THEN
            QMAX = MAX(QMAX,YOBS(I))
          ENDIF
        ENDDO
        IF(QMAX .lt. 1000.) THEN
          DO I=1,NPTS
            YOBS(I) = 300.0 * YOBS(I)
          ENDDO
        ENDIF
*
1200    FORMAT(10F8.0)
C all was OK
        GOTO 999
C----------------------------------------------------------------------
C read a histogram file
C----------------------------------------------------------------------
C  Change the file extension to .EXP
500     CALL NAMEXT(FILINP,FILINP,'EXP')
        CALL GETUNIT(IUEXP)
C   Open experiment file
        call OPNEXP('SHARED','KEEP',FILINP,IUEXP,EFLAG)
        if (EFLAG .ne. 0) then
          write (6,'(/2a)') ' Unable to open file: ',FILINP
          write (6,'(a,I)') ' ERROR CODE',EFLAG
          goto 1888
        endif
 
        temp = ' unable to read title'
        call readexp(IUEXP,'      DESCR ',temp)
        read(temp,'(20a4)') SLID
 
C Get # of histograms in experiment
        TEMP = 'NOT FOUND'
        call readexp(IUEXP,' EXPR  NHST ',temp)
        NHIST = 0
        IF (TEMP .NE. 'NOT FOUND') READ(TEMP,'(I5)') NHIST
 
        WRITE (6,'(A,i5,a)') ' There are',NHIST,
        1       ' histogram(s) on this file'
 
C get the histogram flags
C       1: P    Powder
C          S    Single xtal
C          R    soft constraints (type 'RSN ')
C
C       2: N    Neutron
C          X    X-ray
C
C       3: T    Time of flight
C          C    CW
C
C       4: D    Dummy
C               (blank) use in least squares
C          *    exclude from  least squares
C          F,G,I,N      SCD data
C          R    powder data to be read
C          X    not yet processed by POWPREF
C
105     TEMP = 'NOT FOUND'
        do i=1,nhist,12
          write (key,'(i1)') 1+(I-1)/12
          call readexp(IUEXP,' EXPR  HTYP'//key(1:1),temp)
          IF (TEMP .NE. 'NOT FOUND') then
            READ(TEMP,'(2x,12(a4,1x))') (hhead(j),j=I,MIN(I+11,nhist))
          ELSE
            write (6,'(A)') ' Missing EXP  HTYP'//key(1:1)//' card'
            goto 1888
          ENDIF
          DO J=I,MIN(I+11,nhist)
            HTYPE = hhead(j)
            if (HTYPE(1:1) .eq. 'P') THEN
              if (HTYPE(3:3) .eq. 'T') THEN
                label = 'TOF'
              elseif (HTYPE(3:3) .eq. 'C') THEN
                label = 'CW '
              ELSE
                label = ' ? '
              ENDIF
              if (HTYPE(2:2) .eq. 'N') THEN
                label(4:) = ' Neutron'
              ELSEif (HTYPE(2:2) .eq. 'X') THEN
                label(4:) = ' X-ray  '
              ELSE
                label(4:) = '   ?    '
              ENDIF
              if (HTYPE(4:4) .eq. 'D') THEN
                label(12:) = ' Dummy'
              ELSEif (HTYPE(4:4) .eq. '*') THEN
                label(12:) = ' Excluded'
              ELSEif (HTYPE(4:4) .eq. 'R') THEN
                label(12:) = ' To be read'
              ELSEif (HTYPE(4:4) .eq. 'X') THEN
                label(12:) = ' Not Proc''d'
              ELSEif (HTYPE(4:4) .eq. ' ') THEN
                label(12:) = '  (OK)'
              ELSE
                label(12:) = '   ?'
              ENDIF
            ELSE
              label = 'Not Powder'
            ENDIF
           WRITE (6,'(1x,i4,2x,a,4x,a)') J,HTYPE,label
          ENDDO
        ENDDO
 
 
C Get some histogram information from the exp file
110     IF (NHIST .LE. 0) goto 1888
        IF (NHIST .EQ. 1) THEN
          IHIST = 1
        ELSE
          WRITE (6,'(A)') '$Enter histogram number (<0 to exit): '
          READ (5,'(I)') IHIST
          IF (IHIST .LE. 0) goto 1883
          IF (IHIST .GT. NHIST .or. IHIST .EQ. 0) goto 105
        ENDIF
        HTYPE = hhead(IHIST)
        IF (HTYPE(1:1) .ne. 'P') THEN
          write (6,'(A)') ' Histogram is not powder data'
          GOTO 105
        ENDIF
        IF (HTYPE(4:4) .eq. '*') THEN
          write (6,'(A)') ' WARNING: Histogram is excluded from LS'
        ELSEIF (HTYPE(4:4) .eq. 'D') THEN
          write (6,'(A)') ' WARNING: Histogram is dummy-type'
        ELSEIF (HTYPE(4:4) .ne. ' ') THEN
          write (6,'(2A)') ' Histogram is not appropriate status: ',
        1       HTYPE
          GOTO 105
        ENDIF
        KEY = HSTKEY(IHIST)
        TEMP = 'NOT FOUND'
        call readexp(IUEXP,key,temp)
        IF (TEMP .EQ. 'NOT FOUND') THEN
          type *,'HISTOGRAM NOT PRESENT'
          GOTO 105
        ENDIF
 
!       TEMP = 'NONE'
!       call readexp(IUEXP,key//'  HNAM',temp)
!       WRITE (lo,'(4x,2a)')  'Histogram title: ',temp(:LENCH(TEMP))
!       TEMP = 'NONE'
!       call readexp(IUEXP,key//'  HFIL',temp)
!       WRITE (lo,'(4x,2a)')  'Original data file: ',temp(:LENCH(TEMP))
!       TEMP = 'NONE'
!       call readexp(IUEXP,key//' BANK',temp)
!       WRITE (lo,'(4x,2a)')  'bank # ',temp(:LENCH(TEMP))
 
C Get histogram parameters
        TEMP = 'NOT FOUND'
        call readexp(IUEXP,key//' IRAD ',temp)
        IF (TEMP .NE. 'NOT FOUND') THEN
          READ(TEMP,'(I5)') IWAVE
          IF (IWAVE .EQ. 1) THEN
            iele = 1
          ELSEIF (IWAVE .EQ. 2) THEN
            iele = 2
          ELSEIF (IWAVE .EQ. 3) THEN
            iele = 4
          ELSEIF (IWAVE .EQ. 4) THEN
            iele = 5
          ENDIF
        ELSE
          iele = 6
        ENDIF
 
        TEMP = 'NOT FOUND'
        call readexp(IUEXP,key//' ICONS',temp)
        IF (TEMP .EQ. 'NOT FOUND') THEN
          WRITE(6,'(A)')
        1       ' Warning -- no HST hh CONS record on .EXP file'
        ELSEif (HTYPE(3:3) .ne. 'T') THEN
          READ(TEMP,'(3F10.0,10x,f10.0)') wave(1,6),wave(2,6),zero
        1       ,wave(5,6)
          if (wave(5,6) .gt. 0.0) wave(5,6) = 1./wave(5,6)
          if (wave(1,6) .eq. wave(2,6) .or. wave(5,6) .eq. 0) iele = 6
C convert to degrees
          ZERO = ZERO/100.
        ELSE
          READ(TEMP,'(3F10.0)') DIFC,DIFA,ZERO
C convert to millsec
          DIFC = DIFC/1000.
          DIFA = DIFA/1000.
          ZERO = ZERO/1000.
          wave(1,6) = 0
          wave(2,6) = 0
          wave(5,6) = 0
        ENDIF
 
        IF (ZERO .ne. 0. .and. HTYPE(3:3) .ne. 'T') THEN
          WRITE(6,'(A)') '$Apply zero correction ([Y]/N)? '
          READ (5,'(A1)',err=1,end=883) ians
          if (ians .eq. 'n' .or. ians .eq. 'N') ZERO = 0.0
        ENDIF
 
C Open HISTOGRAM file
        CALL GETUNIT(IUPRF)
        call OPNPRF(IUEXP,ihist,NCHAN,'SHARED',.FALSE.,IUPRF)
        call readexp(IUEXP,key//' CHANS',temp)
        READ(TEMP,'(5I10)') ihold,iclmp,nchans,ichk,mchans
 
        IF (iread .ge. 3 .and. iread .le. 5) THEN
C get the max d value
          TMAX = 180.0
          IF (HTYPE(3:3).EQ.'T' ) THEN
            call readexp(IUEXP,key//'I ITYP',temp)
            READ(temp,'(15X,F10.4)') TMAX
            TMAX = 180.0/TMAX
          ENDIF
C Get BACKGROUND parameters
          DO I=1,12
            BCOF(I) = 0.0
          END DO
          call readexp(IUEXP,key//'BAKGD ',temp)
          READ(temp,'(2I5)') IBAK,NCOF
          j1 = 1
          DO j3 = 1, (NCOF-1)/4+1
            call readexp(IUEXP,key//'BAKGD'//digit(j3),temp)
            READ(temp,'(4E15.6)') (BCOF(I),I=j1,MIN(j1+3,NCOF))
            j1 = j1+4
          ENDDO
        ENDIF
 
C READ FROM HISTOGRAM FILE
        norm = .false.
        NPTS = 0
        DO IREC=1,MIN(MNP,nchans)
          I = READPRF(IUPRF,IREC,ICODE,TTOF,YO,YC,YI,YB,YW,CHWDT,MIN1,MIN2)
          IF (I .NE. 0) GOTO 1888
          IF (HTYPE(3:3) .eq. 'T') THEN
C convert from microsec
            TTOF = TTOF / 1000.
C compute the d-space from time of flight
            IF (DIFA.NE.0.0 ) THEN
              XOBS(NPTS+1) = 2.0*(TTOF-ZERO)/
        1               (DIFC+SQRT(DIFC**2+4.0*DIFA*(TTOF-ZERO)))
            ELSE
              XOBS(NPTS+1) = (TTOF-ZERO)/DIFC
            ENDIF
          ELSE
C convert from centidegrees
            TTOF = TTOF / 100.
            XOBS(NPTS+1) = TTOF-ZERO
          ENDIF
 
          IF (iread .ge. 3 .and. iread .le. 5) THEN
C compute the background
            CALL CALCBAK(HTYPE,TMAX,DIFC,wave(1,6),IBAK,NCOF,
        1       BCOF,TTOF,YB1)
C add the fixed and computed background
            yb = yb + yb1
          ENDIF
          IF (iread .eq. 1) THEN
            YOBS(NPTS+1) = YO
          ELSEIF (iread .eq. 2) THEN
            YOBS(NPTS+1) = YC
          ELSEIF (iread .eq. 3) THEN
            YOBS(NPTS+1) = YB
          ELSEIF (iread .eq. 4) THEN
            YOBS(NPTS+1) = YO - YB
          ELSEIF (iread .eq. 5) THEN
            YOBS(NPTS+1) = YC - YB
          ELSEIF (iread .eq. 6) THEN
            YOBS(NPTS+1) = YO - YC
          ENDIF
 
          IF (IREC .eq. 1 .and. YI .ne. 1. .and. YI .gt. 0) THEN
            WRITE(6,'(A)') '$Normalize by incident spectrum (Y/[N]): '
            read (5,'(A1)',err=1883) ians
            if (ians .eq. 'Y' .or. ians .eq. 'y') norm = .true.
          ENDIF
          IF (norm .and. YI .gt. 0.0) THEN
            YOBS(NPTS+1) = YOBS(NPTS+1)*YI
          ENDIF
          NPTS = NPTS + 1
        ENDDO
        BANG(1) = XOBS(1)
        BANG(2) = XOBS(NPTS)
C the step size may be approximate -- lets hope that VAXCONV can handle this
        BANG(3) = (BANG(2) - BANG(1)) / (NPTS - 1.)
        CLOSE (IUEXP)
        CLOSE (IUPRF)
        GOTO 999
1883    CLOSE (IUEXP)
        CLOSE (IUPRF)
        goto 883
1888    CLOSE (IUEXP)
        CLOSE (IUPRF)
        goto 888
C----------------------------------------------------------------------
777     IERR = 1
        GOTO 999
888     IERR = 2
        GOTO 999
883     IERR = 3
999     RETURN
        END
$EODBHTCOPY**