Index: trunk/Cosmo/RadioBeam/README
===================================================================
--- trunk/Cosmo/RadioBeam/README	(revision 3949)
+++ trunk/Cosmo/RadioBeam/README	(revision 3973)
@@ -26,5 +26,5 @@
 
 A.3/ List of other common files :
-  - fgndsub.h fgndsub.h : radio source, synchrotron subtraction class 
+  - fgndsub.h fgndsub.cc : radio source, synchrotron subtraction class 
   - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
   - specpk.h specpk.cc : 3D power spectrum computation classes  
Index: trunk/Cosmo/RadioBeam/applobe.cc
===================================================================
--- trunk/Cosmo/RadioBeam/applobe.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/applobe.cc	(revision 3973)
@@ -3,6 +3,14 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-  Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName 
-                 [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] 
+  Usage: applobe [-g -t -fib]  Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName 
+                 [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] 
+   o -t : triangle beam shape in k-space
+   o -g : gaussian beam shape in k-space 
+   o -fib: application d'un lobe fixe (independant de la frequence) 
+   o Diameter/Four2DRespTableFile : diametre dish ou nom de fichier PPF reponse 2D
+   o In3DPPFName Out3DPPFName : Noms des fichiers (TArray<r_4> 3D) entree-sortie 
+   o TargetBeamDiam : Correction à un beam fiduciel, independant de la frequence 
+     avec specification de la valeur du diametre pour la frequence la plus basse
+   o ResmapleFactor : reechantillonnage du cube (2 -> surechantillonnage, 0.5 : 1 sur 2)
 ---------------------------------------------------------------  */
 
@@ -48,13 +56,40 @@
 
   InitTim();   // Initializing the CPU timer
+
+  if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
+    cout << "Usage: applobe [-t -g -fib] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n" 
+	 << "               [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
+    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
+      cout << "   -t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
+	   << "   -fib :  Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n" 
+           << "   Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n" 
+           << "   In3DPPFName Out3DPPFName: Input/Output PPF file names  (TArray<r_4> 3D) \n" 
+           << "   TargetBeamDiam: Corrected beam target diameter (D/Lambda at lowest frequency) \n" 
+           << "      DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm)) \n"
+	   << "   ResmapleFactor : Resampling (2-> oversampling, 0.5 : 1/2 undersampling) \n" << endl;    
+    }  
+    return 1;
+  }
+
   Timer tm("applobe");
 
+  // decodage arguments 
+  bool fggaussian=true;  // true -> gaussian beam 
+  bool fixedbeam=false;  // true -> apply freq. independent beam 
+
+  // decodage argument optionnel 
+  bool fgoptarg=true;
+  while (fgoptarg) {
+    string fbo = arg[1];
+    if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
+    else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
+    else if (fbo=="-fib")  { fixedbeam=true; arg++; narg--; }
+    else fgoptarg=false;
+  }
   if (narg < 3) {
-    cout << "Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName " 
-	 << "             [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] \n" << endl;
-    return 1;
-  }
-
-  // decodage arguments 
+    cout << " applobe/error arguments , applobe -h for help " << endl;
+    return 2;
+  }
+
   bool fgresptbl=true;
   double DIAMETRE=100.;
@@ -64,13 +99,13 @@
     DIAMETRE=atof(arg[1]);
   }
-  else resptblname=arg[1];
-
+  else  resptblname=arg[1];
+  
   string inname = arg[2];
   string outname = arg[3];
   bool fgcorrbeam=false;
-  double tbeamDoL=135;
+  double tbeamDiam=50.;
   if (narg>4) {
-    tbeamDoL=atof(arg[4]);
-    if (tbeamDoL>1)  fgcorrbeam=true;
+    tbeamDiam=atof(arg[4]);
+    if (tbeamDiam>1.)  fgcorrbeam=true;
   }
   bool fgresample=false;
@@ -83,5 +118,6 @@
   int rc = 91;
 
-  cout << " ====== applobe :   Input skycube  name= " << inname << " OutName=" << outname;
+  cout << " ====== applobe :   Input skycube  name= " << inname << " OutName=" << outname << endl;
+  cout << ((fggaussian)?" Gaussian ":" Triangular") << ((fixedbeam)?" Fixed (freq.independent)":"") << " beams" << endl;
   bool fginmap=true;
   try {
@@ -113,5 +149,6 @@
     double lambda = conv.getLambda();
     
-    Four2DResponse fresp(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
+    int typbeam=(fggaussian)?1:2;
+    Four2DResponse fresp(typbeam, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
     Four2DResponse* fresp_p=&fresp;
     Four2DRespTable resptbl;
@@ -119,4 +156,5 @@
       cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
       resptbl.readFromPPF(resptblname);
+      resptbl.renormalize(1.);
       fresp_p=&resptbl;
     }
@@ -126,12 +164,16 @@
 
     if (fgcorrbeam) {
-      double DoL = tbeamDoL;
+      double DoL = tbeamDiam/lambda;
       double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
-      int typcb = 2;
+      int typcb = (fggaussian)?1:2;
       //      if (fgresptbl) typcb=22;
       Four2DResponse tbeam(typcb, DoL, DoL );
-      cout << "applobe[3]: calling Correct2RefLobe() with target beam D/Lambda=" << DoL  
-	   << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
+      cout << "applobe[3]: calling Correct2RefLobe() with target beam Diameter=" << tbeamDiam 
+	   << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
       beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq);
+    }
+    else if (fixedbeam) {
+      cout << "applobe[3]: calling ApplyLobe() (freq.independent beam) ... " << endl;
+      beam.ApplyLobe(incube, dx, dy, Freq0MHz);
     }
     else {
Index: trunk/Cosmo/RadioBeam/calcpk.cc
===================================================================
--- trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/calcpk.cc	(revision 3973)
@@ -102,5 +102,5 @@
 
 
-    HProf hp = pkc.ComputePk(0.,256);
+    HProf hp = pkc.ComputePk(0.,HPk_NBin);
     {
       cout << "calcpk[4] : writing profile P(k) to PPF file  " << outname << endl;
Index: trunk/Cosmo/RadioBeam/calcpk2.cc
===================================================================
--- trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3973)
@@ -7,5 +7,5 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile 
+Usage: calcpk2 [-t -g] InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile 
                [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]
 ---------------------------------------------------------------  */
@@ -44,9 +44,10 @@
 {
   if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
-    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
-         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDoL] \n"
+    cout << " Usage: [-t -g] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n" 
+         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
 	 << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
     if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
-      cout << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
+      cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n" 
+	   << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
 	   << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
 	   << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
@@ -54,11 +55,11 @@
 	   << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n" 
 	   << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
-	   << "- CorBeamDoL: Beam correction target Diameter/Lambda \n" 
+	   << "- CorBeamDiam: Beam correction target dish diameter  \n" 
 	   << "    These two parameters are used to correct for beam effect for a \n"
 	   << "    target beam (independent of frequency) defined by D/Lambda \n"
-	   << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5) \n" 
+	   << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm) \n" 
 	   << "    default : no beam correction applied \n " 
 	   << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n" 
-	   << "    default : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
+	   << "    default (0.) : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
 	   << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
 	   << "    foreground subtraction. default is P2 \n" 
@@ -74,4 +75,18 @@
   int rc = 0;
   try {
+    bool fggaussian=true;  // true -> gaussian beam 
+    // decodage argument optionnel 
+    bool fgoptarg=true;
+    while (fgoptarg) {
+      string fbo = arg[1];
+      if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
+      else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
+      else fgoptarg=false;
+    }
+    if (narg < 6) {
+      cout << " calcpk2/error arguments , applobe -h for help " << endl;
+      return 2;
+    }
+
     string inppflss = arg[1];
     r_4 rfaclss = atof(arg[2]);
@@ -102,8 +117,8 @@
       }
     }
-    double tbeamDoL=0.;
+    double tbeamDiam=0.;
     if (narg>8) {
-      tbeamDoL=atof(arg[8]);
-      if (tbeamDoL<1.)  fgcorrbeam=false;
+      tbeamDiam=atof(arg[8]);
+      if (tbeamDiam<1.)  fgcorrbeam=false;
     }
     bool fgclnsrc=true;
@@ -170,4 +185,5 @@
       cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
       resptbl.readFromPPF(resptblname);
+      resptbl.renormalize(1.);
       arep_p=&resptbl;
     }
@@ -175,7 +191,7 @@
 	      << " DoL=" << DIAMETRE/lambda << " ) " << endl;
 
-    double DoL = tbeamDoL;
+    double DoL = tbeamDiam/lambda;
     double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
-    int typcb = 2;
+    int typcb = (fggaussian)?1:2;
     //    if (fgresptbl) typcb=22;
     Four2DResponse tbeam(typcb, DoL, DoL );
@@ -183,6 +199,6 @@
     ForegroundCleaner  cleaner(*arep_p, tbeam, skycube); 
     if (fgcorrbeam) {
-      cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam D/Lambda=" << DoL 
-	   << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
+      cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam  
+	   << " D/Lambda=" << DoL  << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
       cleaner.BeamCorrections();
     }
@@ -218,5 +234,5 @@
     pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
     
-    HProf hp = pkc.ComputePk(0.,256);
+    HProf hp = pkc.ComputePk(0.,HPk_NBin);
 
     tm.Split(" Done ComputePk ");    
Index: trunk/Cosmo/RadioBeam/cubedef.h
===================================================================
--- trunk/Cosmo/RadioBeam/cubedef.h	(revision 3949)
+++ trunk/Cosmo/RadioBeam/cubedef.h	(revision 3973)
@@ -7,12 +7,13 @@
 // Definition tailles et bornes du cube angles (X,Y), frequence (Z)
 
-static sa_size_t NTheta=360;
-static sa_size_t NPhi=360;
+static sa_size_t NTheta=600;
+static sa_size_t NPhi=800;
 static sa_size_t NFreq = 256;
 
+// Carte 30x50 deg, centre en delta = +10 deg , alpha= 150 deg = 10h00
 static double Theta0Degre = 65.;  // -> Delta = +30 deg
-static double Phi0Degre = 135.;   // -> alpha = 9h00
+static double Phi0Degre = 130.;   // -> alpha = 8.66 heures 
 static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
-static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
+static double PhiSizeDegre = 40.;    // Taille de la carte en degre selon delta (axe Y)
 
 /* ---  Parametres avec generation Reza from HASLAM + NVSS 
@@ -34,11 +35,13 @@
 
 // z = 0.6, freq_center=884 MHz , dA_comov ~= 2185 Mpc ,  
-//       5 arcmin -> 3 Mpc , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
+//       5 arcmin -> 3 Mpc (3.17 Mpc) , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
+//       3 arcmin -> 1.9 Mpc, 0.5 MHz -> 2.8 Mpc 
 // Taille de cellule pour le calcul du spectre de puissance 3D 
 static double ComovDA = 2185.;
-static double XCellSizeMpc = 3.;
-static double YCellSizeMpc = 3.;
-static double ZCellSizeMpc = 3.;
-
+static double XCellSizeMpc = 1.9;
+static double YCellSizeMpc = 1.9;
+static double ZCellSizeMpc = 2.8;
+// Nb de bin de HProf pour calcul P(k)
+static int HPk_NBin = 384;
 
 
@@ -52,5 +55,6 @@
 static double sigPLidx2 = 0.15;
 
-static double PLidxSrc = -2.2;  // index de la loi de puissance des sources  
+// Generation de la loi de puissance des sources : 
+static double PLidxSrc = -2.0;  // index de la loi de puissance des sources  
 static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index  
 
Index: trunk/Cosmo/RadioBeam/lobe.cc
===================================================================
--- trunk/Cosmo/RadioBeam/lobe.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/lobe.cc	(revision 3973)
@@ -18,4 +18,35 @@
 }
 
+/* --Methode-- */
+void BeamEffect::ApplyLobe(TArray< TF >& a, double dx, double dy, double f0)
+// dx, dy en radioans, f0, df en MHz
+{
+  Timer tm("BeamEffect::ApplyLobe");
+  FFTWServer ffts(true);                     
+  ffts.setNormalize(true); 
+  
+  H21Conversions conv;
+  conv.setFrequency(f0);
+  fresp_.setLambda(conv.getLambda());
+
+  TArray< complex<TF> > fourAmp;
+  double dkx = DeuxPI/(double)a.SizeX()/dx;
+  double dky = DeuxPI/(double)a.SizeY()/dy;
+  
+  for(sa_size_t kz=0; kz<a.SizeZ(); kz++) {
+    TArray< TF > slice( a(Range::all(), Range::all(), kz) );
+    ffts.FFTForward(slice, fourAmp);
+    ApplyLobeK2D(fresp_, fourAmp, dkx, dky);
+    ffts.FFTBackward(fourAmp, slice, true);
+    if (kz%20==0)  cout << "BeamEffect::ApplyLobe() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
+  }
+  double mean, sigma;
+  TF min, max;
+  a.MinMax(min, max);
+  MeanSigma(a, mean, sigma);
+  cout << " BeamEffect::ApplyLobe() - Result Min=" << min << " Max=" << max 
+       << " Mean=" << mean << " Sigma=" << sigma << endl; 
+  return;
+}
 
 
Index: trunk/Cosmo/RadioBeam/lobe.h
===================================================================
--- trunk/Cosmo/RadioBeam/lobe.h	(revision 3949)
+++ trunk/Cosmo/RadioBeam/lobe.h	(revision 3973)
@@ -27,5 +27,10 @@
   // Definition de l'objet avec la reponse en frequence de l'instrument 
   BeamEffect(Four2DResponse& resp, bool preservefreq0=true);
-  // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement  
+
+  // Applique l'effet d'un lobe fixe (frequency independent) au cube 3D (2 angles, frequence)
+  void ApplyLobe(TArray< TF >& a, double dx, double dy, double freq);  
+
+  // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement 
+  // En faisant varier la reponse de lobe avec la frequence 
   void ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df);  
 
Index: trunk/Cosmo/RadioBeam/mdish.cc
===================================================================
--- trunk/Cosmo/RadioBeam/mdish.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/mdish.cc	(revision 3973)
@@ -24,7 +24,7 @@
   switch (typ_) 
     {
-    case 1:   // Reponse gaussienne parabole diametre D exp[ - 0.5 (lambda  k_g / D )^2 ]
+    case 1:   // Reponse gaussienne parabole diametre D exp[ -2 (lambda  k_g / D )^2 ]
       wk = sqrt(kx*kx+ky*ky)/dx_;
-      wk = 0.5*wk*wk;
+      wk = 2*wk*wk;
       return exp(-wk);
       break;
Index: trunk/Cosmo/RadioBeam/mdish.h
===================================================================
--- trunk/Cosmo/RadioBeam/mdish.h	(revision 3949)
+++ trunk/Cosmo/RadioBeam/mdish.h	(revision 3973)
@@ -21,5 +21,5 @@
 
 // -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (frequences angulaires theta,phi)
-// typ=1 : Reponse gaussienne parabole diametre D exp[ - 0.5  (lambda  k_g / D )^2 ]
+// typ=1 : Reponse gaussienne parabole diametre D exp[ - 2  (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) 
Index: trunk/Cosmo/RadioBeam/plpkn.pic
===================================================================
--- trunk/Cosmo/RadioBeam/plpkn.pic	(revision 3949)
+++ trunk/Cosmo/RadioBeam/plpkn.pic	(revision 3973)
@@ -2,26 +2,79 @@
 ### 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
+###    Janvier 2011 ,  BAORadio/Reza
 ########################################################################
 
-if ( $# < 1 ) then 
-  echo ' Usage: exec plpkn PPFName_pknoise'
-  return 
-endif
-
-echo "---> openppf $1 "
-openppf $1
+echo ' -----> plpknew.pic : opening PPF files ... "
+delobjs *
+# openppf ../../PkNoise/cmvhpkz.ppf 
+# openppf ../pkz0p25.ppf
+openppf ../pkz0p7.ppf
+rename hpkz hpkz0p7
+openppf ../pkz1p0.ppf
+rename hpkz hpkz1p0
+
+defscript pknopen 
+  set f $1
+  set nm $2
+  echo '------ pknopen File= ' $f ' name=' $nm 
+  openppf $f
+  rename h1dnoise h1dn$nm
+  rename h1drep h1dr$nm
+  rename fracmodok fracok$nm
+  mv dtnoise pkn$nm
+  rename h2drep h2dr$nm
+endscript 
+
+pknopen pknf11x11.ppf A121d 
+pknopen pknnan128.ppf B128d 
+pknopen pknconfC.ppf  C129d 
+pknopen pknD50m.ppf   D50m 
+pknopen pknD75m.ppf   D75m 
+pknopen pknD100m.ppf   D100m 
+pknopen pknD200m.ppf   D200m 
+pknopen pknf20x20.ppf E400d
+pknopen pknf4cyl.ppf F4cyl
+pknopen pknf4cylp.ppf F4cylp
+pknopen pknf8cyl.ppf G8cyl
+pknopen pknf8cylp.ppf G8cylp
+
+
+listobjs 
+
+#  A z = 0.7
+z = 0.7
+c = 3.e5
+H = 102
+Da = 2488 
+nu21 = 1.42e9
 
 echo '----> Executing anapkn.pic'
 exec anapkn.pic
 setup5 
-scalewz 0.7 2500 100 1
+# scalewz 0.7 2488 102 1
 # scalewz 0.25 989.13 80.26 1
-y1 = 500*$cct21
-y2 = 9e4*$cct21
-xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
+scalewz 1 3300 120.5 1
+Da = 3300 
+
+
+y1 = 400*$cct21
+y2 = 7e4*$cct21
+xyl = "xylimits=0.005,0.5,$y1,$y2 logx logy minorticks"
+# H = 71.9 km/s/Mpc = 1.0271 x 70 km/s
+h70 = 1.0271
+h70cube = $h70*$h70*$h70
+set kk pow(10.,x)/$h70
+
 defscript plpklss
-  n/plot hpkz.val*${cct21}%$kk ! ! "same notit nsta connectpoints black $xyl line=solid,2"
-  addtext 0.03 2000 '*** P(k)-LSS ***' 'font=helvetica,bolditalic,16 black'
+  scalewz 1 3300 120.5 1
+#  scalewz 0.7 2488 102 1
+  convpk2t21
+  cct21h70 = $cct21*$h70cube
+  n/plot hpkz1p0.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "
+#  n/plot hpkz0p7.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "
+
+  addtext 0.01 4000 '[ P(k)-LSS z=1.0 ]' 'font=helvetica,bolditalic,16 black'
+  setaxelabels 'k_comov  (h70 Mpc^-1)'  'P21(k) mK^2 x (Mpc/h70)^3' 'font=helvetica,bolditalic,16'
+
 endscript  
 
@@ -33,14 +86,8 @@
 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 
+Lambda = 0.21*(1+$z)
 Lam2 = $Lambda*$Lambda 
 
@@ -52,133 +99,177 @@
   Dyol = $Dy/$Lambda 
 
-  FOV = (1.2*1.2*$Lam2/$Dx/$Dy) 
+#   FOV = (1.2*1.2*$Lam2/$Dx/$Dy) 
+  FOV = ($Lam2/$Dx/$Dy) 
+
   FOVDEG = $FOV*$R2D2 
   NPointing = 10000/$FOVDEG
   tinteg = 365*86400/$NPointing
-  PNOISE = $Tsys*$Tsys/$tinteg/$Dxol/$Dyol
+  PNOISE = 2.*$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
+  PNOISE = $PNOISE*1.05
   echo " FOV = $FOV deg^2  NPointing= $NPointing"
   echo " tinteg= $tinteg sec  PNOISE= $PNOISE mK^2"
 endscript 
 
-defscript plnoisedish 
+
+defscript plnoiseDishes
+  Dx = 75
+  Dy = 75
+  calcul 
+  plot2d pknD75m  k/$Da pnoise*$PNOISE/5 fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
   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 navyblue grid cpts nsta notit' 
+  plot2d pknD100m  k/$Da pnoise*$PNOISE/20 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
+
+  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'
+
+
+
+  set lines ( '(a) 75m Dish, 5 beams'  '(b) 100m Dish, 10 beams'  '(b) 100m Dish, 100 beams'  )
+  set cols ( red blueviolet blue )
+  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
+  settitle ' PNoise(k) : 5/10/100 beams/polar @z=1' ' ' 'font=helvetica,bold,16' 
+endscript
+
+defscript plnoiseEFGH
+  Dx = 12*0.9
+  Dy = 0.844248*0.8
+  calcul 
+  plot2d pknF4cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts forestgreen nsta notit' 
+  plot2d pknF4cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts green nsta notit' 
+
+  plot2d pknG8cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violetred nsta notit' 
+  plot2d pknG8cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violet nsta notit' 
+
+  Dx = 5*0.9
+  Dy = 5*0.9
+  calcul 
+  plot2d pknE400d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts turquoise nsta notit'
+
+  set lines (  '(e) 20x20:400xD=5m' '(f) 4Cyl-12mx85m, 400 rec/pol' '(fp) 4Cylp-12mx85m, 400 rec/pol'  '(g) 8Cyl-12mx105m, 960 rec/pol' '(gp) 8Cylp-12mx105m, 960 rec/pol' )
+  set cols ( turquoise forestgreen green violetred violet  )
+  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
+  settitle ' PNoise(k) : Dishes/Cylinders, 400/400/960 recv/pol @z=1' ' ' 'font=helvetica,bold,16' 
+endscript 
+
+defscript plnoiseABCD
+  Dx = 5*0.9
+  Dy = 5*0.9
+  calcul 
+  plot2d pknA121d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts magenta nsta notit'
+  plot2d pknB128d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
+  plot2d pknC129d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts orange nsta notit'
+
+  Dx = 75
+  Dy = 75
+  calcul 
+  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
+
+  set lines (  '(a) 11x11:121xD=5m' '(b) 128xD=5m' '(c) 129xD=5m' '(d) 75m Dish, 100 beams' )
+  set cols ( magenta red orange blueviolet )
+  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
+  settitle ' PNoise(k) : 121,128,129,100 receiver/polar @z=1' ' ' 'font=helvetica,bold,16' 
+endscript
+
+defscript plfracABCD
+  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
+  plot2d pknA121d k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts magenta nsta notit'
+  plot2d pknB128d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts red nsta notit'
+  plot2d pknC129d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts orange nsta notit'
+  plot2d pknD75m  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts violetblue nsta notit'
+endscript
+
+defscript plfracEFGH
+  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
+  plot2d pknF4cyl k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts forestgreen nsta notit'
+  plot2d pknF4cylp k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts green nsta notit'
+  plot2d pknG8cyl k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts gold nsta notit'
+  plot2d pknG8cylp  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts yellow nsta notit'
+  plot2d pknE400d  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts magenta nsta notit'
+
+endscript
+
+defscript pldishes 
+
+  Dx = 50
+  Dy = 50
+  calcul 
+  plot2d pknD50m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts skyblue nsta notit'
+
+  Dx = 75
+  Dy = 75
+  calcul 
+  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'
+
+  Dx = 100
+  Dy = 100
+  calcul 
+  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts navyblue nsta notit'
+
   Dx = 200
   Dy = 200
   calcul 
-  plot2d noiseD2 x/$Da val*$PNOISE nb>10 'same line=solid,2  cpts grey nsta notit' 
-  set lines ( '100mDish' '200mDish'  )
-  set cols ( navyblue grey  )
-  textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.1,0.3,0.7,0.8' 
-  setaxelabels 'k (Mpc^-1) ' 'PNoise(k) mk^2 Mpc^3' 'font=helvetica,bolditalic,16' 
-endscript 
-
-defscript plnoiseA
-  Dx = 5*0.95
-  Dy = 5*0.95
-  calcul 
-  plot2d noisemdf64 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
-#  plot2d noisemds x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
-  plot2d noisemdsB x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts siennared nsta notit'
-  plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violetred nsta notit'
-  Dx = 25*0.3
-  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 (  'FilledA:64x5mD' 'SparseB:72x5mD' 'SparseC:129x5mD' 'Pitts2Cyl(=64C)' 'PerfPitts2Cyl(=64C)' )
-  set cols ( red siennared magenta  forestgreen green )
+  plot2d pknD200m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
+
+  set lines (  'D:Dish50m' 'D:Dish75m' 'D:Dish100m' 'D:Dish200m' )
+  set cols ( skyblue blue navyblue blueviolet )
   textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
-  settitle ' PNoise(k) : Dishes/Cylinders, 64/72/129 channels' ' ' 'font=helvetica,bold,16' 
-endscript 
-
-defscript plnoiseB
-  Dx = 5*0.95
-  Dy = 5*0.95
-  calcul 
-  plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
-  plot2d noisemdf x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts gold nsta notit' 
-  Dx = 10*0.95
-  Dy = 0.5*0.9
-  calcul 
-  plot2d noise3cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit' 
-  plot2d noise3cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts magenta nsta notit' 
-
-  Dx = 12*0.95
-  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 royalblue nsta notit' 
-
-  set lines ( 'SparseC:129x5mD' 'Filled400x5m' '3xCyl10x64m(=384C)' 'Perf3xCyl10x64m(=384C)' )
-  set lines ( $lines '8Cyl12x96m(=1536C)' 'Perf8Cyl12x96m(=1536C)' )
-  set cols ( red gold  violet magenta  blue royalblue )
-  textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.25,0.55,0.6,0.95' 
-  settitle ' PNoise(k) : Dishes/Cylinders, 129/384/400/1536 channels' ' ' 'font=helvetica,bold,16'	
-endscript
-
-defscript xxx
-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' 
-
-
-disp noiseD 'logy nsta' 
-disp noiseD2 'same grey nsta'
-disp noisemdf 'same cyan nsta'
-disp noisemdf64 'same brown nsta'
-disp noisemds 'same red nsta'
-disp noisemdsC 'same orange nsta'
-
-disp noisemdsB 'same yellow nsta'
-disp noisemdsfp 'same yellow nsta'
-# disp noisemdsd7 'same gold nsta'
-disp noisefcyl 'same blue nsta'
-disp noisefcylP 'same skyblue nsta'
-disp noise3cylP 'same magenta nsta'
-disp noise3cyl 'same violet nsta'
-disp noise2cyl 'same forestgreen nsta'
-disp noise2cylP 'same green nsta'
-
-endscript 
-
-defscript AA
-  plnoisedish
-  plpklss
-  plnoiseA
+  settitle ' PNoise(k) : Dish D=50m,75m,100m,200m' ' ' 'font=helvetica,bold,16' 
+endscript
+
+
+
+defscript plnoisenancay
+#  Dx = 5*0.95
+#  Dy = 5*0.95
+#  calcul 
+#  plot2d pknnan24 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
+#  plot2d pknnan25 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts orange nsta notit'
+#  Dx = 3.5*0.95
+#  Dy = 3.5*0.95
+#  calcul 
+#  plot2d pknnan36 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
+#  plot2d pknnan40 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit'
+#  Dx = 7*0.9
+#  Dy = 2*$Lambda*0.9
+#  calcul 
+#  plot2d pknpit2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
+  Dx = 9*0.9
+  Dy = 2*$Lambda*0.9
+  calcul 
+  plot2d pknpit2cylw x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
+  set lines ( '2 Cylinders 20mx9m' )
+  set cols ( red )
+  textdrawer lines cols 'noframe font=helvetica,bold,16 inset=0.4,0.4,0.4,0.4' 
+
+endscript
+
+defscript ABCD
+#  plnoisedish
+#  zone 1 2
+#  newwin 1 1 1000 600
+  plpklss
+  plnoiseABCD
+#  plfracABCD
+endscript 
+
+defscript EFGH
+#  plnoisedish
+#  zone 1 2
+#  newwin 1 1 900 600
+  plpklss
+  plnoiseEFGH
+#  plfracEFGH
+endscript 
+
+defscript DDDD
+#  plnoisedish
+#  zone 1 2
+  plpklss
+  plnoiseDishes
+#  plfracABCD
 endscript 
 
@@ -189,29 +280,62 @@
 endscript 
 
+defscript CC
+  plnoisedish
+  plpklss
+  plnoiseC2
+endscript 
+
+defscript nancay
+  plnoisedish
+  plpklss
+  plnoisenancay
+endscript 
+
 defscript POSCOV
-  disp posspB red
+  openppf hdt_repnan128.ppf
+  rename mdish mdB128d
+  rename h2rep uvB128d
+
+  openppf hdt_repf11x11.ppf
+  rename mdish mdA121d
+  rename h2rep uvA121d
+
+  openppf hdt_repconfC.ppf
+  rename mdish mdC129d
+  rename h2rep uvC129d
+
+  newwin 2 1 800 400
+  zone 2  1
+  nt2d mdB128d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red '
   setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16' 
-  settitle ' Config B dish positions - 72 dishes  ' ' ' 'font=helvetica,bold,16' 
-#  w2ps
-   
-  disp posspC red
+  settitle '(b) 128 D=5m dishes in 8 rows ' ' ' 'font=helvetica,bold,16' 
+  nt2d mdC129d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red'
   setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16' 
-  settitle ' Config C dish positions - 129 dishes  ' ' ' 'font=helvetica,bold,16' 
-#  w2ps
-  disp mfill64 'h2disp=img colbr128 notit'
-  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
-  settitle ' u-v coverage , Filled 8x8 - 64 x 5m Dishes' ' ' 'font=helvetica,bold,16' 
-#  w2ps
-  disp m3cyl 'h2disp=img colbr128 notit'
-  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
-  settitle ' u-v coverage , 3 Cylinders 10mx64m ' ' ' 'font=helvetica,bold,16' 
-#  w2ps
-  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 , Filled 20x20 - 400 x 5m Dishes' ' ' 'font=helvetica,bold,16' 
-#  w2ps
-  disp msparsC 'h2disp=img colbr128 notit'
-  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16' 
-  settitle 'u-v coverage, Sparse-C: 129x5mD Over 80mx80m (Rot~Pi/4)' ' ' 'font=helvetica,bold,16' 
+  settitle '(c) 129 D=5m dishes ' ' ' 'font=helvetica,bold,16' 
+  pssetfilename configab.ps
+  w2ps
+  w2eps configab.eps
+  psclosefile 
+
+  newwin 2 2 800 800 
+  disp uvA121d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
+  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
+  settitle 'u-v coverage, (a) 11x11 D=5m dishes Over 55mx55m' ' ' 'font=helvetica,bold,12' 
+  disp uvB128d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
+  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
+  settitle 'u-v coverage, (b) 8 row of 16xD=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12' 
+  disp uvC129d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
+  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
+  settitle 'u-v coverage, (c) 129 D=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12' 
+  disp h2drD75m 'h2disp=img colbr128 h2dyn=0.02,0.8 notit nsta'
+  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
+  settitle 'u-v coverage, (d) D=75 m dish with 100 beams' ' ' 'font=helvetica,bold,12' 
+
+  pssetfilename uvcovabcd.ps
+  w2ps
+  w2eps uvcovabcd.eps
+  psclosefile 
+
+#  imag uvB128d 'colbr128 stdaxes showcmap=right' 
   
 endscript  
Index: trunk/Cosmo/RadioBeam/srcat2cube.cc
===================================================================
--- trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3949)
+++ trunk/Cosmo/RadioBeam/srcat2cube.cc	(revision 3973)
@@ -278,5 +278,5 @@
       lowoksrccnt++;
     }
-    else slo=5.;
+    else slo=-1.;
     for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
       double rapfreq = pow((freq0+k*dfreq)/infreq, slo);
Index: trunk/Cosmo/RadioBeam/subtractradsrc.cmd
===================================================================
--- trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3949)
+++ trunk/Cosmo/RadioBeam/subtractradsrc.cmd	(revision 3973)
@@ -12,4 +12,7 @@
 csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,3 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
 
+# 1.c/ To run SimLSS with GSM map parametersand 40x40 deg maps (DeltaFreq=500 MHz)
+csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 600,1.9 -y 800,1.9 -z 256,2.8 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
+
 # 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention 
 #  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray) 
@@ -18,5 +21,5 @@
 
 csh> cat > racube.pic
-set f lssz056
+set f lssz060
 readfits ${f}_r.fits
 rename ${f}_r map
@@ -33,21 +36,26 @@
 print lsscube
 # expmeansig lsscube val
-saveppf lsscube lsscube.ppf 
+saveppf lsscube lsscubez060.ppf 
 
 csh> spiapp -term -exec racube.pic 
 
+#### Cube LSS 40x30 deg (3')  @ z=0.6 ( lsscubez060.ppf ) 
+#### -> Size= 122880000 Mean=-7.01664e-05 Sigma=2.53016 Min=-13.7439 Max=14.4648
 
 ## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
 # 2.a/ Synchrotron map from HASLAM 400 MHz map
-csh> ./Objs/syncube syncmap_eq.fits syncube.ppf
+csh> ./Objs/syncube syncmap_eq.fits syncube.ppf syncmap.ppf
 # 2.b/ radio source cube from NVSS catalog
-csh> ./Objs/srcat2cube nvss.fits nvsscube.ppf
+csh> ./Objs/srcat2cube -nvss nvss.fits nvsscube.ppf nvssmap.ppf
+# Or from the north20 catalog :
+csh> ./Objs/srcat2cube -north20 north20cm.fits north20cube.ppf north20map.ppf
+
 # 2.c/ Add the two cubes using the following spiapp script 
 csh> cat > sumcubes.pic
 openppf syncube.ppf
-openppf srcnor3d.ppf
+openppf radsrccube.ppf
 # expmeansig syncube val 
-# expmeansig srcnor3d val
-c++exec TArray<r_4> fgndcube = syncube+srcnor3d; KeepObj(fgndcube);
+# expmeansig nvsscube val
+c++exec TArray<r_4> fgndcube = syncube+radsrccube; KeepObj(fgndcube);
 print fgndcube
 # expmeansig fgndcube val
@@ -56,42 +64,66 @@
 csh> spiapp -term -exec sumcubes.pic 
 
+#### syncube:Mean= 1.8101  Sigma= 0.326538  Min= 0.857019 Max= 3.58987
+#### nvsscube: Mean= 1.95073  Sigma= 1.68515  Min= 0.857019 Max= 428.398
+#### fgndcube=syncube+nvsscube: Mean= 0.140623  Sigma= 1.65068  Min= 0 Max= 426.559
+#### north20: fgndcube_north:
+
+## Step 2.b/ Produce foreground cube from GSM 
+csh> ./Objs/gsm2cube ../Catalogs/GSM/ 1 256 fgndcube_gsm.ppf
+
 ## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
-csh> ./Objs/applobe 50. fgndcube.ppf fgndcube_lobe.ppf
-csh> ./Objs/applobe 50. lsscube.ppf lsscube_lobe.ppf
-## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 130  
-csh> ./Objs/applobe 50. lsscube_lobe.ppf lsscube_corlobe.ppf 130
-csh> ./Objs/applobe 50. fgndcube_lobe.ppf fgndcube_corlobe.ppf 130
+csh> set ddish=55.
+csh> set ddishcor=50.
+csh> ./Objs/applobe $ddish fgndcube.ppf fgndcube_lobe.ppf
+csh> ./Objs/applobe -fib $ddish fgndcube.ppf fgndcube_flobe.ppf
+csh> ./Objs/applobe $ddish lsscube.ppf lsscube_lobe.ppf
+csh> ./Objs/applobe -fib $ddish lsscube.ppf lsscube_flobe.ppf
+## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150  (55 m @ z=0.7 - 820 MHz)  
+csh> ./Objs/applobe $ddish lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
+csh> ./Objs/applobe $ddish fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor 
+
+## Step 3.c/ Apply lobe (Filled 11x11 5m dishes array) effect on foreground cube and LSS cube
+csh> ./Objs/applobe repf11x11.ppf fgndcube.ppf fgndcube_lobe.ppf
+csh> ./Objs/applobe repf11x11.ppf lsscube.ppf lsscube_lobe.ppf
+## Step 3.d/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150 (55 m @ z=0.7 - 820 MHz) 
+csh> ./Objs/applobe repf11x11.ppf lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
+csh> ./Objs/applobe repf11x11.ppf fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor
 
 ### Step 4/ Compute power spectra 
-## mass to temperature converion factor   CT21 ~= 0.2 mK 
+## mass to temperature converion factor   CT21 ~= 0.21 mK for gHI=2% , 0.11 for gHI=1% , 0.13 for gHI=0.008x(1+0.6)
 ## Foreground maps are in temperature 
 ## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
-## Tsys ~ 50 K , DeltaFreq ~ 0.275 MHz , t_obs ~ 1 day ~ 80 000 s.
-## sigma_noise ~ 0.35 mK 
+## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
+## sigma_noise ~ 0.25 mK -> 3 mK  
 # 4.a/ LSS power spectrum  without noise
-csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.2 
+csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.13
 # and with noise  
-csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.35 
+csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.13  3
 #  with the lobe effect 
-csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2  
-csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.2 0.35  
-csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2  
+csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.13
+csh> ./Objs/calcpk lsscube_flobe.ppf lsspkflobe.ppf 0.13
+csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.13 3  
+csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.13  
 
 # 4.b/ Foreground power spectrum 
 csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
 csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
+csh> ./Objs/calcpk fgndcube_flobe.ppf fgndpkflobe.ppf 1000
 csh> ./Objs/calcpk fgndcube_corlobe.ppf fgndpkcorlobe.ppf 1000
 
 # 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
-csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0. P2
+csh> set beamdesc=repf11x11.ppf
+csh> set ddishcor=50.
+csh> set noiselev=1.
+csh> ./Objs/calcpk2 lsscube.ppf 0.13 fgndcube.ppf 1000 subpk.ppf $noiselev $beamdesc 0. 0. P2
 # 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction  
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. P2 
-# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130 
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3. P2 
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpklobe.ppf $noiselev $beamdesc 0. 0. P2 
+# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam= $ddishcor
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2 
 #  Or using a linear fit for foreground subtraction (old version)
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf 0.35 50. 130. 3. P1
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf $noiselev $beamdesc $ddishcor 0. P2 
 # 4.f / Estimate residual noise from Foreground removal : 
-csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3. P2
-csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.
+csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residcorlobe.ppf  $noiselev $beamdesc $ddishcor 0. P2
+csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residnocor.ppf $noiselev $beamdesc 0. 0. P2
 
 ### Step 5 / Check the results using spiapp 
@@ -100,7 +132,9 @@
 openppf fgndpk.ppf 
 openppf fgndpklobe.ppf 
+openppf fgndpkflobe.ppf 
 openppf fgndpkcorlobe.ppf 
 openppf lsspk.ppf 
 openppf lsspklobe.ppf 
+openppf lsspkflobe.ppf 
 openppf lsspklobewn.ppf 
 openppf subpklobe.ppf 
@@ -136,5 +170,5 @@
 
 #  Calcul du volume total en Mpc^3
-set VOL 3*3*3*360*360*256
+set VOL (1.9*1.9*2.8*800*600*256)
 plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
 plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared' 
