Changeset 610 in Sophya for trunk/SophyaLib


Ignore:
Timestamp:
Nov 22, 1999, 12:25:47 AM (26 years ago)
Author:
ansari
Message:

Corrections diverses - programme tgsky.cc (generation aleatoire de ciel)

et skymixer.cc completee - Reza 21/11/99

Location:
trunk/SophyaLib/SkyT
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/SkyT/Makefile

    r601 r610  
    2424myobjs = $(OBJ)radspec.o $(OBJ)radspecvector.o $(OBJ)nupower.o $(OBJ)blackbody.o $(OBJ)specresp.o $(OBJ)specrespvector.o $(OBJ)squarefilt.o  $(OBJ)trianglefilt.o $(OBJ)gaussfilt.o $(OBJ)convtools.o
    2525
    26 myexe = $(EXE)easyTest $(EXE)skymixer
     26myexe = $(EXE)easyTest $(EXE)skymixer $(EXE)tgsky $(EXE)tgrsr
    2727
    2828all : $(myexe)
    … …  
    3636        $(LINK.cc)  $^ $(LIBS) -o $@
    3737
     38$(EXE)tgsky :  $(OBJ)tgsky.o
     39        $(LINK.cc)  $^ $(LIBS) -o $@
     40
     41$(EXE)tgrsr :  $(OBJ)tgrsr.o $(myobjs)
     42        $(LINK.cc)  $^ $(LIBS) -o $@
    3843
    3944
    4045
     46

  • trunk/SophyaLib/SkyT/gaussfilt.cc

    r607 r610  
    11//--------------------------------------------------------------------------
    22// File and Version Information:
    3 //      $Id: gaussfilt.cc,v 1.2 1999-11-20 21:00:48 ansari Exp $
     3//      $Id: gaussfilt.cc,v 1.3 1999-11-21 23:25:45 ansari Exp $
    44//
    55// Description:
    … …  
    5050  else {
    5151    double tmp = (nu-_nu0)/_s;
    52     return(_a * exp(tmp*tmp));
     52    return(_a * exp(-tmp*tmp));
    5353  }
    5454}
    … …  
    7070{
    7171  os << "GaussianFilter::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
    72   os << " T = A * Exp(((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= "
     72  os << " T = A * Exp(-((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= "
    7373     << _s << "  A= " << _a << endl; 
    7474  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;

  • trunk/SophyaLib/SkyT/gaussfilt.h

    r607 r610  
    22//--------------------------------------------------------------------------
    33// File and Version Information:
    4 //      $Id: gaussfilt.h,v 1.2 1999-11-20 21:00:49 ansari Exp $
     4//      $Id: gaussfilt.h,v 1.3 1999-11-21 23:25:45 ansari Exp $
    55//
    66// Description:
    … …  
    2727//              ---------------------
    2828
    29 // Spectral response in the form  A Exp(((nu-nu0)/s)^2)
     29// Spectral response in the form  A Exp(-((nu-nu0)/s)^2)
    3030
    3131class GaussianFilter:public SpectralResponse

  • trunk/SophyaLib/SkyT/nupower.cc

    r607 r610  
    11//--------------------------------------------------------------------------
    22// File and Version Information:
    3 //      $Id: nupower.cc,v 1.2 1999-11-20 21:00:49 ansari Exp $
     3//      $Id: nupower.cc,v 1.3 1999-11-21 23:25:45 ansari Exp $
    44//
    55// Description:
    … …  
    2424// Constructor --
    2525//----------------
    26 PowerLawSpectra::PowerLawSpectra(double a, double b, double nu0, double numin, double numax)
     26PowerLawSpectra::PowerLawSpectra(double a, double b, double nu0, double dnu, double numin, double numax)
    2727  : RadSpectra(numin, numax)
    2828{
    … …  
    3030  _b = b;
    3131  _nu0 = nu0;
     32  _dnu = (dnu > 1.e-19) ? dnu : 1.;
    3233}
    3334
    … …  
    4041PowerLawSpectra::flux(double nu) const
    4142{
    42   if ((nu < _numin) || (nu > _numax)) return(0.);
    43   else return(  _a * pow((nu-_nu0), _b) );
     43  if ((nu< _numin) || (nu > _numax)) return(0.);
     44  double tmp = (nu-_nu0)/_dnu;
     45  if (tmp <= 0.)  return(0.);
     46  else return(  _a * pow(tmp, _b) );
    4447}
    4548
    … …  
    4750PowerLawSpectra::Print(ostream& os) const
    4851{
    49   os << "PowerLawSpectra::Print f = a (nu-nu0)^b " << _a << "(nu-" << _nu0
    50      << ") ^ " << _b << endl;
     52  os << "PowerLawSpectra::Print f = a ((nu-nu0)/dnu)^b " << _a << "((nu-" << _nu0
     53     << ") / " << _dnu << ") ^ " << _b << endl;
    5154  os << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
    5255  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;

  • trunk/SophyaLib/SkyT/nupower.h

    r607 r610  
    22//--------------------------------------------------------------------------
    33// File and Version Information:
    4 //      $Id: nupower.h,v 1.2 1999-11-20 21:00:49 ansari Exp $
     4//      $Id: nupower.h,v 1.3 1999-11-21 23:25:46 ansari Exp $
    55//
    66// Description:
    … …  
    1717#include "convtools.h"
    1818
    19 // Power law spectra  f = a (nu-nu0)^b   
     19// Power law spectra  f = a ((nu-nu0)/dnu)^b   
    2020class PowerLawSpectra : public RadSpectra
    2121{
    2222public:   //Constructor
    2323 
    24   PowerLawSpectra(double a, double b, double nu0, double numin=0., double numax=9.e49);
     24  PowerLawSpectra(double a, double b, double nu0, double dnu, double numin=0., double numax=9.e49);
    2525
    2626  virtual ~PowerLawSpectra();
    … …  
    3232
    3333protected:
    34   double _a, _b, _nu0;
     34  double _a, _b, _nu0, _dnu;
    3535};
    3636

  • trunk/SophyaLib/SkyT/radspecvector.cc

    r607 r610  
    11//--------------------------------------------------------------------------
    22// File and Version Information:
    3 //      $Id: radspecvector.cc,v 1.2 1999-11-20 21:00:50 ansari Exp $
     3//      $Id: radspecvector.cc,v 1.3 1999-11-21 23:25:46 ansari Exp $
    44//
    55// Description:
    … …  
    5757{
    5858  if ( (nu < _numin) || (nu > _numax) ) return(0.);
    59   double value = -99;
     59  double value = 0.;
    6060  int sizeVecOfNu = _vecOfNu.NElts();
    61   double minVal = _vecOfNu(0);
    62   if(nu <= minVal) return _vecOfFDeNu(0);
     61  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
    6362  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
    6463 
    6564  for (int i=1; i<sizeVecOfNu; i++)
    6665    {
    67       if(nu>= minVal && nu < _vecOfNu(i))
     66      if(nu < _vecOfNu(i))
    6867        {
    6968          double up = _vecOfFDeNu(i) ;
    … …  
    7574          return value;
    7675        }
    77       else
    78         {
    79           minVal = _vecOfNu(i);
    80         }       
    8176    }
    82   return value;
     77  return _vecOfFDeNu(sizeVecOfNu-1);
    8378}
    8479

  • trunk/SophyaLib/SkyT/skm.d

    r607 r610  
    1 #  MAPPATH  PathForFITS files
    2 @MAPPATH /exp/eros/Reza/CartesRT
    3 #  SKYMIX  NbComponents  M_Gorski
     1#  SKYMIX  NbComponents  NSide_HealPix
    42@SKYMIX  2  64
    5 #  FILTER  Nu0 Sigma_Nu Tmax NuMin NuMax
    6 @FILTER 143 12 0.84 50. 300.
    73
    8 #  Sky components
    9 #  MAPFITSFILEi FITSfilename [Normalisation]
    10 @MAPFITSFILE1  map_essai.fits
     4#  READMAP  FitsName
     5#  To Read the already prepared (mixed) Output MAP from FITS file
     6#  FitsName should be the complete path to the FITS file
     7READMAP xxx.fits
     8
     9#  Defining the Detection filter pass - band
     10#  GAUSSFILTER or FILTERFITSFILE  card
     11#  GAUSSFILTER -> Gaussian filter
     12#  FILTERFITSFILE -> Filter (nu, T(nu)) from FITS file
     13#  GAUSSFILTER  Nu0 Sigma_Nu Tmax NuMin NuMax
     14#  FILTERFITSFILE FileName NuMin NuMax
     15#       FileName should be a complete path
     16@GAUSSFILTER 143 12 0.84 50. 300.
     17
     18#  MAPPATH  PathForFITS   
     19#  path for input Sky Component FITS files , and EmissionSpectra files
     20MAPPATH /exp/eros/Reza/CartesRT
     21
     22
     23#  ----        Sky components  ---------
     24#      ----     HealPix maps   ----
     25#  MAPFITSFILEi FITSfilename Normalisation
     26@MAPFITSFILE1  map_essai.fits  1.0
    1127@MAPFITSFILE2    comp2.fits    0.25
    1228
     29#      ---   Emission Spectra definition -----
    1330#   SPECTRAFITSFILEi  FITSfilename  Fmin Fmax
    1431@SPECTRAFITSFILE1  spec1.fits 
    15 @SPECTRAFUNC2  3. 4. 5.
     32
     33#   Other possible definition of emission spectra
     34#   BLACKBODYi  temperature
     35BLACKBODY1 2.726
     36#   For power-law spectra  f = a ((nu-nu0)/dnu)^b)
     37#   POWERLAWSPECTRAi  a nu0 dnu b Fmin Fmax
     38POWERLAWSPECTRA1 1. 150. 50. -0.5 100. 500.
    1639
    1740#  Define the Debug level
    18 @DEBUGLEVEL  5
     41@DEBUGLEVEL  0
     42#  Define the Print level
     43@PRINTLEVEL  0

  • trunk/SophyaLib/SkyT/skymixer.cc

    r607 r610  
    2626#include "gaussfilt.h"
    2727
    28 // ------ Function declaration
    29 void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap,
    30                   PixelMap<float>& mapToAdd, RadSpectra& rs, double K=1.);
    31 void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap,
    32                   PixelMap<double>& mapToAdd, RadSpectra& rs, double K=1.);
    33 
     28// -----------------------------------------------------------------
     29// ------------- Function declaration ------------------------------
    3430int CheckCards(DataCards & dc, string & msg);
    3531char * BuildFITSFileName(string const & fname);
     32SpectralResponse * getSpectralResponse(DataCards & dc);
     33RadSpectra * getEmissionSpectra(DataCards & dc, int nc);
    3634void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
    3735void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
    3836
    39 //  ----- Variable globale ------------
     37template <class T>
     38void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap,
     39                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K=1.);
     40template <class T>
     41void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
     42// -----------------------------------------------------------------
     43
     44//  ----- Global (static) variables ------------
     45static bool rdmap = false;    // true -> Read map first
    4046static char mapPath[256];     // Path for input maps
    4147static int  hp_nside = 32;    // HealPix NSide
    4248static int  nskycomp = 0;     // Number of sky components
    43 static int  debuglev = 0;     // debuglevel
     49static int  debuglev = 0;     // Debug Level
     50static int  printlev = 0;     // Print Level
     51static POutPersist * so = NULL;  // Debug PPFOut file
    4452
    4553// -------------------------------------------------------------------------
    … …  
    5765string msg;
    5866int rc = 0;
    59 POutPersist * so = NULL;
     67RadSpectra * es = NULL;
     68SpectralResponse * sr = NULL;
     69double moy, sig;
     70
     71DataCards dc;
     72so = NULL;
    6073
    6174try {
    6275  string dcard = arg[1];
    6376  cout << " Decoding parameters from file " << dcard << endl;
    64   DataCards dc(dcard);
     77  dc.ReadFile(dcard);
    6578
    6679  rc = CheckCards(dc, msg);
    67   if (rc) goto Fin;
    68 
    69   cout << " skymix/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
    70   cout << "  ... MapPath = " << (string)mapPath << "  DbgLev= " << debuglev << endl;
     80  if (rc) {
     81    cerr << " Error condition -> Rc= " << rc << endl;
     82    cerr << " Msg= " << msg << endl;
     83    return(rc);
     84  }
     85}
     86catch (PException exc) {
     87  msg = exc.Msg();
     88  cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
     89  return(90);
     90}   
     91
     92
     93cout << " skymix/Info : NComp = " <<  nskycomp << " SphereGorski_NSide= " << hp_nside << endl;
     94cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev
     95     << "  PrintLev= " << printlev << endl;
    7196
    7297// We create an output persist file for writing debug objects
    73   if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
    74 
    75   SphereGorski<float> outgs(hp_nside);
    76   cout << " Output Gorski Map  created - NbPixels= " << outgs.NbPixels() << endl;
    77   outgs.SetPixels(0.);
     98if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
     99
     100SphereGorski<float> outgs(hp_nside);
     101bool okout = false;
     102
     103try {
     104  if (rdmap) {  // Reading map from FITS file
     105    FitsIoServer fios;
     106    char ifnm[256];
     107    strncpy(ifnm, dc.SParam("READMAP", 1).c_str(), 255);
     108    ifnm[255] = '\0';
     109    cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
     110    fios.load(outgs, ifnm, 2);
     111    cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
     112  }
     113  else {
     114    cout << " Output HealPix Map  created - NbPixels= " << outgs.NbPixels() << endl;
     115    outgs.SetPixels(0.);
     116  }
    78117
    79118  // Decoding detection pass-band filter
    80   double nu0 = dc.DParam("FILTER", 0, 10.);
    81   double s = dc.DParam("FILTER", 1, 1.);
    82   double a = dc.DParam("FILTER", 2, 1.);
    83   double numin = dc.DParam("FILTER", 3, 0.1);
    84   double numax = dc.DParam("FILTER", 4, 9999);
    85   GaussianFilter filt(nu0, s, a, numin, numax);
    86   cout << " Filter decoded - Created " << endl;
    87   cout << filt << endl;
    88 
    89 // FOR debug
    90   if (debuglev > 0)  SpectralResponse2Nt(filt, *so, "filter");
    91 
     119  sr = getSpectralResponse(dc);
    92120  PrtTim(" After FilterCreation ");
    93121
    94   SphereGorski<float> * ings = NULL;  // Our input map 
    95122  FitsIoServer fios; // Our FITS IO Server
    96   char * flnm, buff[64];
     123  char * flnm, buff[90];
    97124  string key;
    98125
     126  SphereGorski<float> ings(hp_nside);  // The input map
     127  double K = 1.;
    99128  // Loop over sky component
    100129  int sk;
    101130  for(sk = 0; sk<nskycomp; sk++) {
    102131    cout << " Processing sky component No " << sk+1 << endl;
    103     if (ings) { delete ings;  ings = NULL; }
    104     ings = new SphereGorski<float>(hp_nside);
    105     sprintf(buff, "MAPFITSFILE%d", sk+1);
    106     key = buff;
     132    sprintf(buff, "%d", sk+1);
     133    key = (string)"MAPFITSFILE" + buff;
    107134    flnm = BuildFITSFileName(dc.SParam(key, 0));
    108135    cout << " Reading Input FITS map " << (string)flnm << endl;
    109     fios.load(*ings, flnm, 1);
    110 
     136    fios.load(ings, flnm, 2);
     137    K = dc.DParam(key, 1, 1.);
    111138    if (debuglev > 4) {  // Writing tne input map to the outppf
    112139      FIO_SphereGorski<float> fiog(ings);
    113140      fiog.Write(*so, key);
    114       } 
    115     }
    116   }   // End of try block
     141      }
     142    // getting Emission spectra 
     143    if (es) { delete es; es = NULL; }
     144    es = getEmissionSpectra(dc, sk);
     145   
     146    if (printlev > 2) {
     147      MeanSig(ings.DataBlock(), moy, sig );
     148      cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
     149    }
     150    addComponent(*sr, outgs, ings, *es, K);
     151    okout = true;
     152    if (printlev > 1) {
     153      MeanSig(outgs.DataBlock(), moy, sig );
     154      cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
     155    }
     156    sprintf(buff, "End of Proc. Comp. %d ", sk+1);
     157    cout << " --------------------------------------- \n" << endl;   
     158    PrtTim(buff);
     159  }   // End of sky component loop
     160}   // End of try block
    117161
    118162catch (PException exc) {
    119163  msg = exc.Msg();
    120164  cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
    121   rc = 90;
     165  rc = 50;
    122166  }   
    123167
    124 
    125 Fin:
    126 if (so) delete so;  //  Closing the debug ppf file
    127 if (rc == 0)  return(0);
    128 cerr << " Error condition -> Rc= " << rc << endl;
    129 cerr << " Msg= " << msg << endl;
    130 return(rc);
     168// Saving the output map in FITS format
     169if (okout) {
     170  FitsIoServer fios;
     171  fios.save(outgs, arg[2]);
     172  cout << "Output Map (SphereGorski<float>) written to FITS file "
     173       << (string)(arg[2]) << endl;
     174  PrtTim("End of WriteFITS ");
     175// Saving the output map in PPF format
     176  if (narg > 3) {
     177   POutPersist s(arg[3]);
     178   FIO_SphereGorski<float> fiog(&outgs) ;
     179   fiog.Write(s);
     180   cout << "Output Map (SphereGorski<float>) written to POutPersist file " 
     181        << (string)(arg[3]) << endl;
     182   PrtTim("End of WritePPF ");
     183   }
     184  }
     185  if (so) delete so;  //  Closing the debug ppf file
     186  return(rc);
    131187}
    132188
    … …  
    135191//   Function to check datacards
    136192{
     193rdmap = false;
    137194mapPath[0] = '\0';
    138195hp_nside = 32;
    139196nskycomp = 0;
    140197debuglev  = 0;
     198printlev = 0;
    141199
    142200int rc = 0;
     201string key, key2;
     202
    143203  //  Cheking datacards
    144   if ( (!dc.HasKey("SKYMIX")) || (!dc.HasKey("FILTER")) ) {
     204  if (dc.NbParam("SKYMIX") < 2)  {
    145205     rc = 71;
    146      msg = "Invalid parameters  - NO @SKYMIX or @FILTER card ";
     206     msg = "Invalid parameters  - Check @SKYMIX card ";
    147207     return(rc);
    148208  }
     209  key = "READMAP";
     210  if (dc.HasKey(key)) {
     211    if (dc.NbParam(key) < 1) {
     212     rc = 72;
     213     msg = "Invalid parameters  - Check @READMAP card ";
     214     return(rc);
     215    }
     216   else rdmap = true;
     217  }
     218
     219//  Checking detection filter specification
     220  key = "GAUSSFILTER";
     221  key2 = "FILTERFITSFILE";
     222  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) {
     223    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
     224    rc = 73;  return(rc);
     225    }
    149226
    150227  // Decoding number of component and pixelisation parameter
    … …  
    155232  if (ncomp < 1)  {
    156233    msg = "Invalid parameters  - Check datacards @SKYMIX ";
    157     rc = 72;
     234    rc = 74;
    158235    return(rc);
    159236  }
    160237
     238  // Checking detection filter specification
    161239//  Checking input FITS file specifications
    162240  int kc;
    163   string key, key2;
    164241  char buff[256];
    165242  bool pb = false;
     243  string key3;
    166244  for(kc=0; kc<ncomp; kc++) {
    167245    sprintf(buff, "MAPFITSFILE%d", kc+1);
    … …  
    173251    sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
    174252    key = buff;
    175     sprintf(buff, "SPECTRAFUNC%d", kc+1);
     253    sprintf(buff, "BLACKBODY%d", kc+1);
    176254    key2 = buff;
    177     if ( (dc.NbParam(key) < 1) && (dc.NbParam(key2) < 2) )  {   
    178       msg = "Missing card or invalid parameters : " + key + " or " + key2;
     255    sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
     256    key3 = buff;
     257    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6))  {   
     258      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3;
    179259      pb = true;  break;
    180260    }
    … …  
    183263
    184264  if (pb)  {
    185     rc = 72;
    186     return(72);
    187   }
    188 
    189 //  Checking detection filter specification
    190   key = "FILTER";
    191   if (dc.NbParam(key) < 3) {
    192     msg = "Missing card or invalid parameters : " + key;
    193     rc = 73;  return(rc);
    194     }
     265    rc = 75;
     266    return(75);
     267  }
     268
    195269
    196270// Initialiazing parameters
    … …  
    206280  key = "DEBUGLEVEL";
    207281  debuglev =  dc.IParam(key, 0, 0);
     282  key = "PRINTLEVEL";
     283  printlev =  dc.IParam(key, 0, 0);
    208284  return(rc);
    209285}
    … …  
    219295
    220296/* Nouvelle-Fonction */
    221 // template <class T>
    222 void addComponent(SpectralResponse&  sr, PixelMap<float>& finalMap,
    223                   PixelMap<float>& mapToAdd, RadSpectra& rs, double K)
     297SpectralResponse * getSpectralResponse(DataCards & dc)
     298{
     299  SpectralResponse * filt = NULL;
     300
     301  string key = "FILTERFITSFILE";
     302  string key2 = "GAUSSFILTER";
     303  string ppfname = "filter";
     304
     305  if (dc.HasKey(key) ) {      // Reading FITS filter file
     306    FitsIoServer fios;
     307    char ifnm[256];
     308    strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
     309    ifnm[255] = '\0';
     310    Matrix mtx(2,10);
     311    fios.load(mtx, ifnm);
     312    double numin = dc.DParam(key, 2, 1.);
     313    double numax = dc.DParam(key, 3, 9999.);
     314    Vector nu(mtx.NCols());
     315    Vector fnu(mtx.NCols());
     316    for(int k=0; k<mtx.NCols(); k++) {
     317      nu(k) = mtx(0, k);
     318      fnu(k) = mtx(1, k);
     319    }
     320  filt = new SpecRespVec(nu, fnu, numin, numax);
     321  ppfname = key;
     322  }
     323  else if (dc.HasKey(key2) ) {   // creating GaussianFilter
     324    double nu0 = dc.DParam(key2, 0, 10.);
     325    double s = dc.DParam(key2, 1, 1.);
     326    double a = dc.DParam(key2, 2, 1.);
     327    double numin = dc.DParam(key2, 3, 0.1);
     328    double numax = dc.DParam(key2, 4, 9999);
     329    filt = new GaussianFilter(nu0, s, a, numin, numax);
     330    ppfname = key2;
     331    }
     332  if (filt == NULL) throw ParmError("datacard error ! No detection filter");
     333  cout << " Filter decoded - Created " << endl;
     334  cout << *filt << endl;
     335
     336// for debug
     337  if (debuglev > 1)  SpectralResponse2Nt(*filt, *so, ppfname);
     338  return(filt);
     339}
     340
     341/* Nouvelle-Fonction */
     342RadSpectra * getEmissionSpectra(DataCards & dc, int nc)
     343{
     344  char numb[16];
     345  sprintf(numb, "%d", nc+1);
     346
     347  string key = (string)"SPECTRAFITSFILE" + numb;
     348  string key2 = (string)"BLACKBODY" + numb;
     349  string key3 = (string)"POWERLAWSPECTRA" + numb;
     350  string ppfname = "espectra";
     351
     352  RadSpectra * rs = NULL;
     353  if (dc.HasKey(key) ) {    // Reading emission spectra from file
     354    char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
     355    cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
     356    FitsIoServer fios;
     357    Matrix mtx(2,10);
     358    fios.load(mtx, ifnm);
     359    double numin = dc.DParam(key, 2, 1.);
     360    double numax = dc.DParam(key, 3, 9999.);
     361    Vector nu(mtx.NCols());
     362    Vector tnu(mtx.NCols());
     363    for(int k=0; k<mtx.NCols(); k++) {
     364      nu(k) = mtx(0, k);
     365      tnu(k) = mtx(1, k);
     366    }
     367  rs = new RadSpectraVec(nu, tnu, numin, numax);
     368  ppfname = key;
     369  }
     370  else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
     371    rs = new BlackBody(dc.DParam(key2, 0, 2.726));
     372    ppfname = key2;
     373  }
     374  else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
     375    double a = dc.DParam(key3, 0, 1.);
     376    double nu0 = dc.DParam(key3, 1, 100.);
     377    double dnu = dc.DParam(key3, 2, 10.);
     378    double b = dc.DParam(key3, 3, 0.);
     379    double numin = dc.DParam(key3, 4, 0.1);
     380    double numax = dc.DParam(key3, 5, 9999);
     381    rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
     382    ppfname = key3;
     383  }
     384  if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
     385  cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
     386  cout << *rs << endl;
     387// for debug
     388  if (debuglev > 2)  RadSpec2Nt(*rs, *so, ppfname);
     389  return(rs);
     390}
     391
     392
     393
     394
     395/* Nouvelle-Fonction */
     396template <class T>
     397void addComponent(SpectralResponse&  sr, PixelMap<T>& finalMap,
     398                  PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
    224399{
    225400  // finalMap = finalMap + coeff* mapToAdd
    … …  
    229404    throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
    230405  double coeff = rs.filteredIntegratedFlux(sr) * K;
     406  if (printlev > 1)
     407    cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
    231408  for(int i=0; i<finalMap.NbPixels(); i++)
    232409    {
    … …  
    236413
    237414/* Nouvelle-Fonction */
    238 void addComponent(SpectralResponse&  sr, PixelMap<double>& finalMap,
    239                   PixelMap<double>& mapToAdd, RadSpectra& rs, double K)
    240 {
    241   // finalMap = finalMap + coeff* mapToAdd
    242   // coeff    =  convolution of sr and rs
    243   // compute the coefficient corresponding to mapToAdd
    244   if (finalMap.NbPixels() != mapToAdd.NbPixels())   
    245     throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
    246   double coeff = rs.filteredIntegratedFlux(sr) * K;
    247   for(int i=0; i<finalMap.NbPixels(); i++)
    248     {
    249       finalMap(i) += coeff * mapToAdd(i);
    250     }
    251 }
    252 
     415template <class T>
     416void  MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig)
     417
     418{
     419  gmoy=0.;
     420  gsig = 0.;
     421  double valok;
     422  for(int k=0; k<(int)dbl.Size(); k++) {
     423    valok = dbl(k);
     424    gmoy += valok;  gsig += valok*valok;
     425  }
     426  gmoy /= (double)dbl.Size();
     427  gsig = gsig/(double)dbl.Size() - gmoy*gmoy;
     428  if (gsig >= 0.) gsig = sqrt(gsig);
     429}
    253430
    254431/* Nouvelle-Fonction */

  • trunk/SophyaLib/SkyT/specrespvector.cc

    r607 r610  
    11//--------------------------------------------------------------------------
    22// File and Version Information:
    3 //      $Id: specrespvector.cc,v 1.2 1999-11-20 21:00:53 ansari Exp $
     3//      $Id: specrespvector.cc,v 1.3 1999-11-21 23:25:47 ansari Exp $
    44//
    55// Description:
    … …  
    5757{
    5858  if ( (nu < _numin) || (nu > _numax) ) return(0.);
    59   double value = -99;
     59  double value = 0.;
    6060  int sizeVecOfNu = _vecOfNu.NElts();
    61   double minVal = _vecOfNu(0);
    62   if(nu <= minVal) return _vecOfFDeNu(0);
     61  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
    6362  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
    6463
    6564  for (int i=1; i<sizeVecOfNu; i++)
    6665    {
    67       if(nu>= minVal && nu < _vecOfNu(i))
     66      if(nu < _vecOfNu(i))
    6867        {
    6968          double up = _vecOfFDeNu(i) ;
    … …  
    7574          return value;
    7675        }
    77       else
    78         {
    79           minVal = _vecOfNu(i);
    80         }       
    8176    }
    8277  return value;
Note: See TracChangeset for help on using the changeset viewer.