source: Sophya/trunk/Cosmo/RadioBeam/pknoise.cc@ 3790

#include "machdefs.h"
#include "sopnamsp.h"
#include <iostream>
#include <string>
#include <math.h>

#include <typeinfo>

#include "mdish.h"
#include "specpk.h"
#include "histinit.h"
// #include "fiosinit.h"
// #include "fitsioserver.h"

#include "randr48.h"      

#include "timing.h"
#include "ctimer.h"

typedef DR48RandGen RandomGenerator ;

// ---------------------------------------------------------------------
// Test main program for computing interferometer noise power spectrum
// R. Ansari  - Avril 2010 
// ---------------------------------------------------------------------

class PkNoiseCalculator
{
public:
  PkNoiseCalculator(Four3DPk& pk3, Four2DResponse& rep, 
                    double s2cut=100., int ngen=1, const char* tit="PkNoise")
    : pkn3d(pk3), frep(rep), S2CUT(s2cut), NGEN(ngen), title(tit) 
  { } 
  inline void SetS2Cut(double s2cut=100.)
  {  S2CUT=s2cut; }
  HProf Compute()
  {
    Timer tm(title.c_str());
    tm.Nop();
    HProf hnd;
    cout << "PkNoiseCalculator::Compute() " << title << "  NGEN=" << NGEN << " S2CUT=" << S2CUT << endl;
    for(int igen=0; igen<NGEN; igen++) {
      pkn3d.ComputeNoiseFourierAmp(frep);
      if (igen==0) hnd = pkn3d.ComputePk(S2CUT);
      else pkn3d.ComputePkCumul(hnd,S2CUT);
    }
    return hnd;
  }

  Four3DPk& pkn3d;
  Four2DResponse& frep;
  double S2CUT;
  int NGEN;
  string title;
};

//-----------------------------------------------------------------------------------
//  Fonctions de creation de configuration d'interfero avec des dishs
//-----------------------------------------------------------------------------------

vector<Dish> CreateFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish=5., double Eta=0.9);
vector<Dish> CreateConfigA(double Ddish=5., double Eta=0.9);
vector<Dish> CreateConfigB(double Ddish=5., double Eta=0.9);
vector<Dish> CreateConfigC(double Ddish=5., double Eta=0.9);
vector<Dish> CreateConfigD(double Ddish=5., double Eta=0.9);


vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW=10., double cylRL=0.5, 
                                   double etaW=0.9, double etaRL=0.9, bool fgscid=true);


//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
int main(int narg, const char* arg[])
{
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  {
    cout<< " Usage: pknoise [OutPPFName NGen S2Cut Lambda] " << endl;
    cout<< "   Default: OutPPFName=pknoise.ppf, NGen=1 " << endl;
    cout<< "            S2CUT=0. , Lambda=0.357  \n" << endl;

    return 1;
  }
  cout << " ==== pknoise.cc program , test of SpectralShape and MassDist2D classes ==== " << endl;
  // make sure SOPHYA modules are initialized 
  SophyaInit();  
  //  FitsIOServerInit();
  InitTim();
  //--- decoding command line arguments 
  string outfile = "pknoise.ppf";  
  if (narg>1) outfile = arg[1];
  if (outfile==".")  outfile = "pknoise.ppf";
  int NMAX = 1;
  if (narg>2) NMAX = atoi(arg[2]);
  double SCut=0.;
  if (narg>3) SCut = atof(arg[3]);
  double LAMBDA=0.357 ;  // 21 cm at z=0.7 
  if (narg>4) LAMBDA = atof(arg[4]);

  //-- end command line arguments
  
  int rc = 1;  
  try {  // exception handling try bloc at top level
    cout << "0/ pknoise.cc: Executing, output file= " << outfile << endl;  
    POutPersist po(outfile);
    cout << " 1.a/ Instanciating object type SpectralShape " << endl;
    SpectralShape spec(2);
    cout << " 1.b/ Wrinting spectral shape vector (name= Pk) to output PPF " << endl; 
    Histo hpk = spec.GetPk(1024);
    po << PPFNameTag("Pk") << hpk;

    double D = 100.;
    double lambda = LAMBDA;
    double Dol = D/lambda;
    cout << " 2.a/ Instanciating Four2DResponse(1/2/3...) " << endl;
    Four2DResponse dishg(1,Dol,Dol);
    Four2DResponse dish(2,Dol,Dol);
    Four2DResponse dish2(2,Dol*2.,Dol*2.);
    Four2DResponse dishsq(3,Dol,Dol/5.);
    cout << " 2.b/ Writing Four2DResponse Histo2D to output ppf " << endl;
    Histo2D hdg = dishg.GetResponse();
    Histo2D hd = dish.GetResponse();
    Histo2D hd2 = dish2.GetResponse();
    Histo2D hdsq = dishsq.GetResponse();
    po << PPFNameTag("dishg") << hdg;
    po << PPFNameTag("dish") << hd;
    po << PPFNameTag("dish2") << hd2;
    po << PPFNameTag("dishsq") << hdsq;

    cout << " 2.c/ Creating MultiDish Filled Array " << endl;
    double Ddish = 5.;
    double Ddish2 = 7.5;
    double Eta=0.95;
    int cnt=0;
    vector<Dish> vdplein = CreateFilledSqConfig(20, Ddish, Eta);
    vector<Dish> vdpl64 = CreateFilledSqConfig(8, Ddish, Eta);

    vector<Dish> vdsparse = CreateConfigA(Ddish, Eta);
    vector<Dish> vdsparseD7 = CreateConfigA(Ddish2, Eta);
    //    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
    vector<Dish> vdsparseB = CreateConfigB(Ddish, Eta);
    vector<Dish> vdsparseC = CreateConfigC(Ddish, Eta);


    double cylW=12.;   // Largeur des cylindres 
    double cylRL=0.5;  // Longeur des elements de reception le long du cylindre
    double etaW=0.95;   // Efficacite de couverture en largeur
    double etaRL=0.9;   // Efficacite de couverture le long du cylindre
    vector<Dish> vcylplein = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
    vector<Dish> vcylplP = CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);

    cylW=10.;
    cylRL=0.5; 
    vector<Dish> v3cyl = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
    vector<Dish> v3cylP = CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
    cylW=25.;
    cylRL=0.5; 
    etaW=0.3;
    etaRL=0.9;
    vector<Dish> v2cyl = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, true);
    vector<Dish> v2cylP = CreateFilledCylConfig(2, 32, cylW, cylRL, etaW, etaRL, false);

    double LMAX = D;
    bool fgnoauto = true;
    int NRX=100;
    int NRY=100;

    MultiDish mdfill(lambda, LMAX, vdplein, fgnoauto);
    mdfill.SetRespHisNBins(NRX,NRY);   
    Histo2D hrfill = mdfill.GetResponse();
    PrtTim("Apres mdfill.GetResponse()");

    MultiDish mdfill64(lambda, LMAX, vdpl64, fgnoauto);
    mdfill64.SetRespHisNBins(NRX,NRY);   
    {
      Histo2D hpos=mdfill64.PosDist(10,10,10.*Ddish);
      po << PPFNameTag("posf64") << hpos;
    }
    Histo2D hrf64 = mdfill64.GetResponse();
    PrtTim("Apres mdfill64.GetResponse()");

    MultiDish mdsparse(lambda, LMAX, vdsparse, fgnoauto); 
    mdsparse.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    mdsparse.SetRespHisNBins(NRX,NRY);
    {
      Histo2D hpos=mdsparse.PosDist(22,22,22.*Ddish);;
      po << PPFNameTag("posspA") << hpos;
    }
    Histo2D hrsp = mdsparse.GetResponse();
    PrtTim("Apres mdsparse.GetResponse()");

    /*
    MultiDish mdsparseD7(lambda, LMAX, vdsparseD7, fgnoauto); 
    mdsparseD7.SetThetaPhiRange(M_PI/4.,16, M_PI/4., 16);
    mdsparseD7.SetRespHisNBins(NRX,NRY);
    Histo2D hrspd7 = mdsparseD7.GetResponse();
    PrtTim("Apres mdsparseD7.GetResponse()");
    */ 
    MultiDish mdsparseB(lambda, LMAX, vdsparseB, fgnoauto); 
    mdsparseB.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    mdsparseB.SetRespHisNBins(NRX,NRY);
    {
      Histo2D hpos=mdsparseB.PosDist(15,15,15.*Ddish);
      po << PPFNameTag("posspB") << hpos;
    }
    Histo2D hrspB = mdsparseB.GetResponse();
    PrtTim("Apres mdsparseB.GetResponse()");                        

    MultiDish mdsparseC(lambda, LMAX, vdsparseC, fgnoauto); 
    mdsparseC.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    mdsparseC.SetRespHisNBins(NRX,NRY);
    {
      Histo2D hpos=mdsparseC.PosDist(20,20,20.*Ddish);
      po << PPFNameTag("posspC") << hpos;
    }
    Histo2D hrspC = mdsparseC.GetResponse();
    PrtTim("Apres mdsparseC.GetResponse()");                        

    MultiDish mdsparseBfp(lambda, LMAX, vdsparseB, fgnoauto); 
    mdsparseBfp.SetRespHisNBins(NRX,NRY);
    Histo2D hrspBfp = mdsparseBfp.GetResponse();
    PrtTim("Apres mdsparseBfp.GetResponse()");                      


    MultiDish mcylfill(lambda, LMAX, vcylplein, fgnoauto);
    mcylfill.SetRespHisNBins(NRX,NRY);
    Histo2D hfcyl = mcylfill.GetResponse();
    PrtTim("Apres mcylfill.GetResponse()");
    MultiDish mcylfillP(lambda, LMAX, vcylplP, fgnoauto);
    mcylfillP.SetRespHisNBins(NRX,NRY);
    Histo2D hfcylP = mcylfillP.GetResponse();
    PrtTim("Apres mcylfillP.GetResponse()");

    MultiDish md3cyl(lambda, LMAX, v3cyl, fgnoauto);
    md3cyl.SetRespHisNBins(NRX,NRY);
    Histo2D h3cyl = md3cyl.GetResponse();
    PrtTim("Apres md3cyl.GetResponse()");
    MultiDish md3cylP(lambda, LMAX, v3cylP, fgnoauto);
    md3cylP.SetRespHisNBins(NRX,NRY);
    Histo2D h3cylP = md3cylP.GetResponse();
    PrtTim("Apres md3cylP.GetResponse()");

    MultiDish md2cyl(lambda, LMAX, v2cyl, fgnoauto);
    md2cyl.SetRespHisNBins(NRX,NRY);
    Histo2D h2cyl = md2cyl.GetResponse();
    PrtTim("Apres md2cyl.GetResponse()");
    MultiDish md2cylP(lambda, LMAX, v2cylP, fgnoauto);
    md2cylP.SetRespHisNBins(NRX,NRY);
    Histo2D h2cylP = md2cylP.GetResponse();
    PrtTim("Apres md2cylP.GetResponse()");

    po << PPFNameTag("mfill") << hrfill;
    po << PPFNameTag("mfill64") << hrf64;
    po << PPFNameTag("mspars") << hrsp;
    //    po << PPFNameTag("msparsd7") << hrspd7;
    po << PPFNameTag("msparsB") << hrspB;
    po << PPFNameTag("msparsC") << hrspC;
    po << PPFNameTag("msparsBfp") << hrspBfp;
    po << PPFNameTag("mcylf") << hfcyl;
    po << PPFNameTag("m3cyl") << h3cyl;
    po << PPFNameTag("m2cyl") << h2cyl;
    po << PPFNameTag("mcylfP") << hfcylP;
    po << PPFNameTag("m3cylP") << h3cylP;
    po << PPFNameTag("m2cylP") << h2cylP;

    PrtTim("Done computing multi-dish response");
 

    Four2DRespTable mdf(hrfill, Dol);
    Four2DRespTable mdf64(hrf64, Dol);
    Four2DRespTable mds(hrsp,   Dol);
    // Four2DRespTable mdsfp(hrspfp,   Dol);
    //    Four2DRespTable mdsd7(hrspd7,   Dol);
    Four2DRespTable mdsB(hrspB,   Dol);
    Four2DRespTable mdsC(hrspC,   Dol);
    Four2DRespTable mdsBfp(hrspBfp,   Dol);

    Four2DRespTable mcylf(hfcyl, Dol);
    Four2DRespTable m3cyl(h3cyl, Dol);
    Four2DRespTable m2cyl(h2cyl, Dol);
    Four2DRespTable mcylfP(hfcylP, Dol);
    Four2DRespTable m3cylP(h3cylP, Dol);
    Four2DRespTable m2cylP(h2cylP, Dol);

    
    cout << " 3.a/ Instanciating object type Four3DPk  " << endl;
    RandomGenerator rg;
    Four3DPk m3d(rg);
    m3d.SetCellSize(2.*DeuxPI, 2.*DeuxPI, 2.*DeuxPI); 
    cout << " 3.b/ call ComputeFourierAmp()  NGEN=" << NMAX << endl;
    HProf hrpk;
    for(int igen=0; igen<NMAX; igen++) {
      m3d.ComputeFourierAmp(spec);
      if (igen==0)  hrpk = m3d.ComputePk();
      else m3d.ComputePkCumul(hrpk);
    }
    PrtTim("md.ComputeFourierAmp() done");
    po << PPFNameTag("recPk") << hrpk;

    cout << " 4/ Computing Noise P(k) using PkNoiseCalculator ..." << endl;
#define NCONFIG 14
    Four2DResponse* f2rep[NCONFIG]={&dish, &dish2, &mdf, &mdf64, &mds, &mdsB, &mdsC, &mdsBfp,
                                    &mcylf, &mcylfP, &m3cyl, &m3cylP, &m2cyl, &m2cylP};
    const char* tits[NCONFIG]={"Dish100m", "Dish200m","F20x20Dish5m","F8x8Dish5m",
                               "S68Dish5m","S72Dish5m","S129CDish5m","S72BDish5mFP",
                               "F8Cyl","F8CylP","F3Cyl","F3CylP","BiCyl","BiCylP"};
    const char* tags[NCONFIG]={"noiseD", "noiseD2","noisemdf","noisemdf64","noisemds",
                               "noisemdsB","noisemdsC","noisemdsBfp",
                               "noisefcyl","noisefcylP","noise3cyl","noise3cylP", "noise2cyl","noise2cylP"};
    vector<int> nbdishes;
    nbdishes.push_back(1);
    nbdishes.push_back(1);
    nbdishes.push_back(vdplein.size());
    nbdishes.push_back(vdpl64.size());
    nbdishes.push_back(vdsparse.size());
    // nbdishes.push_back(vdsparse.size());
    nbdishes.push_back(vdsparseB.size());
    nbdishes.push_back(vdsparseC.size());
    nbdishes.push_back(vdsparseB.size());
    nbdishes.push_back(vcylplein.size());
    nbdishes.push_back(vcylplP.size());
    nbdishes.push_back(v3cyl.size());
    nbdishes.push_back(v3cylP.size());
    nbdishes.push_back(v2cyl.size());
    nbdishes.push_back(v2cylP.size());

    for(int lc=0; lc<NCONFIG; lc++) {
      PkNoiseCalculator pkn(m3d, *(f2rep[lc]), SCut/(double)nbdishes[lc], NMAX, tits[lc]);
      HProf hpn = pkn.Compute();
      po << PPFNameTag(tags[lc]) << hpn;
    } 
    rc = 0;
  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " pknoise.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << "\n...exc.Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {  // catching standard C++ exceptions
    cerr << " pknoise.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " pknoise.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  PrtTim("End-pknoise");
  cout << " ==== End of pknoise.cc program  Rc= " << rc << endl;
  return rc;    
}


//-----------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------
//  Fonctions de creation de configuration d'interfero avec des dishs
//-----------------------------------------------------------------------------------
/* --Fonction -- */
vector<Dish> CreateFilledSqConfig(int nd, double Ddish, double Eta)
{
  vector<Dish>  vd;
  int cnt=0;
  for(int i=0; i<nd; i++) 
    for(int j=0; j<nd; j++) {
      cnt++; 
      vd.push_back(Dish(cnt, i*Ddish, j*Ddish, Eta*Ddish));
    }
  cout << ">>>CreateFilledSqConfig(" << nd << "," << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
  
  return vd;
}

/* --Fonction -- */
vector<Dish> CreateSemiFilledSqConfig(int nd, double Ddish, double Eta)
{
  vector<Dish>  vd;
  int cnt=0;
  int fgst=1;
  for(int i=0; i<nd; i++) {
    fgst = (fgst+1)%2;
    for(int j=0; j<nd; j++) {
      if (j%2==fgst)  continue;
      cnt++; 
      vd.push_back(Dish(cnt, i*Ddish, j*Ddish, Eta*Ddish));
    }
  }
  cout << ">>>CreateSemiFilledSqConfig(" << nd << "," << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
  
  return vd;
}

/* --Fonction -- */
vector<Dish> CreateConfigA(double Ddish, double Eta)
{
  vector<Dish>  vd;
  int cnt=0;
  for(int i=0; i<18; i++) {
    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, i*Ddish, 17.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<17)) {
      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,17.*Ddish,i*Ddish,Eta*Ddish));
    }
  }
  cout << ">>>CreateConfigA(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
  return vd;
}

/* --Fonction -- */
vector<Dish> CreateConfigB(double Ddish, double Eta)
{
  vector<Dish>  vd;
  int cnt=0;
  /*
  for(int i=0; i<13; i++) {
    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, i*Ddish, 12.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<12)) {
      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,12.*Ddish,i*Ddish,Eta*Ddish));
    }
  }
  for(int i=0; i<5; i++) {
    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish,4.*Ddish,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 8.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<4)) {
      cnt++; vd.push_back(Dish(cnt,4.*Ddish,(i+4)*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,8.*Ddish,(i+4)*Ddish,Eta*Ddish));
    }
  }
  */
  for(int i=0; i<11; i++) {
    cnt++; vd.push_back(Dish(cnt, i*Ddish,0.,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, i*Ddish, 10.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<10)) {
      cnt++; vd.push_back(Dish(cnt,0.,i*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,10.*Ddish,i*Ddish,Eta*Ddish));
    }
  }
  for(int i=0; i<7; i++) {
    cnt++; vd.push_back(Dish(cnt, (i+2)*Ddish, 2.*Ddish,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, (i+2)*Ddish, 8.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<6)) {
      cnt++; vd.push_back(Dish(cnt,2.*Ddish,(i+2)*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,8.*Ddish,(i+2)*Ddish,Eta*Ddish));
    }
  }
  for(int i=0; i<3; i++) {
    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 4.*Ddish,Eta*Ddish));
    cnt++; vd.push_back(Dish(cnt, (i+4)*Ddish, 6.*Ddish,Eta*Ddish));
    if ((i>0)&&(i<2)) {
      cnt++; vd.push_back(Dish(cnt,4.*Ddish,(i+4)*Ddish,Eta*Ddish));
      cnt++; vd.push_back(Dish(cnt,6.*Ddish,(i+4)*Ddish,Eta*Ddish));
    }
  }


  cout << ">>>CreateConfigB(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;
  return vd;
}


/* --Fonction -- */
vector<Dish> CreateConfigC(double Ddish, double Eta)
{
  vector<int_4> lesx, lesy;
  
  int max = 16;
  for(int i=0; i<4; i++) 
    for(int j=0; j<4; j++) {
      lesx.push_back(i); lesy.push_back(j);
      lesx.push_back(max-i); lesy.push_back(max-j);
      lesx.push_back(i); lesy.push_back(max-j);
      lesx.push_back(max-i); lesy.push_back(j);
    }
  
  for(int i=5; i<12; i+=2) {
    lesx.push_back(i); lesy.push_back(0);
    lesx.push_back(i); lesy.push_back(max);
    lesx.push_back(0); lesy.push_back(i);
    lesx.push_back(max); lesy.push_back(i);
  }
  
  for(int i=4; i<=12; i+=2) 
    for(int j=4; j<=12; j+=2) {
      lesx.push_back(i); lesy.push_back(j);
    }
  
  for(int i=5; i<=11; i+=2) {
    lesx.push_back(i); lesy.push_back(4);
    lesx.push_back(i); lesy.push_back(max-4);
    lesx.push_back(4); lesy.push_back(i);
    lesx.push_back(max-4); lesy.push_back(i);
  }
  
  EnumeratedSequence esx,esy;
  esx = 2,5;
  esy = 5,2;
  
  for(int k=0; k<esx.Size(); k++) {
    int_4 ix=esx.Value(k);
    int_4 iy=esy.Value(k);
    
    lesx.push_back(ix); lesy.push_back(iy);
    lesx.push_back(max-ix); lesy.push_back(iy);
    lesx.push_back(ix); lesy.push_back(max-iy);
    lesx.push_back(max-ix); lesy.push_back(max-iy);
  }
 cout << "CreateConfigC/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
 
 vector<Dish> vd;
 int cnt=0;
 for(size_t i=0; i<lesx.size(); i++) {
   cnt++; vd.push_back(Dish(cnt, ((double)lesx[i])*Ddish,((double)lesy[i])*Ddish,Eta*Ddish));
 }

 cout << ">>>CreateConfigC(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;

 return vd;

}
/* --Fonction -- */
vector<Dish> CreateConfigD(double Ddish, double Eta)
{
EnumeratedSequence es;
es = 0,1,3,4,6,8,10,12,14,16,18,20,22,24,25,27,28;
vector<int_4> lesx, lesy;
for(int k=0; k<es.Size(); k++) {
  lesx.push_back(es.Value(k));   lesy.push_back(0);
  lesx.push_back(es.Value(k));   lesy.push_back(28);
}
for(int k=1; k<es.Size()-1; k++) {
  lesy.push_back(es.Value(k));   lesx.push_back(0);
  lesy.push_back(es.Value(k));   lesx.push_back(28);
}
for(int k=1; k<=5; k++) {
  lesy.push_back(k);   lesx.push_back(5);
  lesy.push_back(28-k);   lesx.push_back(28-5);
  if (k!=5) { 
    lesx.push_back(k);   lesy.push_back(5);
    lesx.push_back(28-k);   lesy.push_back(28-5);
  }

  lesy.push_back(k);   lesx.push_back(28-5);
  lesy.push_back(28-k);   lesx.push_back(5);
  if (k!=5) { 
    lesx.push_back(28-k);   lesy.push_back(5);
    lesx.push_back(k);   lesy.push_back(28-5);
  }
}

for(int k=6; k<=13; k++) {
  lesy.push_back(k);   lesx.push_back(k);
  lesy.push_back(28-k);   lesx.push_back(28-k);
  lesy.push_back(k);   lesx.push_back(28-k);
  lesy.push_back(28-k);   lesx.push_back(k);
}


lesx.push_back(14); lesy.push_back(14);

EnumeratedSequence esx,esy;
esx = 1,2,4;
esy = 12,11,13;

for(int k=0; k<esx.Size(); k++) {
  int_4 ix=esx.Value(k);
  int_4 iy=esy.Value(k);
  lesx.push_back(ix); lesy.push_back(iy);
  lesx.push_back(28-ix); lesy.push_back(iy);
  lesx.push_back(ix); lesy.push_back(28-iy);
  lesx.push_back(28-ix); lesy.push_back(28-iy);
  
  lesy.push_back(ix); lesx.push_back(iy);
  lesy.push_back(28-ix); lesx.push_back(iy);
  lesy.push_back(ix); lesx.push_back(28-iy);
  lesy.push_back(28-ix); lesx.push_back(28-iy);
 }  
for(int k=5; k<=13; k+=2) {
  lesy.push_back(k);   lesx.push_back(14);
  lesy.push_back(28-k);   lesx.push_back(14);
  lesy.push_back(14);   lesx.push_back(k);
  lesy.push_back(14);   lesx.push_back(28-k);
}

 cout << "CreateConfigB/Debug: -checkSize/lesx=" << lesx.size() << " -Check/lesy=" << lesy.size() << endl;
 
 vector<Dish> vd;
 int cnt=0;
 for(size_t i=0; i<lesx.size(); i++) {
   cnt++; vd.push_back(Dish(cnt, ((double)lesx[i])*Ddish,((double)lesy[i])*Ddish,Eta*Ddish));
 }

 cout << ">>>CreateConfigD(" << Ddish << "," << Eta << ") ---> NDishes=" << vd.size() << endl;

 return vd;
}

/* --Fonction -- */
vector<Dish> CreateFilledCylConfig(int ncyl,  int nRL, double cylW, double cylRL, double etaW, double etaRL, bool fgscid)
{
  vector<Dish>  vd;
  int cnt=0;

  for(int i=0; i<ncyl; i++) 
    for(int j=0; j<nRL; j++) {
      cnt++;
      int rid = (fgscid) ? i+1 : cnt; 
      vd.push_back(Dish(rid, i*cylW, j*cylRL, etaW*cylW, etaRL*cylRL));
    }
  cout << ">>>CreateFilledCylConfig(" << ncyl << "," << nRL << "," << cylW << "," << cylRL << "," 
       << etaW << "," << etaRL << "," << ((fgscid)?" RId=CylNum":"Cnt") 
       << ") ---> NDishes=" << vd.size() << endl;

  return vd;
}