source: Sophya/trunk/Cosmo/RadioBeam/subtractradsrc.cmd@ 3973

#####################################################################################
####  Commands to run the different programs to produce foreground maps 
####  and compute radio-source subtracted P(k)
#####################################################################################

### Cube definition in file cubedef.h 

### Step 1/ Produce an LSS data cube  with appropriate size and redshift using SimLSS
# 1.a/ Run SimLSS
csh> ~/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
# 1.b/ To run SimLSS with GSM map parameters (DeltaFreq=500 MHz)
csh> ~/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

# 1.c/ To run SimLSS with GSM map parametersand 40x40 deg maps (DeltaFreq=500 MHz)
csh> ~/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

# 1.c/ Change the X and Z axis of the cube to adapt it to RadioBeam package convention 
#  SimLSS output : the radial (redshift) direction along X axis of the cube (TArray) 
#  RadioBeam cubes : the radial (redshift) direction along Z axis of the cube (TArray) 
#  Execucte the following script in spiapp :

csh> cat > racube.pic
set f lssz060
readfits ${f}_r.fits
rename ${f}_r map
print map 
c++exec \
  TArray<r_4> omap(map.SizeY(),map.SizeZ(),map.SizeX()-2 ); \
  for(sa_size_t i=0;i<omap.SizeX();i++) \
    for(sa_size_t j=0;j<omap.SizeY();j++) \
      for(sa_size_t k=0;k<omap.SizeZ();k++) \
        omap(i,j,k)=map(k+1,i,j);  \
  KeepObj(omap); 

rename omap lsscube
print lsscube
# expmeansig lsscube val
saveppf lsscube lsscubez060.ppf 

csh> spiapp -term -exec racube.pic 

#### Cube LSS 40x30 deg (3')  @ z=0.6 ( lsscubez060.ppf ) 
#### -> Size= 122880000 Mean=-7.01664e-05 Sigma=2.53016 Min=-13.7439 Max=14.4648

## Step 2/ Produce synchrotron and radio source sky cubes  (cube unit is Temparature- Kelvin) 
# 2.a/ Synchrotron map from HASLAM 400 MHz map
csh> ./Objs/syncube syncmap_eq.fits syncube.ppf syncmap.ppf
# 2.b/ radio source cube from NVSS catalog
csh> ./Objs/srcat2cube -nvss nvss.fits nvsscube.ppf nvssmap.ppf
# Or from the north20 catalog :
csh> ./Objs/srcat2cube -north20 north20cm.fits north20cube.ppf north20map.ppf

# 2.c/ Add the two cubes using the following spiapp script 
csh> cat > sumcubes.pic
openppf syncube.ppf
openppf radsrccube.ppf
# expmeansig syncube val 
# expmeansig nvsscube val
c++exec TArray<r_4> fgndcube = syncube+radsrccube; KeepObj(fgndcube);
print fgndcube
# expmeansig fgndcube val
saveppf fgndcube fgndcube.ppf

csh> spiapp -term -exec sumcubes.pic 

#### syncube:Mean= 1.8101  Sigma= 0.326538  Min= 0.857019 Max= 3.58987
#### nvsscube: Mean= 1.95073  Sigma= 1.68515  Min= 0.857019 Max= 428.398
#### fgndcube=syncube+nvsscube: Mean= 0.140623  Sigma= 1.65068  Min= 0 Max= 426.559
#### north20: fgndcube_north:

## Step 2.b/ Produce foreground cube from GSM 
csh> ./Objs/gsm2cube ../Catalogs/GSM/ 1 256 fgndcube_gsm.ppf

## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
csh> set ddish=55.
csh> set ddishcor=50.
csh> ./Objs/applobe $ddish fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe -fib $ddish fgndcube.ppf fgndcube_flobe.ppf
csh> ./Objs/applobe $ddish lsscube.ppf lsscube_lobe.ppf
csh> ./Objs/applobe -fib $ddish lsscube.ppf lsscube_flobe.ppf
## 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)  
csh> ./Objs/applobe $ddish lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
csh> ./Objs/applobe $ddish fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor 

## Step 3.c/ Apply lobe (Filled 11x11 5m dishes array) effect on foreground cube and LSS cube
csh> ./Objs/applobe repf11x11.ppf fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe repf11x11.ppf lsscube.ppf lsscube_lobe.ppf
## 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) 
csh> ./Objs/applobe repf11x11.ppf lsscube_lobe.ppf lsscube_corlobe.ppf $ddishcor
csh> ./Objs/applobe repf11x11.ppf fgndcube_lobe.ppf fgndcube_corlobe.ppf $ddishcor

### Step 4/ Compute power spectra 
## 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)
## Foreground maps are in temperature 
## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
## Tsys ~ 50 K , DeltaFreq ~ 0.5 MHz , t_obs ~ 1 day ~ 80 000 s.
## sigma_noise ~ 0.25 mK -> 3 mK  
# 4.a/ LSS power spectrum  without noise
csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.13
# and with noise  
csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.13  3
#  with the lobe effect 
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.13
csh> ./Objs/calcpk lsscube_flobe.ppf lsspkflobe.ppf 0.13
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.13 3  
csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.13  

# 4.b/ Foreground power spectrum 
csh> ./Objs/calcpk fgndcube.ppf fgndpk.ppf 1000
csh> ./Objs/calcpk fgndcube_lobe.ppf fgndpklobe.ppf 1000
csh> ./Objs/calcpk fgndcube_flobe.ppf fgndpkflobe.ppf 1000
csh> ./Objs/calcpk fgndcube_corlobe.ppf fgndpkcorlobe.ppf 1000

# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
csh> set beamdesc=repf11x11.ppf
csh> set ddishcor=50.
csh> set noiselev=1.
csh> ./Objs/calcpk2 lsscube.ppf 0.13 fgndcube.ppf 1000 subpk.ppf $noiselev $beamdesc 0. 0. P2
# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect, without beam correction  
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpklobe.ppf $noiselev $beamdesc 0. 0. P2 
# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam= $ddishcor
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf $noiselev $beamdesc $ddishcor 0. P2 
#  Or using a linear fit for foreground subtraction (old version)
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.13 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf $noiselev $beamdesc $ddishcor 0. P2 
# 4.f / Estimate residual noise from Foreground removal : 
csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residcorlobe.ppf  $noiselev $beamdesc $ddishcor 0. P2
csh> ./Objs/calcpk2 lsscube.ppf 0. fgndcube_lobe.ppf 1000 residnocor.ppf $noiselev $beamdesc 0. 0. P2

### Step 5 / Check the results using spiapp 
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
delobjs *
openppf fgndpk.ppf 
openppf fgndpklobe.ppf 
openppf fgndpkflobe.ppf 
openppf fgndpkcorlobe.ppf 
openppf lsspk.ppf 
openppf lsspklobe.ppf 
openppf lsspkflobe.ppf 
openppf lsspklobewn.ppf 
openppf subpklobe.ppf 

openppf subpkcorlobe.ppf 
openppf subpknolss.ppf

# openppf subpknolssnocor.ppf

disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-11,8e-6 gold'
disp lsspklobe 'same nsta orange'
disp lsspklobewn 'same nsta siennared'
settitle ' Pk[LSS] - without normalisation' ' ' 'font=helvetica,bold,16 black' 


disp fgndpk 'logx logy nsta xylimits=0.005,2.,1e-10,1. navyblue'
disp fgndpklobe 'same nsta blue'
disp fgndpkcorlobe 'same nsta skyblue'
disp lsspk 'same nsta gold'
disp lsspklobewn 'same nsta siennared' 
# settitle 'Pk[LSS] , Pk[Foreground] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'
settitle 'Pk[LSS] , Pk[Foreground=GSM] and lobe effect (Dish D=50 m)' ' ' 'font=helvetica,bold,18'

set lines ( 'Pk[Foreground]'  'Pk[fgnd]*Lobe' 'Pk[fgnd]*Lobe/Corrected' 'Pk[LSS]' 'Pk[LSS]*Lobe+Noise' )
set cols ( navyblue blue skyblue gold siennared ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'


disp lsspk 'logx logy nsta xylimits=0.005,2.,4e-9,4e-5 gold'
disp lsspklobewn 'same nsta siennared' 
disp subpkcorlobe 'same nsta red'
disp subpknolss 'same nsta green'

#  Calcul du volume total en Mpc^3
set VOL (1.9*1.9*2.8*800*600*256)
plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 siennared' 
plot2d subpkcorlobe x val*$VOL 1 'same nsta cpts marker=box,5 red' 
plot2d subpklobe x val*$VOL 1 'same nsta cpts marker=box,5 blueviolet' 
plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 

# settitle 'Recovered Pk[LSS] In=LSS+(GSM) (D=50 m)' ' ' 'font=helvetica,bold,18'
settitle 'Recovered Pk[LSS] In=LSS+(Haslam+North20cm) (D=50 m)' ' ' 'font=helvetica,bold,18'
setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[ExtractedLSS]' 'Pk[ExtLSS,NoBeamCor]' 'Pk[residual,NoLSS]' )
set cols ( gold siennared red blueviolet green ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'


plot2d lsspk x val*$VOL 1 'logx logy nsta xylimits=0.01,2.,10.,1e4 cpts marker=box,5 gold'
plot2d lsspklobewn x val*$VOL 1 'same nsta cpts marker=box,5 red' 
plot2d subpknolss x val*$VOL 1 'same nsta cpts marker=box,5 green' 
plot2d subpknolssnocor x val*$VOL 1 'same nsta cpts marker=box,5 magenta' 
setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=helvetica,bolditalic,16' 
settitle 'Recovered Pk[LSS] and residual systematics' ' ' 'font=helvetica,bold,18'
set lines ( 'Pk[LSS]'  'Pk[LSS*lobe+noise]' 'Pk[residual,NoLSS]' 'Pk[residual,NoLSS,NoBeamCorrection]' )
set cols ( gold red green magenta ) 
textdrawer lines cols 'font=helvetica,bold,16 frame'