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