######################################################################
#### Script de trace de P(k) en temperature , effet de bruit, etc ####
###     Fev - Avril 2008 ,  BAORadio/Reza
###  Tout est dans ce script piapp + fichier de donnees PPF cmvhpkz.ppf
## > lancer spiapp 
## > exec anapkn.pic 
## > Faire copier/coller des commandes a la fin, ds script xxtoto
######################################################################

##1 / On ouvre le fichier extrait d'un PPF fait par SimLSS (cmv) 
# delobjs *
openppf cmvhpkz.ppf
set gdz 1.
set k x
set kk pow(10.,x)

setaxesatt 'font=times,bold,16 fixedfontsize'



#  --- Cosmologie WMAP , O_L ~ 0.7 O_m ~ 0.25  H0 = 0.7 x100 km/s/Mpc
Om = 0.25
OL = 0.7 
#  ====> On se met a z = 0.5 , lambda ~ 0.315 , nu ~ 946 MHz  (spectre P(k) du fichier cmvobserv3d.ppf)
z = 0.5
# On considere une fraction de HI (HI/Baryon) = 0.02
gHI = 0.02/0.01

####################################################################
## Script pour calcul du coefficient de conversion P(k)-masse en P(k) temperature en mK^2
defscript convpk2t21 '  --> calcul le coefficient de conversion P(k) en P(k)_temperature mK^2'
  cgeom = (1+$z)*(1+$z)/sqrt($OL+pow((1+$z),3.)*$Om)
  cct21 = 0.054*$cgeom*$gHI
  cct21 = $cct21*$cct21
  echo "----  cct21 (NoFactCroiss)= $cct21"
#  On approxime le facteur de croissance lineaire au carre par 
#  [(1+0.5)/(1+z)]^2  --> erreur ~ 10%
#  cct21 = $cct21*pow(1.5/(1+$z),2)
  echo "---- Conv P(k)->Temperature mK^2 : cct21= $cct21"
endscript

# On execute ce script pour calculer cct21 a z=0.5
convpk2t21

# A z = 0.5 , DA = 1280 Mpc  (ComovDA = 1920 Mpc) , H ~ 88 km/s/Mpc
# A z = 1 , DA = 1700 Mpc  (ComovDA = 3400 Mpc) , H ~ 116 km/s/Mpc
#  10 arcmin -> l=5.58 Mpc , k=2 pi/l= 1.1
#  1 deg = 60 arcmin -> l=55.8 Mpc , k=2 pi/l= 0.11
#  5 deg = 300 arcmin -> l=279 Mpc , k ~ 0.02
Hz05 = 88 
H = 88
DAz05 = 1920 
DA = 1920 

##  ---- VOL1,VOL2,VOL3 
##  VOL0: 10deg x 10deg x Delta_z=0.2 : 335 x 335 x 680 Mpc^3 = 76 10^6 Mpc^3
##  VOL1: 20deg x 20deg x Delta_z=0.2 : VOL1 = 300 10^6 Mpc^3
##  VOL2: 360deg x 10deg x Delta_z=0.2 : VOL2 = 2700 10^6 Mpc^3
##  VOL3: 360deg x 30deg x Delta_z=0.2 : VOL3 = 10800 10^6 Mpc^3
##  VOL4: 40deg x 40deg x Delta_z=0.2 : VOL3 = 1200 10^6 Mpc^3
set angcov ( '20x20 deg^2' '360x10 deg^2' '360x30 deg^2' '40x40 deg^2' )
VOL0 = 76e6
VOL1 = 300e6
VOL2 = 2700e6
VOL3 = 10800e6
VOL4 = 1200e6


## ----- Les antennes , environ 15x15 ~ 225-250 antennes elementaires 
#  On considere des reflecteurs de ~ 7-8 metres de diametres, S~40 m^2 
#  Ou 400 antennes de D~5 metres
#  S_totale ~ 10 000 m^2  , reparties sur ~ 200mx200m , filling factor ~ 0.25
#  FOV instantane ~ 2.5 degres de diametre  (pour D~6-7 m)
#  FOV instantane ~ 4.3 x 4.3 degres^2 (pour D~5 m)
#  Resolution ~ 0.31/200 -> ~ 5 arcmin  
#  VolumePixel ~ 3 Mpc x 3 Mpc x 4 Mpc ~ 36 Mpc^3 
vpixz05 = 50
vpix = 50
#  kappa ~ (1/fill_factor)^2
kappa = 25 

#  Setup sur 200mx200m , 10000 m^2
defscript setup1 
  vpixz05 = 36
  vpix = 36
  kappa = 200*200/10000
  kappa = $kappa*$kappa
  echo "setup-1: VPix=$vpix  Kappa=$kappa"
endscript
#  Setup sur 300mx300m , 10000 m^2
defscript setup2 
  vpixz05 = 16
  vpix = 16
  kappa = 300*300/10000
  kappa = $kappa*$kappa
  echo "setup-2: VPix=$vpix  Kappa=$kappa"
endscript

#  Setup sur 400mx400m , 20000 m^2
defscript setup3
  vpixz05 = 16
  vpix = 16
  kappa = 400*400/20000
  kappa = $kappa*$kappa
  echo "setup-3: VPix=$vpix  Kappa=$kappa"
endscript

#  Setup 64 x D=5m sur 100mx100m , 10000 m^2
defscript setup5
  vpixz05 = 100
  vpix = 100
  kappa = 100*100/1200
  kappa = $kappa*$kappa
  echo "setup-5: VPix=$vpix  Kappa=$kappa"
endscript


## delta k pour le calcul du nombre de modes, correspond a ~ l=5m/lambda / DA
delkz05 = 0.01 
delk = 0.01 


#  --- Temps d'integration  a z=0.5 avec FOV=4x4 degres^2  (D=5m) 
#  120 jours pour VOl0 en ~ 7 passes (7 couvertures elementaires)
#  360 jours pour VOl2 en ~ 16x7 passes (225 couvertures elementaires)
#  360 jours pour VOl1 en ~ 4x7 passes (28 couvertures elementaires)
#  360 jours pour VOl2 en ~ 90x3 passes (270 couvertures elementaires)
#  360 jours pour VOl3 en ~ 90x9 passes (810 couvertures elementaires)

set ndays ( '360 days' '360 days'  '360 days' '360 days' ) 
TINT1 = 1.1e6
TINT2 = 1.1e5
TINT3 = 3.8e4
TINT4 = 1.e5

#  Bruit Tsys = 100 K = 10^5 mK 
TSYS = 1e5
#  Bruit Tsys = 50 K = 5.10^4 mK 
TSYS = 5e4

DELNU = 1.e6

#############################################################
### Setup Nancay multi-beam 
defscript setupnancaymb ' --> Setup Nancay multi-beam '
#  4x22 arcmin @ 1420 MHz (21 cm) 
#  6x33 arcmin^2 -> 3.4 x 18.5 Mpc @ z=0.5 (1920 Mpc) 
  vpixz05 = 250
  vpix = 250
#  kappa ~ (1/fill_factor)^2
  kappa = 1
#  FOV ~ 5 deg^2  (10degx30') @ z=0.5 
  FOV = 5 
  echo "setup-Nancay-MB: VPix=$vpix  Kappa=$kappa"
  TSYS = 200.e3
  TINT1 = 4.e5
  TINT2 = 1.e5
  TINT3 = 0.25e5
  TINT4 = 0.6e4
  set ndays ( '60 days' '60 days'  '60 days' '60 days' ) 
  set vols ( $VOL1 $VOL2 $VOL3 )
  set tints ( $TINT1 $TINT2 $TINT3 ) 
endscript

# Definition tableau de volumes, couleurs ...
set vols ( $VOL1 $VOL2 $VOL3 )
set tints ( $TINT1 $TINT2 $TINT3 )
set cols ( red orange blue )

#############################################################
### Script de scaling des quantites avec z 
defscript scalewz ' Scale with z , Usage scalewz z DA H TimeFac'
  if ( $# < 5 ) then 
    echo ' scalewz/error, Usage scalewz z DA H TimeFac'
  endif
  z = $1
  DA = $2
  H = $3
  timfac = $4
  vfac = ($DA/$DAz05)*($DA/$DAz05)*($H/$Hz05)
  zfac = (1+$z)/1.5
  vpix = $vpixz05*$vfac*$zfac*$zfac
  delk = $delkz05*$zfac
  set vols ( $VOL1 $VOL2 $VOL3 )
  set tints ( $TINT1 $TINT2 $TINT3 )
  for i 0:3
    vols[i] = $vols[i]*$vfac
    tints[i] = $tints[i]*$zfac*$zfac*$timfac
  end  
  convpk2t21
  echo "ScaleZ/Info: z=$z DA=$DA VPix=$vpix DelK=$delk"
  echo "ScaleZ/Info VOLS= " $vols
  echo "ScaleZ/Info TINTS= " $tints
endscript 


#############################################################
### Script de calcul de PNOISE
defscript pnoise ' --> Calcul de PNOISE , appel pnoise tint'
  tint = $1
  PNOISE = $kappa*$TSYS*$TSYS*$vpix/($tint*$DELNU)
  echo " --pnoise tint= $tint --> PNOISE= $PNOISE"
endscript  



#############################################################
### Script de calcul de coefficient de conversion P/k en sigma[P(k)]
defscript csigpk ' --> Calcul de Sigma_P(k) = P/k * $CSPK  ( appel csigpk VolSurvey)'
  vsurv = $1
  CSPK = sqrt(4*Pi*Pi/($delk*$vsurv))
  echo " --csigpk V_survey= $vsurv --> CSPK= $CSPK"
endscript  

#############################################################
### Script de calcul de la valeur de k pour Deltax = 10,50,100 m @ z donne
defscript kvalfdx ' --> Calcul de k fonction de x'
  set lesdx ( 5 10 20 50 100 )
  lambdaz = 0.21*(1+$z)
  set lesk ( 1 1 1 1 1 )
  for i 0:$#lesdx
    lesk[i] = 2*Pi*$lesdx[i]/$DA/$lambdaz
  end
  echo "---kvalfdx/ lesdx=" $lesdx 
  echo "---kvalfdx/ lesk=" $lesk 
endscript  
# On l'execute
kvalfdx

#############################################################
###   Trace de P(k) 
defscript showpkt  '  Affichage P(k) et P(k) en Temp mK^2 '
  zone 1 2 
  y1 = 500
  y2 = 9e4
  xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
  n/plot hpkz.val%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2"
  setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k) ' 'font=times,bolditalic,16' 
  settitle 'Mass power spectrum h=0.7, z=0.5' ' '  'font=times,bolditalic,16' 

  y1 = $y1*$cct21
  y2 = $y2*$cct21
  xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
  n/plot hpkz.val*${cct21}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2"
  setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=times,bolditalic,16' 
  settitle 'H_I sky temperature power spectrum (mK^2) , h=0.7, z=0.5' ' '  'font=times,bolditalic,16' 
endscript
#-----------------------------


#############################################################
###    Script de trace de Sigma_P(k)/P(k)
defscript showspk  ' Trace de Sigma_P(k)/P(k) ( appel showspk n=0..3) '
  xyl = "xylimits=0.01,0.5,1.e-3,0.4 logx logy minorticks"
  set vols ( $VOL1 $VOL2 $VOL3 )
  set cols ( red orange blue )
  csigpk $vols[0]
  n/plot hpkz.$CSPK/$kk%$kk ! ! "notit nsta connectpoints $cols[0] $xyl line=solid,2"
  for i 1:$1
    csigpk $vols[i]
    n/plot hpkz.$CSPK/$kk%$kk ! ! "same nsta connectpoints $cols[i] line=solid,2"
  end

  setaxelabels 'k (Mpc^-1)   h=0.7' 'Sigma_P(k)/P(k)  ' 'font=times,bolditalic,16' 
  settitle 'Relative Sample variance (h=0.7)' ' '  'font=times,bolditalic,16' 
  addline 0.01 0.1 0.4 0.1 'black line=dotted,1 same'

  kvalfdx 
  set kcols ( orange turquoise violet gold orange ) 
  for i 0:$#lesk
    addline $lesk[i] 1.e-3 $lesk[i] 0.4 "$kcols[i] line=dashed,1 same"
    tx = $lesk[i]*0.9
    ty = 0.02
    tatt = 'font=times,bolditalic,16 textdirvertup vertcenter horizcenter'
    addtext $tx $ty " Dist= $lesdx[i] meters " "$kcols[i] $tatt"
  end
endscript 

#############################################################
### Script de trace de P(k) +/- Sigma[P(k) 
defscript showpk  ' Trace de P(k) avec +/-sigma ( appel showpk n=0..3) '
  y1 = 500*$cct21
  y2 = 9e4*$cct21
  xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
  n/plot hpkz.val*${cct21}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2"
#  if ( $# > 1 )  then  
    for i 0:$1 
      csigpk $vols[i]
      n/plot hpkz.val*${cct21}*(1+$CSPK/$kk)%$kk ! ! "same nstat connectpoints $cols[i] line=solid,1"
      n/plot hpkz.val*${cct21}*(1-$CSPK/$kk)%$kk ! ! "same nstat connectpoints $cols[i] line=solid,1"
   end
#  endif
  setaxelabels 'k (Mpc^-1)   h=0.7' 'P(k)    (mK^2 Mpc^3)' 'font=times,bolditalic,16' 
  settitle 'H_I sky temperature power spectrum T21-mk' ' '  'font=times,bolditalic,16' 
endscript 


#############################################################
### Script de trace de PNoise superpose a P(k)
defscript shownoise  ' Trace de PNOISE ( appel shownoise n=0..3) '
  x1 = 0.01 
  x2 = 0.5
  for i 0:$1 
    pnoise $tints[i]
    func $PNOISE $x1 $x2 10 "same nstat connectpoints $cols[i] line=dashed,2"
    tx = 0.1
    ty = $PNOISE*1.1
    addtext $tx $ty "PNoise $angcov[i] $ndays[i] " "font=times,bolditalic,16 $cols[i]"
  end
#  addtext 0.1 1000 ' PNoise (dashed)' 'font=times,bolditalic,16 red'
endscript 


defscript doall
  showpkt 
  w2ps
  w2eps pkt.eps
  sleep 2
  showpk 3 
  shownoise 3
  showspk 3
  w2ps
  w2eps pknoise.eps
endscript 


defscript dopk 
  showpk 3 
  shownoise 3
  if ( $# > 0 )  then
    settitle "$1" ' ' 'font=times,bold,16' 
  endif
  showspk 3
  if ( $# > 1 )  then
    settitle "$2" ' ' 'font=times,bold,16' 
  endif
endscript


defscript xxtoto
## Exemple de commandes
newwin 1 2
setup1
scalewz 0.5 1920 88 1
dopk 'P(k), PNoise - mK^2 , Setup-1 z=0.5'  'Relative Sample Variance, z=0.5'
##
newwin 1 2
setup1
scalewz 1 3400 115 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-1 z=1.0'   'Relative Sample Variance, z=1.0'
##
newwin 1 2
setup1
scalewz 1.5 4540 150 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-1 z=1.5'   'Relative Sample Variance, z=1.5'
##
newwin 1 2
setup1
scalewz 2.0 5420 192 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-1 z=2.0'   'Relative Sample Variance, z=2.0'
##
newwin 1 2
setup1
scalewz 2.5 6120 237 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-1 z=2.5'   'Relative Sample Variance, z=2.5'
##
newwin 1 2
setup2
scalewz 0.5 1920 88 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-2 z=0.5'  'Relative Sample Variance, z=0.5'
##
##
newwin 1 2
setup2
scalewz 1 3400 115 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-2 z=1.0' 'Relative Sample Variance, z=1.0'
##
newwin 1 2
setup3
scalewz 1 3400 115 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-3 z=1.0'   'Relative Sample Variance, z=1.0'
##
newwin 1 2
setup3
scalewz 1.5 4500 150 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Setup-3 z=1.5'  'Relative Sample Variance, z=1.5'

### Setup Nancay multi-beam
newwin 1 2
setupnancaymb
scalewz 0.5 1920 88 1
dopk 'P(k), PNoise - mK^2 Mpc^3, Nancay-MultiBeam  z=0.5'  'Relative Sample Variance, z=0.5'

### Setup pour plpkn @ z=0.7
setup5 
scalewz 0.7 2500 100 1
  y1 = 500*$cct21
  y2 = 9e4*$cct21
  xyl = "xylimits=0.01,0.5,$y1,$y2 logx logy minorticks"
  n/plot hpkz.val*${cct21}%$kk ! ! "same notit nsta connectpoints blueviolet $xyl line=solid,2"

endscript