       PROGRAM VRIPF

* Blame for this program resides with Peter A. Bergbusch, Dept. of Physics,
* University of Regina. 

* Last update: Sept. 28, 2005.

*23456789012345678901234567890123456789012345678901234567890123456789012
       PARAMETER (NMX = 1000, NBN = 1000)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DOUBLE PRECISION LUM(NMX),MASS(NMX),MAG(NMX),IPFW,LFNW,CFNW,MV,
     1                  IPFTH(NMX),IPFTM(NBN),IPFTC(NBN),IPFTB(NBN),
     2                  IPFTT(NBN),IPFTW(NBN),IPFTG(NBN),MBOL(NMX)
       DIMENSION PHIBIN(NBN,3),BINLIM(NBN+1),BINDTH(NBN+1),rnorm(3),
     1           BINCEN(NBN),CSCOEF(4,NMX),CUMPHI(NMX,3),CSC(4,NMX),
     2           PHI(NMX,3),TEFF(NMX),GRAV(NMX),CI(NMX),DTH(NMX),
     3           AGE(10),
     4           DOBS(NMX),X(3),POWER(3),
     5           TM(NMX),TC(NMX),TP(NMX),TD(NMX),DIPF(NBN,3),
     6           CIPF(NBN,3),cumgb(3)
       CHARACTER INFILE*60,OUTFILE*60,COLIND(5)*3,f8*63,
     1           HLINE(6)*65,ALINE*11,BAND(4)*1
       LOGICAL  YES,IPF,LFN,CFN,MORE

       PARAMETER (C1 = 10.60917D0)

       DATA COLIND/'B-V','B-R','B-I','V-R','V-I'/
       DATA BAND/'B','V','R','I'/
       DATA IPF/.FALSE./, LFN/.FALSE./, CFN/.FALSE./


       CALL IOFILE(7,'IN','FILE.ISO',.FALSE.,' ',
     1               INFILE,.TRUE.,'          Input ISOCHRONE File: ')

* Other INPUT files are BVRILO.DAT on unit 10 and BVRIHI.DAT on unit 11. 
* They are opened and closed inside subroutine COLOR.

       DO I = 1, 7
         READ (7,'(A)') HLINE(I)
       END DO


       WRITE(*,'(//T33,''Available Magnitudes & Colour Indices'')')
       WRITE(*,'(/T40,'' 1) B          1) B-V'')')
       WRITE(*,'( T40,'' 2) V          2) B-R'')')
       wRITE(*,'( T40,'' 3) R          3) B-I'')')
       WRITE(*,'( T40,'' 4) I          4) V-R'')')
       WRITE(*,'( T40,''               5) V-I''//)')
       WRITE(*,'(T11,A,$)')
     1       ' Select a Magnitude(1-4) & Colour Index (1-5): '
       READ(*,*) NMG, NCI

       WRITE(*,'(//t11,A,$)')
     1  ' Power law mass spectrum exponent (1+x) '
       WRITE(*,'(/t11,A,$)')
     1  '               Enter three values of x: '
       READ(*,*) (X(I), I=1,3)
       DO I = 1, 3
         POWER(I) = 1.0D0 + X(I)
       END DO
       WRITE(*,'(/)')

* The OUTPUT files which can be opened have the file extensions .IPF on
* unit 8, .LFN on unit 12, and .CFN on unit 14

       IF (YES('Do you want an IPF?')) THEN
         IPF = .TRUE.
         WRITE(*,'(t25,A,$)') 'Input IPF bin width: '
         READ(*,*) IPFW
         CALL PARSE (INFILE,'NAME',OUTFILE,NCH)
         NCHOUT = NCH + 4
         OUTFILE(1:NCHOUT) = OUTFILE(1:NCH)//'.ipf'
         CALL IOFILE(8,'OUT','FILE.IPF',.FALSE.,' ',
     1                  OUTFILE,.FALSE.,' OUTPUT File: ')
         HLINE(1)(2:11)='IPFs      '
         WRITE(8,'(A,2X,A,2X,A,2x,''x ='',sp3f5.1)')
     1             HLINE(1)(1:14),BAND(NMG),COLIND(NCI),X
         DO I = 2, 7
           CALL STR_TRIM(65,HLINE(I),NCH)
           WRITE(8,'(A)') HLINE(I)(1:nch)
         END DO
         DO I = 2, 7
           CALL STR_TRIM(65,HLINE(I),NCH)
         END DO
       END IF

       IF (YES(' Do you want a LF?')) THEN
         LFN = .TRUE.
         WRITE(*,'(t25,A,$)') ' Input LF bin width: '
         READ(*,*) LFNW
         CALL PARSE (INFILE,'NAME',OUTFILE,NCH)
         NCHOUT = NCH + 4
         OUTFILE(1:NCHOUT) = OUTFILE(1:NCH)//'.lfn'
         CALL IOFILE(12,'OUT','FILE.LFN',.FALSE.,' ',
     1                   OUTFILE,.FALSE.,' OUTPUT File: ')
         HLINE(1)(2:11)='LFs       '
         WRITE(12,'(A,2X,A,2X,A,2x,''x ='',sp3f5.1)') 
     1              HLINE(1)(1:14),BAND(NMG),COLIND(NCI),X
         DO I = 2, 7
           CALL STR_TRIM(65,HLINE(I),NCH)
           WRITE(12,'(A)') HLINE(I)(1:nch)
         END DO
       END IF

       IF (YES(' Do you want a CF?')) THEN
         CFN = .TRUE.
         WRITE(*,'(t25,A,$)') ' Input CF bin width: '
         READ(*,*) CFNW
         CALL PARSE (INFILE,'NAME',OUTFILE,NCH)
         NCHOUT = NCH + 4
         OUTFILE(1:NCHOUT) = OUTFILE(1:NCH)//'.cfn'
         CALL IOFILE(14,'OUT','FILE.CFN',.FALSE.,' ',
     1                   OUTFILE,.FALSE.,' OUTPUT File: ')
         HLINE(1)(2:11)='CFs       '
         WRITE(14,'(A,2X,A,2X,A,2x,''x ='',sp3f5.1)') 
     1              HLINE(1)(1:14),BAND(NMG),COLIND(NCI),X
         DO I = 2, 7
           CALL STR_TRIM(65,HLINE(I),NCH)
           WRITE(14,'(A)') HLINE(I)(1:nch)
         END DO
       END IF

       READ(HLINE(1),'(12x,I2)') NISO
       READ(HLINE(2),'(11X,F7.3)') FE

       WRITE(*,'(/,I6,'' ISOCHRONES'',t21,''[Fe/H] = '',f5.2)') 
     2        NISO,FE

       CALL BVRI(0,FE,G0,T0,BMV,VMR,VMI,BC)

*  Process one isochrone at a time

       DO NI = 1, NISO

         BLUE = 99.9

         READ (7,'(/A)') ALINE
         READ (7,'(F6.2,I6)') AGE(NI),NPTS
         WRITE(*,'(t15,'' Age, NPTS: '',F6.2,I6)') AGE(NI),NPTS
         DOBS(1) = 0.0D0

*23456789012345678901234567890123456789012345678901234567890123456789012

         DO NP = 1, NPTS

           READ(7,'(4X,F10.6,F9.6,F13.10,F11.7,D14.7)') LUM(NP),
     1          TEFF(NP),MASS(NP),DTH(NP),PHI0
           MBOL(NP) = 4.75D0 - 2.5D0*LUM(NP)
           DO J = 1, 3
             PHI(NP,J) = PHI0/MASS(NP)**POWER(J)
           END DO
           
           GRAV(NP) = LOG10(MASS(NP)) - C1 + 4.0D0*TEFF(NP) - LUM(NP)

           CALL BVRI(1,FE,GRAV(NP),TEFF(NP),BMV,VMR,VMI,BC)

           MV = MBOL(NP) - BC

           IF (NCI.EQ.1) THEN
             CI(NP) = BMV
           ELSE IF (NCI.EQ.2) THEN
             CI(NP) = BMV + VMR
           ELSE IF (NCI.EQ.3) THEN
             CI(NP) = BMV + VMI
           ELSE IF (NCI.EQ.4) THEN
             CI(NP) = VMR
           ELSE IF (NCI.EQ.5) THEN
             CI(NP) = VMI
           END IF

           IF (NMG.EQ.1) THEN
             MAG(NP) = MV + BMV
           ELSE IF (NMG.EQ.2) THEN
             MAG(NP) = MV
           ELSE IF (NMG.EQ.3) THEN
             MAG(NP) = MV - VMR
           ELSE IF (NMG.EQ.4) THEN
             MAG(NP) = MV - VMI
           END IF
*           write(25,'(i4,4f7.3)') np,mv,bmv,mv-vmi,bmv+vmi

           IF (NP.GT.1) THEN
             DOBS(NP) = DOBS(NPM1) + 
     1       DELTAMAG(CI(NPM1),CI(NP),MAG(NPM1),MAG(NP))
           END IF
           NPM1 = NP

* Take the bluest point on the isochrone as the turnoff point

           IF (CI(NP).LT.BLUE) THEN
             BLUE = CI(NP)
             ITO = NP
           END IF

         END DO
*
* THE NEXT SECTION OF CODE DEALS WITH THE CONSTRUCTION OF BOTH DIFFERENTIAL
* AND CUMULATIVE ISOCHRONE PROBABILITY FUNCTIONS.
*
         IF (IPF) THEN

* Define the reference point for IPFs at 0.05 mag to the red of the turnoff
* on the subgiant branch

           I = ITO
           
           if (ito.lt.npts) then
           BLUE05 = BLUE + 5.0D-2
           DO WHILE (CI(I).LT.BLUE05)
             I = I + 1
           END DO           
           IM1 = I - 1
           D0 = DOBS(IM1) + (BLUE05 - CI(IM1))*
     1          (DOBS(I) - DOBS(IM1))/(CI(I) - CI(IM1))
           else
           d0=dobs(ito)
           end if
           
           DO I = 1, NPTS
             DOBS(I) = DOBS(I) - D0
           END DO 

           CALL MAKE_IPFBINS(DOBS(1),DOBS(NPTS),IPFW,BINCEN,BINLIM,NBIN)

* Now correlate bin limits on the observer's plane with those on the 
* theoretical plane

           CALL AKM_CSPL (NPTS,DOBS,DTH,CSCOEF)
           NBINP1 = NBIN + 1
           DO I = 1, NBINP1
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINLIM(I),IPFTH(I),DFP)
           END DO

           CALL AKM_CSPL (NPTS,DOBS,MAG,CSCOEF)
           DO I = 1, NBIN
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINCEN(I),IPFTM(I),DFP)
           END DO

           CALL AKM_CSPL (NPTS,DOBS,CI,CSCOEF)
           DO I = 1, NBIN
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINCEN(I),IPFTC(I),DFP)
           END DO

           CALL AKM_CSPL (NPTS,DOBS,MBOL,CSCOEF)
           DO I = 1, NBIN
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINCEN(I),IPFTB(I),DFP)
           END DO

           CALL AKM_CSPL (NPTS,DOBS,TEFF,CSCOEF)
           DO I = 1, NBIN
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINCEN(I),IPFTT(I),DFP)
           END DO

           CALL AKM_CSPL (NPTS,DOBS,MASS,CSCOEF)
           DO I = 1, NBIN
             CALL AKM_EVAL (NPTS,DOBS,CSCOEF,BINCEN(I),IPFTW(I),DFP)
             IPFTG(I) = LOG10(IPFTW(I)) - C1 + 4.0D0*IPFTT(I)
     1                  + 4.0D-1*(IPFTB(I) - 4.75D0)
           END DO


*  PATCHES ******************************************************************
*  1. Patch to fix the sign on the "Distance". Since the distance is measured
*  in magnitudes, it should decrease with increasing evolutionary state
*  along an isochrone. That is, bright stars at the RGB tip should have
*  -ve distances.
*****************************************************************************

           DO I = 1, NBIN
             BINCEN(I) = -BINCEN(I)
           END DO

* Now integrate across each bin and write out the results

           DO J = 1, 3
             DO JJ = 1, NPTS
               TP(JJ) = PHI(JJ,J)
             END DO

             CALL AKM_CSPL (NPTS,DTH,TP,CSCOEF)

             DO I = 1, NBIN
               IP1 = I + 1
               CALL AKM_INT (IPFTH(I),IPFTH(IP1),NPTS,DTH,CSCOEF,
     1                       DIPFT)
               CALL AKM_INT (IPFTH(I),IPFTH(NBINP1),NPTS,DTH,CSCOEF,
     1                       CIPF(I,J))
               IF (I.LT.NBIN) THEN
                 DIPF(I,J) = DIPFT/IPFW
               ELSE
                 DIPF(I,J) = DIPFT/(DOBS(NPTS)-BINLIM(NBIN))
               END IF
             END DO
           rnorm(j) = 5.0D0 - LOG10(CIPF(1,J))
           END DO


           WRITE(8,'(/,A)') ALINE
           WRITE(8,'(F6.2,I5)') AGE(NI),NBIN

      f8 = '(9X,A,5X,A,4X,a,5x,a,2(2X,A),5X,A,1X,8X,A,spf5.1,'//
     1     '2(11x,A,F5.1))'
*2345678901234567890123456789012345678901234567890123456789012345678901234567890
           WRITE(8,f8)
     1             BAND(NMG),COLIND(NCI),'Mass','Mbol','log Te','log g',
     2             'd','x =',x(1),'x =',x(2),'x =',x(3)
*     1           '(9X,A,5X,A,4X,a,5x,a,2(2X,A),5X,A,1X,8X,A,spf5.1,
*     1             2(11x,A,F5.1))')

           DO I = 1, NBIN
             DO J=1,3
               DIPF(I,J) = rnorm(j) + LOG10(DIPF(I,J))
               CIPF(I,J) = rnorm(j) + LOG10(CIPF(I,J))
             END DO

* Write out the results to the .ipf file. 
* Each line lists the point number (I),the magnitude in the specified 
* passband (IPFTM), the colour index (IPFTC), the mass (IPFTW), the 
* bolometric magnitude (IPFTB), the effective temperature (IPFTT), 
* the surface gravity log g (IPFTG), the distance along the isochrone
* (BINCEN), and three pairs of the differential (DIPF) and cumulative
* (CIPF) isochrone population functions for the three mass spectrum
* exponents specified in the header.

             WRITE (8,'(I5,2F7.3,F10.7,F7.3,F7.4,F7.3,F8.3,
     1                  3(f10.5,F9.5))') 
     1             I,IPFTM(I),IPFTC(I),IPFTW(I),IPFTB(I),IPFTT(I),
     2             IPFTG(I),BINCEN(I),(DIPF(I,J), CIPF(I,J), J=1,3)
           END DO

* Write the results for the RGB tip

*           WRITE (8,'(5X,2F7.3,F10.7,F7.3,F7.4,F7.3,F8.3)')
*     1       MAG(NPTS),CI(NPTS),MASS(NPTS),MBOL(NPTS)+3.0D-2,
*     2       TEFF(NPTS),GRAV(NPTS),-DOBS(NPTS)

         END IF
*
*  THE NEXT SECTION OF CODE DEALS WITH THE CONSTRUCTION OF LUMINOSITY 
*  FUNCTIONS
*
         IF (LFN) THEN

*  Get the range of magnitudes to set up the bins

           BRIGHT = 9.999d1
           FAINT = -9.999d1

           DO I = 1, NPTS
             BRIGHT = MIN(BRIGHT,MAG(I))
             FAINT = MAX(FAINT,MAG(I))
           END DO

           CALL MAKE_LFBINS(BRIGHT,FAINT,LFNW,BINCEN,BINLIM,NBIN)

*  Set up the spline to interpolate distances with respect to magnitude to
*  find the bin limits in distance that correspond to the magnitude limits
*  of each bin

           CALL AKM_CSPL (NPTS,DTH,MAG,CSC)

*  Integrations are performed for three values of the mass spectrum exponent

           DO J = 1, 3

*  Zero the bins

           DO I = 1, NBIN
             PHIBIN(I,J) = 0.0D0
             CUMPHI(I,J) = 0.0D0
           END DO

*  Construct the spline in the theoretical plane

           DO JJ = 1, NPTS
             TP(JJ) = PHI(JJ,J)
           END DO

           CALL AKM_CSPL (NPTS,DTH,TP,CSCOEF)

*  Now find the first point that lies within the region spanned by the bins.

           MORE = .TRUE.
           K = 0
           DO WHILE (MORE)
             K = K + 1
             IBP = INT ((BINLIM(1) - MAG(K))/LFNW + 1.0D0)
             MORE = IBP.LT.1
           END DO

*  Interpolate in DTH wrt magnitude to establish the faint limit (DTH1) 
*  of the first bin

*23456789012345678901234567890123456789012345678901234567890123456789012

           IF (K.EQ.1) THEN

             DTH1 = DTH(1)

           ELSE

             KM1 = K - 1
             CALL AKM_BISECTION (NPTS,CSC,DTH,DTH(KM1),DTH(K),
     1                           BINLIM(1),DTH1)

             IF (IBP.GT.1) THEN

               DO JB = 2, IBP
                 JBM1 = JB-1
                CALL AKM_BISECTION (NPTS,CSC,DTH,DTH(KM1),DTH(K),
     1                              BINLIM(jB),DTH2)
                 CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
                 CUMPHI(JBM1,J) = CUMPHI(JBM1,J) + DPHI
                 DTH1 = DTH2
               END DO 

             END IF

           END IF

           K = K + 1               

*  Now scan through the rest of the points to find the other bin limits and
*  to accumulate the contents for each magnitude bin.
                    
           DO I = K, NPTS

*  Determine which bin the current point belongs in

             IB = INT ((BINLIM(1) - MAG(I))/LFNW + 1.0D0)
             IB = MIN(IB,NBIN)

*  If the current data point lies in a different bin than the last one, 
*  then interpolate the distance wrt magnitude to establish the upper 
*  limit (DTH2) for the integration of the previous bin. Then integrate 
*  over the interval (DTH1,DTH2) and accumulate the contents of that bin.

             IBDIFF = IB - IBP

             IF (abs(ibdiff).eq.1) THEN

               IF (IBDIFF.EQ.1) THEN
                 BINREF = BINLIM(IB)
                 jb = ib
               ELSE
                 BINREF = BINLIM(IBP)
                 jb = ibp
               END IF

               IM1 = I - 1
              CALL AKM_BISECTION (NPTS,CSC,DTH,DTH(IM1),
     1                            DTH(I),BINref,DTH2)
               CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
               CUMPHI(IBP,J) = CUMPHI(IBP,J) + DPHI
               IBP = IB
               DTH1 = DTH2
               if (i.eq.npts) then
               DTH2 = DTH(I)      
               CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
               CUMPHI(IBP,J) = CUMPHI(IBP,J) + DPHI
               end if

             ELSE if (abs(ibdiff).gt.1) then

               istep = sign(1,ibdiff)
               IM1 = I - 1
               DO JI = 1, IBDIFF, istep
                 JB = IBP + istep
                 CALL AKM_BISECTION (NPTS,CSC,DTH,DTH(IM1),
     1                               DTH(I),BINLIM(JB),DTH2)
                 CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
                 CUMPHI(IBP,J) = CUMPHI(IBP,J) + DPHI
                 IBP = JB
                 DTH1 = DTH2
               END DO
               if (i.eq.npts) then
               DTH2 = DTH(I)      
               CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
               CUMPHI(IBP,J) = CUMPHI(IBP,J) + DPHI
               end if

             ELSE IF (I.EQ.NPTS) THEN

               DTH2 = DTH(I)      
               CALL AKM_INT (DTH1,DTH2,NPTS,DTH,CSCOEF,DPHI)
               CUMPHI(IBP,J) = CUMPHI(IBP,J) + DPHI

             END IF

           END DO

*23456789012345678901234567890123456789012345678901234567890123456789012
           do i=1,nbin-1
              phibin(i,j)=cumphi(i,j)/lfnw
           end do
           phibin(nbin,j)=cumphi(nbin,j)/(binlim(nbin)-bright)

           DO I = NBIN-1, 1, -1
             IP1 = I + 1
             CUMPHI(I,J) = CUMPHI(IP1,J) + CUMPHI(I,J)
           END DO
           rnorm(j) = 5.0D0 - LOG10(Cumphi(1,J))

         END DO

           WRITE(12,'(/,A)') ALINE
           WRITE(12,'(F6.2,I5)') AGE(NI),NBIN
           WRITE(12,'(9X,A,9X,A,SPF5.1,2(11X,A,F5.1))') BAND(NMG),
     1         'x =',x(1),'x =',x(2),'x =',x(3)

           DO I = 1, NBIN
             DO J = 1, 3
               PHIBIN(I,J) = rnorm(j) + LOG10(PHIBIN(I,J))
               CUMPHI(I,J) = rnorm(j) + LOG10(CUMPHI(I,J))
             END DO
             WRITE (12,'(I5,F7.3,3(f10.5,f9.5))') I,BINCEN(I),
     1                    (PHIBIN(I,J), CUMPHI(I,J),J = 1, 3)
           END DO

         END IF
*
*  THE NEXT SECTION OF CODE DEALS WITH THE CONSTRUCTION OF DIFFERENTIAL COLOR
*  FUNCTIONS FOR THE TURNOFF-TO-RGBTIP PHASES OF EVOLUTION AND WITH THE 
*  EVALUATION OF THE SGB COLOR FUNCTION RATIOS.
*
         IF (CFN) THEN

* Construct the differential color function bins first


           BLUE = CI(ITO)
           RED = CI(1)

           call Make_cfbins(.true.,blue,red,cfnw,nms,bincen,binlim)

           RED = CI(NPTS)

           CALL MAKE_CFBINS(.TRUE.,BLUE,RED,CFNW,NBIN,BINCEN,BINLIM)
           NBINP1 = NBIN + 1

*  Correlate the color bin limits with distance on the IPF

           J = 0
           DO I = ITO, NPTS
             J = J + 1 
             TC(J) = CI(I)
             TD(J) = DTH(I)
             TM(J) = MAG(I)
           END DO
           NPPTO = J

           CALL AKM_CSPL(NPPTO,TC,TD,CSCOEF)

           DO I = 1, NBINP1
             CALL AKM_EVAL (NPPTO,TC,CSCOEF,BINLIM(I),BINDTH(I),DFP)
           END DO

*  Zero the bins

           DO J = 1, 3
           DO I = 1, NBIN
             PHIBIN(I,J) = 0.0D0
           END DO

*  Set up the spline for integration

           K = 0
           DO I = ITO, NPTS
             K = K + 1
             TP(K) = PHI(I,J)
           END DO

           CALL AKM_CSPL (NPPTO,TD,TP,CSCOEF)

*  Integrate the spline across each bin
           DO I = 1, NBIN
             IP1 = I + 1
             CALL AKM_INT (BINDTH(I),BINDTH(IP1),NPPTO,TD,CSCOEF,
     1                    PHIBIN(I,J))
             CALL AKM_INT (BINDTH(I),BINDTH(NBINp1),NPPTO,TD,CSCOEF,
     1                    CUMPHI(I,J))
           END DO
             cumgb(j) = cumphi(1,j)
           END DO

           NBINM1 = NBIN - 1

           do j = 1, 3
             PHIBIN(1,J) = LOG10(PHIBIN(1,J)/(binlim(2) - blue))
             CUMPHI(1,J) = LOG10(CUMPHI(1,J))

           DO I = 2, NBINM1
             PHIBIN(I,J) = LOG10(PHIBIN(I,J)/CFNW)
             CUMPHI(I,J) = LOG10(CUMPHI(I,J))
           END DO

           CUMPHI(NBIN,J) = LOG10(PHIBIN(NBIN,J))
           PHIBIN(NBIN,J) = LOG10(PHIBIN(NBIN,J)/(RED-BINLIM(NBIN)))
           end do

           ntot = nms + nbin
           k = ntot
           do i = nbin, 1, -1
             bincen(k) = bincen(i)
             do j = 1, 3
               phibin(k,j) = phibin(i,j)
               cumphi(k,j) = cumphi(i,j)
             end do
           k = k - 1
           end do

* Construct the differential color function bins first

           RED = CI(1)

           CALL MAKE_CFBINS(.TRUE.,BLUE,RED,CFNW,NBIN,BINCEN,BINLIM)
           NBINP1 = NBIN + 1

*  Correlate the color bin limits with distance on the IPF

           J = 0
           DO I = ITO, 1, -1
             J = J + 1 
             TC(J) = CI(I)
             TD(J) = DTH(I)
             TM(J) = MAG(I)
           END DO
           NPPTO = J

           CALL AKM_CSPL(NPPTO,TC,TD,CSCOEF)

           DO I = 1, NBINP1
             CALL AKM_EVAL (NPPTO,TC,CSCOEF,BINLIM(I),BINDTH(I),DFP)
           END DO
           bindth(nbinp1) = dth(1)

*  Zero the bins

           DO J = 1, 3
           DO I = 1, NBIN
             PHIBIN(I,J) = 0.0D0
           END DO

*  Set up the spline for integration

           K = 0
           DO I = ITO, 1, -1
             K = K + 1
             TP(K) = PHI(I,J)
           END DO

           CALL AKM_CSPL (NPPTO,TD,TP,CSCOEF)

*  Integrate the spline across each bin
           DO I = 1 ,NBIN
             IP1 = I + 1
             CALL AKM_INT (BINDTH(I),BINDTH(IP1),NPPTO,TD,CSCOEF,
     1                    PHIBIN(I,J))
             CALL AKM_INT (BINDTH(1),BINDTH(ip1),NPPTO,TD,CSCOEF,
     1                    CUMPHI(I,J))
           END DO
           rnorm(j) = 5.0D0 - log10(Cumphi(nbin,J) + cumgb(j))
           END DO

           NBINM1 = NBIN - 1

           do j = 1, 3
             PHIBIN(1,J) = LOG10(PHIBIN(1,J)/(binlim(2) - blue))
             CUMPHI(1,J) = LOG10(CUMPHI(1,J)+cumgb(j))

           DO I = 2, NBINM1
             PHIBIN(I,J) = LOG10(PHIBIN(I,J)/CFNW)
             CUMPHI(I,J) = LOG10(CUMPHI(I,J)+cumgb(j))
           END DO

           CUMPHI(NBIN,J) = LOG10(cumPHI(NBIN,J)+cumgb(j))
           PHIBIN(NBIN,J) = LOG10(PHIBIN(NBIN,J)/(RED-BINLIM(NBIN)))
           end do

           do j = 1, 3
             do i = 1, ntot
               phibin(i,j) = rnorm(j) + phibin(i,j)
               cumphi(i,j) = rnorm(j) + cumphi(i,j)
             end do
           end do

           WRITE(14,'(/,A,a)') ALINE,'   Nms'
           WRITE(14,'(F6.2,I5,I6)') AGE(NI),Ntot, nms
           WRITE(14,'(8X,A,9X,A,SPF5.1,2(11X,A,F5.1))') COLIND(NCI),
     1         'x =',x(1),'x =',x(2),'x =',x(3)
           k = 0
           do i = nbin, 1, -1
             k = k + 1
             WRITE(14,'(I5,F7.3,3(f10.5,F9.5))') k,BINCEN(I),
     1                    (PHIBIN(i,J),CUMPHI(i,J), J=1,3)
           end do
           do i = nms+1, ntot
             WRITE(14,'(I5,F7.3,3(f10.5,F9.5))') I,BINCEN(I),
     1                    (PHIBIN(I,J),CUMPHI(I,J), J=1,3)
           end do

     
         END IF

       END DO

       CLOSE (UNIT=7,STATUS='KEEP')
       IF (IPF) CLOSE (UNIT=8,STATUS='KEEP')
       IF (LFN) CLOSE (UNIT=12,STATUS='KEEP')
       IF (CFN) CLOSE (UNIT=14,STATUS='KEEP')

       STOP
       END

*23456789012345678901234567890123456789012345678901234567890123456789012

       SUBROUTINE MAKE_IPFBINS(D1,D2,WIDTH,BINCEN,BINLIM,NBIN)

*      1998-May-21
*      This subroutine sets up the binning in "distance" along an isochrone.
*      The distance is defined with respect to some well-defined point on
*      the isochrone, such as a point 0.05 mag to the red of the turnoff
*      point on the subgiant branch.

       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION BINCEN(*),BINLIM(*)

       HWIDTH = WIDTH/2.0d0

       NB1 = INT(ABS(D1 + HWIDTH)/WIDTH)
       NB2 = INT(ABS(D2 - HWIDTH)/WIDTH)
       K = NB1
       DO J = 1, NB1
         BINCEN(J) = -K*WIDTH
         K = K - 1
       END DO

       NBIN = NB1 + 1
       BINCEN(NBIN) = 0.0D0

       DO J = 1, NB2
         NBIN = NBIN + 1
         BINCEN(NBIN) = J*WIDTH
       END DO

       DO J = 1, NBIN
         BINLIM(J) = BINCEN(J) - HWIDTH
       END DO

       BINLIM(NBIN+1) = BINCEN(NBIN) + HWIDTH

       IF (D2.GT.BINLIM(NBIN+1)) THEN
         NBIN = NBIN + 1
         BINCEN(NBIN) = 5.0D-1*(BINLIM(NBIN) + D2)
         BINLIM(NBIN+1) = D2
       END IF

       RETURN
       END

       SUBROUTINE MAKE_LFBINS(BRIGHT,FAINT,WIDTH,BINCEN,BINLIM,NBIN)

*  1998-May-24
*  This subroutine sets up the magnitude bins for the luminosity functions. 
*  The strategy used gives bin limits which are integer multiples of the
*  chosen bin width, except for the brightest bin at the RGB tip, which 
*  has the magnitude of the brightest isochrone point as its upper bin limit.

       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION BINCEN(*),BINLIM(*)

       HWIDTH = 5.0D-1*WIDTH

       IF (FAINT.GT.0.0D0) THEN
         NB1 = INT((FAINT - HWIDTH)/WIDTH)
       ELSE
         NB1 = INT((FAINT + HWIDTH)/WIDTH)
       END IF

       NB1 = MAX(1,NB1)

       BINLIM(1) = NB1*WIDTH + HWIDTH

       I = 1
       FBl = BINLIM(1)
       NBIN = INT((FBL - BRIGHT)/WIDTH) + 1

       DO I = 2, NBIN
         IM1 = I - 1
         BINLIM(I) = FBL - (IM1)*WIDTH
         BINCEN(IM1) = BINLIM(I) + HWIDTH
       END DO

       BINCEN(NBIN) = 5.0D-1*(BINLIM(NBIN) + BRIGHT)
       BINLIM(NBIN+1) = BRIGHT

       RETURN
       END

       SUBROUTINE MAKE_CFBINS(DCF,BLUE,RED,WIDTH,NBIN,BINCEN,BINLIM)

*  1998-May-25
*  This subroutine sets up the bins for the color functions and the color 
*  function ratios. When DCF is .TRUE., the bin limits are set as integer 
*  multiples of the bin width, running from the turnoff point (BLUE) up the 
*  RGB. The last bin at the tip of the RGB has RED as its outer limit.
*  When DCF is .FALSE. three bins are defined across the SGB, as described 
*  in Bergbusch & VandenBerg (1997 AJ,   ,  ) for the construction of bin 
*  ratios.

       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       LOGICAL DCF
       DIMENSION BINCEN(*),BINLIM(*)
       PARAMETER (OFFSET = 5.0D-2)

       HWIDTH = 5.0D-1*WIDTH

       IF (DCF) THEN

         IF (BLUE.LT.0.0D0) THEN
           BINLIM(1) = INT(BLUE/WIDTH)*WIDTH
         ELSE
           BINLIM(1) = (INT(BLUE/WIDTH) + 1)*WIDTH
         END IF

         NBIN = INT((RED - BINLIM(1))/WIDTH)

         DO I = 1, NBIN
           IP1 = I + 1
           BINLIM(IP1) = BINLIM(1) + I*WIDTH
         END DO

         IF (BINLIM(IP1).LT.RED) THEN
           NBIN = IP1
           binlim(nbin+1) = red
         END IF

         if (blue.lt.binlim(1)) then
           nbin = nbin + 1
           do i = nbin+1, 2, -1
             binlim(i) = binlim(i-1)
           end do
           binlim(1) = blue
         end if

       ELSE
         DBMV = RED - BLUE - OFFSET
         WIDTH = DBMV/3.0D0
         BINLIM(1) = BLUE + OFFSET
         BINLIM(2) = BINLIM(1) + WIDTH - 3.0d-1*offset
         BINLIM(3) = BINLIM(1) + 2*WIDTH - 3.0d-1*offset
         BINLIM(4) = RED
       END IF

       DO I = 1, NBIN
         BINCEN(I) = 5.0D-1*(BINLIM(I) + BINLIM(I+1))
       END DO

       RETURN
       END

       REAL FUNCTION DELTAMAG*8(X1,X2,Y1,Y2)

*      1998-MAY-21
*      The color index is the X-coordinate and the magnitude is the 
*      Y-coordinate. The color scale is expanded by a factor of 4 in order
*      to make the structure of the subgiant region easier to see.

       REAL*8 X1,X2,Y1,Y2

       DELTAMAG = SQRT((Y2-Y1)**2 + 1.6D1*(X2-X1)**2)

       RETURN
       END
       
