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

#include "machdefs.h"
#include <iostream.h>
#include <typeinfo>
#include <math.h>

#include "tvector.h"
#include "cimage.h"
#include "sambainit.h"

#include "radspecvector.h"
#include "specrespvector.h"
#include "nupower.h"
#include "gaussfilt.h"

#include "fitsioserver.h"

#include "timing.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
*/

// Test 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[]) 
{
  SophyaInit();
  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);
    }

  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 ");

  FitsIoServer fios;
  cout << "\n Writing Matrix NRows= " << mtx.NRows() << "  NCols= " 
       << mtx.NCols() << " to FITS file: " << (string)(arg[4]) << endl;
  fios.save(mtx, arg[4]);
  PrtTim("End of Matrix->FITS ");

  if (narg > 5) { // Writing ImageR4 fo FITS file
    cout << "wrting ImageR4(" << npt << ", 2) to FITS file: " 
         << (string)(arg[5]) << endl;
    fios.save(img, arg[5]);
    PrtTim("End of ImageR4->FITS ");
  }
  cout <<  " ============ End of tgrsr program ======== " << endl;
  return 0;
}