#####################################################################################
####  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/ 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 lssz056
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 lsscube.ppf 

csh> spiapp -term -exec racube.pic 


## 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
# 2.b/ radio source cube from NVSS catalog
csh> ./Objs/srcat2cube nvss.fits nvsscube.ppf
# 2.c/ Add the two cubes using the following spiapp script 
csh> cat > sumcubes.pic
openppf syncube.ppf
openppf srcnor3d.ppf
# expmeansig syncube val 
# expmeansig srcnor3d val
c++exec TArray<r_4> fgndcube = syncube+srcnor3d; KeepObj(fgndcube);
print fgndcube
# expmeansig fgndcube val
saveppf fgndcube fgndcube.ppf

csh> spiapp -term -exec sumcubes.pic 

## Step 3/ Apply lobe (50 meter diameter array) effect on foreground cube and LSS cube
csh> ./Objs/applobe 50. fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe 50. lsscube.ppf lsscube_lobe.ppf
## Step 3.b/ Correct for the lobe effect by bringing all to the beam of Diam/Lambda = 130  
csh> ./Objs/applobe 50. lsscube_lobe.ppf lsscube_corlobe.ppf 130
csh> ./Objs/applobe 50. fgndcube_lobe.ppf fgndcube_corlobe.ppf 130

### Step 4/ Compute power spectra 
## mass to temperature converion factor   CT21 ~= 0.2 mK 
## Foreground maps are in temperature 
## Noise fluctuations Sigma^2 ~ T_sys^2 / t_obs * DeltaFreq
## Tsys ~ 50 K , DeltaFreq ~ 0.275 MHz , t_obs ~ 1 day ~ 80 000 s.
## sigma_noise ~ 0.35 mK 
# 4.a/ LSS power spectrum  without noise
csh> ./Objs/calcpk lsscube.ppf lsspk.ppf 0.2 
# and with noise  
csh> ./Objs/calcpk lsscube.ppf lsspkwn.ppf 0.2 0.35 
#  with the lobe effect 
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobe.ppf 0.2  
csh> ./Objs/calcpk lsscube_lobe.ppf lsspklobewn.ppf 0.2 0.35  
csh> ./Objs/calcpk lsscube_corlobe.ppf lsspkcorlobe.ppf 0.2  

# 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_corlobe.ppf fgndpkcorlobe.ppf 1000

# 4.c/ Extract LSS P(k) from Foreground+LSS+noise , after cleaning/subtraction without beam 
csh> ./Objs/calcpk2 lsscube.ppf 0.2 fgndcube.ppf 1000 subpk.ppf 0.35 50. 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.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 50. 0. 3. P2 
# 4.e / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting to a beam of Diam/Lambda = 130 
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobe.ppf 0.35 50. 130. 3. P2 
#  Or using a linear fit for foreground subtraction (old version)
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpkcorlobep1.ppf 0.35 50. 130. 3. P1
# 4.f / Estimate residual noise from Foreground removal : 
csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolss.ppf 0.35 50. 130. 3. P2
csh> ./Objs/calcpk2 zerolss.ppf 0. fgndcube_lobe.ppf 1000 subpknolssnocor.ppf 0.35 50. 0. 3.

### Step 5 / Check the results using spiapp 
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
delobjs *
openppf fgndpk.ppf 
openppf fgndpklobe.ppf 
openppf fgndpkcorlobe.ppf 
openppf lsspk.ppf 
openppf lsspklobe.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 3*3*3*360*360*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'