source: Sophya/trunk/SophyaProg/PMixer/tgrsr.cc@ 2250

#include <typeinfo>
#include "pmixer.h"

/*!\ingroup PMixer
 * \file tgrsr.cc
 * \brief \b PROGRAM \b tgrsr  <BR>
 * Program to generate different type of RadSpectra
 * and SpectralResponse in the form of RadSpectraVec and SpecRespVec
 * in FITS format
 */


// ------------- Main program --------------
int main(int narg, char* arg[]) 
{
#if !defined(SunOS) || !defined(__GNUG__)
  // ca se plante sur Sun !
  SophyaInit();
#endif

  InitTim();   // Initializing the CPU timer 

  if ((narg < 5) || ((narg>1) && (strcmp(arg[1],"-h") == 0)) )  {
    cout << " tgrsr : Generation of RadSpectraVec and SpecRespVec in FITS  " << endl;
    cout << " Usage: tgrsr 0/1 Params NPoints FitsFileName [ImageR4FITSName] [ppfname] " << endl;
    cout << " 0 -> GaussianFilter  T = A exp(((nu-nu0)/dnu)^2" << endl;
    cout << " Params : A,Nu0,DNu,MinFreq,MaxFreq " << endl;
    cout << " 1 -> PowerLawSpectra f = a ((nu-nu0)/dnu)^b " << endl;
    cout << " Params : A,Nu0,dnu,b,MinFreq,MaxFreq " << endl;
    exit(0);
    }

try {
  int typ = 0; 
  typ = atoi(arg[1]);
  int k, npt = 100;
  npt = atoi(arg[3]);
  if (npt < 5) npt = 5;
  cout << " ---- tgrsr/Info: Type= " << typ << " NPoints= " << npt << endl;
  Matrix mtx(2, npt);
  ImageR4 img(npt, 2);

  char buff[256];
  if (typ == 0) {  // Generation de SpecRespVec from GaussianFilter
    double a,nu0,dnu,fmin,fmax;
    a = 1.0;
    nu0 = 143.;
    dnu = 10.;
    fmin = 10.;
    fmax = 1000.;
    sscanf(arg[2],"%lg,%lg,%lg,%lg,%lg",&a,&nu0,&dnu,&fmin,&fmax);
    cout << " GaussianFilter  T = A exp(((nu-nu0)/dnu)^2" << endl;
    sprintf(buff,"A=%g Nu0=%g DNu=%g FMin=%g FMax=%g",
            a,nu0,dnu,fmin,fmax);
    cout << "Decoded params: " << (string)buff << endl;
    GaussianFilter sr(nu0, dnu, a, fmin, fmax);
    cout << sr << endl;
    cout << " Creating SpecRespVec() ... " << endl;
    Vector v1(10);
    Vector v2(10);
    SpecRespVec srv(v1, v2, fmin, fmax);
    Vector & nu = srv.getNuVec();
    Vector & tnu = srv.getTNuVec();
    nu.ReSize(npt);
    tnu.ReSize(npt);
    double freq;
    double dfreq = (fmax-fmin)/(npt-3);
    for(k=0; k<npt; k++) {
      freq = fmin+(k-1)*dfreq;
      nu(k) = freq;
      tnu(k) = sr.transmission(freq);
      //      if ((k>10) && (k<25))  cout << " DBG - K = " << k <<
      //        " Nu= " << nu(k) << " Tnu= " << tnu(k) << endl; 
      mtx(0,k) = nu(k);
      mtx(1,k) = tnu(k);
      img(k,0) = nu(k);
      img(k,1) = tnu(k);
    }
    if ((typ == 0) && (narg > 6)) { // Writing to ppf  file
      cout << "writing SpecRespVec(" << npt << ") to PPF file: " 
           << (string)(arg[6]) << endl;
      cout << sr ; 
      cout << srv ; 
      {
        POutPersist pos(arg[6]);
        pos << srv;
      }
      PrtTim("End of ImageR4->PPF ");
    }

  }
  else { // Generation de RadSpectraVec from PowerLawSpectra
    double a,nu0,dnu,b,fmin,fmax;
    a = 1.0;
    nu0 = 143.;
    b = 0.;
    dnu = 10.;
    fmin = 10.;
    fmax = 1000.;
    sscanf(arg[2],"%lg,%lg,%lg,%lg,%lg,%lg",&a,&nu0,&dnu,&b,&fmin,&fmax);
    cout << " PowerLawSpectra f = a ((nu-nu0)/dnu)^b " << endl;
    sprintf(buff,"A=%g Nu0=%g DNu=%g B=%g FMin=%g FMax=%g",
            a,nu0,dnu,b,fmin,fmax);
    cout << "Decoded params: " << (string)buff << endl;
    PowerLawSpectra rs(a,b, nu0, dnu, fmin, fmax);
    cout << rs << endl;
    cout << " Creating RadSpectraVec() ... " << endl;
    Vector v1(10);
    Vector v2(10);
    RadSpectraVec rsv(v1, v2, fmin, fmax);
    Vector & nu = rsv.getNuVec();
    Vector & fnu = rsv.getFNuVec();
    nu.ReSize(npt);
    fnu.ReSize(npt);
    double freq;
    double dfreq = (fmax-fmin)/(npt-3);
    for(k=0; k<npt; k++) {
      freq = fmin+(k-1)*dfreq;
      nu(k) = freq;
      fnu(k) = rs.flux(freq);
      //      if ((k>10) && (k<25))  cout << " DBG2 - K = " << k <<
      //        " Nu= " << nu(k) << " Fnu= " << fnu(k) << endl; 
      mtx(0,k) = nu(k);
      mtx(1,k) = fnu(k);
      img(k,0) = nu(k);
      img(k,1) = fnu(k);
    }
  }

  PrtTim("End of filling ");

  {
  cout << "\n Writing Matrix NRows= " << mtx.NRows() << "  NCols= " 
       << mtx.NCols() << " to FITS file: " << (string)(arg[4]) << endl;
  FitsOutFile fios(arg[4]);
  fios.firstImageOnPrimaryHeader(false);  // use secondary header
  DVList dvl;
  if (typ == 0) 
    dvl["PDMTYPE"] = "FGRSPEC";
  dvl.SetComment("PDMTYPE", "Planck Data Model Type"); 
  dvl["SOPHYVER"] =  SophyaVersion();
  dvl.SetComment("SOPHYVER", "Sophya Version"); 
  dvl["TGRSRKW"] =  "Test Keyword from tgrsr.cc";
  fios.DVListIntoPrimaryHeader(dvl);  
  fios << mtx ; 
  PrtTim("End of Matrix->FITS ");
  }

  if (narg > 5) { // Writing ImageR4 fo FITS file
    cout << "writing ImageR4(" << npt << ", 2) to FITS file: " 
         << (string)(arg[5]) << endl;
    FitsOutFile fios(arg[5]);
    fios.firstImageOnPrimaryHeader(true);  // use primary  header
    fios << img ; 
    PrtTim("End of ImageR4->FITS ");
  }
}
catch (PThrowable & exc) {
  cerr << "tgrsr/Error Catched Exception " << (string)typeid(exc).name() 
       << " - Msg= " << exc.Msg() << endl;
}
catch (...) {
  cerr << "tgrsr/Error some other exception was caught ! " << endl;
}

  cout <<  " ============ End of tgrsr program ======== " << endl;
  return 0;
}