C ***  PROGRAM: INSCAL - INSTRUMENT CALIBRATION PROGRAM
C ***  VERSION: 880921  (YYMMDD)
C ***  AUTHOR:  S.A.HOWARD
C ***           DEPARTMENT OF CERAMIC ENGINEERING
C *             UNIVERISTY OF MISSOURI - ROLLA
C *
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1, WGFILE*32
      COMMON /WSGFID/ WGFILE
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      COMMON /PLTTYP/ KPTYPE
C *
C *   VAX 11/780 ERROR HANDLING ROUTINES
CVAX  EXTERNAL HDLR !THIS IS FOR VAX 780 ERRORS
CVAX  CALL LIB$ESTABLISH( HDLR )
CVAX  I = LIB$SIM_TRAP( SIGARGS, MECHARGS )
C *
      ITTI = 5
      ITTO = 6
      IDFI = 10
      ISCR = 11
      IOPL = 12
      DO 1 I=1,6
1     WGBAK(I) = 0.0
      DO 2 I=1,4
      DO 2 J=1,3
2     WGPRO(I,J) = 0.0
      DO 3 I=1,5
      WGRAD(I,1) = 0.0
      WGRAD(I,2) = 0.0
3     CONTINUE
      DO 4 I=1,3
4     WGCAL(I) = 0.0
      NWGPAR = 0
      TFWBAK = .FALSE.
      TFWPRO = .FALSE.
      TFWWAV = .FALSE.
      TFWCAL = .FALSE.
      KPTYPE = 0
C *
1000  WRITE(ITTO,1001)
1001  FORMAT(' INSCAL - Instrument calibration program  V 880921')
C *
C *  DETERMINE THE DISPOSITION OF THE CALIBRATION FILE
      CALL GETWSG
C *
C **  DETERMINE FUNCTION
C *
1100  WRITE(ITTO,1102)
1102  FORMAT(/,' ***  Choose the function to perform')
      WRITE(ITTO,1110)
1110  FORMAT(/,'   B - Background evaluation',/,
     1         '   C - Determine 2theta calibration curve',/,
     1         '   R - Radiation wavelength determination',/,
     1         '   P - Profile calibration',/,
     1         '   W - Write parameters to calibration data file',/,
     1         '   E - Exit program',//,
     1         ' Choice ? ')
      READ(ITTI,1120) QANS
1120  FORMAT(A1)
C *
C *  DO IT
      IF( QANS.EQ.'B' .OR. QANS.EQ.'b' ) GO TO 2000
      IF( QANS.EQ.'P' .OR. QANS.EQ.'p' ) GO TO 3000
      IF( QANS.EQ.'W' .OR. QANS.EQ.'w' ) GO TO 4000
      IF( QANS.EQ.'R' .OR. QANS.EQ.'r' ) GO TO 5000
      IF( QANS.EQ.'C' .OR. QANS.EQ.'c' ) GO TO 6000
      IF( QANS.EQ.'E' .OR. QANS.EQ.'e' ) GO TO 9000
      WRITE(ITTO,1130)
1130  FORMAT(/,' --Unrecognizable option:  Try again.')
      GO TO 1100
C *
C **  EVALUATE BACKGROUND SHAPE
C *
2000  CALL EVLBAK
      GO TO 1100
C *
C **  PERFORM PROFILE CALIBRATION
C *
3000  CALL EVLWSG
      GO TO 1100
C *
C **  WRITE PARAMETERS TO INSTRUMENT DATA FILE
C *
4000  CALL WRIWSG
      GO TO 1100
C *
C **  DETERMINE RADIATION WAVELENGTHS
C *
5000  CALL EVLRAD
      GO TO 1100
C *
C **  DETERMINE TWO-THETA CALIBRATION
C *
6000  CALL EVLCAL
      GO TO 1100
C *
C *  THATS IT
C *
9000  WRITE(ITTO,9010)
9010  FORMAT(/,' * Normal termination of INSCAL *',/)
CVAX
C      CALL EXIT
CVAX
CIBM
      STOP
CIBM
      END
C ***  GETWSG - READ CURRENT CALIBRATION FILE OR CREATE A NEW FILE
C *
      SUBROUTINE GETWSG
C *
      LOGICAL EXISTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1, QOPT*1, QLID*10
      CHARACTER QDATE*9, QTIME*8, LINE*72, WGFILE*32, FILEXX*20
      CHARACTER FILEID*20, BACKID*10, PROFID*10, WAVEID*10, CALBID*10
      PARAMETER ( FILEID = '*** INS CALIBRATION' )
      PARAMETER ( BACKID = 'BACKGROUND' )
      PARAMETER ( PROFID = 'PROFILES  ' )
      PARAMETER ( WAVEID = 'RADIATION ' )
      PARAMETER ( CALBID = 'DELTA 2THT' )
      COMMON /WSGFID/ WGFILE
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      DIMENSION WCU(5,2)
      EQUIVALENCE ( LINE(1:20), FILEXX )
      DATA WCU/ 1.5405981, 1.544435, 1.533169, 0.00, 0.00,
     1          1.00     , 0.50    , 0.008166, 0.00, 0.00/
C *
C **  DETERMINE FILE NAME AND WHETHER OR NOT IT CURRENTLY EXISTS
C *
2000  WRITE(ITTO,2010)
2010  FORMAT(/,' Instrument calibration data file name ? <WSGDAT> ')
      READ(ITTI,2020) WGFILE
2020  FORMAT(A32)
      IF( WGFILE.EQ.' ' .OR. WGFILE.EQ.'/' ) WGFILE = 'WSGDAT'
      IF( WGFILE(1:3).NE.'WSG' .AND. WGFILE(1:3).NE.'wsg' ) THEN
        WRITE(ITTO,2022)
2022    FORMAT(/,' --Instrument data file name must begin with "WSG":',
     1  ' try again')
        GO TO 2000
      ENDIF
C *
C *  TEST FOR FILE EXISTENCE
CVAX
      INQUIRE(FILE=WGFILE,EXIST=EXISTS)
      IF(EXISTS) GO TO 3000
CIBM  KLUDGE BECAUSE IBM INQUIRE STMT IS WIERD
C      WRITE(ITTO,2110)
C2110  FORMAT(' Does this calibration file currently exist ? <Y> ')
C      READ(ITTI,2203) QANS
C      IF( QANS.NE.'N' .AND. QANS.NE.'n' ) GO TO 3000
CIBM
C *
C *  INSTRUMENT CALIBRATION FILE DOES NOT EXIST... CREATE ONE
2200  WRITE(ITTO,2205)
2205  FORMAT(' -File does not exist: Creating new file')
CVAX
C      OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='NEW',FORM='FORMATTED',
C     1 CARRIAGECONTROL='LIST')
CVAX
CIBM
      OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='NEW',FORM='FORMATTED')
CIBM
      CALL DATE(QDATE)
      CALL TIME(QTIME)
      WRITE(IDFI,2210) FILEID,QDATE,QTIME
2210  FORMAT(1X,A20,27X,'FILED:',1X,A9,1X,A8)
      WRITE(ITTO,2220)
2220  FORMAT(/,' 72 character identifier for the file ? ')
      READ(ITTI,2230) LINE
2230  FORMAT(A72)
      WRITE(IDFI,2235) LINE
2235  FORMAT(1X,A72)
      ZERO = 0.0
      DO 2240 I=1,5
2240  WRITE(IDFI,2250) WAVEID,WCU(I,1),WCU(I,2),ZERO,QDATE,QTIME
      WRITE(IDFI,2250) BACKID,ZERO,ZERO,ZERO,QDATE,QTIME
      WRITE(IDFI,2250) BACKID,ZERO,ZERO,ZERO,QDATE,QTIME
2250  FORMAT(1X,A10,1P,3(1X,E13.6),1X,A9,1X,A8)
      DO 2260 I=1,4
2260  WRITE(IDFI,2250) PROFID,ZERO,ZERO,ZERO,QDATE,QTIME
      WRITE(IDFI,2250) CALBID,ZERO,ZERO,ZERO,QDATE,QTIME
      CLOSE(UNIT=IDFI)
      WRITE(ITTO,2270)
2270  FORMAT(/,' The default radiation wavelengths are those of Cu.',
     1 /,' Use option -R- to change them or perform W determination.')
      GO TO 9000
C *
C *  A FILE IS THERE... DETERMINE IF IT IS A CALIBRATION FILE
3000  OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='OLD',ERR=2200)
      READ(IDFI,2235,END=3020) LINE
      IF(FILEID.EQ.FILEXX) GO TO 3100
3020  WRITE(ITTO,3030)
3030  FORMAT(' --File is not an instrument calibration file')
      CLOSE(UNIT=IDFI,STATUS='KEEP')
C *
C *  AT THIS POINT WE SHOULD CHECK TO SEE IF THE FILE SHOULD BE OVERWRI
      STOP' CAN NOT OVERWRITE THIS FILE... YOU MAY WANT TO DELETE IT'
C *
C *  INSTRUMENT CALIBRATION FILE EXISTS... GET PARAMETERS FROM FILE
3100  WRITE(ITTO,3110)
3110  FORMAT(' -Instrument data file exists:  Reading parameters:',/)
      WRITE(ITTO,2235) LINE
      READ(IDFI,2235,END=3200) LINE
      WRITE(ITTO,2235) LINE
C *
C *  GET RADIATION PARAMETERS
      DO 3140 I=1,5
      READ(IDFI,3150,END=3200) QLID,WGRAD(I,1),WGRAD(I,2),Z1,QDATE,QTIME
      WRITE(ITTO,3150) QLID,WGRAD(I,1),WGRAD(I,2),Z1,QDATE,QTIME
3140  CONTINUE
3150  FORMAT(1X,A10,3(1X,E13.6),1X,A9,1X,A8)
C *
C *  GET THE BACKGROUND PARAMETERS
      READ(IDFI,3150,END=3200) QLID,(WGBAK(J),J=1,3),QDATE,QTIME
      WRITE(ITTO,3150) QLID,(WGBAK(J),J=1,3),QDATE,QTIME
      READ(IDFI,3150,END=3200) QLID,(WGBAK(J),J=4,6),QDATE,QTIME
      WRITE(ITTO,3150) QLID,(WGBAK(J),J=4,6),QDATE,QTIME
C *
C *  GET THE PROFILE PARAMETERS
      DO 3160 I=1,4
      READ(IDFI,3150,END=3200) QLID,(WGPRO(I,J),J=1,3),QDATE,QTIME
      WRITE(ITTO,3150) QLID,(WGPRO(I,J),J=1,3),QDATE,QTIME
3160  CONTINUE
C *
C *  GET THE TWO-THETA CALIBRATION PARAMETERS
      READ(IDFI,3150,END=3200) QLID,(WGCAL(J),J=1,3),QDATE,QTIME
      WRITE(ITTO,3150) QLID,(WGCAL(J),J=1,3),QDATE,QTIME
      CLOSE(UNIT=IDFI)
      GO TO 9000
C *
C *  A PREMATURE EOF HIT WAS ENCOUNTERED... USER MUST RECOVER
3200  CLOSE(UNIT=IDFI)
      WRITE(ITTO,3210)
3210  FORMAT(/,
     1' The end-of-file was encountered when it was not expected.',/,
     1' INSCAL can not recover from this error.  You should check',/,
     1' the file and determine if any of this information can be',/,
     1' saved.  Determine what caused corruption of the calibration',/,
     1' file.  Best advice is to use the existing calibration data',/,
     1' in file WSGPRO and go through the calibration process again.')
CVAX
C      CALL EXIT
CVAX
CIBM
      STOP
CIBM
C *
C *  CALIBRATION FILE READ OR CREATED
9000  RETURN
      END
C ***  EVLBAK - EVALUATE BACKGROUND SHAPE
C *
      SUBROUTINE EVLBAK
C *
      LOGICAL ASCII,LFINIS,LPLOTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QFILE*32, QANS*1, QTITLE*80, QTITLS(3)*20
      CHARACTER QDATE*9, QTIME*8
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /PLTTYP/ KPTYPE
      DIMENSION YDATA(10),KCHANL(3,18),BKFD(3:6),BKF(3:6)
      DIMENSION NPLSET(3),NCPLOT(3),OFFSET(3),NUMPEN(3),MARKER(3)
      DATA QTITLS(1)/'Raw data            '/
      DATA QTITLS(2)/'Points in background'/
      DATA QTITLS(3)/'Refined background  '/
      DATA BKFD(3)/     0.0000 /
      DATA BKFD(4)/  2500.0000 /
      DATA BKFD(5)/     3.0000 /
      DATA BKFD(6)/     2.0000 /
C *
C **  OPEN THE PATTERN FILE AND GET THE PATTERN INFORMATION
C *
C *  GET THE FILE NAME AND TYPE
1000  WRITE(ITTO,1002)
1002  FORMAT(/,' * Background evaluation')
      WRITE(ITTO,1010)
1010  FORMAT(/,' The name of the file containing background ? <EXIT> ')
      READ(ITTI,1020) QFILE
1020  FORMAT(A32)
      IF( QFILE(1:1).EQ.' ' .OR. QFILE(1:1).EQ.'/' ) GO TO 9000
      WRITE(ITTO,1040)
1040  FORMAT(' (F)ormatted or (U)nformatted file ? <F> ')
      READ(ITTI,1050) QANS
1050  FORMAT(A1)
      ASCII = QANS.NE.'U' .AND. QANS.NE.'u'
C *
C *  GET FILE PARAMETERS
1100  IF(ASCII) THEN
        OPEN(UNIT=IDFI,FILE=QFILE,STATUS='OLD',ERR=1110,
     1  FORM='FORMATTED')
        READ(IDFI,1120) QTITLE
1120    FORMAT(A80)
        READ(IDFI,1130) QDATE,QTIME
1130    FORMAT(A9,A8)
        READ(IDFI,1140) FBANG,FEANG,FAINC,FTIME
1140    FORMAT(4F10.4)
       ELSE
        OPEN(UNIT=IDFI,FILE=QFILE,STATUS='OLD',ERR=1110,
     1  FORM='UNFORMATTED')
        READ(IDFI) QTITLE
        READ(IDFI) QDATE,QTIME
        READ(IDFI) FBANG,FEANG,FAINC,FTIME
      ENDIF
      GO TO 1200
C *
C *  AN ERROR OCCURED ON FILE OPENING
1110  WRITE(ITTO,1112)
1112  FORMAT(/,' --The pattern file busy or not there.')
      GO TO 9000
C *
C *  CONFIRM FILE PARAMETERS
1200  WRITE(ITTO,1210) QTITLE(1:72),QDATE,QTIME,FBANG,FEANG,FAINC,FTIME
1210  FORMAT(/,' The run time parameters for this file are:',//,
     1 1X,A72,//,
     1 1X,'Date and time of pattern collection:',2X,A9,2X,A8,//,
     1 1X,'Starting angle   ',F8.4,3X,'Ending angle     ',F8.4,/,
     1 1X,'Angle increment  ',F8.4,3X,'Count time       ',F8.4)
C *
C *  DETERMINE NUMBER OF RECORDS AND POINTS PER RECORD TO READ
1300  NHID = NINT((FEANG-FBANG)/FAINC) + 1
      IF(NHID.GT.40960) THEN
        WRITE(ITTO,1310)
1310    FORMAT(/,' -Background limited to 4096 pts')
        NHID = 40960
      ENDIF
      NHIR = NHID / 10
      IF(NHIR.LT.1) STOP'--TOO FEW POINTS IN PATTERN'
      NHIP = 4096 / NHIR
      IF(NHIP.EQ.0) STOP'--TOO MANY POINTS IN PATTERN'
      IF(NHIP.GT.10) NHIP = 10
      NINTVL = 10 / NHIP
      IF(NINTVL.EQ.10) NINTVL = 1
      NTOGET = 1 + ( NHIP - 1 ) * NINTVL
C *
C *  READ PATTERN SEGMENTS FOR BACKGROUND DATA
1400  NPTS = 0
      TEMP1 = FBANG - FAINC
      TEMP2 = 10.0 * FAINC
      IF(ASCII) THEN
        DO 1410 I=1,NHIR
        READ(IDFI,1420,END=1490) (YDATA(J),J=1,NTOGET,NINTVL)
1420    FORMAT(10F7.0)
        ABREC = TEMP1 + TEMP2 * REAL(I-1)
        DO 1430 J=1,NTOGET,NINTVL
        NPTS    = NPTS + 1
        X(NPTS) = ABREC + REAL(J) * FAINC
        Y(NPTS) = YDATA(J)
1430    CONTINUE
1410    CONTINUE
       ELSE
        DO 1440 I=1,NHIR
        READ(IDFI,END=1490) (YDATA(J),J=1,NTOGET,NINTVL)
        ABREC = TEMP1 + TEMP2 * REAL(I-1)
        DO 1450 J=1,NTOGET,NINTVL
        NPTS    = NPTS + 1
        X(NPTS) = ABREC + REAL(J) * FAINC
        Y(NPTS) = YDATA(J)
1450    CONTINUE
1440    CONTINUE
      ENDIF
      GO TO 1500
1490  WRITE(ITTO,1492)
1492  FORMAT(/,
     1' -The end of the pattern file was prematurely hit.  Check the',
     1' scan angles in',/,
     1'  the file for correctness.  This program will use the',
     1' available data.')
      NHIR = I - 1
C *
C *   CONFIRM NUMBER OF POINTS READ FROM THE FILE
1500  CLOSE(UNIT=IDFI)
      WRITE(ITTO,1510) NHIR,NHIP,NPTS
1510  FORMAT(/,' Records read ',I5,'  Points/record ',I2,
     1'  Total number of points ',I5)
C *
C **  DETERMINE BACKGROUND MODEL TO USE
C *
C *  DETERMINE IF WSGBAK HAS NON-ZERO VALUES.  IF SO THEN EITHER A
C *  CALIBRATION FILE HAS BEEN READ OR BACKGROUND DETERMINATION HAS
C *  ALREADY BEEN PERFORMED.  GET HIGHEST NON-ZERO PARAMETER.
2000  NHIPRM = 0
      DO 2010 I=1,6
2010  IF(WGBAK(I).NE.0.0) NHIPRM = I
C *
C *  GET THE MODEL TO USE
      WRITE(ITTO,2020)
2020  FORMAT(/,' *  Background model',//,
     1' Model    Form',/,
     1'   1      BKG(2theta) = Offset',/,
     1'   2      BKG(2theta) = Slope*2theta + Offset',/,
     1'   3      BKG(2theta) = Pearson VII + Slope*2theta + Offset',//,
     1' Choice ? (1-3) <2> ')
      MODELB = 2
      READ(ITTI,*) MODELB
      IF( MODELB.LT.1 .OR. MODELB.GT.3 ) THEN
        MODELB = 2
        WRITE(ITTO,2030) MODELB
2030    FORMAT(' -Choice out of range: Value set to ',I1)
      ENDIF
      IF(MODELB.EQ.1) THEN
        NVECS = 1
        GO TO 3300
       ELSE IF(MODELB.EQ.2) THEN
        NVECS = 2
        GO TO 3300
       ENDIF
C *
C *  GET CHOICE OF PARAMTERS FOR PEARSON VII PROFILE
C *
2100  NVECS = 6
      IF(NHIPRM.EQ.NVECS) GO TO 2200
2105  WRITE(ITTO,2110)
2110  FORMAT(/,
     1' You have chosen to use a background model with a Pearson VII',
     1' profile included.',/,
     1' Normally the Pearson VII profile would be placed at the origin',
     1' to accomodate',/,
     1' the sharply increasing background due to scattering at low',
     1' angles.  You could',/,
     1' also use the profile to fit an amorphous profile if you should',
     1' so desire.',//,
     1' The default values will be approximately correct to account',
     1' for low angle',/,
     1' scattering.  Adjust the values to suit your needs.',//,
     1' Note the constraint on the angle during refinement is +/- 0.5',
     1' degrees while',/,
     1' the constraints on the other parameters are +/- 50% of their',
     1' values.')
      DO 2120 I=3,6
      BVECS(I) = BKFD(I)
      BKF(I) = 0.0
2120  CONTINUE
2130  WRITE(ITTO,2140) (BVECS(I),I=3,6)
2140  FORMAT(/,' Starting values for the Pearson VII profile:',//,
     1 ' Position   Intensity   FWHM   Exponent',/,
     1 1X,F8.4,3X,F9.2,2X,F5.2,4X,F5.2)
      WRITE(ITTO,2142)
2142  FORMAT(/,
     1 ' New non-zero values to replace above? <use above values> ')
      READ(ITTI,*) (BKF(I),I=3,6)
      NEW = 0
      DO 2150 I=3,6
      IF(BKF(I).EQ.0.0) GO TO 2150
      NEW = 1
      BVECS(I) = BKF(I)
      BMINS(I) = 0.5 * BKF(I)
      BMAXS(I) = 1.5 * BKF(I)
      BKF(I) = 0.0
2150  CONTINUE
      BMINS(3) = BVECS(3) - 0.5
      BMAXS(3) = BVECS(3) + 0.5
      IF(NEW.EQ.1) GO TO 2130
      GO TO 3300
C *
C *  INQUIRE AS TO WHETHER CURRENT VALUES SHOULD BE USED AS STARTING POINT
2200  WRITE(ITTO,2210)
2210  FORMAT(/,' Parameters values for the Pearson VII profile already',
     1' exist.  Either',/,
     1' they were read from the calibration file or you have already',
     1' performed background',/,
     1' refinement.')
      WRITE(ITTO,2140) (BVECS(I),I=3,6)
      WRITE(ITTO,2212)
2212  FORMAT(/,
     1' Do you want to use the current parameters for initial values?',
     1' <Y> ')
      READ(ITTI,1050) QANS
      IF( QANS.EQ.'N' .OR. QANS.EQ.'n' ) GO TO 2105
      DO 2220 I=3,6
      BVECS(I) = WGBAK(I)
      BMINS(I) = 0.5 * BVECS(I)
      BMAXS(I) = 1.5 * BVECS(I)
2220  CONTINUE
      BMINS(3) = BVECS(3) - 0.5
      BMAXS(3) = BVECS(3) + 0.5
C *
C **  GET PLOTTING PARAMETERS AND SET UP INITIAL INFORMATION
C *
3300  WRITE(ITTO,3310)
3310  FORMAT(/,' Plot file name ? (* = NONE) <BPLOT> ')
      READ(ITTI,3320) QFILE
3320  FORMAT(A32)
      LPLOTS = QFILE(1:1).NE.'*'
      IF(.NOT.LPLOTS) GO TO 4000
      IF ( QFILE(1:1).EQ.' ' .OR. QFILE(1:1).EQ.'/' ) QFILE = 'BPLOT'
C *
C *  GET A TITLE FOR THE DATA SET
      WRITE(ITTO,3352)
3352  FORMAT(/,' Enter a title for the data set ? ')
      READ(ITTI,3354) QTITLE
3354  FORMAT(A80)
C *
C *  OPEN THE PLOT FILE
      IF(KPTYPE.EQ.0) THEN
        OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='UNFORMATTED',
     1  ERR=3330)
       ELSE
CVAX
C         OPEN(UNIT=IOPL,FILE=QFILE,STATUS='NEW',FORM='FORMATTED',
C     1   CARRIAGECONTROL='LIST',ERR=3330)
CVAX
CIBM
         OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='FORMATTED',
     1   ERR=3330)
CIBM
      ENDIF
      GO TO 3350
3330  WRITE(ITTO,3340)
3340  FORMAT(' --Can not open specified plot file:  Continuing.')
      LPLOTS = .FALSE.
      GO TO 4000
C *
C *  SET DEFAULT PLOT PARAMETERS
C *  RAW DATA
3350  NPLSET(1) = 1
      NCPLOT(1) = 1
      OFFSET(1) = 0.0
      NUMPEN(1) = 1
      MARKER(1) = 1
C *  BACKGROUND DATA
      NPLSET(2) = 2
      NCPLOT(2) = 1
      OFFSET(2) = 0.0
      NUMPEN(2) = 2
      MARKER(2) = 1
C * REFINED BACKGROUND
      NPLSET(3) = 3
      NCPLOT(3) = 1
      OFFSET(3) = 0.0
      NUMPEN(3) = 3
      MARKER(3) = 1
C *
      NUMSTS = 3
      NPLOTS = 1
      IAXIS = 0
      KYID = 0
      DPI = 2.0
      FS = 200.0
C *
C *  WRITE THE INITIAL INFORMATION TO THE PLOT FILE
3400  IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) QTITLE
        WRITE(IOPL) FBANG, FEANG, FAINC
        WRITE(IOPL) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
        WRITE(IOPL) QTITLS(1),NCPLOT(1),OFFSET(1),NUMPEN(1),MARKER(1)
        WRITE(IOPL) QTITLS(2),NCPLOT(2),OFFSET(2),NUMPEN(2),MARKER(2)
        WRITE(IOPL) QTITLS(3),NCPLOT(3),OFFSET(3),NUMPEN(3),MARKER(3)
       ELSE
        WRITE(IOPL,3410) QTITLE
3410    FORMAT (A80)
        WRITE(IOPL,3420) FBANG,FEANG,FAINC
3420    FORMAT(3F10.4)
        WRITE(IOPL,3430) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
3430    FORMAT(3I9,2F9.2,I9)
        WRITE(IOPL,3440) QTITLS(1),NCPLOT(1),OFFSET(1),NUMPEN(1),
     1  MARKER(1)
3440    FORMAT(A20,I6,F6.2,2I6)
        WRITE(IOPL,3440) QTITLS(2),NCPLOT(2),OFFSET(2),NUMPEN(2),
     1  MARKER(2)
        WRITE(IOPL,3440) QTITLS(3),NCPLOT(3),OFFSET(3),NUMPEN(3),
     1  MARKER(3)
      ENDIF
C *
C *  WRITE THE RAW DATA TO THE PLOT FILE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(1),NPTS,FBANG
          WRITE(IOPL) (X(I),Y(I)/FTIME,I=1,NPTS)
         ELSE
          WRITE(IOPL,3450) NPLSET(1),NPTS,FBANG
3450      FORMAT(2I6,F10.4)
          WRITE(IOPL,3460) (X(I),Y(I)/FTIME,I=1,NPTS)
3460      FORMAT(2F10.4)
        ENDIF
      ENDIF
C *
C **  PERFORM THE OPTIMIZATION
C *
4000  WRITE(ITTO,4010)
4010  FORMAT(/,
     1' Several optimize/reject rounds will be performed with a limit',
     1' on the number of',/,
     1' cycles used during optimization.  After this the optimization',
     1' will be performed',/,
     1' with a larger limit on the number of cycles.  If you are a',
     1' coffee drinker now',/,
     1' would be a good time to get yourself a cup.')
      IR = 1
C *
C *  COMPUTE NECESSARY FACTORS
C *
      YAVG = 0.0
      DO 4012 I=1,NPTS
      YAVG   = YAVG + Y(I)
      WGT(I) = 1.0 / AMAX1( 1.0, Y(I) )
4012  CONTINUE
      YAVG = YAVG / REAL(NPTS)
C *
C *  SET CONSTRAINTS ON THE REFINEMENT PARAMETERS
      BVECS(1) = YAVG
      BMINS(1) = 0.0100
      BMAXS(1) = 2.0 * YAVG
      BVECS(2) = -0.005
      BMINS(2) = -5.000
      BMAXS(2) =  5.000
C *
C *  INITIAL COUNTS IN ANGULAR RANGE
4025  DO 4026 I=1,18
4026  KCHANL(3,I) = 0
      DO 4028 I=1,NPTS
      KCH = INT( X(I) ) / 10 + 1
      KCHANL(3,KCH) = KCHANL(3,KCH) + 1
4028  CONTINUE
C *
C *  USE THE GAUSS-NEWTON
      LFINIS = .FALSE.
4030  ITER = 5
4040  CALL BSOLVE(IR,ICON,RERR,ITER)
      IF(LFINIS) GO TO 6000
C *
C **  ELIMINATE POINTS TSIGMA SIGMA ABOVE COMPUTED LINE
C *
5000  TSIGMA = 2.0
      NDEL = 0
      CALL GNBACK( BVECS, YC )
      DO 5010 ICHK=1,NPTS
      YHI = YC(ICHK) + TSIGMA * SQRT( Y(ICHK) )
      IF(Y(ICHK).LT.YHI) GO TO 5010
      X(ICHK) = 0.0
      NDEL = NDEL + 1
5010  CONTINUE
      IF(NDEL.EQ.0) GO TO 5100
      NTEMP = NPTS
      NPTS = 0
      DO 5020 I=1,NTEMP
      IF(X(I).EQ.0.0) GO TO 5020
      NPTS = NPTS + 1
      X(NPTS)   = X(I)
      Y(NPTS)   = Y(I)
      YC(NPTS)  = YC(I)
      WGT(NPTS) = WGT(I)
5020  CONTINUE
C----------------------------------------TEMP---------------------------
C *
C * WRITE THE CALCULATED BACKGROUND TO THE PLOT FILE
C      IF(LPLOTS) THEN
C        IF(KPTYPE.EQ.0)THEN
C          WRITE(IOPL) NPLSET(3),NPTS,FBANG
C          WRITE(IOPL) (X(I),YC(I)/FTIME,I=1,NPTS)
C         ELSE
C          WRITE(IOPL,6050) NPLSET(3),NPTS,FBANG
C          WRITE(IOPL,3460) (X(I),YC(I)/FTIME,I=1,NPTS)
C        ENDIF
C      ENDIF
C----------------------------------------TEMP---------------------------
C *
5100  NCYCLE = NCYCLE + 1
      WRITE(ITTO,9992) NCYCLE,NDEL,NPTS
9992  FORMAT(/,' Cycle ',I2,'  PTS deleted ',I4,'  Remaining ',I4)
      DO 5103 I=1,18
5103  KCHANL(1,I) = 0
      DO 5105 I=1,NPTS
      KCH = INT( X(I) ) / 10 + 1
      KCHANL(1,KCH) = KCHANL(1,KCH) + 1
5105  CONTINUE
      DO 5106 I = 1,18
      KCHANL(2,I) = KCHANL(3,I) - KCHANL(1,I)
      KCHANL(3,I) = KCHANL(1,I)
5106  CONTINUE
      WRITE(ITTO,5107) ((KCHANL(I,J),J=1,9),I=1,2),
     1 ((KCHANL(I,J),J=10,18),I=1,2)
5107  FORMAT(/,
     1 ' ANGS 000-010 010-020 020-030 030-040 040-050 050-060 060-070',
     1 ' 070-080 080-090',/,
     1 ' REM  ',9(2X,I4,2X),/,
     1 ' DEL  ',9(2X,I4,2X),/,
     1 ' ANGS 090-100 100-110 110-120 120-130 130-140 140-150 150-160',
     1 ' 160-170 170-180',/,
     1 ' REM  ',9(2X,I4,2X),/,
     1 ' DEL  ',9(2X,I4,2X))
      IF(NDEL.GT.0) GO TO 4030
C *
C *  NOW WE NEED ACCURATE VALUES FOR THE PARAMETERS
      LFINIS = .TRUE.
      ITER = 20
      WRITE(ITTO,5110)
5110  FORMAT(/,' Once more through the refinement for best parameters.')
      GO TO 4040
C *
C *  WRITE THE RESULTS
6000  WRITE(ITTO,6010)
6010  FORMAT(/,' The refined parameters are as follows:',/)
      DO 6020 I=1,NVECS
      WGBAK(I) = BVECS(I)
      WRITE(ITTO,6030) I,WGBAK(I)
6030  FORMAT(' Parameter ',I1,': ',E13.6)
6020  CONTINUE
      TFWBAK = .TRUE.
      WRITE(ITTO,6040) RERR,NCYCLE,ICON
6040  FORMAT(/,' Residual error after refinement:   ',F10.6,/,
     1         ' Total number of refinement cycles: ',I10,/,
     1         ' Optimization convergence parameter:',I10)
C *
C *  WRITE THE REMAINING BACKGROUND POINTS TO THE PLOT FILE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(2),NPTS,FBANG
          WRITE(IOPL) (X(I),Y(I)/FTIME,I=1,NPTS)
        ELSE
          WRITE(IOPL,6050) NPLSET(2),NPTS,FBANG
6050      FORMAT(2I6,F10.4)
          WRITE(IOPL,3460) (X(I),Y(I)/FTIME,I=1,NPTS)
        ENDIF
      ENDIF
C *
C * WRITE THE CALCULATED BACKGROUND TO THE PLOT FILE
      IF(LPLOTS) THEN
        CALL GNBACK( BVECS, YC )
        IF(KPTYPE.EQ.0)THEN
          WRITE(IOPL) NPLSET(3),NPTS,FBANG
          WRITE(IOPL) (X(I),YC(I)/FTIME,I=1,NPTS)
         ELSE
          WRITE(IOPL,6050) NPLSET(3),NPTS,FBANG
          WRITE(IOPL,3460) (X(I),YC(I)/FTIME,I=1,NPTS)
        ENDIF
      ENDIF
C *
C *  BACKGROUND EVALUATED
9000  IF(LPLOTS) CLOSE(UNIT=IOPL)
      RETURN
      END
C ***  EVLRAD - EVALUATE RADIATION WAVELENGTHS
C *
      SUBROUTINE EVLRAD
C *
      LOGICAL USECOP,USELPC
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1
      PARAMETER ( DTORD2=8.72664626E-3, RTODT2=114.591559 )
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      DIMENSION WCU(3,2)
      DATA WCU/ 1.5405981, 1.544435, 1.533169,
     1          1.00     , 0.50    , 0.008166/
C *
C *  INITIALIZE VALUES
      DO 10 I=1,5
      WGRAD(I,1) = 0.0
      WGRAD(I,2) = 0.0
10    CONTINUE
      NWAVE = 0
C *
C *  DETERMINE FUNCTION
1000  WRITE(ITTO,1002)
1002  FORMAT(/,' * Radiation wavelength and relative intensity',
     1 ' determination')
1004  WRITE(ITTO,1010)
1010  FORMAT(/,'   S - Set Cu k-alpha as default wavelengths',/,
     1         '   P - Perform Wavelength determination',/,
     1         '   E - Exit this section',//,
     1         ' Choice ? ')
      READ(ITTI,1020) QANS
1020  FORMAT(A1)
C *
C *  DO IT
1200  IF( QANS.EQ.'S' .OR. QANS.EQ.'s' ) GO TO 2000
      IF( QANS.EQ.'P' .OR. QANS.EQ.'p' ) GO TO 3000
      IF( QANS.EQ.'E' .OR. QANS.EQ.'e' ) GO TO 9000
      GO TO 1004
C *
C *  GET THE INFORMATION REGARDING THE RADIATION WAVELENGTHS
2000  DO 2005 I=1,3
      WGRAD(I,1) = WCU(I,1)
      WGRAD(I,2) = WCU(I,2)
2005  CONTINUE
      NWAVE = 3
      WRITE(ITTO,2010)
2010  FORMAT(/,' The standard Cu radiation wavelengths are default.')
2020  WRITE(ITTO,2030)
2030  FORMAT(/,' Wavelength   Relative Intensity')
      DO 2040 I=1,NWAVE
2040  WRITE(ITTO,2050) WGRAD(I,1),WGRAD(I,2)
2050  FORMAT(2X,F8.5,9X,F6.4)
      WRITE(ITTO,2060)
2060  FORMAT(/,' Are these values acceptable ? <Y> ')
      READ(ITTI,1020) QANS
      IF( QANS.NE.'N' .AND. QANS.NE.'n' ) THEN
        TFWWAV = .TRUE.
        GO TO 9000
      ENDIF
C *
C *  GET THE NEW RADIATION WAVELENGTHS
2100  DO 2110 I=1,5
      WGRAD(I,1) = 0.0
      WGRAD(I,2) = 0.0
2110  CONTINUE
      NWAVE = 0
      WRITE(ITTO,2120)
2120  FORMAT(/,' Enter the new wavelengths in (A) and relative',
     1   ' intensities by pairs.',/,' A blank line ends input.')
      DO 2130 I=1,5
      A = 0.0
      B = 0.0
      READ(ITTI,*) A,B
      IF( A.LE.0.0 .OR. B.LE.0.0 ) GO TO 2020
      NWAVE = NWAVE + 1
      WGRAD(NWAVE,1) = A
      WGRAD(NWAVE,2) = B
2130  CONTINUE
      GO TO 2020
C *
C **  READ INSTRUMENT PROFILE FILE AND DETERMINE THE ALPHA-2 WAVELENGTH
C *    AND RELATIVE INTENSITY
C *
3000  NWGPAR = 0
      CALL GETPRO
      IF(NWGPAR.EQ.0) THEN
        WRITE(ITTO,3002)
3002    FORMAT(' --No lines to evaluate: Try again')
        GO TO 1000
      ENDIF
C *
C *  DETERMINE PRIMARY WAVELENGTH
      WRITE(ITTO,3010) WCU(1,1)
3010  FORMAT(/,' The default primary wavelength is Cu k-alpha 1 = ',
     1 F9.6,' A.',/,' New value ? <USE DEFAULT> ')
      A = 0.0
      READ(ITTI,*) A
      USECOP = .NOT. (A.GT.0.0)
      IF(USECOP) THEN
        WGRAD(1,1) = WCU(1,1)
        WGRAD(1,2) = WCU(1,2)
       ELSE
        WGRAD(1,1) = A
        WGRAD(1,2) = 1.0
      ENDIF
C *
C *  GET THE CALCULATED TWO-THETAS FOR THE LINE
      CALL ASS2TH( 1, NWGPAR, WGRAD(1,1) )
C *
C *  GET AN ESTIMATE OF THE TWO-THETA ERROR FOR THIS DIFFRACTION LINE
      DEL2TH = WSGPAR(1,12) - WSGPAR(1,1)
      WRITE(ITTO,1234) WSGPAR(1,1),DEL2TH
1234  FORMAT(/,' Apparent two-theta error at ',F8.4,' degrees is ',F6.4)
C *
C *  BASED ON THE TRUE D-SPACING, CALCULATE THE APPARENT WAVELENGTHS
      TTHT1C = WSGPAR(1,1) + DEL2TH
      A1 = 2.0 * WSGPAR(1,13) * SIN( TTHT1C * DTORD2 )
      A2 = 0.0
      A3 = 0.0
      R21  = 0.0
      R21P = 0.0
      R31  = 0.0
      R31P = 0.0
      IF(NWGPAR.GT.1) THEN
        TTHT2C = WSGPAR(2,1) + DEL2TH
        A2 = 2.0 * WSGPAR(1,13) * SIN( TTHT2C * DTORD2 )
        XX = WSGPAR(2,2) * PLORNZ(TTHT1C) / PLORNZ(TTHT2C)
        R21  = WSGPAR(2,2) / WSGPAR(1,2)
        R21P = XX          / WSGPAR(1,2)
      ENDIF
      IF(NWGPAR.GT.2) THEN
        TTHT3C = WSGPAR(3,1) + DEL2TH
        A3 = 2.0 * WSGPAR(1,13) * SIN( TTHT3C * DTORD2 )
        XX = WSGPAR(3,2) * PLORNZ(TTHT1C) / PLORNZ(TTHT3C)
        R31  = WSGPAR(3,2) / WSGPAR(1,2)
        R31P = XX          / WSGPAR(1,2)
      ENDIF
      WRITE(ITTO,3030) A1,A2,R21,R21P,A3,R31,R31P
3030  FORMAT(/,'  Radiation  Wavelength   Ratio  Ratio w/ L/P',/,
     1         '     A1      ',F10.6,'  1.0000    1.0000',/,
     1         '     A2      ',F10.6,'  ',F6.4,'    ',F6.4,/,
     1         '     A3      ',F10.6,'  ',F6.4,'    ',F6.4)
      WRITE(ITTO,3032)
3032  FORMAT(/,' Use L/P corrected intensities ? <Y> ')
      READ(ITTI,1020) QANS
      USELPC = QANS.NE.'N' .AND. QANS.NE.'n'
C *
C *  PUT THE NEW WAVELENGTHS IN THE DATA ARRAY
      WGRAD(1,1) = A1
      WGRAD(1,2) = 1.0
      IF(A2.GT.0.0) THEN
        WGRAD(2,1) = A2
        WGRAD(2,2) = R21
        IF(USELPC) WGRAD(2,2) = R21P
      ENDIF
      IF(A3.GT.0.0) THEN
        WGRAD(3,1) = A3
        WGRAD(3,2) = R31
        IF(USELPC) WGRAD(3,2) = R31P
      ENDIF
      IF(USECOP) THEN
        WRITE(ITTO,3040) WCU(3,1)
3040    FORMAT(/,' Cu k-alpha 3 wavelength = ',F9.6,' Use it ? <N> ')
        READ(ITTI,1020) QANS
        IF( QANS.NE.'Y' .AND. QANS.NE.'y' ) GO TO 3100
        WGRAD(3,1) = WCU(3,1)
        WGRAD(3,2) = WCU(3,2)
      ENDIF
3100  TFWWAV = .TRUE.
C *
9000  RETURN
      END
C ***  EVLWSG - EVALUATE PROFILES FROM STANDARDS FOR CALIBRATION
C *
      SUBROUTINE EVLWSG
C *
      LOGICAL LPLOTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1
      CHARACTER QDATE*9, QTIME*8, LINE*72
      CHARACTER QTITLS(8)*20, QFILE*32, QTITLE*80
      PARAMETER ( DTOR=1.745329252E-02, RTOD=57.29577951 )
      PARAMETER ( DTORD2=8.72664626E-3, RTODT2=114.591559 )
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      COMMON /PLTTYP/ KPTYPE
      DIMENSION NPLSET(8),NCPLOT(8),OFFSET(8),NUMPEN(8),MARKER(8)
      DATA QTITLS(1)/'FWHM front  obs     '/
      DATA QTITLS(2)/'FWHM front  calc    '/
      DATA QTITLS(3)/'EXP  front  obs     '/
      DATA QTITLS(4)/'EXP  front  calc    '/
      DATA QTITLS(5)/'FWHM back   obs     '/
      DATA QTITLS(6)/'FWHM back   calc    '/
      DATA QTITLS(7)/'EXP  back   obs     '/
      DATA QTITLS(8)/'EXP  back   calc    '/
C *
      NWGPAR = 0
C *
C *  DETERMINE FUNCTION
1000  WRITE(ITTO,1002)
1002  FORMAT(/,' * Profile shape calibration')
1004  WRITE(ITTO,1010)
1010  FORMAT(/,'   R - Read profile data file',/,
     1         '   P - Perform profile analysis',/,
     1         '   E - Exit this section',//,
     1         ' Choice ? ')
      READ(ITTI,1020) QANS
1020  FORMAT(A1)
C *
C *  DO IT
1200  IF( QANS.EQ.'R' .OR. QANS.EQ.'r' ) GO TO 2000
      IF( QANS.EQ.'P' .OR. QANS.EQ.'p' ) GO TO 3000
      IF( QANS.EQ.'E' .OR. QANS.EQ.'e' ) GO TO 9000
      GO TO 1004
C *
C **  READ PROFILE DESCRIPTIORS FROM PROFILE DATA FILE
C *
2000  CALL GETPRO
      GO TO 1000
C *
C **  DETERMINE THE CALIBRATION CURVES
C *
3000  NPTS = NWGPAR
      WRITE(ITTO,3010) NPTS
3010  FORMAT(/,' A total of ',I3,' lines are to be evaluated')
C *
C *  DATA CAN BE PLOTTED.  CHECK FOR PLOTTING BY NAME...
4000  WRITE(ITTO,4010)
4010  FORMAT(/,' Plot file name ? (* = NONE) <PPLOT> ')
      READ(ITTI,4020) QFILE
4020  FORMAT(A32)
      LPLOTS = QFILE(1:1).NE.'*'
      IF(.NOT.LPLOTS) GO TO 3100
      IF( QFILE(1:1).EQ.' ' .OR. QFILE(1:1).EQ.'/' ) QFILE = 'PPLOT'
C *
C *  GET A TITLE FOR THE DATA SET
      WRITE(ITTO,4052)
4052  FORMAT(/,' Enter a title for the data set ? ')
      READ(ITTI,4054) QTITLE
4054  FORMAT(A80)
C *
C *  OPEN THE PLOT FILE
      IF(KPTYPE.EQ.0) THEN
        OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='UNFORMATTED',
     1  ERR=4030)
       ELSE
CVAX
C     OPEN(UNIT=IOPL,FILE=QFILE,STATUS='NEW',FORM='FORMATTED',
C    1  CARRIAGECONTROL='LIST',ERR=4030)
CVAX
CIBM
      OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='FORMATTED',
     1  ERR=4030)
CIBM
      ENDIF
      GO TO 4050
4030  WRITE(ITTO,4040)
4040  FORMAT(' --Can not open specified plot file')
      LPLOTS = .FALSE.
      GO TO 3100
C *
C *  SET DEFAULT PARAMETERS
C *  FWHM FRONT
4050  NPLSET(1) =  1
      NCPLOT(1) = -1
      OFFSET(1) =  0.0
      NUMPEN(1) =  1
      MARKER(1) =  1
      NPLSET(2) =  2
      NCPLOT(2) =  1
      OFFSET(2) =  0.0
      NUMPEN(2) =  1
      MARKER(2) =  1
C *  EXP FRONT
      NPLSET(3) =  3
      NCPLOT(3) = -1
      OFFSET(3) =  0.0
      NUMPEN(3) =  2
      MARKER(3) =  1
      NPLSET(4) =  4
      NCPLOT(4) =  1
      OFFSET(4) =  0.0
      NUMPEN(4) =  2
      MARKER(4) =  1
C *  FWHM BACK
      NPLSET(5) =  5
      NCPLOT(5) = -1
      OFFSET(5) =  0.0
      NUMPEN(5) =  3
      MARKER(5) =  1
      NPLSET(6) =  6
      NCPLOT(6) =  1
      OFFSET(6) =  0.0
      NUMPEN(6) =  3
      MARKER(6) =  1
C *  EXP BACK
      NPLSET(7) =  7
      NCPLOT(7) = -1
      OFFSET(7) =  0.0
      NUMPEN(7) =  4
      MARKER(7) =  1
      NPLSET(8) =  8
      NCPLOT(8) =  1
      OFFSET(8) =  0.0
      NUMPEN(8) =  4
      MARKER(8) =  1
C *
      NUMSTS = 8
      NPLOTS = 1
      IAXIS = 0
      KYID = 1
      DPI = 2.0
      FS = 2.0
C *
C *  DETERMINE THE RANGE FOR THE PLOT FILE OUTPUT
      STRT = 180.0
      FINS =   0.0
      DO 4060 I=1,NPTS
      IF(STRT.GT.WSGPAR(I,1)) STRT = WSGPAR(I,1)
      IF(FINS.LT.WSGPAR(I,1)) FINS = WSGPAR(I,1)
4060  CONTINUE
      STRT = AINT( STRT/10.0 ) * 10.0
      FINS = AINT( FINS/10.0 ) * 10.0 + 10.0
      NPOUT = 100
      FINC = ( FINS - STRT ) / REAL( NPOUT - 1 )
C *
C *  WRITE THE INITIAL INFORMATION TO THE PLOT FILE
4100  IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) QTITLE
        WRITE(IOPL) STRT,FINS,FINC
        WRITE(IOPL) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
        WRITE(IOPL) QTITLS(1),NCPLOT(1),OFFSET(1),NUMPEN(1),MARKER(1)
        WRITE(IOPL) QTITLS(2),NCPLOT(2),OFFSET(2),NUMPEN(2),MARKER(2)
        WRITE(IOPL) QTITLS(3),NCPLOT(3),OFFSET(3),NUMPEN(3),MARKER(3)
        WRITE(IOPL) QTITLS(4),NCPLOT(4),OFFSET(4),NUMPEN(4),MARKER(4)
        WRITE(IOPL) QTITLS(5),NCPLOT(5),OFFSET(5),NUMPEN(5),MARKER(5)
        WRITE(IOPL) QTITLS(6),NCPLOT(6),OFFSET(6),NUMPEN(6),MARKER(6)
        WRITE(IOPL) QTITLS(7),NCPLOT(7),OFFSET(7),NUMPEN(7),MARKER(7)
        WRITE(IOPL) QTITLS(8),NCPLOT(8),OFFSET(8),NUMPEN(8),MARKER(8)
       ELSE
        WRITE(IOPL,4110) QTITLE
4110    FORMAT(A80)
        WRITE(IOPL,4120) STRT,FINS,FINC
4120    FORMAT(3F10.4)
        WRITE(IOPL,4130) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
4130    FORMAT(3I9,2F9.2,I9)
        WRITE(IOPL,4140) QTITLS(1),NCPLOT(1),OFFSET(1),NUMPEN(1),
     1  MARKER(1)
4140    FORMAT(A20,I6,F6.2,2I6)
        WRITE(IOPL,4140) QTITLS(2),NCPLOT(2),OFFSET(2),NUMPEN(2),
     1  MARKER(2)
        WRITE(IOPL,4140) QTITLS(3),NCPLOT(3),OFFSET(3),NUMPEN(3),
     1  MARKER(3)
        WRITE(IOPL,4140) QTITLS(4),NCPLOT(4),OFFSET(4),NUMPEN(4),
     1  MARKER(4)
        WRITE(IOPL,4140) QTITLS(5),NCPLOT(5),OFFSET(5),NUMPEN(5),
     1  MARKER(5)
        WRITE(IOPL,4140) QTITLS(6),NCPLOT(6),OFFSET(6),NUMPEN(6),
     1  MARKER(6)
        WRITE(IOPL,4140) QTITLS(7),NCPLOT(7),OFFSET(7),NUMPEN(7),
     1  MARKER(7)
        WRITE(IOPL,4140) QTITLS(8),NCPLOT(8),OFFSET(8),NUMPEN(8),
     1  MARKER(8)
      ENDIF
C *
C *  MOVE FRONT FWHM TO DATA ARRAYS
3100  DO 3110 I=1,NPTS
      X(I) = TAN( WSGPAR(I,1) * DTORD2 )
      Y(I) = ( WSGPAR(I,3)/2.0 )**2
3110  CONTINUE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(1),NPTS,0.0
          WRITE(IOPL) (WSGPAR(I,1),WSGPAR(I,3),I=1,NPTS)
         ELSE
          WRITE(IOPL,3112) NPLSET(1),NPTS,0.0
3112      FORMAT(2I6,F10.4)
          WRITE(IOPL,3114)  (WSGPAR(I,1),WSGPAR(I,3),I=1,NPTS)
3114      FORMAT(2F10.4)
        ENDIF
      ENDIF
C *
C *  USE A 2ND ORDER FIT
      RERR = PFIT( X, Y, W, NPTS, 2, PYINT, BVECS, 0 )
C *
C *  WRITE THE RESULTS
      WRITE(ITTO,3120)
3120  FORMAT(/,' The refined parameters are as follows:',/)
      DO 3130 I=1,3
      WGPRO(1,I) = BVECS(I)
      WRITE(ITTO,3140) I,WGPRO(1,I)
3140  FORMAT(' Parameter ',I2,': ',E13.6)
3130  CONTINUE
      WRITE(ITTO,3150) RERR
3150  FORMAT(/,' Standard deviation of fit:   ',F10.6)
      IF(.NOT.LPLOTS) GO TO 3200
      DO 3160 I=1,NPOUT
      X(I) = STRT + REAL(I-1) * FINC
      TANTHT = TAN( X(I) * DTORD2 )
      Y(I) = 2.0*SQRT(WGPRO(1,1)+WGPRO(1,2)*TANTHT+WGPRO(1,3)*TANTHT**2)
3160  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(2),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,3112) NPLSET(2),NPOUT,0.0
        WRITE(IOPL,3114) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C *
C *  MOVE FRONT EXPONENT TO DATA ARRAYS
3200  DO 3210 I=1,NPTS
      X(I) = WSGPAR(I,1)
      Y(I) = WSGPAR(I,4)
3210  CONTINUE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(3),NPTS,0.0
          WRITE(IOPL) (X(I),Y(I),I=1,NPTS)
         ELSE
          WRITE(IOPL,3112) NPLSET(3),NPTS,0.0
          WRITE(IOPL,3114) (X(I),Y(I),I=1,NPTS)
        ENDIF
      ENDIF
C *
C *  USE 2ND ORDER FIT
      RERR = PFIT( X, Y, W, NPTS, 2, PYINT, BVECS, 0 )
C *
C *  WRITE THE RESULTS
      WRITE(ITTO,3120)
      DO 3230 I=1,3
      WGPRO(2,I) = BVECS(I)
      WRITE(ITTO,3140) I,WGPRO(2,I)
3230  CONTINUE
      WRITE(ITTO,3150) RERR
      IF(.NOT.LPLOTS) GO TO 3300
      DO 3240 I=1,NPOUT
      X(I) = STRT + REAL(I-1) * FINC
      Y(I) = WGPRO(2,1) + WGPRO(2,2)*X(I) + WGPRO(2,3)*X(I)**2
3240  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(4),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,3112) NPLSET(4),NPOUT,0.0
        WRITE(IOPL,3114) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C *
C *  MOVE BACK FWHM TO DATA ARRAYS
3300  DO 3310 I=1,NPTS
      X(I) = TAN( WSGPAR(I,1) * DTORD2 )
      Y(I) = ( WSGPAR(I,5)/2.0 )**2
3310  CONTINUE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(5),NPTS,0.0
          WRITE(IOPL) (WSGPAR(I,1),WSGPAR(I,5),I=1,NPTS)
         ELSE
          WRITE(IOPL,3112) NPLSET(5),NPTS,0.0
          WRITE(IOPL,3114) (WSGPAR(I,1),WSGPAR(I,5),I=1,NPTS)
        ENDIF
      ENDIF
C *
C *  SET INITIAL VALUES
C *
C *  USE A 2ND ORDER FIT
      RERR = PFIT( X, Y, W, NPTS, 2, PYINT, BVECS, 0 )
C *
C *  WRITE THE RESULTS
      WRITE(ITTO,3120)
      DO 3330 I=1,3
      WGPRO(3,I) = BVECS(I)
      WRITE(ITTO,3140) I,WGPRO(3,I)
3330  CONTINUE
      WRITE(ITTO,3150) RERR
      IF(.NOT.LPLOTS) GO TO 3400
      DO 3340 I=1,NPOUT
      X(I) = STRT + REAL(I-1) * FINC
      TANTHT = TAN( X(I) * DTORD2 )
      Y(I) = 2.0*SQRT(WGPRO(3,1)+WGPRO(3,2)*TANTHT+WGPRO(3,3)*TANTHT**2)
3340  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(6),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,3112) NPLSET(6),NPOUT,0.0
        WRITE(IOPL,3114) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C *
C *  MOVE BACK EXPONENT TO DATA ARRAYS
3400  DO 3410 I=1,NPTS
      X(I) = WSGPAR(I,1)
      Y(I) = WSGPAR(I,6)
3410  CONTINUE
      IF(LPLOTS) THEN
        IF(KPTYPE.EQ.0) THEN
          WRITE(IOPL) NPLSET(7),NPTS,0.0
          WRITE(IOPL) (X(I),Y(I),I=1,NPTS)
         ELSE
          WRITE(IOPL,3112) NPLSET(7),NPTS,0.0
          WRITE(IOPL,3114) (X(I),Y(I),I=1,NPTS)
        ENDIF
      ENDIF
C *
C *  USE A 2ND ORDER FIT
      RERR = PFIT( X, Y, W, NPTS, 2, PYINT, BVECS, 0 )
C *
C *  WRITE THE RESULTS
      WRITE(ITTO,3120)
      DO 3430 I=1,3
      WGPRO(4,I) = BVECS(I)
      WRITE(ITTO,3140) I,WGPRO(4,I)
3430  CONTINUE
      WRITE(ITTO,3150) RERR
      IF(.NOT.LPLOTS) GO TO 3500
      DO 3440 I=1,NPOUT
      X(I) = STRT + REAL(I-1) * FINC
      Y(I) = WGPRO(4,1) + WGPRO(4,2)*X(I) + WGPRO(4,3)*X(I)**2
3440  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(8),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,3112) NPLSET(8),NPOUT,0.0
        WRITE(IOPL,3114) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C *
C *  CALIBRATION COMPLETE
3500  TFWPRO = .TRUE.
      IF(LPLOTS) CLOSE(UNIT=IOPL)
      GO TO 1000
C *
C *  EXIT STAGE RIGHT
9000  RETURN
      END
C ***  EVLCAL - EVALUATE DELTA TWO-THETA CALIBRATION CURVE
C *
      SUBROUTINE EVLCAL
C *
      LOGICAL LPLOTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1, QTITLE*80, QFILE*32, QTITLS(5)*20
      PARAMETER ( DTOR=1.745329252E-02, RTOD=57.29577951 )
      PARAMETER ( DTORD2=8.72664626E-3, RTODT2=114.591559 )
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      COMMON /PLTTYP/ KPTYPE
      DIMENSION NPLSET(5),NCPLOT(5),OFFSET(5),NUMPEN(5),MARKER(5)
      DATA QTITLS(1)/'Two-theta error     '/
      DATA QTITLS(2)/'Polynomial fit      '/
      DATA QTITLS(3)/'Sample Displacement '/
      DATA QTITLS(4)/'Axial Divergence    '/
      DATA QTITLS(5)/'Sample Transparency '/
C *
      NWGPAR = 0
C *
C *  DETERMINE FUNCTION
1000  WRITE(ITTO,1002)
1002  FORMAT(/,' * Two-theta calibration curve determination')
1004  WRITE(ITTO,1010)
1010  FORMAT(/,'   R - Read profile data file',/,
     1         '   P - Perform calibration',/,
     1         '   E - Exit this section',//,
     1         ' Choice ? ')
      READ(ITTI,1020) QANS
1020  FORMAT(A1)
C *
C *  DO IT
1200  IF( QANS.EQ.'R' .OR. QANS.EQ.'r' ) GO TO 2000
      IF( QANS.EQ.'P' .OR. QANS.EQ.'p' ) GO TO 3000
      IF( QANS.EQ.'E' .OR. QANS.EQ.'e' ) GO TO 9000
      GO TO 1004
C *
C **  READ PROFILE DATA FILE
C *
2000  IBEG = NWGPAR + 1
      CALL GETPRO
C *
C *  GET CALCULATED TWO-THETA VALUES FOR NEW LINES
      IEND = NWGPAR
      CALL ASS2TH( IBEG, IEND, WGRAD(1,1) )
      GO TO 1000
C *
C **  DETERMINE DELTA 2THETA CALIBRATION CURVE
C *
C *  GENERATE TWO-THETA , DELTA TWO-THETA ARRAYS
3000  NPTS = 0
      DO 3010 I=1,NWGPAR
      IF(WSGPAR(I,12).EQ.0.0) GO TO 3010
        NPTS = NPTS + 1
        X(NPTS) = WSGPAR(I,1)
        Y(NPTS) = WSGPAR(I,12) - WSGPAR(I,1)
3010  CONTINUE
      WRITE(ITTO,3020) NPTS
3020  FORMAT(/,' Lines for evaluation: ',I3)
      IF(NPTS.EQ.0) THEN
        WRITE(ITTO,3025)
3025    FORMAT(' --No lines to evaluate: try again')
        GO TO 1000
      ENDIF
C *
C *  DATA CAN BE PLOTTED.  CHECK FOR PLOTTING BY NAME...
C *
4000  WRITE(ITTO,4010)
4010  FORMAT(/,' Plot file name ? (* = NONE) <TPLOT> ')
      READ(ITTI,4020) QFILE
4020  FORMAT(A32)
      LPLOTS = QFILE(1:1).NE.'*'
      IF(.NOT.LPLOTS) GO TO 5000
      IF( QFILE(1:1).EQ.' ' .OR. QFILE(1:1).EQ.'/' ) QFILE = 'TPLOT'
C *
C *  GET A TITLE FOR THE DATA SET
      WRITE(ITTO,4052)
4052  FORMAT(/,' Enter a title for the data set ? ')
      READ(ITTI,4054) QTITLE
4054  FORMAT(A80)
C *
C *  OPEN THE PLOT FILE
      IF(KPTYPE.EQ.0) THEN
        OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='UNFORMATTED',
     1  ERR=4030)
       ELSE
        OPEN(UNIT=IOPL,FILE=QFILE,STATUS='UNKNOWN',FORM='FORMATTED',
     1  ERR=4030)
      ENDIF
      GO TO 4050
4030  WRITE(ITTO,4040)
4040  FORMAT(' --Can not open specified plot file')
      LPLOTS = .FALSE.
      GO TO 5000
C *
C *  SET DEFAULT PARAMETERS
C *  TWO-THETA ERROR OBS
4050  NPLSET(1) =  1
      NCPLOT(1) = -1
      OFFSET(1) =  0.0
      NUMPEN(1) =  1
      MARKER(1) =  1
C *  TWO-THETA ERROR POLYNOMIAL
      NPLSET(2) =  2
      NCPLOT(2) =  1
      OFFSET(2) =  0.0
      NUMPEN(2) =  1
      MARKER(2) =  1
C *  TWO-THETA ERROR SPECIMEN DISPLACEMENT
      NPLSET(3) =  3
      NCPLOT(3) =  1
      OFFSET(3) =  0.0
      NUMPEN(3) =  2
      MARKER(3) =  1
C *  TWO-THETA ERROR AXIAL DIVERGENCE
      NPLSET(4) =  4
      NCPLOT(4) =  1
      OFFSET(4) =  0.0
      NUMPEN(4) =  3
      MARKER(4) =  1
C *  TWO-THETA ERROR THICK SPECIMEN TRANSPARENCY
      NPLSET(5) =  5
      NCPLOT(5) =  1
      OFFSET(5) =  0.0
      NUMPEN(5) =  4
      MARKER(5) =  1
C *
      NUMSTS = 5
      NPLOTS = 1
      IAXIS = 0
      KYID = 2
      DPI = 2.0
      FS = 0.2
C *
C *  DETERMINE THE RANGE FOR THE PLOT FILE OUTPUT
      STRT = 180.0
      FINS =   0.0
      DO 4060 I=1,NPTS
      IF(STRT.GT.X(I)) STRT = X(I)
      IF(FINS.LT.X(I)) FINS = X(I)
4060  CONTINUE
      STRT = AINT( STRT/10.0 ) * 10.0
      FINS = AINT( FINS/10.0 ) * 10.0 + 10.0
      NPOUT = 100
      FINC = ( FINS - STRT ) / REAL( NPOUT - 1 )
C *
C *  WRITE THE INITIAL INFORMATION TO THE PLOT FILE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) QTITLE
        WRITE(IOPL) STRT,FINS,FINC
        WRITE(IOPL) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
        DO 4062 I=1,NUMSTS
4062    WRITE(IOPL) QTITLS(I),NCPLOT(I),OFFSET(I),NUMPEN(I),MARKER(I)
        WRITE(IOPL) NPLSET(1),NPTS,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPTS)
       ELSE
        WRITE(IOPL,4110) QTITLE
4110    FORMAT(A80)
        WRITE(IOPL,4120) STRT,FINS,FINC
4120    FORMAT(3F10.4)
        WRITE(IOPL,4130) NUMSTS,NPLOTS,IAXIS,DPI,FS,KYID
4130    FORMAT(3I9,2F9.2,I9)
        DO 4132 I=1,NUMSTS
4132    WRITE(IOPL,4140) QTITLS(I),NCPLOT(I),OFFSET(I),NUMPEN(I),
     1  MARKER(I)
4140    FORMAT(A20,I6,F6.2,2I6)
        WRITE(IOPL,4150) NPLSET(1),NPTS,0.0
4150    FORMAT(2I6,F10.4)
        WRITE(IOPL,4160)  (X(I),Y(I),I=1,NPTS)
4160    FORMAT(2F10.4)
      ENDIF
C---
C-----------------------------------------------------------------------
C---
C---  SECOND ORDER POLYNOMIAL FIT
C---
5000  RERR = PFIT( X, Y, W, NPTS, 2, PYINT, BVECS, 0 )
C *
C *  WRITE THE RESULTS
      WRITE(ITTO,5010)
5010  FORMAT(/,' The coefficients from the 2nd order polynomial fit:',/)
      DO 5020 I=1,3
      WGCAL(I) = BVECS(I)
      WRITE(ITTO,5030) I,WGCAL(I)
5030  FORMAT(' Parameter ',I2,': ',E13.6)
5020  CONTINUE
      WRITE(ITTO,5040) RERR
5040  FORMAT(/,' Standard deviation after refinement:   ',F10.6)
      TFWCAL = .TRUE.
C *
C *  WRITE CORRECTED VALUES
      WRITE(ITTO,5050)
5050  FORMAT(//,' The calibration curve uses the polynomial fit.',//,
     1' Corrected line positions based on the reference lines (*)')
      WRITE(ITTO,5052)
5052  FORMAT(/,
     1' Line  2tht obs     Corr     2tht corr   Delta    2tht calc',/,
     1' ----  ---------  ---------  ---------  --------  ---------')
      DO 5060 I=1,NWGPAR
      TTHT = WSGPAR(I,1)
      CORR = WGCAL(1) + WGCAL(2)*TTHT + WGCAL(3)*TTHT**2
      TTHTC = TTHT + CORR
      IF(WSGPAR(I,12).EQ.0.0) THEN
        QANS = ' '
        WRITE(ITTO,5062) I,TTHT,QANS,CORR,TTHTC
5062    FORMAT(1X,I3,3X,F8.4,A1,2X,F8.4,3X,F8.4,2X,F8.4,3X,F8.4)
       ELSE
        QANS = '*'
        DELTA = WSGPAR(I,12) - TTHTC
        WRITE(ITTO,5062) I,TTHT,QANS,CORR,TTHTC,DELTA,WSGPAR(I,12)
      ENDIF
5060  CONTINUE
C *
C *  GENERATE THE PLOTS
      IF(.NOT.LPLOTS) GO TO 5100
      DO 5070 I=1,NPOUT
      X(I) = STRT + REAL(I-1) * FINC
      Y(I) = WGCAL(1) + WGCAL(2)*X(I) + WGCAL(3)*X(I)**2
5070  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(2),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,4150) NPLSET(2),NPOUT,0.0
        WRITE(IOPL,4160) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C---
C-----------------------------------------------------------------------
C---
C---  SPECIMEN DISPLACEMENT ERROR CORRECTION
C---
5100  NPTS = 0
      DO 5110 I=1,NWGPAR
      IF(WSGPAR(I,12).EQ.0.0) GO TO 5110
        NPTS = NPTS + 1
        X(NPTS) = COS(DTORD2*WSGPAR(I,1))
        Y(NPTS) = WSGPAR(I,12) - WSGPAR(I,1)
5110  CONTINUE
      RERR = PFIT( X, Y, W, NPTS, 1, PYINT, BVECS, 0 )
      WRITE(ITTO,5120) BVECS(1),BVECS(2)
5120  FORMAT(//,' Results based on specimen displacement:',//,
     1          ' Zero point shift = ',F9.6,/,
     1          '            -2s/R = ',F9.6)
      WRITE(ITTO,5040) RERR
C *
C *  GENERATE THE PLOTS
      IF(.NOT.LPLOTS) GO TO 5200
      DO 5130 I=1,NPOUT
      ANGLE = STRT + REAL(I-1) * FINC
      X(I) = ANGLE
      Y(I) = BVECS(1) + BVECS(2) * COS(DTORD2*ANGLE)
5130  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(3),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,4150) NPLSET(3),NPOUT,0.0
        WRITE(IOPL,4160) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C---
C-----------------------------------------------------------------------
C---
C---  AXIAL DIVERGENCE ERROR CORRECTION
C---
5200  NPTS = 0
      DO 5210 I=1,NWGPAR
      IF(WSGPAR(I,12).EQ.0.0) GO TO 5210
        NPTS = NPTS + 1
        X(NPTS) = COS(DTOR*WSGPAR(I,1)) / SIN(DTOR*WSGPAR(I,1))
        Y(NPTS) = WSGPAR(I,12) - WSGPAR(I,1)
5210  CONTINUE
      RERR = PFIT( X, Y, W, NPTS, 1, PYINT, BVECS, 0 )
      WRITE(ITTO,5220) BVECS(1),BVECS(2)
5220  FORMAT(//,' Results based on axial divergence:',//,
     1          ' Zero point shift = ',F9.6,/,
     1          '     -h**2/3*R**2 = ',F9.6)
      WRITE(ITTO,5040) RERR
C *
C *  GENERATE THE PLOTS
      IF(.NOT.LPLOTS) GO TO 5300
      DO 5230 I=1,NPOUT
      ANGLE = STRT + REAL(I-1) * FINC
      X(I) = ANGLE
      Y(I) = BVECS(1) + BVECS(2) * COS(DTOR*ANGLE) / SIN(DTOR*ANGLE)
5230  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(4),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,4150) NPLSET(4),NPOUT,0.0
        WRITE(IOPL,4160) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C---
C-----------------------------------------------------------------------
C---
C---  THICK SPECIMEN TRANSPARENCY
C---
5300  NPTS = 0
      DO 5310 I=1,NWGPAR
      IF(WSGPAR(I,12).EQ.0.0) GO TO 5310
        NPTS = NPTS + 1
        X(NPTS) = SIN(DTOR*WSGPAR(I,1))
        Y(NPTS) = WSGPAR(I,12) - WSGPAR(I,1)
5310  CONTINUE
      RERR = PFIT( X, Y, W, NPTS, 1, PYINT, BVECS, 0 )
      WRITE(ITTO,5320) BVECS(1),BVECS(2)
5320  FORMAT(//,' Results based on thick specimen transparency:',//,
     1          ' Zero point shift = ',F9.6,/,
     1          '           -1/2uR = ',F9.6)
      WRITE(ITTO,5040) RERR
C *
C *  GENERATE THE PLOTS
      IF(.NOT.LPLOTS) GO TO 5400
      DO 5330 I=1,NPOUT
      ANGLE = STRT + REAL(I-1) * FINC
      X(I) = ANGLE
      Y(I) = BVECS(1) + BVECS(2) * SIN(DTOR*ANGLE)
5330  CONTINUE
      IF(KPTYPE.EQ.0) THEN
        WRITE(IOPL) NPLSET(5),NPOUT,0.0
        WRITE(IOPL) (X(I),Y(I),I=1,NPOUT)
       ELSE
        WRITE(IOPL,4150) NPLSET(5),NPOUT,0.0
        WRITE(IOPL,4160) (X(I),Y(I),I=1,NPOUT)
      ENDIF
C---
C-----------------------------------------------------------------------
C---
5400  CLOSE(UNIT=IOPL)
      GO TO 1000
C *
9000  RETURN
      END
C ***  GETPRO - GET PROFILE DATA FROM FILE
C *
      SUBROUTINE GETPRO
C *
      LOGICAL EXISTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1
      CHARACTER LINE*72, WGFILE*32
      CHARACTER FILEID*17, FILEXX*17, LINEID*10, LINEXX*10
      PARAMETER ( FILEID = '*** PROFILE DATA ' , LINEID = '** SET ID ' )
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      EQUIVALENCE ( LINE(2:18) , FILEXX )
      EQUIVALENCE ( LINE(2:11) , LINEXX )
C *
C **  READ FROM PROFILE DATA FILE
C *
1000  WRITE(ITTO,1010)
1010  FORMAT(/,' Profile data file name ? <WSGPRO> ')
      READ(ITTI,1020) WGFILE
1020  FORMAT(A32)
      IF( WGFILE.EQ.' ' .OR. WGFILE.EQ.'/' ) WGFILE = 'WSGPRO'
C *
C *  DETERMINE IF FILE IS A PROFILE DATA FILE
1100  OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='OLD',ERR=1102)
      GO TO 1106
1102  WRITE(ITTO,1104)
1104  FORMAT(/,' --File is busy or not there.')
      GO TO 9000
1106  READ(IDFI,1110,END=1120) LINE
1110  FORMAT(A72)
      IF(FILEID.EQ.FILEXX) GO TO 1310
1120  WRITE(ITTO,1130)
1130  FORMAT(' --File is not a profile data file')
      CLOSE(UNIT=IDFI,STATUS='KEEP')
      GO TO 9000
C *
C *  SCAN FILE FOR PROFILE DATA SETS...  READ IF INDICATED
1310  READ(IDFI,1110,END=1400) LINE
      IF(LINEID.NE.LINEXX) GO TO 1310
      WRITE(ITTO,1315) LINE(13:63)
1315  FORMAT(/,' Data set: ',A51)
      WRITE(ITTO,1320)
1320  FORMAT(' Read this profile data set ? <N> ')
      READ(ITTI,1322) QANS
1322  FORMAT(A1)
      IF( QANS.NE.'Y' .AND. QANS.NE.'y' ) THEN
1323    READ(IDFI,1110,END=1400) LINE
        IF(LINEID.NE.LINEXX) GO TO 1323
        GO TO 1310
      ENDIF
      READ(LINE(35:37),1325) KWSGLN
1325  FORMAT(I3)
      NTEMP = 0
      DO 1330 I=1,KWSGLN
      NWGPAR = NWGPAR + 1
      READ(IDFI,1340,END=1350) WSGPAR(NWGPAR,1),WSGPAR(NWGPAR,2),KPROF,
     1  (WSGPAR(NWGPAR,J),J=3,10)
      WSGPAR(NWGPAR,14) = REAL( KPROF )
1340  FORMAT(1X,F8.4,2X,F10.2,19X,I3,/,
     1  2(6X,E13.6,1X,E13.6,1X,E13.6,1X,E13.6,/))
      NTEMP = NTEMP + 1
1330  CONTINUE
      WRITE(ITTO,1332) NTEMP
1332  FORMAT(/,' The number of lines read from this set = ',I3)
      READ(IDFI,1110,END=1350) LINE
      GO TO 1310
C *
C *  PREMATURE EOF HIT
1350  WRITE(ITTO,1360) NTEMP
1360  FORMAT(' --Premature end-of-file hit:  ',I3,' lines used anyway.')
C *
C *  PROFILE DESCRIPTORS READ FROM FILE
1400  CLOSE(UNIT=IDFI,STATUS='KEEP')
C *
9000  RETURN
      END
C *** ASS2TH - DETERMINE CALCULATED TWO-THETAS FOR LINES
C *
      SUBROUTINE ASS2TH( IBEG, IEND, LAMBDA )
C *
C *  NLINES INDICATES TOTAL NUMBER OF LINES IN DDAT ARRAY.  IF FIRST
C *  VALUE IN DDAT IS 0.0, THEN SECOND IS INDEX INTO DDID ARRAY.
C *
      REAL LAMBDA
      CHARACTER DDID(3)*4
      PARAMETER ( NLINES = 12 )
      PARAMETER ( DELMAX=0.10 )
      PARAMETER ( DTORD2=8.72664626E-3, RTODT2=114.591559 )
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /PROPRM/ NWGPAR,WSGPAR(200,14)
      DIMENSION DDAT(NLINES,2)
      DATA DDID /'SI  ','AG  ','W   '/
      DATA DDAT( 1,1), DDAT( 1,2) /  0.00000,   1.0/
      DATA DDAT( 2,1), DDAT( 2,2) /  3.13552, 100.0/
      DATA DDAT( 3,1), DDAT( 3,2) /  1.92011,  55.0/
      DATA DDAT( 4,1), DDAT( 4,2) /  1.63747,  30.0/
      DATA DDAT( 5,1), DDAT( 5,2) /  1.35772,   6.0/
      DATA DDAT( 6,1), DDAT( 6,2) /  1.24593,  11.0/
      DATA DDAT( 7,1), DDAT( 7,2) /  1.10857,  12.0/
      DATA DDAT( 8,1), DDAT( 8,2) /  1.04517,   6.0/
      DATA DDAT( 9,1), DDAT( 9,2) /  0.96005,   3.0/
      DATA DDAT(10,1), DDAT(10,2) /  0.91799,   7.0/
      DATA DDAT(11,1), DDAT(11,2) /  0.85870,   8.0/
      DATA DDAT(12,1), DDAT(12,2) /  0.82820,   3.0/
C *
C *  SEARCH LINE LIST FOR MATCHES
2000  DO 2010 I=IBEG,IEND
      DOBS = LAMBDA / 2.0 / SIN( WSGPAR(I,1) * DTORD2 )
      WSGPAR(I,11) = DOBS
      WSGPAR(I,12) = 0.0
      WSGPAR(I,13) = 0.0
      ETEST = 1.0E+10
      DO 2020 J=1,NLINES
      IF(DDAT(J,1).EQ.0.0) THEN
        NID = INT(DDAT(J,2))
        GO TO 2020
      ENDIF
      TTHT = RTODT2 * ASIN( LAMBDA / 2.0 / DDAT(J,1) )
      DELTA = ABS( TTHT-WSGPAR(I,1) )
      IF( DELTA.LT.DELMAX .AND. DELTA.LT.ETEST ) THEN
        WSGPAR(I,12) = TTHT
        WSGPAR(I,13) = DDAT(J,1)
        WRITE(ITTO,2030) I,DDID(NID),WSGPAR(I,1),(WSGPAR(I,K),K=11,13)
2030    FORMAT(' Line ',I3,1X,A4,'   2TO ',F8.4,' DO ',F9.6,'   2TC ',
     1  F8.4,' DC ',F9.6)
      ETEST = DELTA
      ENDIF
2020  CONTINUE
2010  CONTINUE
C *
      RETURN
      END
C ***  WRIWSG - WRITE CALIBRATION INFORMATION TO CALIBRATION FILE
C *
      SUBROUTINE WRIWSG
C *
      LOGICAL EXISTS
      LOGICAL TFWBAK,TFWPRO,TFWWAV,TFWCAL
      CHARACTER QANS*1, QOPT*1, WGNAME*32
      CHARACTER QDATE*9, QTIME*8, LINE*72, WGFILE*32, FILEXX*20
      CHARACTER FILEID*20, BACKID*10, PROFID*10, WAVEID*10, CALBID*10
      PARAMETER ( FILEID = '*** INS CALIBRATION' )
      PARAMETER ( BACKID = 'BACKGROUND' )
      PARAMETER ( PROFID = 'PROFILES  ' )
      PARAMETER ( WAVEID = 'RADIATION ' )
      PARAMETER ( CALBID = 'DELTA 2THT' )
      COMMON /WSGFID/ WGFILE
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /WGFLAG/ TFWBAK,TFWPRO,TFWWAV,TFWCAL
      COMMON /WGDATA/ WGBAK(6),WGPRO(4,3),WGRAD(5,2),WGCAL(3)
      EQUIVALENCE ( LINE(1:20), FILEXX )
C *
C *  DETERMINE FUNCTION
1000  WRITE(ITTO,1002)
1002  FORMAT(/,' * Recording of calibration parameters')
1004  WRITE(ITTO,1010)
1010  FORMAT(/,'   B - Write background parameters',/,
     1         '   C - Write 2theta calibration parameters',/,
     1         '   R - Write radiation parameters',/,
     1         '   P - Write profile parameters',/,
     1         '   A - Write all sets of parameters',/,
     1         '   E - Return to main menu',//,
     1         ' Choice ? ')
      READ(ITTI,1020) QOPT
1020  FORMAT(A1)
      IF( QOPT.EQ.'E' .OR. QOPT.EQ.'e' ) GO TO 9000
      IF( QOPT.EQ.'B' .OR. QOPT.EQ.'b' ) GO TO 2000
      IF( QOPT.EQ.'C' .OR. QOPT.EQ.'c' ) GO TO 2000
      IF( QOPT.EQ.'R' .OR. QOPT.EQ.'r' ) GO TO 2000
      IF( QOPT.EQ.'P' .OR. QOPT.EQ.'p' ) GO TO 2000
      IF( QOPT.EQ.'A' .OR. QOPT.EQ.'a' ) GO TO 2000
      GO TO 1004
C *
C **  CONFIRM THE CALIBRATION FILE NAME
C *
2000  WRITE(ITTO,2010) WGFILE
2010  FORMAT(/,' Instrument calibration file name is: ',A32)
C *
C *  OPEN OLD AND A TEMPORARY FILE TO COPY EXISTING DATA
3000  OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='OLD',FORM='FORMATTED')
CVAX
C     OPEN(UNIT=ISCR,STATUS='SCRATCH',FORM='FORMATTED',
C    1CARRIAGECONTROL='LIST')
CVAX
CIBM
      OPEN(UNIT=ISCR,STATUS='SCRATCH',FORM='FORMATTED')
CIBM
C *
C *  COPY OLD FILE HEADER RECORDS
3100  READ(IDFI,3110,END=3150) LINE
3110  FORMAT(1X,A72)
      WRITE(ISCR,3110) LINE
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,3110) LINE
      GO TO 3500
3150  WRITE(ITTO,3160)
3160  FORMAT(' --Premature end of instrument calibration file: Abort.')
      CLOSE(UNIT=IDFI,STATUS='KEEP')
      CLOSE(UNIT=ISCR)
      GO TO 9000
C *
C *  WRITE THE NEW RADIATION WAVELENGTHS IF THEY HAVE BEEN DETERMINED
3500  IF(.NOT.TFWWAV) GO TO 3550
      IF( QOPT.NE.'R' .AND. QOPT.NE.'r' .AND.
     1    QOPT.NE.'A' .AND. QOPT.NE.'a') GO TO 3550
      WRITE(ITTO,3505)
3505  FORMAT(' -Writing radiation parameters')
      CALL DATE( QDATE )
      CALL TIME( QTIME )
      ZERO = 0.0
      DO 3510 I=1,5
3510  READ(IDFI,3110,END=3150) LINE
      DO 3520 I=1,5
3520  WRITE(ISCR,2250) WAVEID,WGRAD(I,1),WGRAD(I,2),ZERO,QDATE,QTIME
2250  FORMAT(1X,A10,1P,3(1X,E13.6),1X,A9,1X,A8)
      GO TO 3200
3550  DO 3560 I=1,5
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,3110) LINE
3560  CONTINUE
C *
C *  WRITE THE NEW BACKGROUND PARAMETERS IF THEY HAVE BEEN DETERMINED
3200  IF(.NOT.TFWBAK) GO TO 3250
      IF( QOPT.NE.'B' .AND. QOPT.NE.'b'  .AND.
     1    QOPT.NE.'A' .AND. QOPT.NE.'a') GO TO 3250
      WRITE(ITTO,3205)
3205  FORMAT(' -Writing background parameters')
      CALL DATE( QDATE )
      CALL TIME( QTIME )
      READ(IDFI,3110,END=3150) LINE
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,2250) BACKID,(WGBAK(I),I=1,3),QDATE,QTIME
      WRITE(ISCR,2250) BACKID,(WGBAK(I),I=4,6),QDATE,QTIME
      GO TO 3300
3250  DO 3260 I=1,2
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,3110) LINE
3260  CONTINUE
C *
C *  WRITE THE NEW PROFILE PARAMETERS IF THEY HAVE BEEN DETERMINED
3300  IF(.NOT.TFWPRO) GO TO 3350
      IF( QOPT.NE.'P' .AND. QOPT.NE.'p'  .AND.
     1    QOPT.NE.'A' .AND. QOPT.NE.'a') GO TO 3350
      WRITE(ITTO,3305)
3305  FORMAT(' -Writing profile parameters')
      CALL DATE( QDATE )
      CALL TIME( QTIME )
      READ(IDFI,3110,END=3150) LINE
      READ(IDFI,3110,END=3150) LINE
      READ(IDFI,3110,END=3150) LINE
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,2250) PROFID,(WGPRO(1,I),I=1,3),QDATE,QTIME
      WRITE(ISCR,2250) PROFID,(WGPRO(2,I),I=1,3),QDATE,QTIME
      WRITE(ISCR,2250) PROFID,(WGPRO(3,I),I=1,3),QDATE,QTIME
      WRITE(ISCR,2250) PROFID,(WGPRO(4,I),I=1,3),QDATE,QTIME
      GO TO 3400
3350  DO 3360 I=1,4
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,3110) LINE
3360  CONTINUE
C *
C *  WRITE THE NEW DELTA 2THETA PARAMETERS IF THEY HAVE BEEN DETERMINED
3400  IF(.NOT.TFWCAL) GO TO 3450
      IF( QOPT.NE.'C' .AND. QOPT.NE.'c' .AND.
     1    QOPT.NE.'A' .AND. QOPT.NE.'a') GO TO 3450
      WRITE(ITTO,3405)
3405  FORMAT(' -Writing delta 2theta calibration parameters')
      CALL DATE( QDATE )
      CALL TIME( QTIME )
      READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,2250) CALBID,(WGCAL(I),I=1,3),QDATE,QTIME
      GO TO 3600
3450  READ(IDFI,3110,END=3150) LINE
      WRITE(ISCR,3110) LINE
C *
C *  COPY THE TEMP FILE BACK TO THE PERMANENT FILE
3600  CLOSE(UNIT=IDFI)
      REWIND (UNIT=ISCR)
CVAX
C     OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='NEW',FORM='FORMATTED',
C    1CARRIAGECONTROL='LIST')
CVAX
CIBM
      OPEN(UNIT=IDFI,FILE=WGFILE,STATUS='UNKNOWN',FORM='FORMATTED')
CIBM
3610  READ(ISCR,3110,END=3699) LINE
      WRITE(IDFI,3110) LINE
      GO TO 3610
3699  CLOSE(UNIT=IDFI,STATUS='KEEP')
      CLOSE(UNIT=ISCR,STATUS='DELETE')
C *
C *  SUBJOB COMPLETE
9000  RETURN
      END
C ***  PLORNZ - LORENTZ/POLARIZATION FACTOR AT X TWO-THETA
C *
      FUNCTION PLORNZ( TWOTHT )
      PARAMETER ( DTOR=1.745329252E-02 )
      PARAMETER ( GC=0.8001525593 )
C *
C *  THIS IS THE LORENTZ-POLARIZATION CORRECTED FOR GR MONOCHROMETER
      ANG = TWOTHT * DTOR
      ANGD2 = ANG / 2.0
      PLORNZ  = (1.+GC*COS(ANG)**2) / (COS(ANGD2)*SIN(ANGD2)**2)
C *
      RETURN
      END
C ***  DATE - DUMMY ROUTINE
C *    LIMIT TO 9 CHARACTERS.  USED ONLY AS A TIME STAMP.
C *
C      SUBROUTINE DATE( QDATE )
C      CHARACTER QDATE*9
C        QDATE = 'DD-MMM-YY'
CIBMPC
C      INTEGER*2 IYEAR,IMONTH,IDAY
C      CHARACTER QMONTH(12)*3
C      DATA QMONTH/ 'Jan','Feb','Mar','Apr','May','Jun',
C     1             'Jul','Aug','Sep','Oct','Nov','Dec'/
C      CALL GETDAT( IYEAR, IMONTH, IDAY )
C      IYEAR = IYEAR - 1900
C      WRITE(QDATE,1) IDAY,IYEAR
C1     FORMAT(I2,'-',3X,'-',I2)
C      QDATE(4:6) = QMONTH(IMONTH)
CIBMPC
C      RETURN
C      END
C ***  TIME - DUMMY ROUTINE
C *    LIMIT TO 8 CHARACTERS.  USED ONLY AS A TIME STAMP.
C *
C      SUBROUTINE TIME( QTIME )
C      CHARACTER QTIME*8
C        QTIME = 'HH:MM:SS'
CIBMPC
C      INTEGER*2 IH,IM,IS,IHS
C      CALL GETTIM( IH, IM, IS, IHS )
C      WRITE(QTIME,1) IH,IM,IS
C1     FORMAT(I2,':',I2,':',I2)
C      IF(QTIME(4:4).EQ.' ') QTIME(4:4) = '0'
C      IF(QTIME(7:7).EQ.' ') QTIME(7:7) = '0'
CIBMPC
C      RETURN
C      END
C ***  CPUSEC - RETURNS CPU EXECUTION TIME IN SECONDS
C *
C *  ARGUMENT IS USED FOR DELTA TIME INTERVALS
C *
      FUNCTION CPUSEC( START )
C *
CVAX  ISTAT = LIB$STAT_TIMER(2,ITIME)
CVAX  TIME = REAL( ITIME ) / 100.0
CVAX  CPUSEC = TIME - START
      CPUSEC = 0.0
C *
      RETURN
      END
C ***  BSOLVE (33) - LEAST-SQUARES OPTIMIZATION ALGORITHM AFTER MARQUARDT
C *
      SUBROUTINE BSOLVE( IR, ICON, RERR, ITER )
      LOGICAL LTOCOM
      DOUBLE PRECISION PHI,PH,OBSTOT,DG,DN,GN,COSG,GAMM,DPSUM
      PARAMETER ( MEMSCR = 20480 )
      PARAMETER ( MEMV = 10 )
      PARAMETER ( MEMVP2 = MEMV + 2, MEMVT2 = MEMV * 2 )
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /SCRTCH/ SCRDAT(MEMSCR)
      DIMENSION AMAT(MEMV,MEMVP2),AKP2(MEMV),ACKP2(MEMV),PRMS(MEMVT2)
      DIMENSION ISENSE(10),ZDIFF(4096),YTMP(4096)
      DATA ISENSE/ 10*0 /
C *
C **  INITIALIZATION
C *
C *  SET ITERATION LIMIT BASED ON VALUE OF ITER
      ITERLM = 99
      IF( ITER.GT.1 .AND. ITER.LT.1000 ) ITERLM = ITER
      ITER = 0
C *
C *  DO NOT EXCEED THE MAXIMUM NUMBER OF VARIABLES
100   WRITE(ITTO,110) NVECS,ITERLM
110   FORMAT(/,' Variables for refinement: ',I3,'  Iteration',
     1' Limit: ',I3,/)
      IF(NVECS.LE.MEMV) GO TO 130
      WRITE(ITTO,120) MEMV
120   FORMAT('0-Number of variables exceeds ',I3,', Refinement skipped')
      ICON = -5
      GO TO 9000
C *
C *  OPEN THE SCRACTH FILE
130   OPEN(UNIT=ISCR,STATUS='SCRATCH',FORM='UNFORMATTED')
C *
C *  INITIALIZE VARIABLES
      K = NVECS
      N = NPTS
      KP1 = K + 1
      KP2 = K + 2
      FNU = 2.0
      FL = 1.00
      TAU = 0.001
      EPS = 0.0001
      PHMIN = 0.0
      ICON = 999
C *
C *  CAN PARAMETERS BE TRANSFERRED TO THE SCRATCH ARRAY?
      LTOCOM = ( K * KP1 ) .LE. MEMSCR
C *
C **  DETERMINE THE DENOMINATOR FOR THE ERROR
C *
1000  OBSTOT = 0.0
      IF(IR.EQ.0) THEN
        DO 1005 I=1,N
1005    OBSTOT = OBSTOT + Y(I)**2
       ELSE
1009    DO 1010 I=1,N
        OBSTOT = OBSTOT + WGT(I) * Y(I)**2
1010    CONTINUE
      ENDIF
C *
C *  COPY ORIGINAL PARAMETERS
1015  DO 1020 J1=1,K
      PRMS(J1) = BVECS(J1)
      J1PK = J1 + K
      PRMS(J1PK) = ABS( BVECS(J1) ) + 1.0E-2
1020  CONTINUE
C *
C *  DETERMINE BEGINNING ERROR
      GO TO 1700
C *
C **  EVALUATE DERIVATIVES
C *
C * GENERATE PATTERN BASED ON LAST PARAMETER VALUES
C1100  CALL GNBACK( PRMS, YC )
C *
C *  CALCULATE THE (WEIGHTED) DIFFERENCE AND STORE IN ZDIFF
1100  IF(IR.EQ.0) THEN
        DO 1104 I=1,N
1104    ZDIFF(I) = DBLE( Y(I) ) - DBLE( YC(I) )
       ELSE
        DO 1106 I=1,N
1106  ZDIFF(I) = DBLE( WGT(I) ) * ( DBLE( Y(I) ) - DBLE( YC(I) ) )
      ENDIF
C *
C *  INIT THE SCRATCH MEMORY
      CALL SCRMEM( FUNK1, JUNK2, JUNK3, JUNK4, 0 )
C *
C *  REPEAT FOR EACH PARAMETER
      DO 1110 J1=1,K
C *
C *  DETERMINE INCREMENTAL PARAMETER SHIFTS
      DO 1120 J2=1,K
1120  PRMS(J2) = BVECS(J2)
      DEN = 0.00001 * AMAX1( ABS(PRMS(J1)), PRMS(J1+K) )
      IF( (PRMS(J1)+DEN) .LE. BMAXS(J1) ) THEN
        PRMS(J1) = PRMS(J1) + DEN
       ELSE
        PRMS(J1) = PRMS(J1) - DEN
        DEN = -DEN
      ENDIF
C *
C *  EVALUATE FUNCTION WITH PARAMETER SHIFTS
1140  CALL GNBACK( PRMS, YTMP )
C *
C *  FIND BEGINNING AND ENDING NONZERO DERIVATIVE POINTERS
      DO 1160 I=1,N
1160  IF(YTMP(I).NE.YC(I)) GO TO 1165
      ISTART = 1
      IFINIS = 1
      GO TO 1190
1165  ISTART = I
      DO 1170 I=N,1,-1
1170  IF(YTMP(I).NE.YC(I)) GO TO 1175
      ISTART = 1
      I = 1
1175  IFINIS = I
C *
C *  CALCULATE DERIVATIVES
1190  DO 1150 J2=ISTART,IFINIS
1150  YTMP(J2) = ( YTMP(J2) - YC(J2) ) / DEN
C *
C *  WRITE THE DERIVATIVES TO THE SCRATCH MEMORY
      CALL SCRMEM( YTMP, ISTART, IFINIS, J1, 1 )
C *
C *  GET LAST MATRIX COLUMN
      DPSUM = 0.0D00
      DO 1225 J2=ISTART,IFINIS
1225  DPSUM = DPSUM + DPROD( YTMP(J2), ZDIFF(J2) )
      AMAT(J1,KP1) = SNGL( DPSUM )
C *
C *  NEXT PARAMETER
1110  CONTINUE
C *
C **  SET UP CORRECTION EQUATIONS
C *
C *  REPEAT FOR ALL PARAMETERS
1200  DO 1210 J1=1,K
C *
C *  RETRIEVE THE CORRECT ELEMENTS FROM THE SCRATCH MEMORY
      CALL SCRMEM( YTMP, ISTAR0, IFINI0, J1, -1 )
C *
C *  DETERMINE THE DIAGONAL COMPONENT
      DPSUM = 0.0D00
      IF(IR.EQ.0) THEN
        DO 1235 J3=ISTAR0,IFINI0
1235    DPSUM = DPSUM + DPROD( YTMP(J3), YTMP(J3) )
       ELSE
        DO 1236 J3=ISTAR0,IFINI0
1236    DPSUM = DPSUM + DBLE(WGT(J3)) * DPROD( YTMP(J3), YTMP(J3) )
      ENDIF
      AMAT(J1,J1) = SNGL( DPSUM )
      IF(J1.EQ.K) GO TO 1210
C *
C *  IF WGTING IS NECESSARY THEN DO IT NOW
      IF(IR.NE.0) THEN
        DO 1238 I=ISTAR0,IFINI0
1238    YTMP(I) = DPROD( YTMP(I), WGT(I) )
      ENDIF
C *
C *  DETERMINE OFF-DIAGONAL ELEMENTS
      DO 1230 J2=J1+1,K
      AMAT(J1,J2) = 0.0
      AMAT(J2,J1) = 0.0
      CALL SCRMEM( YC, ISTAR1, IFINI1, J2, -1 )
C *  DETERMINE OVERLAP REGION OF DERIVATIVE SETS
      I0 = IFINI0 - ISTAR1
      I1 = IFINI1 - ISTAR0
      IF( I0.LT.1 .OR. I1.LT.1 ) GO TO 1230
      ISTAR2 = MAX0( ISTAR0, ISTAR1 )
      IFINI2 = MIN0( IFINI0, IFINI1 )
C *  FIND NUMERICAL VALUES
      DPSUM = 0.0D00
      DO 1250 J3=ISTAR2,IFINI2
1250  DPSUM = DPSUM + DPROD( YTMP(J3), YC(J3) )
      AMAT(J1,J2) = SNGL( DPSUM )
      AMAT(J2,J1) = SNGL( DPSUM )
1230  CONTINUE
1210  CONTINUE
C *
C *  CORRECT FOR NEAR ZERO DIAGONAL ELEMENT
      DO 1260 J1=1,K
      IF(AMAT(J1,J1).GT.1.0E-20) GO TO 1260
      IF(ITER1.EQ.0) WRITE(ITTO,1265) J1,ITER
1265  FORMAT('   Parameter: ',I3,'  Eliminated on cycle: ',I4)
      DO 1270 J3=1,KP1
1270  AMAT(J1,J3) = 0.0
      AMAT(J1,J1) = 1.0
1260  CONTINUE
C *
C *  IF JUST SETTING UP MATRIX THEN EXIT HERE
C      IF(ITER1.EQ.1) GO TO 2320
C *
C *  CALCULATE VECTOR COMPONENT
      GN = 0.0
      DO 1255 J1=1,K
1255  GN = GN + AMAT(J1,KP1)**2
C *
C **  SCALE CORRECTION EQUATIONS
C *
1300  DO 1310 J1=1,K
1310  AKP2(J1) = SQRT( AMAT(J1,J1) )
      DO 1315 J1=1,K
      AMAT(J1,KP1) = AMAT(J1,KP1) / AKP2(J1)
      DO 1315 J2=1,K
      AMAT(J1,J2) = AMAT(J1,J2) / ( AKP2(J1)*AKP2(J2) )
1315  CONTINUE
C *
C *  VERIFY MATRIX
      IF(ISENSE(4).GT.0) THEN
        WRITE(ITTO,1317) (J1,J1=1,K)
1317    FORMAT(/,6X,30(1X,I3))
        DO 1320 J1=1,K
1320    WRITE(ITTO,4) J1,(NINT(AMAT(J1,LL)*100.0),LL=1,K)
4       FORMAT(1X,I3,2X,30(1X,I3))
      ENDIF
C *
C **  DETERMINE LAMBDA AND ADD TO DIAGONAL ELEMENTS
C *
C *  WRITE THE A MATRIX TO THE SCRATCH FILE FOR POSSIBLE LATER USE
1400  IF(LTOCOM) THEN
        NTEMP = 0
        DO 1402 J1=1,K
        DO 1402 J2=1,KP1
        NTEMP = NTEMP + 1
        SCRDAT(NTEMP) = AMAT(J1,J2)
1402    CONTINUE
       ELSE
        REWIND (UNIT=ISCR)
        WRITE(ISCR) ((AMAT(J1,J2),J1=1,K),J2=1,KP1)
      ENDIF
C *
C *  DECREASE THE SIZE OF FL
      FL = FL / FNU
C *
C *  ADD TO DIAGONAL
1415  DO 1425 J1=1,K
1425  AMAT(J1,J1) = AMAT(J1,J1) + FL
C *
C **  SOLVE CORRECTION EQUATIONS
C *
1500  DO 1510 L1=1,K
      L1P1 = L1 + 1
      DO 1515 L2=L1P1,KP1
      IF(AMAT(L1,L1).EQ.0.0) THEN
        WRITE(ITTO,1502)
1502    FORMAT(' -BSOLVE: Avoiding zero division in matrix solution.')
        GO TO 1516
      ENDIF
1515  AMAT(L1,L2) = AMAT(L1,L2) / AMAT(L1,L1)
1516  DO 1520 L3=1,K
      IF( L1-L3 ) 1530,1520,1530
1530  DO 1525 L4=L1P1,KP1
1525  AMAT(L3,L4) = AMAT(L3,L4) - AMAT(L1,L4) * AMAT(L3,L1)
1520  CONTINUE
1510  CONTINUE
C *
C **  SOLVE FOR NEW SOLUTION VECTOR
C *
1600  DN = 0.0
      DG = 0.0
      DO 1610 J1=1,K
      ACKP2(J1) = AMAT(J1,KP1) / AKP2(J1)
      IF(ISENSE(4).EQ.2) WRITE(ITTO,6) J1,ACKP2(J1)
6     FORMAT(' PARAMETER ',I3,'  ADJUSTMENT ',E13.6)
      BPADJ =  BVECS(J1) + ACKP2(J1)
      IF(BPADJ.GT.BMAXS(J1)) WRITE(ITTO,10) J1,ITER,BPADJ,BMAXS(J1)
10    FORMAT('   Param ',I3,' Iter ',I3,': Value ',E13.6,
     1       ' Max limit ',E13.6)
      BPADJ = AMIN1( BMAXS(J1), BPADJ )
      IF(BPADJ.LT.BMINS(J1)) WRITE(ITTO,11) J1,ITER,BPADJ,BMINS(J1)
11    FORMAT('   Param ',I3,' Iter ',I3,': Value ',E13.6,
     1       ' Min limit ',E13.6)
      PRMS(J1) = AMAX1( BMINS(J1),BPADJ )
      DG = DG + ACKP2(J1) * AMAT(J1,KP1) * AKP2(J1)
      DN = DN + ACKP2(J1) * ACKP2(J1)
1610  CONTINUE
      COSG = DG / DSQRT( DN*GN )
      IF(ISENSE(4).EQ.2) WRITE(ITTO,7) COSG
7     FORMAT(' COSG = ',E13.6)
      JGAM = 0
      IF( COSG ) 1615,1620,1620
1615  JGAM = 2
      COSG = -COSG
1620  IF(COSG.GT.1.0) COSG = 1.00
      GAMM = DATAN( DSQRT(1.00-COSG*COSG)/COSG )
      GAMM = GAMM * 180.00 / 3.1415926531615
      IF(JGAM.GT.0) GAMM = 180.0 - GAMM
      IF(ISENSE(4).EQ.2) WRITE(ITTO,5)GAMM
5     FORMAT(' GAMM = ',E13.6)
C *
C **  EVALUATE FUNCTION AND FIND NEW ERROR
C *
C *  GENERATE PROFILES WITH NEW PARAMETERS
1700  CALL GNBACK( PRMS, YC )
C *
C *  EVALUATE SUM-OF-SQUARES ERROR
      PHI = 0.0
      IF(IR.EQ.0) THEN
        DO 1710 J1=1,N
1710    PHI = PHI + ( Y(J1) - YC(J1) )**2
       ELSE
        DO 1760 J1=1,N
1760    PHI = PHI + WGT(J1) * ( Y(J1) - YC(J1) )**2
      ENDIF
1780  PHI = DSQRT( PHI / OBSTOT )
      WRITE(ITTO,2) ITER,PHI,FL
2     FORMAT(' Iter = ',I3,'  Error =',F9.6,'  Weight = ',F11.6)
C++++++++
C      OPEN(UNIT=2,STATUS='NEW')
C      WRITE(2,2) ITER,PHI,RLX
C      CLOSE(UNIT=2)
C++++++++
      IF(PHI.LT.1.0E-10) GO TO 2200
      IF(ITER.EQ.0) GO TO 2120
      IF(PHI.GE.PH) GO TO 2000
C *
C *  EPSILON TEST
C *
C *  TEST FOR CONVERGENCE BASED ON SMALL DEVIATION OF PARAMETERS
1800  ICON = 0
      DO 1810 J1=1,K
      DELTA = ABS( PRMS(J1) - BVECS(J1) )
      IF( DELTA/(TAU+ABS(PRMS(J1))) .GT. EPS ) ICON = ICON + 1
1810  CONTINUE
      IF(ISENSE(4).EQ.2) WRITE(ITTO,3)ICON
3     FORMAT(' ICON = ',I4)
      IF(ICON.NE.0) GO TO 1900
C *
C *  PERFORM GAMMA EPSILON TEST (IS LAMBDA LARGE?)
      IF( FL.GT.1.0 .AND. GAMM.LE.45.0 ) ICON = -3
      IF(ICON.EQ.-3) WRITE(ITTO,1820) ICON,GAMM,FL
1820  FORMAT(' *BSOLVE* GAMM,FL = ',I2,2X,F10.4,1X,E12.6)
      GO TO 2100
C *
C *  GAMMA LAMBDA TEST
C *
C *  HAVE THE TWO SOLUTIONS VECTORS DIVERGED?
1900  IF( FL.GT.1.0 .AND. GAMM.GT.90.0 ) ICON = -2
      IF(ICON.EQ.-2) WRITE(ITTO,1825) ICON,GAMM,FL
1825  FORMAT(' *BSOLVE ABORT* GAMM,FL = ',I2,2X,F10.4,1X,E12.6)
      GO TO 2100
C *
C **  ERROR INCREASED... TAKE APPROPRIATE ACTION
C *
C *  INCREASE FL IF THE LIMIT HAS NOT BEEN HIT
2000  IF(FL.GT.1000.0) THEN
        ICON = -1
        WRITE(ITTO,2010)
2010    FORMAT(' -Maximum limit for diagonal weight hit:  Halting',
     1  ' refinement')
        PHI = PH
        GO TO 2200
      ENDIF
      FL = FL * FNU
C *
C *  RETRIEVE MATRIX
      IF(LTOCOM) THEN
        NTEMP = 0
        DO 2020 J1=1,K
        DO 2020 J2=1,KP1
        NTEMP = NTEMP + 1
        AMAT(J1,J2) = SCRDAT(NTEMP)
2020    CONTINUE
       ELSE
        REWIND (UNIT=ISCR)
        READ(ISCR) ((AMAT(J1,J2),J1=1,K),J2=1,KP1)
      ENDIF
      GO TO 1415
C *
C *  UPDATE VALUES FOR PARAMETERS AND BEST ERROR PH
2100  DO 2110 J2=1,K
2110  BVECS(J2) = PRMS(J2)
2120  PH = PHI
C *
C *  TEST FOR CONVERGENCE
      IF(ICON.LE.0) GO TO 2200
C *
C *  NEXT ITERATION IF THE LIMIT HAS NOT BEEN HIT
      ITER = ITER + 1
      IF(ITER.LE.ITERLM) GO TO 1100
      WRITE(ITTO,2140)
2140  FORMAT(' -Iteration count limit hit:  Halting refinement')
      ITER = ITER - 1
      ICON = -5
C *
C *  RETURN WITH RESULTS
2200  RERR = SNGL( PHI )
C *
C *  REPEAT FOR ESTIMATION OF DEVIATIONS
C      IF(ICON.NE.0) THEN
C         WRITE(ITTO,2325)
C2325     FORMAT(' -Not proper convergence: No ESDs calculated.')
C         GO TO 9000
C      ENDIF
C      ITER1 = 1
C      GO TO 1100
C2320  CALL CHOL( AMAT, MEMV, MEMVP1, K, BVECS )
C *
9000  CLOSE(UNIT=ISCR)
      END
C ***  SCRMEM - HANDLE SCRATCH MEMORY REQUIREMENTS FOR REF. ALGORITHMS
C *
      SUBROUTINE SCRMEM( DATAIN, IBEG, IEND, ISET, ICODE )
C *
C *  ICODE:  0 -> INITIALIZE
C *         +1 -> STORE
C *         -1 -> RETRIEVE
C *
C *  SET PARAMETER MEMSCR TO WHATEVER ROOM IS AVAILABLE/NECESSARY.
C *  SETTING MEMSCR TO 1 HAS ALL DERIVATIVES WRITTEN TO SCRATCH FILE.
C *  VIRTUAL MACHINES.... THIS ONE'S FOR YOU!
C *
      PARAMETER ( MEMV = 10 )
      PARAMETER ( MEMSCR = 20480 )
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      COMMON /SCRTCH/ SCRDAT(MEMSCR)
      DIMENSION DATAIN(2), IPOINT(MEMV,4)
C *
C **  RESET POINTERS ( ICODE = 0 )
C *
1000  IF(ICODE.NE.0) GO TO 2000
      DO 1010 I=1,MEMV
      IPOINT(I,1) = 0
      IPOINT(I,2) = 0
      IPOINT(I,3) = 0
      IPOINT(I,4) = 0
1010  CONTINUE
      IMLAST = 0
      REWIND (UNIT=ISCR)
      IFLAST = 0
      IFCURR = 0
      GO TO 9000
C *
C **  MOVE DATA FROM *DATAIN* TO *SCRDAT* OR *DISK* IF *SCRDAT* IS FULL
C *
2000  IF(ICODE.NE.1) GO TO 3000
C *
C *  CHECK THAT INDICATED DATA SET IS EMPTY
      IF( IPOINT(ISET,1).EQ.0 .AND. IPOINT(ISET,2).EQ.0 ) GO TO 2100
      WRITE(ITTO,2010) ISET
2010  FORMAT(' --Write requested to set ',I3,' but data set exists.')
      STOP
C *
C *  DETERMINE IF THIS DATA SET CAN BE STORED IN *SCRDAT*
2100  NSTORE = IEND - IBEG + 1
      IF(IMLAST+NSTORE.LE.MEMSCR) GO TO 2900
C *
C *  THE CURRENT FILE-SET MUST BE THE SAME AS THE LAST-FILE SET.
C *  WRITES TO FILE MUST BE SEQUENTIAL
      IF(IFLAST.EQ.IFCURR) GO TO 2200
      WRITE(ITTO,2030) ISET,IFLAST,IFCURR
2030  FORMAT(' -Invalid request to write set ',I3,'  IFLAST ',I3,
     1  '  IFCURR ',I3)
      STOP
C *
C *  STORE THE DATA ON *DISK*
2200  IFLAST = IFLAST + 1
      IFCURR = IFLAST
      IPOINT(ISET,1) = -IFCURR
      IPOINT(ISET,2) = -IFCURR
      WRITE(ISCR) IBEG,IEND,(DATAIN(I),I=IBEG,IEND)
      IPOINT(ISET,3) = IBEG
      IPOINT(ISET,4) = IEND
      GO TO 9000
C *
C *  STORE THE REQUESTED DATA IN *SCRDAT*
2900  IPOINT(ISET,1) = IMLAST + 1
      DO 2910 I=IBEG,IEND
      IMLAST = IMLAST + 1
      SCRDAT(IMLAST) = DATAIN(I)
2910  CONTINUE
      IPOINT(ISET,2) = IMLAST
      IPOINT(ISET,3) = IBEG
      IPOINT(ISET,4) = IEND
      GO TO 9000
C *
C **  RETRIEVE DATA FROM *SCRDAT* OR *DISK*
C *
3000  IF(ICODE.NE.-1) GO TO 8000
C *
C *  DETERMINE IF DATA SET EXISTS
      IF( IPOINT(ISET,1).NE.0 .AND. IPOINT(ISET,2).NE.0 ) GO TO 3100
      WRITE(ITTO,3010) ISET
3010  FORMAT(' --Request for data set ',I3,' invalid: No such set.')
      STOP
C *
C *  DETERMINE IF DATA SET IS IN *SCRDAT* OR ON *DISK*
3100  IF( IPOINT(ISET,1).GT.0 .AND. IPOINT(ISET,2).GT.0 ) GO TO 3900
C *
C *  DETERMINE WHICH SET ON *DISK* IS TO BE RETRIEVED
      NDISK = -IPOINT(ISET,1)
C *
C *  IF THE CURRENT-SET IS ONE LESS THAN NDISK THEN NEXT READ IS CORRECT
      IF(IFCURR.EQ.NDISK-1) GO TO 3200
C *
C *  IF IFCURR IS > OR = TO NDISK THEN WE MUST REWIND
      IF(IFCURR.GE.NDISK) REWIND (UNIT=ISCR)
      IFCURR = 0
C *
C *  DETERMINE NUMBER OF SETS TO BYPASS FOR CORRECT SET
      NSKIP = NDISK - IFCURR - 1
      IF(NSKIP.GT.0) THEN
        DO 3120 I=1,NSKIP
        READ(ISCR) AJUNK
        IFCURR = IFCURR + 1
3120    CONTINUE
      ENDIF
C *
C *  READ THE DATA SET
3200  READ(ISCR) IBEG,IEND,(DATAIN(I),I=IBEG,IEND)
      IFCURR = IFCURR + 1
      GO TO 9000
C *
C *  RETRIEVE THE DATA FROM *SCRDAT*
3900  IBEG = IPOINT(ISET,3)
      IEND = IPOINT(ISET,4)
      NPNT = IPOINT(ISET,1)
      DO 3910 I=IBEG,IEND
      DATAIN(I) = SCRDAT(NPNT)
      NPNT = NPNT + 1
3910  CONTINUE
      GO TO 9000
C *
C **  UNKNOWN REQUEST CODE
C *
8000  WRITE(ITTO,8010) ICODE
8010  FORMAT(' --Unknown request code:',I6)
C *
C *  RETURN RETURN RETURN
9000  RETURN
      END
C *** PFIT
C *
C * FUNCTION PFIT FITS AN NTH DEGREE LINEAR LEAST-SQUARES POLYNOMIAL
C * CURVE TO THE M DATA POINTS CONTAINED IN THE X AND Y ARRAYS.
C * RETURNED ARE THE STANDARD DEVIATION OF THE POINTS ABOUT THE
C * REGRESSION CURVE AS THE VALUE OF THE FUNCTION AND THE Y INTERCEPT
C * AND POLYNOMIAL COEFFICIENTS IN THE CALLING ARGUMENTS.  IF THE VALUE
C * OF N IS POSITIVE THEN THE POLYNOMIAL IS EVALUATED SO THE Y INTERCEPT
C * IS REFINED AS A VARIABLE.  HOWEVER, IF THE VALUE OF N IS NEGATIVE
C * THEN THE Y INTERCEPT IS PINNED AT THE VALUE YINTCP AND THE SYSTEM OF
C * SIMULTANEOUS USED TO DETERMINE THE REGRESSION COEFFICIENTS ARE
C * EVALUATED ACCORDINGLY.
C *
C *  CALLING ARGUMENTS:
C *
C *     X      = X DATA ARRAY
C *     Y      = Y DATA ARRAY
C *     W      = OPTIONAL WEIGHTING FOR Y VALUES
C *     M      = NUMBER OF POINTS IN REGRESSION
C *     N      = DEGREE OF POLYNOMIAL FIT.  IF N IS POSITIVE THE Y
C *              INTERCEPT VALUE IS REFINED.  IF N IS NEGATIVE THEN
C *              THE POLYNOMIAL CURVE WILL BE REFINED SUCH THAT
C *              THE VALUE OF YINTCP WILL BE THE Y INTERCEPT.
C *     PYINT  = THE PINNED Y INTERCEPT
C *     IWGT   = 0 - NO WEIGHTING  1 - WEIGHTING
C *
C *  RETURN ARGUMENTS:
C *
C *     PFIT   = STANDARD DEVIATION OF THE POINTS ABOUT THE
C *              REGRESSION CURVE.
C *     PB     = THE REGRESSION COEFFICIENTS, LOW TO HIGH POWER.
C *
      FUNCTION PFIT( X, Y, W, M, N, PYINT, PB, IWGT )
C *
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL PFIT, X(M), Y(M), W(M), PYINT, PB(10)
      DIMENSION C(11,11), SX(20), SYX(10), CYX(10), B(10)
C *
C *  NDEG IS THE DEGREE OF THE POLYNOMIAL TO FIT
      NDEG = IABS( N )
C *
C *  INITIALIZE THE FUNCTION VALUE TO A LARGE NUMBER
      PFIT = 1.0E15
C *
C *  THE NUMBER OF POINTS MUST BE TWO GREATER THAN THE POWER
      IF(M.LT.NDEG+2) RETURN
C *
C *  FIRST THROUGH TENTH DEGREE FIT ONLY
      IF( NDEG.LT.1 .OR. NDEG.GT.10 ) STOP'PFIT: BAD N'
C *
C *  COMPUTE CONSTANTS
      FM = DBLE( M )
      NTWO = 2 * NDEG
      NP1 = NDEG + 1
      NP2 = NDEG + 2
C *
C *  INITIALIZE VARIABLES
      SY = 0.0
      SYY = 0.0
      DO 1 I=1,NDEG
      NPI = NDEG + I
      SX(I) = 0.0
      SX(NPI) = 0.0
1     SYX(I) = 0.0
C *
C *  COMPUTE SUMS AND PRODUCTS WITHOUT WEIGHTS
2000  IF(IWGT.EQ.0) THEN
        DO 2010 I = 1,M
        SY = SY + Y(I)
        SYY = SYY + ( Y(I) * Y(I) )
        DUM = 1.0
        DO 2020 J = 1,NDEG
        DUM = DUM * X(I)
        SX(J) = SX(J) + DUM
2020    SYX(J) = SYX(J) + ( Y(I) * DUM )
        DO 2010 J = NP1,NTWO
        DUM = DUM * X(I)
2010    SX(J) = SX(J) + DUM
       ELSE
C *  COMPUTE SUMS AND PRODUCTS WITH WEIGHTS
        FM = 0.0
        DO 2030 I = 1,M
        FM = FM + W(I)
        SY = SY + ( W(I) * Y(I) )
        SYY = SYY + ( W(I) * Y(I) * Y(I) )
        DUM = W(I)
        DO 2040 J = 1,NDEG
        DUM = DUM * X(I)
        SX(J) = SX(J) + DUM
2040    SYX(J) = SYX(J) + ( Y(I) * DUM )
        DO 2030 J = NP1,NTWO
        DUM = DUM * X(I)
2030    SX(J) = SX(J) + DUM
      ENDIF
C *
C *  EXECUTE ACCORDING TO THE SIGN OF N.  IF POSITIVE THEN REFINE THE
C *  Y INTERCEPT VALUE.  IF NOT THEN NO REFINEMENT IS NECESSARY.
      IF(N.LT.0) GO TO 100
C *
C *  CREATE COEFFICIENT MATRIX
      DO 5 I = 1,NDEG
      CYX(I) = SYX(I) - ( SY * SX(I) / FM )
      C(I,NP1) = CYX(I)
      DO 5 J = 1,NDEG
      IPJ = I + J
5     C(I,J) = SX(IPJ) - ( SX(I) * SX(J) / FM )
C *
C *  CALL SIMUL TO SOLVE SIMULTANEOUS EQUATIONS
      DET = SIMUL( NDEG, C, B, 11 )
C     WRITE(7,300)DET
C300  FORMAT(' THE DETERMINANT EQUALS ',F10.2)
      IF(DET.NE.0.0) GO TO 8
      RETURN
C *
C *  COMPUTE INTERCEPT YINT AND STANDARD DEVIATION SDEV
8     CYY = SYY - SY * SY / FM
      DO 9 I=1,NDEG
      SY = SY - B(I) * SX(I)
9     CYY = CYY - B(I) * CYX(I)
      YINT = SY / FM
      DENOM = M - NDEG - 1
      PFIT = SNGL( DSQRT( DABS( CYY / DENOM ) ) )
      GO TO 200
C *
C *  COMPUTE THE POLYNOMIAL WITH THE PINNED INTERCEPT
C *
C *  CREATE COEFFICIENT MATRIX
100   DO 110 I=1,NDEG
      CYX(I) = SYX(I) - SY * SX(I) / FM
      DO 120 J=1,NDEG
      IPJ = I + J
120   C(I,J) = SX(IPJ)
110   C(I,NP1) = SYX(I) - DBLE(PYINT) * SX(I)
C *
C *  CALL SIMUL TO SOLVE SIMULTANEOUS EQUATIONS
      DET = SIMUL( NDEG, C, B, 11 )
C     WRITE(7,300) DET
      IF(DET.NE.0.0) GO TO 130
      RETURN
C *
C *  RETURN Y INTERCEPT VALUE
130   YINT = DBLE( PYINT )
C *
C *  COMPUTE STANDARD DEVIATION
      CYY = SYY - SY * SY / FM
      DO 150 I=1,NDEG
150   CYY = CYY - B(I) * CYX(I)
      DENOM = M - NDEG - 1
      PFIT = SNGL( DSQRT( DABS( CYY / DENOM ) ) )
C *
C *  TRANSFER THE COEFFICIENTS
C *
200   PB(1) = SNGL( YINT )
      DO 210 I=1,NDEG
210   PB(I+1) = SNGL( B(I) )
C *
      RETURN
      END
C ***  FUNCTION SIMUL
C *
C *  FUNCTION SIMUL SOLVES SIMULTANEOUS LINEAR EQUATIONS WHICH ARE
C *  SOLVED IN PLACE. THE GAUSS-JORDEN COMPLETE ELIMINATION METHOD
C *  IS EMPLOYED WITH THE MAXIMUM PIVOT STRATEGY.  ROW AND COLUMN
C *  SUBSCRIPTS OF SUCCESSIVE PIVOT ELEMENTS ARE SAVED IN ORDER IN
C *  THE IROW AND JCOL ARRAYS RESPECTIVELY.  K IS THE PIVOT COUNTER,
C *  PIVOT THE ALGEBRAIC VALUE OF THE PIVOT ELEMENT, MAX
C *  THE NUMBER OF COLUMNS IN A AND DETER THE DETERMINANT OF THE
C *  COEFFICIENT MATRIX.  THE SOLUTIONS ARE COMPUTED IN THE (N+1)TH
C *  COLUMN OF A AND THEN UNSCRAMBLED AND PUT IN PROPER ORDER IN
C *  X(1)...X(N) USING THE PIVOT SUBSCRIPT INFORMATION AVAILABLE
C *  IN THE IROW AND JCOL ARRAYS,  THE SIGN OF THE DETERMINANT IS
C *  ASJUSTED, IF NECESSARY, BY DETERIMINING IF AN ODD OR EVEN NUMBER
C *  OF PAIRWISE INTERCHANGES IS REQUIRED TO PUT THE ELEMENTS OF THE
C *  JORD ARRAY IN ASCENDING SEQUENCE WHERE JORD(IROW(I))=JCOL(I).
C *  FUNCTION SIMUL COMPUTES THE VALUE OF THE DETERMINANT AND
C *  IS RETURNED AS THE VALUE OF THE FUNCTION.  SHOULD THE POTENTIAL
C *  PIVOT OF LARGEST MAGNITUDE BE SMALLER IN MAGNITUDE THAN EPS,
C *  THE MATRIX IS CONSIDERED TO BE SINGULAR AND A TRUE ZERO IS
C *  RETURNED AS THE VALUE OF THE FUNCTION.
C *
      FUNCTION SIMUL( N, A, X, NRC )
C *
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /LUNS/ ITTI,ITTO,IDFI,ISCR,IOPL
      DIMENSION IROW(30),JCOL(30),JORD(30),A(NRC,NRC),X(N)
      DATA EPS/1.00E-20/
C *
      MAX=N+1
C *
C *  INSURE N IS LESS THAN 30
      IF(N.GT.30) STOP'SIMUL: N TOO LARGE'
C *
C *  BEGIN ELIMINATION PROCEDURE
5     DETER=1.0
      DO 18 K=1,N
      KM1=K-1
C *
C *  SEARCH FOR PIVOT ELEMENT
      PIVOT=0.0
      DO 11 I=1,N
      DO 11 J=1,N
C *  SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN=1,KM1
      DO 8 JSCAN=1,KM1
      IF(I.EQ.IROW(ISCAN))GO TO 11
8     IF(J.EQ.JCOL(JSCAN))GO TO 11
9     IF(DABS(A(I,J)).LE.DABS(PIVOT))GO TO 11
      PIVOT=A(I,J)
      IROW(K)=I
      JCOL(K)=J
11    CONTINUE
C *
C *  INSURE THAT THE SELECTED PIVOT ELEMENT IS LARGER THAN EPS
      IF(DABS(PIVOT).GT.EPS)GO TO 13
      WRITE(ITTO,100)
100   FORMAT(/,' -SIMUL: WARNING - PIVOT ELEMENT TOO SMALL')
      SIMUL=0.0
      RETURN
C *
C *  UPDATE THE DETERMINATE VALUE
13    IROWK=IROW(K)
      JCOLK=JCOL(K)
      DETER=DETER*PIVOT
C *
C *  NORMALIZE PIVOT ROW ELEMENTS
      DO 14 J=1,MAX
14    A(IROWK,J)=A(IROWK,J)/PIVOT
C *
C *  CARRY OUT ELIMINATION AND CARRY OUT INVERSE
      A(IROWK,JCOLK)=1.0/PIVOT
      DO 18 I=1,N
      AIJCK=A(I,JCOLK)
      IF(I.EQ.IROWK)GO TO 18
      A(I,JCOLK)=-AIJCK/PIVOT
      DO 17 J=1,MAX
17    IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J)
18    CONTINUE
C *
C *  ORDER SOLUTION VALUES AND CREATE JORD ARRAY
      DO 20 I=1,N
      IROWI=IROW(I)
      JCOLI=JCOL(I)
      JORD(IROWI)=JCOLI
20    X(JCOLI)=A(IROWI,MAX)
C *
C *  ADJUST THE SIGN OF THE DETERMINANT
      INTCH=0
      NM1=N-1
      DO 22 I=1,NM1
      IP1=I+1
      DO 22 J=IP1,N
      IF(JORD(J).GE.JORD(I))GO TO 22
      JTEMP=JORD(J)
      JORD(J)=JORD(I)
      JORD(I)=JTEMP
      INTCH=INTCH+1
22    CONTINUE
      IF(INTCH/2*2.NE.INTCH) DETER=-DETER
C *
C *  RETURN WITH RESULTS
24    SIMUL = DETER
      RETURN
      END
C ***  GNBACK - GENERATE VALUES FOR BACKGROUND
C *
      SUBROUTINE GNBACK( PRM, YOUT )
      LOGICAL LDONE
      COMMON /XYDATA/ NPTS,X(4096),Y(4096),YC(4096),WGT(4096)
      COMMON /REFPAR/ NVECS,BVECS(10),BMINS(10),BMAXS(10)
      DIMENSION PRM(10), YOUT(4096)
C *
C *  P1 - OFFSET                       P2 - SLOPE FOR LINEAR SEGMENT
C *  P3 - POSITION OF PVII PROFILE     P4 - INTENSITY OF PVII
C *  P5 - FWHM OF PVII                 P6 - EXPONENT OF PVII
C *
      P1 =  PRM(1)
      P2 =  PRM(2)
      P3 =  PRM(3)
      P4 =  PRM(4)
      P5 =  PRM(5)
      P6 =  PRM(6)
C *
      IF(NVECS.EQ.1) THEN
        DO 1000 I=1,NPTS
1000    YOUT(I) = P1
       ELSE IF(NVECS.EQ.2) THEN
        DO 1020 I=1,NPTS
1020    YOUT(I) = P2*X(I) + P1
       ELSE
        TMP = ( 2.0**(1.0/P6) - 1.0 ) / ( (P5/2.0)**2 )
        POW = -P6
        LDONE = .FALSE.
        DO 1040 I=1,NPTS
        YOUT(I) = P2*X(I) + P1
        IF(LDONE) GO TO 1040
        XDIFF = X(I) - P3
        PVII = P4 * ( 1.0 + TMP * (XDIFF*XDIFF) )**POW
        YOUT(I) = YOUT(I) + PVII
        IF( XDIFF.GT.0.0 .AND. PVII.LT.0.01 ) LDONE = .TRUE.
1040    CONTINUE
      ENDIF
C *
      RETURN
      END
