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