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

/*  ------------------------ 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"

#include "cubedef.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;
  }

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

   LocalMap<r_4> outmap(NPhi,NTheta,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,NTheta);
    double tet0 = Angle(Theta0Degre,Angle::Degree).ToRadian();
    double phi0 = Angle(Phi0Degre,Angle::Degree).ToRadian();
    double dtet = Angle(ThetaSizeDegre,Angle::Degree).ToRadian()/(double)NTheta;
    double dphi = Angle(PhiSizeDegre,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);
      }
    }

    double mean, sigma;
    MeanSigma(omap, mean, sigma);
    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << " Sizes:" << endl;
    omap.Show(); 

    //    Sph2Sph(inmap,outmap);
    
    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;
    }

    TArray<r_4> ocube(NPhi,NTheta,NFreq);

    double infreq = 400.; // Frequence carte input en MHz
    double freq0 = Freq0MHz;  // Freq0 du cube de sortie 
    double dfreq = FreqSizeMHz/(double)NFreq;  

    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[4]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << ocube;
    }

    rc = 0;
  }
  catch (PThrowable& exc) {
    cerr << " syncube.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;
}