#####################################################################################
####  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/ 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); 

print omap
saveppf omap 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 nvsscube.ppf
# expmeansig syncube val 
# expmeansig nvsscube val
c++exec TArray<r_4> fgndcube = syncube+nvsscube; 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_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.
# 4.d / 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 subpklobe.ppf 0.35 50. 130. 3.

### Step 5 / Check the results using spiapp 
delobjs *
openppf fgndpk.ppf 
disp fgndpk 'logx logy'
# Adjust the limits
openppf fgndpklobe.ppf 
disp fgndpklobe 'same red'
disp fgndpklobe 'same red nsta'
openppf fgndpkcorlobe.ppf 
disp fgndpkcorlobe 'same magenta nsta'
openppf lsspk.ppf 
disp lsspk 'same gold nsta'
openppf lsspklobe.ppf 
disp lsspklobe 'same green nsta'
openppf subpklobe.ppf 
disp subpklobe 'same blue nsta'