Index: trunk/Cosmo/SimLSS/Makefile
===================================================================
--- trunk/Cosmo/SimLSS/Makefile	(revision 3314)
+++ trunk/Cosmo/SimLSS/Makefile	(revision 3315)
@@ -37,4 +37,11 @@
 LIBR = $(LIB)libcmvsimbao.a
 
+PROGSTEST = \
+        $(EXE)cmvtluc $(EXE)cmvchkwhu $(EXE)hu_sigma8
+
+PROGSTESTOBJ = \
+        $(OBJ)hu_tf_fit.o $(OBJ)hu_power.o \
+        $(OBJ)cmvtluc.o $(OBJ)cmvchkwhu.o $(OBJ)hu_sigma8.o
+
 #----
 all: lib prog
@@ -44,9 +51,12 @@
 prog: $(PROGS)
 
+progtest: $(PROGSTEST)
+
 clean:
 	rm -f *.o
-	rm -f $(PROGS) $(LIBR)
-	rm -f $(PROGSOBJ) $(LIBROBJ)
+	rm -f $(PROGS) $(PROGSOBJ)
+	rm -f $(LIBR) $(LIBROBJ)
 	rm -rf $(OBJ)/CmvBAO_cxxrep/
+	rm -f $(PROGSTEST) $(PROGSTESTOBJ)
 
 ##############################################################################
@@ -72,5 +82,5 @@
 cmvtuniv: $(EXE)cmvtuniv
 	echo $@ " done"
-$(EXE)cmvtuniv: $(OBJ)cmvtuniv.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtuniv: $(OBJ)cmvtuniv.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtuniv.o $(MYLIB) 
 $(OBJ)cmvtuniv.o: cmvtuniv.cc
@@ -80,5 +90,5 @@
 cmvtransf: $(EXE)cmvtransf
 	echo $@ " done"
-$(EXE)cmvtransf: $(OBJ)cmvtransf.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtransf: $(OBJ)cmvtransf.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtransf.o $(MYLIB) 
 $(OBJ)cmvtransf.o: cmvtransf.cc
@@ -88,5 +98,5 @@
 cmvtgrowth: $(EXE)cmvtgrowth
 	echo $@ " done"
-$(EXE)cmvtgrowth: $(OBJ)cmvtgrowth.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtgrowth: $(OBJ)cmvtgrowth.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtgrowth.o $(MYLIB) 
 $(OBJ)cmvtgrowth.o: cmvtgrowth.cc
@@ -96,5 +106,5 @@
 cmvtstpk: $(EXE)cmvtstpk
 	echo $@ " done"
-$(EXE)cmvtstpk: $(OBJ)cmvtstpk.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtstpk: $(OBJ)cmvtstpk.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstpk.o $(MYLIB) 
 $(OBJ)cmvtstpk.o: cmvtstpk.cc
@@ -104,5 +114,5 @@
 cmvtstsch: $(EXE)cmvtstsch
 	echo $@ " done"
-$(EXE)cmvtstsch: $(OBJ)cmvtstsch.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtstsch: $(OBJ)cmvtstsch.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstsch.o $(MYLIB) 
 $(OBJ)cmvtstsch.o: cmvtstsch.cc
@@ -112,5 +122,5 @@
 cmvtstblack: $(EXE)cmvtstblack
 	echo $@ " done"
-$(EXE)cmvtstblack: $(OBJ)cmvtstblack.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtstblack: $(OBJ)cmvtstblack.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstblack.o $(MYLIB) 
 $(OBJ)cmvtstblack.o: cmvtstblack.cc
@@ -120,5 +130,5 @@
 cmvtvarspec: $(EXE)cmvtvarspec
 	echo $@ " done"
-$(EXE)cmvtvarspec: $(OBJ)cmvtvarspec.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtvarspec: $(OBJ)cmvtvarspec.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtvarspec.o $(MYLIB) 
 $(OBJ)cmvtvarspec.o: cmvtvarspec.cc
@@ -128,5 +138,5 @@
 cmvdefsurv: $(EXE)cmvdefsurv
 	echo $@ " done"
-$(EXE)cmvdefsurv: $(OBJ)cmvdefsurv.o $(LIB)libcmvsimbao.a
+$(EXE)cmvdefsurv: $(OBJ)cmvdefsurv.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvdefsurv.o $(MYLIB) 
 $(OBJ)cmvdefsurv.o: cmvdefsurv.cc
@@ -136,5 +146,5 @@
 cmvobserv3d: $(EXE)cmvobserv3d
 	echo $@ " done"
-$(EXE)cmvobserv3d: $(OBJ)cmvobserv3d.o $(LIB)libcmvsimbao.a
+$(EXE)cmvobserv3d: $(OBJ)cmvobserv3d.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvobserv3d.o $(MYLIB) 
 $(OBJ)cmvobserv3d.o: cmvobserv3d.cc
@@ -144,5 +154,5 @@
 cmvtintfun: $(EXE)cmvtintfun
 	echo $@ " done"
-$(EXE)cmvtintfun: $(OBJ)cmvtintfun.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtintfun: $(OBJ)cmvtintfun.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtintfun.o $(MYLIB) 
 $(OBJ)cmvtintfun.o: cmvtintfun.cc
@@ -152,5 +162,5 @@
 cmvconcherr: $(EXE)cmvconcherr
 	echo $@ " done"
-$(EXE)cmvconcherr: $(OBJ)cmvconcherr.o $(LIB)libcmvsimbao.a
+$(EXE)cmvconcherr: $(OBJ)cmvconcherr.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvconcherr.o $(MYLIB) 
 $(OBJ)cmvconcherr.o: cmvconcherr.cc
@@ -161,5 +171,5 @@
 cmvtpoisson: $(EXE)cmvtpoisson
 	echo $@ " done"
-$(EXE)cmvtpoisson: $(OBJ)cmvtpoisson.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtpoisson: $(OBJ)cmvtpoisson.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtpoisson.o $(MYLIB) 
 $(OBJ)cmvtpoisson.o: cmvtpoisson.cc
@@ -169,5 +179,5 @@
 cmvtinterp: $(EXE)cmvtinterp
 	echo $@ " done"
-$(EXE)cmvtinterp: $(OBJ)cmvtinterp.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtinterp: $(OBJ)cmvtinterp.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtinterp.o $(MYLIB) 
 $(OBJ)cmvtinterp.o: cmvtinterp.cc
@@ -177,5 +187,5 @@
 cmvtstagn: $(EXE)cmvtstagn
 	echo $@ " done"
-$(EXE)cmvtstagn: $(OBJ)cmvtstagn.o $(LIB)libcmvsimbao.a
+$(EXE)cmvtstagn: $(OBJ)cmvtstagn.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtstagn.o $(MYLIB) 
 $(OBJ)cmvtstagn.o: cmvtstagn.cc
@@ -185,5 +195,5 @@
 cmvschdist: $(EXE)cmvschdist
 	echo $@ " done"
-$(EXE)cmvschdist: $(OBJ)cmvschdist.o $(LIB)libcmvsimbao.a
+$(EXE)cmvschdist: $(OBJ)cmvschdist.o $(LIBR)
 	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvschdist.o $(MYLIB) 
 $(OBJ)cmvschdist.o: cmvschdist.cc
@@ -191,4 +201,22 @@
 
 
+##############################################################################
+cmvchkwhu: $(EXE)cmvchkwhu
+	echo $@ " done"
+$(EXE)cmvchkwhu: $(OBJ)cmvchkwhu.o $(OBJ)hu_tf_fit.o $(OBJ)hu_power.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvchkwhu.o $(OBJ)hu_tf_fit.o $(OBJ)hu_power.o $(MYLIB) 
+$(OBJ)cmvchkwhu.o: cmvchkwhu.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvchkwhu.cc
+$(OBJ)hu_tf_fit.o: hu_tf_fit.c hu_tf_fit.h
+	$(CCOMPILE) -o $@ hu_tf_fit.c
+$(OBJ)hu_power.o: hu_power.c hu_power.h
+	$(CCOMPILE) -o $@ hu_power.c
+
+hu_sigma8: $(EXE)hu_sigma8
+	echo $@ " done"
+$(EXE)hu_sigma8: $(OBJ)hu_sigma8.o $(OBJ)hu_power.o
+	$(CLINK) -o $@ $(OBJ)hu_sigma8.o $(OBJ)hu_power.o -lc -lm
+$(OBJ)hu_sigma8.o: hu_sigma8.c
+	$(CCOMPILE) -o $@ hu_sigma8.c
 
 ##############################################################################
Index: trunk/Cosmo/SimLSS/cmvchkwhu.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvchkwhu.cc	(revision 3315)
+++ trunk/Cosmo/SimLSS/cmvchkwhu.cc	(revision 3315)
@@ -0,0 +1,157 @@
+// Check cmv / Wayne Hu coding of the transfert functions
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include "ntuple.h"
+
+#include "constcosmo.h"
+#include "pkspectrum.h"
+#include "hu_tf_fit.h"
+#include "hu_power.h"
+
+int main(int narg,char *arg[])
+{
+ double k1 = 1e-6, k2 = 10.; int_4 npt = 10000;
+ double h100 = 0.71, Om0 = 0.267804, Ob0 = 0.0444356, tcmb = T_CMB_Par;
+ // Modele MDM (hu_power.c)  Eisenstein & Hu Apj 511 p5-15, 1999 (ou astro-ph/9710252)
+ double Ohdm = 0.0, Ol0 = 0.7, redshift=0.;
+ int degen_hdm = 0;
+
+ cout<<"k1="<<k1<<"  k2="<<k2<<" npt="<<npt<<endl;
+ cout<<"h100="<<h100<<" Om0="<<Om0<<" Ob0="<<Ob0<<" Tcmb="<<tcmb<<endl;
+ cout<<"Om0*h^2="<<Om0*h100*h100<<" Ob0*h^2="<<Ob0*h100*h100<<endl;
+ double fbaryon = Ob0/Om0;
+ cout<<"fbaryon="<<fbaryon<<endl;
+
+ cout<<endl<<"TFset_parameters"<<endl;
+ TFset_parameters(Om0*h100*h100,fbaryon,tcmb);
+ TFprint_parameters();
+ double soundfithu = TFsound_horizon_fit(Om0,fbaryon,h100)/h100;
+ cout<<"TFsound_horizon_fit = "<<soundfithu<<" Mpc"<<endl;
+
+ cout<<endl<<"TransfertEisenstein with baryon"<<endl;
+ TransfertEisenstein Tfull(h100,Om0-Ob0,Ob0,tcmb,false,1);
+
+ cout<<endl<<"TransfertEisenstein with baryon without oscillation approx.2"<<endl;
+ TransfertEisenstein Tnowig(h100,Om0-Ob0,Ob0,tcmb,false,1);
+ Tnowig.SetNoOscEnv(2);
+
+ cout<<endl<<"TransfertEisenstein without baryon"<<endl;
+ TransfertEisenstein Tnob(h100,Om0,0.,tcmb,true,1);
+
+ cout<<endl<<"TFmdm_set_cosm"<<endl;
+ int rcmdm = TFmdm_set_cosm(Om0,Ob0,Ohdm,degen_hdm,Ol0,h100,redshift);
+ cout<<"  rc="<<rcmdm<<endl;
+
+ cout<<endl<<"Compare kpeak"<<endl;
+ double kpeak = Tfull.KPeak();
+ double kpeakhu = TFk_peak(Om0,fbaryon,h100)*h100;
+ cout<<"kpeak = "<<kpeak<<" "<<kpeakhu<<"  Mpc^-1  -> diff "<<kpeak-kpeakhu<<endl;
+
+ const int n = 12;
+ char *vname[n] = {"k","tf","tfnw","tfnb","tfc","tfb"
+                  ,"hutf","hutfnw","hutfnb","hutfc","hutfb"
+                  ,"hutfmdm"};
+ NTuple nt(n,vname);
+ double xnt[n];
+
+ double kmax[5]={0,0,0,0,0};
+ double tmax[5]={0.,0.,0.,0.,0.};
+ double lnk1 = log10(k1), lnk2=log10(k2), dlnk=(lnk2-lnk1)/npt;
+ for(double lnk=lnk1;lnk<lnk2+dlnk/2.;lnk+=dlnk) {
+   double k = pow(10.,lnk);
+   //if(fabs(lnk-log10(0.0186209))>1e-5) continue;
+   Tfull.SetReturnPart(0); double tf = Tfull(k);
+   double tfnw = Tnowig(k);
+   double tfnb = Tnob(k);
+   Tfull.SetReturnPart(1); double tfc = Tfull(k);
+   Tfull.SetReturnPart(2); double tfb = Tfull(k);
+   double hutfc, hutfb;
+   double hutf = TFfit_onek(k,&hutfb,&hutfc); 
+   double hutfnw = TFnowiggles(Om0,fbaryon,h100,tcmb,k/h100);
+   double hutfnb = TFzerobaryon(Om0,h100,tcmb,k/h100);
+   double hutfmdm = TFmdm_onek_mpc(k);
+
+
+   if(fabs(tf-hutf)>tmax[0])     {kmax[0]=k; tmax[0]=fabs(tf-hutf);}
+   if(fabs(tfnw-hutfnw)>tmax[1]) {kmax[1]=k; tmax[1]=fabs(tfnw-hutfnw);}
+   if(fabs(tfnb-hutfnb)>tmax[2]) {kmax[2]=k; tmax[2]=fabs(tfnb-hutfnb);}
+   if(fabs(tfc-hutfc)>tmax[3])   {kmax[3]=k; tmax[3]=fabs(tfc-hutfc);}
+   if(fabs(tfb-hutfb)>tmax[4])   {kmax[4]=k; tmax[4]=fabs(tfb-hutfb);}
+
+   xnt[0]=k;
+   xnt[1]=tf; xnt[2]=tfnw; xnt[3]=tfnb; xnt[4]=tfc; xnt[5]=tfb;
+   xnt[6]=hutf; xnt[7]=hutfnw; xnt[8]=hutfnb; xnt[9]=hutfc; xnt[10]=hutfb;
+   xnt[11] = hutfmdm;
+   nt.Fill(xnt);
+ }
+
+ cout<<"Maximum difference for:"<<endl;
+ cout<<"tf:   "<<tmax[0]<<"  for k="<<kmax[0]<<endl;
+ cout<<"tfnw: "<<tmax[1]<<"  for k="<<kmax[1]<<endl;
+ cout<<"tfnb: "<<tmax[2]<<"  for k="<<kmax[2]<<endl;
+ cout<<"tfc:  "<<tmax[3]<<"  for k="<<kmax[3]<<endl;
+ cout<<"tfb:  "<<tmax[4]<<"  for k="<<kmax[4]<<endl;
+
+ cout<<">>>> Ecriture"<<endl;
+ string tag = "cmvchkwhu.ppf";
+ POutPersist pos(tag);
+ tag = "nt"; pos.PutObject(nt,tag);
+
+ return 0;
+}
+
+/*
+openppf cmvchkwhu.ppf
+
+set k k
+set k log10(k)
+
+######################
+zone 3 1
+n/plot nt.tf%$k ! ! "nsta connectpoints"
+n/plot nt.hutf%$k ! ! "nsta connectpoints red same"
+n/plot nt.tfc%$k ! ! "nsta connectpoints"
+n/plot nt.hutfc%$k ! ! "nsta connectpoints red same"
+n/plot nt.tfb%$k ! ! "nsta connectpoints"
+n/plot nt.hutfb%$k ! ! "nsta connectpoints red same"
+zone 2 1
+n/plot nt.tfnw%$k ! ! "nsta connectpoints"
+n/plot nt.hutfnw%$k ! ! "nsta connectpoints red same"
+n/plot nt.tfnb%$k ! ! "nsta connectpoints"
+n/plot nt.hutfnb%$k ! ! "nsta connectpoints red same"
+
+######################
+zone 3 1
+n/plot nt.tf-hutf%$k ! ! "nsta connectpoints"
+n/plot nt.tfc-hutfc%$k ! ! "nsta connectpoints"
+n/plot nt.tfb-hutfb%$k ! ! "nsta connectpoints"
+zone 2 1
+n/plot nt.tfnw-hutfnw%$k ! ! "nsta connectpoints"
+n/plot nt.tfnb-hutfnb%$k ! ! "nsta connectpoints"
+
+zone 3 1
+n/plot nt.(tf-hutf)/hutf%$k ! ! "nsta connectpoints"
+n/plot nt.(tfc-hutfc)/hutfc%$k ! ! "nsta connectpoints"
+n/plot nt.(tfb-hutfb)/hutfb%$k ! ! "nsta connectpoints"
+zone 2 1
+n/plot nt.(tfnw-hutfnw)/hutfnw%$k ! ! "nsta connectpoints"
+n/plot nt.(tfnb-hutfnb)/hutfnb%$k ! ! "nsta connectpoints"
+
+######################
+zone
+n/plot nt.tf%$k ! ! "nsta connectpoints"
+n/plot nt.tfnw%$k ! ! "nsta connectpoints red same"
+n/plot nt.tfnb%$k ! ! "nsta connectpoints blue same"
+n/plot nt.hutfmdm%$k ! ! "nsta connectpoints green same"
+
+zone
+n/plot nt.tf/tfnw%$k tfnw>0. ! "nsta connectpoints"
+n/plot nt.hutf/hutfnw%$k hutfnw>0. ! "nsta connectpoints same blue"
+addline -1000 1 1000 1 "red"
+*/
Index: trunk/Cosmo/SimLSS/hu_power.c
===================================================================
--- trunk/Cosmo/SimLSS/hu_power.c	(revision 3315)
+++ trunk/Cosmo/SimLSS/hu_power.c	(revision 3315)
@@ -0,0 +1,274 @@
+/* Fitting Formulae for CDM + Baryon + Massive Neutrino (MDM) cosmologies. */
+/* Daniel J. Eisenstein & Wayne Hu, Institute for Advanced Study */
+
+/* There are two primary routines here, one to set the cosmology, the
+other to construct the transfer function for a single wavenumber k. 
+You should call the former once (per cosmology) and the latter as 
+many times as you want. */
+
+/* TFmdm_set_cosm() -- User passes all the cosmological parameters as
+	arguments; the routine sets up all of the scalar quantites needed 
+	computation of the fitting formula.  The input parameters are: 
+	1) omega_matter -- Density of CDM, baryons, and massive neutrinos,
+				in units of the critical density. 
+	2) omega_baryon -- Density of baryons, in units of critical. 
+	3) omega_hdm    -- Density of massive neutrinos, in units of critical 
+	4) degen_hdm    -- (Int) Number of degenerate massive neutrino species 
+	5) omega_lambda -- Cosmological constant 
+	6) hubble       -- Hubble constant, in units of 100 km/s/Mpc 
+	7) redshift     -- The redshift at which to evaluate */
+
+/* TFmdm_onek_mpc() -- User passes a single wavenumber, in units of Mpc^-1.
+	Routine returns the transfer function from the Eisenstein & Hu
+	fitting formula, based on the cosmology currently held in the 
+	internal variables.  The routine returns T_cb (the CDM+Baryon
+	density-weighted transfer function), although T_cbn (the CDM+
+	Baryon+Neutrino density-weighted transfer function) is stored
+	in the global variable tf_cbnu. */
+
+/* We also supply TFmdm_onek_hmpc(), which is identical to the previous
+	routine, but takes the wavenumber in units of h Mpc^-1. */
+
+/* We hold the internal scalar quantities in global variables, so that
+the user may access them in an external program, via "extern" declarations. */
+
+/* Please note that all internal length scales are in Mpc, not h^-1 Mpc! */
+
+/* -------------------------- Prototypes ----------------------------- */
+
+int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
+	int degen_hdm, double omega_lambda, double hubble, double redshift);
+double TFmdm_onek_mpc(double kk);
+double TFmdm_onek_hmpc(double kk);
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Convenience from Numerical Recipes in C, 2nd edition
+------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
+static double sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+*/
+#define SQR(a) (a*a)
+
+/* ------------------------- Global Variables ------------------------ */
+
+/* The following are set in TFmdm_set_cosm() */
+double   alpha_gamma,	/* sqrt(alpha_nu) */
+	alpha_nu,	/* The small-scale suppression */
+	beta_c,		/* The correction to the log in the small-scale */
+	num_degen_hdm,	/* Number of degenerate massive neutrino species */
+	f_baryon,	/* Baryon fraction */
+	f_bnu,		/* Baryon + Massive Neutrino fraction */
+	f_cb,		/* Baryon + CDM fraction */
+	f_cdm,		/* CDM fraction */
+	f_hdm,		/* Massive Neutrino fraction */
+	growth_k0,	/* D_1(z) -- the growth function as k->0 */
+	growth_to_z0,	/* D_1(z)/D_1(0) -- the growth relative to z=0 */
+	hhubble,	/* Need to pass Hubble constant to TFmdm_onek_hmpc() */
+	k_equality,	/* The comoving wave number of the horizon at equality*/
+	obhh,		/* Omega_baryon * hubble^2 */
+	omega_curv,	/* = 1 - omega_matter - omega_lambda */
+	omega_lambda_z, /* Omega_lambda at the given redshift */
+	omega_matter_z,	/* Omega_matter at the given redshift */
+	omhh,		/* Omega_matter * hubble^2 */
+	onhh,		/* Omega_hdm * hubble^2 */
+	p_c,		/* The correction to the exponent before drag epoch */
+	p_cb,		/* The correction to the exponent after drag epoch */
+	sound_horizon_fit,  /* The sound horizon at the drag epoch */
+	theta_cmb,	/* The temperature of the CMB, in units of 2.7 K */
+	y_drag,		/* Ratio of z_equality to z_drag */
+	z_drag,		/* Redshift of the drag epoch */
+	z_equality;	/* Redshift of matter-radiation equality */
+
+/* The following are set in TFmdm_onek_mpc() */
+double	gamma_eff,	/* Effective \Gamma */
+	growth_cb,	/* Growth factor for CDM+Baryon perturbations */
+	growth_cbnu,	/* Growth factor for CDM+Baryon+Neutrino pert. */
+	max_fs_correction,  /* Correction near maximal free streaming */
+	qq,		/* Wavenumber rescaled by \Gamma */
+	qq_eff,		/* Wavenumber rescaled by effective Gamma */
+	qq_nu,		/* Wavenumber compared to maximal free streaming */
+	tf_master,	/* Master TF */
+	tf_sup,		/* Suppressed TF */
+	y_freestream; 	/* The epoch of free-streaming for a given scale */
+
+/* Finally, TFmdm_onek_mpc() and TFmdm_onek_hmpc() give their answers as */
+double   tf_cb,		/* The transfer function for density-weighted
+			CDM + Baryon perturbations. */
+	tf_cbnu;	/* The transfer function for density-weighted
+			CDM + Baryon + Massive Neutrino perturbations. */
+
+/* By default, these functions return tf_cb */
+
+/* ------------------------- TFmdm_set_cosm() ------------------------ */
+int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
+	int degen_hdm, double omega_lambda, double hubble, double redshift)
+/* This routine takes cosmological parameters and a redshift and sets up
+all the internal scalar quantities needed to compute the transfer function. */
+/* INPUT: omega_matter -- Density of CDM, baryons, and massive neutrinos,
+				in units of the critical density. */
+/* 	  omega_baryon -- Density of baryons, in units of critical. */
+/* 	  omega_hdm    -- Density of massive neutrinos, in units of critical */
+/* 	  degen_hdm    -- (Int) Number of degenerate massive neutrino species */
+/*        omega_lambda -- Cosmological constant */
+/* 	  hubble       -- Hubble constant, in units of 100 km/s/Mpc */
+/*        redshift     -- The redshift at which to evaluate */
+/* OUTPUT: Returns 0 if all is well, 1 if a warning was issued.  Otherwise,
+	sets many global variables for use in TFmdm_onek_mpc() */
+{
+    double z_drag_b1, z_drag_b2, omega_denom;
+    int qwarn;
+    qwarn = 0;
+
+    theta_cmb = 2.728/2.7;	/* Assuming T_cmb = 2.728 K */
+
+    /* Look for strange input */
+    if (omega_baryon<0.0) {
+	fprintf(stderr,
+	  "TFmdm_set_cosm(): Negative omega_baryon set to trace amount.\n");
+	qwarn = 1;
+    }
+    if (omega_hdm<0.0) {
+	fprintf(stderr,
+	  "TFmdm_set_cosm(): Negative omega_hdm set to trace amount.\n");
+	qwarn = 1;
+    }
+    if (hubble<=0.0) {
+	fprintf(stderr,"TFmdm_set_cosm(): Negative Hubble constant illegal.\n");
+	exit(1);  /* Can't recover */
+    } else if (hubble>2.0) {
+	fprintf(stderr,"TFmdm_set_cosm(): Hubble constant should be in units of 100 km/s/Mpc.\n");
+	qwarn = 1;
+    }
+    if (redshift<=-1.0) {
+	fprintf(stderr,"TFmdm_set_cosm(): Redshift < -1 is illegal.\n");
+	exit(1);
+    } else if (redshift>99.0) {
+	fprintf(stderr,
+	  "TFmdm_set_cosm(): Large redshift entered.  TF may be inaccurate.\n");
+	qwarn = 1;
+    }
+    if (degen_hdm<1) degen_hdm=1;
+    num_degen_hdm = (double) degen_hdm;	
+	/* Have to save this for TFmdm_onek_mpc() */
+    /* This routine would crash if baryons or neutrinos were zero, 
+	so don't allow that */
+    if (omega_baryon<=0) omega_baryon=1e-5;
+    if (omega_hdm<=0) omega_hdm=1e-5;
+
+    omega_curv = 1.0-omega_matter-omega_lambda;
+    omhh = omega_matter*SQR(hubble);
+    obhh = omega_baryon*SQR(hubble);
+    onhh = omega_hdm*SQR(hubble);
+    f_baryon = omega_baryon/omega_matter;
+    f_hdm = omega_hdm/omega_matter;
+    f_cdm = 1.0-f_baryon-f_hdm;
+    f_cb = f_cdm+f_baryon;
+    f_bnu = f_baryon+f_hdm;
+
+    /* Compute the equality scale. */
+    z_equality = 25000.0*omhh/SQR(SQR(theta_cmb));	/* Actually 1+z_eq */
+    k_equality = 0.0746*omhh/SQR(theta_cmb);
+
+    /* Compute the drag epoch and sound horizon */
+    z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
+    z_drag_b2 = 0.238*pow(omhh,0.223);
+    z_drag = 1291*pow(omhh,0.251)/(1.0+0.659*pow(omhh,0.828))*
+		(1.0+z_drag_b1*pow(obhh,z_drag_b2));
+    y_drag = z_equality/(1.0+z_drag);
+
+    sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1.0+10.0*pow(obhh,0.75));
+
+    /* Set up for the free-streaming & infall growth function */
+    p_c = 0.25*(5.0-sqrt(1+24.0*f_cdm));
+    p_cb = 0.25*(5.0-sqrt(1+24.0*f_cb));
+
+    omega_denom = omega_lambda+SQR(1.0+redshift)*(omega_curv+
+			omega_matter*(1.0+redshift));
+    omega_lambda_z = omega_lambda/omega_denom;
+    omega_matter_z = omega_matter*SQR(1.0+redshift)*(1.0+redshift)/omega_denom;
+    growth_k0 = z_equality/(1.0+redshift)*2.5*omega_matter_z/
+	    (pow(omega_matter_z,4.0/7.0)-omega_lambda_z+
+	    (1.0+omega_matter_z/2.0)*(1.0+omega_lambda_z/70.0));
+    growth_to_z0 = z_equality*2.5*omega_matter/(pow(omega_matter,4.0/7.0)
+	    -omega_lambda + (1.0+omega_matter/2.0)*(1.0+omega_lambda/70.0));
+    growth_to_z0 = growth_k0/growth_to_z0;	
+    
+    /* Compute small-scale suppression */
+    alpha_nu = f_cdm/f_cb*(5.0-2.*(p_c+p_cb))/(5.-4.*p_cb)*
+	pow(1+y_drag,p_cb-p_c)*
+	(1+f_bnu*(-0.553+0.126*f_bnu*f_bnu))/
+	(1-0.193*sqrt(f_hdm*num_degen_hdm)+0.169*f_hdm*pow(num_degen_hdm,0.2))*
+	(1+(p_c-p_cb)/2*(1+1/(3.-4.*p_c)/(7.-4.*p_cb))/(1+y_drag));
+    alpha_gamma = sqrt(alpha_nu);
+    beta_c = 1/(1-0.949*f_bnu);
+    /* Done setting scalar variables */
+    hhubble = hubble;	/* Need to pass Hubble constant to TFmdm_onek_hmpc() */
+    return qwarn;
+}
+
+/* ---------------------------- TFmdm_onek_mpc() ---------------------- */
+
+double TFmdm_onek_mpc(double kk)
+/* Given a wavenumber in Mpc^-1, return the transfer function for the
+cosmology held in the global variables. */
+/* Input: kk -- Wavenumber in Mpc^-1 */
+/* Output: The following are set as global variables:
+	growth_cb -- the transfer function for density-weighted
+			CDM + Baryon perturbations. 
+ 	growth_cbnu -- the transfer function for density-weighted
+			CDM + Baryon + Massive Neutrino perturbations. */
+/* The function returns growth_cb */
+{
+    double tf_sup_L, tf_sup_C;
+    double temp1, temp2;
+
+    qq = kk/omhh*SQR(theta_cmb);
+
+    /* Compute the scale-dependent growth functions */
+    y_freestream = 17.2*f_hdm*(1+0.488*pow(f_hdm,-7.0/6.0))*
+		SQR(num_degen_hdm*qq/f_hdm);
+    temp1 = pow(growth_k0, 1.0-p_cb);
+    temp2 = pow(growth_k0/(1+y_freestream),0.7);
+    growth_cb = pow(1.0+temp2, p_cb/0.7)*temp1;
+    growth_cbnu = pow(pow(f_cb,0.7/p_cb)+temp2, p_cb/0.7)*temp1;
+
+    /* Compute the master function */
+    gamma_eff =omhh*(alpha_gamma+(1-alpha_gamma)/
+		(1+SQR(SQR(kk*sound_horizon_fit*0.43))));
+    qq_eff = qq*omhh/gamma_eff;
+
+    tf_sup_L = log(2.71828+1.84*beta_c*alpha_gamma*qq_eff);
+    tf_sup_C = 14.4+325/(1+60.5*pow(qq_eff,1.11));
+    tf_sup = tf_sup_L/(tf_sup_L+tf_sup_C*SQR(qq_eff));
+
+    qq_nu = 3.92*qq*sqrt(num_degen_hdm/f_hdm);
+    max_fs_correction = 1+1.2*pow(f_hdm,0.64)*pow(num_degen_hdm,0.3+0.6*f_hdm)/
+		(pow(qq_nu,-1.6)+pow(qq_nu,0.8));
+    tf_master = tf_sup*max_fs_correction;
+
+    /* Now compute the CDM+HDM+baryon transfer functions */
+    tf_cb = tf_master*growth_cb/growth_k0;
+    tf_cbnu = tf_master*growth_cbnu/growth_k0;
+    return tf_cb;
+}
+
+/* ---------------------------- TFmdm_onek_hmpc() ---------------------- */
+
+double TFmdm_onek_hmpc(double kk)
+/* Given a wavenumber in h Mpc^-1, return the transfer function for the
+cosmology held in the global variables. */
+/* Input: kk -- Wavenumber in h Mpc^-1 */
+/* Output: The following are set as global variables:
+	growth_cb -- the transfer function for density-weighted
+			CDM + Baryon perturbations. 
+ 	growth_cbnu -- the transfer function for density-weighted
+			CDM + Baryon + Massive Neutrino perturbations. */
+/* The function returns growth_cb */
+{
+    return TFmdm_onek_mpc(kk*hhubble);
+}
+
+
Index: trunk/Cosmo/SimLSS/hu_power.h
===================================================================
--- trunk/Cosmo/SimLSS/hu_power.h	(revision 3315)
+++ trunk/Cosmo/SimLSS/hu_power.h	(revision 3315)
@@ -0,0 +1,17 @@
+#ifndef HU_POWER_SEEN
+#define HU_POWER_SEEN
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
+	int degen_hdm, double omega_lambda, double hubble, double redshift);
+double TFmdm_onek_mpc(double kk);
+double TFmdm_onek_hmpc(double kk);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // HU_POWER_SEEN
Index: trunk/Cosmo/SimLSS/hu_sigma8.c
===================================================================
--- trunk/Cosmo/SimLSS/hu_sigma8.c	(revision 3315)
+++ trunk/Cosmo/SimLSS/hu_sigma8.c	(revision 3315)
@@ -0,0 +1,526 @@
+/* Code to calculate sigma8 for COBE-normalized CDM and MDM models, 
+using the approximation of Eisenstein & Hu (1997). */
+
+/* This code will calculate sigma_8 and other quantities for CDM and
+MDM models.  Because the fitting form for the transfer function does not
+include the baryon oscillations, the answer will not be accurate for
+large baryon fractions.  The COBE normalization is from the papers
+on Bunn & White (1997) and Hu & White (1997). */
+
+/* Accuracy is generally better than 5%.  Recall that the COBE normalization
+is only 7% accurate (1-sigma).  This program gives low answers compared to
+CMBfast for lower values of Omega_0*h^2, and vice versa for higher values. 
+Values of Omega_0 near 0.2 can produce answers that are 5% low. */
+
+/* The integrator used here is a bit crude, but suffices for rapid
+10^-3 accuracy.  You could get better performance by substituting an
+integrator from a numerical package.  See the end of the file for details. */
+
+/* This program uses the routines TFmdm_set_cosm() and TFmdm_onek_hmpc(),
+available from http://www.sns.ias.edu/~whu/power/power.html in the file
+power.c.  To compile this program, you'll need that file as well, e.g.
+	cc -O sigma8.c power.c -o sigma8 -lm 		*/
+
+/* Lists of versions and bugs:
+v1.0 -- 10/23/97 -- Original release version.
+v1.1 -- 11/13/97 -- We correct a missing factor of omega_m in the
+		    conversion between mass and length scales [lines
+		    32 and 34 of main()].  This only affects the results 
+		    if you explicitly used the -mass option.
+v1.2 -- 04/08/98 -- We correct a typo in the second-order tilt scaling
+		    of the COBE normalization for flat models with tensors.
+v1.3 -- 09/15/98 -- We correct a typo in the tilt scaling of the COBE
+		    normalization for open models with tensors.
+v1.4 -- 12/02/98 -- We correct the COBE normalization for tilted open models,
+		    undoing the change of v1.3 and fixing what was
+		    wrong there to begin with.
+*/
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
+	int degen_hdm, double omega_lambda, double hubble, double redshift);
+double TFmdm_onek_hmpc(double kk);
+
+double sig2tophat(double x);
+double sig2gauss(double x);
+double sig2sharpk(double x);
+double sig2tophatnu(double x);
+double sig2gaussnu(double x);
+double sig2sharpknu(double x);
+
+#define WEIRD  -2538	/* A strange number, unlikely to be entered */
+#define RHOCRIT 2.78e11	/* The critical density, for h=1, in Msun/Mpc^3 */
+
+/* Convenience from Numerical Recipes in C, 2nd edition
+------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
+static double sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+*/
+#define SQR(a) (a*a)
+
+/* Global variable */
+double   omega_m, omega_b, omega_n, 
+	numdegen, lambda, redshift, 
+	hubble, scale, mass, tilt;
+double   deltaH;		/* COBE-normalization */
+double   sig2_scale;	/* For passing a length scale to the integrator */
+int 	qtensors, quiet, qredshift;
+
+/* ========================== User interaction routines ============== */
+void usage()
+/* Print out information about the allowed options and then quit. */
+{
+    fprintf(stderr,"Command line options:\n");
+    fprintf(stderr,"-om or -omega <omega_matter, including CDM, baryons, massive neutrinos>\n");
+    fprintf(stderr,"-ob <omega_baryon> -fb <baryon_fraction>\n");
+    fprintf(stderr,"-on <omega_neutrino> -fn <neutrino_fraction>\n");
+    fprintf(stderr,"-numdegen or -Nnu <number of degenerate neutrino species>, default 1\n");
+    fprintf(stderr,"-lambda <omega_lambda> -flat (choose lambda appropriately)\n");
+    fprintf(stderr,"-h or -hubble <H_0, units of 100 km/s/Mpc>\n");
+    fprintf(stderr,"-H or -Hubble <H_0, units of 1 km/s/Mpc>\n");
+    fprintf(stderr,"-tilt or -n <scalar spectral tilt>, no tensors.\n");
+    fprintf(stderr,"-Tilt or -N <scalar spectral tilt>, power-law tensors.\n");
+    fprintf(stderr,"-z <redshift>, if none given then appropriate value assumed.\n");
+    fprintf(stderr,"The following request the rms fluctuations inside a tophat window function:\n");
+    fprintf(stderr,"-scale <comoving length scale, h^-1 Mpc>\n");
+    fprintf(stderr,"-mass <mass scale, h^-1 Msun>\n");
+    fprintf(stderr,"If neither -scale or -mass is given, sigma_8 and other quantities are computed.\n");
+    fprintf(stderr,"-quiet (print only answer to stdout)\n");
+    exit(1);
+    return;
+}
+
+void crash(char *s)
+{
+    fprintf(stderr,"%s\n", s); exit(1); return;
+}
+
+void warn(char *s)
+{
+    fprintf(stderr,"%s\n", s); return;
+}
+
+
+void parse_command_line(int argc, char *argv[])
+/* Parse the command line, and look for illegal entries */
+{
+    int i;
+    if (argc==1) usage();
+    i = 1;
+    omega_m = 1.0;
+    omega_b = 0.05;
+    omega_n = 0.0;
+    numdegen = 1.0;
+    lambda = 0.0;
+    redshift = 0.0;
+    hubble = 0.5;
+    scale = -1.0;
+    mass = -1.0;
+    qtensors = 0;
+    quiet = 0;
+    tilt = 1.0;
+    qredshift = 0;
+
+    while (i<argc) {
+	if (!strcmp(argv[i],"-om")||!strcmp(argv[i],"-omega")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&omega_m); ++i;
+	} else if (!strcmp(argv[i],"-ob")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&omega_b); ++i;
+	} else if (!strcmp(argv[i],"-fb")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&omega_b); ++i;
+	    omega_b *= -1.0;	/* Signal that a fraction was desired */
+	} else if (!strcmp(argv[i],"-on")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&omega_n); ++i;
+	} else if (!strcmp(argv[i],"-fn")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&omega_n); ++i;
+	    omega_n *= -1.0;	/* Signal that a fraction was desired */
+	} else if (!strcmp(argv[i],"-h")||!strcmp(argv[i],"-hubble")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&hubble); ++i;
+	} else if (!strcmp(argv[i],"-H")||!strcmp(argv[i],"-Hubble")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&hubble); ++i;
+	    hubble /= 100.0;	/* Units were 1 km/s/Mpc */
+	} else if (!strcmp(argv[i],"-lambda")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&lambda); ++i;
+	} else if (!strcmp(argv[i],"-flat")) {
+	    ++i;
+	    lambda = WEIRD;
+	} else if (!strcmp(argv[i],"-numdegen")||!strcmp(argv[i],"-Nnu")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&numdegen); ++i;
+	} else if (!strcmp(argv[i],"-n")||!strcmp(argv[i],"-tilt")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&tilt); ++i;
+	    qtensors = 0;
+	} else if (!strcmp(argv[i],"-N")||!strcmp(argv[i],"-Tilt")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&tilt); ++i;
+	    qtensors = 1;
+	} else if (!strcmp(argv[i],"-z")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&redshift); ++i;
+	    qredshift = 1;
+	} else if (!strcmp(argv[i],"-scale")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&scale); ++i;
+	} else if (!strcmp(argv[i],"-mass")) {
+	    ++i;
+	    if (i>=argc) usage();
+	    sscanf(argv[i],"%lf",&mass); ++i;
+	} else if (!strcmp(argv[i],"-quiet")) {
+	    ++i;
+	    quiet = 1;
+	} else usage();
+    }
+
+    /* Look for illegal entries */
+    if (omega_m<=0.0) crash("Omega must be positive.");
+    if (hubble<=0.0) crash("Hubble constant must be positive.");
+	else if (hubble>10.0) warn("Did you get the Hubble units right?");
+    if (numdegen<0) {
+	warn("Negative number of degenerate species reset to 1.");
+	numdegen = 1.0;
+    }
+    if (redshift<=-1.0) crash("Redshift must exceed -1.");
+    if (scale>0.0 && mass>0.0) {
+	warn("Input length OR mass scale.  Defaulting to length.");
+	mass = -1.0;
+    }
+
+    /* Now look for fractions, if they were input. */
+    if (omega_b<0.0) omega_b = omega_m*(-omega_b);	
+    if (omega_n<0.0) omega_n = omega_m*(-omega_n);	
+    if (lambda== WEIRD) lambda = 1-omega_m;	/* Set to be flat */
+
+    if (omega_b+omega_n>omega_m) crash("Baryons + Neutrinos exceeds Omega_0.");
+	else if (omega_b+omega_n>0.6*omega_m) 
+	    warn("CDM fraction is below range of TF fitting formula.");
+    if (tilt>1.3 || tilt < 0.7) 
+	warn("Tilt out of range of Bunn and White fitting formula.");
+    if (tilt!=1 && (omega_m<0.2 || omega_m>1.0)) 
+	warn("Omega out of range of Bunn and White fitting formula.");
+    if (tilt==1 && (omega_m<0.2 || omega_m>1.6 || lambda<0.0 || lambda>0.8))
+	warn("Omega or Lambda out of range of Bunn and White fitting formula.");
+    if (qtensors!=0&&(tilt>1.0||
+	(fabs(lambda)>1e-5&&fabs(1.0-omega_m-lambda)>1e-5))) 
+	crash("Tensor formulae only hold for tilts < 1 and simple cosmologies.");
+    if (omega_b/omega_m>0.4) warn("Baryon oscillations will be large.  TF fitting form will be inaccurate.");
+    if (omega_n/omega_m>0.4) warn("Neutrino fraction outside of range of TF fitting form.");
+    if (omega_m*hubble*hubble>0.4) warn("Omega_0 h^2 outside of range of TF fitting form.");
+    if (omega_m*hubble*hubble<0.06) warn("Omega_0 h^2 outside of range of TF fitting form.");
+    if (redshift>29) warn("Redshift is higher than range of TF fitting form.");
+
+    return;
+}
+
+/* ==================================== MAIN() ========================= */
+
+void main(int argc, char *argv[])
+{
+    double erfcc(double x);
+    double sigma, sigma8;
+    extern double growth_to_z0, z_equality, alpha_nu, beta_c;
+    double cobenorm();
+    double integrate(double (*func)(double));
+    parse_command_line(argc, argv);
+
+    /* Set all the background quantities */
+    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
+			hubble, redshift);
+    deltaH = cobenorm();
+
+    if (!quiet) {	/* Echo cosmological information to the screen */
+	printf("Omega_0  = %5.3f, Omega_baryon = %5.3f, Omega_neutrino = %5.3f\n",
+		omega_m, omega_b, omega_n);
+	printf("Hubble   = %5.3f, Lambda       = %5.3f, N_neutrinos    = %d\n",
+		hubble, lambda, (int) numdegen);
+	printf("Redshift = %5.3f, Growth = %6.4f relative to z=0\n",
+		redshift, growth_to_z0);
+	printf("Tilt     = %5.3f, COBE normalization: deltaH = %8.3e",
+		tilt, deltaH);
+	if (tilt!=1.0) 
+	    if (qtensors) printf(" with tensors.\n");
+	    else printf(" without tensors.\n");
+	else printf("\n");
+    }
+    /* If we've been asked to compute for a particular scale, do that */
+    if (scale>0 || mass>0) {
+	if (mass<0)
+	    mass = 4.0*M_PI/3.0*omega_m*RHOCRIT*scale*scale*scale;
+	else if (scale<0) 
+	    scale = pow(3.0/4.0/M_PI/RHOCRIT/omega_m*mass, 1.0/3.0);
+	/* Integrate tophat for the given scale */
+	sig2_scale = scale;
+	sigma = sqrt(integrate(sig2tophatnu));
+
+	/* Print out the result */
+	if (!quiet) printf("For length scale %6.4f h^-1 Mpc, mass scale %8.3e h^-1 Msun.\nrms tophat fluctuations are ",
+		sig2_scale, mass);
+	printf("%6.4f\n", sigma);
+    } else {	
+	/* No scale was requested, so we do the standard tasks */
+	/* Calculate sigma_8 */
+	sig2_scale = 8.0;	
+	sigma8 = sqrt(integrate(sig2tophatnu));
+	if (quiet) {
+	    printf("%6.4f\n", sigma8);
+	    return;	/* This is all we're doing */
+	}
+	printf("Sigma_8  = %6.4f\n", sigma8);
+
+	/* Calculate sigma_50 */
+	sig2_scale = 50.0;
+	sigma = sqrt(integrate(sig2tophatnu));
+	printf("Sigma_50 = %6.4f   Sigma_50/sigma_8 = %6.4f\n",
+		sigma, sigma/sigma8);
+
+	/* We default to z=4 for DLA */
+	if (!qredshift) {
+	    redshift = 4.0;		
+	    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
+			hubble, redshift);
+	}
+	sig2_scale = 50.0/100.0/sqrt(1+redshift)/sqrt(2*3.14159)*
+		pow(8.0*3.14159/3.0/sqrt(89.0),1./3.);
+	sigma = sqrt(integrate(sig2gaussnu));
+	printf("Omega_gas = %8.3e at redshift %5.3f due to \n   sigma(50 km/s) = sigma(%5.3f h^-1 Mpc) = %6.4f (DLA, Klypin et al. 1994)\n", 
+	    omega_b*erfcc(1.33/1.414/sigma), redshift, sig2_scale, sigma);
+
+	/* We default to z=3 for Lyman alpha forest */
+	if (!qredshift) {
+	    redshift = 3.0;		
+	    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
+			hubble, redshift);
+	}
+	sig2_scale = 1.0/(34.0*sqrt(omega_m/2.0));	/* Gnedin's k_34 */
+	sigma = sqrt(integrate(sig2gauss));
+	printf("Sigma_Lya = %6.4f at redshift %5.3f (Ly-alpha forest, Gnedin 1997)\n", 
+		sigma, redshift);
+
+	printf("Small-scale: alpha_nu = %6.4f, beta_c = %6.4f\n", 
+		alpha_nu, beta_c);
+    }
+    return;
+}
+
+/* ========================== Integration Kernals ======================= */
+
+double W2tophat(double x)
+/* The square of the tophat window */
+{
+    double j1onx;
+    if (x<0.03) 
+	j1onx = (1.0/3.0-x*x*1.0/30.0);
+    else j1onx = (cos(x)-sin(x)/x)/x/x;
+    return 9.0*j1onx*j1onx;
+}
+
+double sig2tophat(double x)
+/* This is the integrand for the mass fluctuation integral
+sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
+broken and recombined so that one integrates on (0,1] in a new variable x. */
+/* Here, the window is a real-space tophat of radius sig2_scale, the
+latter being passed from outside in units of h^-1 Mpc */
+{
+    double powerspcb(double k);
+    return (powerspcb(x/sig2_scale)*W2tophat(x)+
+ 	    powerspcb(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
+}
+
+double sig2tophatnu(double x)
+/* Same as sig2tophat(), but with P_cbnu */
+{
+    double powerspcbnu(double k);
+    return (powerspcbnu(x/sig2_scale)*W2tophat(x)+
+ 	    powerspcbnu(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
+}
+
+double sig2gauss(double x)
+/* This is the integrand for the mass fluctuation integral
+sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
+broken and recombined so that one integrates on (0,1] in a new variable x. */
+/* Here, the window is a gaussian of radius sig2_scale, the
+latter being passed from outside in units of h^-1 Mpc */
+{
+    double powerspcb(double k);
+    return (powerspcb(x/sig2_scale)*exp(-x*x)+
+ 	    powerspcb(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
+}
+
+double sig2gaussnu(double x)
+/* Same as sig2gauss(), but with P_cbnu */
+{
+    double powerspcbnu(double k);
+    return (powerspcbnu(x/sig2_scale)*exp(-x*x)+
+ 	    powerspcbnu(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
+}
+
+double sig2sharpk(double x)
+/* This is the integrand for the mass fluctuation integral
+sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
+broken and recombined so that one integrates on (0,1] in a new variable x. */
+/* Here, the window is a k-space tophat of radius 1/sig2_scale, the
+latter being passed from outside in units of h^-1 Mpc */
+{
+    double powerspcb(double k);
+    return (powerspcb(x/sig2_scale))/x;
+}
+
+double sig2sharpknu(double x)
+/* Same as sig2sharpk(), but with P_cbnu */
+{
+    double powerspcbnu(double k);
+    return (powerspcbnu(x/sig2_scale))/x;
+}
+
+/* ============================ Power Spectrum and COBE ================ */
+double powerspcbnu(double k)
+/* Returns Delta^2(k), COBE-normalized, based on the approximations of 
+Eisenstein & Hu (1997).  k is in h^-1 Mpc. */
+/* TFmdm_set_cosm() and cobenorm() must have been called before this */
+/* Returns the density-weighted CDM+Baryon+Neutrino power spectrum */
+{
+    extern double tf_cbnu, growth_to_z0;
+    TFmdm_onek_hmpc(k);  /* This sets the value in tf_cbnu */
+    return pow(2997.0*k, tilt+3.0)*SQR(deltaH*tf_cbnu*growth_to_z0);
+}
+
+double powerspcb(double k)
+/* Returns Delta^2(k), COBE-normalized, based on the approximations of 
+Eisenstein & Hu (1997).  k is in h^-1 Mpc. */
+/* TFmdm_set_cosm() and cobenorm() must have been called before this */
+/* Returns the density-weighted CDM+Baryon power spectrum */
+{
+    extern double growth_to_z0;
+    return pow(2997.0*k, tilt+3.0)*SQR(deltaH*TFmdm_onek_hmpc(k)*growth_to_z0); 
+}
+
+double cobenorm()
+/* Return the Bunn & White (1997) fit for delta_H */
+/* Given lambda, omega_m, qtensors, and tilt */
+/* Open model with tensors is from Hu & White */
+{
+    double n;
+    n = tilt-1;
+    if (fabs(omega_m+lambda-1.0)<1e-5) {	/* Flat universe */
+	if (qtensors) 
+	    return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
+		exp(n+1.97*n*n);
+	else return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
+		exp(-0.95*n-0.169*n*n);
+	/* Error -- was 0.0169 in earlier version of code */
+    } else if (fabs(lambda)<1e-5) {	/* No lambda */
+	if (qtensors)
+	    return 1.95e-5*pow(omega_m,-0.35-0.19*log(omega_m)-0.15*n)*
+		exp(+1.02*n+1.7*n*n);
+	 /* Error -- was exp(-1.02n-1.7*n*n) in earlier version (v1-1.3)*/
+	else return 1.95e-5*pow(omega_m, -0.35-0.19*log(omega_m)-0.17*n)*
+		exp(-n-0.14*n*n);
+	 /* Error -- was exp(-n-0.14*n*n) in earlier version (v1-1.2)*/
+	 /* Error -- was exp(n+0.14*n*n) in earlier version (v1.3) */
+    } else return 1e-5*(2.422-1.166*exp(omega_m)+0.800*exp(lambda)
+	+3.780*omega_m-2.267*omega_m*exp(lambda)+0.487*SQR(omega_m)+
+	0.561*lambda+3.329*lambda*exp(omega_m)-8.568*omega_m*lambda+
+	1.080*SQR(lambda));
+}
+
+/* The following routine is from Numerical Recipes in C, 2nd edition,
+by Press, Teulkolsky, Vetterling, and Flannery, (Cambridge Univ. Press) */
+double erfcc(double x)
+{
+	double t,z,ans;
+	z=fabs(x);
+	t=1.0/(1.0+0.5*z);
+	ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
+		t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
+		t*(-0.82215223+t*0.17087277)))))))));
+	return x >= 0.0 ? ans : 2.0-ans;
+}
+
+/* ======================= Integration Engine ======================= */
+
+double midpoint_int(double (*func)(double))
+/* Integrate the function func() on the range (0,1] */
+/* Use midpoint rule, expanding by factor of 3 each time */
+/* Linearly extrapolate to zero stepsize.  Not fancy */
+{
+    int iter, bins, j;
+    double width, xeval, sum[11];
+
+    sum[0] = (*func)(0.5);	/* First step is easy */
+    bins = 1; width = 1.0;
+    for (iter=1; iter<11; iter++) {	/* Now start increasing the bins */
+	bins *= 3;
+	width /= 3.0;
+	sum[iter]=0.0;
+	for (xeval=width/2.0,j=1;j<=bins;xeval+=width,j++) {
+	    if (j%3==2) continue;	/* We've already done this one */
+	    sum[iter]+=(*func)(xeval);
+	}
+	/* Add the new evaluations to the old ones. */
+	sum[iter]=sum[iter]*width+sum[iter-1]/3.0;
+	/* printf("%d iteration, sum = %f\n", iter, sum[iter]); */
+	/* Now decide if we're done */
+	if (iter<4) continue;	/* We certainly want to do some work */
+	/* If error scales as width^2, then linear extrapolation in
+	quantity width^2 is sum[iter]+(sum[iter]-sum[iter-1])/8.0 */
+	if (fabs( (sum[iter]+(sum[iter]-sum[iter-1])/8.0)/
+		  (sum[iter-1]+(sum[iter-1]-sum[iter-2])/8.0)-1.0)<1e-4) 
+	    return sum[iter]+(sum[iter]-sum[iter-1])/8.0;
+	/* if (fabs(sum[iter-1]/sum[iter]-1)<1e-4) return sum[iter];  */
+    }
+    crash("Too many steps in the integration.");
+}
+
+double integrate(double (*func)(double))
+/* Integrate the function func() on the range (0,1] */
+{
+    double midpoint_int(double (*func)(double));
+    return midpoint_int(func);
+}
+
+/* The above is a simple integrator.  You could get faster execution times
+by using a more accurate integrator.  The following is an example from 
+Numerical Recipes, but all that is needed is for the call to integrate()
+to return the value of func() integrated on the range (0,1]. */
+
+#ifdef YES_NR
+#include "qromo.c"
+#include "midpnt.c"
+#include "polint.c"
+#include "nrutil.c"
+
+float integrate(float (*func)(float))
+/* Integrate the function func() on the range (0,1] */
+{
+    float qromo(float (*func)(float), float a, float b,
+	float (*choose)(float(*)(float),float,float,int));
+    float midpnt(float (*func)(float), float a, float b, int n);
+    return qromo(func,0,1,midpnt);
+}
+#endif
+
Index: trunk/Cosmo/SimLSS/hu_tf_fit.c
===================================================================
--- trunk/Cosmo/SimLSS/hu_tf_fit.c	(revision 3315)
+++ trunk/Cosmo/SimLSS/hu_tf_fit.c	(revision 3315)
@@ -0,0 +1,358 @@
+/* The following routines implement all of the fitting formulae in 
+Eisenstein \& Hu (1997) */
+
+/* There are two sets of routines here.  The first set,
+
+	TFfit_hmpc(), TFset_parameters(), and TFfit_onek(),
+
+calculate the transfer function for an arbitrary CDM+baryon universe using
+the fitting formula in Section 3 of the paper.  The second set,
+
+	TFsound_horizon_fit(), TFk_peak(), TFnowiggles(), and TFzerobaryon(),
+
+calculate other quantities given in Section 4 of the paper. */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+void TFset_parameters(double omega0hh, double f_baryon, double Tcmb);
+double TFfit_onek(double k, double *tf_baryon, double *tf_cdm); 
+
+void TFfit_hmpc(double omega0, double f_baryon, double hubble, double Tcmb,
+	int numk, double *k, double *tf_full, double *tf_baryon, double *tf_cdm);
+
+double TFsound_horizon_fit(double omega0, double f_baryon, double hubble);
+double TFk_peak(double omega0, double f_baryon, double hubble);
+double TFnowiggles(double omega0, double f_baryon, double hubble, 
+		double Tcmb, double k_hmpc);
+double TFzerobaryon(double omega0, double hubble, double Tcmb, double k_hmpc);
+
+/* ------------------------ DRIVER ROUTINE --------------------------- */
+/* The following is an example of a driver routine you might use. */
+/* Basically, the driver routine needs to call TFset_parameters() to
+set all the scalar parameters, and then call TFfit_onek() for each 
+wavenumber k you desire. */
+
+/* While the routines use Mpc^-1 units internally, this driver has been
+written to take an array of wavenumbers in units of h Mpc^-1.  On the
+other hand, if you want to use Mpc^-1 externally, you can do this by
+altering the variables you pass to the driver:  
+	omega0 -> omega0*hubble*hubble, hubble -> 1.0 		*/
+
+/* INPUT: omega0 -- the matter density (baryons+CDM) in units of critical 
+	  f_baryon -- the ratio of baryon density to matter density 
+	  hubble -- the Hubble constant, in units of 100 km/s/Mpc
+	  Tcmb -- the CMB temperature in Kelvin. T<=0 uses the COBE value 2.728.
+	  numk -- the length of the following zero-offset array
+	  k[] -- the array of wavevectors k[0..numk-1]  */
+
+/* INPUT/OUTPUT: There are three output arrays of transfer functions. 
+All are zero-offset and, if used, must have storage [0..numk-1] declared
+in the calling program.  However, if you substitute the NULL pointer for
+one or more of the arrays, then that particular transfer function won't
+be outputted. The transfer functions are:
+
+	tf_full[] -- The full fitting formula, eq. (16), for the matter
+			transfer function. 
+	tf_baryon[] -- The baryonic piece of the full fitting formula, eq. 21.
+	tf_cdm[] -- The CDM piece of the full fitting formula, eq. 17. */
+
+/* Again, you can set these pointers to NULL in the function call if
+you don't want a particular output. */
+
+/* Various intermediate scalar quantities are stored in global variables, 
+so that you might more easily access them.  However, this also means that
+you would be better off not simply #include'ing this file in your programs,
+but rather compiling it separately, calling only the driver, and using
+extern declarations to access the intermediate quantities. */
+
+/* ----------------------------- DRIVER ------------------------------- */
+
+void TFfit_hmpc(double omega0, double f_baryon, double hubble, double Tcmb,
+	int numk, double *k, double *tf_full, double *tf_baryon, double *tf_cdm)
+/* Remember: k[0..numk-1] is in units of h Mpc^-1. */
+{
+    int j;
+    double tf_thisk, baryon_piece, cdm_piece;
+
+    TFset_parameters(omega0*hubble*hubble, f_baryon, Tcmb);
+
+    for (j=0;j<numk;j++) {
+    	tf_thisk = TFfit_onek(k[j]*hubble, &baryon_piece, &cdm_piece); 
+	if (tf_full!=NULL) tf_full[j] = tf_thisk;
+	if (tf_baryon!=NULL) tf_baryon[j] = baryon_piece;
+	if (tf_cdm!=NULL) tf_cdm[j] = cdm_piece;
+    }
+    return;
+}
+
+/* ------------------------ FITTING FORMULAE ROUTINES ----------------- */
+
+/* There are two routines here.  TFset_parameters() sets all the scalar
+parameters, while TFfit_onek() calculates the transfer function for a 
+given wavenumber k.  TFfit_onek() may be called many times after a single
+call to TFset_parameters() */
+
+/* Global variables -- We've left many of the intermediate results as 
+global variables in case you wish to access them, e.g. by declaring
+them as extern variables in your main program. */
+/* Note that all internal scales are in Mpc, without any Hubble constants! */
+
+double	omhh,		/* Omega_matter*h^2 */
+	obhh,		/* Omega_baryon*h^2 */
+	theta_cmb,	/* Tcmb in units of 2.7 K */
+	z_equality,	/* Redshift of matter-radiation equality, really 1+z */
+	k_equality,	/* Scale of equality, in Mpc^-1 */
+	z_drag,		/* Redshift of drag epoch */
+	R_drag,		/* Photon-baryon ratio at drag epoch */
+	R_equality,	/* Photon-baryon ratio at equality epoch */
+	sound_horizon,	/* Sound horizon at drag epoch, in Mpc */
+	k_silk,		/* Silk damping scale, in Mpc^-1 */
+	alpha_c,	/* CDM suppression */
+	beta_c,		/* CDM log shift */
+	alpha_b,	/* Baryon suppression */
+	beta_b,		/* Baryon envelope shift */
+	beta_node,	/* Sound horizon shift */
+	k_peak,		/* Fit to wavenumber of first peak, in Mpc^-1 */
+	sound_horizon_fit,	/* Fit to sound horizon, in Mpc */
+	alpha_gamma;	/* Gamma suppression in approximate TF */
+
+/* Convenience from Numerical Recipes in C, 2nd edition
+------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
+static double sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+static double cubearg;
+#define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg)
+static double pow4arg;
+#define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg)
+	Yes, I know the last one isn't optimal; it doesn't appear much */
+#define SQR(a) (a*a)
+#define CUBE(a) (a*a*a)
+#define POW4(a) (a*a*a*a)
+
+void TFset_parameters(double omega0hh, double f_baryon, double Tcmb)
+/* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */
+/* Input: omega0hh -- The density of CDM and baryons, in units of critical dens,
+		multiplied by the square of the Hubble constant, in units
+		of 100 km/s/Mpc */
+/* 	  f_baryon -- The fraction of baryons to CDM */
+/*        Tcmb -- The temperature of the CMB in Kelvin.  Tcmb<=0 forces use
+			of the COBE value of  2.728 K. */
+/* Output: Nothing, but set many global variables used in TFfit_onek(). 
+You can access them yourself, if you want. */
+/* Note: Units are always Mpc, never h^-1 Mpc. */
+{
+    double z_drag_b1, z_drag_b2;
+    double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y;
+
+    if (f_baryon<=0.0 || omega0hh<=0.0) {
+	fprintf(stderr, "TFset_parameters(): Illegal input.\n");
+	exit(1);
+    }
+    omhh = omega0hh;
+    obhh = omhh*f_baryon;
+    if (Tcmb<=0.0) Tcmb=2.728;	/* COBE FIRAS */
+    theta_cmb = Tcmb/2.7;
+
+    z_equality = 2.50e4*omhh/POW4(theta_cmb);  /* Really 1+z */
+    k_equality = 0.0746*omhh/SQR(theta_cmb);
+
+    z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
+    z_drag_b2 = 0.238*pow(omhh,0.223);
+    z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow(omhh,0.828))*
+		(1+z_drag_b1*pow(obhh,z_drag_b2));
+    
+    R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag));
+    R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality);
+
+    sound_horizon = 2./3./k_equality*sqrt(6./R_equality)*
+	    log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality)));
+
+    k_silk = 1.6*pow(obhh,0.52)*pow(omhh,0.73)*(1+pow(10.4*omhh,-0.95));
+
+    alpha_c_a1 = pow(46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532));
+    alpha_c_a2 = pow(12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582));
+    alpha_c = pow(alpha_c_a1,-f_baryon)*
+		pow(alpha_c_a2,-CUBE(f_baryon));
+    
+    beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708));
+    beta_c_b2 = pow(0.395*omhh, -0.0266);
+    beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1));
+
+    y = z_equality/(1+z_drag);
+    alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)));
+    alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G;
+
+    beta_node = 8.41*pow(omhh, 0.435);
+    beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1);
+
+    k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon;
+    sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75));
+
+    alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
+		SQR(f_baryon);
+    
+    return;
+}
+
+double TFfit_onek(double k, double *tf_baryon, double *tf_cdm)
+/* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
+	  *tf_baryon, *tf_cdm -- Input value not used; replaced on output if
+				the input was not NULL. */
+/* Output: Returns the value of the full transfer function fitting formula.
+		This is the form given in Section 3 of Eisenstein & Hu (1997).
+  	  *tf_baryon -- The baryonic contribution to the full fit.
+	  *tf_cdm -- The CDM contribution to the full fit. */
+/* Notes: Units are Mpc, not h^-1 Mpc. */
+{
+    double T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha;
+    double q, xx, xx_tilde, q_eff;
+    double T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full;
+    double T_0_L0, T_0_C0, T_0, gamma_eff; 
+    double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles;
+
+    k = fabs(k);	/* Just define negative k as positive */
+    if (k==0.0) {
+	if (tf_baryon!=NULL) *tf_baryon = 1.0;
+	if (tf_cdm!=NULL) *tf_cdm = 1.0;
+	return 1.0;
+    }
+
+    q = k/13.41/k_equality;
+    xx = k*sound_horizon;
+
+    T_c_ln_beta = log(2.718282+1.8*beta_c*q);
+    T_c_ln_nobeta = log(2.718282+1.8*q);
+    T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08));
+    T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08));
+
+    T_c_f = 1.0/(1.0+POW4(xx/5.4));
+    T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) +
+	    (1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q));
+    
+    s_tilde = sound_horizon*pow(1+CUBE(beta_node/xx),-1./3.);
+    xx_tilde = k*s_tilde;
+
+    T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q));
+    T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1+SQR(xx/5.2))+
+		alpha_b/(1+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4)));
+
+    f_baryon = obhh/omhh;
+    T_full = f_baryon*T_b + (1-f_baryon)*T_c;
+
+    /* Now to store these transfer functions */
+    if (tf_baryon!=NULL) *tf_baryon = T_b;
+    if (tf_cdm!=NULL) *tf_cdm = T_c;
+    return T_full;
+}
+
+/* ======================= Approximate forms =========================== */
+
+double TFsound_horizon_fit(double omega0, double f_baryon, double hubble)
+/* Input: omega0 -- CDM density, in units of critical density
+	  f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
+	  hubble -- Hubble constant, in units of 100 km/s/Mpc */
+/* Output: The approximate value of the sound horizon, in h^-1 Mpc. */
+/* Note: If you prefer to have the answer in  units of Mpc, use hubble -> 1
+and omega0 -> omega0*hubble^2. */ 
+{
+    double omhh, sound_horizon_fit_mpc;
+    omhh = omega0*hubble*hubble;
+    sound_horizon_fit_mpc = 
+	44.5*log(9.83/omhh)/sqrt(1+10.0*pow(omhh*f_baryon,0.75));
+    return sound_horizon_fit_mpc*hubble;
+}
+
+double TFk_peak(double omega0, double f_baryon, double hubble)
+/* Input: omega0 -- CDM density, in units of critical density
+	  f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
+	  hubble -- Hubble constant, in units of 100 km/s/Mpc */
+/* Output: The approximate location of the first baryonic peak, in h Mpc^-1 */
+/* Note: If you prefer to have the answer in  units of Mpc^-1, use hubble -> 1
+and omega0 -> omega0*hubble^2. */ 
+{
+    double omhh, k_peak_mpc;
+    omhh = omega0*hubble*hubble;
+    k_peak_mpc = 2.5*3.14159*(1+0.217*omhh)/TFsound_horizon_fit(omhh,f_baryon,1.0);
+    return k_peak_mpc/hubble;
+}
+
+double TFnowiggles(double omega0, double f_baryon, double hubble, 
+		double Tcmb, double k_hmpc)
+/* Input: omega0 -- CDM density, in units of critical density
+	  f_baryon -- Baryon fraction, the ratio of baryon to CDM density.
+	  hubble -- Hubble constant, in units of 100 km/s/Mpc
+	  Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of
+			COBE FIRAS value of 2.728 K
+	  k_hmpc -- Wavenumber in units of (h Mpc^-1). */
+/* Output: The value of an approximate transfer function that captures the
+non-oscillatory part of a partial baryon transfer function.  In other words,
+the baryon oscillations are left out, but the suppression of power below
+the sound horizon is included. See equations (30) and (31).  */
+/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1
+and omega0 -> omega0*hubble^2. */ 
+{
+    double k, omhh, theta_cmb, k_equality, q, xx, alpha_gamma, gamma_eff;
+    double q_eff, T_nowiggles_L0, T_nowiggles_C0;
+
+    k = k_hmpc*hubble;	/* Convert to Mpc^-1 */
+    omhh = omega0*hubble*hubble;
+    if (Tcmb<=0.0) Tcmb=2.728;	/* COBE FIRAS */
+    theta_cmb = Tcmb/2.7;
+
+    k_equality = 0.0746*omhh/SQR(theta_cmb);
+    q = k/13.41/k_equality;
+    xx = k*TFsound_horizon_fit(omhh, f_baryon, 1.0);
+
+    alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)*
+		SQR(f_baryon);
+    gamma_eff = omhh*(alpha_gamma+(1-alpha_gamma)/(1+POW4(0.43*xx)));
+    q_eff = q*omhh/gamma_eff;
+
+    T_nowiggles_L0 = log(2.0*2.718282+1.8*q_eff);
+    T_nowiggles_C0 = 14.2 + 731.0/(1+62.5*q_eff);
+    return T_nowiggles_L0/(T_nowiggles_L0+T_nowiggles_C0*SQR(q_eff));
+}
+
+/* ======================= Zero Baryon Formula =========================== */
+
+double TFzerobaryon(double omega0, double hubble, double Tcmb, double k_hmpc)
+/* Input: omega0 -- CDM density, in units of critical density
+	  hubble -- Hubble constant, in units of 100 km/s/Mpc
+	  Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of
+			COBE FIRAS value of 2.728 K
+	  k_hmpc -- Wavenumber in units of (h Mpc^-1). */
+/* Output: The value of the transfer function for a zero-baryon universe. */
+/* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1
+and omega0 -> omega0*hubble^2. */ 
+{
+    double k, omhh, theta_cmb, k_equality, q, T_0_L0, T_0_C0;
+
+    k = k_hmpc*hubble;	/* Convert to Mpc^-1 */
+    omhh = omega0*hubble*hubble;
+    if (Tcmb<=0.0) Tcmb=2.728;	/* COBE FIRAS */
+    theta_cmb = Tcmb/2.7;
+
+    k_equality = 0.0746*omhh/SQR(theta_cmb);
+    q = k/13.41/k_equality;
+
+    T_0_L0 = log(2.0*2.718282+1.8*q);
+    T_0_C0 = 14.2 + 731.0/(1+62.5*q);
+    return T_0_L0/(T_0_L0+T_0_C0*q*q);
+}
+
+
+/* Added by cmv 24/08/2007 */
+void TFprint_parameters(void)
+{
+  printf("omhh=%g obhh=%g theta_cmb=%g\n",omhh,obhh,theta_cmb);
+  printf("z_equality=%g k_equality=%g R_equality=%g\n",z_equality,k_equality,R_equality);
+  printf("z_drag=%g R_drag=%g sound_horizon=%g\n",z_drag,R_drag,sound_horizon);
+  printf("k_silk=%g\n",k_silk);
+  printf("alpha_c=%g beta_c=%g\n",alpha_c,beta_c);
+  printf("alpha_b=%g beta_b=%g\n",alpha_b,beta_b);
+  printf("beta_node=%g\n",beta_node);
+  printf("k_peak=%g sound_horizon_fit=%g\n",k_peak,sound_horizon_fit);
+  printf("alpha_gamma=%g\n",alpha_gamma);
+}
Index: trunk/Cosmo/SimLSS/hu_tf_fit.h
===================================================================
--- trunk/Cosmo/SimLSS/hu_tf_fit.h	(revision 3315)
+++ trunk/Cosmo/SimLSS/hu_tf_fit.h	(revision 3315)
@@ -0,0 +1,29 @@
+#ifndef HU_TF_FIT_SEEN
+#define HU_TF_FIT_SEEN
+
+//http://background.uchicago.edu/~whu/transfer/transferpage.html
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void TFset_parameters(double omega0hh, double f_baryon, double Tcmb);
+double TFfit_onek(double k, double *tf_baryon, double *tf_cdm); 
+
+void TFfit_hmpc(double omega0, double f_baryon, double hubble, double Tcmb,
+	int numk, double *k, double *tf_full, double *tf_baryon, double *tf_cdm);
+
+double TFsound_horizon_fit(double omega0, double f_baryon, double hubble);
+double TFk_peak(double omega0, double f_baryon, double hubble);
+double TFnowiggles(double omega0, double f_baryon, double hubble, 
+		double Tcmb, double k_hmpc);
+double TFzerobaryon(double omega0, double hubble, double Tcmb, double k_hmpc);
+
+/* Added by cmv 24/08/2007 */
+void TFprint_parameters(void);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // HU_TF_FIT_SEEN
