Changeset 3973 in Sophya for trunk/Cosmo


Ignore:
Timestamp:
Apr 18, 2011, 5:30:44 PM (14 years ago)
Author:
ansari
Message:

Corrections diverses: choix lobe gaussien/triangle et specif DishDiameter au lieu de DoL ds applobe/calcpk2, possibilite application lobe freq.independante ds applobe, Reza 18/04/2011

Location:
trunk/Cosmo/RadioBeam
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/README

    r3829 r3973  
    2626
    2727A.3/ List of other common files :
    28   - fgndsub.h fgndsub.h : radio source, synchrotron subtraction class
     28  - fgndsub.h fgndsub.cc : radio source, synchrotron subtraction class
    2929  - lobe.h lobe.cc : 2D (k_alpha, k_delta) beam effect on 3D sky cubes
    3030  - specpk.h specpk.cc : 3D power spectrum computation classes 

  • trunk/Cosmo/RadioBeam/applobe.cc

    r3797 r3973  
    33    R. Ansari , C. Magneville - Juin 2010
    44
    5   Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName
    6                  [TargetBeamDoL] [ResmapleFactor=0.5,0.333...]
     5  Usage: applobe [-g -t -fib]  Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName
     6                 [TargetBeamDiam] [ResmapleFactor=0.5,0.333...]
     7   o -t : triangle beam shape in k-space
     8   o -g : gaussian beam shape in k-space
     9   o -fib: application d'un lobe fixe (independant de la frequence)
     10   o Diameter/Four2DRespTableFile : diametre dish ou nom de fichier PPF reponse 2D
     11   o In3DPPFName Out3DPPFName : Noms des fichiers (TArray<r_4> 3D) entree-sortie
     12   o TargetBeamDiam : Correction à un beam fiduciel, independant de la frequence
     13     avec specification de la valeur du diametre pour la frequence la plus basse
     14   o ResmapleFactor : reechantillonnage du cube (2 -> surechantillonnage, 0.5 : 1 sur 2)
    715---------------------------------------------------------------  */
    816
    … …  
    4856
    4957  InitTim();   // Initializing the CPU timer
     58
     59  if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
     60    cout << "Usage: applobe [-t -g -fib] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n"
     61         << "               [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
     62    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
     63      cout << "   -t -g : Triangular / gaussian beam shape (def=gaussian) \n"
     64           << "   -fib :  Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n"
     65           << "   Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n"
     66           << "   In3DPPFName Out3DPPFName: Input/Output PPF file names  (TArray<r_4> 3D) \n"
     67           << "   TargetBeamDiam: Corrected beam target diameter (D/Lambda at lowest frequency) \n"
     68           << "      DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm)) \n"
     69           << "   ResmapleFactor : Resampling (2-> oversampling, 0.5 : 1/2 undersampling) \n" << endl;   
     70    } 
     71    return 1;
     72  }
     73
    5074  Timer tm("applobe");
    5175
     76  // decodage arguments
     77  bool fggaussian=true;  // true -> gaussian beam
     78  bool fixedbeam=false;  // true -> apply freq. independent beam
     79
     80  // decodage argument optionnel
     81  bool fgoptarg=true;
     82  while (fgoptarg) {
     83    string fbo = arg[1];
     84    if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
     85    else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
     86    else if (fbo=="-fib")  { fixedbeam=true; arg++; narg--; }
     87    else fgoptarg=false;
     88  }
    5289  if (narg < 3) {
    53     cout << "Usage: applobe Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName "
    54          << "             [TargetBeamDoL] [ResmapleFactor=0.5,0.333...] \n" << endl;
    55     return 1;
    56   }
    57 
    58   // decodage arguments
     90    cout << " applobe/error arguments , applobe -h for help " << endl;
     91    return 2;
     92  }
     93
    5994  bool fgresptbl=true;
    6095  double DIAMETRE=100.;
    … …  
    6499    DIAMETRE=atof(arg[1]);
    65100  }
    66   else resptblname=arg[1];
    67 
     101  else  resptblname=arg[1];
     102 
    68103  string inname = arg[2];
    69104  string outname = arg[3];
    70105  bool fgcorrbeam=false;
    71   double tbeamDoL=135;
     106  double tbeamDiam=50.;
    72107  if (narg>4) {
    73     tbeamDoL=atof(arg[4]);
    74     if (tbeamDoL>1)  fgcorrbeam=true;
     108    tbeamDiam=atof(arg[4]);
     109    if (tbeamDiam>1.)  fgcorrbeam=true;
    75110  }
    76111  bool fgresample=false;
    … …  
    83118  int rc = 91;
    84119
    85   cout << " ====== applobe :   Input skycube  name= " << inname << " OutName=" << outname;
     120  cout << " ====== applobe :   Input skycube  name= " << inname << " OutName=" << outname << endl;
     121  cout << ((fggaussian)?" Gaussian ":" Triangular") << ((fixedbeam)?" Fixed (freq.independent)":"") << " beams" << endl;
    86122  bool fginmap=true;
    87123  try {
    … …  
    113149    double lambda = conv.getLambda();
    114150   
    115     Four2DResponse fresp(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
     151    int typbeam=(fggaussian)?1:2;
     152    Four2DResponse fresp(typbeam, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
    116153    Four2DResponse* fresp_p=&fresp;
    117154    Four2DRespTable resptbl;
    … …  
    119156      cout << "applobe[2.b]: initializing Four2DRespTable from file" << resptblname << endl;
    120157      resptbl.readFromPPF(resptblname);
     158      resptbl.renormalize(1.);
    121159      fresp_p=&resptbl;
    122160    }
    … …  
    126164
    127165    if (fgcorrbeam) {
    128       double DoL = tbeamDoL;
     166      double DoL = tbeamDiam/lambda;
    129167      double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
    130       int typcb = 2;
     168      int typcb = (fggaussian)?1:2;
    131169      //      if (fgresptbl) typcb=22;
    132170      Four2DResponse tbeam(typcb, DoL, DoL );
    133       cout << "applobe[3]: calling Correct2RefLobe() with target beam D/Lambda=" << DoL
    134            << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
     171      cout << "applobe[3]: calling Correct2RefLobe() with target beam Diameter=" << tbeamDiam
     172           << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
    135173      beam.Correct2RefLobe(tbeam, incube, dx, dy, Freq0MHz, dfreq);
     174    }
     175    else if (fixedbeam) {
     176      cout << "applobe[3]: calling ApplyLobe() (freq.independent beam) ... " << endl;
     177      beam.ApplyLobe(incube, dx, dy, Freq0MHz);
    136178    }
    137179    else {

  • trunk/Cosmo/RadioBeam/calcpk.cc

    r3829 r3973  
    102102
    103103
    104     HProf hp = pkc.ComputePk(0.,256);
     104    HProf hp = pkc.ComputePk(0.,HPk_NBin);
    105105    {
    106106      cout << "calcpk[4] : writing profile P(k) to PPF file  " << outname << endl;

  • trunk/Cosmo/RadioBeam/calcpk2.cc

    r3830 r3973  
    77    R. Ansari , C. Magneville - Juin 2010
    88
    9 Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile
     9Usage: calcpk2 [-t -g] InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile
    1010               [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]
    1111---------------------------------------------------------------  */
    … …  
    4444{
    4545  if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
    46     cout << " Usage: calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
    47          << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDoL] \n"
     46    cout << " Usage: [-t -g] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
     47         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
    4848         << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
    4949    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
    50       cout << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
     50      cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n"
     51           << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
    5152           << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
    5253           << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
    … …  
    5455           << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n"
    5556           << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
    56            << "- CorBeamDoL: Beam correction target Diameter/Lambda \n"
     57           << "- CorBeamDiam: Beam correction target dish diameter \n"
    5758           << "    These two parameters are used to correct for beam effect for a \n"
    5859           << "    target beam (independent of frequency) defined by D/Lambda \n"
    59            << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5) \n"
     60           << "    DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm) \n"
    6061           << "    default : no beam correction applied \n "
    6162           << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n"
    62            << "    default : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
     63           << "    default (0.) : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
    6364           << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
    6465           << "    foreground subtraction. default is P2 \n"
    … …  
    7475  int rc = 0;
    7576  try {
     77    bool fggaussian=true;  // true -> gaussian beam
     78    // decodage argument optionnel
     79    bool fgoptarg=true;
     80    while (fgoptarg) {
     81      string fbo = arg[1];
     82      if (fbo=="-t")  { fggaussian=false; arg++; narg--; }
     83      else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
     84      else fgoptarg=false;
     85    }
     86    if (narg < 6) {
     87      cout << " calcpk2/error arguments , applobe -h for help " << endl;
     88      return 2;
     89    }
     90
    7691    string inppflss = arg[1];
    7792    r_4 rfaclss = atof(arg[2]);
    … …  
    102117      }
    103118    }
    104     double tbeamDoL=0.;
     119    double tbeamDiam=0.;
    105120    if (narg>8) {
    106       tbeamDoL=atof(arg[8]);
    107       if (tbeamDoL<1.)  fgcorrbeam=false;
     121      tbeamDiam=atof(arg[8]);
     122      if (tbeamDiam<1.)  fgcorrbeam=false;
    108123    }
    109124    bool fgclnsrc=true;
    … …  
    170185      cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
    171186      resptbl.readFromPPF(resptblname);
     187      resptbl.renormalize(1.);
    172188      arep_p=&resptbl;
    173189    }
    … …  
    175191              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    176192
    177     double DoL = tbeamDoL;
     193    double DoL = tbeamDiam/lambda;
    178194    double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
    179     int typcb = 2;
     195    int typcb = (fggaussian)?1:2;
    180196    //    if (fgresptbl) typcb=22;
    181197    Four2DResponse tbeam(typcb, DoL, DoL );
    … …  
    183199    ForegroundCleaner  cleaner(*arep_p, tbeam, skycube);
    184200    if (fgcorrbeam) {
    185       cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam D/Lambda=" << DoL
    186            << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
     201      cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam
     202           << " D/Lambda=" << DoL  << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
    187203      cleaner.BeamCorrections();
    188204    }
    … …  
    218234    pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
    219235   
    220     HProf hp = pkc.ComputePk(0.,256);
     236    HProf hp = pkc.ComputePk(0.,HPk_NBin);
    221237
    222238    tm.Split(" Done ComputePk ");   

  • trunk/Cosmo/RadioBeam/cubedef.h

    r3829 r3973  
    77// Definition tailles et bornes du cube angles (X,Y), frequence (Z)
    88
    9 static sa_size_t NTheta=360;
    10 static sa_size_t NPhi=360;
     9static sa_size_t NTheta=600;
     10static sa_size_t NPhi=800;
    1111static sa_size_t NFreq = 256;
    1212
     13// Carte 30x50 deg, centre en delta = +10 deg , alpha= 150 deg = 10h00
    1314static double Theta0Degre = 65.;  // -> Delta = +30 deg
    14 static double Phi0Degre = 135.;   // -> alpha = 9h00
     15static double Phi0Degre = 130.;   // -> alpha = 8.66 heures
    1516static double ThetaSizeDegre = 30.;  // Taille de la carte en degre selon alpha (axe X)
    16 static double PhiSizeDegre = 30.;    // Taille de la carte en degre selon delta (axe Y)
     17static double PhiSizeDegre = 40.;    // Taille de la carte en degre selon delta (axe Y)
    1718
    1819/* ---  Parametres avec generation Reza from HASLAM + NVSS
    … …  
    3435
    3536// z = 0.6, freq_center=884 MHz , dA_comov ~= 2185 Mpc , 
    36 //       5 arcmin -> 3 Mpc , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
     37//       5 arcmin -> 3 Mpc (3.17 Mpc) , 0.550 MHz -> 3 Mpc  (H(ze) ~ 96 km/s/Mpc)
     38//       3 arcmin -> 1.9 Mpc, 0.5 MHz -> 2.8 Mpc
    3739// Taille de cellule pour le calcul du spectre de puissance 3D
    3840static double ComovDA = 2185.;
    39 static double XCellSizeMpc = 3.;
    40 static double YCellSizeMpc = 3.;
    41 static double ZCellSizeMpc = 3.;
    42 
     41static double XCellSizeMpc = 1.9;
     42static double YCellSizeMpc = 1.9;
     43static double ZCellSizeMpc = 2.8;
     44// Nb de bin de HProf pour calcul P(k)
     45static int HPk_NBin = 384;
    4346
    4447
    … …  
    5255static double sigPLidx2 = 0.15;
    5356
    54 static double PLidxSrc = -2.2;  // index de la loi de puissance des sources 
     57// Generation de la loi de puissance des sources :
     58static double PLidxSrc = -2.0;  // index de la loi de puissance des sources 
    5559static double sigPLidxSrc = 0.15;  // Sigma de la variation (gaussienne) de l'index 
    5660

  • trunk/Cosmo/RadioBeam/lobe.cc

    r3797 r3973  
    1818}
    1919
     20/* --Methode-- */
     21void BeamEffect::ApplyLobe(TArray< TF >& a, double dx, double dy, double f0)
     22// dx, dy en radioans, f0, df en MHz
     23{
     24  Timer tm("BeamEffect::ApplyLobe");
     25  FFTWServer ffts(true);                     
     26  ffts.setNormalize(true);
     27 
     28  H21Conversions conv;
     29  conv.setFrequency(f0);
     30  fresp_.setLambda(conv.getLambda());
     31
     32  TArray< complex<TF> > fourAmp;
     33  double dkx = DeuxPI/(double)a.SizeX()/dx;
     34  double dky = DeuxPI/(double)a.SizeY()/dy;
     35 
     36  for(sa_size_t kz=0; kz<a.SizeZ(); kz++) {
     37    TArray< TF > slice( a(Range::all(), Range::all(), kz) );
     38    ffts.FFTForward(slice, fourAmp);
     39    ApplyLobeK2D(fresp_, fourAmp, dkx, dky);
     40    ffts.FFTBackward(fourAmp, slice, true);
     41    if (kz%20==0)  cout << "BeamEffect::ApplyLobe() done kz=" << kz << " / a.SizeZ()=" << a.SizeZ() << endl;
     42  }
     43  double mean, sigma;
     44  TF min, max;
     45  a.MinMax(min, max);
     46  MeanSigma(a, mean, sigma);
     47  cout << " BeamEffect::ApplyLobe() - Result Min=" << min << " Max=" << max
     48       << " Mean=" << mean << " Sigma=" << sigma << endl;
     49  return;
     50}
    2051
    2152

  • trunk/Cosmo/RadioBeam/lobe.h

    r3797 r3973  
    2727  // Definition de l'objet avec la reponse en frequence de l'instrument
    2828  BeamEffect(Four2DResponse& resp, bool preservefreq0=true);
    29   // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement 
     29
     30  // Applique l'effet d'un lobe fixe (frequency independent) au cube 3D (2 angles, frequence)
     31  void ApplyLobe(TArray< TF >& a, double dx, double dy, double freq); 
     32
     33  // Applique l'effet de lobe au cube 3D (2 angles, frequence), pour chaque plan de frequence successivement
     34  // En faisant varier la reponse de lobe avec la frequence
    3035  void ApplyLobe3D(TArray< TF >& a, double dx, double dy, double f0, double df); 
    3136

  • trunk/Cosmo/RadioBeam/mdish.cc

    r3947 r3973  
    2424  switch (typ_)
    2525    {
    26     case 1:   // Reponse gaussienne parabole diametre D exp[ - 0.5 (lambda  k_g / D )^2 ]
     26    case 1:   // Reponse gaussienne parabole diametre D exp[ -2 (lambda  k_g / D )^2 ]
    2727      wk = sqrt(kx*kx+ky*ky)/dx_;
    28       wk = 0.5*wk*wk;
     28      wk = 2*wk*wk;
    2929      return exp(-wk);
    3030      break;

  • trunk/Cosmo/RadioBeam/mdish.h

    r3947 r3973  
    2121
    2222// -- Four2DResponse : Reponse instrumentale ds le plan k_x,k_y (frequences angulaires theta,phi)
    23 // typ=1 : Reponse gaussienne parabole diametre D exp[ - 0.5  (lambda  k_g / D )^2 ]
     23// typ=1 : Reponse gaussienne parabole diametre D exp[ - 2  (lambda  k_g / D )^2 ]
    2424// typ=2 : Reponse parabole diametre D  Triangle <= kmax= 2 pi D / lambda   
    2525// typ=3 : Reponse rectangle Dx x Dy  Triangle (|kx|,|k_y|) <= (2 pi Dx / lambda, 2 pi Dx / lambda)

  • trunk/Cosmo/RadioBeam/plpkn.pic

    r3947 r3973  
    22### Script de trace de PNoise(k) et reponse dans le plan (u,v)
    33###  de l'interferometre a partir du fichier PPF produit par pknoise.cc
    4 ###     Fev - Avril 2010 ,  BAORadio/Reza
     4###    Janvier 2011 ,  BAORadio/Reza
    55########################################################################
    66
    7 if ( $# < 1 ) then
    8   echo ' Usage: exec plpkn PPFName_pknoise'
    9   return
    10 endif
    11 
    12 echo "---> openppf $1 "
    13 openppf $1
     7echo ' -----> plpknew.pic : opening PPF files ... "
     8delobjs *
     9# openppf ../../PkNoise/cmvhpkz.ppf
     10# openppf ../pkz0p25.ppf
     11openppf ../pkz0p7.ppf
     12rename hpkz hpkz0p7
     13openppf ../pkz1p0.ppf
     14rename hpkz hpkz1p0
     15
     16defscript pknopen
     17  set f $1
     18  set nm $2
     19  echo '------ pknopen File= ' $f ' name=' $nm
     20  openppf $f
     21  rename h1dnoise h1dn$nm
     22  rename h1drep h1dr$nm
     23  rename fracmodok fracok$nm
     24  mv dtnoise pkn$nm
     25  rename h2drep h2dr$nm
     26endscript
     27
     28pknopen pknf11x11.ppf A121d
     29pknopen pknnan128.ppf B128d
     30pknopen pknconfC.ppf  C129d
     31pknopen pknD50m.ppf   D50m
     32pknopen pknD75m.ppf   D75m
     33pknopen pknD100m.ppf   D100m
     34pknopen pknD200m.ppf   D200m
     35pknopen pknf20x20.ppf E400d
     36pknopen pknf4cyl.ppf F4cyl
     37pknopen pknf4cylp.ppf F4cylp
     38pknopen pknf8cyl.ppf G8cyl
     39pknopen pknf8cylp.ppf G8cylp
     40
     41
     42listobjs
     43
     44#  A z = 0.7
     45z = 0.7
     46c = 3.e5
     47H = 102
     48Da = 2488
     49nu21 = 1.42e9
    1450
    1551echo '----> Executing anapkn.pic'
    1652exec anapkn.pic
    1753setup5
    18 scalewz 0.7 2500 100 1
     54# scalewz 0.7 2488 102 1
    1955# scalewz 0.25 989.13 80.26 1
    20 y1 = 500*$cct21
    21 y2 = 9e4*$cct21
    22 xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
     56scalewz 1 3300 120.5 1
     57Da = 3300
     58
     59
     60y1 = 400*$cct21
     61y2 = 7e4*$cct21
     62xyl = "xylimits=0.005,0.5,$y1,$y2 logx logy minorticks"
     63# H = 71.9 km/s/Mpc = 1.0271 x 70 km/s
     64h70 = 1.0271
     65h70cube = $h70*$h70*$h70
     66set kk pow(10.,x)/$h70
     67
    2368defscript plpklss
    24   n/plot hpkz.val*${cct21}%$kk ! ! "same notit nsta connectpoints black $xyl line=solid,2"
    25   addtext 0.03 2000 '*** P(k)-LSS ***' 'font=helvetica,bolditalic,16 black'
     69  scalewz 1 3300 120.5 1
     70#  scalewz 0.7 2488 102 1
     71  convpk2t21
     72  cct21h70 = $cct21*$h70cube
     73  n/plot hpkz1p0.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "
     74#  n/plot hpkz0p7.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "
     75
     76  addtext 0.01 4000 '[ P(k)-LSS z=1.0 ]' 'font=helvetica,bolditalic,16 black'
     77  setaxelabels 'k_comov  (h70 Mpc^-1)'  'P21(k) mK^2 x (Mpc/h70)^3' 'font=helvetica,bolditalic,16'
     78
    2679endscript 
    2780
    … …  
    3386Delnu = 1.e6
    3487
    35 #  A z = 0.7
    36 z = 0.7
    37 c = 3.e5
    38 H = 100
    39 Da = 2500
    40 nu21 = 1.42e9
    4188nu = $nu21/(1+$z)
    4289pi23 = 8.*Pi*Pi*Pi
    4390
    44 Lambda = 0.357
     91Lambda = 0.21*(1+$z)
    4592Lam2 = $Lambda*$Lambda
    4693
    … …  
    5299  Dyol = $Dy/$Lambda
    53100
    54   FOV = (1.2*1.2*$Lam2/$Dx/$Dy)
     101#   FOV = (1.2*1.2*$Lam2/$Dx/$Dy)
     102  FOV = ($Lam2/$Dx/$Dy)
     103
    55104  FOVDEG = $FOV*$R2D2
    56105  NPointing = 10000/$FOVDEG
    57106  tinteg = 365*86400/$NPointing
    58   PNOISE = $Tsys*$Tsys/$tinteg/$Dxol/$Dyol
     107  PNOISE = 2.*$Tsys*$Tsys/$tinteg/$Dxol/$Dyol
    59108  PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu
    60109#   PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu/$pi23
    61110  PNOISE = $PNOISE*1.e6
     111  PNOISE = $PNOISE*1.05
    62112  echo " FOV = $FOV deg^2  NPointing= $NPointing"
    63113  echo " tinteg= $tinteg sec  PNOISE= $PNOISE mK^2"
    64114endscript
    65115
    66 defscript plnoisedish
     116
     117defscript plnoiseDishes
     118  Dx = 75
     119  Dy = 75
     120  calcul
     121  plot2d pknD75m  k/$Da pnoise*$PNOISE/5 fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
    67122  Dx = 100
    68123  Dy = 100
    69124  calcul
    70   plot2d noiseD x/$Da val*$PNOISE nb>10 'line=solid,2 logy logx xylimits=0.002,0.8,1,1e5 navyblue grid cpts nsta notit'
     125  plot2d pknD100m  k/$Da pnoise*$PNOISE/20 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
     126
     127  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'
     128
     129
     130
     131  set lines ( '(a) 75m Dish, 5 beams'  '(b) 100m Dish, 10 beams'  '(b) 100m Dish, 100 beams'  )
     132  set cols ( red blueviolet blue )
     133  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35'
     134  settitle ' PNoise(k) : 5/10/100 beams/polar @z=1' ' ' 'font=helvetica,bold,16'
     135endscript
     136
     137defscript plnoiseEFGH
     138  Dx = 12*0.9
     139  Dy = 0.844248*0.8
     140  calcul
     141  plot2d pknF4cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts forestgreen nsta notit'
     142  plot2d pknF4cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts green nsta notit'
     143
     144  plot2d pknG8cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violetred nsta notit'
     145  plot2d pknG8cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violet nsta notit'
     146
     147  Dx = 5*0.9
     148  Dy = 5*0.9
     149  calcul
     150  plot2d pknE400d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts turquoise nsta notit'
     151
     152  set lines (  '(e) 20x20:400xD=5m' '(f) 4Cyl-12mx85m, 400 rec/pol' '(fp) 4Cylp-12mx85m, 400 rec/pol'  '(g) 8Cyl-12mx105m, 960 rec/pol' '(gp) 8Cylp-12mx105m, 960 rec/pol' )
     153  set cols ( turquoise forestgreen green violetred violet  )
     154  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35'
     155  settitle ' PNoise(k) : Dishes/Cylinders, 400/400/960 recv/pol @z=1' ' ' 'font=helvetica,bold,16'
     156endscript
     157
     158defscript plnoiseABCD
     159  Dx = 5*0.9
     160  Dy = 5*0.9
     161  calcul
     162  plot2d pknA121d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts magenta nsta notit'
     163  plot2d pknB128d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
     164  plot2d pknC129d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts orange nsta notit'
     165
     166  Dx = 75
     167  Dy = 75
     168  calcul
     169  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
     170
     171  set lines (  '(a) 11x11:121xD=5m' '(b) 128xD=5m' '(c) 129xD=5m' '(d) 75m Dish, 100 beams' )
     172  set cols ( magenta red orange blueviolet )
     173  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35'
     174  settitle ' PNoise(k) : 121,128,129,100 receiver/polar @z=1' ' ' 'font=helvetica,bold,16'
     175endscript
     176
     177defscript plfracABCD
     178  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
     179  plot2d pknA121d k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts magenta nsta notit'
     180  plot2d pknB128d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts red nsta notit'
     181  plot2d pknC129d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts orange nsta notit'
     182  plot2d pknD75m  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts violetblue nsta notit'
     183endscript
     184
     185defscript plfracEFGH
     186  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
     187  plot2d pknF4cyl k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts forestgreen nsta notit'
     188  plot2d pknF4cylp k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts green nsta notit'
     189  plot2d pknG8cyl k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts gold nsta notit'
     190  plot2d pknG8cylp  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts yellow nsta notit'
     191  plot2d pknE400d  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts magenta nsta notit'
     192
     193endscript
     194
     195defscript pldishes
     196
     197  Dx = 50
     198  Dy = 50
     199  calcul
     200  plot2d pknD50m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts skyblue nsta notit'
     201
     202  Dx = 75
     203  Dy = 75
     204  calcul
     205  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'
     206
     207  Dx = 100
     208  Dy = 100
     209  calcul
     210  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts navyblue nsta notit'
     211
    71212  Dx = 200
    72213  Dy = 200
    73214  calcul
    74   plot2d noiseD2 x/$Da val*$PNOISE nb>10 'same line=solid,2  cpts grey nsta notit'
    75   set lines ( '100mDish' '200mDish'  )
    76   set cols ( navyblue grey  )
    77   textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.1,0.3,0.7,0.8'
    78   setaxelabels 'k (Mpc^-1) ' 'PNoise(k) mk^2 Mpc^3' 'font=helvetica,bolditalic,16'
    79 endscript
    80 
    81 defscript plnoiseA
    82   Dx = 5*0.95
    83   Dy = 5*0.95
    84   calcul
    85   plot2d noisemdf64 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
    86 #  plot2d noisemds x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
    87   plot2d noisemdsB x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts siennared nsta notit'
    88   plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violetred nsta notit'
    89   Dx = 25*0.3
    90   Dy = 0.5*0.9
    91   calcul
    92   plot2d noise2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts forestgreen nsta notit'
    93   plot2d noise2cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
    94   set lines (  'FilledA:64x5mD' 'SparseB:72x5mD' 'SparseC:129x5mD' 'Pitts2Cyl(=64C)' 'PerfPitts2Cyl(=64C)' )
    95   set cols ( red siennared magenta  forestgreen green )
     215  plot2d pknD200m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'
     216
     217  set lines (  'D:Dish50m' 'D:Dish75m' 'D:Dish100m' 'D:Dish200m' )
     218  set cols ( skyblue blue navyblue blueviolet )
    96219  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35'
    97   settitle ' PNoise(k) : Dishes/Cylinders, 64/72/129 channels' ' ' 'font=helvetica,bold,16'
    98 endscript
    99 
    100 defscript plnoiseB
    101   Dx = 5*0.95
    102   Dy = 5*0.95
    103   calcul
    104   plot2d noisemdsC x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
    105   plot2d noisemdf x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts gold nsta notit'
    106   Dx = 10*0.95
    107   Dy = 0.5*0.9
    108   calcul
    109   plot2d noise3cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit'
    110   plot2d noise3cylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts magenta nsta notit'
    111 
    112   Dx = 12*0.95
    113   Dy = 0.5*0.9
    114   calcul
    115   plot2d noisefcyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts blue nsta notit'
    116   plot2d noisefcylP x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts royalblue nsta notit'
    117 
    118   set lines ( 'SparseC:129x5mD' 'Filled400x5m' '3xCyl10x64m(=384C)' 'Perf3xCyl10x64m(=384C)' )
    119   set lines ( $lines '8Cyl12x96m(=1536C)' 'Perf8Cyl12x96m(=1536C)' )
    120   set cols ( red gold  violet magenta  blue royalblue )
    121   textdrawer lines cols 'frame font=helvetica,bold,16 inset=0.25,0.55,0.6,0.95'
    122   settitle ' PNoise(k) : Dishes/Cylinders, 129/384/400/1536 channels' ' ' 'font=helvetica,bold,16'     
    123 endscript
    124 
    125 defscript xxx
    126 newwin
    127 disp mfill 'h2disp=img colbr128 notit'
    128 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    129 settitle ' u-v coverage , 400 x 5m Dishes - No Pointing' ' ' 'font=helvetica,bold,16'
    130 
    131 disp dish 'h2disp=img colbr128 notit'
    132 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    133 settitle ' u-v coverage , 100 m Dish' ' ' 'font=helvetica,bold,16'
    134 
    135 disp msparsfp 'h2disp=img colbr128 notit'
    136 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    137 settitle ' u-v coverage , 63 x 5m Dishes T-config - No Pointing' ' ' 'font=helvetica,bold,16'
    138 
    139 disp mspars 'h2disp=img colbr128 notit'
    140 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    141 settitle ' u-v coverage , 63 x 5m Dishes T-config - Pointing ~Pi/4' ' ' 'font=helvetica,bold,16'
    142 
    143 disp mcylf 'h2disp=img colbr128 notit'
    144 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    145 settitle ' u-v coverage , Filled Cylinder Array 8 Cyl 12mx96m ' ' ' 'font=helvetica,bold,16'
    146 
    147 disp mcylfP 'h2disp=img colbr128 notit'
    148 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    149 settitle ' u-v coverage , Perfect Filled Cylinder Array 8 Cyl 12mx96m ' ' ' 'font=helvetica,bold,16'
    150 
    151 disp m2cyl 'h2disp=img colbr128 notit'
    152 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    153 settitle ' u-v coverage , Pittsburgh 2 Cyl 16mx8m , dist=25m ' ' ' 'font=helvetica,bold,16'
    154 
    155 disp m2cylP 'h2disp=img colbr128 notit'
    156 setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    157 settitle ' u-v coverage , Perfect-Pitts. 2 Cyl 16mx8m , dist=25m ' ' ' 'font=helvetica,bold,16'
    158 
    159 
    160 disp noiseD 'logy nsta'
    161 disp noiseD2 'same grey nsta'
    162 disp noisemdf 'same cyan nsta'
    163 disp noisemdf64 'same brown nsta'
    164 disp noisemds 'same red nsta'
    165 disp noisemdsC 'same orange nsta'
    166 
    167 disp noisemdsB 'same yellow nsta'
    168 disp noisemdsfp 'same yellow nsta'
    169 # disp noisemdsd7 'same gold nsta'
    170 disp noisefcyl 'same blue nsta'
    171 disp noisefcylP 'same skyblue nsta'
    172 disp noise3cylP 'same magenta nsta'
    173 disp noise3cyl 'same violet nsta'
    174 disp noise2cyl 'same forestgreen nsta'
    175 disp noise2cylP 'same green nsta'
    176 
    177 endscript
    178 
    179 defscript AA
    180   plnoisedish
    181   plpklss
    182   plnoiseA
     220  settitle ' PNoise(k) : Dish D=50m,75m,100m,200m' ' ' 'font=helvetica,bold,16'
     221endscript
     222
     223
     224
     225defscript plnoisenancay
     226#  Dx = 5*0.95
     227#  Dy = 5*0.95
     228#  calcul
     229#  plot2d pknnan24 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     230#  plot2d pknnan25 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts orange nsta notit'
     231#  Dx = 3.5*0.95
     232#  Dy = 3.5*0.95
     233#  calcul
     234#  plot2d pknnan36 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     235#  plot2d pknnan40 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit'
     236#  Dx = 7*0.9
     237#  Dy = 2*$Lambda*0.9
     238#  calcul
     239#  plot2d pknpit2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
     240  Dx = 9*0.9
     241  Dy = 2*$Lambda*0.9
     242  calcul
     243  plot2d pknpit2cylw x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
     244  set lines ( '2 Cylinders 20mx9m' )
     245  set cols ( red )
     246  textdrawer lines cols 'noframe font=helvetica,bold,16 inset=0.4,0.4,0.4,0.4'
     247
     248endscript
     249
     250defscript ABCD
     251#  plnoisedish
     252#  zone 1 2
     253#  newwin 1 1 1000 600
     254  plpklss
     255  plnoiseABCD
     256#  plfracABCD
     257endscript
     258
     259defscript EFGH
     260#  plnoisedish
     261#  zone 1 2
     262#  newwin 1 1 900 600
     263  plpklss
     264  plnoiseEFGH
     265#  plfracEFGH
     266endscript
     267
     268defscript DDDD
     269#  plnoisedish
     270#  zone 1 2
     271  plpklss
     272  plnoiseDishes
     273#  plfracABCD
    183274endscript
    184275
    … …  
    189280endscript
    190281
     282defscript CC
     283  plnoisedish
     284  plpklss
     285  plnoiseC2
     286endscript
     287
     288defscript nancay
     289  plnoisedish
     290  plpklss
     291  plnoisenancay
     292endscript
     293
    191294defscript POSCOV
    192   disp posspB red
     295  openppf hdt_repnan128.ppf
     296  rename mdish mdB128d
     297  rename h2rep uvB128d
     298
     299  openppf hdt_repf11x11.ppf
     300  rename mdish mdA121d
     301  rename h2rep uvA121d
     302
     303  openppf hdt_repconfC.ppf
     304  rename mdish mdC129d
     305  rename h2rep uvC129d
     306
     307  newwin 2 1 800 400
     308  zone 2  1
     309  nt2d mdB128d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red '
    193310  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16'
    194   settitle ' Config B dish positions - 72 dishes  ' ' ' 'font=helvetica,bold,16'
    195 #  w2ps
    196    
    197   disp posspC red
     311  settitle '(b) 128 D=5m dishes in 8 rows ' ' ' 'font=helvetica,bold,16'
     312  nt2d mdC129d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red'
    198313  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16'
    199   settitle ' Config C dish positions - 129 dishes  ' ' ' 'font=helvetica,bold,16'
    200 #  w2ps
    201   disp mfill64 'h2disp=img colbr128 notit'
    202   setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    203   settitle ' u-v coverage , Filled 8x8 - 64 x 5m Dishes' ' ' 'font=helvetica,bold,16'
    204 #  w2ps
    205   disp m3cyl 'h2disp=img colbr128 notit'
    206   setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    207   settitle ' u-v coverage , 3 Cylinders 10mx64m ' ' ' 'font=helvetica,bold,16'
    208 #  w2ps
    209   disp mfill 'h2disp=img colbr128 notit'
    210   setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    211   settitle ' u-v coverage , Filled 20x20 - 400 x 5m Dishes' ' ' 'font=helvetica,bold,16'
    212 #  w2ps
    213   disp msparsC 'h2disp=img colbr128 notit'
    214   setaxelabels 'kx (Radian^-1)  k=1000 -> ~21 arcmin ' ' ky (Radian^-1) ' 'font=helvetica,bolditalic,16'
    215   settitle 'u-v coverage, Sparse-C: 129x5mD Over 80mx80m (Rot~Pi/4)' ' ' 'font=helvetica,bold,16'
     314  settitle '(c) 129 D=5m dishes ' ' ' 'font=helvetica,bold,16'
     315  pssetfilename configab.ps
     316  w2ps
     317  w2eps configab.eps
     318  psclosefile
     319
     320  newwin 2 2 800 800
     321  disp uvA121d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
     322  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12'
     323  settitle 'u-v coverage, (a) 11x11 D=5m dishes Over 55mx55m' ' ' 'font=helvetica,bold,12'
     324  disp uvB128d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
     325  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12'
     326  settitle 'u-v coverage, (b) 8 row of 16xD=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12'
     327  disp uvC129d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
     328  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12'
     329  settitle 'u-v coverage, (c) 129 D=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12'
     330  disp h2drD75m 'h2disp=img colbr128 h2dyn=0.02,0.8 notit nsta'
     331  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12'
     332  settitle 'u-v coverage, (d) D=75 m dish with 100 beams' ' ' 'font=helvetica,bold,12'
     333
     334  pssetfilename uvcovabcd.ps
     335  w2ps
     336  w2eps uvcovabcd.eps
     337  psclosefile
     338
     339#  imag uvB128d 'colbr128 stdaxes showcmap=right'
    216340 
    217341endscript 

  • trunk/Cosmo/RadioBeam/srcat2cube.cc

    r3829 r3973  
    278278      lowoksrccnt++;
    279279    }
    280     else slo=5.;
     280    else slo=-1.;
    281281    for (sa_size_t k=0; k<ocube.SizeZ(); k++)  {
    282282      double rapfreq = pow((freq0+k*dfreq)/infreq, slo);

  • trunk/Cosmo/RadioBeam/subtractradsrc.cmd

    r3830 r3973  
    1212csh> ~/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
    1313
     14# 1.c/ To run SimLSS with GSM map parametersand 40x40 deg maps (DeltaFreq=500 MHz)
     15csh> ~/Objs/exe/cmvginit3df -a -1 -2 -C -G 0. -F 0 -x 600,1.9 -y 800,1.9 -z 256,2.8 -Z 0.60 -8 1. -n 10000 -O 0,2 -o lssz060 -T 2
     16
    1417# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention
    1518#  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray)
    … …  
    1821
    1922csh> cat > racube.pic
    20 set f lssz056
     23set f lssz060
    2124readfits ${f}_r.fits
    2225rename ${f}_r map
    … …  
    3336print lsscube
    3437# expmeansig lsscube val
    35 saveppf lsscube lsscube.ppf
     38saveppf lsscube lsscubez060.ppf
    3639
    3740csh> spiapp -term -exec racube.pic
    3841
     42#### Cube LSS 40x30 deg (3')  @ z=0.6 ( lsscubez060.ppf )
     43#### -> Size= 122880000 Mean=-7.01664e-05 Sigma=2.53016 Min=-13.7439 Max=14.4648
    3944
    4045## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin)
    4146# 2.a/ Synchrotron map from HASLAM 400 MHz map
    42 csh> ./Objs/syncube syncmap_eq.fits syncube.ppf
     47csh> ./Objs/syncube syncmap_eq.fits syncube.ppf syncmap.ppf
    4348# 2.b/ radio source cube from NVSS catalog
    44 csh> ./Objs/srcat2cube nvss.fits nvsscube.ppf
     49csh> ./Objs/srcat2cube -nvss nvss.fits nvsscube.ppf nvssmap.ppf
     50# Or from the north20 catalog :
     51csh> ./Objs/srcat2cube -north20 north20cm.fits north20cube.ppf north20map.ppf
     52
    4553# 2.c/ Add the two cubes using the following spiapp script
    4654csh> cat > sumcubes.pic
    4755openppf syncube.ppf
    48 openppf srcnor3d.ppf
     56openppf radsrccube.ppf
    4957# expmeansig syncube val
    50 # expmeansig srcnor3d val
    51 c++exec TArray<r_4> fgndcube = syncube+srcnor3d; KeepObj(fgndcube);
     58# expmeansig nvsscube val
     59c++exec TArray<r_4> fgndcube = syncube+radsrccube; KeepObj(fgndcube);
    5260print fgndcube
    5361# expmeansig fgndcube val
    … …  
    5664csh> spiapp -term -exec sumcubes.pic
    5765
     66#### syncube:Mean= 1.8101  Sigma= 0.326538  Min= 0.857019 Max= 3.58987
     67#### nvsscube: Mean= 1.95073  Sigma= 1.68515  Min= 0.857019 Max= 428.398
     68#### fgndcube=syncube+nvsscube: Mean= 0.140623  Sigma= 1.65068  Min= 0 Max= 426.559
     69#### north20: fgndcube_north:
     70
     71## Step 2.b/ Produce foreground cube from GSM
     72csh> ./Objs/gsm2cube ../Catalogs/GSM/ 1 256 fgndcube_gsm.ppf
     73
    5874## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
    59 csh> ./Objs/applobe 50. fgndcube.ppf fgndcube_lobe.ppf
    60 csh> ./Objs/applobe 50. lsscube.ppf lsscube_lobe.ppf
    61 ## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 130 
    62 csh> ./Objs/applobe 50. lsscube_lobe.ppf lsscube_corlobe.ppf 130
    63 csh> ./Objs/applobe 50. fgndcube_lobe.ppf fgndcube_corlobe.ppf 130
     75csh> set ddish=55.
     76csh> set ddishcor=50.
     77csh> ./Objs/applobe $ddish fgndcube.ppf fgndcube_lobe.ppf
     78csh> ./Objs/applobe -fib $ddish fgndcube.ppf fgndcube_flobe.ppf
     79csh> ./Objs/applobe $ddish lsscube.ppf lsscube_lobe.ppf
     80csh> ./Objs/applobe -fib $ddish lsscube.ppf lsscube_flobe.ppf
     81## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150  (55 m @ z=0.7 - 820 MHz) 
     82csh> ./Objs/applobe $ddish lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
     83csh> ./Objs/applobe $ddish fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor
     84
     85## Step 3.c/ Apply lobe (Filled 11x11 5m dishes array) effect on foreground cube and LSS cube
     86csh> ./Objs/applobe repf11x11.ppf fgndcube.ppf fgndcube_lobe.ppf
     87csh> ./Objs/applobe repf11x11.ppf lsscube.ppf lsscube_lobe.ppf
     88## Step 3.d/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 150 (55 m @ z=0.7 - 820 MHz)
     89csh> ./Objs/applobe repf11x11.ppf lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
     90csh> ./Objs/applobe repf11x11.ppf fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor
    6491
    6592### Step 4/ Compute power spectra
    66 ## mass to temperature converion factor   CT21 ~= 0.2 mK
     93## mass to temperature converion factor   CT21 ~= 0.21 mK for gHI=2% , 0.11 for gHI=1% , 0.13 for gHI=0.008x(1+0.6)
    6794## Foreground maps are in temperature
    6895## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
    69 ## Tsys ~ 50 K , DeltaFreq ~ 0.275 MHz , t_obs ~ 1 day ~ 80 000 s.
    70 ## sigma_noise ~ 0.35 mK
     96## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
     97## sigma_noise ~ 0.25 mK -> 3 mK
    7198# 4.a/ LSS power spectrum  without noise
    72 csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.2
     99csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.13
    73100# and with noise 
    74 csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.35
     101csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.13  3
    75102#  with the lobe effect
    76 csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2 
    77 csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.2 0.35 
    78 csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2 
     103csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.13
     104csh> ./Objs/calcpk lsscube_flobe.ppf lsspkflobe.ppf 0.13
     105csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.13 3 
     106csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.13 
    79107
    80108# 4.b/ Foreground power spectrum
    81109csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
    82110csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
     111csh> ./Objs/calcpk fgndcube_flobe.ppf fgndpkflobe.ppf 1000
    83112csh> ./Objs/calcpk fgndcube_corlobe.ppf fgndpkcorlobe.ppf 1000
    84113
    85114# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam
    86 csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 0. 0. P2
     115csh> set beamdesc=repf11x11.ppf
     116csh> set ddishcor=50.
     117csh> set noiselev=1.
     118csh> ./Objs/calcpk2 lsscube.ppf 0.13 fgndcube.ppf 1000 subpk.ppf $noiselev $beamdesc 0. 0. P2
    87119# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction 
    88 csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. P2
    89 # 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130
    90 csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3. P2
     120csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpklobe.ppf $noiselev $beamdesc 0. 0. P2
     121# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam= $ddishcor
     122csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2
    91123#  Or using a linear fit for foreground subtraction (old version)
    92 csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf 0.35 50. 130. 3. P1
     124csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf $noiselev $beamdesc $ddishcor 0. P2
    93125# 4.f / Estimate residual noise from Foreground removal :
    94 csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3. P2
    95 csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.
     126csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residcorlobe.ppf  $noiselev $beamdesc $ddishcor 0. P2
     127csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residnocor.ppf $noiselev $beamdesc 0. 0. P2
    96128
    97129### Step 5 / Check the results using spiapp
    … …  
    100132openppf fgndpk.ppf
    101133openppf fgndpklobe.ppf
     134openppf fgndpkflobe.ppf
    102135openppf fgndpkcorlobe.ppf
    103136openppf lsspk.ppf
    104137openppf lsspklobe.ppf
     138openppf lsspkflobe.ppf
    105139openppf lsspklobewn.ppf
    106140openppf subpklobe.ppf
    … …  
    136170
    137171#  Calcul du volume total en Mpc^3
    138 set VOL 3*3*3*360*360*256
     172set VOL (1.9*1.9*2.8*800*600*256)
    139173plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
    140174plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared'
Note: See TracChangeset for help on using the changeset viewer.