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**