c------------------------------------------------------------------------c program isolf c Aaron Dotter, Steve Bjork, & Brian Chaboyer | updated 6/24/08 c c Contact aaron.dotter@dartmouth.edu with questions or comments c c Program reads in tabulated isochrone files with any number of c c ages and writes out Luminosity Functions based on user-supplied c c input. c c To compile: g77 -o isolf isolf.f c c The program must be run in the directory where the isochrones c c are located. The program assumes either: c c - a power-law IMF where dN/dM ~ M^x c c - a log-normal IMF where dN/dM ~ exp(-log(M/Mc)^2/2sigma^2) c c c c Output file includes a header similar to that of the input c c isochrone, followed by sequential lines of LF data. Including: c c 1. bin number c c 2. absolute magnitude of the bin in the chosen filter c c 3. log_10 cumulative number of stars, total number is nstars c c (default is nstars=10^5, to change, edit nstars below) c c 4. log_10 number of stars in each bin 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 (if 'IArgC' c c and 'getarg' are supported as in g77) 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 ./isolf fehm10afep0.cmd test.lf 8 4 0.1 1 -2.35 c c will generate a new file called 'test.lf' using empirical colors c c at [Fe/H]=-1, [alpha/Fe]=0, LF in F606W, bin size 0.1 mag, and c c power-law Salpeter IMF (x=-2.35). Note that the lognormal IMF c c requires one additional parameter on the command line (sigma). c c c c Note: in interactive mode it is best if both input and output c c files are in the same directory as the isolf executable. c c------------------------------------------------------------------------c implicit none logical lquiet integer m0,m1,m2,mdl,num,jjj,i,j,k,jj,nn,n0,j0,n1,mgm,nages,max, * iage,npts,ier,kc,junk,nstars,imf,iclr,ncol parameter(max=300,mgm=600,nstars=100000,ncol=10) integer cols(ncol), corr double precision dnt(mgm),dntl,ns,nsl,msave(max),tsave(max),logg, * tau,teff,te(max),logl,mag(max),smt,sml(max),dmag,mb0,sigma, * mb,mtot,dm,dn(mgm,3),mgrid(mgm),mx,mxmag,mlast,mnext,slast, * snext,blue,v(77),v1,m,lognormal,tiny,mk1,mk2 character infile*50,outfile*50,hdr1*60,hdr2*60,hdr3*60, * sep*60,nagestr*16,arg*10,phot(ncol)*12,filter*7 data lquiet/.false./ data tiny/1.0d-12/ data cols/4,4,8,8,8,5,7,7,4,77/ data phot/'BVRI ','uvby ','UBVRIJHKs ', * 'HST/ACS-WFC ','HST/WFPC2 ','SDSS ugriz ', * 'ACS GGCT-Emp','ACS GGCT-Syn','Spitzer IRAC', * 'HST WFC3'/ lquiet=.false. c run in quiet mode - remove comments below if(IArgC().gt.0) then lquiet=.true. call getarg(1,infile) call getarg(2,outfile) call getarg(3,arg) read(arg,'(i2)') iclr call getarg(4,arg) read(arg,'(i2)') m0 call getarg(5,arg) read(arg,'(f10.5)') dmag call getarg(6,arg) read(arg,'(i1)') imf call getarg(7,arg) read(arg,'(f10.5)') mx if(imf.eq.2) call getarg(8,arg) read(arg,'(f10.5)') sigma goto 1 endif c run interactively - skip this section if run from the command line c print header info and get user-supplied input !Input file print *, '-------------------------------------------' print *, ' DSEP Luminosity Function Program' print *, '-------------------------------------------' print *, " Enter isochrone filename:" read (*,'(a)') infile !Output file print *, " Enter output filename:" read (*,'(a)') outfile !Choice of photometric system and filter print *, '-------------------------------------------' print *, ' Choose photometric system: (1-8)' do i=1,ncol write(*,'(a3,i2,a2,a12)') ' [',i,'] ',phot(i) enddo print *, ' Enter choice of phot. system:' read *, iclr !fudge factor for filter names and numbers for HST/WFC3 if(iclr.eq.10) then corr = 5 else corr = 0 endif print *, ' Choose filter to build LF from:' do i=1,cols(iclr)/2 write(*,'(2(8x,i2,4x,a7))') i+corr, filter(iclr,i), * i+cols(iclr)/2+corr, filter(iclr,i+cols(iclr)/2) enddo if(mod(cols(iclr),2).eq.1)then ! if cols(iclr) is odd write(*,'(8x,i2,4x,a7)') cols(iclr)+5, filter(iclr,cols(iclr)) endif !write(*,'(8(4x,i2,4x))') (i,i=1,cols(iclr)), !write(*,'(8(1x,a7,1x))') (filter(iclr,i),i=1,cols(iclr)) print *, " Enter choice of filter:" read *, m0 !Bin size print *, '-------------------------------------------' print *, " Enter width of LF bins (in mag):" read *, dmag !IMF print *, '-------------------------------------------' print *, ' Enter choice of IMF:' print *, ' 1. Power-law 2. Log-normal' read *, imf if(imf.eq.1) then !Power law IMF version: print *, '-------------------------------------------' print *, " IMF power law: dN/dM ~ M^x" print *, " Enter value of x:" read *, mx else if(imf.eq.2) then !Log-normal IMF version print *, '-------------------------------------------' print *, " Log-normal IMF" print *, " Enter value of Mc and Sigma:" read*, mx, sigma endif 1 continue !redundant... but necessary to catch corr in quiet mode !fudge factor for filter names and numbers for HST/WFC3 if(iclr.eq.10) then corr = 5 else corr = 0 endif !correct m0, see above m0 = m0 - corr m1 = 0 m2 = 0 c set up color based on choice of filter if(iclr.eq.7.or.iclr.eq.8)then c ACS GGC includes filters in 3 different photometric systems if(m0.eq.1) then m1=1 m2=2 elseif(m0.eq.2.or.m0.eq.3)then m1=2 m2=3 elseif(m0.eq.4.or.m0.eq.5)then m1=4 m2=5 elseif(m0.eq.6.or.m0.eq.7)then m1=6 m2=7 endif else if(iclr.eq.10)then if(m0.eq.39)then !m0=F814W m1=25 !F606W m2=39 else if(m0.eq.75)then !m0=F160W m1=65 !F110W m2=75 else if(m0.lt.63)then !m0 in UVIS, use F814W m1=m0 m2=39 else !m0 in IR, use F160W m1=m0 m2=75 endif else if(m0.lt.cols(iclr)) then m1=m0 m2=m0+1 else m1=m0-1 m2=m0 endif endif open(unit=1,file=infile,status='old',err=900) open(unit=2,file=outfile,status='unknown',err=901) goto 2 900 if(.not.lquiet) print *,'Input file ',infile,' does not exist.' 901 if(.not.lquiet) print *,'Output file ',outfile,' does not exist.' stop c read composition, etc from input file 2 read(1,'(a16,i2)') nagestr, nages read(1,'(a)') sep read(1,'(a)') hdr1 read(1,'(a)') hdr2 read(1,*) read(1,'(a)') hdr3 read(1,*) c write header to output file write(2,'(a16,i2)') nagestr, nages write(2,'(a60)') sep write(2,'(a60)') hdr1 write(2,'(a60)') hdr2 write(2,'(a60)') sep write(2,'(a60)') hdr3 write(2,'(a60)') sep c loop back here to read in each new age; c calculate smt from nstar, for selected s, to get greater precision do jjj=1,nages read(1,'(5x,f6.3,6x,i3)') tau, npts read(1,*) iage= nint(tau*1.0d3) do i=1,npts if (i.gt.max) stop 'i exceeded dimensions' read(1,'(1x,i3,f10.6,99f8.4)')mdl,smt,logg,teff,logl, * (v(k),k=1,cols(iclr)) tsave(i)=v(m1)-v(m2) msave(i)=v(m0) if(i.eq.1) then v1=v(1) mb0=msave(1) endif sml(i)=log10(smt) enddo if(jjj.lt.nages) read(1,*) if(jjj.lt.nages) read(1,*) c store mags and colors in arrays to be used in calculating the LF do j=1,npts mag(j)=msave(j) te(j)=tsave(j) enddo c locate bluest point blue=te(1) do i=2,npts if(te(i).lt.blue) then blue=te(i) mb=mag(i) endif enddo m=(mb0-mb)/dmag mb=mb+m*dmag c find mag(j) next fainter than mb and set mag(j0)=mb do j=2,npts if(mag(j).lt.mb) go to 150 enddo if(.not.lquiet) write(*,*), 'all isochrone pts for',tau, * ' are fainter than lf ', 'lower bound' go to 100 150 j0=j-1 n0=j-2 n1=3 nn=4 if(n0.eq.0.or.(mag(j)-mag(j0))*(mag(j0)-mag(n0)).le.0.0d0)then n0=j0 n1=2 nn=3 endif if(j+1.gt.npts.or.(mag(j+1)-mag(j))*(mag(j)-mag(n0)).le.0.0d0) * nn=nn-1 call parrot (mag(n0),sml(n0),n1,nn,mb,sml(j0),ier) mag(j0)=mb c define mgrid to be faint boundary of mag intervals mxmag = 1.0d2 do i=j0,npts if(mag(i).lt.mxmag) mxmag=mag(i) enddo mgrid(1)=mb do j=2,mgm mgrid(j)=mgrid(1)-(j-1)*dmag if(mgrid(j).lt.mxmag) go to 170 enddo stop 'exceeded mgrid dimension' 170 jj=j if(iage.eq.0) jj=jj-1 c zero out arrays which will contain sums of stars, etc in each interval c and column do j=1,jj+1 do k=1,3 dn(j,k)=0.0d0 enddo enddo c set up parameters for initial point; c kc is number of output columns in table; normally 1, but 2 or 3 c if mags decrease between turnoff and base of rgb i=j0 j=1 kc=1 mlast=mag(i) slast=sml(i) c follow isochrone, dividing nstar into mag strips, and columns (kc) c check to see that points and grid intervals don't get out of phase 200 if(mlast.lt.mgrid(j+1).or.mlast.gt.mgrid(j)) 1 stop 'error - mag grid out of phase' if(mag(i+1)-mlast) 210,220,230 c next mag is brighter than mlast; check kc and c move up one box in grid if necessary 210 if(kc.eq.2.and.mag(i-1).lt.mag(i)) kc=3 if(mlast.eq.mgrid(j+1)) j=j+1 220 mnext=mag(i+1) snext=sml(i+1) if(mag(i+1).ge.mgrid(j+1)) go to 250 c next mag is not in same mgrid interval; interpolate on grid boundary mnext=mgrid(j+1) n0=i-1 n1=3 nn=4 if(n0.eq.0.or.mag(i-1).le.mag(i)) then n0=i n1=2 nn=3 endif if(i+2.gt.npts.or.mag(i+1).le.mag(i+2)) nn=nn-1 call parrot (mag(n0),sml(n0),n1,nn,mnext,snext,ier) if(ier.ne.0) then if(.not.lquiet) print *, '2. parrot error ', ier stop endif go to 250 c next mag is fainter than mlast; check kc and c move down one box in grid if necessary 230 if(kc.eq.1.and.mag(i).lt.mag(i-1)) kc=2 if(mlast.eq.mgrid(j)) j=j-1 mnext=mag(i+1) snext=sml(i+1) if(mag(i+1).le.mgrid(j)) go to 250 c next mag is not in same mgrid interval mnext=mgrid(j) n0=i-1 n1=3 nn=4 if(n0.eq.0.or.mag(i).le.mag(i-1)) then n0=i n1=2 nn=3 endif if(i+2.gt.npts.or.mag(i+2).le.mag(i+1)) nn=nn-1 call parrot (mag(n0),sml(n0),n1,nn,mnext,snext,ier) if(ier.ne.0) then if(.not.lquiet) print *, '3. parrot error ', ier stop endif c bc allow linear extrapolation c sum fraction of dn in mag strip 250 if(snext.eq.slast) go to 280 if(imf.eq.1.and.mx.eq.-1.0d0) then dm=snext-slast else if(imf.eq.1.and.mx.ne.-1.0d0)then dm=10.d0**((mx+1.0d0)*snext)-10.0d0**((mx+1.0d0)*slast) else if(imf.eq.2) then dm=lognormal(snext,mx,sigma)-lognormal(slast,mx,sigma) else if(imf.eq.3) then !Kroupa,Tout,&Gilmore(1993)IMF if(slast.gt.1.0d0) mk1=10.0d0**(-1.7d0*slast) if(slast.le.1.0d0) mk1=10.0d0**(-1.2d0*slast) if(slast.le.0.5d0) mk1=10.0d0**(-0.3d0*slast) if(snext.gt.1.0d0) mk2=10.0d0**(-1.7d0*snext) if(snext.le.1.0d0) mk2=10.0d0**(-1.2d0*snext) if(snext.le.0.5d0) mk2=10.0d0**(-0.3d0*snext) dm=mk2-mk1 endif dn(j,kc)=dn(j,kc)+dm c move to next segment 280 mlast=mnext slast=snext if(mnext.eq.mag(i+1)) then i=i+1 if(i.eq.npts) go to 300 endif go to 200 c go back to one column if mag decrease spans less than one interval 300 if(kc.gt.1) then do j=2,jj if(dn(j,2).gt.0.0d0.and.dn(j-1,2).eq.0.0d0 * .and.dn(j+1,2).eq.0.0d0) kc=1 if(kc.eq.1) then dn(j,1)=dn(j,1)+dn(j,2)+dn(j,3) dn(j,2)=0.0d0 dn(j,3)=0.0d0 endif enddo endif c write individual LF header junk=jj-2 if(imf.eq.1) write(2,1010) tau, junk, mx if(imf.eq.2) write(2,1012) tau, junk, mx, sigma 1010 format('#AGE=',f6.3,' BINS=',i3,' Power-law IMF: x=',f6.3) 1012 format('#AGE=',f6.3,' BINS=',i3,' Log-normal IMF: Mc=',f5.3, * ' sigma=',f5.3) write(2,1011) filter(iclr,m0) 1011 format('#',1x,'N',2x,a7,6x,'Log10(N)',7x,'Log10(dN)') c find total dn ns = 0.0d0 do j=jj-2, 1, -1 dnt(j)=dn(j,1)+dn(j,2)+dn(j,3) c srb 7/03 calculate also the integrated number of stars n ns = ns + dnt(j) enddo c break into two loops to rescale the total number of stars to a constant mtot=dble(nstars)/ns ns=0.0d0 num=0 do j=jj-2,1,-1 num=num+1 dnt(j)=mtot*dnt(j) + tiny !so that the number is never zero ns = ns + dnt(j) !for the next couple of lines if(ns.ne.0.0d0) nsl = log10(ns) c srb 2/03 added output of log(dn), and the integrated number of stars n if(dnt(j).gt.0.0d0) dntl = log10(dnt(j)) c write table write(2,'(i3,f7.3,2f15.4)') num,mgrid(j),nsl,dntl enddo if(jjj.lt.nages) write(2,*) if(jjj.lt.nages) write(2,*) c read more ages 100 continue enddo c close files and quit close(1) close(2) if(.not.lquiet) then print *, '-------------------------------------------' print *, " LF generation completed." print *, " Output written to ", outfile print *, '-------------------------------------------' endif c All done! end c----------------------------------------------------------------------- subroutine parrot (x,y,j,n,x0,y0,ierr) c parrot interpolates between a set of points in arrays x and y. (the x c array must be monotonic, increasing or decreasing.) the curve produced c by parrot follows the data very closely: polynomial "wiggles" between c points are almost entirely eliminated. the first and second derivatives c are everywhere continuous. c parrot finds one rotated parabola passing through the points j-2, j-1, c and j, with vertex exactly at j-1, and a second parabola passing through c j-1, j, and j+1, with vertex at j, then determines the y values c corresponding to x0 for each curve. the returned value y0 is the c weighted sum of the y values on the two parabolas, with the weights c determined by the distance of the interpolated point from points j and c j-1, respectively. (in the first and last intervals, the result is a c weighted average of a parabola and a straight line.) c c input parameters (unchanged by parrot) : c x(i),y(i) - arrays of known points to be interpolated; x monotonic. c n - largest defined index of x(i) and y(i) c x0 - point for which interpolated value is to be returned c j - x(j) and x(j-1) must bracket x0 c c output parameters : c y0 - returned value corresponding to x0 c ierr - error flag = 0 if no errors; = 1 if x0 is outside the c interval between x(1) and x(n); = 2 if x(i) are not c monotonic; = 3 if x0 is outside the interval between c x(j) and x(j-1). if ierr = 1, y0 is a linear c extrapolation. y0 is not set for ierr > 1. implicit none integer ierr, n, i, j, l, j0, j1, j3, iline double precision x(n),y(n),x0,y0,x1,y1,fy,fx,x3, y3,res(2),d1,d3 double precision r, d, ca, sa, rca, p, p2, q, q3, root,theta, cc double precision alpha, alt, phi, test, w, cw, sw, xp,yp,xp2,cw2 double precision a, b, c, aa, xa, dx, dy, wt1, wt2, xo, qr, save c check for error conditions ierr=0 c ierr=1 ? if((x(n)-x0)*(x0-x(1)).lt.0.0d0) go to 400 c ierr=2 ? do l=2,n-1 if((x(l-1)-x(l))*(x(l)-x(l+1)).le.0.0d0) go to 410 enddo c ierr=3 ? if((n-j)*(j-1).lt.0.) go to 420 if((x(j)-x0)*(x0-x(j-1)).lt.0.0d0) go to 420 do 300 i=1,2 j0=i+j-2 xo=x0-x(j0) c translate coordinates to vertex of each future parabola, with x1 c always chosen to be closest to x0 (+ linear extrapolation past ends) c (since j0 is either j or j-1, j1 is chosen to be the other) j1=j-i+1 x1=x(j1)-x(j0) y1=y(j1)-y(j0) if(j0.gt.1) go to 110 x3=x(1)-x(2) y3=y(1)-y(2) go to 125 110 if(j0.lt.n) go to 120 x3=x(n)-x(n-1) y3=y(n)-y(n-1) go to 125 120 j3=j+3*i-5 x3=x(j3)-x(j0) y3=y(j3)-y(j0) 125 res(i)=0.0d0 if(y1.eq.0.0d0.and.y3.eq.0.0d0) go to 195 c scale x coordinates to approximately same size as y coordinates fx=max(abs(x3),abs(x1)) fy=max(abs(y3),abs(y1)) if(fy.eq.0.0d0) fy=1.0d0 x1=x1/fx x3=x3/fx xo=xo/fx y1=y1/fy y3=y3/fy c determine angle of rotation w such that a parabola of the form y=a*x**2 c passes thru the three points (if theta is the angle from the origin c and (x1,y1) to the new (rotated) y axis; need to solve a cubic equation c in y = tan(theta): y**3 + p*y**2 + 2*r*ca*y - r*sa = 0, where c p = ca*(1-r*ca)/sa, ca = cos(alpha), sa = sin(alpha), and alpha is c the angle between the lines from the origin to (x1,y1) and (x3,y3). c if y=x-p/3, can get equation of form x**3 + a*x + b = 0. d1=x1*x1+y1*y1 d3=x3*x3+y3*y3 r=sqrt(d3/d1) d=sqrt(d1*d3) ca=(x1*x3+y1*y3)/d sa=(x1*y3-x3*y1)/d if(abs(sa).lt.0.001d0) go to 150 rca=r*ca p=ca*(1.0d0-rca)/sa p2=p*p q=2.0d0*rca-p2/3.0d0 r=2.0d0*(p2*p/27.0d0-p*rca/3.0d0)-r*sa q3=q*q*q/27.0d0 root=0.25d0*r*r+q3 if(root.gt.0.0d0) then ! 130,129,128 root=sqrt(root) theta=atan(qr(-0.5d0*r+root)+qr(-0.5d0*r-root)-p/3.0d0) go to 180 else if (root.eq.0.0d0) then c evaluate 2 solutions; find theta (0