Changeset 3829 in Sophya for trunk/Cosmo


Ignore:
Timestamp:
Aug 3, 2010, 6:38:00 PM (15 years ago)
Author:
ansari
Message:

Adaptation programme srcat2cube.cc pour utilisation catalogue North20cm avec indice spectral, Reza 03/08/2010

Location:
trunk/Cosmo/RadioBeam
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/README

    r3825 r3829  
    2121  - applobe.cc : Antenna/array beam effect on 3D sky cube
    2222  - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
    23   - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
     23  - srcat2cube.cc : radio source 3D temperature cube from NVSS or North20cm catalogs
    2424  - gsm2cube.cc : Computes 3D temperature cube from a set of GSM (Global Sky Model) healpix PPF maps 
    2525  - tjyk.cc : test de la classe H21Conversions

  • trunk/Cosmo/RadioBeam/calcpk.cc

    r3825 r3829  
    3636{
    3737  if (narg<3) {
    38     cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
     38    cout << " Usage: calcpk In3DMap_PPFName OutPk_PPFName  [InRenormFactor] [PixNoiseLevel] [Four3D_PPFName]" << endl;
    3939    return 1;
    4040  }
    … …  
    5454    if (narg>4) {
    5555      pixsignoise=atof(arg[4]);
    56       fgaddnoise=true;
     56      if (pixsignoise>1.e-9)  fgaddnoise=true;
    5757    }
     58    bool fgsavef3d=false;
     59    string outf3dname="four3d.ppf";
     60    if (narg>5) {
     61      outf3dname=arg[5];
     62      fgsavef3d=true;
     63    }
     64
    5865
    5966    cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
    … …  
    8289    TArray< complex<r_4> > four3d;
    8390    ffts.FFTForward(inmap, four3d);
     91    cout << " calcpk[2.b] inmap.Norm2()=" << inmap.Norm2() << " four3d.Norm2()=" << four3d.Norm2() << endl;
     92
    8493    tm.Split(" After FFTForward ");
    8594
    … …  
    92101    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
    93102
     103
    94104    HProf hp = pkc.ComputePk(0.,256);
    95105    {
    … …  
    97107      POutPersist po(outname);
    98108      po << hp;
    99       outname = "f3d_" + outname;
    100       cout << "calcpk[4.b] : writing four3d to PPF file  " << outname << endl;
    101       POutPersist pof(outname);
     109    }
     110    if (fgsavef3d)   {
     111      cout << "calcpk[4.b] : writing four3d to PPF file  " << outf3dname << endl;
     112      POutPersist pof(outf3dname);
    102113      pof << four3d;
    103114    }

  • trunk/Cosmo/RadioBeam/cubedef.h

    r3826 r3829  
    1111static sa_size_t NFreq = 256;
    1212
    13 static double Theta0Degre = 75.;  // -> Delta = +30 deg
     13static double Theta0Degre = 65.;  // -> Delta = +30 deg
    1414static double Phi0Degre = 135.;   // -> alpha = 9h00
    1515static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
    … …  
    4545//---  Parametres des lois de puissance en frequence
    4646static double AmpPL1 = 1.;   // amp max PowerLaw 1 (synchrotron 
    47 static double PLidx1 = -2.5;  // index de la loi de puissance synchrotron
     47static double PLidx1 = -2.8;  // index de la loi de puissance synchrotron
    4848static double sigPLidx1 = 0.1;  // Sigma de la variation (gaussienne) de index1
    4949// Amplitude max de la 2eme composante en loi de puissance (tirage plat 0 ... AmpPL2)

  • trunk/Cosmo/RadioBeam/radutil.cc

    r3783 r3829  
    1010  setFrequency(freq);
    1111  setOmegaPix(opix);
     12  setCosmoParam();
    1213}
    1314
    … …  
    2829  return jy*1e-26 / omegapix_ * lambda_ * lambda_ / 2. / k_Boltzman_Cst ;
    2930}
     31
     32
     33double H21Conversions::setCosmoParam(double omegamatter, double omegabaryon, double h100, double fracHI)
     34{
     35  OmegaMatter_=omegamatter;
     36  OmegaBaryons_=omegabaryon;
     37  h100_=h100;
     38  fracHI_=fracHI;
     39  OmegaLambda_=1.-OmegaMatter_;  // Je neglige OmegaRadiation
     40  return OmegaLambda_;
     41}
     42
     43double H21Conversions::Mean21cmTemperature_mK()
     44{
     45  double zz = 1.+z_;
     46  double cc=zz*zz/sqrt(OmegaMatter_*zz*zz*zz+OmegaLambda_);
     47  cc *= ((h100_/0.7)*(OmegaBaryons_/0.04)*(fracHI_/0.1));
     48  return (cc*0.57);
     49}

  • trunk/Cosmo/RadioBeam/radutil.h

    r3783 r3829  
    1818  double toKelvin(double jy);    // Conversion de Jansky en  temperature (Kelvin)
    1919
     20  double Mean21cmTemperature_mK();      // Temperature moyenne de l'emission a 21 cm en mK
     21  inline double T21cm_Kelvin() { return Mean21cmTemperature_mK()/1000.;  }
     22  inline double T21cm_mK() { return Mean21cmTemperature_mK();  }
     23
    2024  void setFrequency(double nu);  // on definit la frequence en MHz
    2125  inline void setRedshift(double z)    // on definit le redshift
    … …  
    3135  { double cf=Angle(1.,Angle::ArcMin).ToRadian(); omegapix_ = opix*cf*cf; }
    3236
     37  // Definition des parametres cosmologiques utiles pour le calcul de la temperature d'emission a 21 cm
     38  // retourne la valeur de OmegaLambda (univers plat)
     39  double setCosmoParam(double omegamatter=0.02581, double omegabaryon=0.0441, double h100=0.719, double fracHI=0.02);
     40  inline void setFracHI(double fracHI=0.02) { fracHI_=fracHI; }
     41
    3342  inline double getRedshift() { return z_; }
    3443  inline double getFrequency() { return freq_; }
    … …  
    4251  { double cf=Angle(1.,Angle::Degree).ToRadian(); return omegapix_/cf/cf; }
    4352
     53
    4454  static double SpeedOfLight_Cst;       // Speed of light  m/sec
    4555  static double Freq021cm_Cst;       // Speed of light  m/sec
    … …  
    5060  double lambda_;
    5161  double omegapix_;
     62
     63  // Parametres cosmologiques pour calcul du coefficient de conversion Mass to T21
     64  double OmegaMatter_;
     65  double OmegaBaryons_;
     66  double OmegaLambda_;
     67  double h100_;
     68  double fracHI_;
    5269};
    5370

  • trunk/Cosmo/RadioBeam/srcat2cube.cc

    r3825 r3829  
    3434
    3535#include "cubedef.h"
     36
     37//----------------
     38int nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& cube);
     39int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& cube);
    3640
    3741//----------------------------------------------------------------------------
    … …  
    4852  Timer tm("srcat2cube");
    4953
    50   if (narg < 3) {
    51     cout << "Usage: srccat2cube NVSS_CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
     54  if (narg < 4) {
     55    cout << "Usage: srccat2cube -nvss/-north20 CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
    5256    return 1;
    5357  }
    … …  
    5559
    5660  // decodage arguments
    57   string outname = arg[2];
    58   string inname = arg[1];
     61 
     62  string copt=arg[1];
     63  string outname=arg[3];
     64  string inname=arg[2];
    5965  int rc = 91;
    6066
    61   cout << " ====== srccat2cube :   Input NVSS catalog name= " << inname << " OutName=" << outname;
     67  cout << " ====== srccat2cube : Input catalog name= " << inname << " OutName=" << outname;
    6268  bool fginmap=true;
    6369  try {
    64     DataTable nvss;
    65     cout << "srccat2cube[1]: reading NVSS catalog from " << inname << endl;
     70    DataTable cat;
     71    cout << "srccat2cube[1]: reading source catalog from " << inname << endl;
    6672    {
    6773      FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
    68       fis >> nvss;
    69     }
    70     cout << nvss;
    71     sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
    72     sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
    73     sa_size_t idxf = nvss.IndexNom("S1_4");
    74     sa_size_t idxmajax = nvss.IndexNom("MajAxis");
    75     sa_size_t idxminax = nvss.IndexNom("MinAxis");
    76 
    77     cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf
    78          << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
    79 
     74      fis >> cat;
     75    }
     76    cout << cat ;
    8077    TArray<r_4> omap(NPhi,NTheta);
    81     double tet0 = Theta0Degre;
    82     double phi0 = Phi0Degre;
    83     double tetmax = tet0+ThetaSizeDegre;
    84     double phimax = phi0+PhiSizeDegre;
    85 
    86     cout << "srccat2cube[2]: projecting sources to map ..." << endl;
    87 
    88     sa_size_t srccnt=0;
    89     sa_size_t extendedsrccnt=0;
    90 
    91     double meanflx=0.;
    92     double flxmin=9.e99;
    93     double flxmax=-9.e99;
    94 
    95     double dtet = ThetaSizeDegre/(double)NTheta;
    96     double dphi = PhiSizeDegre/(double)NPhi;
    97     double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
    98 
    99     for (sa_size_t n=0; n<nvss.NRows(); n++)  {
    100       r_8* pline=nvss.GetLineD(n);
    101       double alpha=pline[idxa];  // alpha en degre
    102       double delta=pline[idxd];  // delta en degre
    103       double flx=pline[idxf]*1.e-3;  // flux en Jy
    104       double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.;  // taille (extension de la source en arcmin
    105       if (srcszarcmin<1.) srcszarcmin=1.;
    106       double tet = 90.-delta;
    107       double phi = alpha;
    108       sa_size_t i = (phi-phi0)/dphi;
    109       sa_size_t j = (tet-tet0)/dtet;
    110       if ((i<0)||(i>=omap.SizeX()))  continue;
    111       if ((j<0)||(j>=omap.SizeY()))  continue;
    112       double srat = (4.*srcszarcmin*srcszarcmin)/(mpixsizarcmin*mpixsizarcmin);
    113       if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
    114         omap(i,j) += flx*srat;
    115       }
    116       else {  // on repartit l'energie de la source dans plusieurs pixels
    117         extendedsrccnt++;
    118         for(int bi=-1;bi<=1;bi++) {
    119           for(sa_size_t bj=-1; bj<=1; bj++) {
    120             sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
    121             sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
    122             if ((ii<0)||(ii>=omap.SizeX()))  continue;
    123             if ((jj<0)||(jj>=omap.SizeY()))  continue;     
    124             if ((bi==0)&&(bj==0))  omap(ii,jj) += flx*srat*0.3;
    125             else omap(ii,jj) += flx*srat*0.7/8.;
    126           }
    127         }
    128       }
    129       srccnt++;   meanflx+=flx;
    130       if (flx<flxmin) flxmin=flx;
    131       if (flx>flxmax) flxmax=flx;
    132     }
    133 
    134     cout << "srccat2cube[3]: Output rectangular map computed " << endl;
    135     meanflx /= (double)srccnt;
    136     cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt
    137          << " -> meanFlx=" << meanflx << " min=" << flxmin
    138          << " max=" << flxmax << " Jy" << endl;
    139 
    140     double mean, sigma;
    141     r_4 mintemp, maxtemp;
    142     omap.MinMax(mintemp, maxtemp);
    143     MeanSigma(omap, mean, sigma);
    144     cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
    145     omap.Show();
    146 
    147     H21Conversions conv;
    148     conv.setRedshift(0.);
    149     conv.setOmegaPixDeg2(dphi*dtet);
    150     cout << "srccat2cube[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad" 
    151          << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
    152     omap *= (r_4)conv.toKelvin(1.);
    153     MeanSigma(omap, mean, sigma);
    154     cout << " After conversion  : Mean=" << mean << " Sigma=" << sigma << "  Kelvin " << endl;
     78    TArray<r_4> ocube(NPhi,NTheta,NFreq);
     79    if (copt=="-nvss")  nvssTocube(cat, omap, ocube);
     80    else north20Tocube(cat, omap, ocube);
    15581   
    156     if (narg > 3) {
    157       string ppfname = arg[3];
    158       cout << " srccat2cube[4]: Saving inmap/outmap tp PPF file-> " << ppfname << endl;
    159       POutPersist po(ppfname);
    160       po << PPFNameTag("omap") << omap;
    161     }
    162 
    163     TArray<r_4> ocube(NPhi,NTheta,NFreq);
    164 
    165     double infreq = 1420.; //  frequence de reference du flux des sources
    166     double freq0 = Freq0MHz;  // Freq0 du cube de sortie
    167     double dfreq = FreqSizeMHz/(double)NFreq; 
    168 
    169     ThSDR48RandGen rg;
    170     for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
    171       for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
    172         double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
    173         for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
    174           double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
    175           ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
    176         }
    177       }
    178     }
    179    
    180     // On sauve le cube de sortie
    181     {
    182       cout << " srccat2cube[5]: Saving output cube to -> " << outname << endl;
     82    {  // On sauve le cube de sortie
     83      cout << " srccat2cube[7]: Saving output cube to -> " << outname << endl;
    18384      POutPersist poc(outname);
    18485      poc << ocube;
    18586    }
    186 
     87   
     88    if (narg > 4) {
     89      string ppfname = arg[4];
     90      cout << " srccat2cube[8]: saving 2D source map to PPF file-> " << ppfname << endl;
     91      POutPersist po(ppfname);
     92      po << omap;
     93    }
    18794    rc = 0;
    18895  }
    … …  
    206113
    207114
     115
     116/* -- Fonction -- */
     117int nvssTocube(DataTable& nvss, TArray<r_4>& omap,  TArray<r_4>& ocube)
     118{
     119  sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
     120  sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
     121  sa_size_t idxf = nvss.IndexNom("S1_4");
     122  sa_size_t idxmajax = nvss.IndexNom("MajAxis");
     123  sa_size_t idxminax = nvss.IndexNom("MinAxis");
     124 
     125  cout << " NVSS catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf
     126       << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
     127 
     128  double tet0 = Theta0Degre;
     129  double phi0 = Phi0Degre;
     130  double tetmax = tet0+ThetaSizeDegre;
     131  double phimax = phi0+PhiSizeDegre;
     132 
     133  cout << "srccat2cube/NVSS[2]: projecting sources to map ..." << endl;
     134 
     135  sa_size_t srccnt=0;
     136  sa_size_t extendedsrccnt=0;
     137 
     138  double meanflx=0.;
     139  double flxmin=9.e99;
     140  double flxmax=-9.e99;
     141 
     142  double dtet = ThetaSizeDegre/(double)NTheta;
     143  double dphi = PhiSizeDegre/(double)NPhi;
     144  double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
     145 
     146  for (sa_size_t n=0; n<nvss.NRows(); n++)  {
     147    r_8* pline=nvss.GetLineD(n);
     148    double alpha=pline[idxa];  // alpha en degre
     149    double delta=pline[idxd];  // delta en degre
     150    double flx=pline[idxf]*1.e-3;  // flux en Jy
     151    double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.;  // taille (extension de la source en arcmin
     152    if (srcszarcmin<1.) srcszarcmin=1.;
     153    double tet = 90.-delta;
     154    double phi = alpha;
     155    sa_size_t i = (phi-phi0)/dphi;
     156    sa_size_t j = (tet-tet0)/dtet;
     157    if ((i<0)||(i>=omap.SizeX()))  continue;
     158    if ((j<0)||(j>=omap.SizeY()))  continue;
     159    if (srcszarcmin<(0.5*mpixsizarcmin)) {   // Toute l'energie dans un seul pixel
     160      omap(i,j) += flx;
     161    }
     162    else {  // on repartit l'energie de la source dans plusieurs pixels
     163      extendedsrccnt++;
     164      for(int bi=-1;bi<=1;bi++) {
     165        for(sa_size_t bj=-1; bj<=1; bj++) {
     166          sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
     167          sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
     168          if ((ii<0)||(ii>=omap.SizeX()))  continue;
     169          if ((jj<0)||(jj>=omap.SizeY()))  continue;       
     170          if ((bi==0)&&(bj==0))  omap(ii,jj) += flx*0.3;
     171          else omap(ii,jj) += flx*0.7/8.;
     172        }
     173      }
     174    }
     175    srccnt++;   meanflx+=flx;
     176    if (flx<flxmin) flxmin=flx;
     177    if (flx>flxmax) flxmax=flx;
     178  }
     179 
     180  cout << "srccat2cube/NVSS[3]: Output rectangular map computed " << endl;
     181  meanflx /= (double)srccnt;
     182  cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt
     183       << " -> meanFlx=" << meanflx << " min=" << flxmin
     184       << " max=" << flxmax << " Jy" << endl;
     185 
     186  double mean, sigma;
     187  r_4 minjy, maxjy;
     188  omap.MinMax(minjy, maxjy);
     189  MeanSigma(omap, mean, sigma);
     190  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky ;  Sizes:" << endl;
     191  omap.Show();
     192 
     193  H21Conversions conv;
     194  conv.setRedshift(0.);
     195  conv.setOmegaPixDeg2(dphi*dtet);
     196  cout << "srccat2cube/NVSS[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad" 
     197       << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
     198  omap *= (r_4)conv.toKelvin(1.);
     199  MeanSigma(omap, mean, sigma);
     200  r_4 minT, maxT;
     201  omap.MinMax(minT, maxT);
     202  cout << " NVSS/ After conversion  : Mean=" << mean << " Sigma=" << sigma << " Min=" << minT
     203       << " Max=" << maxT << "  Kelvin " << endl;
     204
     205  double infreq = 1420.; //  frequence de reference du flux des sources
     206  double freq0 = Freq0MHz;  // Freq0 du cube de sortie
     207  double dfreq = FreqSizeMHz/(double)NFreq; 
     208
     209  ThSDR48RandGen rg;
     210  for (sa_size_t j=0; j<ocube.SizeY(); j++)  {
     211    for (sa_size_t i=0; i<ocube.SizeX(); i++)  {
     212      double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
     213      for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
     214        double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
     215        ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
     216      }
     217    }
     218  }
     219  cout << "srccat2cube/NVSS[5] data cube created from sources " << endl;
     220  ocube.MinMax(minT, maxT);
     221  MeanSigma(ocube, mean, sigma);
     222  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
     223 
     224  return srccnt;
     225}
     226
     227/* -- Fonction -- */
     228int north20Tocube(DataTable& nor, TArray<r_4>& omap,  TArray<r_4>& ocube)
     229{
     230
     231  sa_size_t idxa = nor.IndexNom("ra");
     232  sa_size_t idxd = nor.IndexNom("dec");
     233  sa_size_t idxf = nor.IndexNom("flux");
     234  sa_size_t idxslo = nor.IndexNom("SpLO");
     235  sa_size_t idxshi = nor.IndexNom("SpHI");
     236 
     237  cout << " North20cm catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf
     238       << " SpLO: " << idxslo << " SpHI: " << idxshi << endl;
     239 
     240  double tet0 = Theta0Degre;
     241  double phi0 = Phi0Degre;
     242  double tetmax = tet0+ThetaSizeDegre;
     243  double phimax = phi0+PhiSizeDegre;
     244 
     245  cout << "srccat2cube/North20[2]: projecting sources to map ..." << endl;
     246 
     247  sa_size_t srccnt=0;
     248  sa_size_t lowoksrccnt=0;
     249 
     250  double meanflx=0.;
     251  double flxmin=9.e99;
     252  double flxmax=-9.e99;
     253 
     254  double dtet = ThetaSizeDegre/(double)NTheta;
     255  double dphi = PhiSizeDegre/(double)NPhi;
     256
     257  double infreq = 1420.; //  frequence de reference du flux des sources
     258  double freq0 = Freq0MHz;  // Freq0 du cube de sortie
     259  double dfreq = FreqSizeMHz/(double)NFreq; 
     260 
     261  for (sa_size_t n=0; n<nor.NRows(); n++)  {
     262    r_8* pline=nor.GetLineD(n);
     263    double alpha=pline[idxa];  // alpha en degre
     264    double delta=pline[idxd];  // delta en degre
     265    double flx=pline[idxf]*1.e-3;  // flux en Jy
     266    double tet = 90.-delta;
     267    double phi = alpha*360./24.;
     268    sa_size_t i = (phi-phi0)/dphi;
     269    sa_size_t j = (tet-tet0)/dtet;
     270    if ((i<0)||(i>=omap.SizeX()))  continue;
     271    if ((j<0)||(j>=omap.SizeY()))  continue;
     272    omap(i,j) += flx;
     273    srccnt++;   meanflx+=flx;
     274    if (flx<flxmin) flxmin=flx;
     275    if (flx>flxmax) flxmax=flx;
     276    double slo=pline[idxslo];
     277    if (slo<9.) {  // source detected at 80 cm
     278      lowoksrccnt++;
     279    }
     280    else slo=5.;
     281    for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
     282      double rapfreq = pow((freq0+k*dfreq)/infreq, slo);
     283      ocube(i,j,k) += flx*rapfreq;
     284    }
     285  }
     286 
     287  cout << "srccat2cube/North20[3]: Output rectangular map computed " << endl;
     288  meanflx /= (double)srccnt;
     289  cout << " SrcCount in map: " << srccnt << " SpLowOK=" << lowoksrccnt
     290       << " -> meanFlx=" << meanflx << " min=" << flxmin
     291       << " max=" << flxmax << " Jy" << endl;
     292 
     293  double mean, sigma;
     294  r_4 minjy, maxjy;
     295  omap.MinMax(minjy, maxjy);
     296  MeanSigma(omap, mean, sigma);
     297  cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
     298
     299  ocube.MinMax(minjy, maxjy);
     300  MeanSigma(ocube, mean, sigma);
     301  cout << " Cube : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
     302
     303  H21Conversions conv;
     304  conv.setRedshift(0.);
     305  conv.setOmegaPixDeg2(dphi*dtet);
     306  cout << "srccat2cube/North20[4] H21Conversions,  OmegaPix=" << conv.getOmegaPix() << " srad" 
     307       << " @1400MHz  toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
     308 
     309  // Jansky to Kelvin conversion
     310  for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
     311    conv.setFrequency(freq0+k*dfreq);
     312    //    cout << " DBG* Freq= " << freq0+k*dfreq << " -> toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
     313    ocube(Range::all(), Range::all(), Range(k)) *= (r_4)conv.toKelvin(1.);
     314  }
     315  cout << "srccat2cube/North20[5] data cube in Kelvin computed " << endl;
     316
     317  r_4 minT, maxT;
     318  ocube.MinMax(minT, maxT);
     319  MeanSigma(ocube, mean, sigma);
     320  cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
     321 
     322  return srccnt;
     323}

  • trunk/Cosmo/RadioBeam/subtractradsrc.cmd

    r3825 r3829  
    4646csh> cat > sumcubes.pic
    4747openppf syncube.ppf
    48 openppf nvsscube.ppf
     48openppf srcnor3d.ppf
    4949# expmeansig syncube val
    50 # expmeansig nvsscube val
    51 c++exec TArray<r_4> fgndcube = syncube+nvsscube; KeepObj(fgndcube);
     50# expmeansig srcnor3d val
     51c++exec TArray<r_4> fgndcube = syncube+srcnor3d; KeepObj(fgndcube);
    5252print fgndcube
    5353# expmeansig fgndcube val
    … …  
    8989# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130
    9090csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3.
     91# 4.f / Estimate residual noise from Foreground removal :
     92csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3.
     93csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.
    9194
    9295### Step 5 / Check the results using spiapp
    … …  
    103106
    104107openppf subpknolss.ppf
     108openppf subpknolssnocor.ppf
    105109
    106110disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-11,8e-6 gold'
    107111disp lsspklobe 'same nsta orange'
    108112disp lsspklobewn 'same nsta siennared'
     113settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black'
     114
     115#  Calcul du volume total en Mpc^3
     116set VOL 3*3*3*360*360*256
     117plot2d lsspklobe x val*$VOL 1 'same nsta orange'
     118plot2d lsspklobe x val*$VOL 1 'same nsta orange'
    109119
    110120disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
    … …  
    121131disp lsspklobewn 'same nsta siennared'
    122132disp subpkcorlobe 'same nsta red'
    123 disp subpknolsslobe 'same nsta green'
     133disp subpknolss 'same nsta green'
    124134
    125 settitle 'Recovered Pk[LSS] from LSS+Foreground(GSM) (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
    126 set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[residual,NoLSS]' )
    127 set cols ( gold siennared red green )
     135#  Calcul du volume total en Mpc^3
     136set VOL 3*3*3*360*360*256
     137plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
     138plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared'
     139plot2d subpkcorlobe x val*$VOL 1 'same nsta cpts marker=box,5 red'
     140plot2d subpklobe x val*$VOL 1 'same nsta cpts marker=box,5 blueviolet'
     141plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green'
     142
     143settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
     144setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16'
     145set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[ExtLSS,NoBeamCor]' 'Pk[residual,NoLSS]' )
     146set cols ( gold siennared red blueviolet green )
    128147textdrawer lines cols 'font=helvetica,bold,16 frame'
     148
     149
     150plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
     151plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 red'
     152plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green'
     153plot2d subpknolssnocor x val*$VOL 1 'same nsta cpts marker=box,5 magenta'
     154setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16'
     155settitle 'Recovered Pk[LSS] and residual systematics' ' ' 'font=helvetica,bold,18'
     156set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[residual,NoLSS]' 'Pk[residual,NoLSS,NoBeamCorrection]' )
     157set cols ( gold red green magenta )
     158textdrawer lines cols 'font=helvetica,bold,16 frame'

  • trunk/Cosmo/RadioBeam/tjyk.cc

    r3787 r3829  
    2525       << " MHz , Lambda=" << conv.getLambda() << " m , OmegaPix=" << conv.getOmegaPix() << " srad" << endl;
    2626  cout << "  toKelvin(1 Jy)= " << conv.toKelvin(1.) << "    toJansky(1 K)=" << conv.toJansky(1.) << endl;
     27  cout << "  Mean 21cm emission temperature= " << conv.T21cm_mK() << " milliKelvin " << endl;
    2728  cout << " --------------------------------------------------------------- " << endl;
    2829
Note: See TracChangeset for help on using the changeset viewer.