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

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

  Usage:  repicon [-parname value] configId OutPPFName 
---------------------------------------------------------------  */

#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 
// ---------------------------------------------------------------------

void Usage() 
{
  cout << " Usage: repicon [-parname Value] configId OutPPFName \n" 
       << " configIds: f4x4,f8x8,f20x20, confA,confB,confC,confD, hex12,cross11, f3cyl,f8cyl,f3cylp,f8cylp \n"
       << "   f4x4 , f8x8 , f20x20 Filled array of nxn dishes  \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 \n"
       << "   f4cylw : filled array of 4 perfect of wide cylinders \n"
       << "  [ -parname value] : -renmax -z -prt -D -lmax \n" 
       << "    -renmax MaxValue (default : Do NOT renormalize 2D response value \n"   
       << "    -z redshift (default=0.7) --> determines Lambda \n"
       << "    -D DishDiameter (default=5 m) \n"
       << "    -lmax array extension (default=100 m ) for response calculation kmax \n" 
       << "    -repnxy Nx,Ny (default=200,200) ResponseTable binning \n"
       << "    -beamnxy Nx,Ny (default=200,200) Beam sampling \n"
       << "    -rot ThetaMaxDeg,NTheta,PhiMaxDeg,NPhi (default NO-rotate/pointing -> 23,10,30,15 ) \n" 
       << "    -prt PrtLev,PrtModulo (default=0,10) \n" 
       << endl; 
  return;
}

//      ------------------ MAIN PROGRAM ------------------------------
int main(int narg, const char* arg[])
{
  if (((narg>1)&&(strcmp(arg[1],"-h")==0))||(narg<3))  {
    Usage();
    return 1;
  }
  // make sure SOPHYA modules are initialized 
  SophyaInit();  
  //  FitsIOServerInit();
  InitTim();
  //--- decoding command line arguments 
  string config = "f8x8" ;
  string outfile = "repicon.ppf";  
  bool fgrenorm=false;
  double rmax=1.;
  double z_Redshift=0.7 ;  // 21 cm at z=0.7 -> 0.357 m  
  int prtlev=0;
  int prtmod=10;
  
  double Ddish=5.;
  bool fgDfixed=false;
  double LMAX = 100.;  // taille de la zone, pour calcul kmax 
  bool fgLMAXfixed=false;
  double thetamxdeg=23.;  // 23 degres : angle d'inclinaison de l'orbite terrestre
  double phimxdeg=30.;
  int nteta=10;
  int nphi=15;
  bool fgpoint=false;
  int NRX=200;
  int NRY=200;
  int NBX=200;
  int NBY=200;
  
  int ka=1;
  while (ka<(narg-1)) {
    if (strcmp(arg[ka],"-renmax")==0) {
      rmax=atof(arg[ka+1]);  fgrenorm=true;   ka+=2;
    }
    else if (strcmp(arg[ka],"-z")==0) {
      z_Redshift=atof(arg[ka+1]);  ka+=2;
    }
    else if (strcmp(arg[ka],"-D")==0) {
      Ddish=atof(arg[ka+1]);  fgDfixed=true;  ka+=2;
    }
    else if (strcmp(arg[ka],"-lmax")==0) {
      LMAX=atof(arg[ka+1]);  fgLMAXfixed=true;  ka+=2;
    }
    else if (strcmp(arg[ka],"-repnxy")==0) {
      sscanf(arg[ka+1],"%d,%d",&NRX,&NRY);   ka+=2;
    }
    else if (strcmp(arg[ka],"-beamnxy")==0) {
      sscanf(arg[ka+1],"%d,%d",&NBX,&NBY);   ka+=2;
    }
    else if (strcmp(arg[ka],"-rot")==0) {
      sscanf(arg[ka+1],"%lg,%d,%lg,%d",&thetamxdeg,&nteta,&phimxdeg,&nphi);  fgpoint=true;  ka+=2;
    }
    else if (strcmp(arg[ka],"-prt")==0) {
      sscanf(arg[ka+1],"%d,%d",&prtlev,&prtmod);   ka+=2;
    }
    else break; 
  }

  if ((ka+1)>=narg) {
    cout << " repicon / Argument error " << endl;
    Usage();
    return 2;
  }

  config = arg[ka];
  outfile = arg[ka+1];
  //-- end command line arguments
  
  int rc = 1;  
  try {  // exception handling try bloc at top level
    cout << " ==== repicon.cc program , (u,v) plane response ==== " << endl;

    H21Conversions conv;
    double LAMBDA=0.357 ;  // 21 cm at z=0.7 
    conv.setRedshift(z_Redshift);
    LAMBDA = conv.getLambda();
    
    double Eta=0.95;
    double cylW=12.;   // Largeur des cylindres 
    double cylRL=2*LAMBDA;  // 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> vdishes;
    cout << " repicon[1] : creating dish vector/ MultiDish for configuration: " << config << endl;

    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=2*LAMBDA; 
      vdishes=CreateFilledCylConfig(3, 92, cylW, cylRL, etaW, etaRL, true);
    }
    else if (config=="f3cylp") {
      cylW=10.;   cylRL=2*LAMBDA; 
      vdishes=CreateFilledCylConfig(3, 92, cylW, cylRL, etaW, etaRL, false);
    }
    else if (config=="f4cylw") {
      cylW=25.;   cylRL=3*LAMBDA; 
      vdishes=CreateFilledCylConfig(4, 100, cylW, cylRL, etaW, etaRL, false);
    }
    else if (config=="f8cyl") {
      cylW=12.;   cylRL=2*LAMBDA; 
      vdishes=CreateFilledCylConfig(8, 160, cylW, cylRL, etaW, etaRL, true);
    }
    else if (config=="f8cylp") {
      cylW=12.;   cylRL=2*LAMBDA; 
      vdishes=CreateFilledCylConfig(8, 160, cylW, cylRL, etaW, etaRL, false);
    }
    else if (config=="confA") {
      vdishes=CreateConfigA(Ddish, Eta);
    }
    else if (config=="confB") {
      vdishes=CreateConfigB(Ddish, Eta);
    }
    else if (config=="confC") {
      vdishes=CreateConfigC(Ddish, Eta);
    }
    else if (config=="confD") {
      vdishes=CreateConfigD(Ddish, Eta);
    }
    else if (config=="cross11") {
      if (!fgDfixed) Ddish = 12.;
      double base=20.;
      Eta=0.95;
      if (!fgLMAXfixed) LMAX = 250.;
      vdishes=CreateCrossConfig(Ddish,base,Eta);
    }
    else if (config=="hex12") {
      if (!fgDfixed) Ddish = 12.;
      Eta=0.95;
      if (!fgLMAXfixed) LMAX = 350.;
      vdishes=CreateDoubleHexagonConfig(Ddish);
    }
    else {
      cout << " NON valid configuration option -> exit" << endl;
      return 99;
    }
        
    double Dol = LMAX/LAMBDA;
    bool fgnoauto = true;
    
    cout << " repicon[1] : Lambda=" << LAMBDA << " LMAX= " << LMAX << " NDishes=" << vdishes.size() 
         << " D-Dish=" << Ddish << " m." << endl;
 
    MultiDish mdish(LAMBDA, LMAX, vdishes, fgnoauto);
    mdish.SetPrtLevel(prtlev,prtmod);
    mdish.SetRespHisNBins(NRX,NRY); 
    mdish.SetBeamNSamples(NBX,NBY);
    if (fgpoint)  {
      cout << " repicon[1.b] : activating pointing , ThetaMaxDeg=" << thetamxdeg << " NTheta=" << nteta
           << "  PhiMaxDeg= " << phimxdeg << " NPhi=" << nphi << endl;
      mdish.SetThetaPhiRange(Angle(thetamxdeg,Angle::Degree).ToRadian(),nteta,Angle(phimxdeg,Angle::Degree).ToRadian(), nphi);
    }
    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;
}