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

#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.);

// -------------------------------------------------------------------------
//                             main program 
// -------------------------------------------------------------------------
int main(int narg, char * arg[])
{
  if ((narg < 2) || (strcmp(arg[1], "-h") == 0) ) {
    cout << " skymixer / Error args \n Usage: skymixer parameterFile " << endl;
    exit(0);
  }

InitTim();

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

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

  // 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;
    goto problem;
  }

  int kc;
  string key;
  char buff[256];
  bool pb = false;
  for(kc=0; kc<ncomp; kc++) {
    sprintf(buff, "SCMFILE%d", kc+1);
    key = buff;
    if (dc.NbParam(key) < 1) {   
      msg = "Missing or invalid card : " + key;
      pb = true;  break;
    }
    sprintf(buff, "SCSPEC%d", kc+1);
    key = buff;
    if (dc.NbParam(key) < 1)  {   
      msg = "Missing or invalid card : " + key;
      pb = true;  break;
    }

    }

  if (pb)  {
    rc = 72;
    goto problem;
  }

  cout << " skymix/Info : NComp = " <<  ncomp << " M_SphereGorski= " << mg << endl;

  SphereGorski<float> outgs(mg);
  cout << " Creating Output Gorski Map  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;

  PrtTim(" After FilterCreation ");
  }

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


problem:
if (rc == 0)  return(0);
cerr << " Error condition -> Rc= " << rc << endl;
cerr << " Msg= " << msg << endl;
return(rc);
}

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

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