Context Navigation

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

Visit:

Last change on this file since 1282 was 1282, checked in by ansari, 25 years ago
Print de numero de version - Reza 2/11/2000
File size: 25.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: