cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SHB by Aaron Dotter (dotter@uvic.ca) c c 1st version: Feb 2005 c c 2nd version: May 2006-includes interpolation between mass tracks c c 3rd version: Aug 2006-added [Fe/H] interpolation c c 4th version: Jun 2007-fixed bugs in stellar mass interpolation c c 5th version: Feb 2008-added more color options, minor improvements c c 6th version: Jul 2008-improved RNG seeding method c c 7th version: Jun 2009-initialize RNG with system_clock(iseed) c program SHB c *The program generates a synthetic horizontal branch (SHB) from a c c set of HB tracks and a user-specified mass distribution. c c *The program should be located in the top-level directory so that c c the path to the tracks is 'tracks/[color]/[composition]/hb* . c c *The program requires seeding of a random number generator. c c Currently the seed is obtained from a call to 'system_clock()'. c c If this function is unavailable, another source for the seed must c c be supplied by the user. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none integer nz0,nz1,ns,num,ncol,max,ishb parameter(nz0=2,nz1=8,ns=18,num=100,ncol=8,ishb=11,max=1000000) integer maxn(nz0,ns),iage,h,n,na,nm,ifeh,iafe,ni(nz1,6),ni0,im, * cols(ncol),iseed,i,j,k,ii,jj,kk,iclr,nstar,nz character track*55,outfile*21,band(8)*8,arg*10 double precision a(nz0,ns,num),t(nz0,ns,num),l(nz0,ns,num),rn,ave, * g(nz0,ns,num),m(max),dev,sig,ran2,gasdev,qa(4),fa(4),qm(4), * fm(4),a_ran,l_ran,t_ran,g_ran,f_ran(ncol),feh0(nz1),afe0(6), * mhb,feh,mrgb,fltr(nz0,ns,num,ncol),qz(nz0),fz(nz0),mstar(ns), * ar(nz0,4),tr(nz0,4),gr(nz0,4),lr(nz0,4),fr(nz0,4,ncol) logical lquiet c 'ni' lists the initial mass of each [Fe/H],[alpha/Fe] data ni/7,6,5,5,4,3,3,2, ![a/Fe]=-0.2 * 6,6,5,5,4,3,2,1, ! = 0 * 6,6,5,5,4,3,2,1, ! =+0.2 * 6,6,5,4,4,2,ns,ns, ! =+0.4 * 6,5,5,4,4,2,ns,ns, ! =+0.6 * 5,5,5,4,3,1,ns,ns/ ! =+0.8 data cols/4,6,8,8,8,5,7,7/, fz/nz0*0.0d0/, qz/nz0*0.0d0/ data mstar/4.6d-1,4.7d-1,4.8d-1,4.9d-1,5.0d-1,5.2d-1,5.4d-1, * 5.6d-1,5.8d-1,6.0d-1,6.5d-1,7.0d-1,8.0d-1,9.0d-1,1.0d0, * 1.25d0,1.5d0,1.75d0/ data fa/4*0.0d0/, qa/4*0.0d0/, fm/4*0.0d0/, qm/4*0.0d0/ 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/, lquiet/.false./ C run in quiet mode: get arguments from the cmd line and skip interactive parts if(IArgC().gt.0) then lquiet=.true. call getarg(1,arg) read(arg,'(i1)') iclr call getarg(2,arg) read(arg,'(i1)') iafe call getarg(3,arg) read(arg,'(f10.5)') feh call getarg(4,arg) read(arg,'(i6)') nstar call getarg(5,arg) read(arg,'(f10.5)') mrgb call getarg(6,arg) read(arg,'(f10.5)') mhb call getarg(7,arg) read(arg,'(f10.5)') sig call getarg(8,outfile) goto 1 endif c begin interactive input section print *, "-----------------------------------------------" print *, " Synthetic HB Generator" print *, "-----------------------------------------------" print *, ' Enter choice of color-Teff transformation: [1-8]' 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' read *, iclr if(iclr.lt.1.or.iclr.gt.8) stop 'INVALID CHOICE OF COLORS' 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]' print *, '-----------------------------------------------' print *, ' Available [Fe/H] range is ',feh0(1),'->',feh0(nz1) print *, ' Enter new [Fe/H] to interpolate:' read *, feh if(feh.lt.feh0(1).or.feh.gt.feh0(nz1)) * print *, 'EXTRAPOLATION WARNING - [Fe/H] out of bounds' if(feh.gt.0.0d0.and.iafe.gt.3) * stop 'INVALID COMBINATION OF [Fe/H] and [alpha/Fe]' c set up mass distribution print *, "-----------------------------------------------" print *, "Enter: (1) Number of stars [<10^6]" print *, " (2) Mass of stars at RGB tip" print *, " (3) Average mass on HB" print *, " (4) HB Mass Dispersion " print *, "-----------------------------------------------" read *, nstar, mrgb, mhb, sig print *, "-----------------------------------------------" print *, "Enter output file name:" read *, outfile c end of interactive input section 1 continue c locate input [Fe/H] in feh0 array do i=1,nz1 if(feh.gt.feh0(i))then ifeh=i nz=nz0 elseif(feh.eq.feh0(i))then ifeh=i nz=1 endif enddo do i=1,nz qz(i)=feh0(ifeh+i-(nz+1)/2) enddo c interpolate in [Fe/H] call interp(qz,fz,feh,nz) c if using more than one [Fe/H], use the highest minimum mass ni0=1 do i=1,nz k=ni(ifeh+i-(nz+1)/2,iafe) if(k.gt.ni0) ni0=k enddo c read in HB tracks for all [Fe/H]'s and masses do k=1,nz kk=ifeh+k-1 do j=ni0,ns open(ishb,file=track(kk,iafe,iclr,j),status='old') read(ishb,*) read(ishb,'(39x,8a8)') (band(i),i=1,cols(iclr)) do i=1,num read(ishb,*,end=10) a(k,j,i),t(k,j,i),g(k,j,i),l(k,j,i), * (fltr(k,j,i,h),h=1,cols(iclr)) enddo !i-loop 10 close(ishb) maxn(k,j)=i-1 enddo !j-loop enddo !k-loop c initialize some variables for the mass probability function if(mrgb.gt.mstar(ns)) mrgb=mstar(ns) if(mhb.gt.mrgb) mhb=mrgb c initialize RNG for random mass distribution call system_clock(iseed) rn=ran2(-iseed) c process NSTARs of random mass drawn from a Gaussian distribution c but ensure the masses lie within the specified boundaries i=0 do while (i.le.nstar) rn=gasdev(iseed,mhb,sig) if(rn.le.mrgb.and.rn.ge.mstar(ni0)) then i=i+1 m(i)=rn endif enddo c calculate some statistics regarding the mass distribution c because of the imposed boundaries the actual parameters may c differ from the input values ave=0.0d0 do i=1,nstar ave=ave+m(i) enddo ave=ave/dble(nstar) dev=0.0d0 do i=1,nstar dev=dev+(m(i)-ave)*(m(i)-ave) enddo dev=dsqrt(dev/dble(nstar)) c interactive: print results of mass distribution if(.not.lquiet) then print *, "-----------------------------------------------" print *, " Results: Ave. Mass=",ave print *, " Std. Dev.=",dev print *, "-----------------------------------------------" endif c re-initialize RNG for ages call system_clock(iseed) rn=ran2(-iseed) c open file for writing, write header lines open(unit=ishb,file=outfile,status='unknown') write(ishb,'(a3,i8,a10,f5.2,a10,f4.2,a13,f5.3,a13,f5.3)') * '#N=',nstar,' [Fe/H]=',feh,' [a/Fe]=',afe0(iafe), * ' Ave. Mass=',ave,' Std. Dev.=',dev write(ishb,'(a8,a12,11a8)') '# M/Mo ',' Age (yrs) ','Log Teff', * ' Log g ','Log L/Lo ', (band(i),i=1,cols(iclr)) c loop through mass distribution, select random age and interpolate tracks do i=1,nstar do k=ni0,ns-1 !locate random mass in mstar array if(m(i).ge.mstar(k).and.m(i).lt.mstar(k+1)) im=k enddo nm=2 !get mass interpolation coefficients if(im.lt.ni0+1) nm=2 if(im.gt.ns-2) nm=2 do k=1,nm qm(k)=mstar(im+k-nm/2) enddo call interp(qm,fm,m(i),nm) rn=ran2(iseed) !for random age do k=1,nz !loop over [Fe/H] kk=ifeh+k-(nz+1)/2 do j=1,nm !loop over masses jj=im+j-nm/2 ar(k,j)=a(k,jj,1)+rn*(a(k,jj,maxn(k,jj))-a(k,jj,1)) do n=1,maxn(k,jj)-1 !locate random age in age arrays if(ar(k,j).ge.a(k,jj,n).and.ar(k,j).lt.a(k,jj,n+1)) * iage=n enddo na=2 !age interpolation coefficients do n=1,na ii=iage+n-na/2 qa(n)=a(k,jj,ii) enddo call interp(qa,fa,ar(k,j),na) ar(k,j)=0.0d0 tr(k,j)=0.0d0 gr(k,j)=0.0d0 lr(k,j)=0.0d0 do n=1,cols(iclr) fr(k,j,n)=0.0d0 enddo do n=1,na ! perform age-interpolation ii=iage+n-na/2 ar(k,j) = ar(k,j) + fa(n)*a(k,jj,ii) tr(k,j) = tr(k,j) + fa(n)*t(k,jj,ii) gr(k,j) = gr(k,j) + fa(n)*g(k,jj,ii) lr(k,j) = lr(k,j) + fa(n)*l(k,jj,ii) do h=1,cols(iclr) fr(k,j,h)=fr(k,j,h) + fa(n)*fltr(k,jj,ii,h) enddo !h-loop enddo !n-loop enddo !j-loop enddo !k-loop a_ran=0.0d0 !zero results before interpolation t_ran=0.0d0 g_ran=0.0d0 l_ran=0.0d0 do h=1,cols(iclr) f_ran(h)=0.0d0 enddo do k=1,nz ![Fe/H]-interpolation do j=1,nm !mass-interpolation a_ran=a_ran + fz(k)*fm(j)*ar(k,j) t_ran=t_ran + fz(k)*fm(j)*tr(k,j) g_ran=g_ran + fz(k)*fm(j)*gr(k,j) l_ran=l_ran + fz(k)*fm(j)*lr(k,j) do h=1,cols(iclr) f_ran(h)=f_ran(h) + fz(k)*fm(j)*fr(k,j,h) enddo !h-loop enddo !j-loop enddo !k-loop c write results to file, but only if this is a valid point if(a_ran.eq.a_ran) write(ishb,'(0pf8.4,1p1e12.4,0p11f8.4)') * m(i),a_ran,t_ran,g_ran,l_ran,(f_ran(h),h=1,cols(iclr)) enddo !i-loop c close output file and, if interactive mode, print final message close(ishb) if(.not.lquiet)then print *, " Results written to ", outfile print *, "-----------------------------------------------" endif end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function ran2(idum) c random double precision number on [0,1] integer idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv double precision ran2,am,eps,rnmx parameter(im1=2147483563,im2=2147483399,am=1.0d0/im1,imm1=im1-1, * ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,ir2=3791, * ntab=32,ndiv=1+imm1/ntab,eps=1.2d-7,rnmx=1.0d0-eps) integer idum2,j,k,iv(ntab),iy save iv,iy,idum2 data idum2/123456789/, iv/ntab*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do j=ntab+8,1,-1 k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 if (idum.lt.0) idum=idum+im1 if (j.le.ntab) iv(j)=idum enddo iy=iv(1) endif k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 if (idum.lt.0) idum=idum+im1 k=idum2/iq2 idum2=ia2*(idum2-k*iq2)-k*ir2 if (idum2.lt.0) idum2=idum2+im2 j=1+iy/ndiv iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+imm1 ran2=min(am*iy,rnmx) return end C (C) Copr. 1986-92 Numerical Recipes Software ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function gasdev(idum,mean,var) c adapted from numerical recipes function of same name c to include arbitrary mean and variance for gaussian dist. c price is that we do not bother to save one of the dist. c numbers (gset) c generated since it is presumed that the mean and/or c variance will likely change between calls integer idum double precision gasdev,mean,var,rsq,v1,v2,ran2,var2 var2 = var*var 1 v1=var*(2.0d0*ran2(idum)-1.0d0) v2=var*(2.0d0*ran2(idum)-1.0d0) rsq=v1*v1+v2*v2 if (rsq.ge.var2.or.rsq.eq.0.0d0) goto 1 gasdev=v2*sqrt(-2.0d0*log(rsq/var2)/(rsq/var2)) + mean return end c----------------------------------------------------------------------- subroutine interp(a,b,x,n) c {a} are the tabulated values for use in interpolation c {b} are coefficients of the interpolating polynomial 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 track(feh,afe,clr,mass) c returns the name of track file, not pretty but it works! integer i1, i2, ns, nz, ncol parameter(ns=18,nz=8,ncol=8) integer feh, afe, clr, mass, dirl(ncol), sufl(ncol) character cmass(ns)*3,cfeh(nz)*3,cafe(6)*2,cdir(ncol)*9, * csuf(ncol)*9,track*55,comp*11 data cdir/'bvri','uvby','jc2mass','hst_acs','hst_wfpc2','sdss', * 'gct_emp','gct_syn'/ data csuf/'bvri','uvby','jc2mass','hst_acs','hst_wfpc2','ugriz', * 'cmd','phx'/ data cmass/'046','047','048','049','050','052','054','056','058', * '060','065','070','080','090','100','125','150','175'/ data cfeh/'m25','m20','m15','m10','m05','p00','p03','p05'/ data cafe/'m2','p0','p2','p4','p6','p8'/ data dirl/4,4,7,7,9,4,7,7/ data sufl/4,4,7,7,9,5,3,3/ data comp/'fehxxxafexx'/ do i1=1,len(track) track(i1:i1)=' ' enddo comp(4:6)=cfeh(feh) comp(10:11)=cafe(afe) i1=1 i2=7 track(i1:i2)='tracks/' i1=i2+1 i2=i2+dirl(clr) track(i1:i2)=cdir(clr) i1=i2+1 i2=i1 track(i1:i2)='/' i1=i2+1 i2=i2+len(comp) track(i1:i2)=comp i1=i2+1 i2=i1 track(i1:i2)='/' i1=i2+1 i2=i1+1 track(i1:i2)='hb' i1=i2+1 i2=i2+3 track(i1:i2)=cmass(mass) i1=i2+1 i2=i2+len(comp) track(i1:i2)=comp i1=i2+1 i2=i1 track(i1:i2)='.' i1=i2+1 i2=i2+sufl(clr) track(i1:i2)=csuf(clr) return end c------------------------------------------------------------------------ c------------------------- END SHB PROGRAM ------------------------------ c------------------------------------------------------------------------