program abundance implicit none real*8 M200,alpha,beta_mlcs,beta_salt,gamma,z real*8 n_mlcs,n_salt,np_mlcs,np_salt,n real*8 dinvnorm,p,fsky,pscale,nlim,mlim write(6,*) 'M_200 (wrt mean in h^-1 Mpc)?' read*,M200 write(6,*) 'Redshift z?' read*,z write(6,*) 'Survey Area (sq degrees)?' read*,fsky if (fsky.le.1) then write(6,*) 'Caution: input is in sq degrees' write(6,*) ' do you mean',180d0**2*4d0/3.1415926*fsky end if fsky = fsky*3.1415926/4d0/(180d0)**2 write(6,"(' ',1E15.5,' sky fraction')") fsky write(6,*) 'Parameter CL or P level (p=0.95, 95% CL)' read*,p if (p.lt.0) then write(6,*) ' p level unphysical' stop end if pscale=dinvnorm(p) c LCDM median numbers for MLCS2 and SALT SN fitters as of 2010 n_mlcs=n(dlog10(M200),alpha(z),beta_mlcs(z)) np_mlcs=(1d0-0.035*pscale)*n_mlcs + 0.29*pscale n_salt=n(dlog10(M200),alpha(z),beta_salt(z)) np_salt=(1d0-0.035*pscale)*n_salt + 0.29*pscale c calculate the 100p% parameter and 95% sample CL mass nlim = dlog10(-dlog(0.95d0)/fsky)-0.29*pscale nlim = nlim/(1d0-0.035*pscale) mlim = 1d0/alpha(z)* dlog(1d0-nlim/7.65) write(6,*) write(6,'(1I3,"% CL LCDM Number >M,>z ")') Int(100*p) write(6,"(1F15.5,' =N for MLCSk2')") c fsky*(10d0**np_mlcs) write(6,"(' ',1F5.2,'% CL sample exclusion')") c 100d0*dexp(-fsky*(10d0**np_mlcs)) write(6,"(' ',1E11.5,' Ms/h for 95% CL sample exclusion')") c 10d0**(beta_mlcs(z)+mlim) write(6,"(' ',1F4.1,' local mfn slope at M,z')") c gamma(np_mlcs,z) write(6,"(1F15.5,' =N for SALT2 ')") c fsky*(10d0**np_salt) write(6,"(' ',1F5.2,'% CL sample exclusion')") c 100d0*dexp(-fsky*(10d0**np_salt)) write(6,"(' ',1E11.5,' Ms/h for 95% CL sample exclusion')") c 10d0**(beta_salt(z)+mlim) write(6,"(' ',1F4.1,' local mfn slope at M,z')") c gamma(np_salt,z) c calculate the 100p% parameter and 95% sample CL mass nlim = dlog10(-dlog(0.95d0)/fsky) mlim = 1d0/alpha(z)* dlog(1d0-nlim/7.65) write(6,*) write(6,'(1I3,"% CL LCDM Number >M,>z ")') Int(100*0.5) write(6,"(1F15.5,' =N for MLCSk2')") c fsky*(10d0**n_mlcs) write(6,"(' ',1F5.2,'% CL sample exclusion')") c 100d0*dexp(-fsky*(10d0**n_mlcs)) write(6,"(' ',1E11.5,' Ms/h for 95% CL sample exclusion')") c 10d0**(beta_mlcs(z)+mlim) write(6,"(' ',1F4.1,' local mfn slope at M,z')") c gamma(n_mlcs,z) write(6,"(1F15.5,' =N for SALT2 ')") c fsky*(10d0**n_salt) write(6,"(' ',1F5.2,'% CL sample exclusion')") c 100d0*dexp(-fsky*(10d0**n_salt)) write(6,"(' ',1E11.5,' Ms/h for 95% CL sample exclusion')") c 10d0**(beta_salt(z)+mlim) write(6,"(' ',1F4.1,' local mfn slope at M,z')") c gamma(n_salt,z) end c------------------------------------------------- real*8 function n(lgM,alpha,beta) real*8 alpha,beta,lgM n = dexp(alpha*(lgM-beta)) n = 7.65*(1d0-n) return end c------------------------------------------------- c------------------------------------------------- real*8 function alpha(z) real*8 z alpha = 1.06-0.17*dexp(-1.3*z) return end c------------------------------------------------- c------------------------------------------------- real*8 function beta_mlcs(z) real*8 z beta_mlcs =15.565-0.1*dlog10(7.1+10**(5.25*z)) return end c------------------------------------------------- c------------------------------------------------- real*8 function beta_salt(z) real*8 z beta_salt =15.525-0.1*dlog10(7.1+10**(5.25*z)) return end c------------------------------------------------- c------------------------------------------------- real*8 function gamma(lgN,z) real*8 lgN,z gamma = 7.1 - 1.5*dexp(-3d0*z)-1.1*lgN gamma = dexp(gamma)+2.6+1.5*z**2 gamma = -dlog(gamma) return end c------------------------------------------------- c ren-raw chen, rutgers business school c normal inverse c translate from http://home.online.no/~pjacklam/notes/invnorm c a routine written by john herrero c this returns sqrt[2]*InvErf(2p-1) real*8 function dinvnorm(p) real*8 p,p_low,p_high real*8 a1,a2,a3,a4,a5,a6 real*8 b1,b2,b3,b4,b5 real*8 c1,c2,c3,c4,c5,c6 real*8 d1,d2,d3,d4 real*8 z,q,r a1=-39.6968302866538 a2=220.946098424521 a3=-275.928510446969 a4=138.357751867269 a5=-30.6647980661472 a6=2.50662827745924 b1=-54.4760987982241 b2=161.585836858041 b3=-155.698979859887 b4=66.8013118877197 b5=-13.2806815528857 c1=-0.00778489400243029 c2=-0.322396458041136 c3=-2.40075827716184 c4=-2.54973253934373 c5=4.37466414146497 c6=2.93816398269878 d1=0.00778469570904146 d2=0.32246712907004 d3=2.445134137143 d4=3.75440866190742 p_low=0.02425 p_high=1-p_low if(p.lt.p_low) goto 201 if(p.ge.p_low) goto 301 201 q=dsqrt(-2*dlog(p)) z=(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/ &((((d1*q+d2)*q+d3)*q+d4)*q+1) goto 204 301 if((p.ge.p_low).and.(p.le.p_high)) goto 202 if(p.gt.p_high) goto 302 202 q=p-0.5 r=q*q z=(((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q/ &(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1) goto 204 302 if((p.gt.p_high).and.(p.lt.1)) goto 203 203 q=dsqrt(-2*dlog(1-p)) z=-(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/ &((((d1*q+d2)*q+d3)*q+d4)*q+1) 204 dinvnorm=z return end