Changeset 3825 in Sophya


Ignore:
Timestamp:
Aug 2, 2010, 7:29:48 PM (15 years ago)
Author:
ansari
Message:

modifs et ajout de programme pour traitement cartes GSM (Global Sky Model), Reza 02/08/2010

Location:
trunk/Cosmo/RadioBeam
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/README

    r3793 r3825  
    2222  - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
    2323  - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
     24  - gsm2cube.cc : Computes 3D temperature cube from a set of GSM (Global Sky Model) healpix PPF maps 
    2425  - tjyk.cc : test de la classe H21Conversions
    2526

  • trunk/Cosmo/RadioBeam/calcpk.cc

    r3789 r3825  
    5353    bool fgaddnoise=false;
    5454    if (narg>4) {
    55       pixsignoise=atof(arg[3]);
     55      pixsignoise=atof(arg[4]);
    5656      fgaddnoise=true;
    5757    }

  • trunk/Cosmo/RadioBeam/cubedef.h

    r3796 r3825  
    1616static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
    1717
     18/* ---  Parametres avec generation Reza from HASLAM + NVSS */
    1819static double Freq0MHz = 875.;       // premiere frequence
    1920static double FreqSizeMHz = 70.4;    // Taille de la carte en MHz (axe Z)
    2021
    21 /* z = 0.56, freq_center=910 MHz , dA_comov ~= 2060 Mpc , 
    22        5 arcmin -> 3 Mpc , 0.275 MHz -> 3 Mpc  (H(ze) ~ 94 km/s/Mpc */
     22// z = 0.56, freq_center=910 MHz , dA_comov ~= 2060 Mpc , 
     23//       5 arcmin -> 3 Mpc , 0.500 MHz -> 3 Mpc  (H(ze) ~ 94 km/s/Mpc)
    2324// Taille de cellule pour le calcul du spectre de puissance 3D
    2425static double ComovDA = 2060.;
    … …  
    2627static double YCellSizeMpc = 3.;
    2728static double ZCellSizeMpc = 1.5;
     29
     30/* --- Parametres pour utilisation des cartes GSM 
     31static double Freq0MHz = 820.;       // premiere frequence
     32static double FreqSizeMHz = 128.;    // Taille de la carte en MHz (axe Z)
     33
     34// z = 0.6, freq_center=884 MHz , dA_comov ~= 2185 Mpc , 
     35//       5 arcmin -> 3 Mpc , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
     36// Taille de cellule pour le calcul du spectre de puissance 3D
     37static double ComovDA = 2185.;
     38static double XCellSizeMpc = 3.;
     39static double YCellSizeMpc = 3.;
     40static double ZCellSizeMpc = 3.;
     41*/
     42
    2843
    2944//---  Parametres des lois de puissance en frequence

  • trunk/Cosmo/RadioBeam/makefile

    r3792 r3825  
    66include $(SOPHYABASE)/include/sophyamake.inc
    77
    8 all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe
     8all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube
    99
    1010clean :
    … …  
    2929syncube : Objs/syncube
    3030        echo 'makefile : syncube made'
     31
     32gsm2cube : Objs/gsm2cube
     33        echo 'makefile : gsm2cube made'
    3134
    3235srcat2cube : Objs/srcat2cube
    … …  
    9598        $(CXXCOMPILE) -o Objs/applobe.o applobe.cc
    9699
     100### programme gsm2cube
     101Objs/gsm2cube : Objs/gsm2cube.o $(PKGOLIST)
     102        $(CXXLINK) -o Objs/gsm2cube Objs/gsm2cube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
     103
     104Objs/gsm2cube.o : gsm2cube.cc  $(PKGHLIST)
     105        $(CXXCOMPILE) -o Objs/gsm2cube.o gsm2cube.cc
    97106
    98107### les classes / fonctions

  • trunk/Cosmo/RadioBeam/specpk.cc

    r3787 r3825  
    4646{
    4747  typ_=typ;
     48  SetRenormFac();
    4849}
    4950
    … …  
    5253{
    5354  wk/=DeuxPI;
     55  double retv=1.;
    5456  switch (typ_)
    5557    {
    5658    case 1:
    57       return Pnu1(wk);
     59      retv=Pnu1(wk);
    5860      break;
    5961    case 2:
    60       return Pnu2(wk);
     62      retv=Pnu2(wk);
    6163      break;
    6264    case 3:
    63       return Pnu3(wk);
     65      retv=Pnu3(wk);
    6466      break;
    6567    case 4:
    66       return Pnu4(wk);
     68      retv=Pnu4(wk);
    6769      break;
    6870    default :
    … …  
    8385        }
    8486      }
    85       return csp;
    86       }
    87       break;
    88     }
     87      retv=csp;
     88      }
     89      break;
     90    }
     91  return retv*renorm_fac;
    8992}
    9093// Return a vector representing the power spectrum (for checking)
    … …  
    97100}
    98101
     102double SpectralShape::Sommek2Pk(double kmax, int n)
     103{
     104  double dk=kmax/(double)n;
     105  double s=0.;
     106  for(int i=1; i<=n; i++) {
     107    double ck=(double)i*dk;
     108    s += Value(ck)*ck*ck;
     109  }
     110  return s*dk*4.*M_PI;
     111}
    99112//--------------------------------------------------
    100113// -- Four2DResponse class : test P(k) class

  • trunk/Cosmo/RadioBeam/specpk.h

    r3769 r3825  
    3030  inline  double Value(double wk) { return((*this)(wk)); }
    3131// Return a vector representing the power spectrum (for checking)
    32   Histo GetPk(int n=256);       
     32  Histo GetPk(int n=256);
     33  double Sommek2Pk(double kmax=1000., int n=5000);
     34  inline void SetRenormFac(double f=1.) { renorm_fac=f; }
    3335  int typ_;
     36  double renorm_fac;
    3437};
    3538

  • trunk/Cosmo/RadioBeam/srcat2cube.cc

    r3785 r3825  
    7272    sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
    7373    sa_size_t idxf = nvss.IndexNom("S1_4");
    74     cout << " ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf << endl;
     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;
    7579
    7680    TArray<r_4> omap(NPhi,NTheta);
    … …  
    8387
    8488    sa_size_t srccnt=0;
     89    sa_size_t extendedsrccnt=0;
     90
    8591    double meanflx=0.;
    8692    double flxmin=9.e99;
    … …  
    8995    double dtet = ThetaSizeDegre/(double)NTheta;
    9096    double dphi = PhiSizeDegre/(double)NPhi;
     97    double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
     98
    9199    for (sa_size_t n=0; n<nvss.NRows(); n++)  {
    92100      r_8* pline=nvss.GetLineD(n);
    … …  
    94102      double delta=pline[idxd];  // delta en degre
    95103      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.;
    96106      double tet = 90.-delta;
    97107      double phi = alpha;
    … …  
    100110      if ((i<0)||(i>=omap.SizeX()))  continue;
    101111      if ((j<0)||(j>=omap.SizeY()))  continue;
    102       omap(i,j) += flx;
     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      }
    103129      srccnt++;   meanflx+=flx;
    104130      if (flx<flxmin) flxmin=flx;
    … …  
    108134    cout << "srccat2cube[3]: Output rectangular map computed " << endl;
    109135    meanflx /= (double)srccnt;
    110     cout << " SrcCount in map: " << srccnt << " -> meanFlx=" << meanflx << " min=" << flxmin
     136    cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt
     137         << " -> meanFlx=" << meanflx << " min=" << flxmin
    111138         << " max=" << flxmax << " Jy" << endl;
    112139
    113140    double mean, sigma;
     141    r_4 mintemp, maxtemp;
     142    omap.MinMax(mintemp, maxtemp);
    114143    MeanSigma(omap, mean, sigma);
    115144    cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;

  • trunk/Cosmo/RadioBeam/subtractradsrc.cmd

    r3796 r3825  
    99# 1.a/ Run SimLSS
    1010csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,1.5 -Z 0.56 -8 1. -n 10000 -O 0,2 -o lssz056 -T 2
    11 # 1.b/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention
     11# 1.b/ To run SimLSS with GSM map parameters (DeltaFreq=500 MHz)
     12csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 360,3 -y 360,3 -z 256,3 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
     13
     14# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention
    1215#  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray)
    1316#  RadioBeam cubes : the radial (redshift) direction along Z axis of the cube (TArray)
    … …  
    2730  KeepObj(omap);
    2831
    29 print omap
    30 saveppf omap lsscube.ppf
     32rename omap lsscube
     33print lsscube
     34# expmeansig lsscube val
     35saveppf lsscube lsscube.ppf
    3136
    3237csh> spiapp -term -exec racube.pic
     38
    3339
    3440## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin)
    … …  
    6975#  with the lobe effect
    7076csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2 
     77csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.2 0.35 
    7178csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2 
    7279
    … …  
    7885# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam
    7986csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0.
    80 # 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130
    81 csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 130. 3.
     87# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction 
     88csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3.
     89# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130
     90csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3.
    8291
    8392### Step 5 / Check the results using spiapp
     93setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
    8494delobjs *
    8595openppf fgndpk.ppf
    86 disp fgndpk 'logx logy'
    87 # Adjust the limits
    8896openppf fgndpklobe.ppf
    89 disp fgndpklobe 'same red'
    90 disp fgndpklobe 'same red nsta'
    9197openppf fgndpkcorlobe.ppf
    92 disp fgndpkcorlobe 'same magenta nsta'
    9398openppf lsspk.ppf
    94 disp lsspk 'same gold nsta'
    9599openppf lsspklobe.ppf
    96 disp lsspklobe 'same green nsta'
     100openppf lsspklobewn.ppf
    97101openppf subpklobe.ppf
    98 disp subpklobe 'same blue nsta'
     102openppf subpkcorlobe.ppf
     103
     104openppf subpknolss.ppf
     105
     106disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-11,8e-6 gold'
     107disp lsspklobe 'same nsta orange'
     108disp lsspklobewn 'same nsta siennared'
     109
     110disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
     111disp fgndpklobe 'same nsta blue'
     112disp fgndpkcorlobe 'same nsta skyblue'
     113disp lsspk 'same nsta gold'
     114disp lsspklobewn 'same nsta siennared'
     115settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
     116set lines ( 'Pk[Foreground]'  'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
     117set cols ( navyblue blue skyblue gold siennared )
     118textdrawer lines cols 'font=helvetica,bold,16 frame'
     119
     120disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
     121disp lsspklobewn 'same nsta siennared'
     122disp subpkcorlobe 'same nsta red'
     123disp subpknolsslobe 'same nsta green'
     124
     125settitle 'Recovered Pk[LSS] from LSS+Foreground(GSM) (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
     126set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[residual,NoLSS]' )
     127set cols ( gold siennared red green )
     128textdrawer lines cols 'font=helvetica,bold,16 frame'
Note: See TracChangeset for help on using the changeset viewer.