Index: trunk/Cosmo/RadioBeam/README
===================================================================
--- trunk/Cosmo/RadioBeam/README	(revision 3622)
+++ trunk/Cosmo/RadioBeam/README	(revision 3756)
@@ -3,28 +3,18 @@
   Simulation programs for radio lobe reconstruction, noise estimate etc ... 
  
-  2006-2008  R. Ansari (UPS/LAL-IN2P3)    C. Magneville (DAPNIA/SPP)
+  2006-2010  R. Ansari (UPS/LAL-IN2P3)    C. Magneville (DAPNIA/SPP)
               (C) IN2P3/CNRS              (C)  DAPNIA/CEA
  ===============================================================================
 
-
-A/ 3D reconstruction (alpha,delta,freq) simulation program 
-===========================================================
-
 Note :  Uses SOPHYA class library  (http://www.sophya.org )
 
-List of files and classes :
----------------------------
-
-- treccyl.cc : main program (for tests)
-- multicyl.h .cc : Class MultiCylindres , 3D (freq,AngX,AngY) signal simulation 
-  and reconstruction. Uses class MultiBeamCyl
-- mbeamcyl.h .cc : Class MultiBeamCyl 2D (freq, AngX) signal simulation and reconstruction
-- brsrc.h .cc : Class BRSourceGen for managing list of sources 
-  (frequencies and angular position)
-- makefile 
-
-- recbeam.pic : script with example commands for display with piapp
-
-
+A/ (u-v) plane coverage and Noise-P(k) computation for different interferometer configuration
+============================================================================================== 
+List of files :
+  - specpk.h specpk.cc : 3D power spectrum computation classes  
+  - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations 
+  - pknoise.cc : main program , makefile 
+  - plpkn.pic : spiapp script for plotting pknoise.cc program output 
+   
 B/ P(k) in temperature plot + "cosmic variance" + noise ...
 =============================================================
@@ -40,2 +30,17 @@
 C.2/ pa.pic + pa.cc + pa.h : Lobe synthesis (dipoles sum)
 
+D/ 3D reconstruction (alpha,delta,freq) simulation program 
+===========================================================
+
+List of files and classes :
+---------------------------
+
+- treccyl.cc : main program (for tests)
+- multicyl.h .cc : Class MultiCylindres , 3D (freq,AngX,AngY) signal simulation 
+  and reconstruction. Uses class MultiBeamCyl
+- mbeamcyl.h .cc : Class MultiBeamCyl 2D (freq, AngX) signal simulation and reconstruction
+- brsrc.h .cc : Class BRSourceGen for managing list of sources 
+  (frequencies and angular position)
+- makerecyl 
+
+- recbeam.pic : script with example commands for display with piapp
Index: trunk/Cosmo/RadioBeam/anapkn.pic
===================================================================
--- trunk/Cosmo/RadioBeam/anapkn.pic	(revision 3622)
+++ trunk/Cosmo/RadioBeam/anapkn.pic	(revision 3756)
@@ -54,12 +54,14 @@
 
 ##  ---- VOL1,VOL2,VOL3 
-##  VOL1: 10deg x 10deg x Delta_z=0.2 : 335 x 335 x 680 Mpc^3 = 76 10^6 Mpc^3
-##  VOL2: 20deg x 20deg x Delta_z=0.2 : VOL2 = 300 10^6 Mpc^3
-##  VOL3: 60deg x 60deg x Delta_z=0.2 : VOL4 = 2700 10^6 Mpc^3
+##  VOL0: 10deg x 10deg x Delta_z=0.2 : 335 x 335 x 680 Mpc^3 = 76 10^6 Mpc^3
+##  VOL1: 20deg x 20deg x Delta_z=0.2 : VOL1 = 300 10^6 Mpc^3
+##  VOL2: 360deg x 10deg x Delta_z=0.2 : VOL2 = 2700 10^6 Mpc^3
+##  VOL3: 360deg x 30deg x Delta_z=0.2 : VOL3 = 10800 10^6 Mpc^3
 ##  VOL4: 40deg x 40deg x Delta_z=0.2 : VOL3 = 1200 10^6 Mpc^3
-set angcov ( '10x10 deg^2' '20x20 deg^2' '60x60 deg^2' '40x40 deg^2' )
-VOL1 = 76e6
-VOL2 = 300e6
-VOL3 = 2700e6
+set angcov ( '20x20 deg^2' '360x10 deg^2' '360x30 deg^2' '40x40 deg^2' )
+VOL0 = 76e6
+VOL1 = 300e6
+VOL2 = 2700e6
+VOL3 = 10800e6
 VOL4 = 1200e6
 
@@ -104,4 +106,12 @@
 endscript
 
+#  Setup 64 x D=5m sur 100mx100m , 10000 m^2
+defscript setup5
+  vpixz05 = 100
+  vpix = 100
+  kappa = 100*100/1200
+  kappa = $kappa*$kappa
+  echo "setup-5: VPix=$vpix  Kappa=$kappa"
+endscript
 
 
@@ -112,15 +122,20 @@
 
 #  --- Temps d'integration  a z=0.5 avec FOV=4x4 degres^2  (D=5m) 
-#  120 jours pour VOl1 en ~ 7 passes (7 couvertures elementaires)
-#  120 jours pour VOl2 en ~ 4x7 passes (28 couvertures elementaires)
-#  120 jours pour VOl3 en ~ 16x7 passes (225 couvertures elementaires)
-set ndays ( '120 days' '120 days'  '120 days' '120 days' ) 
-TINT1 = 1.5e6
-TINT2 = 3.7e5
-TINT3 = 4.5e4
+#  120 jours pour VOl0 en ~ 7 passes (7 couvertures elementaires)
+#  360 jours pour VOl2 en ~ 16x7 passes (225 couvertures elementaires)
+#  360 jours pour VOl1 en ~ 4x7 passes (28 couvertures elementaires)
+#  360 jours pour VOl2 en ~ 90x3 passes (270 couvertures elementaires)
+#  360 jours pour VOl3 en ~ 90x9 passes (810 couvertures elementaires)
+
+set ndays ( '360 days' '360 days'  '360 days' '360 days' ) 
+TINT1 = 1.1e6
+TINT2 = 1.1e5
+TINT3 = 3.8e4
 TINT4 = 1.e5
 
 #  Bruit Tsys = 100 K = 10^5 mK 
 TSYS = 1e5
+#  Bruit Tsys = 50 K = 5.10^4 mK 
+TSYS = 5e4
 
 DELNU = 1.e6
@@ -151,5 +166,5 @@
 set vols ( $VOL1 $VOL2 $VOL3 )
 set tints ( $TINT1 $TINT2 $TINT3 )
-set cols ( green blue red )
+set cols ( red orange blue )
 
 #############################################################
@@ -239,5 +254,5 @@
   xyl = "xylimits=0.01,0.5,1.e-3,0.4 logx logy minorticks"
   set vols ( $VOL1 $VOL2 $VOL3 )
-  set cols ( green blue red )
+  set cols ( red orange blue )
   csigpk $vols[0]
   n/plot hpkz.$CSPK/$kk%$kk ! ! "notit nsta connectpoints $cols[0] $xyl line=solid,2"
Index: trunk/Cosmo/RadioBeam/makefile
===================================================================
--- trunk/Cosmo/RadioBeam/makefile	(revision 3622)
+++ trunk/Cosmo/RadioBeam/makefile	(revision 3756)
@@ -1,39 +1,26 @@
+############################################################ 
+##  makefile for interferometer response 
+##  R.Ansari - Avril 2010  
+############################################################ 
+
 include $(SOPHYABASE)/include/sophyamake.inc
 
-all : treccyl
+all : pknoise  
 
-clean :
-	rm Objs/*
-depend :
-	mkdir Objs
+clean : 
+	rm Objs/* 
 
-treccyl : Objs/treccyl
-	echo 'makefile : treccyl made'
+pknoise : Objs/pknoise
+	echo 'makefile : pknoise made'
 
-testFFT : Objs/testFFT.o
-######
+Objs/pknoise : Objs/pknoise.o Objs/specpk.o Objs/mdish.o
+	$(CXXLINK) -o Objs/pknoise Objs/pknoise.o Objs/specpk.o Objs/mdish.o $(SOPHYAEXTSLBLIST)
 
-Objs/treccyl : Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o  Objs/brsrc.o
-	$(CXXLINK) -o Objs/treccyl \
-	Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o Objs/brsrc.o \
-	$(SOPHYAEXTSLBLIST)
+Objs/pknoise.o : pknoise.cc specpk.h  mdish.h 
+	$(CXXCOMPILE) -o Objs/pknoise.o pknoise.cc
 
-Objs/treccyl.o : treccyl.cc  multicyl.h mbeamcyl.h brsrc.h 
-	$(CXXCOMPILE) -o Objs/treccyl.o  treccyl.cc
+Objs/specpk.o : specpk.cc specpk.h  mdish.h 
+	$(CXXCOMPILE) -o Objs/specpk.o specpk.cc
 
-Objs/brsrc.o : brsrc.cc  brsrc.h 
-	$(CXXCOMPILE) -o Objs/brsrc.o  brsrc.cc
-
-Objs/mbeamcyl.o : mbeamcyl.cc  mbeamcyl.h brsrc.h 
-	$(CXXCOMPILE) -o Objs/mbeamcyl.o  mbeamcyl.cc
-
-Objs/multicyl.o : multicyl.cc  multicyl.h mbeamcyl.h brsrc.h 
-	$(CXXCOMPILE) -o Objs/multicyl.o  multicyl.cc
-
-testFFT :  Objs/testFFT.o
-	$(CXXLINK) -o testFFT \
-	Objs/testFFT.o \
-	$(SOPHYAEXTSLBLIST)
-
-Objs/testFFT.o : testFFT.cc
-	$(CXXCOMPILE) -o Objs/testFFT.o testFFT.cc
+Objs/mdish.o : mdish.cc mdish.h 
+	$(CXXCOMPILE) -o Objs/mdish.o mdish.cc
Index: trunk/Cosmo/RadioBeam/makerecyl
===================================================================
--- trunk/Cosmo/RadioBeam/makerecyl	(revision 3756)
+++ trunk/Cosmo/RadioBeam/makerecyl	(revision 3756)
@@ -0,0 +1,39 @@
+include $(SOPHYABASE)/include/sophyamake.inc
+
+all : treccyl
+
+clean :
+	rm Objs/*
+depend :
+	mkdir Objs
+
+treccyl : Objs/treccyl
+	echo 'makefile : treccyl made'
+
+testFFT : Objs/testFFT.o
+######
+
+Objs/treccyl : Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o  Objs/brsrc.o
+	$(CXXLINK) -o Objs/treccyl \
+	Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o Objs/brsrc.o \
+	$(SOPHYAEXTSLBLIST)
+
+Objs/treccyl.o : treccyl.cc  multicyl.h mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/treccyl.o  treccyl.cc
+
+Objs/brsrc.o : brsrc.cc  brsrc.h 
+	$(CXXCOMPILE) -o Objs/brsrc.o  brsrc.cc
+
+Objs/mbeamcyl.o : mbeamcyl.cc  mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/mbeamcyl.o  mbeamcyl.cc
+
+Objs/multicyl.o : multicyl.cc  multicyl.h mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/multicyl.o  multicyl.cc
+
+testFFT :  Objs/testFFT.o
+	$(CXXLINK) -o testFFT \
+	Objs/testFFT.o \
+	$(SOPHYAEXTSLBLIST)
+
+Objs/testFFT.o : testFFT.cc
+	$(CXXCOMPILE) -o Objs/testFFT.o testFFT.cc
Index: trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- trunk/Cosmo/RadioBeam/mdish.cc	(revision 3756)
+++ trunk/Cosmo/RadioBeam/mdish.cc	(revision 3756)
@@ -0,0 +1,309 @@
+//  Classes to compute 2D 
+// R. Ansari - Nov 2008, May 2010 
+
+#include "mdish.h"
+
+
+//--------------------------------------------------
+// -- Four2DResponse class 
+//--------------------------------------------------
+// Constructor
+Four2DResponse::Four2DResponse(int typ, double dx, double dy)
+  : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
+{
+}
+
+// Return the response for the wave vecteor (kx,ky)
+double Four2DResponse::Value(double kx, double ky)
+{
+  double wk,wkx,wky;
+  switch (typ_) 
+    {
+    case 1:   // Reponse gaussienne parabole diametre D exp[ - 0.5 (lambda  k_g / D )^2 ]
+      wk = sqrt(kx*kx+ky*ky)/dx_;
+      wk = 0.5*wk*wk;
+      return exp(-wk);
+      break;
+    case 2:   // Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda   
+      wk = sqrt(kx*kx+ky*ky)/dx_/2./M_PI;
+      return ( (wk<1.)?(1.-wk):0.);
+      break;
+    case 3:   // Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda) 
+      wkx = kx/2./M_PI/dx_; 
+      wky = ky/2./M_PI/dy_; 
+      return ( ((wkx<1.)&&(wky<1.))?((1.-wkx)*(1-wky)):0.);
+      break;
+    default:
+      return 1.;
+    }
+}
+// Return a vector representing the power spectrum (for checking) 
+Histo2D Four2DResponse::GetResponse(int nx, int ny)
+{
+  double kxmx = 1.2*DeuxPI*dx_;
+  double kymx = 1.2*DeuxPI*dy_;
+  if (typ_<3) kymx=kxmx; 
+  Histo2D h2(0.,kxmx,nx,0.,kymx,ny);
+
+  for(int j=0; j<h2.NBinY(); j++) 
+    for(int i=0; i<h2.NBinX(); i++) 
+      h2(i,j) = Value((i+0.5)*h2.WBinX(), (j+0.5)*h2.WBinY());
+  return h2;	
+}
+
+//---------------------------------------------------------------
+// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
+//---------------------------------------------------------------
+Four2DRespTable::Four2DRespTable(Histo2D const & hrep, double d)
+  : Four2DResponse(0,d,d) , hrep_(hrep)
+{
+}
+
+double Four2DRespTable::Value(double kx, double ky)
+{
+  int_4 i,j;
+  if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) || 
+       (ky<=hrep_.YMin())||(ky>=hrep_.YMax()) )  return 0.;
+  hrep_.FindBin(kx, ky, i, j);
+  return hrep_(i, j);
+}
+
+//--- Classe simple pour le calcul de rotation 
+class Rotation {
+public:
+  Rotation(double tet, double phi)
+  {
+// (Teta,Phi) = Direction de visee 
+// Les angles d'Euler correspondants sont Teta, Phi+Pi/2
+// Le Pi/2 vient que les rotations d'euler se font dans l'ordre
+//  Autour de oZ d'angle Phi, autour de oN (nouvel axe X) d'angle Teta
+//  Autour du nouvel axe Z (x3) d'angle Psi  (Psi=0 -> cp=1, sp=0.;
+  double ct = cos(tet);
+  double st = sin(tet);
+  // Le Pi/2 echange les axes X<>Y pour theta=phi=0 !
+  //  double cf = cos(phi+M_PI/2);
+  //  double sf = sin(phi+M_PI/2);
+  double cf = cos(phi);
+  double sf = sin(phi);
+  double cp = 1.; // cos((double)pO);
+  double sp = 0.; // sin((double)pO);
+  RE[0][0] = cf*cp-sf*ct*sp;     RE[0][1] = sf*cp+cf*ct*sp;      RE[0][2] = st*sp;
+  RE[1][0] = -cf*sp-sf*ct*cp;    RE[1][1] = -sf*sp+cf*ct*cp;     RE[1][2] = st*cp;
+  RE[2][0] = sf*st;              RE[2][1] = -cf*st;              RE[2][2] = ct;
+  }
+  inline void Do(double& x, double& y)
+  {
+    double xx=x; 
+    double yy=y;
+    x = RE[0][0]*xx+RE[0][1]*yy;
+    y = RE[1][0]*xx+RE[1][1]*yy;
+  }
+  double RE[3][3];
+};
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 2D 
+// -----------------------------------
+QHis2D::QHis2D()
+  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
+{
+}
+QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
+  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
+{
+  Define(xMin, xMax, nxb, yMin, yMax, nyb);
+}
+void QHis2D::Define(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
+{
+  nx=nxb; ny=nyb;
+  xmin=xMin;  xmax=xMax;  
+  ymin=yMin;  ymax=yMax;  
+  dxb=(xmax-xmin)/(double)nx;
+  dyb=(ymax-ymin)/(double)ny;
+  sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
+  aw.ReSize(2,sz);
+  sumw0=0.;
+  return;
+}
+double QHis2D::Add(r_8 x, r_8 y, r_8 w, bool fgfh)
+{
+  sa_size_t ix = (sa_size_t)((x-xmin)/dxb);
+  sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
+  if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
+  double rw = ((ix==0)&&(jy==0)) ? w : 0.;
+  sumw0 += rw; 
+  if (fgfh) aw(ix,jy) += w;
+  return rw;
+}
+Histo2D QHis2D::Convert()
+{
+  Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
+  for(int_4 j=0; j<ny; j++) 
+    for(int_4 i=0; i<nx; i++)  h2(i,j) = aw(i,j);
+  return h2;
+}
+
+//----------------------------------------------------------------------
+//  -- Pour calculer la reponse ds le plan kx,ky d'un system MultiDish 
+//----------------------------------------------------------------------
+MultiDish::MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto)
+  : lambda_(lambda), dmax_(dmax), dishes_(dishes), fgnoauto_(fgnoauto)
+{ 
+  SetThetaPhiRange();
+  SetRespHisNBins();
+  mcnt_=0;
+}
+
+Histo2D MultiDish::GetResponse()
+{
+  cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
+  double kmx = 1.2*DeuxPI*dmax_/lambda_;
+  /*
+  h2w_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
+  h2cnt_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
+  h2w_.Zero();
+  h2cnt_.Zero();
+  */
+  h2w_.Define(0.,kmx,nx_,0.,kmx,ny_);
+
+  double dold = dishes_[0].D/lambda_;
+  double dolx = dishes_[0].Dx/lambda_;
+  double doly = dishes_[0].Dy/lambda_;
+
+  Four2DResponse rd(2, dold, dold);
+  Four2DResponse rdr(3, dolx, doly);
+
+  if (!dishes_[0].isCircular())  rd = rdr;
+
+  double dtet = thetamax_/(double)ntet_;
+  double dphi = phimax_/(double)ntet_;
+
+  double sumw = 0.;
+  for(int kt=0; kt<ntet_; kt++) 
+    for(int jp=0; jp<nphi_; jp++) 
+      sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi);
+
+  double kx0 = DeuxPI*fabs(dishes_[1].X-dishes_[0].X)/lambda_;
+  double ky0 = DeuxPI*fabs(dishes_[1].Y-dishes_[0].Y)/lambda_;
+  int_4 ib, jb;
+  //  h2w_ /= h2cnt_; 
+  Histo2D h2 = h2w_.Convert();
+  h2.FindBin(kx0, ky0, ib, jb);
+  cout << " ---- MultiDish::GetResponse() VMin=" << h2.VMin() << " VMax= " << h2.VMax() 
+       << " h(0,0)=" << h2(0,0) << " h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
+  //  double fnorm=sqrt((double)dishes_.size())/h2.VMax();
+  double fnorm=1.;
+  if (h2.VMax() > sumw) {
+    fnorm=(double)dishes_.size()/h2.VMax();
+    cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " sumw=" << sumw 
+	 << " Renormalizing x NDishes/sumw  " << fnorm << endl;
+  }
+  else {
+    fnorm=(double)dishes_.size()/h2.VMax();
+    cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " VMax=" << h2.VMax()  
+	 << " Renormalizing x NDishes/h2.VMax()   " << fnorm << endl;
+  }
+  h2 *= fnorm;
+  cout << " ---- MultiDish::GetResponse() APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)=" 
+       << h2(0,0) << endl;
+  return h2;
+}
+
+/*
+double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
+{
+  double xxp = kx0+x;
+  double yyp = ky0+y;
+  double sumw=0.;
+  int_4 ib, jb;
+  h2w_.FindBin(xxp, yyp, ib, jb);
+  if ((ib==0)&&(jb==0))  sumw+=w;
+  if (fgfh) { 
+    h2w_.Add(xxp, yyp, w);
+    h2cnt_.Add(xxp, yyp, 1.);
+  }
+  double xxm=kx0-x;
+  double yym=ky0-y;
+  if (xxm>0.)  {
+    h2w_.FindBin(xxm, yyp, ib, jb);
+    if ((ib==0)&&(jb==0))  sumw+=w;
+    if (fgfh) {
+      h2w_.Add(xxm, yyp, w);
+      h2cnt_.Add(xxm, yyp, 1.);
+    }
+    if (yym>0.) { 
+      h2w_.FindBin(xxm, yym, ib, jb);
+      if ((ib==0)&&(jb==0))  sumw+=w;
+      if (fgfh) { 
+	h2w_.Add(xxm, yym, w);
+	h2cnt_.Add(xxm, yym, 1.);
+      }
+    }
+  }
+  if (yym>0.)  { 
+    h2w_.FindBin(xxp, yym, ib, jb);
+    if ((ib==0)&&(jb==0))  sumw+=w;
+    if (fgfh) {
+      h2w_.Add(xxp, yym, w);
+      h2cnt_.Add(xxp, yym, 1.);
+    }
+  }
+  return sumw; 
+}
+*/
+
+double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
+{
+  double xxp = kx0+x;
+  double yyp = ky0+y;
+  double sumw=0.;
+  sumw += h2w_.Add(xxp, yyp, w, fgfh);
+  double xxm=kx0-x;
+  double yym=ky0-y;
+  if (xxm>0.)  {
+    sumw += h2w_.Add(xxm, yyp, w, fgfh);
+    if (yym>0.)  sumw += h2w_.Add(xxm, yym, w, fgfh);
+  }
+  if (yym>0.)  sumw += h2w_.Add(xxp, yym, w, fgfh);
+  return sumw; 
+}
+
+double MultiDish::CumulResp(Four2DResponse& rd, double theta, double phi)
+{
+  //  cout << " MultiDish::CumulResp()  theta=" << theta << " phi=" << phi << endl;
+
+  double dx = h2w_.WBinX()/5;
+  double dy = h2w_.WBinY()/5;
+  int nbx = DeuxPI*rd.Dx()/dx;
+  int nby = DeuxPI*rd.Dy()/dy;
+  dx = DeuxPI*rd.Dx()/(double)nbx;
+  dy = DeuxPI*rd.Dy()/(double)nby;
+  if (mcnt_==0) 
+    cout << " CumulResp() nbx=" << nbx << " nby=" << nby << " dx=" << dx << " dy=" << dy << endl;
+  mcnt_++;
+
+  double sumw = 0.;
+  Rotation rot(theta, phi);
+
+  for(size_t i=0; i<dishes_.size(); i++) {
+    for(size_t j=i; j<dishes_.size(); j++) {
+      double kx0 = DeuxPI*fabs(dishes_[i].X-dishes_[j].X)/lambda_;
+      double ky0 = DeuxPI*fabs(dishes_[i].Y-dishes_[j].Y)/lambda_;
+      rot.Do(kx0, ky0);
+      if (kx0<0) kx0=-kx0;
+      if (ky0<0) ky0=-ky0;
+      bool fgfh= (!fgnoauto_ || (dishes_[i].ReflectorId()!=dishes_[j].ReflectorId()));
+      for(int ix=0; ix<nbx; ix++) 
+	for(int jy=0; jy<nby; jy++) { 
+	  double x = ix*dx;  
+	  double y = jy*dy;  
+	  sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
+	}
+    }
+    //    if (i%10==0) 
+    //      cout << " MultiDish::CumulResp() done i=" << i << " / imax=" << dishes_.size() 
+    //	   << " theta=" << theta << " phi=" << phi << endl;
+  }
+  return sumw;
+}
+
Index: trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- trunk/Cosmo/RadioBeam/mdish.h	(revision 3756)
+++ trunk/Cosmo/RadioBeam/mdish.h	(revision 3756)
@@ -0,0 +1,126 @@
+//  Classes to compute Multi-Dish or CRT-like radio interferometer response
+//  R. Ansari - Avril-Mai 2010 
+
+#ifndef MDISH_SEEN
+#define MDISH_SEEN
+
+#include "machdefs.h"      // SOPHYA .h  
+#include "sopnamsp.h"      // SOPHYA .h  
+#include <math.h>
+#include <iostream>
+#include <vector>
+#include <string> 
+
+#include "genericfunc.h"   // SOPHYA .h  
+#include "array.h"         // SOPHYA .h  
+#include "histats.h"       // SOPHYA .h 
+
+#define DeuxPI 2.*M_PI 
+
+// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (angles theta,phi)
+// typ=1 : Reponse gaussienne parabole diametre D exp[ - 0.5  (lambda  k_g / D )^2 ]
+// typ=2 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda   
+// typ=3 : Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda) 
+class Four2DResponse  {
+public:
+  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
+  Four2DResponse(int typ, double dx, double dy);
+
+  Four2DResponse(Four2DResponse const& a)
+  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; }
+  Four2DResponse& operator=(Four2DResponse const& a)
+  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }
+
+  // Return the 2D response for wave vector (kx,ky)
+  virtual double Value(double kx, double ky); 
+  inline  double operator()(double kx, double ky) 
+  { return Value(kx, ky); } 
+  virtual Histo2D GetResponse(int nx=256, int ny=256);  
+  inline double D() { return dx_; } ; 
+  inline double Dx() { return dx_; } ; 
+  inline double Dy() { return dy_; } ; 
+
+  int typ_;
+  double dx_, dy_; 
+};
+
+// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi) 
+class Four2DRespTable : public  Four2DResponse {
+public:
+  // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
+  Four2DRespTable(Histo2D const & hrep, double d);
+  // Return the 2D response for wave vector (kx,ky)
+  virtual double Value(double kx, double ky); 
+  Histo2D hrep_;
+};
+
+// Classe toute simple pour representer un element de reception de type dish
+class Dish {
+public:
+  // Circular dish 
+  Dish(int id, double x, double y, double diam)
+    :  id_(id), X(x), Y(y), D(diam), Dx(0.), Dy(0.), fgcirc_(true)   {   } 
+  // Receiver with rectangular type answer in kx,ky plane 
+  Dish(int id, double x, double y, double dx, double dy)
+    :  id_(id), X(x), Y(y), D(0.), Dx(dx), Dy(dy), fgcirc_(false)   {   } 
+  Dish(Dish const& a) 
+    :  id_(a.id_), X(a.X), Y(a.Y), D(a.D), Dx(a.Dx), Dy(a.Dy), fgcirc_(a.fgcirc_)     {   } 
+  inline bool isCircular() { return fgcirc_; }
+  inline int ReflectorId() { return id_; }
+
+  int id_;   // numero de reflecteur
+  double D,X,Y;
+  double Dx, Dy;
+  bool fgcirc_;  // false -> rectangular dish 
+};
+
+// -----------------------------------
+// -- Classe ressemblant a un histo 2D 
+class QHis2D {
+public:
+  QHis2D();
+  QHis2D(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
+  void Define(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
+  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
+  inline double WBinX() { return dxb; }
+  inline double WBinY() { return dyb; }
+  Histo2D Convert();
+  
+  r_8 xmin,xmax,ymin,ymax;
+  r_8 dxb,dyb;
+  sa_size_t nx,ny;
+  TArray<r_8> aw;
+  double sumw0;
+};
+
+// -------------------------------------------------------------------
+// -- Pour calculer la reponse ds le plan kx,ky d'un system MultiDish 
+class MultiDish {
+public:
+  MultiDish(double lambda, double dmax, vector<Dish>& dishes, bool fgnoauto=false);
+
+  inline void SetThetaPhiRange(double thetamax=0., int ntet=1, double phimax=0., int nphi=1)
+  { thetamax_=thetamax; ntet_=ntet; phimax_=phimax; nphi_=nphi; }
+
+  inline void SetRespHisNBins(int nx=128, int ny=128)
+  { nx_=nx; ny_=ny; }
+  Histo2D GetResponse();
+
+  double CumulResp(Four2DResponse& rd, double theta=0., double phi=0.);
+  inline size_t NbDishes() { return dishes_.size(); }
+
+  double AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh);
+
+  double lambda_, dmax_;
+  vector<Dish> dishes_;
+  bool fgnoauto_;
+  double thetamax_, phimax_; 
+  int ntet_,nphi_; 
+  int nx_, ny_;
+  //   Histo2D h2w_, h2cnt_;
+  QHis2D h2w_;
+  int mcnt_;
+};
+
+
+#endif 
Index: trunk/Cosmo/RadioBeam/pknoise.cc
===================================================================
--- trunk/Cosmo/RadioBeam/pknoise.cc	(revision 3756)
+++ trunk/Cosmo/RadioBeam/pknoise.cc	(revision 3756)
@@ -0,0 +1,289 @@
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include <iostream>
+#include <string>
+#include <math.h>
+
+#include <typeinfo>
+
+#include "mdish.h"
+#include "specpk.h"
+#include "histinit.h"
+// #include "fiosinit.h"
+// #include "fitsioserver.h"
+
+#include "randr48.h"      
+
+#include "timing.h"
+#include "ctimer.h"
+
+typedef DR48RandGen RandomGenerator ;
+
+// ---------------------------------------------------------------------
+// Test main program for computing interferometer noise power spectrum
+// R. Ansari  - Avril 2010 
+// ---------------------------------------------------------------------
+
+class PkNoiseCalculator
+{
+public:
+  PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, 
+		    double s2cut=100., int ngen=1, const char* tit="PkNoise")
+    : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
+  { } 
+  inline void SetS2Cut(double s2cut=100.)
+  {  S2CUT=s2cut; }
+  HProf Compute()
+  {
+    Timer tm(title.c_str());
+    tm.Nop();
+    HProf hnd;
+    cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
+    for(int igen=0; igen<NGEN; igen++) {
+      pkn3d.ComputeNoiseFourierAmp(frep);
+      if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
+      else pkn3d.ComputePkCumul(hnd,S2CUT);
+    }
+    return hnd;
+  }
+
+  Four3DPk& pkn3d;
+  Four2DResponse& frep;
+  double S2CUT;
+  int NGEN;
+  string title;
+};
+
+int main(int narg, const char* arg[])
+{
+  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  {
+    cout<< " Usage: pknoise [OutPPFName NGen S2Cut Lambda] " << endl;
+    cout<< "   Default: OutPPFName=pknoise.ppf, NGen=1 " << endl;
+    cout<< "            S2CUT=0. , Lambda=0.357  \n" << endl;
+
+    return 1;
+  }
+  cout << " ==== pknoise.cc program , test of SpectralShape and MassDist2D classes ==== " << endl;
+  // make sure SOPHYA modules are initialized 
+  SophyaInit();  
+  //  FitsIOServerInit();
+  InitTim();
+  //--- decoding command line arguments 
+  string outfile = "pknoise.ppf";  
+  if (narg>1) outfile = arg[1];
+  if (outfile==".")  outfile = "pknoise.ppf";
+  int NMAX = 1;
+  if (narg>2) NMAX = atoi(arg[2]);
+  double SCut=0.;
+  if (narg>3) SCut = atof(arg[3]);
+  double LAMBDA=0.357 ;  // 21 cm at z=0.7 
+  if (narg>4) LAMBDA = atof(arg[4]);
+
+  //-- end command line arguments
+  
+  int rc = 1;  
+  try {  // exception handling try bloc at top level
+    cout << "0/ pknoise.cc: Executing, output file= " << outfile << endl;  
+    POutPersist po(outfile);
+    cout << " 1.a/ Instanciating object type SpectralShape " << endl;
+    SpectralShape spec(2);
+    cout << " 1.b/ Wrinting spectral shape vector (name= Pk) to output PPF " << endl; 
+    Histo hpk = spec.GetPk(1024);
+    po << PPFNameTag("Pk") << hpk;
+
+    double D = 100.;
+    double lambda = LAMBDA;
+    double Dol = D/lambda;
+    cout << " 2.a/ Instanciating Four2DResponse(1/2/3...) " << endl;
+    Four2DResponse dishg(1,Dol,Dol);
+    Four2DResponse dish(2,Dol,Dol);
+    Four2DResponse dish2(2,Dol*2.,Dol*2.);
+    Four2DResponse dishsq(3,Dol,Dol/5.);
+    cout << " 2.b/ Writing Four2DResponse Histo2D to output ppf " << endl;
+    Histo2D hdg = dishg.GetResponse();
+    Histo2D hd = dish.GetResponse();
+    Histo2D hd2 = dish2.GetResponse();
+    Histo2D hdsq = dishsq.GetResponse();
+    po << PPFNameTag("dishg") << hdg;
+    po << PPFNameTag("dish") << hd;
+    po << PPFNameTag("dish2") << hd2;
+    po << PPFNameTag("dishsq") << hdsq;
+
+    cout << " 2.c/ Creating MultiDish Filled Array " << endl;
+    double Ddish = 5.;
+    double Ddish2 = 7.5;
+    double Eta=0.95;
+    int cnt=0;
+    vector<Dish> vdplein;
+    for(int i=0; i<20; i++) 
+      for(int j=0; j<20; j++) {
+	cnt++; 
+	vdplein.push_back(Dish(i*20+j+1, i*Ddish, j*Ddish, Eta*Ddish));
+      }
+    vector<Dish> vdsparse;
+    vector<Dish> vdsparseD7;
+
+    cnt=0;
+    for(int i=0; i<=18; i++) {
+      cnt++; vdsparse.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
+      vdsparseD7.push_back(Dish(cnt, i*Ddish2,0.,Eta*Ddish2));
+    }
+    for(int i=-18; i<=18; i++) {
+      if (i==0)  continue;
+      cnt++; vdsparse.push_back(Dish(cnt, 0.,i*Ddish,Eta*Ddish));
+      vdsparseD7.push_back(Dish(cnt, 0.,i*Ddish2,Eta*Ddish2));
+    }
+    for(int i=0; i<4; i++) {
+	cnt++;  vdsparse.push_back(Dish(cnt, (3+2*i)*Ddish,(3+2*i)*Ddish,Eta*Ddish));
+	vdsparseD7.push_back(Dish(cnt, (3+2*i)*Ddish2,(3+2*i)*Ddish2,Eta*Ddish2));
+	cnt++;	vdsparse.push_back(Dish(cnt, (3+2*i)*Ddish,-(3+2*i)*Ddish,Eta*Ddish));
+	vdsparseD7.push_back(Dish(cnt, (3+2*i)*Ddish2,-(3+2*i)*Ddish2,Eta*Ddish2));
+	/*
+	if ((i>0)||(j>0)) {
+	  cnt++;  vdsparse.push_back(Dish(cnt, (5+3*i)*Ddish,(3+2*j)*Ddish,Eta*Ddish));
+	  cnt++;  vdsparse.push_back(Dish(cnt, (5+3*i)*Ddish,-(3+2*j)*Ddish,Eta*Ddish));
+	}
+	*/
+      }
+
+
+    vector<Dish> vcylplein, vcylplP;
+    cnt=0;
+    double cylW=12.;   // Largeur des cylindres 
+    double cylRL=0.5;  // Longeur des elements de reception le long du cylindre
+    for(int i=0; i<8; i++) 
+      for(int j=0; j<192; j++) {
+	vcylplein.push_back(Dish(i+1, i*cylW, j*cylRL, 0.9*cylW, 0.9*cylRL));
+	cnt++; vcylplP.push_back(Dish(cnt, i*cylW, j*cylRL, 1.*cylW, 1.*cylRL));
+      }
+    vector<Dish> v2cyl, v2cylP;
+    cnt=0;
+    for(int i=0; i<2; i++) 
+      for(int j=0; j<32; j++) {
+	v2cyl.push_back(Dish(i+1, i*25, j*cylRL, 0.9*9., 0.9*cylRL));
+	cnt++; v2cylP.push_back(Dish(cnt, i*25, j*cylRL, 0.9*9., 1.*cylRL));
+      }
+    double LMAX = D;
+    bool fgnoauto = true;
+    int NRX=160;
+    int NRY=160;
+
+    MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
+    mdfill.SetRespHisNBins(NRX,NRY);   
+    Histo2D hrfill = mdfill.GetResponse();
+    PrtTim("Apres mdfill.GetResponse()");
+
+    MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto); 
+    mdsparse.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
+    mdsparse.SetRespHisNBins(NRX,NRY);
+    Histo2D hrsp = mdsparse.GetResponse();
+    PrtTim("Apres mdsparse.GetResponse()");
+
+    MultiDish mdsparsefp(lambda, LMAX, vdsparse, fgnoauto); 
+    mdsparsefp.SetRespHisNBins(NRX,NRY);
+    Histo2D hrspfp = mdsparsefp.GetResponse();
+    PrtTim("Apres mdsparsefp.GetResponse()");
+
+    MultiDish mdsparseD7(lambda, LMAX, vdsparseD7, fgnoauto); 
+    mdsparseD7.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
+    mdsparseD7.SetRespHisNBins(NRX,NRY);
+    Histo2D hrspd7 = mdsparseD7.GetResponse();
+    PrtTim("Apres mdsparseD7.GetResponse()");
+			     
+    MultiDish mcylfill(lambda, LMAX, vcylplein, fgnoauto);
+    mcylfill.SetRespHisNBins(NRX,NRY);
+    Histo2D hfcyl = mcylfill.GetResponse();
+    PrtTim("Apres mcylfill.GetResponse()");
+    MultiDish mcylfillP(lambda, LMAX, vcylplP, fgnoauto);
+    mcylfillP.SetRespHisNBins(NRX,NRY);
+    Histo2D hfcylP = mcylfillP.GetResponse();
+
+    MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
+    md2cyl.SetRespHisNBins(NRX,NRY);
+    Histo2D h2cyl = md2cyl.GetResponse();
+    MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
+    md2cylP.SetRespHisNBins(NRX,NRY);
+    Histo2D h2cylP = md2cylP.GetResponse();
+
+    po << PPFNameTag("mfill") << hrfill;
+    po << PPFNameTag("mspars") << hrsp;
+    po << PPFNameTag("msparsfp") << hrspfp;
+    po << PPFNameTag("msparsd7") << hrspd7;
+    po << PPFNameTag("mcylf") << hfcyl;
+    po << PPFNameTag("m2cyl") << h2cyl;
+    po << PPFNameTag("mcylfP") << hfcylP;
+    po << PPFNameTag("m2cylP") << h2cylP;
+
+    PrtTim("Done computing multi-dish response");
+ 
+
+    Four2DRespTable mdf(hrfill, Dol);
+    Four2DRespTable mds(hrsp,   Dol);
+    Four2DRespTable mdsfp(hrspfp,   Dol);
+    Four2DRespTable mdsd7(hrspd7,   Dol);
+
+    Four2DRespTable mcylf(hfcyl, Dol);
+    Four2DRespTable m2cyl(h2cyl, Dol);
+    Four2DRespTable mcylfP(hfcylP, Dol);
+    Four2DRespTable m2cylP(h2cylP, Dol);
+
+    
+    cout << " 3.a/ Instanciating object type Four3DPk  " << endl;
+    RandomGenerator rg;
+    Four3DPk m3d(rg);
+    m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI); 
+    cout << " 3.b/ call ComputeFourierAmp()  NGEN=" << NMAX << endl;
+    HProf hrpk;
+    for(int igen=0; igen<NMAX; igen++) {
+      m3d.ComputeFourierAmp(spec);
+      if (igen==0)  hrpk = m3d.ComputePk();
+      else m3d.ComputePkCumul(hrpk);
+    }
+    PrtTim("md.ComputeFourierAmp() done");
+    po << PPFNameTag("recPk") << hrpk;
+
+    cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
+#define NCONFIG 10
+    Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mds, &mdsfp, &mdsd7, &mcylf, &mcylfP, &m2cyl, &m2cylP};
+    const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","S63Dish5m","S63Dish5mFP","S63Dish7m",
+			       "F8Cyl","F8CylP","BiCyl","BiCylP"};
+    const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemds","noisemdsfp","noisemdsd7",
+			       "noisefcyl","noisefcylP","noise2cyl","noise2cylP"};
+    vector<int> nbdishes;
+    nbdishes.push_back(1);
+    nbdishes.push_back(1);
+    nbdishes.push_back(vdplein.size());
+    nbdishes.push_back(vdsparse.size());
+    nbdishes.push_back(vdsparse.size());
+    nbdishes.push_back(vdsparseD7.size());
+    nbdishes.push_back(vcylplein.size());
+    nbdishes.push_back(vcylplP.size());
+    nbdishes.push_back(v2cyl.size());
+    nbdishes.push_back(v2cylP.size());
+
+    for(int lc=0; lc<NCONFIG; lc++) {
+      PkNoiseCalculator pkn(m3d, *(f2rep[lc]), SCut/(double)nbdishes[lc], NMAX, tits[lc]);
+      HProf hpn = pkn.Compute();
+      po << PPFNameTag(tags[lc]) << hpn;
+    } 
+    rc = 0;
+  }  // End of try bloc 
+  catch (PThrowable & exc) {  // catching SOPHYA exceptions
+    cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
+         << "\n...exc.Msg= " << exc.Msg() << endl;
+    rc = 99;
+  }
+  catch (std::exception & e) {  // catching standard C++ exceptions
+    cerr << " pknoise.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
+    rc = 98;
+  }
+  catch (...) {  // catching other exceptions
+    cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
+    rc = 97;
+  }
+  PrtTim("End-pknoise");
+  cout << " ==== End of pknoise.cc program  Rc= " << rc << endl;
+  return rc;	
+}
+
+
Index: trunk/Cosmo/RadioBeam/plpkn.pic
===================================================================
--- trunk/Cosmo/RadioBeam/plpkn.pic	(revision 3756)
+++ trunk/Cosmo/RadioBeam/plpkn.pic	(revision 3756)
@@ -0,0 +1,123 @@
+########################################################################
+### Script de trace de PNoise(k) et reponse dans le plan (u,v) 
+###  de l'interferometre a partir du fichier PPF produit par pknoise.cc
+###     Fev - Avril 2010 ,  BAORadio/Reza
+########################################################################
+
+openppf pknoise.ppf
+
+setaxesatt 'minorticks font=helvetica,bold,16 autofontsize'
+disp noiseD 'logy nsta' 
+disp noiseD2 'same grey nsta'
+disp noisemdf 'same orange nsta'
+disp noisemds 'same red nsta'
+disp noisemdsfp 'same yellow nsta'
+disp noisemdsd7 'same gold nsta'
+disp noisefcyl 'same blue nsta'
+disp noisefcylP 'same skyblue nsta'
+disp noise2cyl 'same forestgreen nsta'
+disp noise2cylP 'same green nsta'
+
+Rad2Deg = 180/3.141596
+R2D2 = $Rad2Deg*$Rad2Deg
+Tsys = 50 
+Delnu = 1.e6
+
+#  A z = 0.7
+z = 0.7
+c = 3.e5
+H = 100 
+Da = 2500 
+nu21 = 1.42e9
+nu = $nu21/(1+$z)
+pi23 = 8.*Pi*Pi*Pi
+
+Lambda = 0.357 
+Lam2 = $Lambda*$Lambda 
+
+Dx = 100
+Dy = 100
+
+defscript calcul 
+  Dxol = $Dx/$Lambda 
+  Dyol = $Dy/$Lambda 
+
+  FOV = (1.2*1.2*$Lam2/$Dx/$Dy) 
+  FOVDEG = $FOV*$R2D2 
+  NPointing = 10000/$FOVDEG
+  tinteg = 365*86400/$NPointing
+  PNOISE = $Tsys*$Tsys/$tinteg/$Dxol/$Dyol
+  PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu
+#   PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu/$pi23
+  PNOISE = $PNOISE*1.e6
+  echo " FOV = $FOV deg^2  NPointing= $NPointing"
+  echo " tinteg= $tinteg sec  PNOISE= $PNOISE mK^2"
+endscript 
+
+Dx = 100
+Dy = 100
+calcul 
+plot2d noiseD x/$Da val*$PNOISE nb>10 'line=solid,2 logy logx xylimits=0.002,0.8,1,1e5 grid cpts nsta notit' 
+Dx = 200
+Dy = 200
+calcul 
+plot2d noiseD2 x/$Da val*$PNOISE nb>10 'same line=solid,2  cpts grey nsta notit' 
+Dx = 5*0.95
+Dy = 5*0.95
+calcul 
+plot2d noisemdf x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts orange nsta notit' 
+plot2d noisemds x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
+Dx = 7.55*0.95
+Dy = 7.5*0.95
+calcul 
+plot2d noisemdsd7 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts siennared nsta notit' 
+Dx = 12*0.9
+Dy = 0.5*0.9
+calcul 
+plot2d noisefcyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts blue nsta notit' 
+plot2d noisefcylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts skyblue nsta notit' 
+Dx = 8*0.9
+Dy = 0.5*0.9
+calcul 
+plot2d noise2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts forestgreen nsta notit' 
+plot2d noise2cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit' 
+
+set lines ( '100mDish' '200mDish' '400x5mDish' '63x5mDish' '63x7.5mDish'  '8Cyl12x96m' 'Perf8Cyl12x96m' )
+set lines ( $lines 'Pitts2Cyl' 'PerfPitts2Cyl' ) 
+set cols ( black grey  orange red siennared  blue skyblue forestgreen green )
+textdrawer lines cols 'frame font=helvetica,bold,16' 
+setaxelabels 'k (Mpc^-1) ' 'PNoise(k) mk^2 Mpc^3' 'font=helvetica,bolditalic,16' 
+
+
+newwin
+disp mfill 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , 400 x 5m Dishes - No Pointing' ' ' 'font=helvetica,bold,16' 
+
+disp dish 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , 100 m Dish' ' ' 'font=helvetica,bold,16' 
+
+disp msparsfp 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , 63 x 5m Dishes T-config - No Pointing' ' ' 'font=helvetica,bold,16' 
+
+disp mspars 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , 63 x 5m Dishes T-config - Pointing ~Pi/4' ' ' 'font=helvetica,bold,16' 
+
+disp mcylf 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , Filled Cylinder Array 8 Cyl 12mx96m ' ' ' 'font=helvetica,bold,16' 
+
+disp mcylfP 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , Perfect Filled Cylinder Array 8 Cyl 12mx96m ' ' ' 'font=helvetica,bold,16' 
+
+disp m2cyl 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , Pittsburgh 2 Cyl 16mx8m , dist=25m ' ' ' 'font=helvetica,bold,16' 
+
+disp m2cylP 'h2disp=img colbr128 notit'
+setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
+settitle ' u-v coverage , Perfect-Pitts. 2 Cyl 16mx8m , dist=25m ' ' ' 'font=helvetica,bold,16' 
Index: trunk/Cosmo/RadioBeam/specpk.cc
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.cc	(revision 3756)
+++ trunk/Cosmo/RadioBeam/specpk.cc	(revision 3756)
@@ -0,0 +1,345 @@
+//  Classes to compute 3D power spectrum 
+// R. Ansari - Nov 2008, May 2010 
+
+#include "specpk.h"
+#include "randr48.h"      
+
+//------------------------------------
+// Class SpectralShape 
+// -----------------------------------
+
+double Pnu1(double nu) 
+{
+  return ( sqrt(sqrt(nu)) / ((nu+1.0)/0.2) * 
+	   (1+0.2*cos(2*M_PI*(nu-2.)*0.15)*exp(-nu/50.)) );
+}
+
+double Pnu2(double nu) 
+{
+  if (nu < 1.e-9) return 0.;
+  return ((1.-exp(-nu/0.5))/nu*(1+0.25*cos(2*M_PI*nu*0.1)*exp(-nu/20.)) );
+}
+
+
+double Pnu3(double nu) 
+{
+  return ( log(nu/100.+1)*(1+sin(2*M_PI*nu/300))*exp(-nu/4000) );
+}
+
+
+double Pnu4(double nu) 
+{
+  double x = (nu-0.5)/0.05;
+  double rc = 2*exp(-x*x);
+  x = (nu-3.1)/0.27;
+  rc += exp(-x*x);
+  x = (nu-7.6)/1.4;
+  rc += 0.5*exp(-x*x);
+  return ( rc+2.*exp(-x*x) );
+}
+
+//--------------------------------------------------
+// -- SpectralShape class : test P(k) class
+//--------------------------------------------------
+// Constructor
+SpectralShape::SpectralShape(int typ)
+{
+  typ_=typ;
+}
+
+// Return the spectral power for a given wave number wk 
+double SpectralShape::operator() (double wk)
+{
+  wk/=DeuxPI;
+  switch (typ_) 
+    {
+    case 1:
+      return Pnu1(wk);
+      break;
+    case 2:
+      return Pnu2(wk);
+      break;
+    case 3:
+      return Pnu3(wk);
+      break;
+    case 4:
+      return Pnu4(wk);
+      break;
+    default :
+      {
+  // global shape
+      double csp = pow( (2*sin(sqrt(sqrt(wk/7.)))),2.);
+      if (csp < 0.) return 0.;
+      
+      // Adding some pics
+      double picpos[5] = {75.,150.,225.,300.,375.,};
+      
+      for(int k=0; k<5; k++) {
+  	double x0 = picpos[k];
+  	if ( (wk > x0-25.) && (wk < x0+25.) ) {
+	  double x = (wk-x0);
+	  csp *= (1.+0.5*exp(-(x*x)/(2.*5*5)));
+	  break;
+  	}
+      }
+      return csp;
+      }
+      break;
+    }
+}
+// Return a vector representing the power spectrum (for checking) 
+Histo SpectralShape::GetPk(int n)
+{
+  if (n<16) n = 256;
+  Histo h(0.,1024.*DeuxPI,n);
+  for(int k=0; k<h.NBins(); k++)   h(k) = Value((k+0.5)*h.BinWidth());
+  return h;	
+}
+
+//--------------------------------------------------
+// -- Four2DResponse class : test P(k) class
+
+//---------------------------------------------------------------
+// -- Four3DPk class :  3D fourier amplitudes and power spectrum 
+//---------------------------------------------------------------
+// Constructeur avec Tableau des coeff. de Fourier en argument
+Four3DPk::Four3DPk(TArray< complex<TF> > & fourcoedd, RandomGeneratorInterface& rg)
+  : rg_(rg), fourAmp(fourcoedd)
+{
+  SetPrtLevel();
+  SetCellSize();
+}
+// Constructor
+Four3DPk::Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx, sa_size_t szy, sa_size_t szz)
+  : rg_(rg), fourAmp(szx, szy, szz) 
+{
+  SetPrtLevel();
+  SetCellSize();
+}
+
+
+// Generate mass field Fourier Coefficient
+void Four3DPk::ComputeFourierAmp(SpectralShape& pk)
+{
+  // We generate a random gaussian real field  
+  // fourAmp represent 3-D fourier transform of a real input array. 
+  // The second half of the array along Y and Z contain negative frequencies
+  //  double fnorm = 1./sqrt(2.*fourAmp.Size()); 
+  double fnorm = 1.; 
+  double kxx, kyy, kzz;
+  // sa_size_t is large integer type  
+  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+    kzz =  (kz>fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_; 
+    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+      kyy =  (ky>fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+	double kxx=(double)kx*dkx_;
+	double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
+	double amp = sqrt(pk(wk)*fnorm/2.);      
+	fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   // renormalize fourier coeff usin 
+      }
+    }
+  }
+  if (prtlev_>0)
+    cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
+}
+
+// Generate mass field Fourier Coefficient
+void Four3DPk::ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask)
+{
+  TMatrix<r_4> mask(fourAmp.SizeY(), fourAmp.SizeX());
+  // fourAmp represent 3-D fourier transform of a real input array. 
+  // The second half of the array along Y and Z contain negative frequencies
+  double kxx, kyy, kzz;
+  // sa_size_t is large integer type  
+  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+    kzz =  (kz>fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_; 
+    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+      kyy =  (ky>fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+	double kxx=(double)kx*dkx_;
+	double rep = resp(kxx, kyy);
+	if (crmask&&(kz==0))  mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
+	if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
+	else {
+	  double amp = 1./sqrt(rep)/sqrt(2.);
+	  fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
+	}
+      }
+    }
+  }
+  if (prtlev_>1)  fourAmp.Show();
+  if (crmask) {
+    POutPersist po("mask.ppf");
+    po << mask;
+  }
+  if (prtlev_>0)
+    cout << " Four3DPk::ComputeNoiseFourierAmp() done ..." << endl;
+}
+
+// Compute mass field from its Fourier Coefficient
+TArray<TF>  Four3DPk::ComputeMassDens()
+{
+  TArray<TF> massdens;
+// Backward fourier transform of the fourierAmp array   
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  ffts.FFTBackward(fourAmp, massdens, true);
+  //  cout << " Four3DPk::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" <<  massDens.Size() << endl;
+  cout << " Four3DPk::ComputeMassDens() done NPix=" <<  massdens.Size() << endl;
+  return massdens;
+}
+
+// Compute power spectrum as a function of wave number k 
+// cells with amp^2=re^2+im^2>s2cut are ignored
+// Output : power spectrum (profile histogram)
+HProf Four3DPk::ComputePk(double s2cut, double kmin, double kmax, int nbin)
+{
+  // The second half of the array along Y (matrix rows) contain
+  // negative frequencies
+  //  int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
+  // The profile histogram will contain the mean value of FFT amplitude
+  // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
+  //  if (nbin < 1) nbin = nbh/2;
+  HProf hp(kmin, kmax, nbin);
+  hp.SetErrOpt(false);
+  ComputePkCumul(hp, s2cut);
+  return hp;
+}
+
+// Compute power spectrum as a function of wave number k 
+// Cumul dans hp - cells with amp^2=re^2+im^2>s2cut are ignored
+void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
+{
+
+  // fourAmp represent 3-D fourier transform of a real input array. 
+  // The second half of the array along Y and Z contain negative frequencies
+  double kxx, kyy, kzz;
+  // sa_size_t is large integer type  
+  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
+    kzz =  (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_; 
+    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
+      kyy =  (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_; 
+      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
+	double kxx=(double)kx*dkx_;
+	complex<TF> za = fourAmp(kx, ky, kz);
+	if (za.real()>8.e9) continue;
+	double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
+	double amp2 = za.real()*za.real()+za.imag()*za.imag();
+	if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
+	hp.Add(wk, amp2);
+      }
+    }
+  }
+  return;
+}
+
+
+
+//-----------------------------------------------------
+// -- MassDist2D class :  2D mass distribution 
+//-----------------------------------------------------
+// Constructor
+MassDist2D::MassDist2D(GenericFunc& pk, int size, double meandens) 
+: pkSpec(pk) , sizeA((size>16)?size:16) ,  massDens(sizeA, sizeA), 
+  meanRho(meandens) , fg_fourAmp(false) , fg_massDens(false)
+{
+}
+
+// To the computation job
+void MassDist2D::Compute()
+{
+  ComputeFourierAmp();
+  ComputeMassDens(); 	
+}
+
+// Generate mass field Fourier Coefficient
+void MassDist2D::ComputeFourierAmp()
+{
+  if (fg_fourAmp) return; // job already done
+  // We generate a random gaussian real field  
+  double sigma = 1.;
+// The following line fills the array by gaussian random numbers  
+//--Replaced--  massDens = RandomSequence(RandomSequence::Gaussian, 0., sigma);
+// Can be replaced by 
+  DR48RandGen rg;
+  for(sa_size_t ir=0; ir<massDens.NRows(); ir++) {
+  	for(sa_size_t jc=0; jc<massDens.NCols(); jc++) {
+      massDens(ir, jc) = rg.Gaussian(sigma); 
+  	}
+  }
+// --- End of random filling
+
+  // Compute fourier transform of the random gaussian field -> white noise 
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  ffts.FFTForward(massDens, fourAmp);
+    
+  // fourAmp represent 2-D fourier transform of a real input array. 
+  // The second half of the array along Y (matrix rows) contain
+  // negative frequencies
+//  double fnorm = 1./sqrt(2.*fourAmp.Size()); 
+// PUT smaller value for fnorm and check number of zeros
+  double fnorm = 1.; 
+  // sa_size_t is large integer type  
+  for(sa_size_t ky=0; ky<fourAmp.NRows(); ky++) {
+    double kyy = ky;
+    if (ky > fourAmp.NRows()/2) kyy = fourAmp.NRows()-ky;  // negative frequencies 
+    for(sa_size_t kx=0; kx<fourAmp.NCols(); kx++) {
+      double wk = sqrt((double)(kx*kx+kyy*kyy));
+      double amp = pkSpec(wk)*fnorm;      
+      fourAmp(ky, kx) *= amp;   // renormalize fourier coeff using 
+    }
+  }
+  fg_fourAmp = true;
+  cout << " MassDist2D::ComputeFourierAmp() done ..." << endl;
+}
+
+// Compute mass field from its Fourier Coefficient
+void MassDist2D::ComputeMassDens()
+{
+  if (fg_massDens) return; // job already done
+  if (!fg_fourAmp) ComputeFourierAmp();   // Check fourier amp generation
+
+// Backward fourier transform of the fourierAmp array   
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  ffts.FFTBackward(fourAmp, massDens, true);
+// We consider that massDens represents delta rho/rho 
+// rho = (delta rho/rho + 1) * MeanDensity 
+  massDens += 1.;
+// We remove negative values 
+  sa_size_t npbz = 0;
+  for (sa_size_t i=0; i<massDens.NRows(); i++) 
+    for (sa_size_t j=0; j<massDens.NCols(); j++) 
+      if (massDens(i,j) < 0.) { npbz++; massDens(i,j) = 0.; }
+  massDens *= meanRho;
+  cout << " MassDist2D::ComputeMassDens() done NbNeg=" << npbz << " / NPix=" <<  massDens.Size() << endl;
+}
+
+// Compute power spectrum as a function of wave number k 
+// Output : power spectrum (profile histogram)
+HProf MassDist2D::ReconstructPk(int nbin)
+{
+  // The second half of the array along Y (matrix rows) contain
+  // negative frequencies
+  int nbh = sqrt(2.0)*fourAmp.NCols();
+  // The profile histogram will contain the mean value of FFT amplitude
+  // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
+  if (nbin < 1) nbin = nbh/2;
+  HProf hp(-0.5, nbh-0.5, nbin);
+  hp.SetErrOpt(false);
+
+  for(int ky=0; ky<fourAmp.NRows(); ky++) {
+    double kyy = ky;
+    if (ky > fourAmp.NRows()/2)  kyy = fourAmp.NRows()-ky;  // negative frequencies
+    for(int kx=0; kx<fourAmp.NCols(); kx++) {
+      double wk = sqrt((double)(kx*kx+kyy*kyy));
+      complex<r_8> za = fourAmp(ky, kx);
+      double amp = sqrt(za.real()*za.real()+za.imag()*za.imag());
+      hp.Add(wk, amp);
+    }
+  }
+  return hp;
+}
+
Index: trunk/Cosmo/RadioBeam/specpk.h
===================================================================
--- trunk/Cosmo/RadioBeam/specpk.h	(revision 3756)
+++ trunk/Cosmo/RadioBeam/specpk.h	(revision 3756)
@@ -0,0 +1,107 @@
+//  Class examples to generate mass distribution 
+// R.A. for A. Abate , Nov. 2008
+
+#ifndef SPECPK_SEEN
+#define SPECPK_SEEN
+
+#include "machdefs.h"      
+#include "sopnamsp.h"       
+#include <math.h>
+#include <iostream>
+#include <vector>
+#include <string> 
+
+#include "genericfunc.h"  
+#include "array.h"         
+#include "histats.h"       
+#include "fftwserver.h"    
+#include "randinterf.h"      
+
+#include "mdish.h"       
+
+#define DeuxPI 2.*M_PI 
+
+// -- SpectralShape class : test P(k) class
+class SpectralShape : public GenericFunc {
+public:
+  SpectralShape(int typ);
+// Return the value of power spectrum for wave number wk   
+  virtual double operator() (double wk); 
+  inline  double Value(double wk) { return((*this)(wk)); }
+// Return a vector representing the power spectrum (for checking) 
+  Histo GetPk(int n=256);      	
+  int typ_;
+};
+
+
+#define TF  r_4 
+
+// -- Four3DPk class :  3D fourier amplitudes and power spectrum 
+class Four3DPk {
+public:
+// Constructor
+  Four3DPk(TArray< complex<TF> > & fourcoedd, RandomGeneratorInterface& rg);
+  Four3DPk(RandomGeneratorInterface& rg, sa_size_t szx=128, sa_size_t szy=256, sa_size_t szz=128);
+  inline void SetCellSize(double dkx=DeuxPI, double dky=DeuxPI, double dkz=DeuxPI) 
+  { dkx_=dkx;  dky_=dky;  dkz_=dkz; } 
+  inline int SetPrtLevel(int lev=0) 
+  { int olev=prtlev_; prtlev_=lev; return olev; }
+  void ComputeFourierAmp(SpectralShape& pk);
+  void ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask=false);
+// Return the array size 
+  inline sa_size_t NCells() { return fourAmp.Size(); }
+// Set the cell size/step in Fourier Space
+// Return the fourier amplitude matrix  
+  TArray< complex<TF> > GetFourierAmp()
+    { return fourAmp; }
+// Return the mass density matrix    
+  TArray<TF> ComputeMassDens(); 
+
+// Return the reconstructed power spectrum as a profile histogram   
+//  HProf ComputePk(double s2cut=0., double kmin=-0.5*DeuxPI, double kmax=127.5*DeuxPI, int nbin=256);
+  HProf ComputePk(double s2cut=0., double kmin=-DeuxPI, double kmax=255*DeuxPI, int nbin=256);
+  void  ComputePkCumul(HProf& hp, double s2cut=0.);
+
+protected:
+  // member attribute
+  RandomGeneratorInterface& rg_;
+  TArray< complex<TF> > fourAmp;  // complex array of fourier coefficients
+  double dkx_, dky_, dkz_;
+  int prtlev_;
+};
+
+
+// -- MassDist2D class :  2D mass distribution 
+class MassDist2D {
+public:
+// Constructor
+  MassDist2D(GenericFunc& pk, int size=1024, double meandens=1.);
+// Do the computation 
+  void Compute();
+// Return the array size 
+  inline sa_size_t ArrSize() { return sizeA; }
+// Return the fourier amplitude matrix  
+  TMatrix< complex<r_8> > GetFourierAmp()
+    { if (!fg_fourAmp) ComputeFourierAmp(); return fourAmp; }
+// Return the mass density matrix    
+  Matrix GetMassDens()
+    { if (!fg_massDens) ComputeMassDens(); return massDens; }
+
+// Return the reconstructed power spectrum as a profile histogram   
+  HProf ReconstructPk(int nbin=0);
+protected:
+  void ComputeFourierAmp();
+  void ComputeMassDens();
+  
+// member attribute
+  GenericFunc& pkSpec;   // The spectralShape
+  sa_size_t sizeA;       // 2D array size
+  double meanRho;       // Mean Density
+  bool fg_fourAmp;  // true -> fourAmp computed 
+  TMatrix< complex<r_8> > fourAmp;  // complex array of fourier coefficients
+  bool fg_massDens;  // true -> MassDens computed 
+  TMatrix< r_8 > massDens;  	// real array of d rho/rho 
+};
+
+
+#endif 
