#####################################################################################
####  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 effect on foreground cube and LSS cube
csh> ./Objs/applobe fgndcube.ppf fgndcube_lobe.ppf
csh> ./Objs/applobe lsscube.ppf lsscube_lobe.ppf
## Step 3.b/ Correct for the lobe effect by bringing all to a beam of 30 arcmin 
csh> ./Objs/applobe lsscube_lobe.ppf lsscube_corlobe.ppf 1 30
csh> ./Objs/applobe fgndcube_lobe.ppf fgndcube_corlobe.ppf 1 30

### 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 0. 0.
# 4.d / Extract LSS P(k) from Foreground+LSS+noise and beam effect - correcting for beam to 30 arcmin
csh> ./Objs/calcpk2 lsscube_lobe.ppf 0.2 fgndcube_lobe.ppf 1000 subpklobe.ppf 0.35 30. 3.

### Step 5 / Check the results using spiapp 
openppf fgndpk.ppf
openppf lsspk.ppf
openppf lsspkwn.ppf
openppf lsspklobe.ppf