Changeset 3930 in Sophya for trunk/Cosmo/RadioBeam


Ignore:
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:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/interfconfigs.cc

    r3792 r3930  
    234234}
    235235
    236  cout << "CreateConfigB/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
     236 cout << "CreateConfigD/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
    237237 
    238238 vector<Dish> vd;
    … …  
    265265  return vd;
    266266}
     267
     268/* --Fonction -- */
     269vector<Dish> CreateDoubleHexagonConfig(double Ddish, double radius1, double radius2, double Eta)
     270{
     271  vector<Dish>  vd;
     272  int cnt=0;
     273  double ang=0.;
     274  for(int i=0; i<6; i++) {
     275    double x=cos(ang)*radius1;
     276    double y=sin(ang)*radius1;
     277    cnt++;
     278    vd.push_back(Dish(cnt, x, y, Eta*Ddish));
     279    ang += M_PI/3.;
     280  }
     281  ang=M_PI/6.;
     282  for(int i=0; i<6; i++) {
     283    double x=cos(ang)*radius2;
     284    double y=sin(ang)*radius2;
     285    cnt++;
     286    vd.push_back(Dish(cnt, x, y, Eta*Ddish));
     287    ang += M_PI/3.;
     288  }
     289 cout << ">>>CreateDoubleHexagonConfig(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     290 return vd;
     291}

  • trunk/Cosmo/RadioBeam/interfconfigs.h

    r3792 r3930  
    1616//-----------------------------------------------------------------------------------
    1717
     18//  Filled square array of ndxnd dishes
    1819vector<Dish> CreateFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
     20//  Semi filled square array of ndxnd dishes
    1921vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
     22//  Various semi-filled dish configurations 
    2023vector<Dish> CreateConfigA(double Ddish=5., double Eta=0.9);
    2124vector<Dish> CreateConfigB(double Ddish=5., double Eta=0.9);
    … …  
    2427
    2528
     29// Filled cylinder configuration
    2630vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW=10., double cylRL=0.5,
    2731                                   double etaW=0.9, double etaRL=0.9, bool fgscid=true);
    2832
    29 
     33// ASKPAP like - double hexagonal array of 12 dishes
     34vector<Dish> CreateDoubleHexagonConfig(double Ddish=12., double radius1=60., double radius2=150., double Eta=0.95); 
    3035
    3136#endif

  • trunk/Cosmo/RadioBeam/mdish.cc

    r3796 r3930  
    5454  double kymx = 1.2*DeuxPI*dy_;
    5555  if (typ_<3) kymx=kxmx;
    56   Histo2D h2(0.,kxmx,nx,0.,kymx,ny);
    57 
    58   for(int j=0; j<h2.NBinY(); j++)
    59     for(int i=0; i<h2.NBinX(); i++)
    60       h2(i,j) = Value((i+0.5)*h2.WBinX(), (j+0.5)*h2.WBinY());
     56  Histo2D h2(-kxmx,kxmx,nx,-kymx,kymx,ny);
     57
     58  double xbc,ybc;
     59  for(int_4 j=0; j<h2.NBinY(); j++)
     60    for(int_4 i=0; i<h2.NBinX(); i++) {
     61      h2.BinCenter(i,j,xbc,ybc);
     62      h2(i,j) = Value(xbc,ybc);
     63    }
    6164  return h2;   
    6265}

  • trunk/Cosmo/RadioBeam/mdish.h

    r3797 r3930  
    4747  inline  double operator()(double kx, double ky)
    4848  { return Value(kx, ky); }
    49   virtual Histo2D GetResponse(int nx=256, int ny=256); 
     49  virtual Histo2D GetResponse(int nx=255, int ny=255); 
    5050  inline double D() { return dx_; } ;
    5151  inline double Dx() { return dx_; } ;
    … …  
    105105  inline double Diameter() { return D; }
    106106  inline double DiameterX() { return Dx; }
    107   inline double DiameterY() { return Dx; }
     107  inline double DiameterY() { return Dy; }
    108108
    109109  int id_;   // numero de reflecteur

  • trunk/Cosmo/RadioBeam/pknoise.cc

    r3792 r3930  
     1/*  ------------------------ Projet BAORadio --------------------
     2  Programme de calcul du spectre de puissance de bruit pour
     3  un interferometre (spectre moyenne sur 3D -> P_noise(k) )
     4
     5    R. Ansari , C. Magneville - Juin 2010
     6
     7  Usage:  pknoise Diameter/Four2DRespTableFile OutPPFName
     8                 [NGen=1] [S2Cut=100.] [z_redshift=0.7]
     9---------------------------------------------------------------  */
     10
    111#include "machdefs.h"
    212#include "sopnamsp.h"
    … …  
    1020#include "specpk.h"
    1121#include "interfconfigs.h"
     22#include "radutil.h"
    1223
    1324#include "histinit.h"
    … …  
    2334
    2435// ---------------------------------------------------------------------
    25 // Test main program for computing interferometer noise power spectrum
    26 // R. Ansari  - Avril 2010
     36// Test main
     37// R. Ansari  - Avril-Dec 2010
    2738// ---------------------------------------------------------------------
    28 
    29 class PkNoiseCalculator
    30 {
    31 public:
    32   PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep,
    33                     double s2cut=100., int ngen=1, const char* tit="PkNoise")
    34     : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit)
    35   { }
    36   inline void SetS2Cut(double s2cut=100.)
    37   {  S2CUT=s2cut; }
    38   HProf Compute()
    39   {
    40     Timer tm(title.c_str());
    41     tm.Nop();
    42     HProf hnd;
    43     cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
    44     for(int igen=0; igen<NGEN; igen++) {
    45       pkn3d.ComputeNoiseFourierAmp(frep);
    46       if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
    47       else pkn3d.ComputePkCumul(hnd,S2CUT);
    48     }
    49     return hnd;
    50   }
    51 
    52   Four3DPk& pkn3d;
    53   Four2DResponse& frep;
    54   double S2CUT;
    55   int NGEN;
    56   string title;
    57 };
    5839
    5940
    … …  
    6445int main(int narg, const char* arg[])
    6546{
    66   if ((narg>1)&&(strcmp(arg[1],"-h")==0))  {
    67     cout<< " Usage: pknoise [OutPPFName NGen S2Cut Lambda] " << endl;
    68     cout<< "   Default: OutPPFName=pknoise.ppf, NGen=1 " << endl;
    69     cout<< "            S2CUT=0. , Lambda=0.357  \n" << endl;
    70 
     47  if ( (narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0)) )  {
     48    cout << " Usage:  pknoise Diameter/Four2DRespTableFile OutPPFName [RenormalizeMax]\n"
     49         << "         [NGen=1] [S2Cut=100.] [z_redshift=0.7] [PrtLev,PrtModulo=0,10]" << endl;
    7150    return 1;
    7251  }
    73   cout << " ==== pknoise.cc program , test of SpectralShape and MassDist2D classes ==== " << endl;
     52  cout << " ==== pknoise.cc : interferometer noise power spectrum computation ==== " << endl;
    7453  // make sure SOPHYA modules are initialized
    7554  SophyaInit(); 
    … …  
    7756  InitTim();
    7857  //--- decoding command line arguments
    79   string outfile = "pknoise.ppf"; 
    80   if (narg>1) outfile = arg[1];
     58  string tits="pknoise";
     59  char cbuff[64];
     60  bool fgresptbl=false;
     61  double DIAMETRE=100.;
     62  string resptblname;
     63  if (isdigit(*arg[1])) {
     64    fgresptbl=false;
     65    DIAMETRE=atof(arg[1]);
     66    sprintf(cbuff,"pknoise_Dish(%g m)", DIAMETRE);
     67  }
     68  else {
     69    resptblname=arg[1];
     70    fgresptbl=true;
     71    sprintf(cbuff,"pknoise_RespTblName=%s", arg[1]);
     72  }
     73  tits=cbuff;
     74  string outfile = arg[2]; 
    8175  if (outfile==".")  outfile = "pknoise.ppf";
     76  bool fgrenorm=false;
     77  double rmax=1.;
     78  if (narg>3) {
     79    rmax=atof(arg[3]);
     80    fgrenorm=true;
     81  }
    8282  int NMAX = 1;
    83   if (narg>2) NMAX = atoi(arg[2]);
     83  if (narg>4) NMAX = atoi(arg[4]);
    8484  double SCut=0.;
    85   if (narg>3) SCut = atof(arg[3]);
    86   double LAMBDA=0.357 ;  // 21 cm at z=0.7
    87   if (narg>4) LAMBDA = atof(arg[4]);
    88 
     85  if (narg>5) SCut = atof(arg[5]);
     86  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m 
     87  if (narg>6) z_Redshift = atof(arg[6]);
     88  int prtlev=0;
     89  int prtmod=10;
     90  if (narg>7) sscanf(arg[7],"%d,%d",&prtlev,&prtmod);
    8991  //-- end command line arguments
    9092 
    9193  int rc = 1; 
    9294  try {  // exception handling try bloc at top level
    93     cout << "0/ pknoise.cc: Executing, output file= " << outfile << endl; 
     95    cout << " pknoise[0] : Executing, output file= " << outfile << endl; 
    9496    POutPersist po(outfile);
    95     cout << " 1.a/ Instanciating object type SpectralShape " << endl;
    96     SpectralShape spec(2);
    97     cout << " 1.b/ Wrinting spectral shape vector (name= Pk) to output PPF " << endl;
    98     Histo hpk = spec.GetPk(1024);
    99     po << PPFNameTag("Pk") << hpk;
    10097
    101     double D = 100.;
    102     double lambda = LAMBDA;
    103     double Dol = D/lambda;
    104     cout << " 2.a/ Instanciating Four2DResponse(1/2/3...) " << endl;
    105     Four2DResponse dishg(1,Dol,Dol);
    106     Four2DResponse dish(2,Dol,Dol);
    107     Four2DResponse dish2(2,Dol*2.,Dol*2.);
    108     Four2DResponse dishsq(3,Dol,Dol/5.);
    109     cout << " 2.b/ Writing Four2DResponse Histo2D to output ppf " << endl;
    110     Histo2D hdg = dishg.GetResponse();
    111     Histo2D hd = dish.GetResponse();
    112     Histo2D hd2 = dish2.GetResponse();
    113     Histo2D hdsq = dishsq.GetResponse();
    114     po << PPFNameTag("dishg") << hdg;
    115     po << PPFNameTag("dish") << hd;
    116     po << PPFNameTag("dish2") << hd2;
    117     po << PPFNameTag("dishsq") << hdsq;
     98    H21Conversions conv;
     99    conv.setRedshift(z_Redshift);
     100    double lambda = conv.getLambda();
    118101
    119     cout << " 2.c/ Creating MultiDish Filled Array " << endl;
    120     double Ddish = 5.;
    121     double Ddish2 = 7.5;
    122     double Eta=0.95;
    123     int cnt=0;
    124     vector<Dish> vdplein = CreateFilledSqConfig(20, Ddish, Eta);
    125     vector<Dish> vdpl64 = CreateFilledSqConfig(8, Ddish, Eta);
     102    double Dol = DIAMETRE/lambda;
    126103
    127     vector<Dish> vdsparse = CreateConfigA(Ddish, Eta);
    128     vector<Dish> vdsparseD7 = CreateConfigA(Ddish2, Eta);
    129     //    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
    130     vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
    131     vector<Dish> vdsparseC = CreateConfigC(Ddish, Eta);
    132 
    133 
    134     double cylW=12.;   // Largeur des cylindres
    135     double cylRL=0.5;  // Longeur des elements de reception le long du cylindre
    136     double etaW=0.95;   // Efficacite de couverture en largeur
    137     double etaRL=0.9;   // Efficacite de couverture le long du cylindre
    138     vector<Dish> vcylplein = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
    139     vector<Dish> vcylplP = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);
    140 
    141     cylW=10.;
    142     cylRL=0.5;
    143     vector<Dish> v3cyl = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
    144     vector<Dish> v3cylP = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
    145     cylW=25.;
    146     cylRL=0.5;
    147     etaW=0.3;
    148     etaRL=0.9;
    149     vector<Dish> v2cyl = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, true);
    150     vector<Dish> v2cylP = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
    151 
    152     double LMAX = D;
    153     bool fgnoauto = true;
    154     int NRX=100;
    155     int NRY=100;
    156 
    157     MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
    158     mdfill.SetRespHisNBins(NRX,NRY);   
    159     Histo2D hrfill = mdfill.GetResponse();
    160     PrtTim("Apres mdfill.GetResponse()");
    161 
    162     MultiDish mdfill64(lambda, LMAX, vdpl64, fgnoauto);
    163     mdfill64.SetRespHisNBins(NRX,NRY);   
    164     {
    165       Histo2D hpos=mdfill64.PosDist(10,10,10.*Ddish);
    166       po << PPFNameTag("posf64") << hpos;
     104    Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
     105    Four2DResponse* arep_p=&arep;
     106    Four2DRespTable resptbl;
     107    if (fgresptbl) {
     108      cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
     109      resptbl.readFromPPF(resptblname);
     110      arep_p=&resptbl;
     111      if (fgrenorm) {
     112        cout << " pknoise[1.b] call to resptbl.renormalize(" << rmax << ")";
     113        double omax=resptbl.renormalize(rmax);
     114        cout << " ... Old Max=" << omax << endl;
     115      }
    167116    }
    168     Histo2D hrf64 = mdfill64.GetResponse();
    169     PrtTim("Apres mdfill64.GetResponse()");
    170 
    171     MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto);
    172     mdsparse.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    173     mdsparse.SetRespHisNBins(NRX,NRY);
    174     {
    175       Histo2D hpos=mdsparse.PosDist(22,22,22.*Ddish);;
    176       po << PPFNameTag("posspA") << hpos;
    177     }
    178     Histo2D hrsp = mdsparse.GetResponse();
    179     PrtTim("Apres mdsparse.GetResponse()");
    180 
    181     /*
    182     MultiDish mdsparseD7(lambda, LMAX, vdsparseD7, fgnoauto);
    183     mdsparseD7.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
    184     mdsparseD7.SetRespHisNBins(NRX,NRY);
    185     Histo2D hrspd7 = mdsparseD7.GetResponse();
    186     PrtTim("Apres mdsparseD7.GetResponse()");
    187     */
    188     MultiDish mdsparseB(lambda, LMAX, vdsparseB, fgnoauto);
    189     mdsparseB.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    190     mdsparseB.SetRespHisNBins(NRX,NRY);
    191     {
    192       Histo2D hpos=mdsparseB.PosDist(15,15,15.*Ddish);
    193       po << PPFNameTag("posspB") << hpos;
    194     }
    195     Histo2D hrspB = mdsparseB.GetResponse();
    196     PrtTim("Apres mdsparseB.GetResponse()");                       
    197 
    198     MultiDish mdsparseC(lambda, LMAX, vdsparseC, fgnoauto);
    199     mdsparseC.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    200     mdsparseC.SetRespHisNBins(NRX,NRY);
    201     {
    202       Histo2D hpos=mdsparseC.PosDist(20,20,20.*Ddish);
    203       po << PPFNameTag("posspC") << hpos;
    204     }
    205     Histo2D hrspC = mdsparseC.GetResponse();
    206     PrtTim("Apres mdsparseC.GetResponse()");                       
    207 
    208     MultiDish mdsparseBfp(lambda, LMAX, vdsparseB, fgnoauto);
    209     mdsparseBfp.SetRespHisNBins(NRX,NRY);
    210     Histo2D hrspBfp = mdsparseBfp.GetResponse();
    211     PrtTim("Apres mdsparseBfp.GetResponse()");                     
    212 
    213 
    214     MultiDish mcylfill(lambda, LMAX, vcylplein, fgnoauto);
    215     mcylfill.SetRespHisNBins(NRX,NRY);
    216     Histo2D hfcyl = mcylfill.GetResponse();
    217     PrtTim("Apres mcylfill.GetResponse()");
    218     MultiDish mcylfillP(lambda, LMAX, vcylplP, fgnoauto);
    219     mcylfillP.SetRespHisNBins(NRX,NRY);
    220     Histo2D hfcylP = mcylfillP.GetResponse();
    221     PrtTim("Apres mcylfillP.GetResponse()");
    222 
    223     MultiDish md3cyl(lambda, LMAX, v3cyl, fgnoauto);
    224     md3cyl.SetRespHisNBins(NRX,NRY);
    225     Histo2D h3cyl = md3cyl.GetResponse();
    226     PrtTim("Apres md3cyl.GetResponse()");
    227     MultiDish md3cylP(lambda, LMAX, v3cylP, fgnoauto);
    228     md3cylP.SetRespHisNBins(NRX,NRY);
    229     Histo2D h3cylP = md3cylP.GetResponse();
    230     PrtTim("Apres md3cylP.GetResponse()");
    231 
    232     MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
    233     md2cyl.SetRespHisNBins(NRX,NRY);
    234     Histo2D h2cyl = md2cyl.GetResponse();
    235     PrtTim("Apres md2cyl.GetResponse()");
    236     MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
    237     md2cylP.SetRespHisNBins(NRX,NRY);
    238     Histo2D h2cylP = md2cylP.GetResponse();
    239     PrtTim("Apres md2cylP.GetResponse()");
    240 
    241     po << PPFNameTag("mfill") << hrfill;
    242     po << PPFNameTag("mfill64") << hrf64;
    243     po << PPFNameTag("mspars") << hrsp;
    244     //    po << PPFNameTag("msparsd7") << hrspd7;
    245     po << PPFNameTag("msparsB") << hrspB;
    246     po << PPFNameTag("msparsC") << hrspC;
    247     po << PPFNameTag("msparsBfp") << hrspBfp;
    248     po << PPFNameTag("mcylf") << hfcyl;
    249     po << PPFNameTag("m3cyl") << h3cyl;
    250     po << PPFNameTag("m2cyl") << h2cyl;
    251     po << PPFNameTag("mcylfP") << hfcylP;
    252     po << PPFNameTag("m3cylP") << h3cylP;
    253     po << PPFNameTag("m2cylP") << h2cylP;
    254 
    255     PrtTim("Done computing multi-dish response");
    256  
    257 
    258     Four2DRespTable mdf(hrfill, Dol);
    259     Four2DRespTable mdf64(hrf64, Dol);
    260     Four2DRespTable mds(hrsp,   Dol);
    261     // Four2DRespTable mdsfp(hrspfp,   Dol);
    262     //    Four2DRespTable mdsd7(hrspd7,   Dol);
    263     Four2DRespTable mdsB(hrspB,   Dol);
    264     Four2DRespTable mdsC(hrspC,   Dol);
    265     Four2DRespTable mdsBfp(hrspBfp,   Dol);
    266 
    267     Four2DRespTable mcylf(hfcyl, Dol);
    268     Four2DRespTable m3cyl(h3cyl, Dol);
    269     Four2DRespTable m2cyl(h2cyl, Dol);
    270     Four2DRespTable mcylfP(hfcylP, Dol);
    271     Four2DRespTable m3cylP(h3cylP, Dol);
    272     Four2DRespTable m2cylP(h2cylP, Dol);
     117    else cout << " pknoise[1]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
     118              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    273119
    274120   
    275     cout << " 3.a/ Instanciating object type Four3DPk  " << endl;
     121    cout << " pknoise[2]: Instanciating object type Four3DPk  " << endl;
    276122    RandomGenerator rg;
    277123    Four3DPk m3d(rg);
    278124    m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI);
    279     cout << " 3.b/ call ComputeFourierAmp()  NGEN=" << NMAX << endl;
    280     HProf hrpk;
    281     for(int igen=0; igen<NMAX; igen++) {
    282       m3d.ComputeFourierAmp(spec);
    283       if (igen==0)  hrpk = m3d.ComputePk();
    284       else m3d.ComputePkCumul(hrpk);
    285     }
    286     PrtTim("md.ComputeFourierAmp() done");
    287     po << PPFNameTag("recPk") << hrpk;
     125    m3d.SetPrtLevel(prtlev,prtmod);
    288126
    289     cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    290 #define NCONFIG 14
    291     Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mdf64, &mds, &mdsB, &mdsC, &mdsBfp,
    292                                     &mcylf, &mcylfP, &m3cyl, &m3cylP, &m2cyl, &m2cylP};
    293     const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","F8x8Dish5m",
    294                                "S68Dish5m","S72Dish5m","S129CDish5m","S72BDish5mFP",
    295                                "F8Cyl","F8CylP","F3Cyl","F3CylP","BiCyl","BiCylP"};
    296     const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemdf64","noisemds",
    297                                "noisemdsB","noisemdsC","noisemdsBfp",
    298                                "noisefcyl","noisefcylP","noise3cyl","noise3cylP", "noise2cyl","noise2cylP"};
    299     vector<int> nbdishes;
    300     nbdishes.push_back(1);
    301     nbdishes.push_back(1);
    302     nbdishes.push_back(vdplein.size());
    303     nbdishes.push_back(vdpl64.size());
    304     nbdishes.push_back(vdsparse.size());
    305     // nbdishes.push_back(vdsparse.size());
    306     nbdishes.push_back(vdsparseB.size());
    307     nbdishes.push_back(vdsparseC.size());
    308     nbdishes.push_back(vdsparseB.size());
    309     nbdishes.push_back(vcylplein.size());
    310     nbdishes.push_back(vcylplP.size());
    311     nbdishes.push_back(v3cyl.size());
    312     nbdishes.push_back(v3cylP.size());
    313     nbdishes.push_back(v2cyl.size());
    314     nbdishes.push_back(v2cylP.size());
    315 
    316     for(int lc=0; lc<NCONFIG; lc++) {
    317       PkNoiseCalculator pkn(m3d, *(f2rep[lc]), SCut/(double)nbdishes[lc], NMAX, tits[lc]);
    318       HProf hpn = pkn.Compute();
    319       po << PPFNameTag(tags[lc]) << hpn;
    320     }
     127    cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
     128    PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
     129    pkn.SetPrtLevel(prtlev,prtmod);
     130    HProf hpn = pkn.Compute();
     131    po << hpn;
     132 
    321133    rc = 0;
    322134  }  // End of try bloc

  • trunk/Cosmo/RadioBeam/repicon.cc

    r3796 r3930  
     1/*  ------------------------ Projet BAORadio --------------------
     2  Calcul de la reponse 2D (plan (u,v) d'un interferometre
     3    R. Ansari , C. Magneville - Juin-Dec 2010
     4
     5  Usage:  repicon configId OutPPFName [z_Redshift=0.7] [RenormalizeMax]
     6---------------------------------------------------------------  */
     7
    18#include "machdefs.h"
    29#include "sopnamsp.h"
    … …  
    1017#include "specpk.h"
    1118#include "interfconfigs.h"
    12 
     19#include "radutil.h"
     20
     21#include "ntuple.h"
    1322#include "histinit.h"
    1423// #include "fiosinit.h"
    … …  
    2231typedef DR48RandGen RandomGenerator ;
    2332
     33// pour sauver la reponse mdresp et la config des dishes dans un fichier PPF
     34void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp);
     35
    2436// ---------------------------------------------------------------------
    25 // Test main program for computing interferometer response un (u,v) plane
     37// main program for computing interferometer response un (u,v) plane
    2638// R. Ansari  - Avril-Juin 2010
    2739// ---------------------------------------------------------------------
    2840
    29 
    30 
    31 
    32 //-------------------------------------------------------------------------
    3341//      ------------------ MAIN PROGRAM ------------------------------
    34 //-------------------------------------------------------------------------
    3542int main(int narg, const char* arg[])
    3643{
    3744  if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
    38     cout << " Usage: repicon configId OutPPFName [Lambda=0.357] [RenormalizeMax]"
     45    cout << " Usage: repicon configId OutPPFName [z_redshift=0.7] [RenormalizeMax] \n"
    3946         << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes \n"
    40          << "          confA , confB, confC : semi-filled array of dishes \n"
     47         << "          confA , confB, confC, confD : semi-filled array of dishes \n"
     48         << "          hex12 : ASKAP like double hexagonal array of dishes \n"
    4149         << "          f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders " << endl;
    4250    return 1;
    … …  
    7381    double etaW=0.95;   // Efficacite de couverture en largeur
    7482    double etaRL=0.9;   // Efficacite de couverture le long du cylindre
     83
     84    double D = 100.;    // Taille de la zone
     85
    7586    int cnt=0;
    7687
    … …  
    116127      vdishes=CreateConfigC(Ddish, Eta);
    117128    }
     129    else if (config=="confD") {
     130      fgpoint=true;
     131      vdishes=CreateConfigD(Ddish, Eta);
     132    }
     133    else if (config=="hex12") {
     134      fgpoint=true;
     135      Ddish = 12.;
     136      Eta=0.95;
     137      D=350.;
     138      vdishes=CreateDoubleHexagonConfig();
     139    }
    118140    else {
    119141      cout << " NON valid configuration option -> exit" << endl;
    … …  
    121143    }
    122144
    123     double D = 100.;
    124145    double Dol = D/LAMBDA;
    125 
    126146    double LMAX = D;
    127147    bool fgnoauto = true;
    … …  
    146166    cout << " repicon[3] : saving Four2DRespTable for config " << config << " to PPF file " << outfile << endl;
    147167    mdresp.writeToPPF(outfile);
     168
     169    string outfile2 = "hdt_"+outfile;
     170    cout << " repicon[4] : saving H2D-response, multidish config to PPF file " << outfile2 << endl;
     171    SaveDTVecDishesH2Resp(outfile2, vdishes, mdresp);
    148172
    149173    rc = 0;
    … …  
    167191}
    168192
    169 
     193/*-- Nouvelle-Fonction --*/
     194void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp)
     195{
     196  char* names[5]={"did","posx","posy","diam","diamy"};
     197  NTuple ntvd(5,names,64,false);
     198  r_4 xnt[10];
     199  for(size_t i=0; i<vdishes.size(); i++) {
     200    xnt[0]=vdishes[i].ReflectorId();
     201    xnt[1]=vdishes[i].X;
     202    xnt[2]=vdishes[i].Y;
     203    if (vdishes[i].isCircular()) {
     204      xnt[3]=vdishes[i].Diameter();
     205      xnt[4]=0.;
     206    }
     207    else {
     208      xnt[3]=vdishes[i].DiameterX();
     209      xnt[4]=vdishes[i].DiameterY();     
     210    }
     211    ntvd.Fill(xnt);
     212  }
     213  Histo2D h2rep=mdresp.GetResponse();
     214  POutPersist po(outfile);
     215  po << PPFNameTag("mdish") << ntvd;
     216  po << PPFNameTag("h2rep") << h2rep;
     217  return;
     218}
     219

  • trunk/Cosmo/RadioBeam/specpk.cc

    r3825 r3930  
    1 //  Classes to compute 3D power spectrum
    2 // R. Ansari - Nov 2008, May 2010
     1
     2/*  ------------------------ Projet BAORadio --------------------
     3    Classes to compute 3D power spectrum and noise power spectrum
     4    R. Ansari - Nov 2008 ... Dec 2010
     5---------------------------------------------------------------  */
    36
    47#include "specpk.h"
    58#include "randr48.h"     
     9#include "ctimer.h"     
    610
    711//------------------------------------
    … …  
    154158    }
    155159  }
    156   if (prtlev_>0)
     160  if (prtlev_>2)
    157161    cout << " Four3DPk::ComputeFourierAmp() done ..." << endl;
    158162}
    … …  
    182186    }
    183187  }
    184   if (prtlev_>1)  fourAmp.Show();
     188  if (prtlev_>2)  fourAmp.Show();
    185189  if (crmask) {
    186190    POutPersist po("mask.ppf");
    … …  
    233237void Four3DPk::ComputePkCumul(HProf& hp, double s2cut)
    234238{
    235 
     239  uint_8 nmodeok=0;
    236240  // fourAmp represent 3-D fourier transform of a real input array.
    237241  // The second half of the array along Y and Z contain negative frequencies
    … …  
    251255        if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
    252256        hp.Add(wk, amp2);
     257        nmodeok++;
    253258      }
    254259    }
    255260  }
     261  if (prtlev_>1) {
     262    cout << " Four3DPk::ComputePkCumul/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
     263         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
     264  }
    256265  return;
    257266}
    258267
     268//-----------------------------------------------------
     269// -- MassDist2D class :  2D mass distribution
     270// --- PkNoiseCalculator : Classe de calcul du spectre de bruit PNoise(k)
     271// determine par une reponse 2D de l'instrument
     272//-----------------------------------------------------
     273PkNoiseCalculator::PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut, int ngen,
     274                                     const char* tit)
     275  : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit)
     276{
     277  SetPrtLevel();
     278}
     279
     280HProf PkNoiseCalculator::Compute()
     281{
     282  Timer tm(title.c_str());
     283  tm.Nop();
     284  HProf hnd;
     285  cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
     286  for(int igen=0; igen<NGEN; igen++) {
     287    pkn3d.ComputeNoiseFourierAmp(frep);
     288    if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
     289    else pkn3d.ComputePkCumul(hnd,S2CUT);
     290    if ((prtlev_>0)&&(((igen+1)%prtmodulo_)==0))
     291      cout << " PkNoiseCalculator::Compute() - done igen=" << igen << " / MaxNGen=" << NGEN << endl;
     292  }
     293  return hnd;
     294}
    259295
    260296

  • trunk/Cosmo/RadioBeam/specpk.h

    r3825 r3930  
    4848  inline void SetCellSize(double dkx=DeuxPI, double dky=DeuxPI, double dkz=DeuxPI)
    4949  { dkx_=dkx;  dky_=dky;  dkz_=dkz; }
    50   inline int SetPrtLevel(int lev=0)
    51   { int olev=prtlev_; prtlev_=lev; return olev; }
     50  inline int SetPrtLevel(int lev=0, int prtmod=10)
     51  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
    5252  void ComputeFourierAmp(SpectralShape& pk);
    5353  void ComputeNoiseFourierAmp(Four2DResponse& resp, bool crmask=false);
    … …  
    7171  double dkx_, dky_, dkz_;
    7272  int prtlev_;
     73  int prtmodulo_;
    7374};
     75
     76// --- PkNoiseCalculator :
     77// - Classe de calcul du spectre de bruit PNoise(k) determine par une reponse
     78//   2D de l'instrument
     79class PkNoiseCalculator
     80{
     81public:
     82  PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, double s2cut=100., int ngen=1, const char* tit="PkNoise");
     83 
     84  inline void SetS2Cut(double s2cut=100.)
     85  {  S2CUT=s2cut; }
     86  inline double GetS2Cut() { return S2CUT; }
     87  HProf Compute();
     88  inline int SetPrtLevel(int lev=0, int prtmod=10)
     89  { int olev=prtlev_; prtlev_=lev; prtmodulo_=prtmod; return olev; }
     90
     91protected:
     92  Four3DPk& pkn3d;
     93  Four2DResponse& frep;
     94  double S2CUT;
     95  int NGEN;
     96  string title;
     97  int prtlev_;
     98  int prtmodulo_;
     99};
     100
    74101
    75102
Note: See TracChangeset for help on using the changeset viewer.