SUBROUTINE OZPPF1(Y,X,N,C1,C2,C3,C4,F,Q0,Q1,TIME,FTIME,
     *          EXPFIT,STATUS)
*
* Module Number: 13.6.1
*
* Module Name: OZPPF1
*
* 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: OZPPF1
*
* Keywords of accessed files and tables:
* --------------------------------------
*       NONE
*
*
* Subroutines Called:
* -------------------
* CDBS:
*       NONE
* SDAS:
*       umsput
*
* History:
* --------
* version         date            Author          Description
*    1          4/20/87         D. Lindler      Designed and coded
*    2          JAN 88          D. Lindler      New SDAS I/O and standards
*    3          10/20/88	D. Lindler	Renamed from ZPPF1 to OZPPF1
*-------------------------------------------------------------------------
*
* Input parameters
*
*       Y       - vector of observed count rates (real*4)
*       X       - vector of expected count rates (real*4)
*       N       - number of points in X and Y   (integer)
*       C1,C2   - count rate range for computing q0 (real)
*       C3,C4   - count rate range for computing q1, F (real)
*       F       - linear time constant threshold (real*8)
*
* Output parameters
*
*       Q0, Q1, - paired pulse coef. (real)
*       TIME    - vector of time coef. for each data point (real)
*       FTIME   - vector of fitted time coef. (real)
*       EXPFIT  - vector of fitted expected values (real)
*       STATUS  - error status (integer)
*
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
        REAL X(1),Y(1),TIME(1),FTIME(1),EXPFIT(1)
        DOUBLE PRECISION C1,C2,C3,C4,Q0,Q1,F
        INTEGER N
        INTEGER STATUS
C
C LOCAL VARIABLES
C
        INTEGER I,M,ISTAT
        CHARACTER*130   CONTXT
        DOUBLE PRECISION SUM,SUMYT,SUMYY,Y1,XY,DEL
C
C------------------------------------------------------------------------
C
C COMPUTE TIME COEF. FOR EACH POINT IN X AND Y
C
        DO 10 I=1,N
                XY=X(I)*Y(I)
                TIME(I)=0.0
                IF(XY.GT.0.0) TIME(I)=(X(I)-Y(I))/XY
10      CONTINUE
C
C COMPUTE AVERAGE TIME COEF FOR POINTS BETWEEN C1 AND C2 (I.E. Q0)
C
        M=0
        SUM=0.0
        DO 20 I=1,N
                IF((Y(I).GT.C1).AND.(Y(I).LT.C2))THEN
                        SUM=SUM+TIME(I)
                        M=M+1
                ENDIF
20      CONTINUE
        IF(M.LT.1)THEN
                CONTXT='NO OBSERVED COUNT RATES BETWEEN C1 AND C2'
                GO TO 999
        ENDIF
        Q0=SUM/M
C
C PERFORM LINEAR FIT TO POINTS BETWEEN C3 AND C4
C
C   SUBTRACT Q0 FROM TIME AND F FROM X TO MOVE DATA TO BE FIT
C       TO AN ORIGIN OF (0,0).  The fitted line must pass through
C       This origin.  That is why lista is selected only to fit
C       the second parameter.
C
        SUMYT=0.0
        SUMYY=0.0
        DO 70 I=1,N
            IF((Y(I).GT.C3).AND.(Y(I).LT.C4))THEN
                Y1=(Y(I)-F)
                SUMYT=SUMYT+(TIME(I)-Q0)*Y1
                SUMYY=SUMYY+Y1*Y1
            ENDIF
70      CONTINUE
        IF(SUMYY.GT.0)THEN
                Q1=SUMYT/SUMYY
           ELSE
                Q1=0.0D0
                CALL UMSPUT('WARNING: No points in range c3 to c4, '//
     *                  'Q1 set to 0.0',STDOUT,0,ISTAT)
        ENDIF
C
C COMPUTE FITTED TIME AND FITTED OBSERVED COUNTS
C
        DO 100 I=1,N
                FTIME(I)=Q0
                IF(Y(I).GT.F)FTIME(I)=Q0+(Y(I)-F)*Q1
                DEL=1-FTIME(I)*Y(I)
                IF(DEL.LT.1.0E-10)DEL=1.0E-10
                EXPFIT(I)=Y(I)/DEL
100     CONTINUE
C
C DONE
C
        STATUS=0
        GO TO 1000
999     STATUS=1
        CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
1000    RETURN
        END