C
C  PROGRAM COBRAG TO CARRY OUT A COPLAND-BRAGG ANALYSIS FOR
C  INTERNAL STANDARD METHOD QUANTATATIVE ANALYSIS
C
C  THE BINARY INPUT FILE ON LUN 30 MUST BE PREVIOUSLY PREPARED
C  AND WRITTEN BY PROGRAM QUANT
C
C          WRITTEN BY
C
C  R. L. SNYDER   NYS COLLEGE OF CERAMICS
C                 ALFRED UNIVERSITY
C                 ALFRED, NEW YORK
C
C and
C
C  C.R. HUBBARD   NATIONAL BUREAU OF STANDARDS
C                 A221 MATL
C                 WASHINGTON D.C. 20234
C
C
C  DIMENSIONAL INFORMATION:
C  THE MAXIMUM NUMBER OF UNKNOWN PHASES = 10 = N
C  THE NUMBER OF COLUMNS IN THE DETERMINATE IS N+1 OR 11
C  THE MAXIMUM NUMBER OF ROWS IN THE MATRIX = NDW = 100
C  THE DIMENSIONS ON THE WORKING STORAGE ARRAYS WS AND IP
C  ARE LARGE OVERESTIMATES OF THE ACTUAL SPACE REQUIRED
C  FOR MORE INFORMATION SEE THE COMMENTS IN LSEI.
C
	PARAMETER (NDW=100, N=10)
	COMMON/IO/ IRD,ITYP,IPRT,ICOBRG,LIMPRT
	DIMENSION A(NDW,N+1), PRGOPT(20), X(N), WS(1000), IP(500), SI(N)
	DIMENSION V(NDW)
	CHARACTER*100 TITLE
	CHARACTER*20 NAM
	CHARACTER*12  LABEL(NDW,4)
	CHARACTER*6 PHASE(N)
	CHARACTER*8 DAT,TIM
	CHARACTER*1 NOYES
       CHARACTER FILNAM*40
C
C-----LOGICAL UNIT# ASSIGNMENTS
C     IRD IS FOR CRT TERMINAL INPUT OR CARD READER INPUT
C     ITYP IS FOR CRT TERMINAL OUTPUT (REQUESTING USER RESPONSES)
C     IPRT IS FOR 132 COLUMN PRINTER OUTPUT (FINAL RESULTS)
C     ICOBRG IS FOR COPLAND-BRAGG BINARY FILE (SEE QUANT AND COBRAG)
C
	IRD  = 5
	ITYP = 5
	IPRT = 5
	CALL DATIM(DAT,TIM)
C
C  SAY HELLO
C
	WRITE(ITYP,1000) DAT,TIM
1000	FORMAT(/21X,'NBS*QUANT84-COBRAG',5X,'Revision Date: 29 FEB ',
     & '1984'/
     * /10X,'Copland-Bragg internal standard analysis allowing for '
     * 'overlapped lines'/
     * /25X,'Written BY R.L. Snyder and C.R. Hubbard '/
     * /29X,'Run on ',A8,5X,' AT ',A8//)
C
C  SELECT OUTPUT DESTINATION
C
        WRITE(ITYP,1)
1	FORMAT(/' Select Output Mode Desired:  (Output is 132 characte',
     1'rs wide)'/'       1 	- causes output to be put on the inter',
     2'active terminal.'/'       newname 	- causes output to go ',
     3'to the named file.'/'       <CR> 	- causes output to go ',
     4'to a file named COBRA.LIS.'/' Enter the output destination: '$)
       READ(IRD,2) FILNAM
2	FORMAT(A40)
       IF(FILNAM(1:1) .EQ. '1') GO TO 3
       IPRT = 6
       IF(FILNAM(1:1) .EQ. ' ') THEN
		OPEN(UNIT=IPRT,FILE='COBRAG.LIS',STATUS='NEW',
     1 ACCESS='SEQUENTIAL',FORM='FORMATTED')
       		GO TO 3
       ENDIF
	OPEN(UNIT=IPRT,FILE=FILNAM,STATUS='NEW',ACCESS='SEQUENTIAL',
     1 FORM='FORMATTED')
3	ICOBRG= 30
	IF(IPRT.NE.ITYP) WRITE(IPRT,1000)
C-----LIMPRT CONTROLS AMOUNT OF PRINT OUT:
C           =0 FOR FULL LISTING (MATRIX, SOLUTION, FINAL COMPARISON)
C           =1 FOR NORMAL LISTING (SOLUTION AND FINAL COMPARISON)
C           =2 FOR BRIEFEST LISTING (SOLUTION ONLY)
	WRITE(ITYP,11)
11	FORMAT(/' Select a print control option (I1):'/
     & ' 0 = Full listing (input matrix, solution and final comparison'/
     & ' 1 = Normal listing (solution and final comparison'/
     & ' 2 = Brief listing (solution only)')
	READ(IRD,12) LIMPRT
12	FORMAT(I1)
	IF(LIMPRT.LT.0 .OR. LIMPRT.GT.2) LIMPRT=0
C
C  ATTACH THE BINARY INPUT FILE FROM QUANT TO LOGICAL UNIT=ICOBRG
C
	WRITE(ITYP,20)
20	FORMAT(/' Input the name of the Copland Bragg file output by',
     1' QUANT.'/' (default extension = .COP): '$)
	READ(IRD,1007)NAM
1007	FORMAT(A20)
       DO 22 I=1,20
       IF(NAM(I:I) .EQ. ' ')GO TO 24
22	CONTINUE
C  SAFETY EXTENSION ASSUMES 6 LETTERS
       I = 7
24	NAM(I:I+4) = '.COP '
C
C  READ THE HEADER INFORMATION FROM THE BINARY QUANT FILE
C
	OPEN(UNIT=ICOBRG,ERR=50,FILE=NAM,STATUS='OLD',ACCESS='SEQUENTIAL',
     1 FORM='UNFORMATTED')
	GO TO 100
C
C  ERROR ON OPENING THE INPUT BINARY FILE FROM QUANT
C
50	WRITE(ITYP,1010)NAM
1010	FORMAT(/' Error on opening file ',A20,' CALL EXIT')
	GO TO 999
100	READ(ICOBRG)TITLE
	READ(ICOBRG) NUNK,(PHASE(I),I=1,NUNK)
	WRITE(IPRT,1020)NAM,TITLE
1020	FORMAT(/' Analysis of file: ',A20/' ',A100/)
	WRITE(IPRT,1025) NUNK,(PHASE(I),I=1,NUNK)
1025	FORMAT(/' The number of unknown phases = ',I5/
     * /' Phase Names: ',10(1X,A6))
C-----TEST IF ARRAY LIMITS WERE EXCEEDED
	IF(NUNK.LE.N) GOTO 180
	WRITE(IPRT,1026) N
1026	FORMAT(/' ***Number of phases exceeded array limit N=',I3,
     * ' ***CALL EXIT***')
	GOTO 999
C
C-----READ IN THE GENERALIZED COPELAND-BRAGG MATRIX  (A(I,J),J=1,NUNK)
C                 OBSERVED INTENSITY RATIO            A(I,NUNK+1)
C                 VARIANCE OF INTENSITY RATIO         V(I)
C                 PHASE NAME AND HKL'S               (LABEL(I,J),J=1,4)
C
180	DO 200 I=1,NDW
	READ(ICOBRG,END=250)(A(I,J),J=1,NUNK+1),V(I),(LABEL(I,J),J=1,4)
200	CONTINUE
C
C-----TOO MANY EQUATIONS INPUT FOR PROGRAM DIMENSION NDW
C
	WRITE(IPRT,1030) NDW
1030	FORMAT(/' ***Analysis limited to NDW=',I4,' lines*** ')
C
C-----CONTINUE SETTING UP FOR LEAST-SQUARES REFINEMENT
C
250	IROW = I - 1
	MA = IROW
C
C  SET THE LSEI OPTION TO COMPUTE THE COVARIENCE MATRIX
C
	PRGOPT(1) = 4
	PRGOPT(2) = 1
	PRGOPT(3) = 1
	PRGOPT(4) = 1
C
C  SHOULD WE APPLY BOUNDRY CONSTRAINTS ON THE SUM OF THE WT. FRACTIONS?
C
	WRITE(ITYP,1040)
1040	FORMAT(/' Enter the following number:'/
     1'  0.0 - Do not apply constraints to the sum of the wt. fractions'
     2/'  #.# - Constrain solution so that the sum of the wt. fractions'
     3/'        equal the input number (usually 1.0)')
	READ(IRD,1050)XSUM
1050	FORMAT(G15.7)
	ME = 0.0
	IF(XSUM.EQ.0.0)GO TO 280
C
C  IMPOSE THE CONSTRAINT ON THE SUM
C  FIRST PUSH ALL THE MATRIX ELEMENTS DOWN BY ONE TO MAKE ROOM)
C
	ME = 1
	IROW = IROW + 1
	DO 260 I=IROW,2,-1
	DO 260 J=1,NUNK+1
260	A(I,J) = A(I-1,J)
C
C  IMPOSE THE CONSTRAINT BY MAKING THE FIRST ROW ALL 1'S
C
	DO 270 J=1,NUNK
270	A(1,J) = 1.0
	A(1,NUNK+1) = XSUM
	WRITE(IPRT,272) XSUM
272	FORMAT(/' Solution to be constrained such that sum of wt. '
     * 'fractions = ',F6.3)
C
C  SHOULD INEQUALITY BOUNDRY CONSTRAINTS BE APPLIED?
C
280	WRITE(ITYP,1060)
1060	FORMAT(/' Do you wish to apply the following constraints?'/
     1'  Each X(I) .GE. 0.0  and Each X(I) .LE. 1.0'/
     2'  (answer yes or no)')
	READ(IRD,1070)NOYES
1070	FORMAT(A1)
	MG = 0
	IF(NOYES.EQ.'N')GO TO 295
C
C  IMPOSE THE INEQUALITY CONSTRAINTS
C
	MG = 2 * NUNK
	DO 290 I=1,NUNK
	IROW = IROW + 1
	DO 285 J=1,NUNK
	A(IROW,J) = 0.0
	A(IROW+1,J) = 0.0
285	CONTINUE
C  SET THE DIAGONAL ELEMENTS
	A(IROW,I) = 1.0
	A(IROW,NUNK+1) = 0.0
	IROW = IROW + 1
	A(IROW,I) = -1.0
	A(IROW,NUNK+1) = -1.0
290	CONTINUE
	WRITE(IPRT,292)
292	FORMAT(/' Solution to be constrained such that each wt fraction'
     * ' .GE.0.0  AND  .LE.1.0')
C
C  CALL THE LEAST SQUARES ROUTINE
C
295	IF(LIMPRT.NE.0)GO TO 298
C
C  IN DEBUG MODE LIST ENTIRE MATRIX
C
	WRITE(IPRT,1090)ME,MA,MG,NUNK,NDW
1090	FORMAT(/' Number of equalities = ',I5/
     1' Number of least squares equations = ',I5/
     2' Number of inequalities = ',I5/
     3' Number of unknown phases = ',I4/
     4' The ndw dimension = ',I5/)
	M = ME + MA + MG
	DO 296 I=1,M
296	WRITE(IPRT,1100)I,(A(I,J),J=1,NUNK+1)
1100	FORMAT(' Row:',I3,4X,10(1X,F10.6))
298	CALL LSEI(A,NDW,ME,MA,MG,NUNK,PRGOPT,X,R1,R2,MODE,WS,IP)
C
C  CHECK MODE FOR NORMAL COMPLETION
C
	IF(MODE.GT.4)MODE = 4
	MODE = MODE + 1
	GO TO(300,320,340,360,380),MODE
300	WRITE(IPRT,310)
310	FORMAT(/' Successful regression - all constraints both',
     1' equalities and inequalities were satisfied')
	GO TO 400
320	WRITE(IPRT,330)
330	FORMAT(/' Equality constraints are contradictory so',
     1' a generalized inverse solution was used')
	GO TO 400
340	WRITE(IPRT,350)
350	FORMAT(/' The inequality constraints are contradictory')
	GO TO 400
360	WRITE(IPRT,370)
370	FORMAT(/' Both the inequality and equality constraints',
     1' are contradictory')
	GO TO 400
380	WRITE(IPRT,390)
390	FORMAT(/' Insufficient work space in arrays - ',
     1'increase the size of parameter MDW')
	GO TO 999
C
C  OUTPUT THE RESULTS OF THE LEAST SQUARES
C
400	WRITE(IPRT,440)R2,R1
440	FORMAT(/' The residual error from least squares equations=',
     * F10.6/' ','The residual error from equality constraints   =',
     * F10.6//
     * '0',5X,'************************************'/
     * ' ',5X,'Phase    Original Sample   Least Sqs'/
     * ' ',5X,'Name     Weight Fraction     e.s.d. '/
     * ' ',5X,'************************************')
	DO 450 I=1,NUNK
	SIG = SQRT(A(I,I))
	WRITE(IPRT,410)PHASE(I),X(I),SIG
410	FORMAT(' ',5X,A6,F12.3,3X,F12.3)
450	CONTINUE
C
C-----OUTPUT THE CORRELATION MATRIX
C
	WRITE(IPRT,460) (PHASE(I),I=1,NUNK)
460	FORMAT(/' The full varience covarience matrix:'/
     * /7X, 10(6X,A6))
	DO 480 I=1,NUNK
480	SI(I) = SQRT(A(I,I))
	DO 500 I=1,NUNK
	DO 490 J=1,NUNK
490	A(I,J) = A(I,J)/(SI(I)*SI(J))
	WRITE(IPRT,510) PHASE(I),(A(I,J),J=1,NUNK)
510	FORMAT(' ',A6,10(3X,F9.4))
500	CONTINUE
C
C-----CALCULATE AND PRINT THE INTENSITY CALCULATED FROM THE
C-----WEIGHT FRACTIONS
C
	IF(LIMPRT.EQ.2) GOTO 999
	WRITE(IPRT,520)
520	FORMAT(/' Line',3(2X,'Phase   Line'),2X,'Standrd Line',4x,
     * 'Observed  Calculated    Obs-Calc      Sigma    Obs-Calc'/
     *       ' ','  # ',4(2X,'Name   h k l'),3X,
     * 'Int Ratio   Int Ratio    Int Ratio  Int Ratio     /Sigma')
C
C-----BEFOR THE CALC INTENSITIES CAN BE CALCULATED THE MATRIX ELEMENTS
C-----MUST BE REREAD. (THEY WERE REPLACED BY LSEI)
C
	REWIND ICOBRG
	READ(ICOBRG) TITLE
	READ(ICOBRG) NUNK,(PHASE(I),I=1,NUNK)
	DO 528 I = 1,MA
	READ(ICOBRG,END=530) (A(I,J),J=1,NUNK+1),V(I),(LABEL(I,J),J=1,4)
528	CONTINUE
C
C-----BEGIN CALCULATION OF INTENSITY LINE BY LINE
C
530	DO 600 I=1+ME,IROW-MG
	XNT = 0.0
	DO 540 J = 1,NUNK
540	XNT = XNT + A(I,J)*X(J)
	DIFF = A(I,NUNK+1) - XNT
	SIGMA = SQRT(V(I))
	PERCNT = DIFF/SIGMA
	WRITE(IPRT,560) I,(LABEL(I,K),K=1,4),A(I,NUNK+1),XNT,DIFF,SIGMA,
     & PERCNT
560	FORMAT(' ',I3,1X,4(2X,A12),4F12.3,F11.1)
600	CONTINUE
999	CLOSE(UNIT=ICOBRG,STATUS='KEEP')
	CALL EXIT
	END
