Changeset 3787 in Sophya


Ignore:
Timestamp:
Jun 25, 2010, 12:00:30 PM (15 years ago)
Author:
ansari
Message:

Ajout classes / programmes de calcul d'effet de lobe sur les radio-sources, Reza 25/06/2010

Location:
trunk/Cosmo/RadioBeam
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/README

    r3785 r3787  
    1616  - plpkn.pic : spiapp script for plotting pknoise.cc program output
    1717A.2/ power spectrum computation and foreground/radio source
     18  - subtractradsrc.cmd : Explanation and commands to compute various power spectra and foreground subtraction
    1819  - calcpk.cc : power spectrum P(k) computation from 3D mass or temperature map
    1920  - calcpk2.cc : LSS deltaT/T P(k) computation from LSS+foreground cubes after cleaning
     21  - applobe.cc : Antenna/array beam effect on 3D sky cube
    2022  - syncube.cc : Synchrotron 3D temperature cube from HASLAM 400 MHz map
    2123  - srcat2cube.cc : radio source 3D temperature cube from NVSS catalog
    … …  
    2426A.3/ List of other common files :
    2527  - specpk.h specpk.cc : 3D power spectrum computation classes 
     28  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
    2629  - mdish.h mdish.cc : Classes for representing Multi-Dish interferometer configurations
    2730  - radutil.h radutil.cc : utilitaire de conversion a 21 cm

  • trunk/Cosmo/RadioBeam/calcpk.cc

    r3783 r3787  
    1919
    2020#include "qhist.h"
     21#include "lobe.h"
     22
    2123#include "histinit.h"
    2224#include "fftwserver.h"
    … …  
    3335{
    3436  if (narg<3) {
    35     cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] " << endl;
     37    cout << " Usage: calcpk In3DMap OutPk  [InRenormFactor] [PixNoiseLevel] " << endl;
    3638    return 1;
    3739  }
    … …  
    4143    string inppfname = arg[1];
    4244    string outname = arg[2];
    43     r_4 renfact=1.;
     45    double renfact=1.;
    4446    bool fgrenfact=false;
    4547    if (narg>3) {
    … …  
    4749      fgrenfact=true;
    4850    }
     51    double pixsignoise = 0.;
     52    bool fgaddnoise=false;
     53    if (narg>4) {
     54      pixsignoise=atof(arg[3]);
     55      fgaddnoise=true;
     56    }
     57
    4958    cout << "calcpk[1] : reading 3D map from file " << inppfname << endl;
    5059    TArray<r_4> inmap;
    … …  
    5362      pin >> inmap;
    5463    }
    55     if (fgrenfact) inmap *= renfact;
     64    if (fgrenfact) {
     65      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
     66      inmap *= renfact;
     67    }
    5668    double mean, sigma;
    5769    MeanSigma(inmap, mean, sigma);
    … …  
    5971    inmap.Show();
    6072    cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
     73    if (fgaddnoise) {
     74      BeamEffect::AddNoise(inmap, pixsignoise);
     75    }
     76
    6177    tm.Split(" After read ");
    6278   

  • trunk/Cosmo/RadioBeam/calcpk2.cc

    r3784 r3787  
    77    R. Ansari , C. Magneville - Juin 2010
    88
    9 Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk
     9Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPk [PixNoiseLevel]
    1010---------------------------------------------------------------  */
    1111
    … …  
    2222
    2323#include "qhist.h"
     24#include "lobe.h"
     25#include "cubedef.h"
     26
    2427#include "histinit.h"
    2528#include "fftwserver.h"
    … …  
    3033typedef ThSDR48RandGen RandomGenerator ;
    3134//--- Declaration des fonctions
    32 TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx);
     35TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq,
     36                            TArray<r_8>& synctemp, TArray<r_8>& specidx);
    3337
    3438//-------------------------------------------------------------------------
    … …  
    3943{
    4044  if (narg<6) {
    41     cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk " << endl;
     45    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync  OutPk [PixNoiseLevel] " << endl;
    4246    return 1;
    4347  }
    … …  
    5054    r_4 rfacsync = atof(arg[4]);
    5155    string outname = arg[5];
     56
     57    double pixsignoise = 0.;
     58    bool fgaddnoise=false;
     59    if (narg>6) {
     60      pixsignoise=atof(arg[6]);
     61      fgaddnoise=true;
     62    }
     63
     64
    5265    TArray<r_4> inmaplss, inmapsync;
    5366    const char * tits[2]={"LSS", "Sync"};
    … …  
    6881      cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
    6982    }
     83
     84    if (fgaddnoise) {
     85      cout << " calcpk2: adding noise to LSS input cube ... " << endl;
     86      BeamEffect::AddNoise(inmaplss, pixsignoise);
     87    }
     88
    7089    tm.Split(" After read ");
    7190    TArray<r_8> synctemp, specidx;
    7291    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
    73     TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, synctemp, specidx);
     92    double freq0 = Freq0MHz;
     93    double dfreq = FreqSizeMHz/(double)NFreq;
     94    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, freq0, dfreq, synctemp, specidx);
    7495    double mean, sigma;
    7596    MeanSigma(inmap, mean, sigma);
    … …  
    117138
    118139/* --Fonction-- */
    119 TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx)
     140TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, double freq0, double dfreq,
     141                            TArray<r_8>& synctemp, TArray<r_8>& specidx)
     142// Inputs : maplss, mapsyc, freq0, dfreq
     143// Outputs : synctemp, specidx  (reconstructed foreground temperature and spectral index
     144// Return_Array : foreground subtracted LSS signal
    120145{
    121146  bool smo;
    … …  
    131156  Vector vlnf(maplss.SizeZ());
    132157  int nprt = 0;
    133   double freq0=840.;  // Frequence premier index en k (MHz)
    134   double dfreq=1.;   // largeur en frequence de chaque plan (Mhz) 
     158  // double freq0 : Frequence premier index en k (MHz)
     159  // double dfreq :   // largeur en frequence de chaque plan (Mhz) 
    135160  for(sa_size_t i=0; i<maplss.SizeX(); i++)
    136161    for(sa_size_t j=0; j<maplss.SizeY(); j++)  {

  • trunk/Cosmo/RadioBeam/cubedef.h

    r3785 r3787  
    66
    77// Definition tailles et bornes du cube angles (X,Y), frequence (Z)
    8 static sa_size_t NTheta=256;
    9 static sa_size_t NPhi=256;
    10 static sa_size_t NFreq = 128;
    118
    12 static double Theta0Degre = 60.;  // -> Delta = +30 deg
    13 static double Phi0Degre = 120.;   // -> alpha = 8h00
    14 static double ThetaSizeDegre = 60.;  // Taille de la carte en degre selon alpha (axe X)
    15 static double PhiSizeDegre = 60.;    // Taille de la carte en degre selon delta (axe Y)
     9static sa_size_t NTheta=360;
     10static sa_size_t NPhi=360;
     11static sa_size_t NFreq = 200;
     12
     13static double Theta0Degre = 75.;  // -> Delta = +30 deg
     14static double Phi0Degre = 135.;   // -> alpha = 9h00
     15static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
     16static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
    1617
    1718static double Freq0MHz = 840.;       // premiere frequence
    18 static double FreqSizeMHz = 128.;    // Taille de la carte en MHz (axe Z)
     19static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
    1920
    2021
    … …  
    3132static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index 
    3233
     34/*
     35static sa_size_t NTheta=256;
     36static sa_size_t NPhi=256;
     37static sa_size_t NFreq = 128;
     38
     39static double Theta0Degre = 60.;  // -> Delta = +30 deg
     40static double Phi0Degre = 120.;   // -> alpha = 8h00
     41static double ThetaSizeDegre = 60.;  // Taille de la carte en degre selon alpha (axe X)
     42static double PhiSizeDegre = 60.;    // Taille de la carte en degre selon delta (axe Y)
     43static double Freq0MHz = 840.;       // premiere frequence
     44static double FreqSizeMHz = 100.;    // Taille de la carte en MHz (axe Z)
     45*/
     46
    3347#endif

  • trunk/Cosmo/RadioBeam/makefile

    r3785 r3787  
    66include $(SOPHYABASE)/include/sophyamake.inc
    77
    8 all : pknoise calcpk calcpk2 syncube srcat2cube tjyk
     8all : pknoise calcpk calcpk2 syncube srcat2cube tjyk applobe
    99
    1010clean :
    1111        rm Objs/*
     12
     13PKGOLIST =  Objs/lobe.o Objs/specpk.o Objs/mdish.o Objs/qhist.o Objs/radutil.o
     14PKGHLIST =  lobe.h specpk.h mdish.h qhist.h radutil.h cubedef.h
    1215
    1316### les executables
    … …  
    3033        echo 'makefile : tjyk made'
    3134
     35applobe : Objs/applobe
     36        echo 'makefile : applobe made'
     37
    3238### programme pknoise
    33 Objs/pknoise : Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o
    34         $(CXXLINK) -o Objs/pknoise Objs/pknoise.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
     39Objs/pknoise : Objs/pknoise.o $(PKGOLIST)
     40        $(CXXLINK) -o Objs/pknoise Objs/pknoise.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
    3541
    36 Objs/pknoise.o : pknoise.cc specpk.h  mdish.h qhist.h
     42Objs/pknoise.o : pknoise.cc $(PKGHLIST)
    3743        $(CXXCOMPILE) -o Objs/pknoise.o pknoise.cc
    3844
    3945### programme calcpk
    40 Objs/calcpk : Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o
    41         $(CXXLINK) -o Objs/calcpk Objs/calcpk.o Objs/specpk.o Objs/mdish.o Objs/qhist.o  $(SOPHYAEXTSLBLIST)
     46Objs/calcpk : Objs/calcpk.o $(PKGOLIST)
     47        $(CXXLINK) -o Objs/calcpk Objs/calcpk.o $(PKGOLIST)  $(SOPHYAEXTSLBLIST)
    4248
    43 Objs/calcpk.o : calcpk.cc specpk.h  mdish.h qhist.h
     49Objs/calcpk.o : calcpk.cc $(PKGHLIST)
    4450        $(CXXCOMPILE) -o Objs/calcpk.o calcpk.cc
    4551
    4652### programme calcpk2
    47 Objs/calcpk2 : Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o
    48         $(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o Objs/specpk.o Objs/mdish.o Objs/qhist.o $(SOPHYAEXTSLBLIST)
     53Objs/calcpk2 : Objs/calcpk2.o $(PKGOLIST)
     54        $(CXXLINK) -o Objs/calcpk2 Objs/calcpk2.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
    4955
    50 Objs/calcpk2.o : calcpk2.cc specpk.h  mdish.h qhist.h
     56Objs/calcpk2.o : calcpk2.cc $(PKGHLIST)
    5157        $(CXXCOMPILE) -o Objs/calcpk2.o calcpk2.cc
    5258
    5359### programme syncube
    54 Objs/syncube : Objs/syncube.o
    55         $(CXXLINK) -o Objs/syncube Objs/syncube.o $(SOPHYAEXTSLBLIST)
     60Objs/syncube : Objs/syncube.o $(PKGOLIST)
     61        $(CXXLINK) -o Objs/syncube Objs/syncube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
    5662
    5763Objs/syncube.o : syncube.cc
    … …  
    5965
    6066### programme srcat2cube
    61 Objs/srcat2cube : Objs/srcat2cube.o
    62         $(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o  Objs/radutil.o $(SOPHYAEXTSLBLIST)
     67Objs/srcat2cube : Objs/srcat2cube.o $(PKGOLIST)
     68        $(CXXLINK) -o Objs/srcat2cube Objs/srcat2cube.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
    6369
    64 Objs/srcat2cube.o : srcat2cube.cc radutil.h
     70Objs/srcat2cube.o : srcat2cube.cc $(PKGHLIST)
    6571        $(CXXCOMPILE) -o Objs/srcat2cube.o srcat2cube.cc
    6672
    6773### programme tjyk
    68 Objs/tjyk : Objs/tjyk.o Objs/radutil.o
    69         $(CXXLINK) -o Objs/tjyk Objs/tjyk.o Objs/radutil.o $(SOPHYAEXTSLBLIST)
     74Objs/tjyk : Objs/tjyk.o $(PKGOLIST)
     75        $(CXXLINK) -o Objs/tjyk Objs/tjyk.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
    7076
    71 Objs/tjyk.o : tjyk.cc radutil.h
     77Objs/tjyk.o : tjyk.cc $(PKGHLIST)
    7278        $(CXXCOMPILE) -o Objs/tjyk.o tjyk.cc
     79
     80### programme applobe
     81Objs/applobe : Objs/applobe.o $(PKGOLIST)
     82        $(CXXLINK) -o Objs/applobe Objs/applobe.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
     83
     84Objs/applobe.o : applobe.cc $(PKGHLIST)
     85        $(CXXCOMPILE) -o Objs/applobe.o applobe.cc
     86
    7387
    7488### les classes / fonctions
    … …  
    7993        $(CXXCOMPILE) -o Objs/mdish.o mdish.cc
    8094
     95Objs/lobe.o : lobe.cc lobe.h mdish.h qhist.h
     96        $(CXXCOMPILE) -o Objs/lobe.o lobe.cc
     97
    8198Objs/qhist.o : qhist.cc qhist.h
    8299        $(CXXCOMPILE) -o Objs/qhist.o qhist.cc

  • trunk/Cosmo/RadioBeam/mdish.cc

    r3783 r3787  
    1212  : typ_(typ), dx_((dx>1.e-3)?dx:1.), dy_((dy>1.e-3)?dy:1.)
    1313{
     14  setLambdaRef();
     15  setLambda();
    1416}
    1517
    … …  
    1719double Four2DResponse::Value(double kx, double ky)
    1820{
     21  kx *= lambda_ratio_;
     22  ky *= lambda_ratio_;
    1923  double wk,wkx,wky;
    2024  switch (typ_)
    … …  
    6266double Four2DRespTable::Value(double kx, double ky)
    6367{
     68  kx *= lambda_ratio_;
     69  ky *= lambda_ratio_;
    6470  int_4 i,j;
    6571  if ( (kx<=hrep_.XMin())||(kx>=hrep_.XMax()) ||

  • trunk/Cosmo/RadioBeam/mdish.h

    r3783 r3787  
    3434  { typ_ = a.typ_; dx_=a.dx_; dy_=a.dy_; return (*this); }
    3535
     36  inline void setLambdaRef(double lambda=1.)
     37  { lambdaref_ = lambda; }
     38  inline void setLambda(double lambda=1.)
     39  { lambda_ = lambda;  lambda_ratio_ = lambdaref_/lambda_; }
     40 
    3641  // Return the 2D response for wave vector (kx,ky)
    3742  virtual double Value(double kx, double ky);
    … …  
    4550  int typ_;
    4651  double dx_, dy_;
     52  double lambdaref_, lambda_;
     53  double lambda_ratio_;     // lambdaref_/lambda_;
     54
    4755};
    4856

  • trunk/Cosmo/RadioBeam/specpk.cc

    r3783 r3787  
    151151  // fourAmp represent 3-D fourier transform of a real input array.
    152152  // The second half of the array along Y and Z contain negative frequencies
    153   double kxx, kyy, kzz;
     153  double kxx, kyy, kzz, rep, amp;
    154154  // sa_size_t is large integer type 
    155155  for(sa_size_t kz=0; kz<fourAmp.SizeZ(); kz++) {
    … …  
    158158      kyy =  (ky>fourAmp.SizeY()/2) ? -(double)(fourAmp.SizeY()-ky)*dky_ : (double)ky*dky_;
    159159      for(sa_size_t kx=0; kx<fourAmp.SizeX(); kx++) {
    160         double kxx=(double)kx*dkx_;
    161         double rep = resp(kxx, kyy);
     160        kxx=(double)kx*dkx_;
     161        rep = resp(kxx, kyy);
    162162        if (crmask&&(kz==0))  mask(ky,kx)=((rep<1.e-8)?9.e9:(1./rep));
    163163        if (rep<1.e-8)  fourAmp(kx, ky, kz) = complex<TF>(9.e9,0.);
    164164        else {
    165           double amp = 1./sqrt(rep)/sqrt(2.);
     165          amp = 1./sqrt(rep)/sqrt(2.);
    166166          fourAmp(kx, ky, kz) = complex<TF>(rg_.Gaussian(amp), rg_.Gaussian(amp));   
    167167        }

  • trunk/Cosmo/RadioBeam/syncube.cc

    r3785 r3787  
    115115    double mean, sigma;
    116116    MeanSigma(omap, mean, sigma);
    117     cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << "  Jy - Sizes:" << endl;
     117    cout << "syncube[2] omap : Mean=" << mean << " Sigma=" << sigma << " Sizes:" << endl;
    118118    omap.Show();
    119119

  • trunk/Cosmo/RadioBeam/tjyk.cc

    r3785 r3787  
    66int main(int narg, char* arg[])
    77{
     8  if (narg<2) {
     9    cout << " Usage: tjyk redshift [resol_arcmin] " << endl;
     10    return 1;
     11  }
    812
    913  H21Conversions conv;
Note: See TracChangeset for help on using the changeset viewer.