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

#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 "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 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 [Lambda=0.357] [RenormalizeMax]" 
         << " configs: f4x4 , f8x8 , f20x20 Filled array of nxn dishes \n" 
         << "          confA , confB, confC : semi-filled array of dishes \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
    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 {
      cout << " NON valid configuration option -> exit" << endl;
      return 99;
    }

    double D = 100.;
    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) 
      mdish.SetThetaPhiRange(M_PI/6.,12, M_PI/6., 12);
    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);

    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;    
}