C----------------------------------------------------------------------
	SUBROUTINE READIPNS(dname,pname,filform,nfil,STATUS)
C----------------------------------------------------------------------
C Read in file in ARGONNE/IPNS GLASS format -- RDF/SOQ 
C                              TOF format -- ASCII TOF format
C                              TOF RIETVELD format -- *LS.PLT format
C----------------------------------------------------------------------
	LOGICAL status
	INTEGER nfil
	character*(*) dname,pname,filform
	CHARACTER*4 thisform/'IPNS'/
	CHARACTER*20 testform

	REAL PI
	PARAMETER (PI=3.141592654)

	include 'shared$common'
	include '[-.shared]flabel.common'
C
	INTEGER ln,npt
	real YOBS(10),YFIT(10),DV(10)

	LOGICAL exist

C Assume that the read was not successful
	STATUS = .false.

C check the file format qualifier, if one exists
	IF (filform .ne. ' ') THEN
	  I = LENCH(filform)
	  IF (filform(:I) .ne. thisform(:I)) RETURN
	ENDIF
C---------------------------------------------------------------------
C Is there a datafile name?
	IF (dname .eq. ' ') RETURN
	IF (filform(:7) .eq. 'IPNSSOF') GOTO 311
	IF (filform(:7) .eq. 'IPNSTOF') GOTO 211
	IF (filform(:8) .eq. 'IPNSRIET') GOTO 411

C---------------------------------------------------------------------
C   Does the file exist as an RDF file?
	INQUIRE(FILE=dname,EXIST=exist,defaultfile='.RDF')
	IF (.not. exist) goto 311
C try to start reading it
	OPEN(UNIT=IN,TYPE='OLD',ERR=311,FORM='UNFORMATTED',
	1	FILE=dname,defaultfile='.RDF',readonly,shared)
	READ(IN,err=311) ATITLE(nfil)
	READ(IN,err=311)	! HISTRY
	read (IN,err=311) nmin,nmax,delr,rho
C is the header reasonable?
	IF (nmax .le. nmin) goto 311
	IF (delr .gt. 0.1 .or. delr .le. 0.0) goto 311
	IF (rho  .gt. 10000. .or. rho .lt. 0.0) goto 311
	npt = nmax - nmin + 1
	IF (filform .eq. 'IPNSRHO') THEN
	  ityp = 1
	ELSEIF (filform .eq. 'IPNSGOFR') THEN
	  ityp = 2
	ELSEIF (filform .eq. 'IPNSNOFR') THEN
	  ityp = 3
	ELSEif (prompt) then 
	  write(lo,'(a)') ' Choose the variable: ','   1 rho(r)',
	2	'   2  g(r)','   3  n(r)','   4  S(Q)',
	3	'$(Default=1): '
	  read (li,'(I)',END=999) ityp
	  if (ityp .eq. '4') goto 311
	ELSE
	  ityp = 1
	ENDIF	
	MONO(nfil) = .true.
	WAVE(nfil) = 1. 
	WKA1(nfil) = 1. 
	WKA2(nfil) = 1. 
	if (ityp .eq. 2) then 
	  READ (IN,err=311) (BX(i,nfil),y1,B(i,nfil),y2,I=1,NPT)
	  FYTYP(nfil) = 7
	elseif (ityp .eq. 3) then 
	  READ (IN,err=311) (BX(i,nfil),y1,y2,B(i,nfil),y1,y2,I=1,NPT)
	  FYTYP(nfil) = 8
	else
	  ityp = 1
	  READ (IN,err=311) (BX(i,nfil),B(i,nfil),y1,y2,I=1,NPT)
	  FYTYP(nfil) = 6
	endif
	START(nfil) = BX(1,nfil)
	FINISH(nfil) = BX(npt,nfil)
C check for even spacing between points
	DELTA(nfil) = (BX(npt,nfil) - BX(1,nfil))/(npt-1.)
	if (DELTA(nfil) .ne. (BX(2,nfil)-BX(1,nfil))) DELTA(nfil) = 0.0
	NBKGP(nfil) = 0
	NPKS(nfil) = 0
	DO I=1,npt
	  BESD(I,nfil)= 1.0
	ENDDO
	FXTYP(nfil) = XTYP_R
	INQUIRE(UNIT=IN,NAME=cfilnam(nfil))
	close(unit=IN)
	ln = lench(cfilnam(nfil))
	write (LO,'(A)') ' File '//cfilnam(nfil)(:ln)
	1		//' read in ARGONNE RDF format'
	IF (ityp .eq. 1) THEN
	  filform = 'IPNSRHO'
	ELSEIF (ityp .eq. 2) THEN
	  filform =  'IPNSGOFR'
	ELSEIF (ityp .eq. 3) THEN
	  filform = 'IPNSNOFR'
	ENDIF
	NPOINT(nfil) = npt 
	STATUS = .true.
	RETURN
C-------------------------------------------------------------------------
c   Does the file exist as an SOQ file?
311	close(unit=IN)
	INQUIRE(FILE=dname,EXIST=exist,defaultfile='.soq')
	IF (.not. exist) goto 211
	OPEN(UNIT=IN,TYPE='OLD',ERR=211,FORM='UNFORMATTED',
	1	FILE=dname,defaultfile='.soq',readonly,shared)
	READ(IN,err=211) ATITLE(nfil)
	READ(IN,err=211)	! HISTRY
	read (IN,err=211) nmin,nmax,delr,rho,dummy
C is the header reasonable?
	IF (nmax .le. nmin) goto 211
	IF (delr .gt. 0.1 .or. delr .lt. 0.0) goto 211
	IF (rho  .gt. 10000. .or. rho .lt. 0.0) goto 211
	IF (filform .eq. 'IPNSSOFQ') THEN
	  ityp = 1
	ELSEIF (filform .eq. 'IPNSSOFD') THEN
	  ityp = 2
	ELSEif (prompt) then 
!	  ln = lench(cfilnam(nfil))
!	  write(lo,'(a)') '$Read file '//cfilnam(nfil)(:ln)//
!	1	' in IPNS GLASS SOQ format? ([Y]/N): '
!	  read (li,'(A)',END=999) ians
!	  if (ians .eq. 'N' .or. ians .eq. 'n') goto 211
	  write(lo,'(a)') ' Choose the x-axis units: ','   1 Q',
	2	'   2  d-space','$(Default=1): '
	  read (li,'(I)',END=999) ityp
	ELSE
	  ityp = 1
	ENDIF	
	MONO(nfil) = .true.
	WAVE(nfil) = 1. 
	WKA1(nfil) = 1. 
	WKA2(nfil) = 1. 

	npt = nmax - nmin + 1
	READ (IN,err=211) (BX(i,nfil),B(i,nfil),BESD(i,nfil),I=npt,1,-1)
	if (ityp .ne. 2) then 
	  ityp = 1
	  do i=1,npt
	    BX(I,nfil) = 2.0*3.141592654/BX(I,nfil)
	  enddo
	  FXTYP(nfil) = XTYP_Q
	ELSE
	  FXTYP(nfil) = XTYP_D
	ENDIF
	START(nfil) = BX(1,nfil)
	FINISH(nfil) = BX(npt,nfil)
	DELTA(nfil) = 0.0
	NBKGP(nfil) =0
	NPKS(nfil) = 0
	DO I=1,npt
	  BESD(I,nfil)= 1.0
	ENDDO
	FYTYP(nfil) = 5
	INQUIRE(UNIT=IN,NAME=cfilnam(nfil))
	ln = lench(cfilnam(nfil))
	close(unit=IN)
	write (LO,'(A)') ' File '//cfilnam(nfil)(:ln)
	1		//' read in ARGONNE SOQ format'
	IF (ityp .eq. 1) THEN
	  filform = 'IPNSSOFQ'
	ELSEIF (ityp .eq. 2) THEN
	  filform = 'IPNSSOFD'
	ENDIF
	NPOINT(nfil) = npt 
	STATUS = .true.
	RETURN
C----------------------------------------------------------------------
211	close(unit=IN)
C   Does the file exist as an ASCII TOF raw data file?
	INQUIRE(FILE=dname,EXIST=exist,defaultfile='.TOF')
	IF (.not. exist) goto 411
	OPEN(UNIT=IN,FILE=dname,defaultfile='.TOF',TYPE='OLD',
	1	ERR=411,FORM='FORMATTED',readonly,shared)
	READ(IN,'(A55,I4,3X,I2,5X,I4)',err=411) LAB1,npt
	if (npt .le. 0) goto 411
	READ(IN,'(3f8.3,3g12.6,2i5)',err=411) 
	1	angle,dist,dists,difc,difa,zp,ICW,ITMIN
	DT=ICW
	TMIN=ITMIN
	ityp = 1
	IF (filform .eq. 'IPNSTOFT') THEN
	  ityp = 1
	ELSEIF (filform .eq. 'IPNSTOFL') THEN
	  ityp = 2
	ELSEIF (filform .eq. 'IPNSTOFQ') THEN
	  IF (angle .ne. 0.0 .and. angle .ne. 180.) ityp = 3
	ELSEIF (filform .eq. 'IPNSTOFD') THEN
	  IF (angle .ne. 0.0 .and. angle .ne. 180.) ityp = 4
	ELSEif (prompt) then 
!	  ln = lench(cfilnam(nfil))
!	  write(lo,'(a)') '$Read TOF data file '//cfilnam(nfil)(:ln)//
!	1	'? ([Y]/N): '
!	  read (li,'(A)',END=999) ians
!	  if (ians .eq. 'N' .or. ians .eq. 'n') goto 411
232	  IF (angle .eq. 0.0 .or. angle .eq. 180.) then
	    write(lo,'(a)') ' Choose the x-axis: ','   1 TOF (ms)',
	2	'   2  wavelength',
	3	'$(Default=1): '
	    read (li,'(I)',END=999) ityp
	    if (ityp .gt. 2) goto 232
	  ELSE
	    write(lo,'(a)') ' Choose the x-axis: ','   1 TOF (ms)',
	2	'   2  wavelength','   3  Q','   4  D-space',
	3	'$(Default=1): '
	    read (li,'(I)',END=999) ityp
	    if (ityp .gt. 4) goto 232
	  ENDIF
	ENDIF	
	npt = MIN(M$PTS,npt)
	READ(IN,'(10f8.0)',err=411) (B(I,nfil),I=1,npt)
	MONO(nfil) = .true.
	WAVE(nfil) = 1. 
	WKA1(nfil) = 1. 
	WKA2(nfil) = 1. 

	IF (ityp .eq. 0 .OR. ITYP .EQ. 1) THEN
	  ityp = 1
          DO I=1,npt
            BX(i,nfil) = ((I-1.0)*dt+tmin)/1000.
          ENDDO
	  FXTYP(nfil) = XTYP_TOF
	ELSEIF (DIFA .eq. 0) then 
C compute WAVELENGTH directly from TOF and distance
          DO I=1,npt
            BX(i,nfil) = 3956.05*((I-1.0)*dt+tmin)*1.0E-6/DIST
	    IF (ityp .eq. 4) then 
              BX(i,nfil) = 0.5*BX(i,nfil)/SIN(PI*angle/360.)
	    ELSEIF (ityp .eq. 3) then 
              BX(i,nfil) = 4.*PI*SIN(PI*angle/360.)/BX(i,nfil)
	    ENDIF
          ENDDO
	  FXTYP(nfil) = XTYP_LAM
	  IF (ityp .eq. 4) then 
	     FXTYP(nfil) = XTYP_D
	  ELSEIF (ityp .eq. 3) then 
	     FXTYP(nfil) = XTYP_Q
	  ENDIF
	ELSE
C compute d-space from calibration constants
          DO I=1,npt
	    tof = (I-1.0)*dt+tmin
C TOF => D-SPACE
	    CALL QUADTD(TOF,DIFA,DIFC,ZP,BX(i,nfil))
	    IF (ityp .eq. 2) then 
C D-SPACE ==> WAVELENGTH
              BX(i,nfil) = 2.0 * BX(i,nfil) * SIN(PI*angle/360.)
	    ELSEIF (ityp .eq. 3) then 
C D-SPACE ==> Q
              BX(i,nfil) = 2.*PI/BX(i,nfil)
	    ENDIF
          ENDDO
	  TYPE *,'X AXIS IS CALIBRATED'
	  FXTYP(nfil) = XTYP_D
	  IF (ityp .eq. 2) then 
	    FXTYP(nfil) = XTYP_LAM
	  ELSEIF (ityp .eq. 3) then 
	    FXTYP(nfil) = XTYP_Q
	  ENDIF
	ENDIF

	START(nfil) = BX(1,nfil)
	FINISH(nfil) = BX(npt,nfil)
	DELTA(nfil) = 0.0
	NBKGP(nfil) = 0
	NPKS(nfil) = 0
	DO I=1,npt
	  BESD(I,nfil) = max(sqrt(b(i,nfil)),1.)
	ENDDO
	FYTYP(nfil) = 1
	INQUIRE(UNIT=IN,NAME=cfilnam(nfil))
	close(unit=IN)
	ln = lench(cfilnam(nfil))
	write (LO,'(A)') ' File '//cfilnam(nfil)(:ln)
	1		//' read in ARGONNE TOF format'
	IF (ityp .eq. 1) THEN
	  filform = 'IPNSTOFT'
	ELSEIF (ityp .eq. 2) THEN
	  filform =  'IPNSTOFL'
	ELSEIF (ityp .eq. 3) THEN
	  filform =  'IPNSTOFQ'
	ELSEIF (ityp .eq. 4) THEN
	  filform =  'IPNSTOFD'
	ENDIF
	NPOINT(nfil) = npt
	STATUS = .true.
	RETURN
C----------------------------------------------------------------------
C----------------------------------------------------------------------
411	CLOSE(unit=IN)
C   Does the file exist as a Rietveld difference file?
	INQUIRE(FILE=dname,EXIST=exist,defaultfile='.PLT')
	IF (.not. exist) goto 511
	OPEN(UNIT=IN,FILE=dname,defaultfile='.PLT',TYPE='OLD',
	1	ERR=511,FORM='FORMATTED',readonly,shared)
	NPKS(nfil) = 0
C   READ D SPACINGS OF BRAGG REFLECTIONS FROM ALL SAMPLE PHASES IN
C   SEQUENCE FROM UNIT IN
	K=0
	DO I=1,10000
	  J=K+1
	  K=10*I
	  READ(IN,'(8F10.5)',ERR=511) DV
	  DO L=J,K
           IF(DV(L-J+1).EQ.0.0) GO TO 415
	    IF (L .LE. m$pks) THEN 
	      TT(L,nfil) = DV(L-J+1)
	      PKHT(L,nfil) = 100.
	      NPKS(nfil) = L
	    ENDIF
	  ENDDO
	ENDDO
C   READ MINIMUM (IMIN) AND MAXIMUM (IMAX) TIME CHANNEL NUMBERS AND
C   NUMBER OF DATA POINTS (npt), AND NUMBER OF BRAGG REFLECTIONS (NPHZ%)
C   FOR EACH SAMPLE PHASE FROM UNIT IN
415	READ(IN,'(8I10)',err=511) IMIN,IMAX,npt,NPHZ1,NPHZ2,NPHZ3,NPHZ4
	ityp = 1
	IF (filform(9:9) .eq. 'T') THEN
	  ityp = 1
	ELSEIF (filform(9:9) .eq. 'L') THEN
	  ityp = 2
	ELSEIF (filform(9:9) .eq. 'Q') THEN
	  ityp = 3
	ELSEIF (filform(9:9) .eq. 'D') THEN
	  ityp = 4
	ELSEIF (prompt) then 
!	  ln = lench(cfilnam(nfil))
!	  write(lo,'(a)') '$Read TOF Rietveld file '//cfilnam(nfil)(:ln)//
!	1	'? ([Y]/N): '
!	  read (li,'(A)',END=999) ians
!	  if (ians .eq. 'N' .or. ians .eq. 'n') goto 511
332	  write(lo,'(a)') ' Choose the x-axis: ','   1 TOF (ms)',
	2	'   2  wavelength','   3  Q','   4  D-space',
	3	'$(Default=1): '
	  read (li,'(I)',END=999) ityp
	  if (ityp .gt. 4) goto 332
	ENDIF

	JTYP = 1
	IF (filform(10:10) .eq. 'O') THEN
	  Jtyp = 1
	ELSEIF (filform(10:10) .eq. 'C') THEN
	  Jtyp = 2
	ELSEIF (filform(10:10) .eq. 'D') THEN
	  Jtyp = 3
	ELSEIF (prompt) then 
333	  write(lo,'(a)') ' Choose the y-axis: ','   1  Obs',
	2	'   2  calc','   3  Obs-calc',
	3	'$(Default=1): '
	  read (li,'(I)',END=999) Jtyp
	  if (jtyp .gt. 3) goto 333
	ENDIF
C   READ YOBS (YOBS) AND YCALC (YFIT) VALUES, BOTH SANS BACKGROUND, FROM
C   UNIT IN
	MONO(nfil) = .true.
	WAVE(nfil) = 1. 
	WKA1(nfil) = 1. 
	WKA2(nfil) = 1. 
	K=0
	DO I=1,10000
	  J=K+1
	  K=10*I
	  READ(IN,'(10F8.0)',err=511) YOBS,YFIT
	  DO L=J,K
	    IF(YOBS(L-j+1).EQ.-100001.0
	1	.AND. YFIT(L-J+1).EQ.-100001.0) GO TO 420
	    IF (L .LE. m$pts) THEN 
	      IF (jtyp .eq. 2) then 
		B(L,nfil) = YFIT(L-J+1)
	      ELSEIF (jtyp .eq. 2) then 
		B(L,nfil) = YOBS(L-J+1) - YFIT(L-J+1)
	      ELSE
		B(L,nfil) = YOBS(L-J+1) 
	      ENDIF
	    ENDIF
	  ENDDO
	ENDDO
C   READ DIFC, DIFA, ZERO POINT (ZP), CHANNEL WIDTH (DT), TIME
C   DELAY (TMIN) AND MINIMUM (YMIN) AND MAXIMUM (YMAX) COUNTS FROM UNIT IN
420	READ(IN,'(8F10.0)',err=511) DIFC,DIFA,ZP,DT,TMIN,YMIN,YMAX
C compute d-space from calibration constants
	DO I=1,npt
	  TOF = (IMIN + I - 1.)*DT + TMIN
	  IF (ityp .eq. 0 .OR. ityp .eq. 1) then
	    BX(I,nfil) = TOF/1000.
	  ELSE
	    CALL QUADTD(TOF,DIFA,DIFC,ZP,BX(i,nfil))
	    IF (ityp .eq. 2) then 
              BX(i,nfil) = 2.0 * BX(i,nfil) * SIN(PI*angle/360.)
	    ELSEIF (ityp .eq. 3) then 
              BX(i,nfil) = 2.*PI/BX(i,nfil)
	    ENDIF
	  ENDIF
       ENDDO
	FXTYP(nfil) = XTYP_TOF
	IF (ityp .eq. 4) then 
	  FXTYP(nfil) = XTYP_D
	ELSEIF (ityp .eq. 2) then 
	  FXTYP(nfil) = XTYP_LAM
	ELSEIF (ityp .eq. 3) then 
	  FXTYP(nfil) = XTYP_Q
	ENDIF
	NPOINT(nfil) = npt-1
	START(nfil) = BX(1,nfil)
	FINISH(nfil) = BX(npt,nfil)
	DELTA(nfil) = 0.0
	NBKGP(nfil) = 0
	YMIN = B(1,NFIL)
	DO I=1,npt
	  YMIN = MIN(YMIN,B(I,NFIL))
	  BESD(I,nfil) = 1.
	ENDDO
	if (prompt) then 
	  write(lo,'(a,1pg13.4,A)') '$Offset to add to Y values (default=',
	1	max(0.,-ymin),') ? ' 
	  read (li,'(q,f13.4)',ERR=999) i,yoff
	  IF (i .eq. 0) yoff = max(0.,-ymin)
	ELSE
	  yoff = max(0.,-ymin)
	ENDIF
	if (YOFF .ne. 0) then 
	  DO I=1,npt
	    B(I,NFIL) = B(I,NFIL) + YOFF
	  ENDDO
	ENDIF
	FYTYP(nfil) = 3
	INQUIRE(UNIT=IN,NAME=cfilnam(nfil))
	CLOSE(unit=IN)
	ln = lench(cfilnam(nfil))
	write (LO,'(A)') ' File '//cfilnam(nfil)(:ln)
	1		//' read in ARGONNE TOF RIETVELD format'
	filform = 'IPNSRIET'
	IF (ityp .eq. 1) THEN
	  filform(9:9) =  'T'
	ELSEIF (ityp .eq. 2) THEN
	  filform(9:9) =  'L'
	ELSEIF (ityp .eq. 3) THEN
	  filform(9:9) =  'Q'
	ELSEIF (ityp .eq. 4) THEN
	  filform(9:9) =  'D'
	ENDIF
	IF (Jtyp .eq. 1) THEN
	  filform(10:10) =  'O'
	ELSEIF (Jtyp .eq. 2) THEN
	  filform(10:10) =  'C'
	ELSEIF (Jtyp .eq. 3) THEN
	  filform(10:10) =  'D'
	ENDIF
	STATUS = .true.
	RETURN
C----------------------------------------------------------------------
511	CLOSE(unit=IN)
	RETURN
C----------------------------------------------------------------------
999	CLOSE(unit=IN)
	STATUS = .false.
	RETURN
	END
C----------------------------------------------------------------------
      SUBROUTINE QUADTD(D,A,B,C,E)
C
C  ROUTINE UTILIZING THE ANALYTICAL VERSION OF THE QUADRATIC
C  FORMULA TO DETERMINE E USING THE POSITIVE RESULT OF THE SQUARE
C  ROOT OF Q (SEE BELOW)
C
C  IF A=0, ROUTE TO STATEMENT 10
      IF(A.EQ.0) GO TO 10
      TA=2.0*A
      FAC=2.0*TA*(C-D)
      BSQ=B*B
      Q=SQRT(BSQ-FAC)
      E=(-B+Q)/TA
      GO TO 15
   10 CONTINUE
      E=(D-C)/B
   15 CONTINUE
      RETURN
      END
