source: Sophya/trunk/Cosmo/RadioBeam/gsm2cube.cc@ 4065

/*  ------------------------ 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("gsm2cube");

  if (narg < 5) {
    cout << "Usage: gsm2cube GSMPPFDirectoryFileName GSMMapNumStart GSMMapNumEnd  Out3DPPFName  \n" 
         << "  GSMPPFDirectoryFileName in form of Directory/gsmJJJJJ "  << endl;
    return 1;
  }

  // decodage arguments 
  string inpath = arg[1];
  int numapstart = atoi(arg[2]);
  int numapend = atoi(arg[3]);
  string outname = arg[4];

  int rc = 91;

  cout << " ====== gsm2cube :   Input GSM map path= " << inpath << " NumStart=" << numapstart 
       << " End=" << numapend << " OutName=" << outname << endl;
  try {
    int nfreqgsm = numapend-numapstart+1;
    if (nfreqgsm!=NFreq) {
      cout << " gsm2cube/ERROR: (nfreqgsm=" << nfreqgsm << ") <> NFreq (=" << NFreq << ") -> exit 9" << endl;
      return 9;
    }

    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;

    TArray<r_4> ocube(NPhi,NTheta,NFreq);
    double mjd2000 = MJDfrDate(1,1,2000);   // Modified Julian Date pour 1 Janvier 2000

    cout << "gsm2cube[1] Loop over input GSM spherical maps ... " << endl;

    for(sa_size_t kf=0; kf<NFreq; kf++) {
      char mapname[512];
      sprintf(mapname,"%s_%d.ppf",inpath.c_str(),(int)(kf+numapstart));
      cout << kf << "- Reading GSM map from file " << mapname << " ... " ;
      PInPersist pin(mapname);
      SphereHEALPix<r_4> gsm;
      pin >> gsm;
      cout << " Freq=" << gsm.Info()["FMHz"] << " MHz " << endl;
      for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
        for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
          double theta = j*dtet+tet0;
          double phi = i*dphi+phi0;
          double ra, dec;  // Angle horaire (right ascension) [0...24], dec [-90 90]
          ra=Angle(phi).ToDegree()/360.*24.;
          dec=90-Angle(theta).ToDegree();
          double glong, glat;  // Longitude et latitude galactique - en degres [0...360], [-90...90]
          EqtoGal(mjd2000,ra,dec,&glong,&glat);
          phi=Angle(glong,Angle::Degree).ToRadian();
          theta=Angle(90-glat,Angle::Degree).ToRadian();
          ocube(i,j,kf)=gsm(theta,phi);
        }
      }
    }

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

    
    // On sauve le cube de sortie
    {
      cout << " gsm2cube[3]: Saving output cube to -> " << outname << endl;
      POutPersist poc(outname);
      poc << ocube;
    }
    rc = 0;
  }
  catch (PThrowable& exc) {
    cerr << " gsm2cube.cc catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n gsm2cube.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " gsm2cube.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

  cout << ">>>> gsm2cube ------- FIN ----------- Rc=" << rc << endl;
  return rc;
}