Context Navigation

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

Visit:

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

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

Download in other formats: