SUBROUTINE OZPPF
*
* Module Number: 13.6
*
* Module Name: OZPPF
*
* Keyphrase:
* ----------
*       Compute HRS paired pulse coefficients
*
* Description:
* ------------
*       The paired pulse equation is given by:
*
*               y = x/(1-tx)
*       where
*               x is the true count rate
*               y is the observed count rate
*               t is a time coef. which depends on the observed count
*                       rate by the following relation.
*
*                       t = q0  for y <= F
*                       t = q0 + q1*(x-F) for y>F
*
*       This routine computes the values for q0, q1 and F as follows.
*
*       1) The time constant for each pair of input x(i), y(i) is computed
*               by      t(i)=(x-y)/xy
*       2) The constant q0 is computed by the average of t(i) where
*               c1<y(i)<c2.  c1 and c2 are count rates supplied by
*               the user.
*       3) q1 is computed as the slope of the line fit through t(i)
*               versus y(i) for all points where c3<y(i)<c4.  c3 and
*               c4 are user supplied coefficients.
*
* Fortran Name: OZPPF
*
* Keywords of accessed files and tables:
* --------------------------------------
*       OBSERVED        input   Observed count rate file
*       EXPECTED        input   Expected count rate file
*       TABLE           output  table of fit results
*       PPTAB           output  table of paired pulse coef. for RSDP
*       C1              input   parameter giving first count rate in the
*                               region to average for computing q0
*       C2              input   paramter giving the last count rate in
*                               the region to avergage for computing q0
*       C3,C4           input   The count rates specifiing the region to
*                               fit for computing q1.
*       F               input   Threshold for linear time function
*       EPS             input   Allowed error copied to output table
*       THRESHOLD       input   Threshold copied to output table
*       ITER            input   iterations copied to output table
*
*
* Subroutines Called:
* -------------------
* CDBS:
*       OZPPF1
* SDAS:
*       prolog, epilog, getpar, copyhead, opentbl, putrow, getdat,
*       puterr
*
* History:
* --------
* version         date            Author          Description
*    1          4/20/87         D. Lindler      Designed and coded
*    2          10/20/88	D. Lindler	Renamed to OZPPF
*-------------------------------------------------------------------------
C     INCLUDE FILE FOR THE IRAF77 FORTRAN INTERFACE TO THE IRAF VOS
C
C
C     FILE I/O ACCESS MODES
C
      INTEGER   RDONLY
      PARAMETER (RDONLY = 1)
      INTEGER   RDWRIT
      PARAMETER (RDWRIT = 2)
      INTEGER   WRONLY
      PARAMETER (WRONLY = 3)
      INTEGER   APPEND
      PARAMETER (APPEND = 4)
C
C     CODES FOR DATA TYPES
C
      INTEGER   TYBOOL
      PARAMETER (TYBOOL = 1)
      INTEGER   TYCHAR
      PARAMETER (TYCHAR = 2)
      INTEGER   TYINT
      PARAMETER (TYINT = 4)
      INTEGER   TYREAL
      PARAMETER (TYREAL = 6)
      INTEGER   TYDOUB
      PARAMETER (TYDOUB = 7)
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C     UHDAS HEADER PARM TYPES -- CB, DAO, 5-SEP-87
C
      INTEGER GENHDR
      PARAMETER (GENHDR = 0)
      INTEGER IMSPEC
      PARAMETER (IMSPEC = 1)
C
C     THIS SECTION IS FOR PARAMETERS RELEVANT TO TABLE I/O.
C
C     THESE MAY BE SET BY UTPPTI AND/OR READ BY UTPGTI:
C
C                                       LENGTH OF ROW (UNIT = SIZE OF REAL)
      INTEGER   TBRLEN
      PARAMETER (TBRLEN = 1)
C                                       INCREASE ROW LENGTH
      INTEGER   TBIRLN
      PARAMETER (TBIRLN = 2)
C                                       NUMBER OF ROWS TO ALLOCATE
      INTEGER   TBALLR
      PARAMETER (TBALLR = 3)
C                                       INCREASE ALLOC NUM OF ROWS
      INTEGER   TBIALR
      PARAMETER (TBIALR = 4)
C                                       WHICH TYPE OF TABLE? (ROW OR COLUMN)
      INTEGER   TBWTYP
      PARAMETER (TBWTYP = 5)
C                                       MAXIMUM NUMBER OF USER PARAMETERS
      INTEGER   TBMXPR
      PARAMETER (TBMXPR = 6)
C                                       MAXIMUM NUMBER OF COLUMNS
      INTEGER   TBMXCL
      PARAMETER (TBMXCL = 7)
C                                       TYPE = ROW-ORDERED TABLE
      INTEGER   TBTYPR
      PARAMETER (TBTYPR = 11)
C                                       TYPE = COLUMN-ORDERED TABLE
      INTEGER   TBTYPC
      PARAMETER (TBTYPC = 12)
C
C     THESE MAY BE READ BY UTPGTI BUT MAY NOT BE SET:
C
C                                       NUMBER OF ROWS WRITTEN TO
      INTEGER   TBNROW
      PARAMETER (TBNROW = 21)
C
C     END IRAF77.INC
C
C ERROR PROCESSING PARAMETERS
C
        CHARACTER*130   CONTXT
        INTEGER STATUS,ISTAT,ISTATS(20)
C
C OUTPUT PAIRED PULSE COEF. TABLE
C
        CHARACTER*64 PPTAB
        INTEGER IDOUT,COLIDS(9),CTYPE(9),NROWS
        CHARACTER*12 COLNAM(9),CFORM(9),CUNITS(9)
C
C OUTPUT RESULT TABLE
C
        CHARACTER*64 TABLE
        CHARACTER*15 CNAME2(5)
        INTEGER CTYPE2(5)
C
C INPUT FILE I/O
C
        CHARACTER*130 OBSER,EXPECT,NAME,NAME2
        INTEGER IDTEMP,IDTMP2,IDIN,IDIN2,NAXIS,DTYPE,DIMEN(8)
        INTEGER NS,NPOINT
C
C INPUT KEYWORD PARAMETERS
C
        DOUBLE PRECISION C1,C2,C3,C4,THRESH,F
        INTEGER         IDET,DET1,DET2
C
C FIT PARAMETERS
C
        REAL    OBS(50000),EX(50000),TIME(50000),
     *                  FTIME(50000),EXPFIT(50000)
        DOUBLE PRECISION Q0,Q1
C
C OTHER LOCAL VARIABLES
C
        INTEGER I
C
C  DATA DECLARATIONS
C
        DATA COLNAM/'INSTRUMENT','DETECTOR','TAU1','EPSILON','THRESHOLD'
     *          ,'Q0','Q1','F','ITERATIONS'/
        DATA CTYPE/-3,TYINT,6*TYDOUB,TYINT/
        DATA CNAME2/'OBSERVED','EXPECTED','TIME','TIME_FIT',
     *                  'EXPECTED_FIT'/
        DATA CTYPE2/5*TYREAL/
C--------------------------------------------------------------------------
C
C READ CL PARAMETERS
C
        CALL UCLGST('observed',OBSER,ISTATS(1))
        CALL UCLGST('expected',EXPECT,ISTATS(2))
        CALL UCLGST('table',TABLE,ISTATS(3))
        CALL UCLGST('pptab',PPTAB,ISTATS(4))
        CALL UCLGSD('f',F,ISTATS(5))
        CALL UCLGSD('c1',C1,ISTATS(6))
        CALL UCLGSD('c2',C2,ISTATS(7))
        CALL UCLGSD('c3',C3,ISTATS(8))
        CALL UCLGSD('c4',C4,ISTATS(9))
        CALL UCLGSD('threshold',THRESH,ISTATS(10))
        CALL UCLGSI('detector',IDET,ISTATS(11))
        DO 10 I=1,11
                IF(ISTATS(I).NE.0)THEN
                        CONTXT='Error reading CL parameter'
                        GO TO 999
                ENDIF
10      CONTINUE
C
C READ INPUT DATA ---------------------------------------------------------
C
        NPOINT=0
C                                    --->TOTAL NUMBER OF DATA POINTS READ
C
C OPEN TEMPLATES
C
        CALL UIMOTP(OBSER,IDTEMP,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error opening observed rate file template '//
     *                                  OBSER
                GO TO 999
        ENDIF
        CALL UIMOTP(EXPECT,IDTMP2,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error opening expected rate file template '//
     *                          EXPECT
                GO TO 999
        ENDIF
C
C GET NEXT FILE NAME
C
20      CALL UIMXTP(IDTEMP,NAME,ISTAT)
        IF(ISTAT.LT.0)GO TO 100
        IF(ISTAT.NE.0)THEN
                CONTXT='Error getting filename from template '//OBSER
                GO TO 999
        ENDIF
C
C GET NEXT EXPECTED COUNT RATE FILE NAME
C
        CALL UIMXTP(IDTMP2,NAME2,ISTAT)
        IF(ISTAT.LT.0)THEN
                CONTXT='Insufficent expected rate files specified'
                GO TO 999
        ENDIF
        IF(ISTAT.NE.0)THEN
                CONTXT='Error getting filename from template '//EXPECT
                GO TO 999
        ENDIF
C
C OPEN INPUT FILES
C
        CALL UIMOPN(NAME,RDONLY,IDIN,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error opening input file '//NAME
                GO TO 999
        ENDIF
        CALL UIMOPN(NAME2,RDONLY,IDIN2,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error opening input file '//NAME2
                GO TO 999
        ENDIF
C
C READ IMAGE INFO
C
        CALL UIMGID(IDIN,DTYPE,NAXIS,DIMEN,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error reading input file '//NAME
                GO TO 999
        ENDIF
C
C CHECK FOR VALID DATA
C
        IF(NAXIS.NE.1)THEN
                CONTXT='Input data must be one dimensional'
                GO TO 999
        ENDIF
        NS=DIMEN(1)
C
C CHECK TO SEE IF WE STILL HAVE ROOM
C
        IF((NS+NPOINT).GT.50000)THEN
                CONTXT='Total number of input data points can not '//
     *                  'exceed 50000'
                GO TO 999
        ENDIF
C
C READ EXPECTED COUNT RATE FILE FILE INFO
C
        CALL UIMGID(IDIN2,DTYPE,NAXIS,DIMEN,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error reading file '//NAME2
                GO TO 999
        ENDIF
C
C CHECK FOR VALID DATA
C
        IF((NAXIS.NE.1).OR.(DIMEN(1).NE.NS))THEN
                CONTXT='Expected count rate file '//
     *                  'must be same size as observed'
                GO TO 999
        ENDIF
C
C READ DATA
C
        CALL UIGL1R(IDIN,OBS(NPOINT),ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error reading input file '//NAME
                GO TO 999
        ENDIF
        CALL UIGL1R(IDIN2,EX(NPOINT),ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error reading input file '//NAME2
                GO TO 999
        ENDIF
        NPOINT=NPOINT+NS
C
C CLOSE IMAGES
C
        CALL UIMCLO(IDIN,ISTAT)
        CALL UIMCLO(IDIN2,ISTAT)
C
C GO GET NEXT IMAGE
C
        GO TO 20
100     CONTINUE
C
C DONE READING INPUT DATA ------------------------------------------------
        IF(NPOINT.EQ.0)THEN
               CONTXT='No input files found'
               GO TO 999
        ENDIF
C
C
C COMPUTE PAIRED PULSE COEF. AND FITS
C
        CALL OZPPF1(OBS,EX,NPOINT,C1,C2,C3,C4,F,
     *                  Q0,Q1,TIME,FTIME,EXPFIT,STATUS)
        IF(STATUS.NE.0)THEN
                CONTXT='No output files generated'
                GO TO 999
        ENDIF
C
C WRITE FIT RESULT TABLE ---------------------------------------------------
C
C
C OPEN OUTPUT TABLE
C
        CALL UTTINN(TABLE,IDOUT,ISTATS(1))
        CALL UTPPTI(IDOUT,TBWTYP,TBTYPC,ISTATS(2))
        CALL UTPPTI(IDOUT,TBALLR,NPOINT,ISTATS(3))
        CALL UTCDEF(IDOUT,CNAME2,CUNITS,CFORM,CTYPE2,5,COLIDS,ISTATS(4))
        CALL UTTCRE(IDOUT,ISTATS(5))
        DO 200 I=1,5
                IF(ISTATS(I).NE.0)THEN
                        CONTXT='Error creating output table '//TABLE
                        GO TO 999
                ENDIF
200     CONTINUE
C
C COPY RESULTS TO TABLE
C
        CALL UTCPTR(IDOUT,COLIDS(1),1,NPOINT,OBS,ISTATS(1))
        CALL UTCPTR(IDOUT,COLIDS(2),1,NPOINT,EX,ISTATS(2))
        CALL UTCPTR(IDOUT,COLIDS(3),1,NPOINT,TIME,ISTATS(3))
        CALL UTCPTR(IDOUT,COLIDS(4),1,NPOINT,FTIME,ISTATS(4))
        CALL UTCPTR(IDOUT,COLIDS(5),1,NPOINT,EXPFIT,ISTATS(5))
        DO 210 I=1,3
                IF(ISTATS(I).NE.0)THEN
                        CONTXT='Error writing to output table '//TABLE
                        GO TO 999
                ENDIF
210     CONTINUE
        CALL UTTCLO(IDOUT,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error closing output table '//TABLE
                GO TO 999
        ENDIF
C
C WRITE PAIRED PULSE COEFFICIENT TABLE ------------------------------------
C
C
C OPEN OUTPUT TABLE
C
        CALL UTTINN(PPTAB,IDOUT,ISTATS(1))
        CALL UTPPTI(IDOUT,TBRLEN,15,ISTATS(2))
        CALL UTPPTI(IDOUT,TBMXCL,9,ISTATS(3))
        CALL UTCDEF(IDOUT,COLNAM,CUNITS,CFORM,CTYPE,9,COLIDS,ISTATS(4))
        CALL UTTCRE(IDOUT,ISTATS(5))
        DO 300 I=1,5
                IF(ISTATS(I).NE.0)THEN
                        CONTXT='Error creating output table '//PPTAB
                        GO TO 999
                ENDIF
300     CONTINUE
C
C COPY RESULTS TO TABLE, ONE ROW FOR EACH DETECTOR IF IDET=0
C       OTHERWISE ONE ROW FOR THE SEPCIFIED IDET
C
        IF(IDET.EQ.0)THEN
                DET1=1
                DET2=2
           ELSE
                DET1=IDET
                DET2=IDET
        ENDIF
        NROWS=0
        DO 700 IDET=DET1,DET2
            NROWS=NROWS+1
            CALL UTRPTT(IDOUT,COLIDS(1),1,NROWS,'HRS',ISTATS(1))
            CALL UTRPTI(IDOUT,COLIDS(2),1,NROWS,IDET,ISTATS(2))
            CALL UTRPTD(IDOUT,COLIDS(3),1,NROWS,0.0D0,ISTATS(3))
            CALL UTRPTD(IDOUT,COLIDS(4),1,NROWS,0.0D0,ISTATS(4))
            CALL UTRPTD(IDOUT,COLIDS(5),1,NROWS,THRESH,ISTATS(5))
            CALL UTRPTD(IDOUT,COLIDS(6),1,NROWS,Q0,ISTATS(6))
            CALL UTRPTD(IDOUT,COLIDS(7),1,NROWS,Q1,ISTATS(7))
            CALL UTRPTD(IDOUT,COLIDS(8),1,NROWS,F,ISTATS(8))
            CALL UTRPTI(IDOUT,COLIDS(9),1,NROWS,0,ISTATS(9))
            DO 310 I=1,9
                IF(ISTATS(I).NE.0)THEN
                        CONTXT='Error writing to output table'//PPTAB
                        GO TO 999
                ENDIF
310         CONTINUE
700     CONTINUE
        CALL UTTCLO(IDOUT,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='Error closing output table '//TABLE
                GO TO 999
        ENDIF
C
C DONE
C
        GO TO 1000
999     CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
1000    RETURN
        END