Changeset 3930 in Sophya for trunk/Cosmo/RadioBeam/pknoise.cc


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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.