| 1 | #include "machdefs.h"
|
|---|
| 2 | #include <stdlib.h>
|
|---|
| 3 | #include <stdio.h>
|
|---|
| 4 |
|
|---|
| 5 | #include <iostream.h>
|
|---|
| 6 | #include <math.h>
|
|---|
| 7 |
|
|---|
| 8 | #include <string>
|
|---|
| 9 | #include <vector>
|
|---|
| 10 |
|
|---|
| 11 | #include "timing.h"
|
|---|
| 12 | #include "sambainit.h"
|
|---|
| 13 | #include "pexceptions.h"
|
|---|
| 14 | #include "datacards.h"
|
|---|
| 15 | #include "fitsioserver.h"
|
|---|
| 16 |
|
|---|
| 17 | #include "spherehealpix.h"
|
|---|
| 18 | #include "fiospherehealpix.h"
|
|---|
| 19 |
|
|---|
| 20 | #include "radspecvector.h"
|
|---|
| 21 | #include "blackbody.h"
|
|---|
| 22 | #include "nupower.h"
|
|---|
| 23 |
|
|---|
| 24 | #include "squarefilt.h"
|
|---|
| 25 | #include "trianglefilt.h"
|
|---|
| 26 | #include "specrespvector.h"
|
|---|
| 27 | #include "gaussfilt.h"
|
|---|
| 28 |
|
|---|
| 29 | // -----------------------------------------------------------------
|
|---|
| 30 | // ------------- Function declaration ------------------------------
|
|---|
| 31 | int CheckCards(DataCards & dc, string & msg);
|
|---|
| 32 | char * BuildFITSFileName(string const & fname);
|
|---|
| 33 | SpectralResponse * getSpectralResponse(DataCards & dc);
|
|---|
| 34 | RadSpectra * getEmissionSpectra(DataCards & dc, int nc);
|
|---|
| 35 | void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
|
|---|
| 36 | void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
|
|---|
| 37 |
|
|---|
| 38 | template <class T>
|
|---|
| 39 | void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
|
|---|
| 40 | PixelMap<T>& mapToAdd, RadSpectra& rs, double K=1.);
|
|---|
| 41 | template <class T>
|
|---|
| 42 | void MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
|
|---|
| 43 | // -----------------------------------------------------------------
|
|---|
| 44 |
|
|---|
| 45 | // ----- Global (static) variables ------------
|
|---|
| 46 | static bool rdmap = false; // true -> Read map first
|
|---|
| 47 | static char mapPath[256]; // Path for input maps
|
|---|
| 48 | static int hp_nside = 32; // HealPix NSide
|
|---|
| 49 | static int nskycomp = 0; // Number of sky components
|
|---|
| 50 | static int debuglev = 0; // Debug Level
|
|---|
| 51 | static int printlev = 0; // Print Level
|
|---|
| 52 | static POutPersist * so = NULL; // Debug PPFOut file
|
|---|
| 53 |
|
|---|
| 54 | // -------------------------------------------------------------------------
|
|---|
| 55 | // main program
|
|---|
| 56 | // -------------------------------------------------------------------------
|
|---|
| 57 | int main(int narg, char * arg[])
|
|---|
| 58 | {
|
|---|
| 59 | if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
|
|---|
| 60 | cout << " Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
|
|---|
| 61 | exit(0);
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 | InitTim();
|
|---|
| 65 |
|
|---|
| 66 | string msg;
|
|---|
| 67 | int rc = 0;
|
|---|
| 68 | RadSpectra * es = NULL;
|
|---|
| 69 | SpectralResponse * sr = NULL;
|
|---|
| 70 | double moy, sig;
|
|---|
| 71 |
|
|---|
| 72 | DataCards dc;
|
|---|
| 73 | so = NULL;
|
|---|
| 74 |
|
|---|
| 75 | try {
|
|---|
| 76 | string dcard = arg[1];
|
|---|
| 77 | cout << " Decoding parameters from file " << dcard << endl;
|
|---|
| 78 | dc.ReadFile(dcard);
|
|---|
| 79 |
|
|---|
| 80 | rc = CheckCards(dc, msg);
|
|---|
| 81 | if (rc) {
|
|---|
| 82 | cerr << " Error condition -> Rc= " << rc << endl;
|
|---|
| 83 | cerr << " Msg= " << msg << endl;
|
|---|
| 84 | return(rc);
|
|---|
| 85 | }
|
|---|
| 86 | }
|
|---|
| 87 | catch (PException exc) {
|
|---|
| 88 | msg = exc.Msg();
|
|---|
| 89 | cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
|
|---|
| 90 | return(90);
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 |
|
|---|
| 94 | cout << " skymix/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
|
|---|
| 95 | cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
|
|---|
| 96 | << " PrintLev= " << printlev << endl;
|
|---|
| 97 |
|
|---|
| 98 | // We create an output persist file for writing debug objects
|
|---|
| 99 | if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
|
|---|
| 100 |
|
|---|
| 101 | SphereHEALPix<float> outgs(hp_nside);
|
|---|
| 102 | bool okout = false;
|
|---|
| 103 |
|
|---|
| 104 | try {
|
|---|
| 105 | if (rdmap) { // Reading map from FITS file
|
|---|
| 106 | FitsIoServer fios;
|
|---|
| 107 | char ifnm[256];
|
|---|
| 108 | strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
|
|---|
| 109 | ifnm[255] = '\0';
|
|---|
| 110 | cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
|
|---|
| 111 | fios.load(outgs, ifnm, 2);
|
|---|
| 112 | cout << " Output HealPIx Map read - NbPixels= " << outgs.NbPixels() << endl;
|
|---|
| 113 | if (printlev > 0) {
|
|---|
| 114 | MeanSig(outgs.DataBlock(), moy, sig );
|
|---|
| 115 | cout << " MeanSig for outpout map - Mean= " << moy << " Sigma= " << sig << endl;
|
|---|
| 116 | }
|
|---|
| 117 | }
|
|---|
| 118 | else {
|
|---|
| 119 | cout << " Output HealPix Map created - NbPixels= " << outgs.NbPixels() << endl;
|
|---|
| 120 | outgs.SetPixels(0.);
|
|---|
| 121 | }
|
|---|
| 122 |
|
|---|
| 123 | // Decoding detection pass-band filter
|
|---|
| 124 | sr = getSpectralResponse(dc);
|
|---|
| 125 | PrtTim(" After FilterCreation ");
|
|---|
| 126 |
|
|---|
| 127 | FitsIoServer fios; // Our FITS IO Server
|
|---|
| 128 | char * flnm, buff[90];
|
|---|
| 129 | string key;
|
|---|
| 130 |
|
|---|
| 131 | SphereHEALPix<float> ings(hp_nside); // The input map
|
|---|
| 132 | double K = 1.;
|
|---|
| 133 | // Loop over sky component
|
|---|
| 134 | int sk;
|
|---|
| 135 | for(sk = 0; sk<nskycomp; sk++) {
|
|---|
| 136 | cout << " Processing sky component No " << sk+1 << endl;
|
|---|
| 137 | sprintf(buff, "%d", sk+1);
|
|---|
| 138 | key = (string)"MAPFITSFILE" + buff;
|
|---|
| 139 | flnm = BuildFITSFileName(dc.SParam(key, 0));
|
|---|
| 140 | cout << " Reading Input FITS map " << (string)flnm << endl;
|
|---|
| 141 | fios.load(ings, flnm, 2);
|
|---|
| 142 | K = dc.DParam(key, 1, 1.);
|
|---|
| 143 | if (debuglev > 4) { // Writing tne input map to the outppf
|
|---|
| 144 | FIO_SphereHEALPix<float> fiog(ings);
|
|---|
| 145 | fiog.Write(*so, key);
|
|---|
| 146 | }
|
|---|
| 147 | // getting Emission spectra
|
|---|
| 148 | if (es) { delete es; es = NULL; }
|
|---|
| 149 | es = getEmissionSpectra(dc, sk);
|
|---|
| 150 |
|
|---|
| 151 | if (printlev > 2) {
|
|---|
| 152 | MeanSig(ings.DataBlock(), moy, sig );
|
|---|
| 153 | cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
|
|---|
| 154 | }
|
|---|
| 155 | addComponent(*sr, outgs, ings, *es, K);
|
|---|
| 156 | okout = true;
|
|---|
| 157 | if (printlev > 1) {
|
|---|
| 158 | MeanSig(outgs.DataBlock(), moy, sig );
|
|---|
| 159 | cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
|
|---|
| 160 | }
|
|---|
| 161 | sprintf(buff, "End of Proc. Comp. %d ", sk+1);
|
|---|
| 162 | cout << " --------------------------------------- \n" << endl;
|
|---|
| 163 | PrtTim(buff);
|
|---|
| 164 | } // End of sky component loop
|
|---|
| 165 | } // End of try block
|
|---|
| 166 |
|
|---|
| 167 | catch (PException exc) {
|
|---|
| 168 | msg = exc.Msg();
|
|---|
| 169 | cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
|
|---|
| 170 | rc = 50;
|
|---|
| 171 | }
|
|---|
| 172 |
|
|---|
| 173 | // Saving the output map in FITS format
|
|---|
| 174 | if (okout) {
|
|---|
| 175 | FitsIoServer fios;
|
|---|
| 176 | fios.save(outgs, arg[2]);
|
|---|
| 177 | cout << "Output Map (SphereHEALPix<float>) written to FITS file "
|
|---|
| 178 | << (string)(arg[2]) << endl;
|
|---|
| 179 | PrtTim("End of WriteFITS ");
|
|---|
| 180 | // Saving the output map in PPF format
|
|---|
| 181 | if (narg > 3) {
|
|---|
| 182 | POutPersist s(arg[3]);
|
|---|
| 183 | FIO_SphereHEALPix<float> fiog(&outgs) ;
|
|---|
| 184 | fiog.Write(s);
|
|---|
| 185 | cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
|
|---|
| 186 | << (string)(arg[3]) << endl;
|
|---|
| 187 | PrtTim("End of WritePPF ");
|
|---|
| 188 | }
|
|---|
| 189 | }
|
|---|
| 190 | if (so) delete so; // Closing the debug ppf file
|
|---|
| 191 | return(rc);
|
|---|
| 192 | }
|
|---|
| 193 |
|
|---|
| 194 | /* Nouvelle-Fonction */
|
|---|
| 195 | int CheckCards(DataCards & dc, string & msg)
|
|---|
| 196 | // Function to check datacards
|
|---|
| 197 | {
|
|---|
| 198 | rdmap = false;
|
|---|
| 199 | mapPath[0] = '\0';
|
|---|
| 200 | hp_nside = 32;
|
|---|
| 201 | nskycomp = 0;
|
|---|
| 202 | debuglev = 0;
|
|---|
| 203 | printlev = 0;
|
|---|
| 204 |
|
|---|
| 205 | int rc = 0;
|
|---|
| 206 | string key, key2;
|
|---|
| 207 |
|
|---|
| 208 | // Cheking datacards
|
|---|
| 209 | if (dc.NbParam("SKYMIX") < 2) {
|
|---|
| 210 | rc = 71;
|
|---|
| 211 | msg = "Invalid parameters - Check @SKYMIX card ";
|
|---|
| 212 | return(rc);
|
|---|
| 213 | }
|
|---|
| 214 | key = "READMAP";
|
|---|
| 215 | if (dc.HasKey(key)) {
|
|---|
| 216 | if (dc.NbParam(key) < 1) {
|
|---|
| 217 | rc = 72;
|
|---|
| 218 | msg = "Invalid parameters - Check @READMAP card ";
|
|---|
| 219 | return(rc);
|
|---|
| 220 | }
|
|---|
| 221 | else rdmap = true;
|
|---|
| 222 | }
|
|---|
| 223 |
|
|---|
| 224 | // Checking detection filter specification
|
|---|
| 225 | key = "GAUSSFILTER";
|
|---|
| 226 | key2 = "FILTERFITSFILE";
|
|---|
| 227 | if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) {
|
|---|
| 228 | msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
|
|---|
| 229 | rc = 73; return(rc);
|
|---|
| 230 | }
|
|---|
| 231 |
|
|---|
| 232 | // Decoding number of component and pixelisation parameter
|
|---|
| 233 | int mg = 32;
|
|---|
| 234 | int ncomp = 0;
|
|---|
| 235 | ncomp = dc.IParam("SKYMIX", 0, 0);
|
|---|
| 236 | mg = dc.IParam("SKYMIX", 1, 32);
|
|---|
| 237 | if (ncomp < 1) {
|
|---|
| 238 | msg = "Invalid parameters - Check datacards @SKYMIX ";
|
|---|
| 239 | rc = 74;
|
|---|
| 240 | return(rc);
|
|---|
| 241 | }
|
|---|
| 242 |
|
|---|
| 243 | // Checking detection filter specification
|
|---|
| 244 | // Checking input FITS file specifications
|
|---|
| 245 | int kc;
|
|---|
| 246 | char buff[256];
|
|---|
| 247 | bool pb = false;
|
|---|
| 248 | string key3;
|
|---|
| 249 | for(kc=0; kc<ncomp; kc++) {
|
|---|
| 250 | sprintf(buff, "MAPFITSFILE%d", kc+1);
|
|---|
| 251 | key = buff;
|
|---|
| 252 | if (dc.NbParam(key) < 1) {
|
|---|
| 253 | msg = "Missing or invalid card : " + key;
|
|---|
| 254 | pb = true; break;
|
|---|
| 255 | }
|
|---|
| 256 | sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
|
|---|
| 257 | key = buff;
|
|---|
| 258 | sprintf(buff, "BLACKBODY%d", kc+1);
|
|---|
| 259 | key2 = buff;
|
|---|
| 260 | sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
|
|---|
| 261 | key3 = buff;
|
|---|
| 262 | if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6)) {
|
|---|
| 263 | msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3;
|
|---|
| 264 | pb = true; break;
|
|---|
| 265 | }
|
|---|
| 266 |
|
|---|
| 267 | }
|
|---|
| 268 |
|
|---|
| 269 | if (pb) {
|
|---|
| 270 | rc = 75;
|
|---|
| 271 | return(75);
|
|---|
| 272 | }
|
|---|
| 273 |
|
|---|
| 274 |
|
|---|
| 275 | // Initialiazing parameters
|
|---|
| 276 | rc = 0;
|
|---|
| 277 | msg = "OK";
|
|---|
| 278 | nskycomp = ncomp;
|
|---|
| 279 | hp_nside = mg;
|
|---|
| 280 |
|
|---|
| 281 | // Checking for PATH definition card
|
|---|
| 282 | key = "MAPPATH";
|
|---|
| 283 | if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
|
|---|
| 284 | mapPath[255] = '\0';
|
|---|
| 285 | key = "DEBUGLEVEL";
|
|---|
| 286 | debuglev = dc.IParam(key, 0, 0);
|
|---|
| 287 | key = "PRINTLEVEL";
|
|---|
| 288 | printlev = dc.IParam(key, 0, 0);
|
|---|
| 289 | return(rc);
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | static char buff_flnm[1024]; // Mal protege !
|
|---|
| 293 | /* Nouvelle-Fonction */
|
|---|
| 294 | char* BuildFITSFileName(string const & fname)
|
|---|
| 295 | {
|
|---|
| 296 | if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
|
|---|
| 297 | else sprintf(buff_flnm, "%s", fname.c_str());
|
|---|
| 298 | return(buff_flnm);
|
|---|
| 299 | }
|
|---|
| 300 |
|
|---|
| 301 | /* Nouvelle-Fonction */
|
|---|
| 302 | SpectralResponse * getSpectralResponse(DataCards & dc)
|
|---|
| 303 | {
|
|---|
| 304 | SpectralResponse * filt = NULL;
|
|---|
| 305 |
|
|---|
| 306 | string key = "FILTERFITSFILE";
|
|---|
| 307 | string key2 = "GAUSSFILTER";
|
|---|
| 308 | string ppfname = "filter";
|
|---|
| 309 |
|
|---|
| 310 | if (dc.HasKey(key) ) { // Reading FITS filter file
|
|---|
| 311 | FitsIoServer fios;
|
|---|
| 312 | char ifnm[256];
|
|---|
| 313 | strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
|
|---|
| 314 | ifnm[255] = '\0';
|
|---|
| 315 | Matrix mtx(2,10);
|
|---|
| 316 | fios.load(mtx, ifnm);
|
|---|
| 317 | double numin = dc.DParam(key, 2, 1.);
|
|---|
| 318 | double numax = dc.DParam(key, 3, 9999.);
|
|---|
| 319 | Vector nu(mtx.NCols());
|
|---|
| 320 | Vector fnu(mtx.NCols());
|
|---|
| 321 | for(int k=0; k<mtx.NCols(); k++) {
|
|---|
| 322 | nu(k) = mtx(0, k);
|
|---|
| 323 | fnu(k) = mtx(1, k);
|
|---|
| 324 | }
|
|---|
| 325 | filt = new SpecRespVec(nu, fnu, numin, numax);
|
|---|
| 326 | ppfname = key;
|
|---|
| 327 | }
|
|---|
| 328 | else if (dc.HasKey(key2) ) { // creating GaussianFilter
|
|---|
| 329 | double nu0 = dc.DParam(key2, 0, 10.);
|
|---|
| 330 | double s = dc.DParam(key2, 1, 1.);
|
|---|
| 331 | double a = dc.DParam(key2, 2, 1.);
|
|---|
| 332 | double numin = dc.DParam(key2, 3, 0.1);
|
|---|
| 333 | double numax = dc.DParam(key2, 4, 9999);
|
|---|
| 334 | filt = new GaussianFilter(nu0, s, a, numin, numax);
|
|---|
| 335 | ppfname = key2;
|
|---|
| 336 | }
|
|---|
| 337 | if (filt == NULL) throw ParmError("datacard error ! No detection filter");
|
|---|
| 338 | cout << " Filter decoded - Created " << endl;
|
|---|
| 339 | cout << *filt << endl;
|
|---|
| 340 |
|
|---|
| 341 | // for debug
|
|---|
| 342 | if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
|
|---|
| 343 | return(filt);
|
|---|
| 344 | }
|
|---|
| 345 |
|
|---|
| 346 | /* Nouvelle-Fonction */
|
|---|
| 347 | RadSpectra * getEmissionSpectra(DataCards & dc, int nc)
|
|---|
| 348 | {
|
|---|
| 349 | char numb[16];
|
|---|
| 350 | sprintf(numb, "%d", nc+1);
|
|---|
| 351 |
|
|---|
| 352 | string key = (string)"SPECTRAFITSFILE" + numb;
|
|---|
| 353 | string key2 = (string)"BLACKBODY" + numb;
|
|---|
| 354 | string key3 = (string)"POWERLAWSPECTRA" + numb;
|
|---|
| 355 | string ppfname = "espectra";
|
|---|
| 356 |
|
|---|
| 357 | RadSpectra * rs = NULL;
|
|---|
| 358 | if (dc.HasKey(key) ) { // Reading emission spectra from file
|
|---|
| 359 | char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
|
|---|
| 360 | cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
|
|---|
| 361 | FitsIoServer fios;
|
|---|
| 362 | Matrix mtx(2,10);
|
|---|
| 363 | fios.load(mtx, ifnm);
|
|---|
| 364 | double numin = dc.DParam(key, 2, 1.);
|
|---|
| 365 | double numax = dc.DParam(key, 3, 9999.);
|
|---|
| 366 | Vector nu(mtx.NCols());
|
|---|
| 367 | Vector tnu(mtx.NCols());
|
|---|
| 368 | for(int k=0; k<mtx.NCols(); k++) {
|
|---|
| 369 | nu(k) = mtx(0, k);
|
|---|
| 370 | tnu(k) = mtx(1, k);
|
|---|
| 371 | }
|
|---|
| 372 | rs = new RadSpectraVec(nu, tnu, numin, numax);
|
|---|
| 373 | ppfname = key;
|
|---|
| 374 | }
|
|---|
| 375 | else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
|
|---|
| 376 | rs = new BlackBody(dc.DParam(key2, 0, 2.726));
|
|---|
| 377 | ppfname = key2;
|
|---|
| 378 | }
|
|---|
| 379 | else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
|
|---|
| 380 | double a = dc.DParam(key3, 0, 1.);
|
|---|
| 381 | double nu0 = dc.DParam(key3, 1, 100.);
|
|---|
| 382 | double dnu = dc.DParam(key3, 2, 10.);
|
|---|
| 383 | double b = dc.DParam(key3, 3, 0.);
|
|---|
| 384 | double numin = dc.DParam(key3, 4, 0.1);
|
|---|
| 385 | double numax = dc.DParam(key3, 5, 9999);
|
|---|
| 386 | rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
|
|---|
| 387 | ppfname = key3;
|
|---|
| 388 | }
|
|---|
| 389 | if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
|
|---|
| 390 | cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
|
|---|
| 391 | cout << *rs << endl;
|
|---|
| 392 | // for debug
|
|---|
| 393 | if (debuglev > 2) RadSpec2Nt(*rs, *so, ppfname);
|
|---|
| 394 | return(rs);
|
|---|
| 395 | }
|
|---|
| 396 |
|
|---|
| 397 |
|
|---|
| 398 |
|
|---|
| 399 |
|
|---|
| 400 | /* Nouvelle-Fonction */
|
|---|
| 401 | template <class T>
|
|---|
| 402 | void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
|
|---|
| 403 | PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
|
|---|
| 404 | {
|
|---|
| 405 | // finalMap = finalMap + coeff* mapToAdd
|
|---|
| 406 | // coeff = convolution of sr and rs
|
|---|
| 407 | // compute the coefficient corresponding to mapToAdd
|
|---|
| 408 | if (finalMap.NbPixels() != mapToAdd.NbPixels())
|
|---|
| 409 | throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
|
|---|
| 410 | double coeff = rs.filteredIntegratedFlux(sr) * K;
|
|---|
| 411 | if (printlev > 1)
|
|---|
| 412 | cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
|
|---|
| 413 | for(int i=0; i<finalMap.NbPixels(); i++)
|
|---|
| 414 | {
|
|---|
| 415 | finalMap(i) += coeff * mapToAdd(i);
|
|---|
| 416 | }
|
|---|
| 417 | }
|
|---|
| 418 |
|
|---|
| 419 | /* Nouvelle-Fonction */
|
|---|
| 420 | template <class T>
|
|---|
| 421 | void MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig)
|
|---|
| 422 |
|
|---|
| 423 | {
|
|---|
| 424 | gmoy=0.;
|
|---|
| 425 | gsig = 0.;
|
|---|
| 426 | double valok;
|
|---|
| 427 | for(int k=0; k<(int)dbl.Size(); k++) {
|
|---|
| 428 | valok = dbl(k);
|
|---|
| 429 | gmoy += valok; gsig += valok*valok;
|
|---|
| 430 | }
|
|---|
| 431 | gmoy /= (double)dbl.Size();
|
|---|
| 432 | gsig = gsig/(double)dbl.Size() - gmoy*gmoy;
|
|---|
| 433 | if (gsig >= 0.) gsig = sqrt(gsig);
|
|---|
| 434 | }
|
|---|
| 435 |
|
|---|
| 436 | /* Nouvelle-Fonction */
|
|---|
| 437 | void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
|
|---|
| 438 | {
|
|---|
| 439 | char *ntn[2] = {"nu","fnu"};
|
|---|
| 440 | NTuple nt(2,ntn); // Creation NTuple (AVEC new )
|
|---|
| 441 | float xnt[2];
|
|---|
| 442 | double nu;
|
|---|
| 443 | double numin = rs.minFreq();
|
|---|
| 444 | double numax = rs.maxFreq();
|
|---|
| 445 | int nmax = 500;
|
|---|
| 446 | double dnu = (numax-numin)/nmax;
|
|---|
| 447 | for(int k=0; k<nmax; k++) {
|
|---|
| 448 | nu = numin+k*dnu;
|
|---|
| 449 | xnt[0] = nu;
|
|---|
| 450 | xnt[1] = rs.flux(nu);
|
|---|
| 451 | nt.Fill(xnt);
|
|---|
| 452 | }
|
|---|
| 453 | ObjFileIO<NTuple> oiont(nt);
|
|---|
| 454 | oiont.Write(so, name);
|
|---|
| 455 | return;
|
|---|
| 456 | }
|
|---|
| 457 |
|
|---|
| 458 | /* Nouvelle-Fonction */
|
|---|
| 459 | void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
|
|---|
| 460 | {
|
|---|
| 461 | char *ntn[2] = {"nu","tnu"};
|
|---|
| 462 | NTuple nt(2,ntn); // Creation NTuple (AVEC new )
|
|---|
| 463 | float xnt[2];
|
|---|
| 464 | double nu;
|
|---|
| 465 | double numin = sr.minFreq();
|
|---|
| 466 | double numax = sr.maxFreq();
|
|---|
| 467 | int nmax = 500;
|
|---|
| 468 | double dnu = (numax-numin)/nmax;
|
|---|
| 469 | for(int k=0; k<nmax; k++) {
|
|---|
| 470 | nu = numin+k*dnu;
|
|---|
| 471 | xnt[0] = nu;
|
|---|
| 472 | xnt[1] = sr.transmission(nu);
|
|---|
| 473 | nt.Fill(xnt);
|
|---|
| 474 | }
|
|---|
| 475 | ObjFileIO<NTuple> oiont(nt);
|
|---|
| 476 | oiont.Write(so, name);
|
|---|
| 477 | return;
|
|---|
| 478 | }
|
|---|
| 479 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
