Changeset 3947 in Sophya for trunk/Cosmo


Ignore:
Timestamp:
Feb 14, 2011, 12:58:29 AM (15 years ago)
Author:
ansari
Message:

Amelioration et modifs diverses lors de la redaction du papier, Reza 13/02/2011

Location:
trunk/Cosmo/RadioBeam
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/anapkn.pic

    r3769 r3947  
    2525z = 0.5
    2626# On considere une fraction de HI (HI/Baryon) = 0.02
    27 fHI = 0.020/0.1
     27gHI = 0.02/0.01
    2828
    2929####################################################################
    … …  
    3131defscript convpk2t21 '  --> calcul le coefficient de conversion P(k) en P(k)_temperature mK^2'
    3232  cgeom = (1+$z)*(1+$z)/sqrt($OL+pow((1+$z),3.)*$Om)
    33   cct21 = 0.57*$cgeom*$fHI
     33  cct21 = 0.054*$cgeom*$gHI
    3434  cct21 = $cct21*$cct21
    3535  echo "----  cct21 (NoFactCroiss)= $cct21"
    3636#  On approxime le facteur de croissance lineaire au carre par
    3737#  [(1+0.5)/(1+z)]^2  --> erreur ~ 10%
    38   cct21 = $cct21*pow(1.5/(1+$z),2)
     38#  cct21 = $cct21*pow(1.5/(1+$z),2)
    3939  echo "---- Conv P(k)->Temperature mK^2 : cct21= $cct21"
    4040endscript

  • trunk/Cosmo/RadioBeam/interfconfigs.cc

    r3931 r3947  
    168168
    169169 return vd;
    170 
    171 }
     170}
     171
    172172/* --Fonction -- */
    173173vector<Dish> CreateConfigD(double Ddish, double Eta)
    … …  
    243243
    244244 cout << ">>>CreateConfigD(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
    245 
    246245 return vd;
     246}
     247
     248/* --Fonction -- */
     249vector<Dish> CreateConfigNancay12(double Ddish, double Eta)
     250{
     251  EnumeratedSequence esx,esy;
     252  esx = 0,1,3;
     253  esy = 0,1,2,3;
     254
     255  vector<Dish> vd;
     256  int cnt=0;
     257
     258  for(int ix=0; ix<esx.Size(); ix++) {
     259    for(int iy=0; iy<esy.Size(); iy++) {
     260      cnt++;
     261      vd.push_back(Dish(cnt, ((double)esx.Value(ix))*Ddish,((double)esy.Value(iy))*Ddish,Eta*Ddish));
     262    }
     263  }
     264  cout << ">>>CreateConfigNancay12(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     265  return vd;
     266}
     267
     268/* --Fonction -- */
     269vector<Dish> CreateConfigNancay24(double Ddish, double Eta)
     270{
     271  EnumeratedSequence esx,esy;
     272  esx = 0,1,3,5;
     273  esy = 0,1,2,3,4,5;
     274
     275  vector<Dish> vd;
     276  int cnt=0;
     277
     278  for(int ix=0; ix<esx.Size(); ix++) {
     279    for(int iy=0; iy<esy.Size(); iy++) {
     280      cnt++;
     281      vd.push_back(Dish(cnt, ((double)esx.Value(ix))*Ddish,((double)esy.Value(iy))*Ddish,Eta*Ddish));
     282    }
     283  }
     284  cout << ">>>CreateConfigNancay24(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     285  return vd;
     286}
     287
     288/* --Fonction -- */
     289vector<Dish> CreateConfigNancay36(double Ddish, double Eta)
     290{
     291  EnumeratedSequence esx,esy;
     292  esx = 0,1,2,3,4,5;
     293  esy = 0,1,2,3,4,5;
     294
     295  vector<Dish> vd;
     296  int cnt=0;
     297
     298  for(int ix=0; ix<esx.Size(); ix++) {
     299    for(int iy=0; iy<esy.Size(); iy++) {
     300      cnt++;
     301      vd.push_back(Dish(cnt, ((double)esx.Value(ix))*Ddish,((double)esy.Value(iy))*Ddish,Eta*Ddish));
     302    }
     303  }
     304  cout << ">>>CreateConfigNancay36(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     305  return vd;
     306}
     307
     308/* --Fonction -- */
     309vector<Dish> CreateConfigNancay40(double Ddish, double Eta)
     310{
     311  EnumeratedSequence esx,esy;
     312  esx = 0,1,3,5,7;
     313  esy = 0,1,2,3,4,5,6,7;
     314
     315  vector<Dish> vd;
     316  int cnt=0;
     317
     318  for(int ix=0; ix<esx.Size(); ix++) {
     319    for(int iy=0; iy<esy.Size(); iy++) {
     320      cnt++;
     321      vd.push_back(Dish(cnt, ((double)esx.Value(ix))*Ddish,((double)esy.Value(iy))*Ddish,Eta*Ddish));
     322    }
     323  }
     324  cout << ">>>CreateConfigNancay40(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     325  return vd;
     326}
     327
     328/* --Fonction -- */
     329vector<Dish> CreateConfigNancay128(double Ddish, double Eta)
     330{
     331  EnumeratedSequence esx,esy;
     332  esx = 0,1,3,4,7,10,13,15;
     333  esy = 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
     334
     335  vector<Dish> vd;
     336  int cnt=0;
     337
     338  for(int ix=0; ix<esx.Size(); ix++) {
     339    for(int iy=0; iy<esy.Size(); iy++) {
     340      cnt++;
     341      vd.push_back(Dish(cnt, ((double)esx.Value(ix))*Ddish,((double)esy.Value(iy))*Ddish,Eta*Ddish));
     342    }
     343  }
     344  cout << ">>>CreateConfigNancay128(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     345  return vd;
    247346}
    248347

  • trunk/Cosmo/RadioBeam/interfconfigs.h

    r3931 r3947  
    2626vector<Dish> CreateConfigD(double Ddish=5., double Eta=0.9);
    2727
     28// Configuration a 12 dishs - arrange en cylindres (4+4+4)
     29vector<Dish> CreateConfigNancay12(double Ddish=5., double Eta=0.9);
     30// Configuration a 24 dishs - arrange en cylindres (6+6+0+6+0+6)
     31vector<Dish> CreateConfigNancay24(double Ddish=5., double Eta=0.9);
     32// Configuration a 36 cylindres (6+6+6+6+6+6)
     33vector<Dish> CreateConfigNancay36(double Ddish=5., double Eta=0.9);
     34// Configuration a 40 cylindres (8+8+0+8+0+8+0+8)
     35vector<Dish> CreateConfigNancay40(double Ddish=5., double Eta=0.9);
     36
     37// Configuration a 128 dishs - arrange en 8 cylindres de 15 dishs (CC0CC00C00C00C0C)
     38vector<Dish> CreateConfigNancay128(double Ddish=5., double Eta=0.9);
    2839
    2940// Filled cylinder configuration

  • trunk/Cosmo/RadioBeam/mdish.cc

    r3934 r3947  
    6565}
    6666
     67HProf Four2DResponse::GetProjNoiseLevel(int nbin, bool fgnorm1)
     68{
     69  Histo2D h2w = GetResponse(2*nbin, 2*nbin);
     70  r_8 vmin=h2w.VMin();
     71  r_8 vmax=h2w.VMax();
     72  double seuil=vmax/10000.;
     73  if (seuil<1.e-6) seuil=1.e-6;
     74  r_8 facnorm=((fgnorm1)?vmax:1.);
     75  cout << " Four2DResponse::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax
     76       << " facnorm=" << facnorm << " seuil=" << seuil << endl; 
     77  double kmax=2.*M_PI*D();
     78  HProf hp(0,kmax,nbin);
     79  double x,y;
     80  for(sa_size_t j=0; j<h2w.NBinY(); j++) {
     81    for(sa_size_t i=0; i<h2w.NBinX(); i++) {
     82      h2w.BinCenter(i,j,x,y);
     83      double yw=h2w(i,j);
     84      if (yw<seuil) continue;
     85      hp.Add(sqrt(x*x+y*y),facnorm/yw);
     86    }
     87  }
     88  return hp;
     89}
     90
     91HProf Four2DResponse::GetProjResponse(int nbin, bool fgnorm1)
     92{
     93  Histo2D h2w = GetResponse(2*nbin, 2*nbin);
     94  r_8 vmin=h2w.VMin();
     95  r_8 vmax=h2w.VMax();
     96  r_8 facnorm=((fgnorm1)?(1./vmax):1.);
     97  cout << " Four2DResponse::GetProjResponse Min,Max=" << vmin << " , " << vmax
     98       << " facnorm=" << facnorm << endl;
     99  double kmax=2.*M_PI*D();
     100  HProf hp(0,kmax,nbin);
     101  double x,y;
     102  for(sa_size_t j=0; j<h2w.NBinY(); j++) {
     103    for(sa_size_t i=0; i<h2w.NBinX(); i++) {
     104      h2w.BinCenter(i,j,x,y);
     105      hp.Add(sqrt(x*x+y*y),h2w(i,j)*facnorm);
     106    }
     107  }
     108  return hp;
     109}
     110
    67111//---------------------------------------------------------------
    68112// -- Four2DRespTable : Reponse tabulee instrumentale ds le plan k_x,k_y (angles theta,phi)
    … …  
    185229  SetBeamNSamples();
    186230  SetPrtLevel();
     231  fgcomputedone_=false;
    187232  mcnt_=0;
    188233}
    189234
    190 Histo2D MultiDish::GetResponse()
    191 {
    192   cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
     235void MultiDish::ComputeResponse()
     236{
     237  cout << " MultiDish::ComputeResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
    193238  double kmx = 1.2*DeuxPI*dmax_/lambda_;
    194239  double dkmx = kmx/(double)nx_;
    … …  
    198243  h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
    199244  h2w_.SetZeroBin(0.,0.);
     245  kmax_=kmx;
    200246
    201247  double dold = dishes_[0].Diameter()/lambda_;
    … …  
    210256  double dtet = thetamax_/(double)ntet_;
    211257  double dphi = phimax_/(double)nphi_;
    212   cout << " MultiDish::GetResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_
     258  cout << " MultiDish::ComputeResponse() - ThetaMax=" << thetamax_ << " NTheta=" << ntet_
    213259       << " PhiMax=" <<  phimax_ << " NPhi=" << nphi_ << endl;
    214260
    … …  
    224270      sumw += CumulResp(rd, theta, -(phi+M_PI));
    225271    }
    226     if (prtlev_>0) cout << "  MultiDish::GetResponse() done ktheta=" << kt << " / MaxNTheta= " << ntet_ << endl;
    227   }
     272    if (prtlev_>0)
     273      cout << "  MultiDish::ComputeResponse() done ktheta=" << kt << " / MaxNTheta= "
     274           << ntet_ << endl;
     275  }
     276  r_8 rmin,rmax;
     277  h2w_.GetMinMax(rmin,rmax);
     278  cout << "  MultiDish::ComputeResponse() finished : Rep_min,max=" << rmin << "," << rmax << " sumW0="
     279       << sumw << " ?=" << h2w_.SumWBinZero() << endl;
     280  fgcomputedone_=true;
     281  return;
     282}
     283
     284Histo2D MultiDish::GetResponse()
     285{
     286  if (!fgcomputedone_) ComputeResponse();
     287
    228288  double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
    229289  double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
    … …  
    236296    ib=jb=0;
    237297  }
     298  double sumw=h2w_.SumWBinZero();
    238299  double vmax=h2.VMax();
    239300  cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax 
    … …  
    257318  return h2;
    258319}
     320
     321HProf MultiDish::GetProjNoiseLevel(int nbin, bool fgnorm1)
     322{
     323  r_8 vmin,vmax;
     324  h2w_.GetMinMax(vmin,vmax);
     325  double seuil=vmax/10000.;
     326  if (seuil<1.e-6) seuil=1.e-6;
     327  r_8 facnorm=((fgnorm1)?vmax:1.);
     328  cout << " MultiDish::GetProjNoiseLevel Min,Max=" << vmin << " , " << vmax
     329       << " facnorm=" << facnorm << " seuil=" << seuil << endl; 
     330  HProf hp(0,kmax_,nbin);
     331  for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
     332    double y=h2w_.Y(j);
     333    for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
     334      double x=h2w_.X(i);
     335      double yw=h2w_(i,j);
     336      if (yw<seuil) continue;
     337      hp.Add(sqrt(x*x+y*y),facnorm/yw);
     338    }
     339  }
     340  return hp;
     341}
     342
     343HProf MultiDish::GetProjResponse(int nbin, bool fgnorm1)
     344{
     345  r_8 vmin,vmax;
     346  h2w_.GetMinMax(vmin,vmax);
     347  r_8 facnorm=((fgnorm1)?(1./vmax):1.);
     348  cout << " MultiDish::GetProjResponse Min,Max=" << vmin << " , " << vmax
     349       << " facnorm=" << facnorm << endl; 
     350  HProf hp(0,kmax_,nbin);
     351  for(sa_size_t j=0; j<h2w_.NBinY(); j++) {
     352    double y=h2w_.Y(j);
     353    for(sa_size_t i=0; i<h2w_.NBinX(); i++) {
     354      double x=h2w_.X(i);
     355      hp.Add(sqrt(x*x+y*y),h2w_(i,j)*facnorm);
     356    }
     357  }
     358  return hp;
     359}
     360
    259361
    260362Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
    … …  
    312414      double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
    313415      double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
     416      double pgain=dishes_[i].Gain()*dishes_[j].Gain();
    314417      rot.Do(kx0, ky0);
    315418      //    if (kx0<0) kx0=-kx0;
    … …  
    321424          double y = jy*dy;
    322425          if ((ix>0)&&(jy>0)) {
    323             sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
    324             if (j!=i) sumw += AddToHisto(-kx0, -ky0, x, y, rd(x,y), fgfh);
     426            sumw += AddToHisto(kx0, ky0, x, y, rd(x,y)*pgain, fgfh);
     427            if (j!=i) sumw += AddToHisto(-kx0, -ky0, x, y, rd(x,y)*pgain, fgfh);
    325428          }
    326429          else {
    327430            if ((ix==0)&&(jy==0)) {
    328               sumw += h2w_.Add(kx0, ky0, rd(0.,0.), fgfh);
    329               if (j!=i) sumw += h2w_.Add(-kx0, -ky0, rd(0.,0.), fgfh);
     431              sumw += h2w_.Add(kx0, ky0, rd(0.,0.)*pgain, fgfh);
     432              if (j!=i) sumw += h2w_.Add(-kx0, -ky0, rd(0.,0.)*pgain, fgfh);
    330433            }
    331434            else {
    332               double w = rd(x,y);
     435              double w = rd(x,y)*pgain;
    333436              if (ix==0) {
    334437                sumw += h2w_.Add(kx0, ky0+y, w, fgfh);

  • trunk/Cosmo/RadioBeam/mdish.h

    r3933 r3947  
    4242  inline void setLambda(double lambda=1.)
    4343  { lambda_ = lambda;   lambda_ratio_ = lambda_/lambdaref_; }   
     44
     45  inline double getLambdaRef()   { return lambdaref_; }
     46  inline double getLambda()  { return lambda_; }
    4447 
    4548  // Return the 2D response for wave vector (kx,ky)
    … …  
    4851  { return Value(kx, ky); }
    4952  virtual Histo2D GetResponse(int nx=255, int ny=255); 
     53  // Retourne le niveau moyen du bruit projete 1D en fonction (sqrt(u^2+v^2)
     54  HProf GetProjNoiseLevel(int nbin=128, bool fgnorm1=true);
     55  // Retourne la reponse moyenne projetee 1D en fonction (sqrt(u^2+v^2)
     56  HProf GetProjResponse(int nbin=128, bool fgnorm1=true);
     57
    5058  inline double D() { return dx_; } ;
    5159  inline double Dx() { return dx_; } ;
    … …  
    95103  // Circular dish
    96104  Dish(int id, double x, double y, double diam)
    97     :  id_(id), X(x), Y(y), D(diam), Dx(D), Dy(D), fgcirc_(true)   {   }
     105    :  id_(id), X(x), Y(y), D(diam), Dx(D), Dy(D), fgcirc_(true), gain_(1.)   {   }
    98106  // Receiver with rectangular type answer in kx,ky plane
    99107  Dish(int id, double x, double y, double dx, double dy)
    100     :  id_(id), X(x), Y(y), D(sqrt(dx*dy)), Dx(dx), Dy(dy), fgcirc_(false)   {   }
     108    :  id_(id), X(x), Y(y), D(sqrt(dx*dy)), Dx(dx), Dy(dy), fgcirc_(false), gain_(1.)    {   }
    101109  Dish(Dish const& a)
    102     :  id_(a.id_), X(a.X), Y(a.Y), D(a.D), Dx(a.Dx), Dy(a.Dy), fgcirc_(a.fgcirc_)     {   }
     110    :  id_(a.id_), X(a.X), Y(a.Y), D(a.D), Dx(a.Dx), Dy(a.Dy), fgcirc_(a.fgcirc_), gain_(a.gain_)  {   }
     111  inline void setGain(double gain) { gain_=gain; return; }
    103112  inline bool isCircular() { return fgcirc_; }
    104113  inline int ReflectorId() { return id_; }
    … …  
    106115  inline double DiameterX() { return Dx; }
    107116  inline double DiameterY() { return Dy; }
     117  inline double Gain() { return gain_; }
    108118
    109119  int id_;   // numero de reflecteur
    … …  
    111121  double Dx, Dy;
    112122  bool fgcirc_;  // false -> rectangular dish
     123  double gain_;
    113124};
    114125
    … …  
    132143  { beamnx_=nx; beamny_=ny; }
    133144
     145  // Calcul la reponse ds le plan 2D  (u,v) = (kx,ky)
     146  void ComputeResponse();
     147  // Retourne la reponse 2D ds le plan (u,v) = (kx,ky) sous forme d'histo 2D
    134148  Histo2D GetResponse();
     149
     150  // Retourne le niveau moyen du bruit projete 1D en fonction (sqrt(u^2+v^2)
     151  HProf GetProjNoiseLevel(int nbin=128, bool fgnorm1=true);
     152  // Retourne la reponse moyenne projetee 1D en fonction (sqrt(u^2+v^2)
     153  HProf GetProjResponse(int nbin=128, bool fgnorm1=true);
    135154
    136155  double CumulResp(Four2DResponse& rd, double theta=0., double phi=0.);
    … …  
    143162  double AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh);
    144163
    145   double lambda_, dmax_;
     164  double lambda_, dmax_, kmax_;
    146165  vector<Dish> dishes_;
    147166  bool fgnoauto_;
    … …  
    153172  //   Histo2D h2w_, h2cnt_;
    154173  QHis2D h2w_;
     174  bool fgcomputedone_;
    155175  int mcnt_;
    156176  int prtlev_,prtmodulo_;

  • trunk/Cosmo/RadioBeam/pknoise.cc

    r3931 r3947  
    4444       << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
    4545       << "    -scut SCutValue (default=100.) \n"
    46        << "    -ngen NGen (default=1) number of noise fourier amp generations \n"
     46       << "    -ngen NGen (default=0) number of noise fourier amp generations \n"
     47       << "       NGen=0 -> Call ComputeNoisePk(), else generate Fourier Amplitudes (random) \n"
    4748       << "    -z redshift (default=0.7) \n"
    4849       << "    -prt PrtLev,PrtModulo (default=0,10) " << endl;
    … …  
    7475  bool fgrenorm=false;
    7576  double rmax=1.;
    76   int NMAX = 1;
     77  int NMAX = 0;
    7778  double SCut=0.;
    7879  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m 
    … …  
    141142      cout << "pknoise[1]: initializing Four2DRespTable from file" << resptblname << endl;
    142143      resptbl.readFromPPF(resptblname);
     144      cout << "pknoise[1.b] Four2DRespTable.LambdaRef=" << resptbl.getLambdaRef() << " setLambda("
     145           << lambda << ")" << endl;
     146      resptbl.setLambda(lambda);
    143147      arep_p=&resptbl;
    144148      if (fgrenorm) {
    145         cout << " pknoise[1.b] call to resptbl.renormalize(" << rmax << ")";
     149        cout << "pknoise[1.c] call to resptbl.renormalize(" << rmax << ")";
    146150        double omax=resptbl.renormalize(rmax);
    147151        cout << " ... Old Max=" << omax << endl;
    … …  
    159163
    160164    cout << " pknoise[3]: Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    161     PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
    162     pkn.SetPrtLevel(prtlev,prtmod);
    163     HProf hpn = pkn.Compute();
    164     po << hpn;
    165  
     165    if (NMAX>0) {
     166      PkNoiseCalculator pkn(m3d, *(arep_p), SCut, NMAX, tits.c_str());
     167      pkn.SetPrtLevel(prtlev,prtmod);
     168      HProf hpn = pkn.Compute();
     169      po << hpn;
     170    }
     171    else {
     172      Histo fracmodok;
     173      DataTable dtnoise;
     174      HProf hpn = m3d.ComputeNoisePk(*(arep_p),fracmodok,dtnoise,SCut);
     175      po << hpn;
     176      string outfile2 = "x"+outfile;
     177      POutPersist po2(outfile2);
     178      HProf h1dnoise=arep_p->GetProjNoiseLevel();
     179      HProf h1drep=arep_p->GetProjResponse();
     180      po2 << PPFNameTag("dtnoise") << dtnoise;
     181      po2 << PPFNameTag("hpnoise") << hpn;
     182      po2 << PPFNameTag("fracmodok") << fracmodok;
     183      po2 << PPFNameTag("h1dnoise") << h1dnoise;
     184      po2 << PPFNameTag("h1drep") << h1drep;
     185    }
    166186    rc = 0;
    167187  }  // End of try bloc

  • trunk/Cosmo/RadioBeam/plpkn.pic

    r3769 r3947  
    1717setup5
    1818scalewz 0.7 2500 100 1
     19# scalewz 0.25 989.13 80.26 1
    1920y1 = 500*$cct21
    2021y2 = 9e4*$cct21

  • trunk/Cosmo/RadioBeam/qhist.h

    r3783 r3947  
    4343  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
    4444  void SetZeroBin(r_8 x=0, r_8 y=0);
    45   inline double WBinX() { return dxb; }
    46   inline double WBinY() { return dyb; }
     45  inline r_8 WBinX() { return dxb; }
     46  inline r_8 WBinY() { return dyb; }
     47  inline sa_size_t NBinX() { return nx; }
     48  inline sa_size_t NBinY() { return ny; }
     49  inline r_8 X(sa_size_t i) { return (dxb*((double)i+0.5)+xmin); }
     50  inline r_8 Y(sa_size_t j) { return (dyb*((double)j+0.5)+ymin); }
     51  inline r_8 operator()(sa_size_t i, sa_size_t j) { return aw(i,j); }
     52
     53  inline void GetMinMax(r_8& min, r_8& max) { return aw.MinMax(min,max); }
     54
     55  inline double SumWBinZero() { return sumw0; }
     56  void Renormalize(r_8 mfac=1.) { aw*=mfac; }
     57
    4758  Histo2D Convert();
    4859 
     60protected:
    4961  r_8 xmin,xmax,ymin,ymax;
    5062  r_8 dxb,dyb;

  • trunk/Cosmo/RadioBeam/radutil.cc

    r3829 r3947  
    4545  double zz = 1.+z_;
    4646  double cc=zz*zz/sqrt(OmegaMatter_*zz*zz*zz+OmegaLambda_);
    47   cc *= ((h100_/0.7)*(OmegaBaryons_/0.04)*(fracHI_/0.1));
    48   return (cc*0.57);
     47  cc *= ((h100_/0.7)*(OmegaBaryons_/0.044)*(fracHI_/0.01));
     48  return (cc*0.054);
    4949}

  • trunk/Cosmo/RadioBeam/radutil.h

    r3829 r3947  
    3737  // Definition des parametres cosmologiques utiles pour le calcul de la temperature d'emission a 21 cm
    3838  // retourne la valeur de OmegaLambda (univers plat)
    39   double setCosmoParam(double omegamatter=0.02581, double omegabaryon=0.0441, double h100=0.719, double fracHI=0.02);
     39  double setCosmoParam(double omegamatter=0.2581, double omegabaryon=0.0441, double h100=0.719, double fracHI=0.02);
    4040  inline void setFracHI(double fracHI=0.02) { fracHI_=fracHI; }
    4141

  • trunk/Cosmo/RadioBeam/repicon.cc

    r3936 r3947  
    3232
    3333// pour sauver la reponse mdresp et la config des dishes dans un fichier PPF
    34 void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp);
     34void SaveDTVecDishesH2Resp(POutPersist& po, vector<Dish>& vdishes, Four2DRespTable& mdresp);
    3535
    3636// ---------------------------------------------------------------------
    37 // main program for computing interferometer response un (u,v) plane
     37// main program for computing interferometer response in (u,v) plane
    3838// R. Ansari  - Avril-Juin 2010
    3939// ---------------------------------------------------------------------
    … …  
    4242{
    4343  cout << " Usage: repicon [-parname Value] configId OutPPFName \n"
    44        << " configIds: f4x4,f8x8,f20x20, confA,confB,confC,confD, hex12,cross11, f3cyl,f8cyl,f3cylp,f8cylp \n"
    45        << "   f4x4 , f8x8 , f20x20 Filled array of nxn dishes  \n"
     44       << " configIds: f4x4,f8x8,f11x11,f20x20, confA,confB,confC,confD, hex12,cross11, \n"
     45       <<  "           f3cyl,f8cyl,f3cylp,f8cylp, nan12,nan24,nan36,nan40,nan128 \n"
     46       << "   f4x4 , f8x8 , f11x11 , f20x20 Filled array of nxn dishes  \n"
    4647       << "   confA , confB, confC, confD : semi-filled array of dishes \n"
    4748       << "   hex12,cross11 : ASKAP like double hexagonal (12xD=12m), cross config (11xD=12m) \n"
     49       << "   nan12,nan24,nan36,nan40,nan128: 3,4,6,5,8 cylinder like configurations with 12,24,36,40,128 dishes\n"
    4850       << "   f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders \n"
    4951       << "   f4cylw : filled array of 4 perfect of wide cylinders \n"
    50        << "  [ -parname value] : -renmax -z -prt -D -lmax \n"
     52       << "   pit2cyl,pit2cylw : 2 cylinders 15mx7m, 20mx10 for z=0.3 at Pittsburg (perfect) \n"
     53       << "   pit4cyl : 4 cylinders 32mx8m for z=0.3-0.7 at Pittsburg (perfect) \n"
     54       << "  [ -parname value] : -renmax -z -prt -D -lmax -eta -autocor \n"
    5155       << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
    5256       << "    -z redshift (default=0.7) --> determines Lambda \n"
    5357       << "    -D DishDiameter (default=5 m) \n"
     58       << "    -eta fill_factor (default=0.90) \n"
    5459       << "    -lmax array extension (default=100 m ) for response calculation kmax \n"
    5560       << "    -repnxy Nx,Ny (default=200,200) ResponseTable binning \n"
    5661       << "    -beamnxy Nx,Ny (default=200,200) Beam sampling \n"
    5762       << "    -rot ThetaMaxDeg,NTheta,PhiMaxDeg,NPhi (default NO-rotate/pointing -> 23,10,30,15 ) \n"
     63       << "    -autocor : keep antenna auto-correlation signals (default NO) \n"
    5864       << "    -prt PrtLev,PrtModulo (default=0,10) \n"
    5965       << endl;
    … …  
    8288 
    8389  double Ddish=5.;
     90  double Eta=0.9;
    8491  bool fgDfixed=false;
    8592  double LMAX = 100.;  // taille de la zone, pour calcul kmax
    … …  
    94101  int NBX=200;
    95102  int NBY=200;
     103  bool fgnoauto=true;  // do not keep antenna autocorrelation signal
    96104 
    97105  int ka=1;
    … …  
    106114      Ddish=atof(arg[ka+1]);  fgDfixed=true;  ka+=2;
    107115    }
     116    else if (strcmp(arg[ka],"-eta")==0) {
     117      Eta=atof(arg[ka+1]);  ka+=2;
     118    }
    108119    else if (strcmp(arg[ka],"-lmax")==0) {
    109120      LMAX=atof(arg[ka+1]);  fgLMAXfixed=true;  ka+=2;
    … …  
    117128    else if (strcmp(arg[ka],"-rot")==0) {
    118129      sscanf(arg[ka+1],"%lg,%d,%lg,%d",&thetamxdeg,&nteta,&phimxdeg,&nphi);  fgpoint=true;  ka+=2;
     130    }
     131    else if (strcmp(arg[ka],"-autocor")==0) {
     132      fgnoauto=false;   ka+=1;
    119133    }
    120134    else if (strcmp(arg[ka],"-prt")==0) {
    … …  
    143157    LAMBDA = conv.getLambda();
    144158   
    145     double Eta=0.95;
    146159    double cylW=12.;   // Largeur des cylindres
    147160    double cylRL=2*LAMBDA;  // Longeur des elements de reception le long du cylindre
    … …  
    150163
    151164    vector<Dish> vdishes;
    152     cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config << endl;
     165    cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config
     166         << " Ddish=" << Ddish << " m. Eta=" << Eta << endl;
    153167
    154168    if (config=="f4x4") {
    … …  
    157171    else if (config=="f8x8") {
    158172      vdishes=CreateFilledSqConfig(8,Ddish, Eta);
     173    }
     174    else if (config=="f11x11") {
     175      vdishes=CreateFilledSqConfig(11,Ddish, Eta);
    159176    }
    160177    else if (config=="f20x20") {
    … …  
    172189      cylW=25.;   cylRL=3*LAMBDA;
    173190      vdishes=CreateFilledCylConfig(4, 100, cylW, cylRL, etaW, etaRL, false);
     191    }
     192    else if (config=="pit2cyl") {
     193      etaW=0.9; etaRL=0.9;
     194      cylW=7.;   cylRL=2*LAMBDA;
     195      vdishes=CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
     196    }
     197    else if (config=="pit2cylw") {
     198      etaW=0.9; etaRL=0.9;
     199      cylW=9.;   cylRL=2*LAMBDA;
     200      vdishes=CreateFilledCylConfig(2, 40, cylW, cylRL, etaW, etaRL, false);
     201    }
     202    else if (config=="pit4cyl") {
     203      etaW=0.9; etaRL=0.9;
     204      cylW=8.;   cylRL=2*LAMBDA;
     205      vdishes=CreateFilledCylConfig(4, 64, cylW, cylRL, etaW, etaRL, false);
    174206    }
    175207    else if (config=="f8cyl") {
    … …  
    206238      vdishes=CreateDoubleHexagonConfig(Ddish);
    207239    }
     240    else if (config=="nan12") {
     241      vdishes=CreateConfigNancay12(Ddish, Eta);
     242    }
     243    else if (config=="nan24") {
     244      vdishes=CreateConfigNancay24(Ddish, Eta);
     245    }
     246    else if (config=="nan36") {
     247      vdishes=CreateConfigNancay36(Ddish, Eta);
     248    }
     249    else if (config=="nan40") {
     250      vdishes=CreateConfigNancay40(Ddish, Eta);
     251    }
     252    else if (config=="nan128") {
     253      vdishes=CreateConfigNancay128(Ddish, Eta);
     254    }
    208255    else {
    209256      cout << " NON valid configuration option -> exit" << endl;
    … …  
    212259       
    213260    double Dol = LMAX/LAMBDA;
    214     bool fgnoauto = true;
    215261   
    216262    cout << " repicon[1] : Lambda=" << LAMBDA << " LMAX= " << LMAX << " NDishes=" << vdishes.size()
    217          << " D-Dish=" << Ddish << " m." << endl;
     263         << " D-Dish=" << Ddish << " m. " << ((fgnoauto)?" NO-":"With-") << "AutoCorrelation" << endl;
    218264 
    219265    MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
    … …  
    230276    Histo2D hrep = mdish.GetResponse();
    231277    PrtTim("Apres mdish.GetResponse()");
     278    HProf h1dnoise = mdish.GetProjNoiseLevel();
     279    HProf h1drep = mdish.GetProjResponse();
    232280
    233281    Four2DRespTable mdresp(hrep, Dol, LAMBDA);
    … …  
    241289
    242290    string outfile2 = "hdt_"+outfile;
    243     cout << " repicon[4] : saving H2D-response, multidish config to PPF file " << outfile2 << endl;
    244     SaveDTVecDishesH2Resp(outfile2, vdishes, mdresp);
     291    POutPersist po2(outfile2);
     292    cout << " repicon[4] : saving H1D/H2D-response, multidish config to PPF file " << outfile2 << endl;
     293    po2 << PPFNameTag("h1dnoise") << h1dnoise;
     294    po2 << PPFNameTag("h1drep") << h1drep;
     295    SaveDTVecDishesH2Resp(po2, vdishes, mdresp);
    245296
    246297    rc = 0;
    … …  
    265316
    266317/*-- Nouvelle-Fonction --*/
    267 void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp)
     318void SaveDTVecDishesH2Resp(POutPersist& po, vector<Dish>& vdishes, Four2DRespTable& mdresp)
    268319{
    269320  char* names[5]={"did","posx","posy","diam","diamy"};
    … …  
    285336  }
    286337  Histo2D h2rep=mdresp.GetResponse();
    287   POutPersist po(outfile);
    288338  po << PPFNameTag("mdish") << ntvd;
    289339  po << PPFNameTag("h2rep") << h2rep;

  • trunk/Cosmo/RadioBeam/specpk.cc

    r3931 r3947  
    266266}
    267267
     268// Compute noise power spectrum as a function of wave number k
     269// No random generation, computes profile of noise sigma^2
     270// cells with amp^2=re^2+im^2>s2cut are ignored
     271// Output : noise power spectrum (profile histogram)
     272HProf Four3DPk::ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt,
     273                               double s2cut, int nbin, double kmin, double kmax)
     274{
     275  // The second half of the array along Y (matrix rows) contain
     276  // negative frequencies
     277  //  int nbh = sqrt(fourAmp.SizeX()*fourAmp.SizeX()+fourAmp.SizeY()*fourAmp.SizeY()/4.+fourAmp.SizeZ()*fourAmp.SizeY()/4.);
     278  // The profile histogram will contain the mean value of noise sigma^2
     279  // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
     280  //  if (nbin < 1) nbin = nbh/2;
     281  if ((kmax<0.)||(kmax<kmin)) {
     282    kmin=0.;
     283    double maxx=fourAmp.SizeX()*dkx_;
     284    double maxy=fourAmp.SizeY()*dky_/2;
     285    double maxz=fourAmp.SizeZ()*dkz_/2;
     286    kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
     287  }
     288  if (nbin<2) nbin=128;
     289  HProf hp(kmin, kmax, nbin);
     290  hp.SetErrOpt(false);
     291  Histo hmcnt(kmin, kmax, nbin);
     292  Histo hmcntok(kmin, kmax, nbin);
     293
     294  uint_8 nmodeok=0;
     295  // fourAmp represent 3-D fourier transform of a real input array.
     296  // The second half of the array along Y and Z contain negative frequencies
     297  double kxx, kyy, kzz;
     298  // sa_size_t is large integer type 
     299  // We ignore 0th term in all frequency directions ...
     300  for(sa_size_t kz=1; kz<fourAmp.SizeZ(); kz++) {
     301    kzz =  (kz > fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
     302    for(sa_size_t ky=1; ky<fourAmp.SizeY(); ky++) {
     303      kyy =  (ky > fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
     304      for(sa_size_t kx=1; kx<fourAmp.SizeX(); kx++) {  // ignore the 0th coefficient (constant term)
     305        double kxx=(double)kx*dkx_;
     306        double wk = sqrt(kxx*kxx+kyy*kyy+kzz*kzz);
     307        double amp2 = 1./resp(kxx, kyy);
     308        hmcnt.Add(wk);
     309        if ((s2cut>1.e-9)&&(amp2>s2cut))  continue;
     310        hmcntok.Add(wk);
     311        hp.Add(wk, amp2);
     312        nmodeok++;
     313      }
     314    }
     315  }
     316  if ((prtlev_>1)||((prtlev_>0)&&(s2cut>1.e-9))) {
     317    cout << " Four3DPk::ComputeNoisePk/Info : NModeOK=" << nmodeok << " / NMode=" << fourAmp.Size()
     318         << " -> " << 100.*(double)nmodeok/(double)fourAmp.Size() << "%" << endl;
     319  }
     320
     321  fracmodok=hmcntok/hmcnt;
     322
     323  char* nomcol[5] = {"k","pnoise","nmode","nmodok","fracmodok"};
     324  dt.Clear();
     325  dt.AddDoubleColumn(nomcol[0]);   
     326  dt.AddDoubleColumn(nomcol[1]);   
     327  dt.AddIntegerColumn(nomcol[2]);   
     328  dt.AddIntegerColumn(nomcol[3]);   
     329  dt.AddFloatColumn(nomcol[4]);   
     330  DataTableRow  dtr = dt.EmptyRow();
     331  for(int_4 ib=0; ib<hp.NBins(); ib++) {
     332    dtr[0]=hp.BinCenter(ib);
     333    dtr[1]=hp(ib);
     334    dtr[2]=hmcnt(ib);
     335    dtr[3]=hmcntok(ib);
     336    dtr[4]=fracmodok(ib);
     337    dt.AddRow(dtr);
     338  }
     339  return hp;
     340}
     341
    268342//-----------------------------------------------------
    269343// -- MassDist2D class :  2D mass distribution

  • trunk/Cosmo/RadioBeam/specpk.h

    r3930 r3947  
    6464  HProf ComputePk(double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
    6565  void  ComputePkCumul(HProf& hp, double s2cut=0.);
     66
     67  HProf ComputeNoisePk(Four2DResponse& resp, Histo& fracmodok, DataTable& dt,
     68                       double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
    6669
    6770protected:
Note: See TracChangeset for help on using the changeset viewer.