source: Sophya/trunk/Cosmo/RadioBeam/syncube.cc@ 3783

/*  ------------------------ Projet BAORadio -------------------- 
  Programme d'extraction d'une partie de carte synchrotron 
  (HASLAM @ 400 MHz) et fabrication d'un cube 3D (angles,fre) 
    R. Ansari , C. Magneville - Juin 2010 

  Usage: syncube InFitsName Out3DPPFName [Out2DMapName] 
---------------------------------------------------------------  */

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

#include "tvector.h"
#include "srandgen.h"
#include "fioarr.h"
#include "sopemtx.h"
#include "pexceptions.h"

#include "randr48.h"

#include "tvector.h"        // Pour l'utilisation des classes TArray, TMatrix , TVector
#include "matharr.h"

/*
#include "spherehealpix.h"   // Pour les cartes spheriques pixelisees au format HEALPix 
#include "spherethetaphi.h"   // Pour les cartes spheriques pixelisees au format Theta-Phi 
#include "localmap.h"         // Pour les cartes locales 
#include "mapoperation.h"         // Pour les cartes locales 
*/
#include "skymap.h"
#include "mapoperation.h"         // Pour les cartes locales 

#include "fitsspherehealpix.h"   // Pour les I/O fits de HEALPix
#include "fitsspherethetaphi.h"   // Pour les I/O fits de SphereThetaPhi
#include "fitslocalmap.h"         // Pour les I/O fits de LocalMap<T> 

#include "xastropack.h"     // Pour faire les conversions de coordonnees celestes

// Pour l'initialisation des modules 
#include "tarrinit.h"
#include "skymapinit.h"     
#include "fiosinit.h"     


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

//----------------------------------------------------------------------------
//----------------------------------------------------------------------------   
int main(int narg, char* arg[])
{
  // Sophya modules initialization
  TArrayInitiator  _inia;
  SkyMapInitiator  _inis;
  FitsIOServerInitiator  _inif;
  //------- AU LIEU DE ------>  SophyaInit();  

  InitTim();   // Initializing the CPU timer
  Timer tm("syncube");

  if (narg < 3) {
    cout << "Usage: syncube InFitsName Out3DPPFName [Out2DMapName] \n" << endl;
    return 1;
  }

  //---  Parametres des lois de puissance en frequence 
  double AmpPL1 = 1.;   // amp max PowerLaw 2 
  double PLidx1 = -2.5;  // index de la loi de puissance synchrotron 
  double sigPLidx1 = 0.1;  // Sigma de la variation (gaussienne) de index1 
  // Amplitude max de la 2eme composante en loi de puissance (tirage plat 0 ... AmpPL2)
  double AmpPL2 = 0.1;   // amp max PowerLaw 2 
  double PLidx2 = -3.2; 
  double sigPLidx2 = 0.15;

  // decodage arguments 
  string outname = arg[2];
  string inname = arg[1];
  int rc = 91;

  cout << " ====== syncube :   Input map name= " << inname << " OutName=" << outname;
  bool fginmap=true;
  try {
    SphereHEALPix<r_4> inmap0;
    if (fginmap) {  // Lecture eventuelle de la carte en entree 
      cout << " syncube[0]: Reading input map : " << inname << endl;
      FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
      fis >> inmap0;
      cout << inmap0;
    }
    // On fait une carte en doublant la resolution 
    int nside = inmap0.SizeIndex() ;  // Parametre de pixelisation HEALPix - definit la resolution
    SphereHEALPix<r_4> inmap(2*nside);
    for(int_4 kk=0; kk<inmap.NbPixels(); kk++) {
      double theta, phi;  // Theta, Phi en radians
      inmap.PixThetaPhi(kk, theta, phi);
      inmap(kk) = inmap0(theta, phi);
    }    
    //    Sph2Sph(inmap0, inmap);
    cout << "syncube[1]: Input resolution doubled : " << inmap << endl;

    sa_size_t NTet=256;
    sa_size_t NPhi=256;
    LocalMap<r_4> outmap(NPhi,NTet,60.,60.,110.,150.);
    for(int_4 kk=0; kk<outmap.NbPixels(); kk++) {
      double theta, phi;  // Theta, Phi en radians
      outmap.PixThetaPhi(kk, theta, phi);
      outmap(kk) = inmap(theta, phi);
    }    

    TArray<r_4> omap(NPhi,NTet);
    double tet0 = Angle(60.,Angle::Degree).ToRadian();
    double phi0 = Angle(120.,Angle::Degree).ToRadian();
    double dtet = Angle(60.,Angle::Degree).ToRadian()/(double)NTet;
    double dphi = Angle(60.,Angle::Degree).ToRadian()/(double)NPhi;
    for (sa_size_t j=0; j<omap.SizeY(); j++)  {
      double theta = j*dtet+tet0;
      for (sa_size_t i=0; i<omap.SizeX(); i++)  {
        double phi = i*dphi+phi0;
        omap(i,j) = inmap(theta, phi);
      }
    }
    //    Sph2Sph(inmap,outmap);
    
    cout << "syncube[2]: Output rectangular map computed " << outmap << endl;
    if (narg > 3) {
      string ppfname = arg[3];
      cout << " syncube[3]: Saving inmap/outmap tp PPF file-> " << ppfname << endl;
      POutPersist po(ppfname);
      po << PPFNameTag("inmap") << inmap;
      po << PPFNameTag("outmap") << outmap;
      po << PPFNameTag("omap") << omap;
    }

    sa_size_t NFreq = 128;
    TArray<r_4> ocube(NPhi,NTet,NFreq);

    double infreq = 400.; // Frequence carte input en MHz
    double freq0 = 840.;  // Freq0 du cube de sortie 
    double dfreq = 1.;  

    ThSDR48RandGen rg;
    for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
      for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
        double freqexpo = rg.Gaussian(sigPLidx1,PLidx1);
        for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
          double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
          ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
        }
        // On ajoute une autre composante avec un indice spectral different
        freqexpo = rg.Gaussian(sigPLidx2,PLidx2); 
        double famp = rg.Flat01()*AmpPL2;
        for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
          double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
          ocube(i,j,k) += famp*omap(i,j)*rapfreq;
        }

      }
    }
    
    // On sauve le cube de sortie
    {
      cout << " syncube[2]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << ocube;
    }

    rc = 0;
  }
  catch (PThrowable& exc) {
    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n syncube.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " syncube.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

  cout << ">>>> syncube[9] ------- FIN ----------- Rc=" << rc << endl;
  return rc;
}