source: Sophya/trunk/SophyaLib/SkyT/skymixer.cc@ 607

#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>

#include <iostream.h>
#include <math.h>

#include <string>
#include <vector>

#include "timing.h"
#include "sambainit.h"
#include "pexceptions.h"
#include "datacards.h"
#include "fitsioserver.h"

#include "spheregorski.h"

#include "radspecvector.h"
#include "blackbody.h"
#include "nupower.h"

#include "squarefilt.h"
#include "trianglefilt.h"
#include "specrespvector.h"
#include "gaussfilt.h"

// ------ Function declaration 
void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K=1.);
void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K=1.);

int CheckCards(DataCards & dc, string & msg);
char * BuildFITSFileName(string const & fname);
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);

//  ----- Variable globale ------------
static char mapPath[256];     // Path for input maps
static int  hp_nside = 32;    // HealPix NSide
static int  nskycomp = 0;     // Number of sky components
static int  debuglev = 0;     // debuglevel

// -------------------------------------------------------------------------
//                             main program 
// -------------------------------------------------------------------------
int main(int narg, char * arg[])
{
  if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
    cout << "  Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
    exit(0);
  }

InitTim();

string msg;
int rc = 0;
POutPersist * so = NULL;

try {
  string dcard = arg[1];
  cout << " Decoding parameters from file " << dcard << endl;
  DataCards dc(dcard);

  rc = CheckCards(dc, msg);
  if (rc) goto Fin;

  cout << " skymix/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
  cout << "  ... MapPath = " << (string)mapPath << "  DbgLev= " << debuglev << endl;

// We create an output persist file for writing debug objects
  if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");

  SphereGorski<float> outgs(hp_nside);
  cout << " Output Gorski Map  created - NbPixels= " << outgs.NbPixels() << endl;
  outgs.SetPixels(0.);

  // Decoding detection pass-band filter 
  double nu0 = dc.DParam("FILTER", 0, 10.);
  double s = dc.DParam("FILTER", 1, 1.);
  double a = dc.DParam("FILTER", 2, 1.);
  double numin = dc.DParam("FILTER", 3, 0.1);
  double numax = dc.DParam("FILTER", 4, 9999);
  GaussianFilter filt(nu0, s, a, numin, numax);
  cout << " Filter decoded - Created " << endl;
  cout << filt << endl;

// FOR debug
  if (debuglev > 0)  SpectralResponse2Nt(filt, *so, "filter");

  PrtTim(" After FilterCreation ");

  SphereGorski<float> * ings = NULL;  // Our input map  
  FitsIoServer fios; // Our FITS IO Server
  char * flnm, buff[64];
  string key;

  // Loop over sky component 
  int sk;
  for(sk = 0; sk<nskycomp; sk++) {
    cout << " Processing sky component No " << sk+1 << endl;
    if (ings) { delete ings;  ings = NULL; }
    ings = new SphereGorski<float>(hp_nside);
    sprintf(buff, "MAPFITSFILE%d", sk+1);
    key = buff;
    flnm = BuildFITSFileName(dc.SParam(key, 0));
    cout << " Reading Input FITS map " << (string)flnm << endl;
    fios.load(*ings, flnm, 1);

    if (debuglev > 4) {  // Writing tne input map to the outppf
      FIO_SphereGorski<float> fiog(ings);
      fiog.Write(*so, key);
      }  
    }
  }   // End of try block

catch (PException exc) {
  msg = exc.Msg();
  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
  rc = 90;
  }   


Fin:
if (so) delete so;  //  Closing the debug ppf file 
if (rc == 0)  return(0);
cerr << " Error condition -> Rc= " << rc << endl;
cerr << " Msg= " << msg << endl;
return(rc);
}

/* Nouvelle-Fonction */
int CheckCards(DataCards & dc, string & msg)
//   Function to check datacards
{
mapPath[0] = '\0';
hp_nside = 32;
nskycomp = 0;
debuglev  = 0;

int rc = 0;
  //  Cheking datacards
  if ( (!dc.HasKey("SKYMIX")) || (!dc.HasKey("FILTER")) ) {
     rc = 71; 
     msg = "Invalid parameters  - NO @SKYMIX or @FILTER card ";
     return(rc);
  }

  // Decoding number of component and pixelisation parameter
  int mg = 32;
  int ncomp = 0;
  ncomp = dc.IParam("SKYMIX", 0, 0);
  mg = dc.IParam("SKYMIX", 1, 32);
  if (ncomp < 1)  {
    msg = "Invalid parameters  - Check datacards @SKYMIX ";
    rc = 72;
    return(rc);
  }

//  Checking input FITS file specifications 
  int kc;
  string key, key2;
  char buff[256];
  bool pb = false;
  for(kc=0; kc<ncomp; kc++) {
    sprintf(buff, "MAPFITSFILE%d", kc+1);
    key = buff;
    if (dc.NbParam(key) < 1) {   
      msg = "Missing or invalid card : " + key;
      pb = true;  break;
    }
    sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
    key = buff;
    sprintf(buff, "SPECTRAFUNC%d", kc+1);
    key2 = buff;
    if ( (dc.NbParam(key) < 1) && (dc.NbParam(key2) < 2) )  {   
      msg = "Missing card or invalid parameters : " + key + " or " + key2;
      pb = true;  break;
    }

    }

  if (pb)  {
    rc = 72;
    return(72);
  }

//  Checking detection filter specification 
  key = "FILTER";
  if (dc.NbParam(key) < 3) { 
    msg = "Missing card or invalid parameters : " + key;
    rc = 73;  return(rc); 
    }

// Initialiazing parameters
  rc = 0;
  msg = "OK";
  nskycomp = ncomp;
  hp_nside = mg;  

// Checking for PATH definition card
  key = "MAPPATH";
  if (dc.NbParam(key) < 3)  strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
  mapPath[255] = '\0';
  key = "DEBUGLEVEL";
  debuglev =  dc.IParam(key, 0, 0);
  return(rc);
}

static char buff_flnm[1024];   // Mal protege !
/* Nouvelle-Fonction */
char* BuildFITSFileName(string const & fname)
{
if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
else sprintf(buff_flnm, "%s", fname.c_str());
return(buff_flnm);
}

/* Nouvelle-Fonction */
// template <class T> 
void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap, 
                  PixelMap<float>& mapToAdd, RadSpectra& rs, double K)
{
  // finalMap = finalMap + coeff* mapToAdd
  // coeff    =  convolution of sr and rs 
  // compute the coefficient corresponding to mapToAdd
  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
  double coeff = rs.filteredIntegratedFlux(sr) * K;
  for(int i=0; i<finalMap.NbPixels(); i++)
    {
      finalMap(i) += coeff * mapToAdd(i);
    }
}

/* Nouvelle-Fonction */
void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap, 
                  PixelMap<double>& mapToAdd, RadSpectra& rs, double K)
{
  // finalMap = finalMap + coeff* mapToAdd
  // coeff    =  convolution of sr and rs 
  // compute the coefficient corresponding to mapToAdd
  if (finalMap.NbPixels() != mapToAdd.NbPixels())   
    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
  double coeff = rs.filteredIntegratedFlux(sr) * K;
  for(int i=0; i<finalMap.NbPixels(); i++)
    {
      finalMap(i) += coeff * mapToAdd(i);
    }
}


/* Nouvelle-Fonction */
void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
{
  char *ntn[2] = {"nu","fnu"};
  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
  float xnt[2];
  double nu;
  double numin = rs.minFreq();
  double numax = rs.maxFreq();
  int nmax = 500;
  double dnu = (numax-numin)/nmax;
  for(int k=0; k<nmax; k++) {
    nu = numin+k*dnu;
    xnt[0] = nu;
    xnt[1] = rs.flux(nu);
    nt.Fill(xnt);
  }
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, name);
  return;
}

/* Nouvelle-Fonction */
void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
{
  char *ntn[2] = {"nu","tnu"};
  NTuple nt(2,ntn);  // Creation NTuple (AVEC new )
  float xnt[2];
  double nu;
  double numin = sr.minFreq();
  double numax = sr.maxFreq();
  int nmax = 500;
  double dnu = (numax-numin)/nmax;
  for(int k=0; k<nmax; k++) {
    nu = numin+k*dnu;
    xnt[0] = nu;
    xnt[1] = sr.transmission(nu);
    nt.Fill(xnt);
  }
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, name);
  return;
}