Index: /trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.cc	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/mdish.cc	(revision 4030)
@@ -48,5 +48,5 @@
     }
 }
-// Return a vector representing the power spectrum (for checking) 
+// Return a 2 D histrogram as the response(kx,ky)
 Histo2D Four2DResponse::GetResponse(int nx, int ny)
 {
@@ -61,4 +61,21 @@
       h2.BinCenter(i,j,xbc,ybc);
       h2(i,j) = Value(xbc,ybc);
+    }
+  return h2;	
+}
+
+// Return a 2 D histrogram as the response(u=kx/2 pi, v=ky/2 pi)
+Histo2D Four2DResponse::GetResponseUV(int nx, int ny)
+{
+  double kxmx = 1.2*dx_;
+  double kymx = 1.2*dy_;
+  if (typ_<3) kymx=kxmx; 
+  Histo2D h2(-kxmx,kxmx,nx,-kymx,kymx,ny);
+
+  double xbc,ybc;
+  for(int_4 j=0; j<h2.NBinY(); j++) 
+    for(int_4 i=0; i<h2.NBinX(); i++) {
+      h2.BinCenter(i,j,xbc,ybc);
+      h2(i,j) = Value(xbc*DeuxPI,ybc*DeuxPI);
     }
   return h2;	
Index: /trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- /trunk/Cosmo/RadioBeam/mdish.h	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/mdish.h	(revision 4030)
@@ -50,5 +50,9 @@
   inline  double operator()(double kx, double ky) 
   { return Value(kx, ky); } 
+  // Retourne l'histo de reponse en fonction de k_x,k_y
   virtual Histo2D GetResponse(int nx=255, int ny=255);  
+  // Retourne l'histo de reponse en fonction de u=k_x/2 pi, v=k_y/2 pi
+  virtual Histo2D GetResponseUV(int nx=255, int ny=255);  
+
   // Retourne le niveau moyen du bruit projete 1D en fonction (sqrt(u^2+v^2)
   HProf GetProjNoiseLevel(int nbin=128, bool fgnorm1=true);
Index: /trunk/Cosmo/RadioBeam/pknoise.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/pknoise.cc	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/pknoise.cc	(revision 4030)
@@ -50,4 +50,5 @@
        << "    -cszmpc sx,sy,sz : 3D real space box cell size in Mpc (default=3x3x1.5) \n"   
        << "    -z redshift : redshift  (default=1.0) \n"
+       << "    -pnf : Apply frequency dependent noise factor=ComovDa^2*(1+z)^2/Hz (default= NON) \n"
        << "    -scut SCutValue (default= -100.) \n"
        << "       if SCutValue<0. ==> SCut=MinNoisePower*(-SCutValue) \n"
@@ -91,4 +92,5 @@
   int box3dsz[3]={512,512,256};
   double cellsz[3]={3.,3.,3.};
+  bool fgnoisefreq=false;
   int NBINPK=256;
   int prtlev=0;
@@ -122,4 +124,7 @@
     else if (strcmp(arg[ka],"-nbin")==0) {
       NBINPK=atoi(arg[ka+1]);    ka+=2;
+    }
+    else if (strcmp(arg[ka],"-pnf")==0) {
+      fgnoisefreq=true;   ka++;
     }
     else if (strcmp(arg[ka],"-prt")==0) {
@@ -191,4 +196,5 @@
 	      << " DoL=" << DIAMETRE/lambda << " ) " << endl;
     Histo2D h2drep = arep_p->GetResponse();
+    Histo2D h2drepuv = arep_p->GetResponseUV();
     double repmax= h2drep.VMax();
     if (fgautoscut) {
@@ -210,4 +216,7 @@
     double angscale=1.;
     Vector angscales(box3dsz[2]);
+    double pnoise0=comov_dA_z*comov_dA_z*(1.+z_Redshift)*(1.+z_Redshift)/H_z;
+    Vector pnoisevec(box3dsz[2]);
+    pnoisevec=1.;
     angscales=angscale;
     if (NMAX<0)  { box3dsz[0]=256; box3dsz[1]=256;  box3dsz[2]=128; }
@@ -222,8 +231,16 @@
     cout << " pknoise[2.c]: Freq0=" << freq0 << " dFreq=" << dfreq << " freq(z=" << z_Redshift << ")=" 
 	 << conv.getFrequency() << " zmin=" << 1420.4/freq0-1. << " zmax=" << 1420.4/(freq0+box3dsz[2]*dfreq)-1. << endl;
-    cout << " pknoise[2.d]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
+    if(fgnoisefreq) {
+      for(int kz=0; kz<box3dsz[2]; kz++) {
+	double zreds=1420.4/(freq0+(double)kz*dfreq)-1.;
+	double comda=redshift2da(zreds);
+	pnoisevec(kz)=comda*comda*(1.+zreds)*(1.+zreds)/redshift2hz(zreds)/pnoise0;
+      }
+      cout << " pknoise[2.d]: PnoiseFactor:" << pnoisevec(0) << " ... " << pnoisevec(pnoisevec.Size()-1) << endl;
+    }
+    cout << " pknoise[2.e]: Four3DPk m3d(rg," << box3dsz[0]/2 << "," << box3dsz[1] << "," 
 	 << box3dsz[2] << ")" << endl;
     Four3DPk m3d(rg,box3dsz[0]/2,box3dsz[1],box3dsz[2]);
-    cout << " pknoise[2.d]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
+    cout << " pknoise[2.f]: m3d.SetCellSize(" << dkxmpc << "," << dkympc << "," << dkzmpc 
 	 << ") cell size (Mpc) : " << cellsz[0] << "x" << cellsz[1] << "x" << cellsz[2] << endl;
     m3d.SetCellSize(dkxmpc, dkympc, dkzmpc);
@@ -235,4 +252,5 @@
       pkn.SetFreqRange(freq0, dfreq);
       pkn.SetAngScaleConversion(angscales);
+      pkn.SetPNoiseFactor(pnoisevec);
       pkn.SetPrtLevel(prtlev,prtmod);
       HProf hpn = pkn.Compute(NBINPK);
@@ -255,5 +273,6 @@
     po << PPFNameTag("h1drep") << h1drep;
     po << PPFNameTag("h2drep") << h2drep;
-  
+    po << PPFNameTag("h2drepuv") << h2drepuv;
+
     rc = 0;
   }  // End of try bloc 
Index: /trunk/Cosmo/RadioBeam/sensfgnd21cm.tex
===================================================================
--- /trunk/Cosmo/RadioBeam/sensfgnd21cm.tex	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/sensfgnd21cm.tex	(revision 4030)
@@ -59,5 +59,12 @@
 %% Definition fonction de transfer 
 \newcommand{\TrF}{\mathbf{T}}
-
+%% Definition (u,v) , ...
+\def\uv{\mathrm{u,v}} 
+\def\uvu{\mathrm{u}} 
+\def\uvv{\mathrm{v}} 
+\def\dudv{\mathrm{d u d v}} 
+
+% Commande pour marquer les changements du papiers pour le referee
+\def\changemark{\bf }
 
 %%% Definition pour la section sur les param DE par C.Y
@@ -257,5 +264,6 @@
 temperature, with a total integration time of 86400 seconds (1 day) over a frequency band of 1 MHz. 
 The width of this frequency band is well adapted  to detection of \HI source with an intrinsic velocity 
-dispersion of few 100 km/s. These detection limits should be compared with the expected 21 cm brightness
+dispersion of few 100 km/s.
+These detection limits should be compared with the expected 21 cm brightness
 $S_{21}$ of compact sources which can be computed using the expression below (e.g.\cite{binney.98}) :
 \begin{equation}
@@ -264,5 +272,8 @@
 \end{equation}
  where $ M_{H_I} $ is the neutral hydrogen mass, $\dlum(z)$ is the luminosity distance and $\sigma_v$ 
-is the source velocity dispersion. 
+is the source velocity dispersion.  
+{\changemark The 1 MHz bandwidth mentioned above is only used for computing the
+galaxy detection thresholds and does not determine the total survey bandwidth or frequency resolution
+of an intensity mapping survey.}
 % {\color{red} Faut-il developper le calcul en annexe ? } 
  
@@ -270,5 +281,5 @@
 compact objects with a total \HI \, mass of $10^{10} M_\odot$ and an intrinsic velocity dispersion of 
 $200 \, \mathrm{km/s}$. The luminosity distance is computed for the standard 
-WMAP \LCDM universe. $10^9 - 10^{10} M_\odot$ of neutral gas mass 
+WMAP \LCDM universe \citep{komatsu.11}. $10^9 - 10^{10} M_\odot$ of neutral gas mass 
 is typical for large galaxies \citep{lah.09}. It is clear that detection of \HI sources at cosmological distances
 would require collecting area in the range of $10^6 \mathrm{m^2}$. 
@@ -334,6 +345,7 @@
 
 \subsection{ \HI power spectrum and BAO}
-In the absence of any foreground or background radiation, the brightness temperature 
-for a given direction and wavelength $\TTlam$ would be proportional to 
+In the absence of any foreground or background radiation 
+{\changemark and assuming high spin temperature $\kb T_{spin} \gg h \nu_{21}$},
+the brightness temperature for a given direction and wavelength $\TTlam$ would be proportional to 
 the local \HI number density $\etaHI(\vec{\Theta},z)$ through the relation:
 \begin{equation}
@@ -394,4 +406,10 @@
 \end{equation}
 where $k$ is the comoving wave vector and $ \dang(z) $ is the angular diameter distance.
+{ \changemark The matter power spectrum $P(k)$ has been measured using 
+galaxy surveys, for example by SDSS and 2dF at low redshift $z \lesssim 0.3$
+(\cite{cole.05}, \cite{tegmark.04}). The 21 cm brightness power spectra $P_{T_{21}}(k)$
+shown here are comparable to one from the galaxy surveys, once the mean 21 cm 
+temperature conversion factor $\left( \bar{T}_{21}(z)  \right)^2$ and redshift evolution 
+have been accounted for. }
 % It should be noted that the maximum transverse $k^{comov} $ sensitivity range 
 % for an instrument corresponds approximately to half of its angular resolution. 
@@ -474,5 +492,5 @@
 the common receivers axis, the visibilty would be written as the 2D Fourier transform
 of the product of the sky intensity and the receiver beam, for the angular frequency 
-\mbox{$(u,v)_{12} = 2 \pi( \frac{\Delta x}{\lambda} ,  \frac{\Delta y}{\lambda} )$}:
+\mbox{$(\uv)_{12} = ( \frac{\Delta x}{\lambda} ,  \frac{\Delta y}{\lambda} )$}:
 \begin{equation}
 \vis(\lambda) \simeq  \iint d\alpha d\beta \, \, I(\alpha, \beta) \,  L(\alpha, \beta) 
@@ -483,11 +501,11 @@
 two $(\alpha, \beta)$ angular planes. 
 
-Furthermore, we introduce the conjugate Fourier variables $(u,v)$ and the Fourier transforms
+Furthermore, we introduce the conjugate Fourier variables $(\uv)$ and the Fourier transforms
 of the sky intensity and the receiver beam:
 \begin{center}
 \begin{tabular}{ccc}
-$(\alpha, \beta)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & $(u,v)$ \\
-$I(\alpha, \beta, \lambda)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & ${\cal I}(u,v, \lambda)$ \\
-$L(\alpha, \beta, \lambda)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & ${\cal L}(u,v, \lambda)$ \\
+$(\alpha, \beta)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & $(\uv)$ \\
+$I(\alpha, \beta, \lambda)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & ${\cal I}(\uv, \lambda)$ \\
+$L(\alpha, \beta, \lambda)$ & \hspace{2mm} $\longrightarrow $ \hspace{2mm} & ${\cal L}(\uv, \lambda)$ \\
 \end{tabular} 
 \end{center}
@@ -495,44 +513,54 @@
 The visibility can then be interpreted as the weighted sum of the sky intensity, in an angular 
 wave number domain located around 
-$(u, v)_{12}=2 \pi( \frac{\Delta x}{\lambda} ,  \frac{\Delta y}{\lambda} )$. The weight function is 
+$(\uv)_{12}=2 \pi( \frac{\Delta x}{\lambda} ,  \frac{\Delta y}{\lambda} )$. The weight function is 
 given by the receiver beam Fourier transform. 
 \begin{equation}
-\vis(\lambda) \simeq  \iint d u d v \, \, {\cal I}(u,v, \lambda) \, {\cal L}(u - 2 \pi \frac{\Delta x}{\lambda} , v - 2 \pi \frac{\Delta y}{\lambda} , \lambda)
+\vis(\lambda) \simeq  \iint \dudv \, \, {\cal I}(\uv, \lambda) \, {\cal L}(\uvu - \frac{\Delta x}{\lambda} , \uvv - \frac{\Delta y}{\lambda} , \lambda)
 \end{equation} 
 
 A single receiver instrument would measure the total power integrated in a spot centered around the 
-origin in the $(u,v)$ or the angular wave mode plane. The shape of the spot depends on the receiver 
+origin in the $(\uv)$ or the angular wave mode plane. The shape of the spot depends on the receiver 
 beam pattern, but its extent would be $\sim 2 \pi D / \lambda$, where $D$ is the receiver physical 
 size. 
 
 The correlation signal from a pair of receivers would measure the integrated signal on a similar 
-spot, located around the central angular wave mode  $(u, v)_{12}$ determined by the relative 
+spot, located around the central angular wave mode  $(\uv)_{12}$ determined by the relative 
 position of the two receivers (see figure \ref{figuvplane}).
 In an interferometer with multiple receivers, the area covered by different receiver pairs in the 
-$(u,v)$ plane might overlap and some pairs might measure the same area (same base lines). 
+$(\uv)$ plane might overlap and some pairs might measure the same area (same base lines). 
 Several beams can be formed using different combination of the correlations from a set of  
 antenna pairs.  
 
-An instrument can thus be characterized by its $(u,v)$ plane coverage or response 
-${\cal R}(u,v,\lambda)$. For a single dish with a single receiver in the focal plane, 
+An instrument can thus be characterized by its $(\uv)$ plane coverage or response 
+${\cal R}(\uv,\lambda)$. For a single dish with a single receiver in the focal plane, 
 the instrument response is simply the Fourier transform of the beam. 
 For a single dish with multiple receivers, either as a Focal Plane Array (FPA) or 
 a multi-horn system, each beam (b) will have its own response 
-${\cal R}_b(u,v,\lambda)$. 
+${\cal R}_b(\uv,\lambda)$. 
 For an interferometer, we can compute a raw instrument response 
-${\cal R}_{raw}(u,v,\lambda)$ which corresponds to $(u,v)$ plane coverage by all 
+${\cal R}_{raw}(\uv,\lambda)$ which corresponds to $(u,v)$ plane coverage by all 
 receiver pairs  with uniform weighting. 
 Obviously, different weighting schemes can be used, changing 
-the effective beam shape and thus the response ${\cal R}_{w}(u,v,\lambda)$
+the effective beam shape and thus the response ${\cal R}_{w}(\uv,\lambda)$
 and the noise behaviour. If the same Fourier angular frequency mode is measured 
 by several receiver pairs, the raw instrument response might then be larger 
 that unity. This non normalized instrument response is used to compute the projected 
 noise power spectrum in the following section (\ref{instrumnoise}). 
-We can also define a  normalized instrument response, ${\cal R}_{norm}(u,v,\lambda) \lesssim 1$ as:
-\begin{equation}
-{\cal R}_{norm}(u,v,\lambda) = {\cal R}(u,v,\lambda) / \mathrm{Max_{(u,v)}} \left[ {\cal R}(u,v,\lambda) \right] 
+We can also define a  normalized instrument response, ${\cal R}_{norm}(\uv,\lambda) \lesssim 1$ as:
+\begin{equation}
+{\cal R}_{norm}(\uv,\lambda) = {\cal R}(\uv,\lambda) / \mathrm{Max_{(\uv)}} \left[ {\cal R}(\uv,\lambda) \right] 
 \end{equation} 
 This normalized  instrument response can be used to compute the effective instrument beam, 
 in particular in section \ref{recsec}.  
+
+{\changemark Detection of the reionisation at 21 cm band has been an active field 
+in the last decade and several groups 
+(\cite{rottgering.06}, \cite{bowman.07}, \cite{lonsdale.09}, \cite{parsons.09})  have built 
+instruments to detect reionisation signal around 100 MHz. 
+Several authors have studied the instrumental noise 
+and statistical uncertainties when measuring the reionisation signal power spectrum;
+the methods presented here to compute the instrument response 
+and sensitivities are similar to the ones developed in these publications 
+(\cite{morales.04}, \cite{bowman.06}, \cite{mcquinn.06}). }
 
 \begin{figure}
@@ -543,9 +571,9 @@
 }
 \vspace*{-15mm}
-\caption{Schematic view of the $(u,v)$ plane coverage by interferometric measurement.}
+\caption{Schematic view of the $(\uv)$ plane coverage by interferometric measurement.}
 \label{figuvplane}
 \end{figure}
 
-\subsection{Noise power spectrum}
+\subsection{Noise power spectrum computation}
 \label{instrumnoise}
 Let's consider a total power measurement using a receiver at wavelength $\lambda$, over a frequency 
@@ -555,5 +583,5 @@
 corresponds also to the noise for the visibility $\vis$ measured from two identical receivers, with uncorrelated 
 noise. If the receiver has an effective area $A \simeq \pi D^2/4$ or $A \simeq D_x D_y$, the measurement 
-corresponds to the integration of power over a spot in the angular frequency plane with an area $\sim A/\lambda^2$. The noise spectral density, in the angular frequencies plane (per unit area of angular frequencies  $\frac{\delta u}{ 2 \pi} \times \frac{\delta v}{2 \pi}$), corresponding to a visibility 
+corresponds to the integration of power over a spot in the angular frequency plane with an area $\sim A/\lambda^2$. The noise spectral density, in the angular frequencies plane (per unit area of angular frequencies  $\delta \uvu \times \uvv$), corresponding to a visibility 
 measurement from a pair of receivers can be written as:
 \begin{eqnarray}
@@ -565,15 +593,16 @@
 
 The sky temperature measurement can thus be characterized by the noise spectral power density in 
-the angular frequencies plane $P_{noise}^{(u,v)} \simeq \frac{\sigma_{noise}^2}{A / \lambda^2}$, in $\mathrm{Kelvin^2}$  
-per unit area of angular frequencies  $\frac{\delta u}{ 2 \pi} \times \frac{\delta v}{2 \pi}$:
+the angular frequencies plane $P_{noise}^{(\uv)} \simeq \frac{\sigma_{noise}^2}{A / \lambda^2}$, in $\mathrm{Kelvin^2}$  
+per unit area of angular frequencies  $\delta \uvu \times \delta \uvv$:
 We can characterize the sky temperature measurement with a radio instrument by the noise 
-spectral power density in the angular frequencies plane $P_{noise}(u,v)$ in units of $\mathrm{Kelvin^2}$  
-per unit area of angular frequencies  $\frac{\delta u}{ 2 \pi} \times \frac{\delta v}{2 \pi}$. 
+spectral power density in the angular frequencies plane $P_{noise}(\uv)$ in units of $\mathrm{Kelvin^2}$  
+per unit area of angular frequencies  $\delta \uvu \times \delta \uvv$. 
 For an interferometer made of identical receiver elements, several ($n$) receiver pairs 
-might have the same baseline. The noise power density in the corresponding $(u,v)$ plane area 
+might have the same baseline. The noise power density in the corresponding $(\uv)$ plane area 
 is then reduced by a factor $1/n$. More generally, we can write the instrument  noise 
 spectral power density using the instrument response defined in section \ref{instrumresp} :
 \begin{equation}
-P_{noise}(u,v) = \frac{ P_{noise}^{\mathrm{pair}} } { {\cal R}_{raw}(u,v,\lambda) } 
+P_{noise}(\uv) = \frac{ P_{noise}^{\mathrm{pair}} } { {\cal R}_{raw}(\uv,\lambda) } 
+\label{eq:pnoiseuv}
 \end{equation} 
 
@@ -581,14 +610,45 @@
 $P(k)$ is computed, angles are translated into comoving transverse distances, 
 and frequencies or wavelengths into comoving radial distance, using the following relations:
+{ \changemark
 \begin{eqnarray}
-\delta \alpha , \beta & \rightarrow & \delta \ell_\perp = (1+z) \, \dang(z) \, \delta \alpha,\beta  \\
+\alpha , \beta & \rightarrow &  \ell_\perp = l_x, l_y = (1+z) \, \dang(z) \,  \alpha,\beta  \\
+\uv & \rightarrow & k_\perp = k_x, k_y = 2 \pi \frac{ \uvu , \uvv }{  (1+z) \, \dang(z)  }   \label{eq:uvkxky} \\
 \delta \nu & \rightarrow & \delta \ell_\parallel = (1+z) \frac{c}{H(z)} \frac{\delta \nu}{\nu} 
   = (1+z) \frac{\lambda}{H(z)} \delta \nu \\
-\delta u , \delta v & \rightarrow & \delta k_\perp = \frac{ \delta u \, , \, \delta v }{  (1+z) \, \dang(z)  } \\
-\frac{1}{\delta \nu} & \rightarrow & \delta k_\parallel = \frac{H(z)}{c} \frac{1}{(1+z)} \, \frac{\nu}{\delta \nu}
+% \delta \uvu , \delta \uvv & \rightarrow & \delta k_\perp = 2 \pi \frac{ \delta \uvu \, , \, \delta \uvv }{  (1+z) \, \dang(z)  } \\
+\frac{1}{\delta \nu} & \rightarrow & \delta k_\parallel = 2 \pi \, \frac{H(z)}{c} \frac{1}{(1+z)} \, \frac{\nu}{\delta \nu}
  =  \frac{H(z)}{c} \frac{1}{(1+z)^2} \, \frac{\nu_{21}}{\delta \nu}
 \end{eqnarray}
-
-If we consider a uniform noise spectral density in the $(u,v)$ plane corresponding to the 
+}
+{ \changemark
+A brightness measurement at a  point $(\uv,\lambda)$, covering 
+the 3D spot  $(\delta \uvu, \delta \uvv, \delta \nu)$, would correspond 
+to cosmological power spectrum measurement at a transverse wave mode $(k_x,k_y)$
+defined by the equation \ref{eq:uvkxky}, measured at a redshift given by the observation frequency.
+The measurement noise spectral density given by the equation \ref{eq:pnoisepairD} can then be 
+translated into a 3D noise power spectrum, per unit of spatial frequencies 
+$ \frac{\delta k_x \times \delta k_y \times \delta k_z}{8 \pi^3} $ (units: $ \mathrm{K^2 \times Mpc^3}$) :  
+
+\begin{eqnarray}
+(\uv , \lambda) & \rightarrow & k_x(\uvu),k_y(\uvv), z(\lambda) \\
+P_{noise}(k_x,k_y, z)  & = & P_{noise}(\uv) 
+ \frac{ 8 \pi^3 \delta \uvu \times \delta \uvv }{\delta k_x \times \delta k_y \times \delta k_z} \\
+P_{noise}(k_x,k_y, z) & = &   \left( 2 \, \frac{\Tsys^2}{t_{int} \, \nu_{21} } \, \frac{\lambda^2}{D^2}  \right) 
+   \, \frac{1}{{\cal R}_{raw}} \, \dang^2(z) \frac{c}{H(z)} \, (1+z)^4  
+\label{eq:pnoisekxkz} 
+\end{eqnarray} 
+
+It is worthwhile to notice that the cosmological 3D noise power spectrum does not depend anymore on the 
+individual measurement bandwidth.
+In the following paragraph, we will first consider an ideal instrument 
+with uniform $(\uv)$ coverage
+in order to establish the general noise power spectrum behaviour for cosmological 21 cm surveys.
+The numerical method used to compute the 3D noise power spectrum is then presented in section 
+\ref{pnoisemeth}.
+}
+ 
+\subsubsection{Uniform $(\uv)$ coverage}
+
+If we consider a uniform noise spectral density in the $(\uv)$ plane corresponding to the 
 equation \ref{eq:pnoisepairD} above,  the three dimensional projected noise spectral density 
 can then be written as: 
@@ -644,10 +704,9 @@
 
 For a single dish of diameter $D$, or an interferometric instrument with maximal extent $D$, 
-observations provide information up to $u_{max},v_{max} \lesssim 2 \pi D / \lambda $. This value of 
-$u_{max},v_{max}$ would be mapped to a maximum transverse cosmological wave number
-$k^{\perp}_{max}$:
-\begin{equation}
-k^{\perp}  =  \frac{(u,v)}{(1+z) \dang}  \hspace{8mm} 
-k^{\perp}_{max}  \lesssim  \frac{2 \pi}{\dang \, (1+z)^2} \frac{D}{\lambda_{21}} 
+observations provide information up to $\uvu_{max},\uvv_{max} \lesssim D / \lambda $. This value of 
+$\uvu_{max},\uvv_{max}$ would be mapped to a maximum transverse cosmological wave number
+$k_{\perp}^{max}$:
+\begin{equation}
+k_{\perp}^{max}  \lesssim  \frac{2 \pi}{\dang \, (1+z)^2} \frac{D}{\lambda_{21}} 
 \label{kperpmax}
 \end{equation}   
@@ -657,9 +716,9 @@
 beams and a system noise temperature $\Tsys = 50 \mathrm{K}$.
 The survey is supposed to cover a quarter of sky $\Omega_{tot} = \pi \, \mathrm{srad}$, in one 
-year. The maximum comoving wave number $k_{max}$  is also shown as a function 
+year. The maximum comoving wave number $k^{max}$  is also shown as a function 
 of redshift, for an instrument with $D=100 \, \mathrm{m}$ maximum extent. In order 
 to take into account the radial component of $\vec{k}$ and the increase of 
-the instrument noise level with $k^{\perp}$, we have taken the effective $k_{ max}  $ 
-as half of the maximum transverse $k^{\perp} _{max}$ of \mbox{eq. \ref{kperpmax}}: 
+the instrument noise level with $k_{\perp}$, we have taken the effective $k_{ max} $ 
+as half of the maximum transverse $k_{\perp} ^{max}$ of \mbox{eq. \ref{kperpmax}}: 
 \begin{equation}
 k_{max} (z) = \frac{\pi}{\dang \, (1+z)^2} \frac{D=100 \mathrm{m}}{\lambda_{21}} 
@@ -678,10 +737,60 @@
 \label{pnkmaxfz}
 \end{figure}
-
+ 
+\subsubsection{3D noise power spectrum computation}  
+\label{pnoisemeth}
+{ \changemark
+We describe here the numerical method used to compute the 3D noise power spectrum, for a given instrument 
+response, as presented in section \ref{instrumnoise}. The noise power spectrum is a good indicator to compare 
+sensitivities for different instrument configurations, albeit the noise realization for a real instrument would not be 
+isotropic. 
+\begin{itemize}
+\item We start by a 3D Fourier coefficient grid, with the two first coordinates being the transverse angular wave modes,
+and the third being the frequency $(k_x,k_y,\nu)$. The grid is positioned at the mean redshift $z_0$ for which 
+we want to compute $P_{noise}(k)$. For the results at redshift \mbox{$z_0=1$} discussed in section  \ref{instrumnoise}, 
+the grid cell size and dimensions have been chosen to represent a box in the universe  
+\mbox{$\sim 1500 \times 1500 \times 750 \mathrm{Mpc^3}$}, 
+with $3\times3\times3 \mathrm{Mpc^3}$ cells.
+This correspond to an angular wedge $\sim 25^\circ \times 25^\circ \times (\Delta z \simeq 0.3)$. Given 
+the small angular extent, we have neglected the curvature of redshift shells. 
+\item For each redshift shell $z(\nu)$, we compute a Gaussian noise realization. $(k_x,k_y)$ is 
+converted to the $\uv$ coordinates using the equation \ref{eq:uvkxky}, and the 
+angular diameter distance $\dang(z)$ for \LCDM model with WMAP parameters. 
+The noise variance is taken proportional to $P_{noise}$ :
+\begin{equation}
+\sigma_{re}^2=\sigma_{im}^2 \propto \frac{1}{{\cal R}_{raw}(\uv,\lambda)} \, \dang^2(z) \frac{c}{H(z)} \, (1+z)^4 
+\end{equation}
+\item an FFT is then performed in the frequency or redshift direction to obtain the noise Fourier 
+complex coefficients $n(k_x,k_y,k_z)$ and the power spectrum is estimated as :
+\begin{equation}
+\tilde{P}_{noise}(k) = < | n(k_x,k_y,k_z) |^2 >  \hspace{2mm}  \mathrm{for} \hspace{2mm} 
+  \sqrt{k_x^2+k_y^2+k_z^2} = k 
+\end{equation}
+Noise samples corresponding to small instrument response, typically less than 1\% of the 
+maximum instrument response are ignored when calculating  $\tilde{P}_{noise}(k)$. 
+However, we require to have a significant fraction, typically 20\% to 50\% of the possible modes 
+$(k_x,k_y,k_z)$ measured for a given $k$ value. 
+\item the above steps are repeated a number of time to decrease the statistical fluctuations 
+due to the random generations. The averaged computed noise power spectrum is normalized using 
+equation \ref{eq:pnoisekxkz}  and the instrument and survey parameters ($\Tsys \ldots$).
+\end{itemize} 
+
+It should be noted that it is possible to obtain a  good approximation of noise 
+power spectrum shape, neglecting the redshift or frequency dependence of the 
+instrument response function and $\dang(z)$ for a small redshift interval around $z_0$,
+using a fixed  instrument response ${\cal R}(u,v,\lambda(z_0))$ and 
+a constant the radial distance $\dang(z_0)*(1+z_0)$. 
+\begin{equation}
+\tilde{P}_{noise}(k) = < | n(k_x,k_y,k_z) |^2 > \simeq < P_{noise}(u,v) , k_z > 
+\end{equation}
+The approximate power spectrum obtain through this simplified and much faster 
+method is shown as dashed curves on figure \ref{figpnoisea2g} for few instrument 
+configurations. 
+}
 
 \subsection{Instrument configurations and noise power spectrum}
 \label{instrumnoise}
-We have numerically computed the instrument response ${\cal R}(u,v,\lambda)$ 
-with uniform weights in the $(u,v)$ plane for several instrument configurations:
+We have numerically computed the instrument response ${\cal R}(\uv,\lambda)$ 
+with uniform weights in the $(\uv)$ plane for several instrument configurations:
 \begin{itemize}
 \item[{\bf a} :] A packed array of $n=121 \, D_{dish}=5 \, \mathrm{m}$ dishes, arranged in 
@@ -734,5 +843,5 @@
 \end{figure}
 
-We have used simple triangular shaped dish response in the $(u,v)$ plane. 
+We have used simple triangular shaped dish response in the $(\uv)$ plane. 
 However, we have introduced a filling factor or illumination efficiency 
 $\eta$, relating the effective dish diameter $D_{ill}$ to the 
@@ -740,5 +849,5 @@
 as $\eta^2$ or $\eta_x \eta_y$. 
 \begin{eqnarray}
-{\cal L}_\circ (u,v,\lambda) & = & \bigwedge_{[\pm 2 \pi D^{ill}/ \lambda]}(\sqrt{u^2+v^2})  \\
+{\cal L}_\circ (\uv,\lambda) & = & \bigwedge_{[\pm  D^{ill}/ \lambda]}(\sqrt{u^2+v^2})  \\
  L_\circ (\alpha,\beta,\lambda) & = & \left[ \frac{ \sin (\pi (D^{ill}/\lambda) \sin \theta ) }{\pi (D^{ill}/\lambda) \sin \theta} \right]^2 
 \hspace{4mm} \theta=\sqrt{\alpha^2+\beta^2}
@@ -751,5 +860,13 @@
 rectangular shaped antenna. The illumination efficiency factor has been taken 
 equal to $\eta_x = 0.9$ in the direction of the cylinder width, and $\eta_y = 0.8$ 
-along the cylinder length. It should be noted that the small angle approximation 
+along the cylinder length. {\changemark  We have used double triangular shaped
+response function in the $(\uv)$ for each of the receiver elements along the cylinder:
+\begin{equation}
+ {\cal L}_\Box(\uv,\lambda)  = 
+\bigwedge_{[\pm \eta_x D_x / \lambda]} (\uvu ) \times
+\bigwedge_{[\pm \eta_y D_y / \lambda ]} (\uvv ) 
+\end{equation}
+}
+It should be noted that the small angle approximation 
 used here for the expression of visibilities is not valid for the receivers along 
 the cylinder axis. However, some preliminary numerical checks indicate that 
@@ -758,12 +875,13 @@
 is the case for a transit type telescope.
 
-\begin{equation}
- {\cal L}_\Box(u,v,\lambda)  = 
-\bigwedge_{[\pm 2 \pi D^{ill}_x / \lambda]} (u ) \times
-\bigwedge_{[\pm 2 \pi D^{ill}_y / \lambda ]} (v ) 
-\end{equation}
 Figure \ref{figuvcovabcd} shows the instrument response ${\cal R}(u,v,\lambda)$ 
 for the four configurations (a,b,c,d) with $\sim 100$ receivers per 
-polarisation. The resulting projected noise spectral power density is shown in figure 
+polarisation. 
+
+{\changemark Using the numerical method sketched in section \ref{pnoisemeth}, we have 
+computed the 3D noise power spectrum for each of the eight instrument configurations discussed
+here, with a system noise temperature $\Tsys = 50 \mathrm{K}$, for a one year survey 
+of a quarter of sky $\Omega_{tot} = \pi \, \mathrm{srad}$ at a mean redshift $z_0=1, \nu_0=710 \mathrm{MHz}$.}
+The resulting noise spectral power densities are shown in figure 
 \ref{figpnoisea2g}. The increase of $P_{noise}(k)$ at low $k^{comov} \lesssim 0.02$ 
 is due to the fact that we have ignored all auto-correlation measurements.  
@@ -774,24 +892,24 @@
 \centering
 \mbox{
-\hspace*{-10mm}
-\includegraphics[width=0.90\textwidth]{Figs/uvcovabcd.pdf}
+% \hspace*{-10mm}
+\includegraphics[width=\textwidth]{Figs/uvcovabcd.pdf}
 }
-\caption{(u,v) plane coverage (raw instrument response ${\cal R}(u,v,\lambda)$ 
+\caption{$(\uv)$ plane coverage (raw instrument response ${\cal R}(\uv,\lambda)$ 
 for four configurations.
 (a) 121 $D_{dish}=5$ meter diameter dishes arranged in a compact, square array 
 of $11 \times 11$, (b) 128 dishes arranged in 8 row of 16 dishes each (fig. \ref{figconfbc}), 
 (c) 129 dishes arranged as shown in figure \ref{figconfbc} , (d) single D=75 meter diameter, with 100 beams. 
-(color scale : black $<1$, blue, green, yellow, red $\gtrsim 80$) }
+The common color scale (1 \ldots 80) is shown on the right. }
 \label{figuvcovabcd}
 \end{figure*}
  
 \begin{figure*}
-\vspace*{-25mm}
+\vspace*{-10mm}
 \centering
 \mbox{
-\hspace*{-20mm}
-\includegraphics[width=1.15\textwidth]{Figs/pkna2h.pdf}
+% \hspace*{-5mm}
+\includegraphics[width=\textwidth]{Figs/pkna2h.pdf}
 }
-\vspace*{-40mm}
+\vspace*{-20mm}
 \caption{P(k) LSS power  and noise power spectrum for several interferometer 
 configurations ((a),(b),(c),(d),(e),(f),(g)) with 121, 128, 129, 400 and 960 receivers.}
@@ -810,9 +928,13 @@
 it has been suggested that the smooth frequency dependence of the synchrotron
 emissions can be used to separate the faint LSS signal from the Galactic and radio source 
-emissions.
+emissions. {\changemark Discussion of contribution of different sources 
+to foregrounds for measurement of reionization through redshifted 21 cm,
+as well foreground subtraction using their smooth frequency dependence can 
+be found in (\cite{shaver.99}, \cite{matteo.02},\cite{oh.03}) }
 However, any real radio instrument has a beam shape which changes with 
 frequency: this instrumental effect significantly increases the difficulty and complexity of this component separation 
 technique. The effect of frequency dependent beam shape is some time referred to as {\em 
-mode mixing}. See for example \citep{morales.06}, \citep{bowman.07}. 
+mode mixing}.  {\changemark Effect of frequency dependent beam shape for foreground subtraction and 
+its application to MWA has been discussed in \citep{morales.06}  \citep{bowman.09}.} 
 
 In this section, we present a short description of the foreground emissions and 
@@ -879,5 +1001,7 @@
 varying spectral index. We assign a power law index $\beta = -2.8  \pm 0.15$ to each sky direction.
 $\beta$ has a gaussian distribution centered at -2.8 and with standard 
-deviation $\sigma_\beta = 0.15$.
+deviation $\sigma_\beta = 0.15$. {\changemark The 
+diffuse radio background spectral index has been measured  for example by 
+\citep{platania.98} or \cite{rogers.08} }
 The synchrotron contribution to the sky temperature for each cell is then 
 obtained  through the formula:
@@ -906,6 +1030,6 @@
 of the power spectrum $P(k)$ at $z=0$ computed according to the parametrization of 
 \citep{eisenhu.98}. We have used the standard cosmological parameters, 
- $H_0=71 \, \mathrm{km/s/Mpc}$, $\Omega_m=0.27$, $\Omega_b=0.044$,
-$\Omega_\lambda=0.73$ and $w=-1$.
+ $H_0=71 \, \mathrm{km/s/Mpc}$, $\Omega_m=0.264$, $\Omega_b=0.045$,
+$\Omega_\lambda=0.73$ and $w=-1$ \citep{komatsu.11}.
 An inverse FFT was then performed to compute the matter density fluctuations $\delta \rho / \rho$ 
 in the linear regime, 
@@ -966,4 +1090,20 @@
 to the factor $\sim 10^3$ of the sky brightness temperature fluctuations ($\sim$ K),
 compared to the mK LSS signal.  
+
+{ \changemark Contrary to most similar studies, where it is assumed that bright sources 
+can be nearly perfectly subtracted, our aim was to compute also their 
+effect in the foreground subtraction process. 
+The GSM model lacks the angular resolution needed to compute 
+correctly the effect of bright compact sources for 21 cm LSS observations and 
+the mode mixing due to frequency dependent instrument, while Model-II 
+provides a reasonable description of these compact sources. Our simulated 
+sky cubes have an angular resolution $3'\times3'$, well below the typical 
+$15'$ resolution of the instrument configuration considered here. 
+However, Model-II might lack spatial structures at large scales, above a degree, 
+compared to Model-I/GSM, and the frequency variations as a simple power law
+might not be realistic enough. The differences for the two sky models can be seen
+in their power spectra shown in figure  \ref{pkgsmlss}. We hope that by using 
+these two models, we have explored some of the systematic uncertainties 
+related to foreground subtraction.}
 
 It should also be noted that in section 3, we presented the different instrument 
@@ -1057,5 +1197,5 @@
 The virtual target instrument has a beam width larger than the worst real instrument beam,
 i.e at the lowest observed frequency.  
- \item For each sky direction $(\alpha, \delta)$, a power law $T = T_0 \left( \frac{\nu}{\nu_0} \right)^b$
+\item For each sky direction $(\alpha, \delta)$, a power law $T = T_0 \left( \frac{\nu}{\nu_0} \right)^b$
  is fitted to the beam-corrected brightness temperature. The fit is done through a linear $\chi^2$ fit in 
 the $\lgd ( T ) , \lgd (\nu)$ plane and we show here the results for a pure power law (P1) 
@@ -1066,5 +1206,14 @@
 \end{eqnarray*}
 where $b$ is the power law index and  $T_0 = 10^a$ is the brightness temperature at the 
-reference frequency $\nu_0$:
+reference frequency $\nu_0$.
+
+{\changemark Interferometers have poor response at small $(\uv)$ corresponding to baselines 
+smaller than interferometer element size. The $(0,0)$ mode, corresponding the mean temperature 
+can not be measured with an interferometer. We have used a simple trick to make the power law
+fitting procedure to work: we have set the mean value of the temperature for 
+each frequency plane to a power law with an index close to the synchrotron index
+and we have checked that results are not too sensitive to the arbitrarily fixed mean temperature 
+power law parameters. } 
+
 \item The difference between the beam-corrected sky temperature and the fitted power law
 $(T_0(\alpha, \delta), b(\alpha, \delta))$ is our extracted 21 cm LSS signal.
@@ -1608,6 +1757,8 @@
 \bibitem[Abdalla \& Rawlings (2005)]{abdalla.05} Abdalla, F.B. \& Rawlings, S.  2005,  \mnras, 360,  27     
 
+% reference DETF - DE eq.state parameter figure of merit 
 \bibitem[Albrecht et al. (2006)]{DETF}  Albrecht, A., Bernstein, G., Cahn, R. {\it et al.} (Dark Energy Task Force) 2006, arXiv:astro-ph/0609591 
 
+% Papier sensibilite/reconstruction CRT (cylindres) ansari et al 2008 
 \bibitem[Ansari et al. (2008)]{ansari.08} Ansari R., J.-M. Le Goff, C. Magneville, M. Moniez, N. Palanque-Delabrouille, J. Rich, 
     V. Ruhlmann-Kleider, \& C. Y\`eche , 2008 , arXiv:0807.3614
@@ -1625,9 +1776,11 @@
 % Galactic astronomy, emission HI d'une galaxie
 \bibitem[Binney \& Merrifield (1998)]{binney.98} Binney J. \& Merrifield M. , 1998 {\it Galactic Astronomy} Princeton University Press 
+% 21cm reionisation P(k) estimation and sensitivities 
+\bibitem[Bowman et al. (2006)]{bowman.06} Bowman, J.D.,  Morales, M.F., Hewitt, J.N. 2006, \apj, 638, 20-26
 % MWA description 
 \bibitem[Bowman et al. (2007)]{bowman.07} Bowman, J. D., Barnes, D.G., Briggs, F.H. et al 2007, \aj, 133, 1505-1518  
 
 %% Soustraction avant plans ds MWA 
-\bibitem[Bowman et al. (2009)]{bowman.07} Bowman, J. D., Morales, M., Hewitt, J.N., 2009, \apj, 695, 183-199  
+\bibitem[Bowman et al. (2009)]{bowman.09} Bowman, J. D., Morales, M., Hewitt, J.N., 2009, \apj, 695, 183-199  
 
 %  Intensity mapping/HSHS 
@@ -1641,6 +1794,9 @@
 Taylor, G. B., \& Broderick, J. J. 1998, AJ, 115, 1693
 
-%  Parametrisation P(k) 
-\bibitem[Eisenstein \& Hu  (1998)]{eisenhu.98}  Eisenstein D. \& Hu W. 1998, ApJ 496:605-614 (astro-ph/9709112)
+% Effet des radio-sources sur le signal 21 cm reionisation 
+\bibitem[Di Matteo  et al. (2002)]{matteo.02}  Di Matteo, T., Perna R., Abel T., Rees M.J. 2002, \apj, 564, 576-580
+
+%  Parametrisation P(k)   - (astro-ph/9709112) 
+\bibitem[Eisenstein \& Hu  (1998)]{eisenhu.98}  Eisenstein D. \& Hu W. 1998, \apj 496, 605-614 
 
 % SDSS first BAO observation 
@@ -1661,6 +1817,10 @@
 \bibitem[Jackson (2004)]{jackson.04} Jackson, C.A. 2004, \na, 48, 1187  
 
+% WMAP 7 years cosmological parameters 
+\bibitem[Komatsu et al. (2011)]{komatsu.11}   E. Komatsu, K. M. Smith, J. Dunkley {\it et al.}  2011, \apjs, 192, p. 18 \\
+\mbox{\tt http://lambda.gsfc.nasa.gov/product/map/current/params/lcdm\_sz\_lens\_wmap7.cfm} 
+
 % HI mass in galaxies 
-\bibitem[Lah et al. (2009)]{lah.09}  Philip Lah, Michael B. Pracy, Jayaram N. Chengalur et al.  2009,  \mnras, 399, 1447 
+\bibitem[Lah et al. (2009)]{lah.09}  Philip Lah, Michael B. Pracy, Jayaram N. Chengalur {\it et al.}  2009,  \mnras, 399, 1447 
 % ( astro-ph/0907.1416) 
 
@@ -1672,4 +1832,6 @@
 \bibitem[Larson et al. (2011)]{larson.11}  Larson, D.,  {\it et al.}  (WMAP) 2011, \apjs, 192, 16 
 
+%% Description MWA 
+\bibitem[Lonsdale et al. (2009)]{lonsdale.09}  Lonsdale C.J., Cappallo R.J., Morales M.F.  {\it et al.}  2009, arXiv:0903.1828
 % LSST Science book
 \bibitem[LSST.Science]{lsst.science} 
@@ -1688,6 +1850,16 @@
 \bibitem[Mauskopf et al. (2000)]{mauskopf.00} Mauskopf, P. D., Ade, P. A. R., de Bernardis, P. {\it et al.}  2000, \apjl, 536,59
 
+%% PNoise and cosmological parameters with reionization
+\bibitem[McQuinn et al. (2006)]{mcquinn.06}  McQuinn M., Zahn O., Zaldarriaga M., Hernquist L.,  Furlanetto S.R.
+2006, \apj 653, 815-834 
+
+%  Papier sur la mesure de sensibilite P(k)_reionisation 
+\bibitem[Morales \& Hewitt (2004)]{morales.04}  Morales M. \& Hewitt J., 2004, \apj, 615, 7-18  
+
 %  Papier sur le traitement des observations radio / mode mixing 
 \bibitem[Morales et al. (2006)]{morales.06} Morales, M., Bowman, J.D., Hewitt, J.N., 2006, \apj, 648, 767-773
+
+%% Foreground removal using smooth frequency dependence 
+\bibitem[Oh \& Mack (2003)]{oh.03} Oh S.P. \& Mack K.J., 2003, \mnras, 346, 871-877
 
 %  Global Sky Model Paper 
@@ -1695,4 +1867,10 @@
 \mnras, 388, 247-260 
 
+%%  Description+ resultats PAPER 
+\bibitem[Parsons et al. (2009)]{parsons.09}  Parsons A.R.,Backer D.C.,Bradley  R.F. {\it et al.}  2009, arXiv:0904.2334
+
+% Synchrotron index =-2.8 in the freq range 1.4-7.5 GHz 
+\bibitem[Platania et al. (1998)]{platania.98}   Platania P., Bensadoun M., Bersanelli M. {\it al.} 1998,  \apj 505, 473-483
+
 % Original CRT HSHS paper (Moriond Cosmo 2006 Proceedings) 
 \bibitem[Peterson et al. (2006)]{peterson.06} Peterson, J.B.,  Bandura, K., \& Pen, U.-L. 2006, arXiv:0606104  
@@ -1703,4 +1881,7 @@
 % SDSS BAO 2010  - arXiv:0907.1660 
 \bibitem[Percival et al. (2010)]{percival.10}   Percival, W.J., Reid, B.A., Eisenstein, D.J. {\it et al.},  2010, \mnras, 401, 2148-2168   
+
+% Radio spectral index between 100-200 MHz 
+\bibitem[Rogers \& Bowman (2008)]{rogers.08} Rogers, A.E.E. \& Bowman, J. D. 2008, \aj 136, 641-648
 
 %% LOFAR description 
@@ -1710,4 +1891,7 @@
 %% SDSS-3
 \bibitem[SDSS-III(2008)]{sdss3} SDSS-III 2008, http://www.sdss3.org/collaboration/description.pdf 
+
+% Reionisation: Can the reionization epoch be detected as a global signature in the cosmic background?
+\bibitem[Shaver  et al. (1999))]{shaver.99} Shaver P.A., Windhorst R. A., Madau P., de Bruyn A.G. \aap, 345, 380-390
 
 %  Frank H. Briggs, Matthew Colless, Roberto De Propris, Shaun Ferris, Brian P. Schmidt, Bradley E. Tucker
@@ -1720,4 +1904,7 @@
 % Papier 21cm-BAO Fermilab  ( arXiv:0910.5007)
 \bibitem[Seo et al (2010)]{seo.10} Seo, H.J. Dodelson, S., Marriner, J. et al,  2010, \apj, 721, 164-173 
+
+% Mesure P(k) par SDSS 
+\bibitem[Tegmark et al.  (2004)]{tegmark.04}  Tegmark M., Blanton  M.R, Strauss M.A. et al. 2004, \apj, 606, 702-740
 
 % FFT telescope 
Index: /trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/specpk.cc	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/specpk.cc	(revision 4030)
@@ -209,11 +209,11 @@
 
 // Generate mass field Fourier Coefficient
-void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales)
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales, Vector& noisevec)
 // angscale is a multiplicative factor converting transverse k (wave number) values to angular wave numbers 
 // typically = ComovRadialDistance 
 {
   uint_8 nmodeok=0;
-  if (angscales.Size() != fourAmp.SizeZ()) 
-    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales.Size()!=fourAmp.SizeZ()");
+  if ((angscales.Size() != fourAmp.SizeZ())||(noisevec.Size() != fourAmp.SizeZ()))
+    throw SzMismatchError("ComputeNoiseFourierAmp(): angscales/noisevec.Size()!=fourAmp.SizeZ()");
   H21Conversions conv;
   // fourAmp represent 3-D fourier transform of a real input array. 
@@ -225,6 +225,8 @@
     resp.setLambda(conv.getLambda());
     double angsc=angscales(kz);
+    double noisepow=noisevec(kz);
     if (prtlev_>2)  
-      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc << endl;
+      cout << " Four3DPk::ComputeNoiseFourierAmp(...) - freq=" << f0+kz*df << " -> AngSc=" << angsc 
+	   << " NoisePow=" << noisepow << endl;
     for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
       kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
@@ -234,5 +236,5 @@
 	if (rep<1.e-19)  rep=1.e-19;
 	else  nmodeok++;
-	amp = 1./sqrt(rep)/sqrt(2.);
+	amp = sqrt(noisepow/rep/2.);
 	fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
       }
@@ -445,8 +447,9 @@
 PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen, 
 				     const char* tit)
-  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ())
+  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit), angscales_(pk3.SizeZ()), pnoisefac_(pk3.SizeZ())
 {
   SetFreqRange();
   SetAngScaleConversion();
+  SetPNoiseFactor();
   SetPrtLevel();
 }
@@ -461,5 +464,5 @@
        << " ... " << angscales_(angscales_.Size()-1) << endl;
   for(int igen=0; igen<NGEN; igen++) {
-    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_);
+    pkn3d.ComputeNoiseFourierAmp(frep, freq0_, dfreq_, angscales_, pnoisefac_);
     if (igen==0) hnd = pkn3d.ComputePk(S2CUT,nbin,kmin,kmax,true);
     else pkn3d.ComputePkCumul();
Index: /trunk/Cosmo/RadioBeam/specpk.h
===================================================================
--- /trunk/Cosmo/RadioBeam/specpk.h	(revision 4029)
+++ /trunk/Cosmo/RadioBeam/specpk.h	(revision 4030)
@@ -56,5 +56,5 @@
 // typically = ComovRadialDistance 
   void ComputeNoiseFourierAmp(Four2DResponse& resp, double angscale=1., bool crmask=false);
-  void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales);
+  void ComputeNoiseFourierAmp(Four2DResponse& resp, double f0, double df, Vector& angscales, Vector& noisp);
 
 // Return the array size 
@@ -112,4 +112,8 @@
   inline void SetAngScaleConversion(Vector& angscs)
   { angscales_=angscs; }
+  inline void SetPNoiseFactor(double pnoisef=1.)
+  { pnoisefac_=pnoisef; }
+  inline void SetPNoiseFactor(Vector& pnoisefac)
+  { pnoisefac_=pnoisefac; }
   inline void SetS2Cut(double s2cut=100.)
   {  S2CUT=s2cut; }
@@ -124,4 +128,5 @@
   double freq0_,dfreq_;
   Vector angscales_; 
+  Vector pnoisefac_; 
   double S2CUT;
   int NGEN;
