c------------------------------------------------------------------------c program iso_interp_feh c Aaron Dotter c c Contact dotter@uvic.ca with questions or comments c c 2008 orig: includes full range of photometric systems c c 2008 July: bug fix in setting [a/Fe] for Y=0.33,0.4 (AS) c c 2008 Sept: minor fixes, set default interpolation to cubic (JR) c c 2008 Oct: bug fix in setting -0.5<[Fe/H]<0 with Y=0.33,0.4 (PG) c c 2008 Nov: perform interpolation on [Fe/H] instead of Z/X gives c c smooth results over the entire range of [Fe/H] (TP) c c c c Program interpolates amongst existing isochrones with any number c c of ages (<50) and writes the new isochrones in the same format. c c c c To compile: g77 -o iso_interp_feh iso_interp_feh.f c c c c The program must be able to access the isochrone grids via the c c following path: ./isochrones/[bvri/jc2mass/hst_acs/etc.] c c Interpolation is performed on [Fe/H] at constant Y and [a/Fe]. c c The order or degree of interpolation is set by the integer NISO. c c For example, NISO=2 gives linear interpolation while NISO=4 gives c c cubic interpolation. Default is cubic but reverts to linear near c c the edges of the [Fe/H] grid. c c c c Output from this program can be split into individual ages with c c the isolf_split program. c c c c This program may be run quietly from the command line c c (if 'IArgC' and 'getarg' are supported) by uncommenting the lines c c under 'run in quiet mode' below. To run from the command line c c simply supply the arguments in the same order as requested by c c the interactive mode, for example: c c c c ./iso_interp_feh 1 1 2 -1.3 test.iso c c c c will generate a new file called 'test.iso' using empirical BVRI c c colors with Y=0.25+1.5Z, [alpha/Fe]=0, [Fe/H]=-1.3. c c------------------------------------------------------------------------c implicit none integer i,j,k,l,max,isonum,nl,ifeh,ieep,neweep,fn,ncol integer niso,newimin,newimax,nz,iafe,nint,nafe,iclr,iy,nz1,nz2 parameter(niso=4,max=50,nl=350,nz1=8,nz2=6,isonum=1,nafe=6,ncol=9) integer neeps(niso,max),eeps(niso,max,nl),nages(niso),ioff(niso), * cols(ncol) character astr*5,estr*6,nagestr*19,hdr1*53,sep*60,hdr4*128 character outfile*25,phot(niso)*65,iso_file*45,arg*10 logical lquiet double precision age(niso,max),mass(niso,max,nl),teff(niso,max,nl) * ,gl(niso,max,nl),lum(niso,max,nl),color(niso,max,nl,ncol), * newfeh,qf(niso),ff(niso),newm,newt,newg,newl,newclr(ncol), * y(niso),z(niso),zeff(niso),feh(niso),afe(niso),newzeff, * newz,newy,cln,newfeh0,feh0(nz1),mixl,afe0(nafe),newxeff data sep/'#----------------------------------------------------'/ data feh0/-2.5d0,-2.0d0,-1.5d0,-1.0d0,-0.5d0,0.0d0,0.3d0,0.5d0/ data afe0/-0.2d0,0.0d0,0.2d0,0.4d0,0.6d0,0.8d0/ data cols/4,6,8,8,8,5,7,7,4/ data lquiet/.false./ C run in quiet mode if(IArgC().gt.0) then lquiet=.true. call getarg(1,arg) read(arg,'(i1)') iclr call getarg(2,arg) read(arg,'(i1)') iy if(iy.eq.1) then nz=nz1 else nz=nz2 endif call getarg(3,arg) read(arg,'(i1)') iafe call getarg(4,arg) read(arg,'(f10.5)') newfeh call getarg(5,outfile) goto 1 endif C run interactively C obtain user input: print *, '-----------------------------------------------------' print *, ' Isochrone [Fe/H] Interpolation Program' print *, '-----------------------------------------------------' print *, ' Enter choice of color-Teff transformation: [1-9]' print *, ' [1] BVRI :Vandenberg & Clem 2003' print *, ' [2] uvby :Clem et al 2004' print *, ' [3] UBVRIJHKs :Synthetic + Bessell(1990)/2MASS' print *, ' [4] HST/ACS-WFC :Synthetic + Sirianni et al 2005' print *, ' [5] HST/WFPC2 :Synthetic + Sirianni et al 2005' print *, ' [6] SDSS ugriz :Synthetic + SDSS' print *, ' [7] ACS GGCT-Emp:Vandenberg & Clem + Sirianni' print *, ' [8] ACS GGCT-Syn:Synthetic + Bessell/Sirianni' print *, ' [9] Spitzer IRAC:Synthetic' read *, iclr if(iclr.lt.1.or.iclr.gt.ncol) stop 'INVALID CHOICE OF COLOR' print *, '-----------------------------------------------' print *, ' Enter choice of Y: [1-3]' print *, ' Y= [1] 0.245 + 1.5 Z' print *, ' [2] 0.33 ([Fe/H]<=0 only)' print *, ' [3] 0.40 ([Fe/H]<=0 only)' read *, iy if(iy.lt.1.or.iy.gt.3) stop 'INVALID CHOICE OF Y' if(iy.eq.1) then print *, '-----------------------------------------------' print *, ' Enter choice of [a/Fe]: [1-6]' print *, ' [a/Fe]= [1] -0.2' print *, ' [2] 0.0' print *, ' [3] +0.2' print *, ' [4] +0.4 ([Fe/H]<=0 only)' print *, ' [5] +0.6 ([Fe/H]<=0 only)' print *, ' [6] +0.8 ([Fe/H]<=0 only)' read *, iafe if(iafe.lt.1.or.iafe.gt.6) stop 'INVALID CHOICE OF [a/Fe]' else print *, '-----------------------------------------------' print *, ' Enter choice of [a/Fe]: [1-2]' print *, ' [a/Fe]= [1] 0.0' print *, ' [2] +0.4' read *, iafe if(iafe.eq.2) then iafe=4 else if(iafe.eq.1) then iafe=2 else stop 'INVALID CHOICE OF [a/Fe]' endif endif if(iy.eq.1) then nz=nz1 else nz=nz2 endif print *, '-----------------------------------------------' print *, ' Available [Fe/H] range is ',feh0(1),'->',feh0(nz) print *, ' Enter new [Fe/H] to interpolate:' read *, newfeh if(newfeh.lt.feh0(1).or.newfeh.gt.feh0(nz)) * print *, 'EXTRAPOLATION WARNING - [Fe/H] out of bounds' if(newfeh.gt.0.0d0.and.iafe.gt.3) * stop 'INVALID COMBINATION OF [Fe/H] and [alpha/Fe]' print *, '-----------------------------------------------' print *, ' Enter new filename:' read *, outfile print *, '-----------------------------------------------' 1 continue C locate new [Fe/H] in grid call locate(feh0,nz,newfeh,ifeh) C nint is the number of points to be used for interpolation nint=niso C revert to linear interpolation if near the edge of the grid if(ifeh.lt.nint/2.or.ifeh.gt.nz-(nint/2)) nint=2 C this section of code reads in the isochrone files do j=1,nint fn=j+ifeh-nint/2 C read header open(isonum,file=iso_file(fn,iy,iafe,iclr),status='old') read(isonum,'(a16,i2)') nagestr,nages(j) read(isonum,*) read(isonum,'(a)') hdr1 read(isonum,'(1x,f7.4,f8.4,1p2e11.4,0p2f7.2)') * mixl,y(j),z(j),zeff(j),feh(j),afe(j) read(isonum,*) read(isonum,'(25x,a65)') phot(j) read(isonum,*) if(.not.lquiet) print *, "READ IN [Fe/H]=",feh(j) C read lines of data do i=1,nages(j) read(isonum,'(a5,f6.3,a6,i3)') astr,age(j,i),estr,neeps(j,i) read(isonum,'(a)') hdr4 do k=1,neeps(j,i) !Mass,Teff,log g,Lum,colors of each EEP read(isonum,*) eeps(j,i,k),mass(j,i,k),teff(j,i,k), * gl(j,i,k),lum(j,i,k),(color(j,i,k,l),l=1,cols(iclr)) enddo if(i.lt.nages(j)) read(isonum,*) if(i.lt.nages(j)) read(isonum,*) enddo close(isonum) enddo C finished reading in isochrones C make sure all files contain the same number and range ages: do j=2,nint if(nages(j).ne.nages(j-1)) stop 'INCORRECT NUMBER OF AGES!' if(phot(j).ne.phot(j-1))stop 'COLOR-TEFF NOT SAME IN ALL FILES' do i=1,nages(j) if(age(j,i).ne.age(j-1,i)) stop 'AGES DO NOT MATCH!' enddo enddo C interpolation in [Fe/H] gives smoother results do j=1,nint qf(j)=feh(j) enddo call interp(qf,ff,newfeh,nint) C set new X, Y, Z, and calculate [Fe/H] from them C if result differs from desired result by >5% print warning message newy=0.0d0 newz=0.0d0 do j=1,nint newy=newy+ff(j)*y(j) newz=newz+ff(j)*z(j) enddo c adjust X and Z if [a/Fe] non-zero by introducing Xeff and Zeff c calculate Zeff using relation of Salaris, Chieffi & Straniero C (1993, ApJ, 414, 580) cln=log(1.0d1) newzeff=newz/(6.38d-1*exp(cln*afe0(iafe)) + 3.62d-1) newxeff=1.0d0-newy-newzeff newfeh0=log10(newzeff/(newxeff*2.29d-2)) if(.not.lquiet .and. * abs(newfeh-newfeh0)/abs(newfeh).gt.5.0d-2) then print *, '***WARNING - POOR X,Z INTERPOLATION***' print *, '[Fe/H](Input) - [Fe/H](Interp)=',newfeh-newfeh0 endif C open new file and write header information open(isonum,file=outfile) write(isonum,'(a16,i2)') nagestr,nages(1) write(isonum,'(a)') sep write(isonum,'(a)') hdr1 write(isonum,'(a1,f7.4,f8.4,1p2e11.4,0p2f7.2)') * '#',mixl,newy,newz,newzeff,newfeh,afe0(iafe) write(isonum,'(a)') sep write(isonum,'(a25,a65)') "#**PHOTOMETRIC SYSTEM**: ",phot(1) write(isonum,'(a)') sep C loop through ages do j=1,nages(1) newimax=eeps(1,j,neeps(1,j)) newimin=eeps(1,j,1) do k=2,nint if(newimax.gt.eeps(k,j,neeps(k,j))) * newimax=eeps(k,j,neeps(k,j)) if(newimin.lt.eeps(k,j,1)) newimin=eeps(k,j,1) enddo neweep=newimax-newimin+1 do k=1,nint ioff(k)=newimin-eeps(k,j,1) enddo write(isonum,'(a5,f6.3,a6,i3)') astr,age(1,j),estr,neweep write(isonum,'(a)') hdr4 C interpolate and write results for each EEP do i=1,neweep do k=1,nint do l=k,nint C make sure EEPs are constant across different files if(.not.lquiet.and. * eeps(k,j,i+ioff(k)).ne.eeps(l,j,i+ioff(l))) * print *, "EEPS NOT EQUAL" enddo enddo C zero new quantities newm=0.0d0 newt=0.0d0 newg=0.0d0 newl=0.0d0 do k=1,cols(iclr) newclr(k)=0.0d0 enddo C perform interpolation do k=1,nint newm=newm+ff(k)*mass(k,j,i+ioff(k)) newt=newt+ff(k)*teff(k,j,i+ioff(k)) newg=newg+ff(k)* gl(k,j,i+ioff(k)) newl=newl+ff(k)* lum(k,j,i+ioff(k)) do l=1,cols(iclr) newclr(l)=newclr(l)+ff(k)*color(k,j,i+ioff(k),l) enddo enddo C write new isochrone(s) to file ieep=i-1+newimin write(isonum,'(i4,f10.6,11f8.4)')ieep,newm,newt,newg,newl, * (newclr(k),k=1,cols(iclr)) enddo if(j.lt.nages(1)) write(isonum,*) if(j.lt.nages(1)) write(isonum,*) enddo close(isonum) C all done! end C------------------------------------------------------------------------ C finds the index, j, of x in array xx(n). assuming xx is sorted, x lies C between xx(j) and xx(j+1) -- from NUMERICAL RECIPES IN FORTRAN 77 subroutine locate(xx,n,x,j) implicit none integer n, j, jl, ju, jm double precision xx(n), x jl=0 ju=n+1 10 if(ju-jl.gt.1)then jm=(ju+jl)/2 if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then jl=jm else ju=jm endif goto 10 endif j=jl return end c------------------------------------------------------------------------ subroutine interp(a,b,x,n) c {b} are coefficients of the interpolating polynomial c {a} are the tabulated values for use in interpolation c x is the abscissa to be interpolated c n is the number of points to be used, interpolating polynomial c has order n-1 implicit none integer i,j,n double precision a(n),b(n),x do i=1,n b(i)=1.0d0 do j=1,n if(j.ne.i) b(i)=b(i)*(x-a(j))/(a(i)-a(j)) enddo enddo return end c------------------------------------------------------------------------ function iso_file(fn,iy,iafe,iclr) c there are currently 30 isochrone files, this function returns c the name of the proper file based on [Fe/H] and [a/Fe] implicit none integer fn, iafe, nz, iclr, i1, i2, dirl(9), sufl(9), iy parameter(nz=8) character iso_file*45, feh(nz)*6,afe(6)*5,suf(9)*9,dir(9)*9 data feh/'fehm25','fehm20','fehm15','fehm10','fehm05','fehp00', * 'fehp03','fehp05'/ data afe/'afem2','afep0','afep2','afep4','afep6','afep8'/ data dir/'bvri','uvby','jc2mass','hst_acs','hst_wfpc2', * 'sdss','gct_emp','gct_syn','spitzer'/ data suf/'bvri','uvby','jc2mass','hst_acs','hst_wfpc2', * 'ugriz','cmd','phx','irac'/ data dirl/4,4,7,7,9,4,7,7,7/ data sufl/4,4,7,7,9,5,3,3,4/ do i1=1,len(iso_file) iso_file(i1:i1)=' ' enddo i1=1 i2=11 iso_file(i1:i2)='isochrones/' i1=i2+1 i2=i2+dirl(iclr) iso_file(i1:i2)=dir(iclr) i1=i2+1 i2=i2+1 iso_file(i1:i2)='/' i1=i2+1 i2=i2+6 iso_file(i1:i2)=feh(fn) i1=i2+1 i2=i2+5 iso_file(i1:i2)=afe(iafe) if(iy.eq.2) then i1=i2+1 i2=i2+3 iso_file(i1:i2)='y33' else if(iy.eq.3) then i1=i2+1 i2=i2+3 iso_file(i1:i2)='y40' endif i1=i2+1 i2=i2+1 iso_file(i1:i2)='.' i1=i2+1 i2=i2+sufl(iclr) iso_file(i1:i2)=suf(iclr) return end c------------------------------------------------------------------------ c------------------------------------------------------------------------