Context Navigation

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

Visit:

Last change on this file since 2682 was 2615, checked in by cmv, 21 years ago
using namespace sophya enleve de machdefs.h, nouveau sopnamsp.h cmv 10/09/2004
File size: 31.6 KB

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

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

Download in other formats: