PROGRAM VRISO * Last modified February 3, 2006 * New features include implementation of the subroutine ISOOPT that allows * the user to vary the interpolation relations continuously from purely * linear to purely Akima spline. It is possible to select a hardwired * set of isochrone ages, to read a set of ages from a file (e.g. VRISO.OPT), * or to enter them interactively. * This program is intended for use with the grids of evolutionary tracks * computed by Don VandenBerg at the University of Victoria. PARAMETER (NPM=3000, NTM = 50, NEM = 10) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION T(NPM,NTM),T0(NPM),A(NPM,NTM),A0(NPM),D(NPM,NTM), 1 DI(NPM),Z(NPM,NTM),DLDA(NPM,NTM),DTDA(NPM,NTM) REAL*8 L(NPM,NTM),L0(NPM),MASS(NTM),ICAGE(NTM),MI(NPM), 2 ICM,ICL,ICT,ICD,ICDL,ICDT,iam,ilm,itm INTEGER LEEP(NEM,NTM),NP(NTM) CHARACTER INFILE*60,OUTFIL*40,OLINE(NPM)*57,HLINE(7)*65,CLEEP*35 LOGICAL INTER,MORE,first common /opts/ iam,ilm,itm PARAMETER (CX = 2.302585092994045684D0) CALL IOFILE(8,'IN','FILE.EEP',.FALSE.,' ', 1 INFILE,.TRUE.,' EEP Tracks file:') CALL PARSE(INFILE,'NAME',OUTFIL,NCH) NCHOUT = NCH + 4 OUTFIL(1:NCHOUT) = OUTFIL(1:NCH)//'.iso' * Read header lines from the .EEP file DO I = 1, 7 READ(8,'(A)') HLINE(I) END DO * Extract the number of tracks from the first header line READ(HLINE(1),'(11x,I3)') NT *23456789012345678901234567890123456789012345678901234567890123456789012 WRITE(*,'(/T31,''ZAMS'',T45,'' Primary EEP Models'')') WRITE(*,'(T13,''Track'',T22,''Mass'',T31,''Ages'', 1 T37,''ZAMS'',T42,''MSTO'',T47,''BLHK'',T52,''HZGP'',T57,''BRGB'', 2 T62,''GBPS'',T67,''RGBT'',/)') * Read in the .EEP file track by track DO I = NT, 1, -1 * Read in the mass, number of points, and EEP locations for each track READ(8,'(//F6.3,I6,40X,5I3,2I4)') 1 MASS(I),NP(I),(LEEP(J,I),J=1,5),(LEEP(J,I),J=6,7) do j = 1, 35 cleep(j:j) = ' ' end do do j = 1, 7 if (leep(j,i).gt.0) then ku = 5*j kl = ku-4 write(cleep(kl:ku),'(i5)') leep(j,i) end if end do READ(8,'(4X,F10.6,F9.6,D11.4,f7.4)') L0(1),T0(1),AGE,Z(1,I) A0(1) = LOG10(AGE) D(1,I) = age DO N = 2, NP(I) READ(8,'(4X,F10.6,F9.6,D11.4,f7.4)') L0(N),T0(N),DA,Z(N,I) AGE = AGE + DA A0(N) = LOG10(AGE) D(N,I) = AGE END DO DO N = 1, NP(I) DLDA(N,I) = DIFRED(A0,L0,N,NP(I)) DTDA(N,I) = DIFRED(A0,T0,N,NP(I)) L(N,I) = L0(N) T(N,I) = T0(N) A(N,I) = A0(N) END DO AGE = 10**A(1,I) WRITE(*,'(T14,I3,T21,F5.3,T29,F6.3,T36,A35)') 1 NT-I+1,MASS(I),AGE,CLEEP END DO CLOSE(UNIT=8,STATUS='KEEP') CALL EEP_TRACKS(NT,NP,L,T,A,Z,D,DLDA,DTDA,LEEP,KNT) NPMAX = NP(1) DO I = 1, NT NPMAX = MAX(NPMAX,NP(I)) END DO WRITE(*,'(//)') CALL IOFILE(10,'OUT',OUTFIL(1:nchout),.TRUE., 1 'Default Isochrone file ', 1 OUTFIL,.FALSE.,' Isochrone file:') WRITE(*,'(//)') CALL ISOOPT(ICAGE,NIC) IF (NIC.LT.100) THEN WRITE(10,'('' ISOCHRONES '',I2)') NIC ELSE WRITE(*,'('' More than 99 isochrones! '')') STOP END IF DO I = 2, 7 J = 65 MORE = .TRUE. DO WHILE (MORE) IF (HLINE(I)(J:J).EQ.' ') THEN J = J - 1 ELSE WRITE(10,'(A)') HLINE(I)(1:J) MORE = .FALSE. END IF END DO END DO WRITE(*,'(/)') DO K = 1, NIC WRITE(*,'(T20,''Interpolating '',F5.2, 1 '' Gyr Isochrone.'')') ICAGE(K) ICPT = 0 AGE = LOG10(ICAGE(K)) PMASS = 0.0D0 first=.true. DO I = 1, NPMAX CALL ISO_INTERP (I,KNT,AGE,NP,NT,L,T,A,DLDA,DTDA, 1 MASS,ICM,ICL,ICT,ICDL,ICDT,DADM,DLDM,DTDM,INTER,first) IF (INTER.AND.ICM.GT.PMASS) THEN PMASS = ICM IF (ICPT.EQ.0) THEN PICL = ICL PICT = ICT ICD = 0.0D0 ELSE dICD = DELTAD(PICT,ICT,PICL,ICL) END IF if (icpt.eq.0) then ICPT = 1 di(1) = icd MI(1) = ICM DLDM = DLDM - ICDL*DADM DTDM = DTDM - ICDT*DADM * PARAMETER (XSCALE=1.0D2, YSCALE=1.25D1) DMDD = 1.0D0/SQRT(1.0D2*DTDM*DTDM + 1.5625D0*DLDM*DLDM) *23456789012345678901234567890123456789012345678901234567890123456789012 WRITE(OLINE(1), '(f10.6,f9.6,f13.10,f11.7,1PD14.7)') 1 ICL,ICT,ICM,ICD,DMDD else if (dicd.ge.2.0d-2) then ICPT = ICPT + 1 PICT = ICT PICL = ICL icd = icd + dicd DI(ICPT) = ICD MI(ICPT) = ICM DLDM = DLDM - ICDL*DADM DTDM = DTDM - ICDT*DADM * PARAMETER (XSCALE=1.0D2, YSCALE=1.25D1) DMDD = 1.0D0/SQRT(1.0D2*DTDM*DTDM + 1.5625D0*DLDM*DLDM) *23456789012345678901234567890123456789012345678901234567890123456789012 WRITE(OLINE(ICPT), '(f10.6,f9.6,f13.10,f11.7,1PD14.7)') 1 ICL,ICT,ICM,ICD,DMDD end if END IF END DO WRITE(10,'(/'' Age Npts'')') WRITE(10,'(F6.2,I6)') ICAGE(K),ICPT DO I1 = 1, ICPT Dmdd = DIFRED(di,mi,i1,icpt) WRITE(10,'(I4,A,1pd14.7)') I1,OLINE(I1),dmdd END DO END DO CLOSE(UNIT=10,STATUS='KEEP') CLOSE (UNIT=12,STATUS='KEEP') STOP END SUBROUTINE EEP_TRACKS(NTRK,NP,L,T,A,X,D,DLDA,DTDA,LEEP,K) PARAMETER (NPM=3000, NTM = 50, NEM = 10) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION T(NPM,*),A(NPM,*),D(NPM,*),X(NPM,*),DLDA(NPM,*), 1 DTDA(NPM,*),TEFF(NPM,NTM),AGE(NPM,NTM),XX(NPM,NTM),R(NEM), 2 DD(NPM,NTM),DLDA2(NPM,NTM),DTDA2(NPM,NTM),F(NEM),FF(NEM) REAL*8 L(NPM,*),LUM(NPM,NTM) INTEGER LEEP(NEM,*),NP(*),NP2(NTM),NEEP(NTM),N2EEP(NEM) LOGICAL BLHK, RGBT *23456789012345678901234567890123456789012345678901234567890123456789012 * NTRK input number of tracks in the grid * NP input/output array of number of points per track * L input/output array of log L * T input/output array of log Teff * A input/output array of log age * X input/output array of central hydrogen content * D input/output array of age * DLDA input/output array of d(log L)/d(log t) * DTDA input/output array of d(log Teff)/(d log t) * LEEP input/output array of track point indices that identify EEPs * K output index of lowest mass track that has the MSTO identified * An underlying assumption is that if an EEP doesn't exist on the track * of the most massive star, then it will not be found on the tracks of * lower mass stars. (For example, the Blue Hook only appears in the * tracks of more massive stars.) The following 13 lines of code search * for identified EEPs on the track for the most massive star in the * grid, that is the one with the index NTRK. If, for example, the BLHK * EEP was not identified, it will be eliminated from the interpolation * scheme. Similarly, any subsequent EEPs that were not identified will * be eliminated. * I=2 MSTO; I=3 BLHK; I=4 HZGP; I=5 BRGB; I=6 GBPS; I=7 RGBT MAXEEP = 7 * Track KBH is the highest mass track in which there is no hook KBH = NTRK IF (LEEP(3,NTRK).NE.0) THEN BLHK = .TRUE. DO WHILE (LEEP(3,KBH).GT.0) KBH = KBH - 1 END DO END IF * When there is a blue hook, the interpolation scheme becomes degenerate for * the grids which have the MSTO point. The BLHK EEP connects to the MSTO EEP * in the highest mass track without the hook identified. if (blhk) then do i = 1, kbh if (leep(2,i).gt.0) leep(3,i) = leep(2,i) end do end if * In grids extended to higher masses (say 2.0 Msol) the higher mass tracks * may terminate at or near the RGB pause. To ensure that interpolations are * performed to the RGB tip in lower mass tracks, it is necessary to identify * the first track containing the RGBT point. In the higher mass models, the * GBPS EEP connects to the first RGBT EEP identified. * do i = 1, ntrk * write(*,'(8i6)') i,(leep(j,i),j=1,7) * end do rgbt = leep(7,ntrk).eq.0 * write(*,'(i6)') leep(7,ntrk) krgbt = ntrk do while (rgbt) krgbt = krgbt - 1 rgbt = leep(7,krgbt).eq.0 * write(*,'(2i6,l3)') krgbt,leep(7,krgbt),rgbt end do if (krgbt.lt.ntrk) then krgbt = krgbt + 1 do i = krgbt, ntrk leep(7,i) = leep(6,i) * write(*,'(3i6)') i, leep(6,i),leep(7,i) end do end if * Now eliminate any primary EEPs that are not found in the highest mass track. I = 2 DO WHILE (I.LE.MAXEEP) IF (LEEP(I,NTRK).EQ.0) THEN MAXEEP = MAXEEP - 1 DO J = I, MAXEEP DO K = 1, NTRK LEEP(J,K) = LEEP(J+1,K) END DO END DO END IF I = I + 1 END DO * write(*,'(//t33,'' Primary EEP Table'',/)') * do i = 1, ntrk * write(*,'(t24,7i5)') (leep(j,i),j=1,maxeep) * end do * At this point, the first MAXEEP elements for model number NTRK in the * array LEEP are non-zero, and for the other models, only trailing * entries are set to zero. The next 7 lines of code count the number * of non-zero LEEP entries for each model star. The array NEEP stores * number of valid EEPs for each model. DO I = 1, NTRK J = maxeep DO WHILE (LEEP(J,I).EQ.0) J = J - 1 END DO NEEP(I) = J END DO * Identify the low mass tracks that have only the ZAMS EEP * K is the index of the lowest mass model that was evolved to the MSTO K = 1 DO WHILE (NEEP(K).EQ.1.AND.K.LT.NTRK) K = K + 1 END DO * Set the number of secondary EEP's on each phase interval. n2eep(1) = 800 * r(1) = -6.931471805599453d-1 r(1) = -2.310490601866351d-1 if (blhk) then n2eep(2) = 120 r(2) = -6.931471805599453d-1 n2eep(3) = 400 r(3) = -6.931471805599453d-1 n2eep(4) = 220 r(4) = -2.87682073d-1 n2eep(5) = 200 r(5) = -6.931471805599453d-1 n2eep(6) = 300 r(6) = -6.931471805599453d-1 else n2eep(2) = 400 r(2) = -6.931471805599453d-1 n2eep(3) = 220 r(3) = -2.87682073d-1 n2eep(4) = 200 r(4) = -6.931471805599453d-1 n2eep(5) = 300 r(5) = -6.931471805599453d-1 end if * Work out the scaling for the point distribution. As it is currently set up * the scaling is based on a single ratio (0.5) of the last to the first time * within a region bracketed by successive primary EEPs. do i = 1, neep(ntrk)-1 f(i) = exp(r(i)/(n2eep(i)-1)) ff(i) = 1.0d0 do kj = 2, n2eep(i) ff(i) = ff(i) + f(i)**(kj-1) end do end do * Interpolate track parameters to secondary EEP's on all tracks DO NT = K, NTRK NP2(NT) = 1 KEEP = LEEP(1,NT) DD(1,NT) = D(KEEP,NT) AGE(1,NT) = A(KEEP,NT) LUM(1,NT) = L(KEEP,NT) TEFF(1,NT) = T(KEEP,NT) DLDA2(1,NT) = DLDA(KEEP,NT) DTDA2(1,NT) = DTDA(KEEP,NT) XX(1,NT) = X(KEEP,NT) DO I = 1, NEEP(NT)-1 I1 = LEEP(I,NT) I2 = LEEP(I+1,NT) IF (I1.NE.I2) THEN dt0 = (D(I2,NT) - D(I1,NT))/ff(i) delta = 1.0d0 DO I3 = 1, N2EEP(I) NP2(NT) = NP2(NT) + 1 dist = d(i1,nt) + dt0*delta delta = f(i)**(i3) + delta I4 = KEEP DO WHILE (DIST.GT.D(I4,NT).AND.I4.LE.NP(NT)) I4 = I4 + 1 END DO KEEP = I4 - 1 FACTOR = (DIST - D(KEEP,NT))/(D(I4,NT) - D(KEEP,NT)) DD(NP2(NT),NT) = DIST AGE(NP2(NT),NT) = log10(dist) LUM(NP2(NT),NT) = L(KEEP,NT) + 1 FACTOR* (L(I4,NT) - L(KEEP,NT)) TEFF(NP2(NT),NT) = T(KEEP,NT) + 1 FACTOR* (T(I4,NT) - T(KEEP,NT)) DLDA2(NP2(NT),NT) = DLDA(KEEP,NT) + 1 FACTOR*(DLDA(I4,NT) -DLDA(KEEP,NT)) DTDA2(NP2(NT),NT) = DTDA(KEEP,NT) + 1 FACTOR*(DTDA(I4,NT) -DTDA(KEEP,NT)) XX(NP2(NT),NT) = X(KEEP,NT) + 1 FACTOR*(X(I4,NT) -X(KEEP,NT)) END DO ELSE DO I3 = 1, N2EEP(I) NP2(NT) = NP2(NT) + 1 DD(NP2(NT),NT) = D(I1,NT) AGE(NP2(NT),NT) = A(I1,NT) LUM(NP2(NT),NT) = L(I1,NT) TEFF(NP2(NT),NT) = T(I1,NT) DLDA2(NP2(NT),NT) = DLDA(I1,NT) DTDA2(NP2(NT),NT) = DTDA(I1,NT) XX(NP2(NT),NT) = X(I1,NT) END DO END IF END DO END DO * Extrapolate secondary EEP's to low mass tracks using Xcen as the criterion * for setting the 2nd primary EEP KM1 = K - 1 DO NT = KM1, 1, -1 * First find the interval on the track of the next highest mass that * brackets the Xcen of the last point of the current track. XCEN = X(NP(NT),NT) * XCEN = 0.9*X(NP(NT),NT) J = 1 NTP1 = NT + 1 DO WHILE (XCEN.LE.XX(J,NTP1).and.j.le.np2(ntp1)) J = J + 1 END DO * There are np2(nt) secondary EEPs on the (k-1)th track. The (j-1)th point in * the kth EEP track corresponds to the last useful point on the (k-1)th track, * so there will be j-1 2ndary EEPs on this track. NP2(NT) = J - 1 * Set up the scaling for the point distribution ff(1) = 1.0d0 do kj = 2, np2(nt)-1 ff(1) = ff(1) + f(1)**(kj-1) end do * Interpolate Age with respect to Xcen within the current track to establish * the time interval between the first (primary) EEP and the last secondary EEP. XCEN = XX(NP2(NT),NTP1) J = 1 DO WHILE (XCEN.LT.X(J,NT)) J = J + 1 END DO JM1 = J - 1 FACTOR = (XCEN - X(JM1,NT))/(X(J,NT) - X(JM1,NT)) * Calculate the difference in age between the last secondary EEP and the * ZAMS point on the current low mass track. dt0 = D(Jm1,NT) - d(1,nt) + factor*(d(j,nt) - d(jm1,nt)) * Now get the size of the first time step dt0 = dt0/ff(1) * The first point on the track is a primary EEP DD(1,NT) = D(1,NT) AGE(1,NT) = A(1,NT) LUM(1,NT) = L(1,NT) TEFF(1,NT) = T(1,NT) DLDA2(1,NT) = DLDA(1,NT) DTDA2(1,NT) = DTDA(1,NT) XX(1,NT) = X(1,NT) * Interpolate the interior secondary EEPs for this track with respect * to "distance" along the track KEEP = 1 DELTA = 1.0d0 DO I3 = 2, NP2(NT) DIST = d(1,nt) + dt0*delta delta = f(1)**(i3-1) + delta I4 = KEEP DO WHILE (DIST.GT.D(I4,NT).AND.I4.LE.NP(NT)) I4 = I4 + 1 END DO KEEP = I4 - 1 FACTOR = (DIST - D(KEEP,NT))/(D(I4,NT) - D(KEEP,NT)) DD(I3,NT) = DIST AGE(I3,NT) = log10(dist) LUM(I3,NT) = L(KEEP,NT) + 1 FACTOR*(L(I4,NT) - L(KEEP,NT)) TEFF(I3,NT) = T(KEEP,NT) + 1 FACTOR*(T(I4,NT) - T(KEEP,NT)) DLDA2(I3,NT) = DLDA(KEEP,NT) + 1 FACTOR*(DLDA(I4,NT) -DLDA(KEEP,NT)) DTDA2(I3,NT) = DTDA(KEEP,NT) + 1 FACTOR*(DTDA(I4,NT) -DTDA(KEEP,NT)) XX(I3,NT) = X(KEEP,NT) + 1 FACTOR*(X(I4,NT) -X(KEEP,NT)) END DO END DO * Overwrite the input arrays with the contents of the working arrays; that * is, the original tracks are replaced by the EEP tracks. DO NT = 1, NTRK NP(NT) = NP2(NT) DO I = 1, NP2(NT) D(I,NT) = DD(I,NT) A(I,NT) = AGE(I,NT) L(I,NT) = LUM(I,NT) T(I,NT) = TEFF(I,NT) DLDA(I,NT) = DLDA2(I,NT) DTDA(I,NT) = DTDA2(I,NT) X(I,NT) = XX(I,NT) END DO LEEP(1,NT) = 1 DO I = 2, NEEP(NT) LEEP(I,NT) = LEEP(I-1,NT) + N2EEP(I-1) END DO END DO * write(*,'(//t33,'' Secondary EEP Table'',/)') * do i = 1, ntrk * if (i.lt.k) leep(2,i)=np(i) * write(*,'(t24,7i5)') (leep(j,i),j=1,maxeep) * end do RETURN END SUBROUTINE ISO_INTERP(ID,KNT,AGE,NPT,NTRK,L,T,A,DLDA,DTDA, 1 MASS,ICM,ICL,ICT,ICDL,ICDT,DADM,DLDM,DTDM,INTER,first) PARAMETER (NTM = 50, NPM = 3000) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(NPM,NTM),T(NPM,NTM),DLDA(NPM,NTM),DTDA(NPM,NTM), 1 DL(NTM),DT(NTM),DA(NTM),DDL(NTM),DDT(NTM),DM(NTM), 2 CSCOEF(4,NTM) REAL*8 L(NPM,NTM),MASS(NTM),ICA,ICM,ICL,ICT,ICDL,ICDT REAL*8 IAM,ILM,ITM * REAL*8 W1(5,4),W2(6),C(4) INTEGER NPT(NTM) LOGICAL INTER, FIRST common /opts/ iam,ilm,itm NT = 1 DO WHILE (NPT(NT).LT.ID) NT = NT + 1 IF (NT.GT.NTRK) STOP 'NT > NTRK IN ISOCHRONE_INTERPOLATION' END DO KNT1 = NT NT = NTRK DO WHILE(NPT(NT).LT.ID) NT = NT - 1 IF (NT.LT.1) STOP 'NT < 1 IN ISOCHRONE_INTERPOLATION' END DO KNT2 = NT IF (KNT1.EQ.KNT2) THEN INTER = .FALSE. WRITE(6,'('' KNT1 = KNT2'')') RETURN ELSE IF (A(ID,KNT2).LE.AGE.AND.AGE.LE.A(ID,KNT1)) THEN INTER = .TRUE. NP = 0 DO I = KNT1, KNT2 NP = NP + 1 DM(NP) = MASS(I) DL(NP) = L(ID,I) DT(NP) = T(ID,I) DA(NP) = A(ID,I) DDL(np) = DLDA(ID,I) DDT(NP) = DTDA(ID,I) END DO * Interpolation to the isochrone age in the Age-Mass relation gives the * mass (icm) at this isochrone point. The slope of the relation at this * point allows us to evaluate the derivative d(log t)/dM (dadm) in the * direction of interpolation. CALL AKM_CSPL(NP,Da,Dm,CSCOEF) call akm_eval(np,da,cscoef,age,icm,dadm) dadm = 1.0d0/dadm i = 1 do while (da(i).gt.age) i = i + 1 end do im1 = i - 1 factor = (dm(i) - dm(im1))/(da(i) - da(im1)) icm = iam*icm + (1-iam)*(dm(im1) + factor*(age - da(im1))) dadm = iam*dadm + (1-iam)/factor * Now interpolate between the tracks to get the luminosity (icl), the * temperature (ict) along an evolutionary track of mass icm at the * isochrone point. CALL AKM_CSPL(NP,DM,DL,CSCOEF) CALL AKM_EVAL(NP,DM,CSCOEF,ICM,ICL,DLDM) factor = (dl(i) - dl(im1))/(dm(i) - dm(im1)) icl = ilm*icl + (1-ilm)*(dl(im1) + factor*(icm - dm(im1))) dldm = ilm*dldm + (1-ilm)*factor CALL AKM_CSPL(NP,DM,DT,CSCOEF) CALL AKM_EVAL(NP,DM,CSCOEF,ICM,ICT,DTDM) factor = (dt(i) - dt(im1))/(dm(i) - dm(im1)) ict = itm*ict + (1-itm)*(dt(im1) + factor*(icm - dm(im1))) dtdm = itm*dtdm + (1-itm)*factor CALL AKM_CSPL(NP,DM,DDL,CSCOEF) CALL AKM_EVAL(NP,DM,CSCOEF,ICM,ICDL,ICA) factor = (ddl(i) - ddl(im1))/(dm(i) - dm(im1)) icdl = ilm*icdl + (1-ilm)*(ddl(im1) + factor*(icm - dm(im1))) CALL AKM_CSPL(NP,DM,DDT,CSCOEF) CALL AKM_EVAL(NP,DM,CSCOEF,ICM,ICDT,ICA) factor = (ddt(i) - ddt(im1))/(dm(i) - dm(im1)) icdt = itm*icdt + (1-itm)*(ddt(im1) + factor*(icm - dm(im1))) RETURN ELSE INTER = .FALSE. RETURN END IF END REAL FUNCTION DELTAD*8(X1,X2,Y1,Y2) * 10 APR. 1987 REAL*8 X1,X2,Y1,Y2,XSCALE,YSCALE PARAMETER (XSCALE=1.0D1, YSCALE=1.25D0) DELTAD = SQRT((YSCALE*(Y2-Y1))**2 + (XSCALE*(X2-X1))**2) RETURN END SUBROUTINE ISOOPT(AGE,NA) IMPLICIT DOUBLE PRECISION (A-H, O-Z) REAL*8 IAM,ILM,ITM DIMENSION AT0(38), AGE(50), A(3,10) CHARACTER INFILE*40 LOGICAL FIRST, MORE COMMON /OPTS/ IAM,ILM,ITM DATA AT0/1.0D-2,2.0D-2,5.0D-2,1.0D-1,2.0D-1,5.0D-1,7.5D-1,1.0D0, 1 1.2D0,1.4D0,1.6D0,1.8D0,2.0D0,2.2D0,2.4D0,2.6D0,2.7D0, 2 2.8D0,2.9D0,3.0D0,3.1D0,3.2D0,3.3D0,3.4D0,3.5D0,3.6D0, 3 3.8D0,4.2D0,4.6D0,5.0D0,6.0D0,7.0D0,8.0D0,10.0D0,12.0D0, 4 14.0D0,16.0D0,18.0D0/ IERR = 0 CALL IOFILE(15,'IN','VRISO.OPT',.TRUE.,'Default Options File ', 1 INFILE,.FALSE.,' Options File: ') READ (15,'(/3F5.2)') IAM, ILM, ITM READ (15,'(I2)') ITEST NL = 0 FIRST = .TRUE. DO I = 1, 10 READ (15,'(F5.2,F6.2,F5.2)') (A(J,I),J=1,3) IF (A(1,I).EQ.0.AND.FIRST) THEN NL = I-1 FIRST = .FALSE. END IF END DO CLOSE (UNIT=15,STATUS='KEEP') write(*,'(//t25,''ITERPOLATION OPTIONS''/)') write(*,'(t15,''AM ='',f5.2, 1 '' Age-Mass interpolant; 1=spline, 0=linear''/ 2 t15,''LM ='',f5.2, 3 '' Luminosity-Mass interpolant; 1=spline, 0=linear''/ 4 t15,''TM ='',f5.2, 5 '' Temperature-Mass interpolant; 1=spline, 0=linear'')') 6 iam,ilm,itm IF (ITEST.EQ.1) THEN NA = 38 DO I = 1, 38 AGE(I) = AT0(I) END DO ELSE IF (ITEST.EQ.0) THEN NA = 0 DO I = 1, NL NR = NINT((A(2,I)-A(1,I))/A(3,I)) NA = NA + 1 IF (NA.EQ.1) THEN AGE(NA) = A(1,I) ELSE NA = NA - 1 END IF DO J = 1, NR NA = NA + 1 AGE(NA) = A(1,I) + J*A(3,I) END DO END DO ELSE WRITE(*,'(//)') NA = 0 MORE = .TRUE. DO WHILE (MORE) WRITE(*,'(T20,''Isochrone Age(Gyr) Type E to exit: '',$)') READ(*,*,IOSTAT=IERR) A0 IF (IERR.EQ.0.AND.A0.GT.0.0D0) THEN NA = NA + 1 AGE(NA) = A0 ELSE MORE = .FALSE. END IF END DO END IF RETURN END REAL FUNCTION DIFRED*8 (X,F,I,N) * Function subroutine DIFRED is a double precision subroutine that uses * evaluates the derivative of a function at the ith point using five point * differentiation IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),F(N),H(5),FX(5) IF (I.LT.3) THEN I1 = 1 ELSE IF (I.GT.(N-2)) THEN I1 = N - 4 ELSE I1 = I - 2 END IF DO K = 1,5 H(K) = X(I1) - X(I) FX(K) = F(I1) I1 = I1 + 1 END DO SHI2FI = 0.0D0 SHIFI = 0.0D0 SFI = 0.0D0 SHI = 0.0D0 SHI2 = 0.0D0 SHI3 = 0.0D0 SHI4 = 0.0D0 DO K = 1,5 SFI = SFI + FX(K) SHI = SHI + H(K) SHIFI = SHIFI + H(K)*FX(K) TEMP = H(K)*H(K) SHI2 = SHI2 + TEMP SHI2FI = SHI2FI + FX(K)*TEMP TEMP = TEMP*H(K) SHI3 = SHI3 + TEMP TEMP = TEMP*H(K) SHI4 = SHI4 + TEMP END DO TOP = SHI4*(5.0D0*SHIFI - SHI*SFI) + 1 SHI2*(SHI2FI*SHI + SHI3*SFI - SHI2*SHIFI) - 2 5.0D0*SHI3*SHI2FI BOT = SHI2*(5.0D0*SHI4 + 2.0D0*SHI3*SHI - SHI2*SHI2) - 1 SHI*SHI*SHI4 - 5.0D0*SHI3*SHI3 IF (BOT.LT.1.0D-38) BOT = 1.0D-38 DIFRED = TOP/BOT RETURN END