C     ****************************************************************
C     *                                                              *
C     *                    FORTRAN-77 PROGRAM                        *
C     *                                                              *
C     *   SENSITIVITY ANALYSIS UNDER PARAMETERIZATION OF P0 AND P1   *
C     *                                                              *
C     *                      SOURCE PROGRAM                          *
C     *                                                              *
C     *                       MARCH 4, 1997                          *
C     *                                                              *
C     ****************************************************************
C
C
C
C     This FORTRAN program implements the method of sensitivity analysis
C     of Lin, Psaty and Kronmal under the parameterization of P1 and P0.
C
C     'sen.dat' and 'sen.out' are the input and output file names.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION G(100),P1(100),P0(100),A(100)
C
      OPEN (5,FILE='sen.dat')
      OPEN (6,FILE='sen.out')
C
C     In the input file:
C
C     RP, RL, RU are the point estimate, and lower and upper confidence 
C       limits.
C     Gamma is the odds ratio or hazard ratio of disease associated with
C       the unmeasured confounder.
C     P1 is the prevalence of the unmeasured confounder among the exposed.
C     P0 is the prevalence of the unmeasured confounder among the unexposed.
C     NG,NP1 and NP0 are the numbers of values for Gamma, P1 and P0.
C     G(I), P1(I) and P0(I) are the ith values of G, P1 and P0. 
C
C     The format of the output is similar to Tables 2 and 3 of the paper.
C
      READ (5,*) RP,RL,RU
      READ (5,*) NG,NP1,NP0
      READ (5,*) (G(K),K=1,NG)
      READ (5,*) (P1(I),I=1,NP1)
      READ (5,*) (P0(I),I=1,NP0)
C
      DO 50 K=1,NG
      WRITE (6,900) G(K),(P0(J),J=1,NP0)
      WRITE (6,990)
      DO 20 I=1,NP1
         DO 10 J=1,NP0
            A(J)=(1.D0+(G(K)-1.D0)*P1(I))/(1.D0+(G(K)-1.D0)*P0(J))
   10    CONTINUE
         WRITE (6,900) P1(I),(RP/A(J),J=1,NP0)
         WRITE (6,920) (RL/A(J),RU/A(J),J=1,NP0)
   20 CONTINUE
      WRITE (6,990)
   50 CONTINUE
C
  900 FORMAT (F6.2,20F14.2)
  920 FORMAT (10X,'(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5..2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2,')',1X,
     &            '(',F5.2,',',F5.2,')',1X,'(',F5.2,',',F5.2)
  990 FORMAT (' ')
C
      STOP
      END