Changeset 3769 in Sophya for trunk/Cosmo


Ignore:
Timestamp:
May 7, 2010, 6:44:43 PM (15 years ago)
Author:
ansari
Message:

Corrections/amelioration du programme de calcul de la sensibilite interfero, Reza 07/05/2010

Location:
trunk/Cosmo/RadioBeam
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/anapkn.pic

    r3756 r3769  
    99
    1010##1 / On ouvre le fichier extrait d'un PPF fait par SimLSS (cmv)
    11 delobjs *
     11# delobjs *
    1212openppf cmvhpkz.ppf
    1313set gdz 1.
    … …  
    392392dopk 'P(k), PNoise - mK^2 Mpc^3, Nancay-MultiBeam  z=0.5'  'Relative Sample Variance, z=0.5'
    393393
    394 endscript
     394### Setup pour plpkn @ z=0.7
     395setup5
     396scalewz 0.7 2500 100 1
     397  y1 = 500*$cct21
     398  y2 = 9e4*$cct21
     399  xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
     400  n/plot hpkz.val*${cct21}%$kk ! ! "same notit nsta connectpoints blueviolet $xyl line=solid,2"
     401
     402endscript
     403

  • trunk/Cosmo/RadioBeam/mdish.cc

    r3756 r3769  
    108108  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
    109109{
     110  ixb0 = jyb0 = 0;
    110111}
    111112QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
    … …  
    123124  sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
    124125  aw.ReSize(2,sz);
     126  SetZeroBin();
    125127  sumw0=0.;
    126128  return;
    … …  
    131133  sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
    132134  if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
    133   double rw = ((ix==0)&&(jy==0)) ? w : 0.;
     135  double rw = ((ix==ixb0)&&(jy==jyb0)) ? w : 0.;
    134136  sumw0 += rw;
    135137  if (fgfh) aw(ix,jy) += w;
    136138  return rw;
    137139}
     140void QHis2D::SetZeroBin(r_8 x, r_8 y)
     141{
     142  ixb0 = (sa_size_t)((x-xmin)/dxb);
     143  jyb0 = (sa_size_t)((y-ymin)/dyb);
     144}
    138145Histo2D QHis2D::Convert()
    139146{
     147  int_4 imn,jmn,imx,jmx;
     148  r_8 min = aw(0,0);
     149  r_8 max = aw(0,0);
     150  imn=jmn=imx=jmx=0;
    140151  Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
    141152  for(int_4 j=0; j<ny; j++)
    142     for(int_4 i=0; i<nx; i++)  h2(i,j) = aw(i,j);
     153    for(int_4 i=0; i<nx; i++) {
     154      h2(i,j) = aw(i,j);
     155      if (aw(i,j)>max) {
     156        imx=i;  jmx=j;  max=aw(i,j);
     157      }
     158      if (aw(i,j)<min) {
     159        imn=i;  jmn=j;  min=aw(i,j);
     160      }
     161    }
     162  cout << "QHis2D::Convert()/Info: Nx,Ny=" << nx << "," << ny << " SumW=" << sumw0
     163       << "\n ... Max:" << imx << "," << jmx << " ->" << max
     164       << " @" << imx*dxb+xmin << "," << jmx*dyb+ymin
     165       << "\n ...Min:" << imn << "," << jmn << " ->" << min
     166       << " @" << imn*dxb+xmin << "," << jmn*dyb+ymin << endl;
    143167  return h2;
    144168}
    … …  
    159183  cout << " MultiDish::GetResponse() - NDishes=" << dishes_.size() << " nx=" << nx_ << " ny=" << ny_ << endl;
    160184  double kmx = 1.2*DeuxPI*dmax_/lambda_;
    161   /*
    162   h2w_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
    163   h2cnt_= Histo2D(0.,kmx,nx_,0.,kmx,ny_);
    164   h2w_.Zero();
    165   h2cnt_.Zero();
    166   */
    167   h2w_.Define(0.,kmx,nx_,0.,kmx,ny_);
     185  double dkmx = kmx/(double)nx_;
     186  double dkmy = kmx/(double)ny_;
     187  double kmxx = ((double)nx_+0.5)*dkmx;
     188  double kmxy = ((double)ny_+0.5)*dkmy;
     189  h2w_.Define(-kmxx,kmxx,2*nx_+1,-kmxy,kmxy,2*ny_+1);
     190  h2w_.SetZeroBin(0.,0.);
    168191
    169192  double dold = dishes_[0].D/lambda_;
    … …  
    184207      sumw += CumulResp(rd, (double)kt*dtet, (double)jp*dphi);
    185208
    186   double kx0 = DeuxPI*fabs(dishes_[1].X-dishes_[0].X)/lambda_;
    187   double ky0 = DeuxPI*fabs(dishes_[1].Y-dishes_[0].Y)/lambda_;
    188   int_4 ib, jb;
     209  double kx1 = DeuxPI*(dishes_[0].DiameterX())/lambda_;
     210  double ky1 = DeuxPI*(dishes_[0].DiameterY())/lambda_;
     211  int_4 ib,jb;
    189212  //  h2w_ /= h2cnt_;
    190213  Histo2D h2 = h2w_.Convert();
    191   h2.FindBin(kx0, ky0, ib, jb);
    192   cout << " ---- MultiDish::GetResponse() VMin=" << h2.VMin() << " VMax= " << h2.VMax()
    193        << " h(0,0)=" << h2(0,0) << " h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
     214  h2.FindBin(kx1, ky1, ib, jb);
     215  if ((kx1<0)||(ky1<0)||(kx1>=h2.NBinX())||(ky1>=h2.NBinY())) {
     216    cout << " MultiDish::GetResponse[1]/ERROR kx1,ky1=" << kx1 <<","<< ky1 << " --> ib,jb=" << ib <<","<< jb << endl;
     217    ib=jb=0;
     218  }
     219  double vmax=h2.VMax();
     220  cout << " MultiDish::GetResponse[1] VMin=" << h2.VMin() << " VMax= " << vmax 
     221       << " h(0,0)=" << h2(0,0) << " kx1,ky1->h(" << ib <<"," << jb << ")=" << h2(ib,jb) <<endl;
    194222  //  double fnorm=sqrt((double)dishes_.size())/h2.VMax();
    195223  double fnorm=1.;
    196   if (h2.VMax() > sumw) {
     224  if (vmax > sumw) {
    197225    fnorm=(double)dishes_.size()/h2.VMax();
    198     cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " sumw=" << sumw
    199          << " Renormalizing x NDishes/sumw  " << fnorm << endl;
     226    cout << " MultiDish::GetResponse[2]/Warning h2.VMax()=" << vmax << " >  sumw=" << sumw << endl; 
     227    cout << "   ... NDishes=" << dishes_.size() << " sumw=" << sumw
     228         << " Renormalizing x NDishes/VMax  " << fnorm << endl;
    200229  }
    201230  else {
    202     fnorm=(double)dishes_.size()/h2.VMax();
    203     cout << " ---- MultiDish::GetResponse() NDishes=" << dishes_.size() << " VMax=" << h2.VMax() 
    204          << " Renormalizing x NDishes/h2.VMax()   " << fnorm << endl;
     231    fnorm=(double)dishes_.size()/sumw;
     232    cout << " MultiDish::GetResponse[3] NDishes=" << dishes_.size() << " sumw=" << sumw 
     233         << " Renormalizing x NDishes/sumw   " << fnorm << endl;
    205234  }
    206235  h2 *= fnorm;
    207   cout << " ---- MultiDish::GetResponse() APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)="
     236  cout << " ---- MultiDish::GetResponse/[4] APRES VMin=" << h2.VMin() << " VMax= " << h2.VMax() << " h(0,0)="
    208237       << h2(0,0) << endl;
    209238  return h2;
    210239}
    211240
    212 /*
    213 double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
    214 {
    215   double xxp = kx0+x;
    216   double yyp = ky0+y;
    217   double sumw=0.;
    218   int_4 ib, jb;
    219   h2w_.FindBin(xxp, yyp, ib, jb);
    220   if ((ib==0)&&(jb==0))  sumw+=w;
    221   if (fgfh) {
    222     h2w_.Add(xxp, yyp, w);
    223     h2cnt_.Add(xxp, yyp, 1.);
    224   }
    225   double xxm=kx0-x;
    226   double yym=ky0-y;
    227   if (xxm>0.)  {
    228     h2w_.FindBin(xxm, yyp, ib, jb);
    229     if ((ib==0)&&(jb==0))  sumw+=w;
    230     if (fgfh) {
    231       h2w_.Add(xxm, yyp, w);
    232       h2cnt_.Add(xxm, yyp, 1.);
    233     }
    234     if (yym>0.) {
    235       h2w_.FindBin(xxm, yym, ib, jb);
    236       if ((ib==0)&&(jb==0))  sumw+=w;
    237       if (fgfh) {
    238         h2w_.Add(xxm, yym, w);
    239         h2cnt_.Add(xxm, yym, 1.);
    240       }
    241     }
    242   }
    243   if (yym>0.)  {
    244     h2w_.FindBin(xxp, yym, ib, jb);
    245     if ((ib==0)&&(jb==0))  sumw+=w;
    246     if (fgfh) {
    247       h2w_.Add(xxp, yym, w);
    248       h2cnt_.Add(xxp, yym, 1.);
    249     }
    250   }
    251   return sumw;
    252 }
    253 */
     241Histo2D MultiDish::PosDist(int nx, int ny, double dmax)
     242{
     243  if (dmax<1e-3)  dmax=nx*dishes_[0].Diameter();
     244  double dd = dmax/nx/2.;
     245  Histo2D hpos(-dd,dmax+dd,nx+1,-dd,dmax+dd,ny+1);
     246  for(size_t i=0; i<NbDishes(); i++) {
     247    hpos.Add(dishes_[i].X, dishes_[i].Y);
     248  }
     249  return hpos;
     250}
    254251
    255252double MultiDish::AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh)
    … …  
    261258  double xxm=kx0-x;
    262259  double yym=ky0-y;
    263   if (xxm>0.)  {
    264     sumw += h2w_.Add(xxm, yyp, w, fgfh);
    265     if (yym>0.)  sumw += h2w_.Add(xxm, yym, w, fgfh);
    266   }
    267   if (yym>0.)  sumw += h2w_.Add(xxp, yym, w, fgfh);
     260  //  if (xxm>0.)  {
     261  sumw += h2w_.Add(xxm, yyp, w, fgfh);
     262  // if (yym>0.) 
     263  sumw += h2w_.Add(xxm, yym, w, fgfh);
     264  //  }
     265  // if (yym>0.) 
     266  sumw += h2w_.Add(xxp, yym, w, fgfh);
    268267  return sumw;
    269268}
    … …  
    275274  double dx = h2w_.WBinX()/5;
    276275  double dy = h2w_.WBinY()/5;
    277   int nbx = DeuxPI*rd.Dx()/dx;
    278   int nby = DeuxPI*rd.Dy()/dy;
     276  int nbx = DeuxPI*rd.Dx()/dx+1;
     277  int nby = DeuxPI*rd.Dy()/dy+1;
    279278  dx = DeuxPI*rd.Dx()/(double)nbx;
    280279  dy = DeuxPI*rd.Dy()/(double)nby;
    … …  
    287286
    288287  for(size_t i=0; i<dishes_.size(); i++) {
    289     for(size_t j=i; j<dishes_.size(); j++) {
    290       double kx0 = DeuxPI*fabs(dishes_[i].X-dishes_[j].X)/lambda_;
    291       double ky0 = DeuxPI*fabs(dishes_[i].Y-dishes_[j].Y)/lambda_;
     288    for(size_t j=0; j<dishes_.size(); j++) {
     289      double kx0 = DeuxPI*(dishes_[i].X-dishes_[j].X)/lambda_;
     290      double ky0 = DeuxPI*(dishes_[i].Y-dishes_[j].Y)/lambda_;
    292291      rot.Do(kx0, ky0);
    293       if (kx0<0) kx0=-kx0;
    294       if (ky0<0) ky0=-ky0;
     292      //    if (kx0<0) kx0=-kx0;
     293      //    if (ky0<0) ky0=-ky0;
    295294      bool fgfh= (!fgnoauto_ || (dishes_[i].ReflectorId()!=dishes_[j].ReflectorId()));
    296295      for(int ix=0; ix<nbx; ix++)
    297296        for(int jy=0; jy<nby; jy++) {
    298297          double x = ix*dx; 
    299           double y = jy*dy; 
    300           sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
     298          double y = jy*dy;
     299          if ((ix>0)&&(jy>0)) {
     300            sumw += AddToHisto(kx0, ky0, x, y, rd(x,y), fgfh);
     301          }
     302          else {
     303            if ((ix==0)&&(jy==0))
     304              sumw += h2w_.Add(kx0, ky0, rd(0.,0.), fgfh);
     305            else {
     306              double w = rd(x,y);
     307              if (ix==0) {
     308                sumw += h2w_.Add(kx0, ky0+y, w, fgfh);
     309                sumw += h2w_.Add(kx0, ky0-y, w, fgfh);
     310              }
     311              else {
     312                sumw += h2w_.Add(kx0+x, ky0, w, fgfh);
     313                sumw += h2w_.Add(kx0-x, ky0, w, fgfh);
     314              }
     315            }
     316            //   
     317          }
    301318        }
    302     }
    303319    //    if (i%10==0)
    304320    //      cout << " MultiDish::CumulResp() done i=" << i << " / imax=" << dishes_.size()
    305321    //     << " theta=" << theta << " phi=" << phi << endl;
     322    }
    306323  }
    307324  return sumw;

  • trunk/Cosmo/RadioBeam/mdish.h

    r3756 r3769  
    6060  // Circular dish
    6161  Dish(int id, double x, double y, double diam)
    62     :  id_(id), X(x), Y(y), D(diam), Dx(0.), Dy(0.), fgcirc_(true)   {   }
     62    :  id_(id), X(x), Y(y), D(diam), Dx(D), Dy(D), fgcirc_(true)   {   }
    6363  // Receiver with rectangular type answer in kx,ky plane
    6464  Dish(int id, double x, double y, double dx, double dy)
    65     :  id_(id), X(x), Y(y), D(0.), Dx(dx), Dy(dy), fgcirc_(false)   {   }
     65    :  id_(id), X(x), Y(y), D(sqrt(dx*dy)), Dx(dx), Dy(dy), fgcirc_(false)   {   }
    6666  Dish(Dish const& a)
    6767    :  id_(a.id_), X(a.X), Y(a.Y), D(a.D), Dx(a.Dx), Dy(a.Dy), fgcirc_(a.fgcirc_)     {   }
    6868  inline bool isCircular() { return fgcirc_; }
    6969  inline int ReflectorId() { return id_; }
     70  inline double Diameter() { return D; }
     71  inline double DiameterX() { return Dx; }
     72  inline double DiameterY() { return Dx; }
    7073
    7174  int id_;   // numero de reflecteur
    … …  
    8386  void Define(r_8 xMin,r_8 xMax,int_4 nxBin,r_8 yMin,r_8 yMax,int_4 nyBin);
    8487  double Add(r_8 x, r_8 y, r_8 w, bool fgfh);
     88  void SetZeroBin(r_8 x=0, r_8 y=0);
    8589  inline double WBinX() { return dxb; }
    8690  inline double WBinY() { return dyb; }
    … …  
    9094  r_8 dxb,dyb;
    9195  sa_size_t nx,ny;
     96  sa_size_t ixb0, jyb0;
    9297  TArray<r_8> aw;
    9398  double sumw0;
    … …  
    109114  double CumulResp(Four2DResponse& rd, double theta=0., double phi=0.);
    110115  inline size_t NbDishes() { return dishes_.size(); }
     116  inline Dish&  operator[](size_t k) { return dishes_[k]; }
    111117
     118  virtual Histo2D PosDist(int nx=30, int ny=30, double dmax=0.); 
     119
     120protected:
    112121  double AddToHisto(double kx0, double ky0, double x, double y, double w, bool fgfh);
    113122

  • trunk/Cosmo/RadioBeam/pknoise.cc

    r3756 r3769  
    5555};
    5656
     57//-----------------------------------------------------------------------------------
     58//  Fonctions de creation de configuration d'interfero avec des dishs
     59//-----------------------------------------------------------------------------------
     60
     61vector<Dish> CreateFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
     62vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
     63vector<Dish> CreateConfigA(double Ddish=5., double Eta=0.9);
     64vector<Dish> CreateConfigB(double Ddish=5., double Eta=0.9);
     65vector<Dish> CreateConfigC(double Ddish=5., double Eta=0.9);
     66vector<Dish> CreateConfigD(double Ddish=5., double Eta=0.9);
     67
     68
     69vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW=10., double cylRL=0.5,
     70                                   double etaW=0.9, double etaRL=0.9, bool fgscid=true);
     71
     72
     73//-------------------------------------------------------------------------
     74//      ------------------ MAIN PROGRAM ------------------------------
     75//-------------------------------------------------------------------------
    5776int main(int narg, const char* arg[])
    5877{
    … …  
    115134    double Eta=0.95;
    116135    int cnt=0;
    117     vector<Dish> vdplein;
    118     for(int i=0; i<20; i++)
    119       for(int j=0; j<20; j++) {
    120         cnt++;
    121         vdplein.push_back(Dish(i*20+j+1, i*Ddish, j*Ddish, Eta*Ddish));
    122       }
    123     vector<Dish> vdsparse;
    124     vector<Dish> vdsparseD7;
    125 
    126     cnt=0;
    127     for(int i=0; i<=18; i++) {
    128       cnt++; vdsparse.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
    129       vdsparseD7.push_back(Dish(cnt, i*Ddish2,0.,Eta*Ddish2));
    130     }
    131     for(int i=-18; i<=18; i++) {
    132       if (i==0)  continue;
    133       cnt++; vdsparse.push_back(Dish(cnt, 0.,i*Ddish,Eta*Ddish));
    134       vdsparseD7.push_back(Dish(cnt, 0.,i*Ddish2,Eta*Ddish2));
    135     }
    136     for(int i=0; i<4; i++) {
    137         cnt++;  vdsparse.push_back(Dish(cnt, (3+2*i)*Ddish,(3+2*i)*Ddish,Eta*Ddish));
    138         vdsparseD7.push_back(Dish(cnt, (3+2*i)*Ddish2,(3+2*i)*Ddish2,Eta*Ddish2));
    139         cnt++;  vdsparse.push_back(Dish(cnt, (3+2*i)*Ddish,-(3+2*i)*Ddish,Eta*Ddish));
    140         vdsparseD7.push_back(Dish(cnt, (3+2*i)*Ddish2,-(3+2*i)*Ddish2,Eta*Ddish2));
    141         /*
    142         if ((i>0)||(j>0)) {
    143           cnt++;  vdsparse.push_back(Dish(cnt, (5+3*i)*Ddish,(3+2*j)*Ddish,Eta*Ddish));
    144           cnt++;  vdsparse.push_back(Dish(cnt, (5+3*i)*Ddish,-(3+2*j)*Ddish,Eta*Ddish));
    145         }
    146         */
    147       }
    148 
    149 
    150     vector<Dish> vcylplein, vcylplP;
    151     cnt=0;
     136    vector<Dish> vdplein = CreateFilledSqConfig(20, Ddish, Eta);
     137    vector<Dish> vdpl64 = CreateFilledSqConfig(8, Ddish, Eta);
     138
     139    vector<Dish> vdsparse = CreateConfigA(Ddish, Eta);
     140    vector<Dish> vdsparseD7 = CreateConfigA(Ddish2, Eta);
     141    //    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
     142    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
     143    vector<Dish> vdsparseC = CreateConfigC(Ddish, Eta);
     144
     145
    152146    double cylW=12.;   // Largeur des cylindres
    153147    double cylRL=0.5;  // Longeur des elements de reception le long du cylindre
    154     for(int i=0; i<8; i++)
    155       for(int j=0; j<192; j++) {
    156         vcylplein.push_back(Dish(i+1, i*cylW, j*cylRL, 0.9*cylW, 0.9*cylRL));
    157         cnt++; vcylplP.push_back(Dish(cnt, i*cylW, j*cylRL, 1.*cylW, 1.*cylRL));
    158       }
    159     vector<Dish> v2cyl, v2cylP;
    160     cnt=0;
    161     for(int i=0; i<2; i++)
    162       for(int j=0; j<32; j++) {
    163         v2cyl.push_back(Dish(i+1, i*25, j*cylRL, 0.9*9., 0.9*cylRL));
    164         cnt++; v2cylP.push_back(Dish(cnt, i*25, j*cylRL, 0.9*9., 1.*cylRL));
    165       }
     148    double etaW=0.95;   // Efficacite de couverture en largeur
     149    double etaRL=0.9;   // Efficacite de couverture le long du cylindre
     150    vector<Dish> vcylplein = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
     151    vector<Dish> vcylplP = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);
     152
     153    cylW=10.;
     154    cylRL=0.5;
     155    vector<Dish> v3cyl = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
     156    vector<Dish> v3cylP = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
     157    cylW=25.;
     158    cylRL=0.5;
     159    etaW=0.3;
     160    etaRL=0.9;
     161    vector<Dish> v2cyl = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, true);
     162    vector<Dish> v2cylP = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);
     163
    166164    double LMAX = D;
    167165    bool fgnoauto = true;
    168     int NRX=160;
    169     int NRY=160;
     166    int NRX=100;
     167    int NRY=100;
    170168
    171169    MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
    … …  
    174172    PrtTim("Apres mdfill.GetResponse()");
    175173
     174    MultiDish mdfill64(lambda, LMAX, vdpl64, fgnoauto);
     175    mdfill64.SetRespHisNBins(NRX,NRY);   
     176    {
     177      Histo2D hpos=mdfill64.PosDist(10,10,10.*Ddish);
     178      po << PPFNameTag("posf64") << hpos;
     179    }
     180    Histo2D hrf64 = mdfill64.GetResponse();
     181    PrtTim("Apres mdfill64.GetResponse()");
     182
    176183    MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto);
    177     mdsparse.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
     184    mdsparse.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    178185    mdsparse.SetRespHisNBins(NRX,NRY);
     186    {
     187      Histo2D hpos=mdsparse.PosDist(22,22,22.*Ddish);;
     188      po << PPFNameTag("posspA") << hpos;
     189    }
    179190    Histo2D hrsp = mdsparse.GetResponse();
    180191    PrtTim("Apres mdsparse.GetResponse()");
    181192
    182     MultiDish mdsparsefp(lambda, LMAX, vdsparse, fgnoauto);
    183     mdsparsefp.SetRespHisNBins(NRX,NRY);
    184     Histo2D hrspfp = mdsparsefp.GetResponse();
    185     PrtTim("Apres mdsparsefp.GetResponse()");
    186 
     193    /*
    187194    MultiDish mdsparseD7(lambda, LMAX, vdsparseD7, fgnoauto);
    188195    mdsparseD7.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
    … …  
    190197    Histo2D hrspd7 = mdsparseD7.GetResponse();
    191198    PrtTim("Apres mdsparseD7.GetResponse()");
    192                              
     199    */
     200    MultiDish mdsparseB(lambda, LMAX, vdsparseB, fgnoauto);
     201    mdsparseB.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
     202    mdsparseB.SetRespHisNBins(NRX,NRY);
     203    {
     204      Histo2D hpos=mdsparseB.PosDist(15,15,15.*Ddish);
     205      po << PPFNameTag("posspB") << hpos;
     206    }
     207    Histo2D hrspB = mdsparseB.GetResponse();
     208    PrtTim("Apres mdsparseB.GetResponse()");                       
     209
     210    MultiDish mdsparseC(lambda, LMAX, vdsparseC, fgnoauto);
     211    mdsparseC.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
     212    mdsparseC.SetRespHisNBins(NRX,NRY);
     213    {
     214      Histo2D hpos=mdsparseC.PosDist(20,20,20.*Ddish);
     215      po << PPFNameTag("posspC") << hpos;
     216    }
     217    Histo2D hrspC = mdsparseC.GetResponse();
     218    PrtTim("Apres mdsparseC.GetResponse()");                       
     219
     220    MultiDish mdsparseBfp(lambda, LMAX, vdsparseB, fgnoauto);
     221    mdsparseBfp.SetRespHisNBins(NRX,NRY);
     222    Histo2D hrspBfp = mdsparseBfp.GetResponse();
     223    PrtTim("Apres mdsparseBfp.GetResponse()");                     
     224
     225
    193226    MultiDish mcylfill(lambda, LMAX, vcylplein, fgnoauto);
    194227    mcylfill.SetRespHisNBins(NRX,NRY);
    … …  
    198231    mcylfillP.SetRespHisNBins(NRX,NRY);
    199232    Histo2D hfcylP = mcylfillP.GetResponse();
     233    PrtTim("Apres mcylfillP.GetResponse()");
     234
     235    MultiDish md3cyl(lambda, LMAX, v3cyl, fgnoauto);
     236    md3cyl.SetRespHisNBins(NRX,NRY);
     237    Histo2D h3cyl = md3cyl.GetResponse();
     238    PrtTim("Apres md3cyl.GetResponse()");
     239    MultiDish md3cylP(lambda, LMAX, v3cylP, fgnoauto);
     240    md3cylP.SetRespHisNBins(NRX,NRY);
     241    Histo2D h3cylP = md3cylP.GetResponse();
     242    PrtTim("Apres md3cylP.GetResponse()");
    200243
    201244    MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
    202245    md2cyl.SetRespHisNBins(NRX,NRY);
    203246    Histo2D h2cyl = md2cyl.GetResponse();
     247    PrtTim("Apres md2cyl.GetResponse()");
    204248    MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
    205249    md2cylP.SetRespHisNBins(NRX,NRY);
    206250    Histo2D h2cylP = md2cylP.GetResponse();
     251    PrtTim("Apres md2cylP.GetResponse()");
    207252
    208253    po << PPFNameTag("mfill") << hrfill;
     254    po << PPFNameTag("mfill64") << hrf64;
    209255    po << PPFNameTag("mspars") << hrsp;
    210     po << PPFNameTag("msparsfp") << hrspfp;
    211     po << PPFNameTag("msparsd7") << hrspd7;
     256    //    po << PPFNameTag("msparsd7") << hrspd7;
     257    po << PPFNameTag("msparsB") << hrspB;
     258    po << PPFNameTag("msparsC") << hrspC;
     259    po << PPFNameTag("msparsBfp") << hrspBfp;
    212260    po << PPFNameTag("mcylf") << hfcyl;
     261    po << PPFNameTag("m3cyl") << h3cyl;
    213262    po << PPFNameTag("m2cyl") << h2cyl;
    214263    po << PPFNameTag("mcylfP") << hfcylP;
     264    po << PPFNameTag("m3cylP") << h3cylP;
    215265    po << PPFNameTag("m2cylP") << h2cylP;
    216266
    … …  
    219269
    220270    Four2DRespTable mdf(hrfill, Dol);
     271    Four2DRespTable mdf64(hrf64, Dol);
    221272    Four2DRespTable mds(hrsp,   Dol);
    222     Four2DRespTable mdsfp(hrspfp,   Dol);
    223     Four2DRespTable mdsd7(hrspd7,   Dol);
     273    // Four2DRespTable mdsfp(hrspfp,   Dol);
     274    //    Four2DRespTable mdsd7(hrspd7,   Dol);
     275    Four2DRespTable mdsB(hrspB,   Dol);
     276    Four2DRespTable mdsC(hrspC,   Dol);
     277    Four2DRespTable mdsBfp(hrspBfp,   Dol);
    224278
    225279    Four2DRespTable mcylf(hfcyl, Dol);
     280    Four2DRespTable m3cyl(h3cyl, Dol);
    226281    Four2DRespTable m2cyl(h2cyl, Dol);
    227282    Four2DRespTable mcylfP(hfcylP, Dol);
     283    Four2DRespTable m3cylP(h3cylP, Dol);
    228284    Four2DRespTable m2cylP(h2cylP, Dol);
    229285
    … …  
    244300
    245301    cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
    246 #define NCONFIG 10
    247     Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mds, &mdsfp, &mdsd7, &mcylf, &mcylfP, &m2cyl, &m2cylP};
    248     const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","S63Dish5m","S63Dish5mFP","S63Dish7m",
    249                                "F8Cyl","F8CylP","BiCyl","BiCylP"};
    250     const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemds","noisemdsfp","noisemdsd7",
    251                                "noisefcyl","noisefcylP","noise2cyl","noise2cylP"};
     302#define NCONFIG 14
     303    Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mdf64, &mds, &mdsB, &mdsC, &mdsBfp,
     304                                    &mcylf, &mcylfP, &m3cyl, &m3cylP, &m2cyl, &m2cylP};
     305    const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","F8x8Dish5m",
     306                               "S68Dish5m","S72Dish5m","S129CDish5m","S72BDish5mFP",
     307                               "F8Cyl","F8CylP","F3Cyl","F3CylP","BiCyl","BiCylP"};
     308    const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemdf64","noisemds",
     309                               "noisemdsB","noisemdsC","noisemdsBfp",
     310                               "noisefcyl","noisefcylP","noise3cyl","noise3cylP", "noise2cyl","noise2cylP"};
    252311    vector<int> nbdishes;
    253312    nbdishes.push_back(1);
    254313    nbdishes.push_back(1);
    255314    nbdishes.push_back(vdplein.size());
     315    nbdishes.push_back(vdpl64.size());
    256316    nbdishes.push_back(vdsparse.size());
    257     nbdishes.push_back(vdsparse.size());
    258     nbdishes.push_back(vdsparseD7.size());
     317    // nbdishes.push_back(vdsparse.size());
     318    nbdishes.push_back(vdsparseB.size());
     319    nbdishes.push_back(vdsparseC.size());
     320    nbdishes.push_back(vdsparseB.size());
    259321    nbdishes.push_back(vcylplein.size());
    260322    nbdishes.push_back(vcylplP.size());
     323    nbdishes.push_back(v3cyl.size());
     324    nbdishes.push_back(v3cylP.size());
    261325    nbdishes.push_back(v2cyl.size());
    262326    nbdishes.push_back(v2cylP.size());
    … …  
    288352
    289353
     354//-----------------------------------------------------------------------------------
     355//-----------------------------------------------------------------------------------
     356//  Fonctions de creation de configuration d'interfero avec des dishs
     357//-----------------------------------------------------------------------------------
     358/* --Fonction -- */
     359vector<Dish> CreateFilledSqConfig(int nd, double Ddish, double Eta)
     360{
     361  vector<Dish>  vd;
     362  int cnt=0;
     363  for(int i=0; i<nd; i++)
     364    for(int j=0; j<nd; j++) {
     365      cnt++;
     366      vd.push_back(Dish(cnt, i*Ddish, j*Ddish, Eta*Ddish));
     367    }
     368  cout << ">>>CreateFilledSqConfig(" << nd << "," << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     369 
     370  return vd;
     371}
     372
     373/* --Fonction -- */
     374vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish, double Eta)
     375{
     376  vector<Dish>  vd;
     377  int cnt=0;
     378  int fgst=1;
     379  for(int i=0; i<nd; i++) {
     380    fgst = (fgst+1)%2;
     381    for(int j=0; j<nd; j++) {
     382      if (j%2==fgst)  continue;
     383      cnt++;
     384      vd.push_back(Dish(cnt, i*Ddish, j*Ddish, Eta*Ddish));
     385    }
     386  }
     387  cout << ">>>CreateSemiFilledSqConfig(" << nd << "," << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     388 
     389  return vd;
     390}
     391
     392/* --Fonction -- */
     393vector<Dish> CreateConfigA(double Ddish, double Eta)
     394{
     395  vector<Dish>  vd;
     396  int cnt=0;
     397  for(int i=0; i<18; i++) {
     398    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
     399    cnt++; vd.push_back(Dish(cnt, i*Ddish, 17.*Ddish,Eta*Ddish));
     400    if ((i>0)&&(i<17)) {
     401      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
     402      cnt++; vd.push_back(Dish(cnt,17.*Ddish,i*Ddish,Eta*Ddish));
     403    }
     404  }
     405  cout << ">>>CreateConfigA(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     406  return vd;
     407}
     408
     409/* --Fonction -- */
     410vector<Dish> CreateConfigB(double Ddish, double Eta)
     411{
     412  vector<Dish>  vd;
     413  int cnt=0;
     414  /*
     415  for(int i=0; i<13; i++) {
     416    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
     417    cnt++; vd.push_back(Dish(cnt, i*Ddish, 12.*Ddish,Eta*Ddish));
     418    if ((i>0)&&(i<12)) {
     419      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
     420      cnt++; vd.push_back(Dish(cnt,12.*Ddish,i*Ddish,Eta*Ddish));
     421    }
     422  }
     423  for(int i=0; i<5; i++) {
     424    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish,4.*Ddish,Eta*Ddish));
     425    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 8.*Ddish,Eta*Ddish));
     426    if ((i>0)&&(i<4)) {
     427      cnt++; vd.push_back(Dish(cnt,4.*Ddish,(i+4)*Ddish,Eta*Ddish));
     428      cnt++; vd.push_back(Dish(cnt,8.*Ddish,(i+4)*Ddish,Eta*Ddish));
     429    }
     430  }
     431  */
     432  for(int i=0; i<11; i++) {
     433    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
     434    cnt++; vd.push_back(Dish(cnt, i*Ddish, 10.*Ddish,Eta*Ddish));
     435    if ((i>0)&&(i<10)) {
     436      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
     437      cnt++; vd.push_back(Dish(cnt,10.*Ddish,i*Ddish,Eta*Ddish));
     438    }
     439  }
     440  for(int i=0; i<7; i++) {
     441    cnt++; vd.push_back(Dish(cnt, (i+2)*Ddish, 2.*Ddish,Eta*Ddish));
     442    cnt++; vd.push_back(Dish(cnt, (i+2)*Ddish, 8.*Ddish,Eta*Ddish));
     443    if ((i>0)&&(i<6)) {
     444      cnt++; vd.push_back(Dish(cnt,2.*Ddish,(i+2)*Ddish,Eta*Ddish));
     445      cnt++; vd.push_back(Dish(cnt,8.*Ddish,(i+2)*Ddish,Eta*Ddish));
     446    }
     447  }
     448  for(int i=0; i<3; i++) {
     449    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 4.*Ddish,Eta*Ddish));
     450    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 6.*Ddish,Eta*Ddish));
     451    if ((i>0)&&(i<2)) {
     452      cnt++; vd.push_back(Dish(cnt,4.*Ddish,(i+4)*Ddish,Eta*Ddish));
     453      cnt++; vd.push_back(Dish(cnt,6.*Ddish,(i+4)*Ddish,Eta*Ddish));
     454    }
     455  }
     456
     457
     458  cout << ">>>CreateConfigB(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     459  return vd;
     460}
     461
     462
     463/* --Fonction -- */
     464vector<Dish> CreateConfigC(double Ddish, double Eta)
     465{
     466  vector<int_4> lesx, lesy;
     467 
     468  int max = 16;
     469  for(int i=0; i<4; i++)
     470    for(int j=0; j<4; j++) {
     471      lesx.push_back(i); lesy.push_back(j);
     472      lesx.push_back(max-i); lesy.push_back(max-j);
     473      lesx.push_back(i); lesy.push_back(max-j);
     474      lesx.push_back(max-i); lesy.push_back(j);
     475    }
     476 
     477  for(int i=5; i<12; i+=2) {
     478    lesx.push_back(i); lesy.push_back(0);
     479    lesx.push_back(i); lesy.push_back(max);
     480    lesx.push_back(0); lesy.push_back(i);
     481    lesx.push_back(max); lesy.push_back(i);
     482  }
     483 
     484  for(int i=4; i<=12; i+=2)
     485    for(int j=4; j<=12; j+=2) {
     486      lesx.push_back(i); lesy.push_back(j);
     487    }
     488 
     489  for(int i=5; i<=11; i+=2) {
     490    lesx.push_back(i); lesy.push_back(4);
     491    lesx.push_back(i); lesy.push_back(max-4);
     492    lesx.push_back(4); lesy.push_back(i);
     493    lesx.push_back(max-4); lesy.push_back(i);
     494  }
     495 
     496  EnumeratedSequence esx,esy;
     497  esx = 2,5;
     498  esy = 5,2;
     499 
     500  for(int k=0; k<esx.Size(); k++) {
     501    int_4 ix=esx.Value(k);
     502    int_4 iy=esy.Value(k);
     503   
     504    lesx.push_back(ix); lesy.push_back(iy);
     505    lesx.push_back(max-ix); lesy.push_back(iy);
     506    lesx.push_back(ix); lesy.push_back(max-iy);
     507    lesx.push_back(max-ix); lesy.push_back(max-iy);
     508  }
     509 cout << "CreateConfigC/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
     510 
     511 vector<Dish> vd;
     512 int cnt=0;
     513 for(size_t i=0; i<lesx.size(); i++) {
     514   cnt++; vd.push_back(Dish(cnt, ((double)lesx[i])*Ddish,((double)lesy[i])*Ddish,Eta*Ddish));
     515 }
     516
     517 cout << ">>>CreateConfigC(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     518
     519 return vd;
     520
     521}
     522/* --Fonction -- */
     523vector<Dish> CreateConfigD(double Ddish, double Eta)
     524{
     525EnumeratedSequence es;
     526es = 0,1,3,4,6,8,10,12,14,16,18,20,22,24,25,27,28;
     527vector<int_4> lesx, lesy;
     528for(int k=0; k<es.Size(); k++) {
     529  lesx.push_back(es.Value(k));   lesy.push_back(0);
     530  lesx.push_back(es.Value(k));   lesy.push_back(28);
     531}
     532for(int k=1; k<es.Size()-1; k++) {
     533  lesy.push_back(es.Value(k));   lesx.push_back(0);
     534  lesy.push_back(es.Value(k));   lesx.push_back(28);
     535}
     536for(int k=1; k<=5; k++) {
     537  lesy.push_back(k);   lesx.push_back(5);
     538  lesy.push_back(28-k);   lesx.push_back(28-5);
     539  if (k!=5) {
     540    lesx.push_back(k);   lesy.push_back(5);
     541    lesx.push_back(28-k);   lesy.push_back(28-5);
     542  }
     543
     544  lesy.push_back(k);   lesx.push_back(28-5);
     545  lesy.push_back(28-k);   lesx.push_back(5);
     546  if (k!=5) {
     547    lesx.push_back(28-k);   lesy.push_back(5);
     548    lesx.push_back(k);   lesy.push_back(28-5);
     549  }
     550}
     551
     552for(int k=6; k<=13; k++) {
     553  lesy.push_back(k);   lesx.push_back(k);
     554  lesy.push_back(28-k);   lesx.push_back(28-k);
     555  lesy.push_back(k);   lesx.push_back(28-k);
     556  lesy.push_back(28-k);   lesx.push_back(k);
     557}
     558
     559
     560lesx.push_back(14); lesy.push_back(14);
     561
     562EnumeratedSequence esx,esy;
     563esx = 1,2,4;
     564esy = 12,11,13;
     565
     566for(int k=0; k<esx.Size(); k++) {
     567  int_4 ix=esx.Value(k);
     568  int_4 iy=esy.Value(k);
     569  lesx.push_back(ix); lesy.push_back(iy);
     570  lesx.push_back(28-ix); lesy.push_back(iy);
     571  lesx.push_back(ix); lesy.push_back(28-iy);
     572  lesx.push_back(28-ix); lesy.push_back(28-iy);
     573 
     574  lesy.push_back(ix); lesx.push_back(iy);
     575  lesy.push_back(28-ix); lesx.push_back(iy);
     576  lesy.push_back(ix); lesx.push_back(28-iy);
     577  lesy.push_back(28-ix); lesx.push_back(28-iy);
     578 } 
     579for(int k=5; k<=13; k+=2) {
     580  lesy.push_back(k);   lesx.push_back(14);
     581  lesy.push_back(28-k);   lesx.push_back(14);
     582  lesy.push_back(14);   lesx.push_back(k);
     583  lesy.push_back(14);   lesx.push_back(28-k);
     584}
     585
     586 cout << "CreateConfigB/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
     587 
     588 vector<Dish> vd;
     589 int cnt=0;
     590 for(size_t i=0; i<lesx.size(); i++) {
     591   cnt++; vd.push_back(Dish(cnt, ((double)lesx[i])*Ddish,((double)lesy[i])*Ddish,Eta*Ddish));
     592 }
     593
     594 cout << ">>>CreateConfigD(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
     595
     596 return vd;
     597}
     598
     599/* --Fonction -- */
     600vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW, double cylRL, double etaW, double etaRL, bool fgscid)
     601{
     602  vector<Dish>  vd;
     603  int cnt=0;
     604
     605  for(int i=0; i<ncyl; i++)
     606    for(int j=0; j<nRL; j++) {
     607      cnt++;
     608      int rid = (fgscid) ? i+1 : cnt;
     609      vd.push_back(Dish(rid, i*cylW, j*cylRL, etaW*cylW, etaRL*cylRL));
     610    }
     611  cout << ">>>CreateFilledCylConfig(" << ncyl << "," << nRL << "," << cylW << "," << cylRL << ","
     612       << etaW << "," << etaRL << "," << ((fgscid)?" RId=CylNum":"Cnt")
     613       << ") ---> NDishes=" << vd.size() << endl;
     614
     615  return vd;
     616}

  • trunk/Cosmo/RadioBeam/plpkn.pic

    r3756 r3769  
    55########################################################################
    66
    7 openppf pknoise.ppf
     7if ( $# < 1 ) then
     8  echo ' Usage: exec plpkn PPFName_pknoise'
     9  return
     10endif
     11
     12echo "---> openppf $1 "
     13openppf $1
     14
     15echo '----> Executing anapkn.pic'
     16exec anapkn.pic
     17setup5
     18scalewz 0.7 2500 100 1
     19y1 = 500*$cct21
     20y2 = 9e4*$cct21
     21xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
     22defscript plpklss
     23  n/plot hpkz.val*${cct21}%$kk ! ! "same notit nsta connectpoints black $xyl line=solid,2"
     24  addtext 0.03 2000 '*** P(k)-LSS ***' 'font=helvetica,bolditalic,16 black'
     25endscript 
    826
    927setaxesatt 'minorticks font=helvetica,bold,16 autofontsize'
    10 disp noiseD 'logy nsta'
    11 disp noiseD2 'same grey nsta'
    12 disp noisemdf 'same orange nsta'
    13 disp noisemds 'same red nsta'
    14 disp noisemdsfp 'same yellow nsta'
    15 disp noisemdsd7 'same gold nsta'
    16 disp noisefcyl 'same blue nsta'
    17 disp noisefcylP 'same skyblue nsta'
    18 disp noise2cyl 'same forestgreen nsta'
    19 disp noise2cylP 'same green nsta'
    2028
    2129Rad2Deg = 180/3.141596
    … …  
    5563endscript
    5664
    57 Dx = 100
    58 Dy = 100
    59 calcul
    60 plot2d noiseD x/$Da val*$PNOISE nb>10 'line=solid,2 logy logx xylimits=0.002,0.8,1,1e5 grid cpts nsta notit'
    61 Dx = 200
    62 Dy = 200
    63 calcul
    64 plot2d noiseD2 x/$Da val*$PNOISE nb>10 'same line=solid,2  cpts grey nsta notit'
    65 Dx = 5*0.95
    66 Dy = 5*0.95
    67 calcul
    68 plot2d noisemdf x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts orange nsta notit'
    69 plot2d noisemds x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
    70 Dx = 7.55*0.95
    71 Dy = 7.5*0.95
    72 calcul
    73 plot2d noisemdsd7 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts siennared nsta notit'
    74 Dx = 12*0.9
    75 Dy = 0.5*0.9
    76 calcul
    77 plot2d noisefcyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts blue nsta notit'
    78 plot2d noisefcylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts skyblue nsta notit'
    79 Dx = 8*0.9
    80 Dy = 0.5*0.9
    81 calcul
    82 plot2d noise2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts forestgreen nsta notit'
    83 plot2d noise2cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
    84 
    85 set lines ( '100mDish' '200mDish' '400x5mDish' '63x5mDish' '63x7.5mDish'  '8Cyl12x96m' 'Perf8Cyl12x96m' )
    86 set lines ( $lines 'Pitts2Cyl' 'PerfPitts2Cyl' )
    87 set cols ( black grey  orange red siennared  blue skyblue forestgreen green )
    88 textdrawer lines cols 'frame font=helvetica,bold,16'
    89 setaxelabels 'k (Mpc^-1) ' 'PNoise(k) mk^2 Mpc^3' 'font=helvetica,bolditalic,16'
    90 
    91 
     65defscript plnoisedish
     66  Dx = 100
     67  Dy = 100
     68  calcul
     69  plot2d noiseD x/$Da val*$PNOISE nb>10 'line=solid,2 logy logx xylimits=0.002,0.8,1,1e5 navyblue grid cpts nsta notit'
     70  Dx = 200
     71  Dy = 200
     72  calcul
     73  plot2d noiseD2 x/$Da val*$PNOISE nb>10 'same line=solid,2  cpts grey nsta notit'
     74  set lines ( '100mDish' '200mDish'  )
     75  set cols ( navyblue grey  )
     76  textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.1,0.3,0.7,0.8'
     77  setaxelabels 'k (Mpc^-1) ' 'PNoise(k) mk^2 Mpc^3' 'font=helvetica,bolditalic,16'
     78endscript
     79
     80defscript plnoiseA
     81  Dx = 5*0.95
     82  Dy = 5*0.95
     83  calcul
     84  plot2d noisemdf64 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     85#  plot2d noisemds x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     86  plot2d noisemdsB x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts siennared nsta notit'
     87  plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violetred nsta notit'
     88  Dx = 25*0.3
     89  Dy = 0.5*0.9
     90  calcul
     91  plot2d noise2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts forestgreen nsta notit'
     92  plot2d noise2cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
     93  set lines (  'FilledA:64x5mD' 'SparseB:72x5mD' 'SparseC:129x5mD' 'Pitts2Cyl(=64C)' 'PerfPitts2Cyl(=64C)' )
     94  set cols ( red siennared magenta  forestgreen green )
     95  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35'
     96  settitle ' PNoise(k) : Dishes/Cylinders, 64/72/129 channels' ' ' 'font=helvetica,bold,16'
     97endscript
     98
     99defscript plnoiseB
     100  Dx = 5*0.95
     101  Dy = 5*0.95
     102  calcul
     103  plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     104  plot2d noisemdf x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts gold nsta notit'
     105  Dx = 10*0.95
     106  Dy = 0.5*0.9
     107  calcul
     108  plot2d noise3cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit'
     109  plot2d noise3cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts magenta nsta notit'
     110
     111  Dx = 12*0.95
     112  Dy = 0.5*0.9
     113  calcul
     114  plot2d noisefcyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts blue nsta notit'
     115  plot2d noisefcylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts royalblue nsta notit'
     116
     117  set lines ( 'SparseC:129x5mD' 'Filled400x5m' '3xCyl10x64m(=384C)' 'Perf3xCyl10x64m(=384C)' )
     118  set lines ( $lines '8Cyl12x96m(=1536C)' 'Perf8Cyl12x96m(=1536C)' )
     119  set cols ( red gold  violet magenta  blue royalblue )
     120  textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.25,0.55,0.6,0.95'
     121  settitle ' PNoise(k) : Dishes/Cylinders, 129/384/400/1536 channels' ' ' 'font=helvetica,bold,16'     
     122endscript
     123
     124defscript xxx
    92125newwin
    93126disp mfill 'h2disp=img colbr128 notit'
    … …  
    122155setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    123156settitle ' u-v coverage , Perfect-Pitts. 2 Cyl 16mx8m , dist=25m ' ' ' 'font=helvetica,bold,16'
     157
     158
     159disp noiseD 'logy nsta'
     160disp noiseD2 'same grey nsta'
     161disp noisemdf 'same cyan nsta'
     162disp noisemdf64 'same brown nsta'
     163disp noisemds 'same red nsta'
     164disp noisemdsC 'same orange nsta'
     165
     166disp noisemdsB 'same yellow nsta'
     167disp noisemdsfp 'same yellow nsta'
     168# disp noisemdsd7 'same gold nsta'
     169disp noisefcyl 'same blue nsta'
     170disp noisefcylP 'same skyblue nsta'
     171disp noise3cylP 'same magenta nsta'
     172disp noise3cyl 'same violet nsta'
     173disp noise2cyl 'same forestgreen nsta'
     174disp noise2cylP 'same green nsta'
     175
     176endscript
     177
     178defscript AA
     179  plnoisedish
     180  plpklss
     181  plnoiseA
     182endscript
     183
     184defscript BB
     185  plnoisedish
     186  plpklss
     187  plnoiseB
     188endscript
     189
     190defscript POSCOV
     191  disp posspB red
     192  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16'
     193  settitle ' Config B dish positions - 72 dishes  ' ' ' 'font=helvetica,bold,16'
     194#  w2ps
     195   
     196  disp posspC red
     197  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16'
     198  settitle ' Config C dish positions - 129 dishes  ' ' ' 'font=helvetica,bold,16'
     199#  w2ps
     200  disp mfill64 'h2disp=img colbr128 notit'
     201  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
     202  settitle ' u-v coverage , Filled 8x8 - 64 x 5m Dishes' ' ' 'font=helvetica,bold,16'
     203#  w2ps
     204  disp m3cyl 'h2disp=img colbr128 notit'
     205  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
     206  settitle ' u-v coverage , 3 Cylinders 10mx64m ' ' ' 'font=helvetica,bold,16'
     207#  w2ps
     208  disp mfill 'h2disp=img colbr128 notit'
     209  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
     210  settitle ' u-v coverage , Filled 20x20 - 400 x 5m Dishes' ' ' 'font=helvetica,bold,16'
     211#  w2ps
     212  disp msparsC 'h2disp=img colbr128 notit'
     213  setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
     214  settitle 'u-v coverage, Sparse-C: 129x5mD Over 80mx80m (Rot~Pi/4)' ' ' 'font=helvetica,bold,16'
     215 
     216endscript 

  • trunk/Cosmo/RadioBeam/specpk.cc

    r3756 r3769  
    154154  // sa_size_t is large integer type 
    155155  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
    156     kzz =  (kz>fourAmp.SizeZ()/2) ? (double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
     156    kzz =  (kz>fourAmp.SizeZ()/2) ? -(double)(fourAmp.SizeZ()-kz)*dkz_ : (double)kz*dkz_;
    157157    for(sa_size_t ky=0; ky<fourAmp.SizeY(); ky++) {
    158       kyy =  (ky>fourAmp.SizeY()/2) ? (double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
     158      kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
    159159      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
    160160        double kxx=(double)kx*dkx_;
    … …  
    194194// cells with amp^2=re^2+im^2>s2cut are ignored
    195195// Output : power spectrum (profile histogram)
    196 HProf Four3DPk::ComputePk(double s2cut, double kmin, double kmax, int nbin)
     196HProf Four3DPk::ComputePk(double s2cut, int nbin, double kmin, double kmax)
    197197{
    198198  // The second half of the array along Y (matrix rows) contain
    … …  
    202202  // as a function of wave-number k = sqrt((double)(kx*kx+ky*ky))
    203203  //  if (nbin < 1) nbin = nbh/2;
     204  if ((kmax<0.)||(kmax<kmin)) {
     205    kmin=0.;
     206    double maxx=fourAmp.SizeX()*dkx_;
     207    double maxy=fourAmp.SizeY()*dky_/2;
     208    double maxz=fourAmp.SizeZ()*dkz_/2;
     209    kmax=sqrt(maxx*maxx+maxy*maxy+maxz*maxz);
     210  }
     211  if (nbin<2) nbin=128;
    204212  HProf hp(kmin, kmax, nbin);
    205213  hp.SetErrOpt(false);

  • trunk/Cosmo/RadioBeam/specpk.h

    r3756 r3769  
    5959
    6060// Return the reconstructed power spectrum as a profile histogram   
    61 //  HProf ComputePk(double s2cut=0., double kmin=-0.5*DeuxPI, double kmax=127.5*DeuxPI, int nbin=256);
    62   HProf ComputePk(double s2cut=0., double kmin=-DeuxPI, double kmax=255*DeuxPI, int nbin=256);
     61  HProf ComputePk(double s2cut=0., int nbin=256, double kmin=0., double kmax=-1.);
    6362  void  ComputePkCumul(HProf& hp, double s2cut=0.);
    6463
Note: See TracChangeset for help on using the changeset viewer.