Changeset 3789 in Sophya for trunk/Cosmo

Timestamp:

Jun 27, 2010, 4:57:36 PM (15 years ago)

Author:

ansari

Message:

Corrections diverses, Reza 27/06/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• applobe.cc (modified) (4 diffs)
• calcpk.cc (modified) (3 diffs)
• calcpk2.cc (modified) (4 diffs)
• cubedef.h (modified) (3 diffs)
• fgndsub.cc (modified) (5 diffs)
• fgndsub.h (modified) (1 diff)
• lobe.cc (modified) (4 diffs)
• makefile (modified) (1 diff)
• mdish.cc (modified) (3 diffs)
• mdish.h (modified) (2 diffs)
• subtractradsrc.cmd (modified) (6 diffs)

trunk/Cosmo/RadioBeam/applobe.cc

@@ -3 +3 @@
     R. Ansari , C. Magneville - Juin 2010 
 
-  Usage: applobe In3DPPFName Out3DPPFName
+  Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamArcmin]
 ---------------------------------------------------------------  */
 
@@ -50 +50 @@
 
   if (narg < 3) {
-    cout << "Usage: applobe In3DPPFName Out3DPPFName \n" << endl;
+    cout << "Usage: applobe In3DPPFName Out3DPPFName [ResmapleFactor=0.5,0.333...] [TargetBeamArcmin]\n" << endl;
     return 1;
   }
@@ -57 +57 @@
   string outname = arg[2];
   string inname = arg[1];
+  bool fgresample=false;
+  double fresamp=1.;
+  if (narg>3) {
+    fresamp=atof(arg[3]);
+    if ((fabs(fresamp)>1.e-2)&&(fabs(fresamp-1.)>1.e-2))  fgresample=true; 
+  }
+  bool fgcorrbeam=false;
+  double tbeamarcmin=15.;
+  if (narg>4) {
+    tbeamarcmin=atof(arg[4]);
+    if (tbeamarcmin>1.e-3)  fgcorrbeam=true;
+  }
+
   int rc = 91;
 
@@ -88 +101 @@
     conv.setFrequency(Freq0MHz);
     double lambda = conv.getLambda();
-    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
+    Four2DResponse fresp(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda, lambda);    
     BeamEffect beam(fresp);
 
-    cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
-    beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
+    if (fgcorrbeam) {
+      double DoL = 1.22/ArcminToRadian(tbeamarcmin);
+      Four2DResponse tbeam(2, DoL, DoL );
+      cout << "applobe[3]: calling Correct2RefLobe() with target beam:" << tbeamarcmin 
+           << " arcmin -> D/Lambda=" << DoL << endl;
+      beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq);
+    }
+    else {
+      cout << "applobe[3]: calling ApplyLobe3D() ... " << endl;
+      beam.ApplyLobe3D(incube, dx, dy, Freq0MHz, dfreq);
+    }
+    TArray< r_4 > outcube;
+    if (fgresample) {
+      cout << "applobe[4]: calling ReSample(incube," << fresamp << "," << ",1.) ... " << endl;
+      outcube.Share(beam.ReSample(incube, fresamp, fresamp, 1.));
+    }
+    else outcube.Share(incube);
 
-    cout << "applobe[4]: calling ReSample(incube, 0.5, 0.5, 1.) ... " << endl;
-    TArray< r_4 > outcube = beam.ReSample(incube, 0.5, 0.5, 1.);
-    
     outcube.Show(); 
     MeanSigma(outcube, mean, sigma);

trunk/Cosmo/RadioBeam/calcpk.cc

@@ -20 +20 @@
 #include "qhist.h" 
 #include "lobe.h" 
+#include "cubedef.h" 
 
 #include "histinit.h"
@@ -86 +87 @@
     RandomGenerator rg;
     Four3DPk pkc(four3d, rg);
-    
+    double dkxmpc = DeuxPI/(double)inmap.SizeX()/XCellSizeMpc;
+    double dkympc = DeuxPI/(double)inmap.SizeY()/YCellSizeMpc;
+    double dkzmpc = DeuxPI/(double)inmap.SizeZ()/ZCellSizeMpc;
+    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
+
     HProf hp = pkc.ComputePk(0.,256);
     {
@@ -92 +97 @@
       POutPersist po(outname);
       po << hp;
+      outname = "f3d_" + outname;
+      cout << "calcpk[4.b] : writing four3d to PPF file  " << outname << endl;
+      POutPersist pof(outname);
+      pof << four3d;
     }
     tm.Split(" Done CalcP(k) ");    

trunk/Cosmo/RadioBeam/calcpk2.cc

@@ -61 +61 @@
     if (narg>6) {
       pixsignoise=atof(arg[6]);
-      fgaddnoise=true;
+      if (pixsignoise>1.e-6)  fgaddnoise=true;
     }
     
     bool fgcorrbeam=true;
-    double tbeamarcmin=30.;
+    double tbeamarcmin=15.;
     if (narg>7) {
       tbeamarcmin=atof(arg[7]);
@@ -119 +119 @@
     conv.setFrequency(Freq0MHz);
     double lambda = conv.getLambda();
-    Four2DResponse arep(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda);
+    Four2DResponse arep(2, InterfArrayDiametre/lambda, InterfArrayDiametre/lambda, lambda);
     double DoL = 1.22/ArcminToRadian(tbeamarcmin);
     Four2DResponse tbeam(2, DoL, DoL );
@@ -128 +128 @@
       cleaner.BeamCorrections();
     }
+    cout << " calcpk2[3.b] : calling cleaner.CleanNegatives() ... " << endl;
+    cleaner.CleanNegatives();
     if (fgclnsrc) {
-      cout << "calcpk2[3.b] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
+      cout << "calcpk2[3.c] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
       cleaner.CleanPointSources(nsigsrc);
     }
@@ -150 +152 @@
     RandomGenerator rg;
     Four3DPk pkc(four3d, rg);
+    double dkxmpc = DeuxPI/(double)exlss.SizeX()/XCellSizeMpc;
+    double dkympc = DeuxPI/(double)exlss.SizeY()/YCellSizeMpc;
+    double dkzmpc = DeuxPI/(double)exlss.SizeZ()/ZCellSizeMpc;
+    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
     
     HProf hp = pkc.ComputePk(0.,256);

trunk/Cosmo/RadioBeam/cubedef.h

@@ -9 +9 @@
 static sa_size_t NTheta=360;
 static sa_size_t NPhi=360;
-static sa_size_t NFreq = 200;
+static sa_size_t NFreq = 256;
 
 static double Theta0Degre = 75.;  // -> Delta = +30 deg
@@ -16 +16 @@
 static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
 
-static double Freq0MHz = 840.;       // premiere frequence 
-static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
+static double Freq0MHz = 875.;       // premiere frequence 
+static double FreqSizeMHz = 70.4;    // Taille de la carte en MHz (axe Z)
 
+/* z = 0.56, freq_center=910 MHz , dA_comov ~= 2060 Mpc ,  
+       5 arcmin -> 3 Mpc , 0.275 MHz -> 3 Mpc  (H(ze) ~ 94 km/s/Mpc */
+// Taille de cellule pour le calcul du spectre de puissance 3D 
+static double ComovDA = 2060.;
+static double XCellSizeMpc = 3.;
+static double YCellSizeMpc = 3.;
+static double ZCellSizeMpc = 1.5;
 
 //---  Parametres des lois de puissance en frequence 
@@ -33 +40 @@
 
 // Parametres pour la reponse en Fourier de l'instrument :
-static double InterfArrayDiametre = 64.;
+static double InterfArrayDiametre = 50.;  // Taille du reseau en metres 
 
 /*

trunk/Cosmo/RadioBeam/fgndsub.cc

@@ -33 +33 @@
   beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_);
   cout << " ForegroundCleaner::BeamCorrections() done " << endl;
+}
+
+/* --Methode-- */
+int ForegroundCleaner::CleanNegatives(TF seuil)
+{
+  sa_size_t nneg = 0.;
+  for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++) 
+    for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++) 
+      for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) 
+        if (skycube_(kx, ky, kz) < seuil)  {
+          nneg++; skycube_(kx, ky, kz)=seuil;
+        }
+  cout << " ForegroundCleaner::CleanNegatives " << nneg << " sky-pixels <" << seuil << " changed to" << seuil << endl;
+  return (int)nneg;
 }
 
@@ -82 +96 @@
 // Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
 // Return_Array : foreground subtracted LSS signal 
-{
   sa_size_t sz[5];   sz[0]=skycube_.SizeX();  sz[1]=skycube_.SizeY();
   synctemp.SetSize(2, sz); 
@@ -92 +105 @@
   // double freq0 : Frequence premier index en k (MHz)
   // double dfreq :   // largeur en frequence de chaque plan (Mhz)  
+  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  
+    vlnf(k)=log((double)k*dfreq_+freq0_);
+
+  sa_size_t nbinfini=0;
+  sa_size_t nbbad=0;
+  sa_size_t imodprt=skycube_.SizeX()/6;
+  sa_size_t jmodprt=skycube_.SizeY()/6;
   for(sa_size_t i=0; i<skycube_.SizeX(); i++) 
     for(sa_size_t j=0; j<skycube_.SizeY(); j++)  {
@@ -97 +117 @@
       s1=sx=sx2=sy=sxy=0.;
       for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
-        double lnf=log((double)k*dfreq_+freq0_);
-        vlnf(k)=lnf;
+        double lnf = vlnf(k);
         double ttot=(r_8)(skycube_(i,j,k));
         if (ttot < 1.e-5) continue;
@@ -108 +127 @@
       double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
       double T0 = exp(beta+alpha*vlnf(0));
-      if ((i%16==0)&&(j%27==0)) 
+
+      bool fgnan = false;
+      if (!isfinite(alpha)||(!isfinite(beta))) {
+        cout << "extractLSSCube[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        alpha=beta=-999.;
+        fgnan = true;  nbinfini++;
+      }
+      else {
+        double axp1 = beta+alpha*vlnf(0);
+        double axp2 = beta+alpha*vlnf(vlnf.Size()-1);
+        
+        if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
+        cout << "extractLSSCube[" << i << "," << j << "] BAD alpha=" << alpha 
+             << " beta=" << beta << " T0=" << T0 << " - (mapsync="
+             << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
+        fgnan = true;  nbbad++;
+        }
+      }
+      if ((i%imodprt==0)&&(j%jmodprt==0)) 
         cout << "extractLSSCube[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
-             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,125) << ")" << endl; 
+             << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl; 
       synctemp(i,j) = T0;
       specidx(i,j) = alpha;
-      for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
-        r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
-        omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+      if (fgnan) {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++)   omap(i,j,k) = 0.;
+      }
+      else {
+        for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
+          r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
+          omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
+        }
       }
     }
+  cout << " ForegroundCleaner::extractLSSCube() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
   return omap;
 }
 
-}

trunk/Cosmo/RadioBeam/fgndsub.h

@@ -27 +27 @@
   ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube);
   void BeamCorrections();
+  int CleanNegatives(TF seuil=1.e-6);
   int CleanPointSources(double nsigmas=5.);
   TArray< TF > extractLSSCube(TArray< TF >& synctemp, TArray< TF >& specidx);

trunk/Cosmo/RadioBeam/lobe.cc

@@ -6 +6 @@
 
 
-#include "fftwserver.h"    
+#include "fftwserver.h" 
+#include "matharr.h"    
 #include "ctimer.h"    
 
@@ -38 +39 @@
     conv.setFrequency(f0+kz*df);
     fresp_.setLambda(conv.getLambda());
+    //    cout << " DEBUG*" << kz << " lambda=" << conv.getLambda() << " lambda_ratio_=" << fresp_.lambda_ratio_ << endl;
     ApplyLobeK2D(fresp_, fourAmp, dkx, dky);
     ffts.FFTBackward(fourAmp, slice, true);
     if (kz%20==0)  cout << "BeamEffect::ApplyLobe3D() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
   }
+  double mean, sigma;
+  TF min, max;
+  a.MinMax(min, max);
+  MeanSigma(a, mean, sigma);
+  cout << " BeamEffect::ApplyLobe3D() - Result Min=" << min << " Max=" << max 
+       << " Mean=" << mean << " Sigma=" << sigma << endl; 
   return;
 }
@@ -64 +72 @@
     conv.setFrequency(f0+kz*df);
     fresp_.setLambda(conv.getLambda());
-    Four2DRespRatio rratio(fresp_, rep);
+    Four2DRespRatio rratio(rep, fresp_);
     ApplyLobeK2D(rratio, fourAmp, dkx, dky);
     ffts.FFTBackward(fourAmp, slice, true);
     if (kz%20==0)  cout << "BeamEffect::Correct2RefLobe() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl; 
   }
+  double mean, sigma;
+  TF min, max;
+  a.MinMax(min, max);
+  MeanSigma(a, mean, sigma);
+  cout << " BeamEffect::Correct2RefLobe() - Result Min=" << min << " Max=" << max 
+       << " Mean=" << mean << " Sigma=" << sigma << endl; 
   return;
 }
@@ -74 +88 @@
 /* --Methode-- */
 void BeamEffect::ApplyLobeK2D(Four2DResponse& rep, TArray< complex<TF> >& fourAmp, double dkx, double dky)
-//  dx, dy en radians, lambda en metres 
 {
   double kxx, kyy;

trunk/Cosmo/RadioBeam/makefile

@@ -61 +61 @@
         $(CXXLINK) -o Objs/syncube Objs/syncube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
 
-Objs/syncube.o : syncube.cc 
+Objs/syncube.o : syncube.cc  $(PKGHLIST)
         $(CXXCOMPILE) -o Objs/syncube.o syncube.cc
 

trunk/Cosmo/RadioBeam/mdish.cc

@@ -9 +9 @@
 //--------------------------------------------------
 // Constructor
-Four2DResponse::Four2DResponse(int typ, double dx, double dy)
+Four2DResponse::Four2DResponse(int typ, double dx, double dy, double lambda)
   : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
 {
-  setLambdaRef();
-  setLambda();
+  setLambdaRef(lambda);
+  setLambda(lambda);
 }
 
@@ -78 +78 @@
 // -- Four2DRespRatio : rapport de la reponse entre deux objets Four2DResponse
 //---------------------------------------------------------------
-Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b)
-  : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b) 
+Four2DRespRatio::Four2DRespRatio(Four2DResponse& a, Four2DResponse& b, double divzthr)
+  : Four2DResponse(0, a.D(), a.D()), a_(a), b_(b), divzthr_(divzthr)
 {
 }
@@ -87 +87 @@
   double ra = a_.Value(kx,ky);
   double rb = b_.Value(kx,ky);
-  if (rb<1.e-19) rb = 1.e-19;
+  if (rb<divzthr_) {
+    if (ra<rb)  return 0.;
+    else rb=divzthr_;
+  }
   return (ra/rb);
 }

trunk/Cosmo/RadioBeam/mdish.h

@@ -27 +27 @@
 public:
   // On donne dx=D/lambda=Dx/lambda , dy=Dy/lambda 
-  Four2DResponse(int typ, double dx, double dy);
+  Four2DResponse(int typ, double dx, double dy, double lambda=1.);
 
   Four2DResponse(Four2DResponse const& a)
-  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; }
+  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; lambdaref_=a.lambdaref_;
+    lambda_=a.lambda_; lambda_ratio_=a.lambda_ratio_; }
   Four2DResponse& operator=(Four2DResponse const& a)
-  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }
+  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; lambdaref_=a.lambdaref_;
+    lambda_=a.lambda_; lambda_ratio_=a.lambda_ratio_;  return (*this); }
 
   inline void setLambdaRef(double lambda=1.)
   { lambdaref_ = lambda; }
   inline void setLambda(double lambda=1.)
-  { lambda_ = lambda;  lambda_ratio_ = lambdaref_/lambda_; }
+  { lambda_ = lambda;   lambda_ratio_ = lambda_/lambdaref_; }   
   
   // Return the 2D response for wave vector (kx,ky)
@@ -70 +72 @@
 class Four2DRespRatio : public  Four2DResponse {
 public:
-  Four2DRespRatio(Four2DResponse& a, Four2DResponse& b);
+  Four2DRespRatio(Four2DResponse& a, Four2DResponse& b, double divzthr=1.e-1);
   // Return the ratio a.Value(kx,ky) / b.Value(kx, ky) - with protection against divide by zero 
   virtual double Value(double kx, double ky); 
   Four2DResponse& a_;
   Four2DResponse& b_;
+  double divzthr_;
 };
 

trunk/Cosmo/RadioBeam/subtractradsrc.cmd

@@ -8 +8 @@
 ### Step 1/ Produce an LSS data cube  with appropriate size and redshift using SimLSS
 # 1.a/ Run SimLSS
-csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,4 -y 360,4 -z 128,5 -Z 0.6 -8 1. -n 10000 -O 0,2 -o map3dz06B -T 2
+csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,1.5 -Z 0.56 -8 1. -n 10000 -O 0,2 -o lssz056 -T 2
 # 1.b/ 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) 
@@ -15 +15 @@
 
 csh> cat > racube.pic
-set f map3dz06B
+set f lssz056
 readfits ${f}_r.fits
 rename ${f}_r map
@@ -30 +30 @@
 saveppf omap lsscube.ppf 
 
-csh> spiapp -term -exec lsscube.ppf
+csh> spiapp -term -exec racube.pic 
 
 ## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
@@ -41 +41 @@
 openppf syncube.ppf
 openppf nvsscube.ppf
-expmeansig syncube val 
-expmeansig nvsscube val
+# expmeansig syncube val 
+# expmeansig nvsscube val
 c++exec TArray<r_4> fgndcube = syncube+nvsscube; KeepObj(fgndcube);
 print fgndcube
-expmeansig fgndcube val
+# expmeansig fgndcube val
 saveppf fgndcube fgndcube.ppf
 
@@ -53 +53 @@
 csh> ./Objs/applobe fgndcube.ppf fgndcube_lobe.ppf
 csh> ./Objs/applobe lsscube.ppf lsscube_lobe.ppf
+## Step 3.b/ Correct for the lobe effect by bringing all to a beam of 30 arcmin 
+csh> ./Objs/applobe lsscube_lobe.ppf lsscube_corlobe.ppf 1 30
+csh> ./Objs/applobe fgndcube_lobe.ppf fgndcube_corlobe.ppf 1 30
 
 ### Step 4/ Compute power spectra 
@@ -58 +61 @@
 ## Foreground maps are in temperature 
 ## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
-## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
-## sigma_noise ~ 0.25 mK 
+## Tsys ~ 50 K , DeltaFreq ~ 0.275 MHz , t_obs ~ 1 day ~ 80 000 s.
+## sigma_noise ~ 0.35 mK 
 # 4.a/ LSS power spectrum  without noise
 csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.2 
 # and with noise  
-csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.25 
+csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.35 
 #  with the lobe effect 
 csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2  
+csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2  
+
 # 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_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.25
-# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect 
-csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.25
+csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 0. 0.
+# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting for beam to 30 arcmin
+csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 30. 3.
 
 ### Step 5 / Check the results using spiapp