C***** PROGRAM TO MATCH ELEMENTAL PEAKS FROM XRF ********************
C
C       COPYRIGHT 1982 GENERAL ELECTRIC CO.
C       WRITTEN MARCH 1982 M.F.GARBAUSKAS
C
C        PROGRAM USES PART OF NIH DATA BASE PACKED INTO
C        DIRECT ACCESS FILE 'XF1:LAMBDA.NIH' BY PROGRAM
C        'PACKNIH'
C
C        CUTOFF CONSTRAINTS:
C         MINIMUM RELATIVE INTENSITY: 0.06
C         MAXIMUM NUMBER OF LINES: 75
C         MAXIMUM WAVELENGTH: 25 ANGSTROMS
C
C
C        LOGICAL UNIT ASSIGNMENTS:
C           LSPEC:  FILE FOR INPUT ENERGIES FROM 'SPECPLOT'
C           LNIH:  'POWD:[PUBLIC.BASE.XRF]LAMBDA.NIH'  ELEMENTAL DATA BASE
C           LRES:  RESULTS FILE
C           LUNI:  INPUT
C           LUNP:  PROMPTS
C           LUNO:  OUTPUT
C           LSCR:  SCRATCH FILE 'QUAL.TMP'
C           LOUT:  SCRATCH FILE USED FOR OUTPUT
C
C
C         EXTERNAL SUBROUTINES CALLED:
C          MINMAX   FIND   OUTPUT   BSORT   LIMIT  ZSYM
C	GET2D
C
C      MODIFIED 1-FEB-85 JDW TO READ DIFFRAC-11 *.DIF FILES
C                  ON UNIT LSPEC
C      MODIFIED 13-SEP-87 RLS TO RUN ON VAX
C
C      TASKBUILD OPTIONS:
C         MAXBUF=1100
C         ACTFIL=5
C         ASG=XF1:6
C         DL1:[1,1]FORLIB/DL
C******************************************************************
	DIMENSION THETA(100),WAVE(100),HEIGHT(100)
	BYTE DSN(30),DUMMY(80),FILNAM(30)
	DATA LSPEC/1/,LNIH/2/,LRES/3/,LUNI/5/,LUNP/5/,LUNO/5/
	DATA LSCR/4/,LOUT/6/
C	ASIN(X) = ATAN( X / SQRT( 1.0 - X*X ) )
100	IQ=0
	WRITE(LUNP,9000)
9000	FORMAT('$ Enter filename of DIFFRAC 11 d-I file>')
	READ(LUNI,9010,END=999)IQ,DSN
9010	FORMAT(Q,30A1)
	DSN(IQ+1)=0
	IF(IQ.EQ.0)GO TO 100
C	 OPEN(UNIT=LSPEC,NAME=DSN,TYPE='OLD',SHARED,READONLY)
C--JDW---open binary file
c	CALL FILOPN(1)
c	CLOSE(UNIT=1)
C	CALL SCOPY(DSN,FILNAM)
C	CALL CONCAT(DSN,'.DIF',DSN)
C	CALL DEVUIC(DSN,1)
	OPEN(UNIT=LSPEC,NAME=DSN,TYPE='OLD',SHARED,READONLY,
     &   FORM='UNFORMATTED',ACCESS='DIRECT',RECORDSIZE=1,
     &   ASSOCIATEVARIABLE=LSPAS)
C
150	IQ=0
C	WRITE(LUNP,9020)
C9020	FORMAT('$ ENTER FILENAME FOR RESULTS>')
C	READ(LUNI,9010,END=999)IQ,DSN
C	DSN(IQ+1)=0
C	IF(IQ.EQ.0)GO TO 150
C	OPEN(UNIT=LRES,NAME=DSN,TYPE='NEW')
C
c	CALL CONCAT(FILNAM,'.LST',FILNAM)
D	WRITE(LUNP,841) LRES,FILNAM
D841	FORMAT(I3,30A1,2X,30A1)
	OPEN(UNIT=LRES,NAME=FILNAM,TYPE='NEW')
C   
C---READ HEADER INFORMATION AND TRANSFER TO RESULTS FILE
C       GET 2D FOR CRYSTAL
C  
C	READ(LSPEC,9030,END=999)DUMMY
C
	LSPAS=1
	DUMMY(1)=' '
	READ(LSPEC'LSPAS) (DUMMY(I),I=2,5)	! Sample ID
	READ(LSPEC'LSPAS) (DUMMY(I),I=6,9)
	READ(LSPEC'LSPAS) (DUMMY(I),I=10,13)
	READ(LSPEC'LSPAS) (DUMMY(I),I=14,17)
	DUMMY(18)=' '
	DUMMY(19)=' '
	READ(LSPEC'LSPAS) (DUMMY(I),I=20,23)	! Date of run
	READ(LSPEC'LSPAS) (DUMMY(I),I=24,27)
	READ(LSPEC'LSPAS) (DUMMY(I),I=28,31)
	DUMMY(32)=' '
	DUMMY(33)=' '
	READ(LSPEC'LSPAS) (DUMMY(I),I=34,37)	! Time of run
	READ(LSPEC'LSPAS) (DUMMY(I),I=38,41)
	DUMMY(42)=' '
	DUMMY(43)=' '
	READ(LSPEC'LSPAS) (DUMMY(I),I=44,47)	! Operator's name
	READ(LSPEC'LSPAS) (DUMMY(I),I=48,51)
	READ(LSPEC'LSPAS) (DUMMY(I),I=52,55)
	DO 800, I=56,80
800	DUMMY(I)=' '
C
9030	FORMAT(80A1)
	WRITE(LUNO,9030)DUMMY
	WRITE(LRES,9030)DUMMY
C
C	READ(LSPEC,9035,END=999)(DUMMY(I),I=1,62),TWOD
C9035	FORMAT(62A1,F6.3)
C	WRITE(LUNO,9035)(DUMMY(I),I=1,62),TWOD
C	WRITE(LRES,9035)(DUMMY(I),I=1,62),TWOD
C	READ(LSPEC,9030,END=999)DUMMY
C	WRITE(LUNO,9030)DUMMY
C	WRITE(LRES,9030)DUMMY
	CALL GET2D(LUNI,LUNP,TWOD)
9035	FORMAT(43(' '),'TWO D OF CRYSTAL = ',F6.3)
	WRITE(LUNO,9035) TWOD
	WRITE(LRES,9035) TWOD
	WRITE(LUNO,9030) (' ',I=1,80)
	WRITE(LRES,9030) (' ',I=1,80)
	READ(LSPEC'19) WAVLEN
	READ(LSPEC'21) NPKS
C
	WRITE(LUNP,9300)
9300	FORMAT('$ ENTER ERROR IN 2-THETA (DEFAULT:.10)>')
	READ(LUNI,9310,END=999)D2T
9310	FORMAT(F6.0)
	IF(D2T.EQ.0)D2T=.1
	WRITE(LRES,9311)D2T
9311	FORMAT(1X,'ERROR(DEG)=',F6.3)
	WRITE(LRES,9315)
9315	FORMAT('0',2X,'N',3X,'2THETA',3X,'LAMBDA',3X,'HEIGHT'
     %	 ,5X,'ELEMENTAL LINES MATCHED')
C	DO 280 I=1,3
C	READ(LSPEC,9030,END=999)DUMMY
C280	CONTINUE
C
C---READ DATA FROM INPUT FILE    NL=LINE COUNTER
C
	NL=0
C	DO 290 I=1,100
C285	READ(LSPEC,9140,END=300)THETA(I),WAVE(I),HEIGHT(I)
C9140	FORMAT(4X,2F9.0,9X,F9.0)
C	IF(WAVE(I).EQ.0)GO TO 285
C
	LSPAS = 41
	DO 290, I=1,NPKS
	READ(LSPEC'LSPAS) DSPACE
	THETA(I) = 2. * ASIN( WAVLEN / DSPACE / 2.0 ) * 180.0 / 3.14159
	WAVE(I) = TWOD * WAVLEN / DSPACE / 2.0   
	READ(LSPEC'LSPAS,ERR=300) HEIGHT(I)  ! Intensity
	LSPAS = LSPAS + 4
D	WRITE(LUNP,830) I,NL,THETA(I),WAVE(I),HEIGHT(I)
D830	FORMAT(' I=',I3,'  NL=',I3,'  2TH=',F6.2,'  WAVE=',F6.2,
D    &   '  HEIGHT=',F6.2)
C
	NL=NL+1
290	CONTINUE
	CLOSE(UNIT=LSPEC)
300	OPEN(UNIT=LNIH,NAME='POWD:[PUBLIC.BASE.XRF]LAMBDA.NIH',
     %  ACCESS='DIRECT',
     %	RECORDSIZE=301,TYPE='OLD',READONLY,SHARED)
C  OPEN SCRATCH FILE
	OPEN(UNIT=LSCR,NAME='QUAL.TMP',FORM='UNFORMATTED',
     %  TYPE='SCRATCH',ACCESS='SEQUENTIAL')
C
C  CALCULATE ERROR IN WAVELENGTH (ASSUME COS(THETA)=1.0)
C
	ERR=0.008727*TWOD*D2T
C
	CALL MINMAX(WAVE,WMIN,WMAX,NL)
	WMIN=WMIN-ERR
	WMAX=WMAX+ERR
C
C--GET ORDER INFORMATION
C
	WRITE(LUNP,9400)
9400	FORMAT('$ MAXIMUM ORDER TO BE CONSIDERED [ D:2  MAX:9]>')
	READ(LUNI,9410,END=999)MAXORD
9410	FORMAT(I1)
	IF(MAXORD.EQ.0)MAXORD=2
	WRITE(LUNP,9420)
9420	FORMAT('$ ENTER KV ON X-RAY TUBE [ D:40. ] >')
	READ(LUNI,9423,END=999)AKV
9423	FORMAT(F5.0)
	IF(AKV.EQ.0)AKV=40.
	NORD=1
	M=0
500	CALL FIND(WAVE,MATCH,ERR,WMIN,WMAX,NL,LNIH,LSCR,NORD,AKV)
	WRITE(LUNO,9500)MATCH,NORD
9500    FORMAT('  ** ',I3,' POSSIBLE ELEMENTS IDENTIFIED - ORDER: ',I2)
	NORD=NORD+1
	M=M+MATCH
	IF(NORD.GT.MAXORD)GO TO 520
C
C  ORDERS 2 AND 6 NOT ALLOWED FOR GERMANIUM CRYSTAL
C
	IF(TWOD.NE.6.532)GO TO 500
	IF(NORD.NE.2.AND.NORD.NE.6)GO TO 500
	WRITE(LUNO,9510)NORD
9510	FORMAT('  ** ORDER ',I3,' NOT ALLOWED FOR GE CRYSTAL')
	NORD=NORD+1
C
	IF(NORD.LE.MAXORD)GO TO 500
C
C---WRITE RESULTS OF MATCH -- MAX OF 12 ELEMENTS PER MEASURED LINE
C
520	REWIND LSCR
	CALL OUTPT(THETA,WAVE,HEIGHT,NL,LSCR,LRES,M,ERR,LUNO,LOUT)
	CLOSE(UNIT=LRES,DISP='PRINT')
999	STOP
	END
	SUBROUTINE BSORT(W,INT,LINE,NLINE)
C
C  BUBBLE SORT BASED UPON INTENSITY
C
	REAL*8 LINE(75),TTEMP
	DIMENSION W(75)
	INTEGER INT(75)
	N=NLINE-1
100	NSWAP=0
	DO 300 K=1,N
	IF(INT(K).GE.INT(K+1))GO TO 300
	ITEMP=INT(K+1)
	INT(K+1)=INT(K)
	INT(K)=ITEMP
	TEMP=W(K+1)
	W(K+1)=W(K)
	W(K)=TEMP
	TTEMP=LINE(K+1)
	LINE(K+1)=LINE(K)
	LINE(K)=TTEMP
	NSWAP=1
300	CONTINUE
	IF(NSWAP.NE.0)GO TO 100
	RETURN
	END
	SUBROUTINE FIND(WAVE,MATCH,ERR,WMIN,WMAX,NL,LNIH,LSCR,NORD,AKV)
C
C---ROUTINE MATCHES STRONGEST LINE OF AN ELEMENT WITHIN THE MEASURED
C      RANGE TO THE MEASURED LINES
C   NO ATTEMPT IS MADE TO LOOK FOR OTHER LINES OF THAT ELEMENT
C   MATCHED ELEMENT INFO IS WRITTEN TO SCRATCH FILE FOR LATER RETRIEVAL
C---
	REAL*8 LINE(75)
	DIMENSION WAVE(100),W(75),INT(75)
C
	NREC=1
	MATCH=0
C
C  CALCULATE CUTOFF WAVELENGTH FOR ACCELERATING VOLTAGE USED
C
	CUTOFF=12.396/AKV
C
C  CALUCLUATE WAVELENGTH RANGE TO BE CONSIDERED FOR THIS ORDER
C
	WLOW=WMIN/NORD
	IF(CUTOFF.GT.WLOW)WLOW=CUTOFF
	WHIGH=WMAX/NORD
	IF(WHIGH.LT.WLOW)RETURN
C
C  START MATCH
C
100	READ(LNIH'NREC,ERR=900)IZ,W,INT,LINE
	CALL LIMIT(W,INT,LINE,WLOW,WHIGH,NLINE)
	IF(NLINE.EQ.0)GO TO 150
	CALL BSORT(W,INT,LINE,NLINE)
	STRONG=W(1)
	STRONG=STRONG*NORD
	DO 120 N=1,NL
	DWAVE=ABS(WAVE(N)-STRONG)
	IF(DWAVE.LE.ERR)GO TO 350
120	CONTINUE
	GO TO 150
350	MATCH=MATCH+1
	WRITE(LSCR)IZ,W,INT,LINE,NLINE,NORD
150	NREC=NREC+1
	IF(NREC.LE.91)GO TO 100
900	RETURN
	END
	SUBROUTINE GET2D(LUNI,LUNP,TWOD)
C  Modified 1-FEB for Alfred Univ.
100	WRITE(LUNP,9000)
9000	FORMAT(1X,' ENTER CRYSTAL: 1 - LiF 200'/17X,'2 - RAP 100'/
     %	 17X,'3 - PET'/17X,'4 - EDDT 020'/17X,'5 - Ge'/17X,'6 - ADP'
     %	 /17X,'7 - Other')   ! JDW
	WRITE(LUNP,9010)
9010	FORMAT('$ >')
	READ(LUNI,*,END=999)ICRY
	IF(ICRY.LT.1.OR.ICRY.GT.7)GO TO 100
	IF(ICRY.EQ.1)TWOD=4.027
	IF(ICRY.EQ.2)TWOD= 26.121    ! JDW
	IF(ICRY.EQ.3)TWOD=8.742
	IF(ICRY.EQ.4)TWOD= 8.808     ! JDW
	IF(ICRY.EQ.5)TWOD=6.532
	IF(ICRY.EQ.6)TWOD=10.648
	IF(ICRY.NE.7)GO TO 800
	WRITE(LUNO,9200)
9200	FORMAT('$ ENTER 2D FOR CRYSTAL>')
	READ(LUNI,*,END=999)TWOD
800	RETURN
999	STOP
	END
	SUBROUTINE LIMIT(W,INT,LINE,WMIN,WMAX,NLINE)
C
C  ROUTINE DETERMINES ELEMENTS OF W THAT ARE WITHIN WMIN,WMAX RANGE
C    AND RETURNS WITH THESE LINES IN THE FIRST NLINE ELEMENTS OF W
C
	REAL*8 LINE(75),L(75)
	DIMENSION W(75),W2(75)
	INTEGER INT(75),I2(75)
	N=0
	DO 200 K=1,75
	IF(W(K).LT.WMIN)GO TO 200
	IF(W(K).GT.WMAX)GO TO 200
	N=N+1
	W2(N)=W(K)
	I2(N)=INT(K)
	L(N)=LINE(K)
200	CONTINUE
	NLINE=N
	DO 300 K=1,NLINE
	W(K)=W2(K)
	INT(K)=I2(K)
	LINE(K)=L(K)
300	CONTINUE
	RETURN
	END
	SUBROUTINE MINMAX(WAVE,WMIN,WMAX,NL)
C*** DETERMINE THE MINIMUM AND MAXIMUM WAVELENGTH OF THE MEASURED PATTERN
	DIMENSION WAVE(1)
	WMIN=WAVE(1)
	WMAX=WAVE(1)
	DO 100 I=1,NL
	IF(WAVE(I).LT.WMIN)WMIN=WAVE(I)
	IF(WAVE(I).GT.WMAX)WMAX=WAVE(I)
100	CONTINUE
	RETURN
	END
	SUBROUTINE OUTPT(THETA,WAVE,HEIGHT,NL,LSCR,LRES,MATCH,ERR
     %	,LUNO,LOUT)
	REAL*8 LINE(75),OUTLIN(12)
	DIMENSION THETA(100),WAVE(100),HEIGHT(100),W(75)
	INTEGER INT(75),OUTSYM(12)
	BYTE ORDER(12)
crsx	OPEN(UNIT=LOUT,ACCESS='DIRECT',TYPE='SCRATCH',RECORDSIZE=34)
	OPEN(UNIT=LOUT,ACCESS='DIRECT',TYPE='SCRATCH',RECORDSIZE=156)
	DO 100 I=1,12
	ORDER(I)=0
	OUTSYM(I)=0
	OUTLIN(I)=0
100	CONTINUE
	DO 120 I=1,NL
	II=I
120	WRITE(LOUT'II)ORDER,OUTSYM,OUTLIN
	DO 300 I=1,MATCH
	READ(LSCR)IZ,W,INT,LINE,NLINE,NORD
	ISYM='  '
	CALL ZSYM(IZ,ISYM,LUNO)
	DO 200 K=1,NLINE
	W1=W(K)*NORD
	DO 150 L=1,NL
	LPT=L
	DELW=ABS(W1-WAVE(L))
	IF(DELW.LE.ERR)GO TO 170
150	CONTINUE
	GO TO 200
170	READ(LOUT'LPT)ORDER,OUTSYM,OUTLIN
	DO 180 L=1,12
	L2=L
	IF(OUTLIN(L).EQ.0)GO TO 190
180	CONTINUE
	GO TO 200
190	IF(ISYM.EQ.OUTSYM(L2-1))GO TO 200
	OUTLIN(L2)=LINE(K)
	ORDER(L2)=NORD+"60
	IF(NORD.EQ.1)ORDER(L2)=' '
	OUTSYM(L2)=ISYM
	WRITE(LOUT'LPT)ORDER,OUTSYM,OUTLIN
200	CONTINUE
300	CONTINUE
	DO 500 K=1,NL
	KK=K
	READ(LOUT'KK)ORDER,OUTSYM,OUTLIN
C
	DO 350, J=12,1,-1
	N=J
	IF (ORDER(J).NE.0) GOTO 360
350	CONTINUE
C
360	WRITE(LRES,9500)K,THETA(K),WAVE(K),HEIGHT(K),(ORDER(J),
     %	 OUTSYM(J),OUTLIN(J),J=1,N)
9500	FORMAT('0',I3,1X,F8.3,1X,F8.3,1X,F8.0,3X,4(A1,A2,1X,A8)/
     %	 34X,4(A1,A2,1X,A8)/34X,4(A1,A2,1X,A8))
500	CONTINUE
	CLOSE(UNIT=LOUT)
	RETURN
	END
      SUBROUTINE ZSYM(IZ,ELEM,LUTI)
C
C   PROGRAM TO CONVERT ELEMENT SYMBOL TO ATOMIC NUMBER IF ELEM IS BLANK
C   OR CONVERT ATOMIC NUMBER IF IZ IS ZERO
C
      INTEGER*2 IZ,ELEM,D(104)
C
      IF(IZ.EQ.0.AND.ELEM.EQ.'  ') GO TO 999
      IF(IZ.GT.104) GO TO 990
      D(1)='H '
      D(2)='HE'
      D(3)='LI'
      D(4)='BE'
      D(5)='B '
      D(6)='C '
      D(7)='N '
      D(8)='O '
      D(9)='F '
      D(10)='NE'
      D(11)='NA'
      D(12)='MG'
      D(13)='AL'
      D(14)='SI'
      D(15)='P '
      D(16)='S '
      D(17)='CL'
      D(18)='AR'
      D(19)='K '
      D(20)='CA'
      D(21)='SC'
      D(22)='TI'
      D(23)='V '
      D(24)='CR'
      D(25)='MN'
      D(26)='FE'
      D(27)='CO'
      D(28)='NI'
      D(29)='CU'
      D(30)='ZN'
      D(31)='GA'
      D(32)='GE'
      D(33)='AS'
      D(34)='SE'
      D(35)='BR'
      D(36)='KR'
      D(37)='RB'
      D(38)='SR'
      D(39)='Y '
      D(40)='ZR'
      D(41)='NB'
      D(42)='MO'
      D(43)='TC'
      D(44)='RU'
      D(45)='RH'
      D(46)='PD'
      D(47)='AG'
      D(48)='CD'
      D(49)='IN'
      D(50)='SN'
      D(51)='SB'
      D(52)='TE'
      D(53)='I '
      D(54)='XE'
      D(55)='CS'
      D(56)='BA'
      D(57)='LA'
      D(58)='CE'
      D(59)='PR'
      D(60)='ND'
      D(61)='PM'
      D(62)='SM'
      D(63)='EU'
      D(64)='GD'
      D(65)='TB'
      D(66)='DY'
      D(67)='HO'
      D(68)='ER'
      D(69)='TM'
      D(70)='YB'
      D(71)='LU'
      D(72)='HF'
      D(73)='TA'
      D(74)='W '
      D(75)='RE'
      D(76)='OS'
      D(77)='IR'
      D(78)='PT'
      D(79)='AU'
      D(80)='HG'
      D(81)='TL'
      D(82)='PB'
      D(83)='BI'
      D(84)='PO'
      D(85)='AT'
      D(86)='RN'
      D(87)='FR'
      D(88)='RA'
      D(89)='AC'
      D(90)='TH'
      D(91)='PA'
      D(92)='U '
      D(93)='NP'
      D(94)='PU'
      D(95)='AM'
      D(96)='CM'
      D(97)='BK'
      D(98)='CF'
      D(99)='ES'
      D(100)='FM'
      D(101)='MD'
      D(102)='NO'
      D(103)='LW'
      D(104)='LN'  ! ALL RARE EARTHS ALA JCPDS
C
      IF(ELEM.NE.'  ') GO TO 100
      ELEM=D(IZ)
      RETURN
  100 DO 200 I=1,104
  200 IF(ELEM.EQ.D(I)) IZ=I
      IF(ELEM.EQ.'D ')IZ=1
      RETURN
  999 WRITE(LUTI,300)
  300 FORMAT(1X,'WHICH ONE DO YOU WANT Z OR ELEM SHOULD I GUESS ?')
      STOP
  990 WRITE(LUTI,400)
  400 FORMAT(1X,'ATOMIC NUMBER OUT OF RANGE')
      STOP
      END
