Changeset 904 in Sophya for trunk/SophyaProg/PMixer/skymixer.cc


Ignore:
Timestamp:
Apr 13, 2000, 12:04:40 PM (25 years ago)
Author:
ansari
Message:

Sophie: adding the dipole

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaProg/PMixer/skymixer.cc

    r899 r904  
    5959void integratedMap(SpectralResponse&  sr,   
    6060                   SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
    61   //template <class T>
    62   //void integratedMap(SpectralResponse&  sr,   
    63   //                          PixelMap<T>& betaMap,
    64   //                          double normFreq,
    65   //                          PixelMap<T>& outMap);
     61
     62//
    6663template <class T>
    6764void addComponentBeta(SphereHEALPix<T>& finalMap,
    6865                      SphereHEALPix<T>& mapToAdd,
    6966                      SphereHEALPix<T>& intBetaMap, double K);
     67//
     68template <class T>
     69void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap,
     70               double theta,double phi,double amp,double temp);
    7071//
    7172// -----------------------------------------------------------------
    … …  
    8990    exit(0);
    9091  }
    91 
    92 InitTim();
    93 
    94 string msg;
    95 int rc = 0;
    96 RadSpectra * es = NULL;
    97 SpectralResponse * sr = NULL;
    98 double moy, sig;
    99 
    100 DataCards dc;
    101 so = NULL;
    102 
    103  try {
    104    string dcard = arg[1];
    105    if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
    106    dc.ReadFile(dcard);
    107  
    108    rc = CheckCards(dc, msg);
    109    if (rc) {
    110      cerr << " Error condition -> Rc= " << rc << endl;
    111      cerr << " Msg= " << msg << endl;
    112      return(rc);
    113    }
    114  }
    115  catch (PException exc) {
    116    msg = exc.Msg();
    117    cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
    118    return(90);
    119  }   
    120 
    121 
    122 cout << " skymix/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
    123 cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev
    124      << "  PrintLev= " << printlev << endl;
    125 
    126 // We create an output persist file for writing debug objects
    127 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
    128 
    129  SphereHEALPix<float> outgs(hp_nside);
    130  try{
    131    if (rdmap) {  // Reading map from FITS file
    132      char ifnm[256];
    133      strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
    134      ifnm[255] = '\0';
    135      cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
    136      {
    137        FITS_SphereHEALPix<float> fios(ifnm);
    138        outgs = (SphereHEALPix<float>)fios;
    139      }
    140      if(printlev>0)
    141        cout << " Output HealPIx Map read - NbPixels= " <<
    142          outgs.NbPixels() << endl;
     92 
     93  InitTim();
     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;
     103 
     104  try {
     105    string dcard = arg[1];
     106    if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
     107    dc.ReadFile(dcard);
     108   
     109    rc = CheckCards(dc, msg);
     110    if (rc) {
     111      cerr << " Error condition -> Rc= " << rc << endl;
     112      cerr << " Msg= " << msg << endl;
     113      return(rc);
     114    }
     115  }
     116  catch (PException exc) {
     117    msg = exc.Msg();
     118    cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
     119    return(90);
     120  }   
     121 
     122 
     123  cout << " skymix/Info : NComp = " <<  nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
     124  cout << "  ... MapPath = " << (string)mapPath << "  DebugLev= " << debuglev
     125       << "  PrintLev= " << printlev << endl;
     126 
     127  // We create an output persist file for writing debug objects
     128  if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
     129 
     130  SphereHEALPix<float> outgs(hp_nside);
     131  try{
     132    if (rdmap) {  // Reading map from FITS file
     133      char ifnm[256];
     134      strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
     135      ifnm[255] = '\0';
     136      cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
     137      {
     138        FITS_SphereHEALPix<float> fios(ifnm);
     139        outgs = (SphereHEALPix<float>)fios;
     140      }
     141      if(printlev>0)
     142        cout << " Output HealPIx Map read - NbPixels= " <<
     143          outgs.NbPixels() << endl;
    143144      if (printlev > 0) {
    144145        MeanSig(outgs.DataBlock(), moy, sig );
    … …  
    161162    string key;
    162163   
    163     //    SphereHEALPix<float> ings(hp_nside);  // The input map
    164164    double K = 1.;
    165165    double freqOfMap = 1.;
    … …  
    170170      cout << " Processing sky component No " << sk+1 << endl;
    171171     
     172
    172173      sprintf(buff, "%d", sk+1);
    173       key = (string)"MAPFITSFILE" + buff;
    174       flnm = BuildFITSFileName(dc.SParam(key, 0));
    175      
    176       K = dc.DParam(key, 1, 1.);
    177      
    178      
    179       cout << " Reading Input FITS map " << (string)flnm << endl;
    180       SphereHEALPix<float> ings(hp_nside);
    181       {
    182         FITS_SphereHEALPix<float> fiosIn(flnm);
    183         ings = (SphereHEALPix<float>)fiosIn;
    184       }
    185       if (debuglev > 4) {  // Writing the input map to the outppf
    186         FIO_SphereHEALPix<float> fiog(ings);
    187         fiog.Write(*so, key);
    188       }
    189       if (printlev > 2) {
    190         MeanSig(ings.DataBlock(), moy, sig );
    191         cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
    192       }
    193       bool mapDependentOfFreq = false;
    194       key = (string)"BETAFITSFILE"+ buff;
     174      key = (string)"DIPOLE" + buff;
     175      // check for an eventual dipole
    195176      if(dc.HasKey(key))
    196177        {
    197           mapDependentOfFreq = true;
    198         }
    199 
    200       // getting Emission spectra 
    201       if(!mapDependentOfFreq)
    202         {
    203178          if (es) { delete es; es = NULL; }
    204           es = getEmissionSpectra(dc, sk);
    205           addComponent(*sr, outgs, ings, *es, K);
     179          double temp = -10.;
     180          double theta=  dc.DParam(key,1,1.);
     181          double phi  =  dc.DParam(key,2,1.);
     182          double amp  =  dc.DParam(key,3,1.);
     183          if(dc.NbParam(key)>3)
     184            {
     185              temp =  dc.DParam(key,4,1.);
     186            }
     187          cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
     188          addDipole(*sr, outgs,theta,phi,amp,temp);
    206189        }
    207190      else
    208191        {
    209           key = (string)"BETAFITSFILE"+ buff;
    210           //SphereHEALPix<float> betaMap;
     192          sprintf(buff, "%d", sk+1);
     193          key = (string)"MAPFITSFILE" + buff;
    211194          flnm = BuildFITSFileName(dc.SParam(key, 0));
    212           double normFreq = dc.DParam(key, 1, 1.);
    213           if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
    214           int nSideForInt = dc.DParam(key, 2, 1.);
    215           if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
    216           cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
    217           SphereHEALPix<float> betaMap(hp_nside);
     195         
     196          K = dc.DParam(key, 1, 1.);
     197         
     198          cout << " Reading Input FITS map " << (string)flnm << endl;
     199          SphereHEALPix<float> ings(hp_nside);
    218200          {
    219             FITS_SphereHEALPix<float> fiosBM(flnm);
    220             betaMap = (SphereHEALPix<float>)fiosBM;
     201            FITS_SphereHEALPix<float> fiosIn(flnm);
     202            ings = (SphereHEALPix<float>)fiosIn;
    221203          }
    222204          if (debuglev > 4) {  // Writing the input map to the outppf
    223             FIO_SphereHEALPix<float> fiogs(betaMap);
    224             fiogs.Write(*so, key);
     205            FIO_SphereHEALPix<float> fiog(ings);
     206            fiog.Write(*so, key);
    225207          }
    226 
    227           if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
    228           bool bydefault = true;
    229           if(!bydefault)
    230             addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
     208          if (printlev > 2) {
     209            MeanSig(ings.DataBlock(), moy, sig );
     210            cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
     211          }
     212          bool mapDependentOfFreq = false;
     213          key = (string)"BETAFITSFILE"+ buff;
     214          if(dc.HasKey(key))
     215            {
     216              mapDependentOfFreq = true;
     217            }
     218         
     219          // getting Emission spectra 
     220          if(!mapDependentOfFreq)
     221            {
     222              if (es) { delete es; es = NULL; }
     223              es = getEmissionSpectra(dc, sk);
     224              addComponent(*sr, outgs, ings, *es, K);
     225            }
    231226          else
    232227            {
    233               // integrate the betamap over the SpectralResponse
    234               SphereHEALPix<float> intBetaMap(nSideForInt);
    235               integratedMap(*sr, betaMap, normFreq, intBetaMap);
     228              key = (string)"BETAFITSFILE"+ buff;
     229              //SphereHEALPix<float> betaMap;
     230              flnm = BuildFITSFileName(dc.SParam(key, 0));
     231              double normFreq = dc.DParam(key, 1, 1.);
     232              if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
     233              int nSideForInt = dc.DParam(key, 2, 1.);
     234              if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
     235              cout << "....BetaFits...  Reading Beta FITS map " << (string)flnm << endl;
     236              SphereHEALPix<float> betaMap(hp_nside);
     237              {
     238                FITS_SphereHEALPix<float> fiosBM(flnm);
     239                betaMap = (SphereHEALPix<float>)fiosBM;
     240              }
    236241              if (debuglev > 4) {  // Writing the input map to the outppf
    237                 FIO_SphereHEALPix<float> fiogs2(intBetaMap);
    238                 fiogs2.Write(*so, "INTBETAMAP");
     242                FIO_SphereHEALPix<float> fiogs(betaMap);
     243                fiogs.Write(*so, key);
    239244              }
    240 
    241               betaMap.Resize(8);
    242               MeanSig(intBetaMap.DataBlock(), moy, sig );
    243               if (printlev > 4) cout << "....BetaFits... MeanSig for intBetaMap - Mean= " << moy << " Sigma= " << sig << endl;
    244               // add the integrated beta map
    245               addComponentBeta(outgs,ings,intBetaMap, K);
     245             
     246              if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
     247              bool bydefault = true;
     248              if(!bydefault)
     249                addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
     250              else
     251                {
     252                  // integrate the betamap over the SpectralResponse
     253                  SphereHEALPix<float> intBetaMap(nSideForInt);
     254                  integratedMap(*sr, betaMap, normFreq, intBetaMap);
     255                  if (debuglev > 4) {  // Writing the input map to the outppf
     256                    FIO_SphereHEALPix<float> fiogs2(intBetaMap);
     257                    fiogs2.Write(*so, "INTBETAMAP");
     258                  }
     259                 
     260                  betaMap.Resize(8);
     261                  MeanSig(intBetaMap.DataBlock(), moy, sig );
     262                  if (printlev > 4)
     263                    cout << "....BetaFits... MeanSig for intBetaMap - Mean= "
     264                         << moy << " Sigma= " << sig << endl;
     265                  // add the integrated beta map
     266                  addComponentBeta(outgs,ings,intBetaMap, K);
     267                }
    246268            }
    247         }
    248      
    249       MeanSig(outgs.DataBlock(), moy, sig );
    250       cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
    251       cout << "-------------------------------------------------" << endl;
    252      
    253       sprintf(buff, "End of Processing Component %d ", sk+1);
    254       PrtTim(buff);
     269         
     270          MeanSig(outgs.DataBlock(), moy, sig );
     271          cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
     272          cout << "-------------------------------------------------" << endl;
     273         
     274          sprintf(buff, "End of Processing Component %d ", sk+1);
     275          PrtTim(buff);
     276        }
     277    }
     278  }
     279  catch(PException exc)
     280    {
     281      cout << "catched PException" << endl;
     282      msg = exc.Msg();
     283      cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
     284      rc = 50;
     285      return(50);
    255286    }
    256  }
    257  catch(PException exc)
    258    {
    259      cout << "catched PException" << endl;
    260      msg = exc.Msg();
    261      cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
    262      rc = 50;
    263      return(50);
    264    }
    265  
     287 
    266288 // Saving the output map in FITS format
    267289 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
    … …  
    297319
    298320int rc = 0;
    299 string key, key2;
     321string key, key2,key3;
    300322
    301323  //  Cheking datacards
    … …  
    318340  key = "GAUSSFILTER";
    319341  key2 = "FILTERFITSFILE";
    320   if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) {
    321     msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
     342  key3 = "DIPOLE";
     343  if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) {
     344    msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
    322345    rc = 73;  return(rc);
    323346    }
    … …  
    339362  char buff[256];
    340363  bool pb = false;
    341   string key3;
    342364  string key4;
    343365  string key5;
     366  string key6;
    344367  for(kc=0; kc<ncomp; kc++) {
    345368    sprintf(buff, "MAPFITSFILE%d", kc+1);
    346369    key = buff;
    347     if (dc.NbParam(key) < 1 ) {   
    348       msg = "Missing or invalid card : " + key + " " + key2 ;
     370    sprintf(buff, "DIPOLE%d", kc+1);
     371    key3 = buff;
     372    if (dc.NbParam(key) < 1 && dc.NbParam(key3)<1) {   
     373      msg = "Missing or invalid card : " + key + " " + key2 + " " + key3;
    349374      pb = true;  break;
    350375    }
    … …  
    359384    sprintf(buff, "DIPOLE%d", kc+1);
    360385    key5 = buff;
     386    sprintf(buff, "DERIVBB%d", kc+1);
     387    key6 = buff;
    361388    if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
    362          && (dc.NbParam(key5)<1))  {   
     389         &&  (dc.NbParam(key6)<1) &&  (dc.NbParam(key5)<3))  {   
    363390      msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
    364391      pb = true;  break;
    … …  
    455482  string key = (string)"SPECTRAFITSFILE" + numb;
    456483  string key2 = (string)"BLACKBODY" + numb;
    457   string key5 = (string)"DIPOLE" + numb;
     484  string key5 = (string)"DERIVBB" + numb;
    458485  string key3 = (string)"POWERLAWSPECTRA" + numb;
    459486  string ppfname = "espectra";
    … …  
    482509  }
    483510  else if (dc.HasKey(key5) ) { // Creating Dipole
    484     rs = new DerivBlackBody(dc.DParam(key5, 0, 2.726));
     511    rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
    485512    ppfname = key5;
    486513  }
    … …  
    506533
    507534
     535template <class T>
     536void addDipole(SpectralResponse&  sr,  PixelMap<T>& finalMap,
     537               double theta,double phi,double amp,double temp)
     538{
     539  DerivBlackBody dbb;
     540  if(temp>0) dbb.setTemperature(temp);
     541  double coeff = dbb.filteredIntegratedFlux(sr) * amp;
     542  UnitVector vd(theta,phi);
     543  UnitVector vc(theta,phi);
     544 
     545  for(int i=0; i<finalMap.NbPixels(); i++)
     546    {
     547      double thetar,phir;
     548      finalMap.PixThetaPhi(i,thetar,phir);     
     549      vc.SetThetaPhi(thetar,  phir);
     550      finalMap(i) += vd.Psc(vc)*coeff;
     551    }
     552  if (debuglev > 4) {  // Writing the input map to the outppf
     553    SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
     554    for(int i=0; i<finalMap.NbPixels(); i++)
     555      {
     556        double thetar,phir;
     557        finalMap.PixThetaPhi(i,thetar,phir);     
     558        vc.SetThetaPhi(thetar,  phir);
     559        ings(i) = vd.Psc(vc);
     560      }
     561    FIO_SphereHEALPix<float> fiog(ings);
     562    fiog.Write(*so, "dipole");
     563    cout << "Debug the dipole map....saved in debug file !" << endl;
     564  }
     565}
    508566/* Nouvelle-Fonction */
    509567template <class T>
Note: See TracChangeset for help on using the changeset viewer.