Changeset 3930 in Sophya for trunk/Cosmo

Timestamp:

Dec 23, 2010, 12:49:22 AM (15 years ago)

Author:

ansari

Message:

ajout config double hexagon + programme calcul pknoise a partir du resultat de repicon.cc, Reza 22/12/2010

Location:

trunk/Cosmo/RadioBeam

Files:

• interfconfigs.cc (modified) (2 diffs)
• interfconfigs.h (modified) (2 diffs)
• mdish.cc (modified) (1 diff)
• mdish.h (modified) (2 diffs)
• pknoise.cc (modified) (5 diffs)
• pknoisev1.cc (added)
• repicon.cc (modified) (8 diffs)
• specpk.cc (modified) (5 diffs)
• specpk.h (modified) (2 diffs)

trunk/Cosmo/RadioBeam/interfconfigs.cc

@@ -234 +234 @@
 }
 
- cout << "CreateConfigB/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
+ cout << "CreateConfigD/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
  
  vector<Dish> vd;
@@ -265 +265 @@
   return vd;
 }
+
+/* --Fonction -- */
+vector<Dish> CreateDoubleHexagonConfig(double Ddish, double radius1, double radius2, double Eta)
+{
+  vector<Dish>  vd;
+  int cnt=0;
+  double ang=0.;
+  for(int i=0; i<6; i++) {
+    double x=cos(ang)*radius1;
+    double y=sin(ang)*radius1;
+    cnt++;
+    vd.push_back(Dish(cnt, x, y, Eta*Ddish));
+    ang += M_PI/3.;
+  }
+  ang=M_PI/6.;
+  for(int i=0; i<6; i++) {
+    double x=cos(ang)*radius2;
+    double y=sin(ang)*radius2;
+    cnt++;
+    vd.push_back(Dish(cnt, x, y, Eta*Ddish));
+    ang += M_PI/3.;
+  }
+ cout << ">>>CreateDoubleHexagonConfig(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
+ return vd;
+}

trunk/Cosmo/RadioBeam/interfconfigs.h

@@ -16 +16 @@
 //-----------------------------------------------------------------------------------
 
+//  Filled square array of ndxnd dishes 
 vector<Dish> CreateFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
+//  Semi filled square array of ndxnd dishes 
 vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
+//  Various semi-filled dish configurations  
 vector<Dish> CreateConfigA(double Ddish=5., double Eta=0.9);
 vector<Dish> CreateConfigB(double Ddish=5., double Eta=0.9);
@@ -24 +27 @@
 
 
+// Filled cylinder configuration 
 vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW=10., double cylRL=0.5, 
                                    double etaW=0.9, double etaRL=0.9, bool fgscid=true);
 
-
+// ASKPAP like - double hexagonal array of 12 dishes 
+vector<Dish> CreateDoubleHexagonConfig(double Ddish=12., double radius1=60., double radius2=150., double Eta=0.95);  
 
 #endif

trunk/Cosmo/RadioBeam/mdish.cc

@@ -54 +54 @@
   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());
+  Histo2D h2(-kxmx,kxmx,nx,-kymx,kymx,ny);
+
+  double xbc,ybc;
+  for(int_4 j=0; j<h2.NBinY(); j++) 
+    for(int_4 i=0; i<h2.NBinX(); i++) {
+      h2.BinCenter(i,j,xbc,ybc);
+      h2(i,j) = Value(xbc,ybc);
+    }
   return h2;    
 }

trunk/Cosmo/RadioBeam/mdish.h

@@ -47 +47 @@
   inline  double operator()(double kx, double ky) 
   { return Value(kx, ky); } 
-  virtual Histo2D GetResponse(int nx=256, int ny=256);  
+  virtual Histo2D GetResponse(int nx=255, int ny=255);  
   inline double D() { return dx_; } ; 
   inline double Dx() { return dx_; } ; 
@@ -105 +105 @@
   inline double Diameter() { return D; }
   inline double DiameterX() { return Dx; }
-  inline double DiameterY() { return Dx; }
+  inline double DiameterY() { return Dy; }
 
   int id_;   // numero de reflecteur

trunk/Cosmo/RadioBeam/pknoise.cc

@@ +1 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de calcul du spectre de puissance de bruit pour 
+  un interferometre (spectre moyenne sur 3D -> P_noise(k) )
+
+    R. Ansari , C. Magneville - Juin 2010 
+
+  Usage:  pknoise Diameter/Four2DRespTableFile OutPPFName 
+                 [NGen=1] [S2Cut=100.] [z_redshift=0.7]
+---------------------------------------------------------------  */
+
 #include "machdefs.h"
 #include "sopnamsp.h"
@@ -10 +20 @@
 #include "specpk.h"
 #include "interfconfigs.h"
+#include "radutil.h"
 
 #include "histinit.h"
@@ -23 +34 @@
 
 // ---------------------------------------------------------------------
-// Test main program for computing interferometer noise power spectrum
-// R. Ansari  - Avril 2010 
+// Test main 
+// R. Ansari  - Avril-Dec 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;
-};
 
 
@@ -64 +45 @@
 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;
-
+  if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) )  {
+    cout << " Usage:  pknoise Diameter/Four2DRespTableFile OutPPFName [RenormalizeMax]\n" 
+         << "         [NGen=1] [S2Cut=100.] [z_redshift=0.7] [PrtLev,PrtModulo=0,10]" << endl;
     return 1;
   }
-  cout << " ==== pknoise.cc program , test of SpectralShape and MassDist2D classes ==== " << endl;
+  cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
   // make sure SOPHYA modules are initialized 
   SophyaInit();  
@@ -77 +56 @@
   InitTim();
   //--- decoding command line arguments 
-  string outfile = "pknoise.ppf";  
-  if (narg>1) outfile = arg[1];
+  string tits="pknoise";
+  char cbuff[64];
+  bool fgresptbl=false;
+  double DIAMETRE=100.;
+  string resptblname;
+  if (isdigit(*arg[1])) {
+    fgresptbl=false;
+    DIAMETRE=atof(arg[1]);
+    sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
+  }
+  else { 
+    resptblname=arg[1];
+    fgresptbl=true;
+    sprintf(cbuff,"pknoise_RespTblName=%s", arg[1]);
+  }
+  tits=cbuff;
+  string outfile = arg[2];  
   if (outfile==".")  outfile = "pknoise.ppf";
+  bool fgrenorm=false;
+  double rmax=1.;
+  if (narg>3) {
+    rmax=atof(arg[3]);
+    fgrenorm=true;
+  }
   int NMAX = 1;
-  if (narg>2) NMAX = atoi(arg[2]);
+  if (narg>4) NMAX = atoi(arg[4]);
   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]);
-
+  if (narg>5) SCut = atof(arg[5]);
+  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
+  if (narg>6) z_Redshift = atof(arg[6]);
+  int prtlev=0;
+  int prtmod=10;
+  if (narg>7) sscanf(arg[7],"%d,%d",&prtlev,&prtmod);
   //-- end command line arguments
   
   int rc = 1;  
   try {  // exception handling try bloc at top level
-    cout << "0/ pknoise.cc: Executing, output file= " << outfile << endl;  
+    cout << " pknoise[0] : 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;
+    H21Conversions conv;
+    conv.setRedshift(z_Redshift);
+    double lambda = conv.getLambda();
 
-    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 = CreateFilledSqConfig(20, Ddish, Eta);
-    vector<Dish> vdpl64 = CreateFilledSqConfig(8, Ddish, Eta);
+    double Dol = DIAMETRE/lambda;
 
-    vector<Dish> vdsparse = CreateConfigA(Ddish, Eta);
-    vector<Dish> vdsparseD7 = CreateConfigA(Ddish2, Eta);
-    //    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
-    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
-    vector<Dish> vdsparseC = CreateConfigC(Ddish, Eta);
-
-
-    double cylW=12.;   // Largeur des cylindres 
-    double cylRL=0.5;  // Longeur des elements de reception le long du cylindre
-    double etaW=0.95;   // Efficacite de couverture en largeur
-    double etaRL=0.9;   // Efficacite de couverture le long du cylindre
-    vector<Dish> vcylplein = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
-    vector<Dish> vcylplP = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);
-
-    cylW=10.;
-    cylRL=0.5; 
-    vector<Dish> v3cyl = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
-    vector<Dish> v3cylP = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
-    cylW=25.;
-    cylRL=0.5; 
-    etaW=0.3;
-    etaRL=0.9;
-    vector<Dish> v2cyl = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, true);
-    vector<Dish> v2cylP = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
-
-    double LMAX = D;
-    bool fgnoauto = true;
-    int NRX=100;
-    int NRY=100;
-
-    MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
-    mdfill.SetRespHisNBins(NRX,NRY);   
-    Histo2D hrfill = mdfill.GetResponse();
-    PrtTim("Apres mdfill.GetResponse()");
-
-    MultiDish mdfill64(lambda, LMAX, vdpl64, fgnoauto);
-    mdfill64.SetRespHisNBins(NRX,NRY);   
-    {
-      Histo2D hpos=mdfill64.PosDist(10,10,10.*Ddish);
-      po << PPFNameTag("posf64") << hpos;
+    Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
+    Four2DResponse* arep_p=&arep;
+    Four2DRespTable resptbl;
+    if (fgresptbl) {
+      cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
+      resptbl.readFromPPF(resptblname);
+      arep_p=&resptbl;
+      if (fgrenorm) {
+        cout << " pknoise[1.b] call to resptbl.renormalize(" << rmax << ")"; 
+        double omax=resptbl.renormalize(rmax);
+        cout << " ... Old Max=" << omax << endl;
+      }
     }
-    Histo2D hrf64 = mdfill64.GetResponse();
-    PrtTim("Apres mdfill64.GetResponse()");
-
-    MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto); 
-    mdsparse.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
-    mdsparse.SetRespHisNBins(NRX,NRY);
-    {
-      Histo2D hpos=mdsparse.PosDist(22,22,22.*Ddish);;
-      po << PPFNameTag("posspA") << hpos;
-    }
-    Histo2D hrsp = mdsparse.GetResponse();
-    PrtTim("Apres mdsparse.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 mdsparseB(lambda, LMAX, vdsparseB, fgnoauto); 
-    mdsparseB.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
-    mdsparseB.SetRespHisNBins(NRX,NRY);
-    {
-      Histo2D hpos=mdsparseB.PosDist(15,15,15.*Ddish);
-      po << PPFNameTag("posspB") << hpos;
-    }
-    Histo2D hrspB = mdsparseB.GetResponse();
-    PrtTim("Apres mdsparseB.GetResponse()");                        
-
-    MultiDish mdsparseC(lambda, LMAX, vdsparseC, fgnoauto); 
-    mdsparseC.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
-    mdsparseC.SetRespHisNBins(NRX,NRY);
-    {
-      Histo2D hpos=mdsparseC.PosDist(20,20,20.*Ddish);
-      po << PPFNameTag("posspC") << hpos;
-    }
-    Histo2D hrspC = mdsparseC.GetResponse();
-    PrtTim("Apres mdsparseC.GetResponse()");                        
-
-    MultiDish mdsparseBfp(lambda, LMAX, vdsparseB, fgnoauto); 
-    mdsparseBfp.SetRespHisNBins(NRX,NRY);
-    Histo2D hrspBfp = mdsparseBfp.GetResponse();
-    PrtTim("Apres mdsparseBfp.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();
-    PrtTim("Apres mcylfillP.GetResponse()");
-
-    MultiDish md3cyl(lambda, LMAX, v3cyl, fgnoauto);
-    md3cyl.SetRespHisNBins(NRX,NRY);
-    Histo2D h3cyl = md3cyl.GetResponse();
-    PrtTim("Apres md3cyl.GetResponse()");
-    MultiDish md3cylP(lambda, LMAX, v3cylP, fgnoauto);
-    md3cylP.SetRespHisNBins(NRX,NRY);
-    Histo2D h3cylP = md3cylP.GetResponse();
-    PrtTim("Apres md3cylP.GetResponse()");
-
-    MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
-    md2cyl.SetRespHisNBins(NRX,NRY);
-    Histo2D h2cyl = md2cyl.GetResponse();
-    PrtTim("Apres md2cyl.GetResponse()");
-    MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
-    md2cylP.SetRespHisNBins(NRX,NRY);
-    Histo2D h2cylP = md2cylP.GetResponse();
-    PrtTim("Apres md2cylP.GetResponse()");
-
-    po << PPFNameTag("mfill") << hrfill;
-    po << PPFNameTag("mfill64") << hrf64;
-    po << PPFNameTag("mspars") << hrsp;
-    //    po << PPFNameTag("msparsd7") << hrspd7;
-    po << PPFNameTag("msparsB") << hrspB;
-    po << PPFNameTag("msparsC") << hrspC;
-    po << PPFNameTag("msparsBfp") << hrspBfp;
-    po << PPFNameTag("mcylf") << hfcyl;
-    po << PPFNameTag("m3cyl") << h3cyl;
-    po << PPFNameTag("m2cyl") << h2cyl;
-    po << PPFNameTag("mcylfP") << hfcylP;
-    po << PPFNameTag("m3cylP") << h3cylP;
-    po << PPFNameTag("m2cylP") << h2cylP;
-
-    PrtTim("Done computing multi-dish response");
- 
-
-    Four2DRespTable mdf(hrfill, Dol);
-    Four2DRespTable mdf64(hrf64, Dol);
-    Four2DRespTable mds(hrsp,   Dol);
-    // Four2DRespTable mdsfp(hrspfp,   Dol);
-    //    Four2DRespTable mdsd7(hrspd7,   Dol);
-    Four2DRespTable mdsB(hrspB,   Dol);
-    Four2DRespTable mdsC(hrspC,   Dol);
-    Four2DRespTable mdsBfp(hrspBfp,   Dol);
-
-    Four2DRespTable mcylf(hfcyl, Dol);
-    Four2DRespTable m3cyl(h3cyl, Dol);
-    Four2DRespTable m2cyl(h2cyl, Dol);
-    Four2DRespTable mcylfP(hfcylP, Dol);
-    Four2DRespTable m3cylP(h3cylP, Dol);
-    Four2DRespTable m2cylP(h2cylP, Dol);
+    else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
+              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
 
     
-    cout << " 3.a/ Instanciating object type Four3DPk  " << endl;
+    cout << " pknoise[2]: 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;
+    m3d.SetPrtLevel(prtlev,prtmod);
 
-    cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
-#define NCONFIG 14
-    Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mdf64, &mds, &mdsB, &mdsC, &mdsBfp,
-                                    &mcylf, &mcylfP, &m3cyl, &m3cylP, &m2cyl, &m2cylP};
-    const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","F8x8Dish5m",
-                               "S68Dish5m","S72Dish5m","S129CDish5m","S72BDish5mFP",
-                               "F8Cyl","F8CylP","F3Cyl","F3CylP","BiCyl","BiCylP"};
-    const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemdf64","noisemds",
-                               "noisemdsB","noisemdsC","noisemdsBfp",
-                               "noisefcyl","noisefcylP","noise3cyl","noise3cylP", "noise2cyl","noise2cylP"};
-    vector<int> nbdishes;
-    nbdishes.push_back(1);
-    nbdishes.push_back(1);
-    nbdishes.push_back(vdplein.size());
-    nbdishes.push_back(vdpl64.size());
-    nbdishes.push_back(vdsparse.size());
-    // nbdishes.push_back(vdsparse.size());
-    nbdishes.push_back(vdsparseB.size());
-    nbdishes.push_back(vdsparseC.size());
-    nbdishes.push_back(vdsparseB.size());
-    nbdishes.push_back(vcylplein.size());
-    nbdishes.push_back(vcylplP.size());
-    nbdishes.push_back(v3cyl.size());
-    nbdishes.push_back(v3cylP.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;
-    } 
+    cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
+    PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
+    pkn.SetPrtLevel(prtlev,prtmod);
+    HProf hpn = pkn.Compute();
+    po << hpn;
+ 
     rc = 0;
   }  // End of try bloc 

trunk/Cosmo/RadioBeam/repicon.cc

@@ +1 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Calcul de la reponse 2D (plan (u,v) d'un interferometre 
+    R. Ansari , C. Magneville - Juin-Dec 2010 
+
+  Usage:  repicon configId OutPPFName [z_Redshift=0.7] [RenormalizeMax]
+---------------------------------------------------------------  */
+
 #include "machdefs.h"
 #include "sopnamsp.h"
@@ -10 +17 @@
 #include "specpk.h"
 #include "interfconfigs.h"
-
+#include "radutil.h"
+
+#include "ntuple.h"
 #include "histinit.h"
 // #include "fiosinit.h"
@@ -22 +31 @@
 typedef DR48RandGen RandomGenerator ;
 
+// pour sauver la reponse mdresp et la config des dishes dans un fichier PPF 
+void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp);
+
 // ---------------------------------------------------------------------
-// Test main program for computing interferometer response un (u,v) plane
+// main program for computing interferometer response un (u,v) plane
 // R. Ansari  - Avril-Juin 2010 
 // ---------------------------------------------------------------------
 
-
-
-
-//-------------------------------------------------------------------------
 //      ------------------ MAIN PROGRAM ------------------------------
-//-------------------------------------------------------------------------
 int main(int narg, const char* arg[])
 {
   if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
-    cout << " Usage: repicon configId OutPPFName [Lambda=0.357] [RenormalizeMax]" 
+    cout << " Usage: repicon configId OutPPFName [z_redshift=0.7] [RenormalizeMax] \n" 
          << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes \n" 
-         << "          confA , confB, confC : semi-filled array of dishes \n"
+         << "          confA , confB, confC, confD : semi-filled array of dishes \n"
+         << "          hex12 : ASKAP like double hexagonal array of dishes \n"
          << "          f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders " << endl; 
     return 1;
@@ -73 +81 @@
     double etaW=0.95;   // Efficacite de couverture en largeur
     double etaRL=0.9;   // Efficacite de couverture le long du cylindre
+
+    double D = 100.;    // Taille de la zone
+
     int cnt=0;
 
@@ -116 +127 @@
       vdishes=CreateConfigC(Ddish, Eta);
     }
+    else if (config=="confD") {
+      fgpoint=true;
+      vdishes=CreateConfigD(Ddish, Eta);
+    }
+    else if (config=="hex12") {
+      fgpoint=true;
+      Ddish = 12.;
+      Eta=0.95;
+      D=350.;
+      vdishes=CreateDoubleHexagonConfig();
+    }
     else {
       cout << " NON valid configuration option -> exit" << endl;
@@ -121 +143 @@
     }
 
-    double D = 100.;
     double Dol = D/LAMBDA;
-
     double LMAX = D;
     bool fgnoauto = true;
@@ -146 +166 @@
     cout << " repicon[3] : saving Four2DRespTable for config " << config << " to PPF file " << outfile << endl;
     mdresp.writeToPPF(outfile);
+
+    string outfile2 = "hdt_"+outfile;
+    cout << " repicon[4] : saving H2D-response, multidish config to PPF file " << outfile2 << endl;
+    SaveDTVecDishesH2Resp(outfile2, vdishes, mdresp);
 
     rc = 0;
@@ -167 +191 @@
 }
 
-
+/*-- Nouvelle-Fonction --*/
+void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp)
+{
+  char* names[5]={"did","posx","posy","diam","diamy"};
+  NTuple ntvd(5,names,64,false);
+  r_4 xnt[10];
+  for(size_t i=0; i<vdishes.size(); i++) {
+    xnt[0]=vdishes[i].ReflectorId();
+    xnt[1]=vdishes[i].X;
+    xnt[2]=vdishes[i].Y;
+    if (vdishes[i].isCircular()) {
+      xnt[3]=vdishes[i].Diameter();
+      xnt[4]=0.;
+    }
+    else {
+      xnt[3]=vdishes[i].DiameterX();
+      xnt[4]=vdishes[i].DiameterY();      
+    }
+    ntvd.Fill(xnt);
+  }
+  Histo2D h2rep=mdresp.GetResponse();
+  POutPersist po(outfile);
+  po << PPFNameTag("mdish") << ntvd;
+  po << PPFNameTag("h2rep") << h2rep;
+  return;
+}
+

trunk/Cosmo/RadioBeam/specpk.cc

@@ -1 @@
-//  Classes to compute 3D power spectrum 
-// R. Ansari - Nov 2008, May 2010 
+
+/*  ------------------------ Projet BAORadio -------------------- 
+    Classes to compute 3D power spectrum and noise power spectrum
+    R. Ansari - Nov 2008 ... Dec 2010 
+---------------------------------------------------------------  */
 
 #include "specpk.h"
 #include "randr48.h"      
+#include "ctimer.h"      
 
 //------------------------------------
@@ -154 +158 @@
     }
   }
-  if (prtlev_>0)
+  if (prtlev_>2)
     cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
 }
@@ -182 +186 @@
     }
   }
-  if (prtlev_>1)  fourAmp.Show();
+  if (prtlev_>2)  fourAmp.Show();
   if (crmask) {
     POutPersist po("mask.ppf");
@@ -233 +237 @@
 void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
 {
-
+  uint_8 nmodeok=0;
   // fourAmp represent 3-D fourier transform of a real input array. 
   // The second half of the array along Y and Z contain negative frequencies
@@ -251 +255 @@
         if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
         hp.Add(wk, amp2);
+        nmodeok++;
       }
     }
   }
+  if (prtlev_>1) {
+    cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size() 
+         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
+  }
   return;
 }
 
+//-----------------------------------------------------
+// -- MassDist2D class :  2D mass distribution 
+// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k) 
+// determine par une reponse 2D de l'instrument
+//-----------------------------------------------------
+PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen, 
+                                     const char* tit)
+  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
+{
+  SetPrtLevel();
+}
+
+HProf PkNoiseCalculator::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);
+    if ((prtlev_>0)&&(((igen+1)%prtmodulo_)==0)) 
+      cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
+  }
+  return hnd;
+}
 
 

trunk/Cosmo/RadioBeam/specpk.h

@@ -48 +48 @@
   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; }
+  inline int SetPrtLevel(int lev=0, int prtmod=10) 
+  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
   void ComputeFourierAmp(SpectralShape& pk);
   void ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask=false);
@@ -71 +71 @@
   double dkx_, dky_, dkz_;
   int prtlev_;
+  int prtmodulo_;
 };
+
+// --- PkNoiseCalculator : 
+// - Classe de calcul du spectre de bruit PNoise(k) determine par une reponse
+//   2D de l'instrument 
+class PkNoiseCalculator
+{
+public:
+  PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut=100., int ngen=1, const char* tit="PkNoise");
+  
+  inline void SetS2Cut(double s2cut=100.)
+  {  S2CUT=s2cut; }
+  inline double GetS2Cut() { return S2CUT; }
+  HProf Compute();
+  inline int SetPrtLevel(int lev=0, int prtmod=10) 
+  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
+
+protected:
+  Four3DPk& pkn3d;
+  Four2DResponse& frep;
+  double S2CUT;
+  int NGEN;
+  string title;
+  int prtlev_;
+  int prtmodulo_;
+};
+