source: Sophya/trunk/Cosmo/RadioBeam/repicon.cc@ 3931

/*  ------------------------ Projet BAORadio -------------------- 
  Calcul de la reponse 2D (plan (u,v) d'un interferometre 
    R. Ansari , C. Magneville - Juin-Dec 2010 

  Usage:  repicon configId OutPPFName [z_Redshift=0.7] [RenormalizeMax]
---------------------------------------------------------------  */

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

#include <typeinfo>

#include "mdish.h"
#include "specpk.h"
#include "interfconfigs.h"
#include "radutil.h"

#include "ntuple.h"
#include "histinit.h"
// #include "fiosinit.h"
// #include "fitsioserver.h"

#include "randr48.h"      

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

typedef DR48RandGen RandomGenerator ;

// pour sauver la reponse mdresp et la config des dishes dans un fichier PPF 
void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp);

// ---------------------------------------------------------------------
// main program for computing interferometer response un (u,v) plane
// R. Ansari  - Avril-Juin 2010 
// ---------------------------------------------------------------------

//      ------------------ MAIN PROGRAM ------------------------------
int main(int narg, const char* arg[])
{
  if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
    cout << " Usage: repicon configId OutPPFName [z_redshift=0.7] [RenormalizeMax] \n" 
         << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes (D=5m) \n" 
         << "          confA , confB, confC, confD : semi-filled array of dishes \n"
         << "          hex12,cross11 : ASKAP like double hexagonal (12xD=12m), cross config (11xD=12m) \n"
         << "          f3cyl, f8cyl , f3cylp, f8cylp : filled array of non perfect/perfect of n cylinders " << endl; 
    return 1;
  }
  cout << " ==== repicon.cc program , (u,v) plane response ==== " << endl;
  // make sure SOPHYA modules are initialized 
  SophyaInit();  
  //  FitsIOServerInit();
  InitTim();
  //--- decoding command line arguments 
  string config = "f8x8" ;
  if (narg>1) config = arg[1];
  string outfile = "repicon.ppf";  
  if (narg>2) outfile = arg[2];
  if (outfile==".")  outfile = "repicon.ppf";
  double LAMBDA=0.357 ;  // 21 cm at z=0.7 
  if (narg>3) LAMBDA = atof(arg[3]);
  bool fgrenorm=false;
  double rmax=1.;
  if (narg>4) {
    rmax=atof(arg[4]);
    fgrenorm=true;
  }
  //-- end command line arguments
  
  int rc = 1;  
  try {  // exception handling try bloc at top level

    double Ddish = 5.;
    double Ddish2 = 7.5;
    double Eta=0.95;
    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

    double D = 100.;    // Taille de la zone

    int cnt=0;

    vector<Dish> vdishes;
    cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config << endl;

    bool fgpoint=false;
    if (config=="f4x4") {
      vdishes=CreateFilledSqConfig(4,Ddish, Eta);
    }
    else if (config=="f8x8") {
      vdishes=CreateFilledSqConfig(8,Ddish, Eta);
    }
    else if (config=="f20x20") {
      vdishes=CreateFilledSqConfig(20,Ddish, Eta);
    }
    else if (config=="f3cyl") {
      cylW=10.;   cylRL=0.5; 
      vdishes=CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, true);
    }
    else if (config=="f3cylp") {
      cylW=10.;   cylRL=0.5; 
      vdishes=CreateFilledCylConfig(3, 128, cylW, cylRL, etaW, etaRL, false);
    }
    else if (config=="f8cyl") {
      cylW=12.;   cylRL=0.5; 
      vdishes=CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, true);
    }
    else if (config=="f8cylp") {
      cylW=12.;   cylRL=0.5; 
      vdishes=CreateFilledCylConfig(8, 192, cylW, cylRL, etaW, etaRL, false);
    }
    else if (config=="confA") {
      fgpoint=true;
      vdishes=CreateConfigA(Ddish, Eta);
    }
    else if (config=="confB") {
      fgpoint=true;
      vdishes=CreateConfigB(Ddish, Eta);
    }
    else if (config=="confC") {
      fgpoint=true;
      vdishes=CreateConfigC(Ddish, Eta);
    }
    else if (config=="confD") {
      fgpoint=true;
      vdishes=CreateConfigD(Ddish, Eta);
    }
    else if (config=="cross11") {
      fgpoint=true;
      Ddish = 12.;
      double base=20.;
      Eta=0.95;
      D=250.;
      vdishes=CreateCrossConfig(Ddish,base,Eta);
    }
    else if (config=="hex12") {
      fgpoint=true;
      Ddish = 12.;
      Eta=0.95;
      D=350.;
      vdishes=CreateDoubleHexagonConfig(Ddish);
    }

    else {
      cout << " NON valid configuration option -> exit" << endl;
      return 99;
    }

    double Dol = D/LAMBDA;
    double LMAX = D;
    bool fgnoauto = true;
    int NRX=200;
    int NRY=200;

    MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
    mdish.SetRespHisNBins(NRX,NRY); 
    if (fgpoint)  // 23 degres : angle d'inclinaison de l'orbite terrestre
      mdish.SetThetaPhiRange(Angle(23.,Angle::Degree).ToRadian(),10, M_PI/6., 15);
    cout << " repicon[2] : calling mdish.GetResponse() ..."<< endl;

    Histo2D hrep = mdish.GetResponse();
    PrtTim("Apres mdish.GetResponse()");

    Four2DRespTable mdresp(hrep, Dol, LAMBDA);
    if (fgrenorm) {
      cout << " repicon[2.b] call to mdresp.renormalize(" << rmax << ")"; 
      double omax=mdresp.renormalize(rmax);
      cout << " Old Max=" << omax << endl;
    }
    cout << " repicon[3] : saving Four2DRespTable for config " << config << " to PPF file " << outfile << endl;
    mdresp.writeToPPF(outfile);

    string outfile2 = "hdt_"+outfile;
    cout << " repicon[4] : saving H2D-response, multidish config to PPF file " << outfile2 << endl;
    SaveDTVecDishesH2Resp(outfile2, vdishes, mdresp);

    rc = 0;
  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " repicon.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 << " repicon.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " repicon.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  PrtTim("End-repicon");
  cout << " ==== End of repicon.cc program  Rc= " << rc << endl;
  return rc;    
}

/*-- Nouvelle-Fonction --*/
void SaveDTVecDishesH2Resp(string& outfile, vector<Dish>& vdishes, Four2DRespTable& mdresp)
{
  char* names[5]={"did","posx","posy","diam","diamy"};
  NTuple ntvd(5,names,64,false);
  r_4 xnt[10];
  for(size_t i=0; i<vdishes.size(); i++) {
    xnt[0]=vdishes[i].ReflectorId();
    xnt[1]=vdishes[i].X;
    xnt[2]=vdishes[i].Y;
    if (vdishes[i].isCircular()) {
      xnt[3]=vdishes[i].Diameter();
      xnt[4]=0.;
    }
    else {
      xnt[3]=vdishes[i].DiameterX();
      xnt[4]=vdishes[i].DiameterY();      
    }
    ntvd.Fill(xnt);
  }
  Histo2D h2rep=mdresp.GetResponse();
  POutPersist po(outfile);
  po << PPFNameTag("mdish") << ntvd;
  po << PPFNameTag("h2rep") << h2rep;
  return;
}