source: Sophya/trunk/SophyaProg/PMixer/skymixer.cc@ 1295

Last change on this file since 1295 was 1295, checked in by ansari, 25 years ago

Transport de mot-cle FITS depuis HDU 2 (et pas 1) ds skymixer , Reza 6/11/2000
File size: 27.0 KB
RevLine 
[931]1#include "pmixer.h"
[761]2
[928]3/*!
[931]4 * \defgroup PMixer PMixer module
5 * This module contains programs which:
6 * <UL>
7 * <LI> add several sky components, taking into account their
8 * radiation spectra and convoluting them with a given filter
9 * response : skymixer
10 * <LI> create a map with point source : extractRS
11 * <LI> generate sky components, radiation spectra and spectral
12 * response (small generator of maps) : tgsky and tgrsr
13 * </UL>
14 * A detailed description may be found at:
15 */
[928]16/*!
17 * \ingroup PMixer
18 * \file skymixer.cc
19 *\brief \b PROGRAM \b skyMixer <BR>
20 * add several sky components, taking into account their
21 *radiation spectra and convoluting them with a given filter
22 *response
[931]23 */
[928]24
[761]25// -----------------------------------------------------------------
26// ------------- Function declaration ------------------------------
27int CheckCards(DataCards & dc, string & msg);
28char * BuildFITSFileName(string const & fname);
29SpectralResponse * getSpectralResponse(DataCards & dc);
30RadSpectra * getEmissionSpectra(DataCards & dc, int nc);
[1252]31void SKM_MergeFITSKeywords(char * flnm);
[1295]32void SKM_MergeFITSKeywords2(DVList & dvl);
[761]33void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
34void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
[1252]35
[876]36// to add different sky components and corresponding tools
37//----------------------------------------------------------
38template <class T>
39void addComponent(SpectralResponse& sr,
40 PixelMap<T>& finalMap,
41 PixelMap<T>& mapToAdd,
42 RadSpectra& rs, double K=1.);
43//
[761]44template <class T>
[880]45void addComponentBeta(SphereHEALPix<T>& finalMap,
46 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
47 SphereHEALPix<T>& betaMap, double normFreq, double K);
[876]48//
[761]49template <class T>
[876]50void integratedMap(SpectralResponse& sr,
[880]51 SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
[904]52
53//
[876]54template <class T>
[880]55void addComponentBeta(SphereHEALPix<T>& finalMap,
56 SphereHEALPix<T>& mapToAdd,
57 SphereHEALPix<T>& intBetaMap, double K);
[876]58//
[904]59template <class T>
60void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
61 double theta,double phi,double amp,double temp);
62//
[761]63// -----------------------------------------------------------------
64
65// ----- Global (static) variables ------------
66static bool rdmap = false; // true -> Read map first
67static char mapPath[256]; // Path for input maps
68static int hp_nside = 32; // HealPix NSide
69static int nskycomp = 0; // Number of sky components
70static int debuglev = 0; // Debug Level
71static int printlev = 0; // Print Level
[1252]72static POutPersist * so = NULL; // Debug PPFOut file
73static DVList * dvl_fitskw = NULL; // Global DVList for all FITS Keywords
[761]74
[1252]75// --------- SkyMixer Version --------------
76static double skm_version = 1.4;
77
[761]78// -------------------------------------------------------------------------
79// main program
80// -------------------------------------------------------------------------
81int main(int narg, char * arg[])
82{
83 if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
84 cout << " Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
85 exit(0);
86 }
[886]87
[904]88 InitTim();
[1282]89
90 cout << " skymixer Version= " << skm_version << " SophyaVersion= "
91 << SophyaVersion() << endl;
[904]92
93 string msg;
94 int rc = 0;
95 RadSpectra * es = NULL;
96 SpectralResponse * sr = NULL;
97 double moy, sig;
98
99 DataCards dc;
100 so = NULL;
[1280]101 // DVList for merging all FITS keywords and datacards
[1252]102 dvl_fitskw = new DVList;
[904]103
104 try {
105 string dcard = arg[1];
106 if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
107 dc.ReadFile(dcard);
108
109 rc = CheckCards(dc, msg);
110 if (rc) {
111 cerr << " Error condition -> Rc= " << rc << endl;
112 cerr << " Msg= " << msg << endl;
113 return(rc);
114 }
115 }
116 catch (PException exc) {
117 msg = exc.Msg();
118 cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
119 return(90);
120 }
121
122
123 cout << " skymix/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
124 cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
125 << " PrintLev= " << printlev << endl;
126
127 // We create an output persist file for writing debug objects
128 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
129
130 SphereHEALPix<float> outgs(hp_nside);
131 try{
132 if (rdmap) { // Reading map from FITS file
133 char ifnm[256];
134 strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
135 ifnm[255] = '\0';
136 cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
137 {
[937]138 FITS_SphereHEALPix<float> fios(&outgs);
[1062]139 fios.Read(ifnm,2);
[1252]140 // Getting FITS keywords in primary header
[1295]141 // SKM_MergeFITSKeywords(ifnm);
142 // Getting FITS keywords in data object header
143 SKM_MergeFITSKeywords2(outgs.Info());
[904]144 }
145 if(printlev>0)
146 cout << " Output HealPIx Map read - NbPixels= " <<
147 outgs.NbPixels() << endl;
[886]148 if (printlev > 0) {
149 MeanSig(outgs.DataBlock(), moy, sig );
[899]150 cout << " MeanSig for outpout map - Mean= " <<
151 moy << " Sigma= " << sig << endl;
[886]152 }
[761]153 }
[886]154 else {
[899]155 if(printlev>0)
156 cout << " Output HealPix Map created - NbPixels= " <<
157 outgs.NbPixels() << endl;
[886]158 outgs.SetPixels(0.);
159 }
[876]160
[886]161 // Decoding detection pass-band filter
162 sr = getSpectralResponse(dc);
163 PrtTim(" After FilterCreation ");
[876]164
[886]165 char * flnm, buff[90];
166 string key;
[876]167
[886]168 double K = 1.;
169 double freqOfMap = 1.;
170 // Loop over sky component
171 int sk;
172 for(sk = 0; sk<nskycomp; sk++) {
[899]173 cout << "------------------------------------" << endl;
[886]174 cout << " Processing sky component No " << sk+1 << endl;
175
[904]176
[886]177 sprintf(buff, "%d", sk+1);
[904]178 key = (string)"DIPOLE" + buff;
179 // check for an eventual dipole
[899]180 if(dc.HasKey(key))
181 {
[886]182 if (es) { delete es; es = NULL; }
[904]183 double temp = -10.;
184 double theta= dc.DParam(key,1,1.);
185 double phi = dc.DParam(key,2,1.);
186 double amp = dc.DParam(key,3,1.);
187 if(dc.NbParam(key)>3)
188 {
189 temp = dc.DParam(key,4,1.);
190 }
191 cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
192 addDipole(*sr, outgs,theta,phi,amp,temp);
[886]193 }
194 else
195 {
[904]196 sprintf(buff, "%d", sk+1);
197 key = (string)"MAPFITSFILE" + buff;
[886]198 flnm = BuildFITSFileName(dc.SParam(key, 0));
[904]199
200 K = dc.DParam(key, 1, 1.);
201
202 cout << " Reading Input FITS map " << (string)flnm << endl;
203 SphereHEALPix<float> ings(hp_nside);
[899]204 {
[937]205 FITS_SphereHEALPix<float> fiosIn(&ings);
[1062]206 fiosIn.Read(flnm,2);
[1252]207 // Getting FITS keywords in primary header
[1295]208 // SKM_MergeFITSKeywords(flnm);
209 // Getting FITS keywords in data object header
210 SKM_MergeFITSKeywords2(ings.Info());
[899]211 }
212 if (debuglev > 4) { // Writing the input map to the outppf
[904]213 FIO_SphereHEALPix<float> fiog(ings);
214 fiog.Write(*so, key);
[899]215 }
[904]216 if (printlev > 2) {
217 MeanSig(ings.DataBlock(), moy, sig );
218 cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
219 }
220 bool mapDependentOfFreq = false;
221 key = (string)"BETAFITSFILE"+ buff;
222 if(dc.HasKey(key))
223 {
224 mapDependentOfFreq = true;
225 }
226
227 // getting Emission spectra
228 if(!mapDependentOfFreq)
229 {
230 if (es) { delete es; es = NULL; }
231 es = getEmissionSpectra(dc, sk);
232 addComponent(*sr, outgs, ings, *es, K);
233 }
[886]234 else
235 {
[904]236 key = (string)"BETAFITSFILE"+ buff;
237 //SphereHEALPix<float> betaMap;
238 flnm = BuildFITSFileName(dc.SParam(key, 0));
239 double normFreq = dc.DParam(key, 1, 1.);
240 if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
241 int nSideForInt = dc.DParam(key, 2, 1.);
242 if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
243 cout << "....BetaFits... Reading Beta FITS map " << (string)flnm << endl;
244 SphereHEALPix<float> betaMap(hp_nside);
245 {
[937]246 FITS_SphereHEALPix<float> fiosBM(&betaMap);
[1062]247 fiosBM.Read(flnm,2);
[1252]248 // Getting FITS keywords in primary header
[1295]249 // SKM_MergeFITSKeywords(flnm);
250 // Getting FITS keywords in data object header
251 SKM_MergeFITSKeywords2(betaMap.Info());
[904]252 }
[937]253 if (printlev > 2) {
254 MeanSig(betaMap.DataBlock(), moy, sig );
255 cout << " MeanSig for Beta map - Mean= " << moy << " Sigma= " << sig << endl;
256 }
[899]257 if (debuglev > 4) { // Writing the input map to the outppf
[904]258 FIO_SphereHEALPix<float> fiogs(betaMap);
259 fiogs.Write(*so, key);
[899]260 }
[904]261
262 if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
263 bool bydefault = true;
264 if(!bydefault)
265 addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
266 else
267 {
268 // integrate the betamap over the SpectralResponse
269 SphereHEALPix<float> intBetaMap(nSideForInt);
270 integratedMap(*sr, betaMap, normFreq, intBetaMap);
271 if (debuglev > 4) { // Writing the input map to the outppf
272 FIO_SphereHEALPix<float> fiogs2(intBetaMap);
273 fiogs2.Write(*so, "INTBETAMAP");
274 }
275
276 betaMap.Resize(8);
[937]277 if (printlev > 4)
278 {
279 MeanSig(intBetaMap.DataBlock(), moy, sig );
280 cout << "....BetaFits... MeanSig for intBetaMap - Mean= "
281 << moy << " Sigma= " << sig << endl;
282 }
[904]283 // add the integrated beta map
284 addComponentBeta(outgs,ings,intBetaMap, K);
285 }
[886]286 }
[904]287
288 MeanSig(outgs.DataBlock(), moy, sig );
289 cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
290 cout << "-------------------------------------------------" << endl;
291
292 sprintf(buff, "End of Processing Component %d ", sk+1);
293 PrtTim(buff);
294 }
295 }
296 }
297 catch(PException exc)
298 {
299 cout << "catched PException" << endl;
300 msg = exc.Msg();
301 cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
302 rc = 50;
303 return(50);
[886]304 }
[904]305
[986]306 try {
[886]307 // Saving the output map in FITS format
[986]308 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
309 << (string)(arg[2]) << endl;
310 {
[1239]311 FitsOutFile fios(arg[2]);
312 fios.firstImageOnPrimaryHeader(false); // Use secondary header
[1252]313 DVList& dvl = (*dvl_fitskw);
314 dvl["PDMTYPE"] = "COMPMAP";
315 dvl.SetComment("PDMTYPE", "Planck Data Model Type");
316 dvl["SOPHYVER"] = SophyaVersion();
317 dvl.SetComment("SOPHYVER", "Sophya Version number");
318 dvl["SKYMVER"] = skm_version;
319 dvl.SetComment("SKYMVER", "skymixer Version number");
[1295]320 // fios.DVListIntoPrimaryHeader(dvl);
321 outgs.Info() = dvl;
[1239]322 fios << outgs ;
[986]323 }
324 PrtTim("End of WriteFITS ");
325 // Saving the output map in PPF format
326 if (narg > 3) {
327 POutPersist s(arg[3]);
328 FIO_SphereHEALPix<float> fiog(&outgs) ;
329 fiog.Write(s);
330 cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
331 << (string)(arg[3]) << endl;
332 PrtTim("End of WritePPF ");
333 }
[899]334 }
[986]335 catch(PException exc)
336 {
337 cout << "catched PException (2)" << endl;
338 msg = exc.Msg();
339 cerr << " !!!! skymixer(2) - Catched exception - Msg= " << exc.Msg() << endl;
340 rc = 55;
341 return(55);
342 }
343
[886]344 if (so) delete so; // Closing the debug ppf file
345 return(rc);
[761]346}
347
348/* Nouvelle-Fonction */
349int CheckCards(DataCards & dc, string & msg)
350// Function to check datacards
351{
352rdmap = false;
353mapPath[0] = '\0';
354hp_nside = 32;
355nskycomp = 0;
356debuglev = 0;
357printlev = 0;
358
359int rc = 0;
[904]360string key, key2,key3;
[761]361
[1280]362 DVList & dvl = *dvl_fitskw;
[761]363 // Cheking datacards
364 if (dc.NbParam("SKYMIX") < 2) {
365 rc = 71;
366 msg = "Invalid parameters - Check @SKYMIX card ";
367 return(rc);
368 }
[1280]369 dvl.SetS("SKYMIX", dc.GetParams("SKYMIX"));
[761]370 key = "READMAP";
371 if (dc.HasKey(key)) {
372 if (dc.NbParam(key) < 1) {
373 rc = 72;
374 msg = "Invalid parameters - Check @READMAP card ";
375 return(rc);
376 }
[1280]377 else rdmap = true;
378 dvl.SetS(key, dc.GetParams(key));
[761]379 }
380
381// Checking detection filter specification
382 key = "GAUSSFILTER";
383 key2 = "FILTERFITSFILE";
[904]384 key3 = "DIPOLE";
385 if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) {
386 msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
[761]387 rc = 73; return(rc);
388 }
389
[1280]390 if (dc.HasKey(key)) dvl.SetS("GAUSFILT", dc.GetParams(key));
391 if (dc.HasKey(key2)) dvl.SetS("FILTFILE", dc.GetParams(key2));
392 if (dc.HasKey(key3)) dvl.SetS("DIPOLE", dc.GetParams(key3));
393
[761]394 // Decoding number of component and pixelisation parameter
395 int mg = 32;
396 int ncomp = 0;
397 ncomp = dc.IParam("SKYMIX", 0, 0);
398 mg = dc.IParam("SKYMIX", 1, 32);
399 if (ncomp < 1) {
400 msg = "Invalid parameters - Check datacards @SKYMIX ";
401 rc = 74;
402 return(rc);
403 }
404
405 // Checking detection filter specification
[876]406 // Checking input FITS file specifications
[761]407 int kc;
408 char buff[256];
409 bool pb = false;
[876]410 string key4;
411 string key5;
[904]412 string key6;
[761]413 for(kc=0; kc<ncomp; kc++) {
414 sprintf(buff, "MAPFITSFILE%d", kc+1);
415 key = buff;
[904]416 sprintf(buff, "DIPOLE%d", kc+1);
417 key3 = buff;
418 if (dc.NbParam(key) < 1 && dc.NbParam(key3)<1) {
419 msg = "Missing or invalid card : " + key + " " + key2 + " " + key3;
[761]420 pb = true; break;
421 }
[1280]422
423 sprintf(buff, "MAPFIL%d", kc+1);
424 if (dc.HasKey(key)) dvl.SetS(buff, dc.GetParams(key));
425 sprintf(buff, "DIPOLE%d", kc+1);
426 if (dc.HasKey(key3)) dvl.SetS(buff, dc.GetParams(key3));
427
[761]428 sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
429 key = buff;
430 sprintf(buff, "BLACKBODY%d", kc+1);
431 key2 = buff;
432 sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
433 key3 = buff;
[876]434 sprintf(buff, "BETAFITSFILE%d", kc+1);
435 key4 = buff;
436 sprintf(buff, "DIPOLE%d", kc+1);
437 key5 = buff;
[904]438 sprintf(buff, "DERIVBB%d", kc+1);
439 key6 = buff;
[876]440 if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
[904]441 && (dc.NbParam(key6)<1) && (dc.NbParam(key5)<3)) {
[876]442 msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
[761]443 pb = true; break;
444 }
445
[1280]446 sprintf(buff, "SPECTF%d", kc+1);
447 if (dc.HasKey(key)) dvl.SetS(buff, dc.GetParams(key));
448 sprintf(buff, "BLACKB%d", kc+1);
449 if (dc.HasKey(key2)) dvl.SetS(buff, dc.GetParams(key2));
450 sprintf(buff, "PLAWSP%d", kc+1);
451 if (dc.HasKey(key3)) dvl.SetS(buff, dc.GetParams(key3));
452 sprintf(buff, "BETAFI%d", kc+1);
453 if (dc.HasKey(key4)) dvl.SetS(buff, dc.GetParams(key4));
454 sprintf(buff, "DIPOLE%d", kc+1);
455 if (dc.HasKey(key5)) dvl.SetS(buff, dc.GetParams(key5));
456 sprintf(buff, "DERIBB%d", kc+1);
457 if (dc.HasKey(key6)) dvl.SetS(buff, dc.GetParams(key6));
458
[761]459 }
460
461 if (pb) {
462 rc = 75;
463 return(75);
464 }
465
466
467// Initialiazing parameters
468 rc = 0;
469 msg = "OK";
470 nskycomp = ncomp;
471 hp_nside = mg;
472
[1280]473
[761]474// Checking for PATH definition card
475 key = "MAPPATH";
476 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
477 mapPath[255] = '\0';
478 key = "DEBUGLEVEL";
479 debuglev = dc.IParam(key, 0, 0);
480 key = "PRINTLEVEL";
481 printlev = dc.IParam(key, 0, 0);
482 return(rc);
483}
484
485static char buff_flnm[1024]; // Mal protege !
486/* Nouvelle-Fonction */
487char* BuildFITSFileName(string const & fname)
488{
489if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
490else sprintf(buff_flnm, "%s", fname.c_str());
491return(buff_flnm);
492}
493
494/* Nouvelle-Fonction */
495SpectralResponse * getSpectralResponse(DataCards & dc)
496{
497 SpectralResponse * filt = NULL;
498
499 string key = "FILTERFITSFILE";
500 string key2 = "GAUSSFILTER";
501 string ppfname = "filter";
502
503 if (dc.HasKey(key) ) { // Reading FITS filter file
504 char ifnm[256];
505 strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
506 ifnm[255] = '\0';
[1120]507 Matrix mtx;
[1239]508 FitsInFile fiis(ifnm);
509 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
510 fiis >> mtx ;
[1252]511 // Getting FITS keywords in primary header
[1295]512 // SKM_MergeFITSKeywords(ifnm);
513 // Getting FITS keywords in data object header
514 SKM_MergeFITSKeywords2(mtx.Info());
[761]515 double numin = dc.DParam(key, 2, 1.);
516 double numax = dc.DParam(key, 3, 9999.);
517 Vector nu(mtx.NCols());
518 Vector fnu(mtx.NCols());
519 for(int k=0; k<mtx.NCols(); k++) {
520 nu(k) = mtx(0, k);
521 fnu(k) = mtx(1, k);
522 }
523 filt = new SpecRespVec(nu, fnu, numin, numax);
524 ppfname = key;
525 }
526 else if (dc.HasKey(key2) ) { // creating GaussianFilter
527 double nu0 = dc.DParam(key2, 0, 10.);
528 double s = dc.DParam(key2, 1, 1.);
529 double a = dc.DParam(key2, 2, 1.);
530 double numin = dc.DParam(key2, 3, 0.1);
531 double numax = dc.DParam(key2, 4, 9999);
532 filt = new GaussianFilter(nu0, s, a, numin, numax);
533 ppfname = key2;
534 }
535 if (filt == NULL) throw ParmError("datacard error ! No detection filter");
[899]536 if(printlev>0)
537 {
538 cout << endl;
539 cout << " Filter decoded - Created " << endl;
540 cout << *filt << endl;
541 }
542 // for debug
[761]543 if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
544 return(filt);
545}
546
547/* Nouvelle-Fonction */
548RadSpectra * getEmissionSpectra(DataCards & dc, int nc)
549{
550 char numb[16];
551 sprintf(numb, "%d", nc+1);
552
553 string key = (string)"SPECTRAFITSFILE" + numb;
554 string key2 = (string)"BLACKBODY" + numb;
[904]555 string key5 = (string)"DERIVBB" + numb;
[761]556 string key3 = (string)"POWERLAWSPECTRA" + numb;
557 string ppfname = "espectra";
558
559 RadSpectra * rs = NULL;
560 if (dc.HasKey(key) ) { // Reading emission spectra from file
561 char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
562 cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
[1120]563 Matrix mtx;
[1239]564 FitsInFile fiis(ifnm);
565 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
566 fiis >> mtx ;
[1252]567 // Getting FITS keywords in primary header
[1295]568 // SKM_MergeFITSKeywords(ifnm);
569 // Getting FITS keywords in data object header
570 SKM_MergeFITSKeywords2(mtx.Info());
[761]571 double numin = dc.DParam(key, 2, 1.);
572 double numax = dc.DParam(key, 3, 9999.);
573 Vector nu(mtx.NCols());
574 Vector tnu(mtx.NCols());
575 for(int k=0; k<mtx.NCols(); k++) {
576 nu(k) = mtx(0, k);
577 tnu(k) = mtx(1, k);
578 }
579 rs = new RadSpectraVec(nu, tnu, numin, numax);
580 ppfname = key;
581 }
582 else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
583 rs = new BlackBody(dc.DParam(key2, 0, 2.726));
584 ppfname = key2;
585 }
[876]586 else if (dc.HasKey(key5) ) { // Creating Dipole
[904]587 rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
[876]588 ppfname = key5;
589 }
[761]590 else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
591 double a = dc.DParam(key3, 0, 1.);
592 double nu0 = dc.DParam(key3, 1, 100.);
593 double dnu = dc.DParam(key3, 2, 10.);
594 double b = dc.DParam(key3, 3, 0.);
595 double numin = dc.DParam(key3, 4, 0.1);
596 double numax = dc.DParam(key3, 5, 9999);
597 rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
598 ppfname = key3;
599 }
600 if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
601 cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
602 cout << *rs << endl;
603// for debug
604 if (debuglev > 2) RadSpec2Nt(*rs, *so, ppfname);
605 return(rs);
606}
607
608
609
610
[904]611template <class T>
612void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
613 double theta,double phi,double amp,double temp)
614{
615 DerivBlackBody dbb;
616 if(temp>0) dbb.setTemperature(temp);
617 double coeff = dbb.filteredIntegratedFlux(sr) * amp;
618 UnitVector vd(theta,phi);
619 UnitVector vc(theta,phi);
620
621 for(int i=0; i<finalMap.NbPixels(); i++)
622 {
623 double thetar,phir;
624 finalMap.PixThetaPhi(i,thetar,phir);
625 vc.SetThetaPhi(thetar, phir);
626 finalMap(i) += vd.Psc(vc)*coeff;
627 }
628 if (debuglev > 4) { // Writing the input map to the outppf
629 SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
630 for(int i=0; i<finalMap.NbPixels(); i++)
631 {
632 double thetar,phir;
633 finalMap.PixThetaPhi(i,thetar,phir);
634 vc.SetThetaPhi(thetar, phir);
635 ings(i) = vd.Psc(vc);
636 }
637 FIO_SphereHEALPix<float> fiog(ings);
638 fiog.Write(*so, "dipole");
639 cout << "Debug the dipole map....saved in debug file !" << endl;
640 }
641}
[761]642/* Nouvelle-Fonction */
643template <class T>
644void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
645 PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
646{
647 // finalMap = finalMap + coeff* mapToAdd
648 // coeff = convolution of sr and rs
649 // compute the coefficient corresponding to mapToAdd
650 if (finalMap.NbPixels() != mapToAdd.NbPixels())
651 throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
652 double coeff = rs.filteredIntegratedFlux(sr) * K;
653 if (printlev > 1)
654 cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
655 for(int i=0; i<finalMap.NbPixels(); i++)
656 {
657 finalMap(i) += coeff * mapToAdd(i);
658 }
659}
[876]660/* Nouvelle-Fonction */
661template <class T>
[880]662void addComponentBeta(SphereHEALPix<T>& finalMap,
663 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
664 SphereHEALPix<T>& betaMap, double normFreq, double K)
[876]665{
666 // finalMap = finalMap + coeff* mapToAdd
667 // coeff = convolution of sr and rs
668 // compute the coefficient corresponding to mapToAdd
669 // betaMap is the map of (beta(theta,phi))
[761]670
[876]671 int nbpix = finalMap.NbPixels();
672 if (nbpix != mapToAdd.NbPixels())
673 throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
674 if (printlev > 1)
675 {
676 cout << "addComponentBeta - Coeff= " << K << endl;
677 cout << "nb pixels: " << finalMap.NbPixels() << endl;
678 }
[885]679 SphereHEALPix<T> bigBetaMap(sqrt((double)nbpix/12));
[876]680 if(nbpix != betaMap.NbPixels())
681 {
682 Sph2Sph(betaMap,bigBetaMap);
683 }
684 for(int i=0; i<finalMap.NbPixels(); i++)
685 {
686 // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
687 RadSpectra* rs = new PowerLawSpectra
688 (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
689 double coeff = rs->filteredIntegratedFlux(sr);
690 finalMap(i) += coeff*K*mapToAdd(i);
691 }
692}
693
694template <class T>
695void integratedMap(SpectralResponse& sr,
[880]696 SphereHEALPix<T>& betaMap,
[876]697 double normFreq,
[880]698 SphereHEALPix<T>& intBetaMap)
[876]699{
[899]700 PowerLawSpectra rs(1.,-2., 0., normFreq);
701
702 if(betaMap.NbPixels()!=intBetaMap.NbPixels())
[876]703 {
[899]704 Sph2Sph(betaMap,intBetaMap);
705 for(int i=0; i<intBetaMap.NbPixels(); i++)
706 {
707 rs.setExp(-intBetaMap(i));
708 double coeff = rs.filteredIntegratedFlux(sr);
709 intBetaMap(i) = coeff;
710 }
[876]711 }
[899]712 else
713 {
714 for(int i=0; i<intBetaMap.NbPixels(); i++)
715 {
716 rs.setExp(-betaMap(i));
717 double coeff = rs.filteredIntegratedFlux(sr);
718 intBetaMap(i) = coeff;
719 }
720 }
[876]721}
722
723template <class T>
[880]724void addComponentBeta(SphereHEALPix<T>& finalMap,
725 SphereHEALPix<T>& mapToAdd,SphereHEALPix<T>& intBetaMap, double K)
[876]726{
727 // finalMap = finalMap + coeff* mapToAdd
728 // coeff = convolution of sr and rs
729 // compute the coefficient corresponding to mapToAdd
730 // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
731 // the spectralResponse
732 // different from addComponentBeta(PixelMap<T>& finalMap,
733 // PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
734 // since it permits to use a intBetaMap with a different number of pixels than
735 // the other maps
736
737 int nbpix = finalMap.NbPixels();
738 if (nbpix != mapToAdd.NbPixels())
739 throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
740 double coeff = K;
741
742 if(nbpix != intBetaMap.NbPixels())
743 {
744 for(int i=0; i<finalMap.NbPixels();i++)
745 {
[899]746 double teta,phi,val;
747 finalMap.PixThetaPhi(i, teta, phi);
748 int pixel = intBetaMap.PixIndexSph(teta,phi);
749 val = intBetaMap.PixVal(pixel);
750 finalMap(i) += coeff*mapToAdd(i)*val;
[876]751 }
752 }
753 else
754 {
755 for(int i=0; i<finalMap.NbPixels();i++)
756 {
757 finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
758 }
759 }
760 if (printlev > 1)
761 {
762 cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= " << K << endl;
763 }
764}
765
766
767
768
769
770/* Nouvelle-Fonction */
[1252]771void SKM_MergeFITSKeywords(char * flnm)
772{
773 DVList dvl;
774 FitsFile::FitsExtensionType typeOfExtension;
775 int naxis;
776 vector<int> naxisn;
777 FitsFile::FitsDataType dataType;
778 FitsInFile::GetBlockType(flnm, 1, typeOfExtension, naxis, naxisn, dataType, dvl);
779 // Cleaning the keywords
780#define SZexlst 21
781 char *exlst[SZexlst]=
782 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
783 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
784 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
785 "XTENSION","INSTRUME","TELESCOP"};
786 char kwex[32];
787 int i,l;
788 for (i=0; i<SZexlst; i++) {
789 strncpy(kwex, exlst[i], 32);
790 l = strlen(kwex)-1;
791 if ((kwex[l] != '*') && (kwex[l] != '#')) {
792 dvl.DeleteKey(kwex);
793 }
794 else {
795 bool fgd = (kwex[l] == '#') ? true : false;
796 list<string> lstsup;
797 kwex[l] = '\0';
798 DVList::ValList::const_iterator it;
799 for (it = dvl.Begin(); it != dvl.End(); it++) {
800 if ((*it).first.substr(0,l) != kwex) continue;
801 if (fgd && !isdigit((*it).first[l])) continue;
802 lstsup.push_back((*it).first);
803 }
804 list<string>::iterator it2;
805 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
806 dvl.DeleteKey(*it2);
807 }
808 }
809 dvl_fitskw->Merge(dvl);
810}
[761]811
812/* Nouvelle-Fonction */
[1295]813void SKM_MergeFITSKeywords2(DVList & dvl)
814{
815 // Cleaning the keywords
816#define SZexlst 21
817 char *exlst[SZexlst]=
818 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
819 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
820 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
821 "XTENSION","INSTRUME","TELESCOP"};
822 char kwex[32];
823 int i,l;
824 for (i=0; i<SZexlst; i++) {
825 strncpy(kwex, exlst[i], 32);
826 l = strlen(kwex)-1;
827 if ((kwex[l] != '*') && (kwex[l] != '#')) {
828 dvl.DeleteKey(kwex);
829 }
830 else {
831 bool fgd = (kwex[l] == '#') ? true : false;
832 list<string> lstsup;
833 kwex[l] = '\0';
834 DVList::ValList::const_iterator it;
835 for (it = dvl.Begin(); it != dvl.End(); it++) {
836 if ((*it).first.substr(0,l) != kwex) continue;
837 if (fgd && !isdigit((*it).first[l])) continue;
838 lstsup.push_back((*it).first);
839 }
840 list<string>::iterator it2;
841 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
842 dvl.DeleteKey(*it2);
843 }
844 }
845 dvl_fitskw->Merge(dvl);
846}
847
848/* Nouvelle-Fonction */
[761]849void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
850{
851 char *ntn[2] = {"nu","fnu"};
852 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
853 float xnt[2];
854 double nu;
855 double numin = rs.minFreq();
856 double numax = rs.maxFreq();
857 int nmax = 500;
858 double dnu = (numax-numin)/nmax;
859 for(int k=0; k<nmax; k++) {
860 nu = numin+k*dnu;
861 xnt[0] = nu;
862 xnt[1] = rs.flux(nu);
863 nt.Fill(xnt);
864 }
865 ObjFileIO<NTuple> oiont(nt);
866 oiont.Write(so, name);
867 return;
868}
869
870/* Nouvelle-Fonction */
871void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
872{
873 char *ntn[2] = {"nu","tnu"};
874 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
875 float xnt[2];
876 double nu;
877 double numin = sr.minFreq();
878 double numax = sr.maxFreq();
879 int nmax = 500;
880 double dnu = (numax-numin)/nmax;
881 for(int k=0; k<nmax; k++) {
882 nu = numin+k*dnu;
883 xnt[0] = nu;
884 xnt[1] = sr.transmission(nu);
885 nt.Fill(xnt);
886 }
887 ObjFileIO<NTuple> oiont(nt);
888 oiont.Write(so, name);
889 return;
890}
891
Note: See TracBrowser for help on using the repository browser.