Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3488)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3489)
@@ -1,9 +1,13 @@
 ===============================================================================
   Projet BAORadio / HSHS
-  Simulation program for radio lobe reconstruction 
+  Simulation programs for radio lobe reconstruction, noise estimate etc ... 
  
-  2006-2007  R. Ansari (UPS/LAL-IN2P3)    C. Magneville (DAPNIA/SPP)
+  2006-2008  R. Ansari (UPS/LAL-IN2P3)    C. Magneville (DAPNIA/SPP)
               (C) IN2P3/CNRS              (C)  DAPNIA/CEA
  ===============================================================================
+
+
+A/ 3D reconstruction (alpha,delta,freq) simulation program 
+===========================================================
 
 Note :  Uses SOPHYA class library  (http://www.sophya.org )
@@ -21,2 +25,17 @@
 
 - recbeam.pic : script with example commands for display with piapp
+
+
+B/ P(k) in temperature plot + "cosmic variance" + noise ...
+=============================================================
+Note : to be executed in spiapp 
+List of Files : anapkn.pic + cmvhpkz.ppf
+
+C/ Interferometer configuration / noise estimate 
+=================================================
+Note : to be executed in spiapp 
+
+C.1/ formnoise.pic + vnpairs.cc + vnpairs.h 
+
+C.2/ pa.pic + pa.cc + pa.h : Lobe synthesis (dipoles sum)
+
Index: /trunk/Cosmo/RadioBeam/anapkn.pic
===================================================================
--- /trunk/Cosmo/RadioBeam/anapkn.pic	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/anapkn.pic	(revision 3489)
@@ -0,0 +1,379 @@
+######################################################################
+#### 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
+fHI = 0.020/0.1 
+
+####################################################################
+## 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.57*$cgeom*$fHI
+  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 
+##  VOL1: 10deg x 10deg x Delta_z=0.2 : 335 x 335 x 680 Mpc^3 = 76 10^6 Mpc^3
+##  VOL2: 20deg x 20deg x Delta_z=0.2 : VOL2 = 300 10^6 Mpc^3
+##  VOL3: 60deg x 60deg x Delta_z=0.2 : VOL4 = 2700 10^6 Mpc^3
+##  VOL4: 40deg x 40deg x Delta_z=0.2 : VOL3 = 1200 10^6 Mpc^3
+set angcov ( '10x10 deg^2' '20x20 deg^2' '60x60 deg^2' '40x40 deg^2' )
+VOL1 = 76e6
+VOL2 = 300e6
+VOL3 = 2700e6
+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
+
+
+
+## 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 VOl1 en ~ 7 passes (7 couvertures elementaires)
+#  120 jours pour VOl2 en ~ 4x7 passes (28 couvertures elementaires)
+#  120 jours pour VOl3 en ~ 16x7 passes (225 couvertures elementaires)
+set ndays ( '120 days' '120 days'  '120 days' '120 days' ) 
+TINT1 = 1.5e6
+TINT2 = 3.7e5
+TINT3 = 4.5e4
+TINT4 = 1.e5
+
+#  Bruit Tsys = 100 K = 10^5 mK 
+TSYS = 1e5
+
+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 ( green blue red )
+
+#############################################################
+### 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 ( green blue red )
+  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'
+
+endscript
Index: /trunk/Cosmo/RadioBeam/formnoise.pic
===================================================================
--- /trunk/Cosmo/RadioBeam/formnoise.pic	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/formnoise.pic	(revision 3489)
@@ -0,0 +1,81 @@
+#######################################################
+##### Script de trace de PNoise pour differentes configurations 
+####### BAORadio --  R. Ansari - Avril 2008 
+### ===> Fichiers formnoise.pic + vnpairs.cc + vnpairs.h
+
+## --- Script d'execution 
+defscript go 
+  delobjs *
+  c++execfrf vnpairs.cc vnpairs.h
+endscript 
+## ---- Script prjkt H2NPairsName ProfName Nbin 
+defscript prjkt 'Usage: prjkt H2NPairsName ProfName Nbin' 
+  if ( $# < 3 ) then
+     echo 'Usage: prjkt H2NPairsName ProfName Nbin'
+     return
+  endif 
+  set h2n $1
+  set prn k_theta_$2
+  set nbin $3
+  # Longueur d'onde a z ~ 1 
+  LAMBDA = 0.42
+  # Facteur de distance transverse a z ~ 1 (Mpc)
+  DCT = 3400.
+  newprof $prn 0. 3000 $nbin 
+  projprof $prn $h2n 2*M_PI*sqrt(x*x+y*y)/$LAMBDA val 
+  set prn k_Mpc_$2
+  newprof $prn 0. 1. $nbin 
+  projprof $prn $h2n 2*M_PI*sqrt(x*x+y*y)/$LAMBDA/$DCT val 
+  disp $prn 
+endscript 
+
+## --- Script de representation de config antennes
+defscript drantennes
+  zone 1 
+  nt2d recp xr yr - - - - 'red marker=circle,13 nsta '
+  settitle 'Configuration antennes - D=5m' ' ' 'font=times,bolditalic,18'
+  setaxelabels 'x (metres)' 'y (metres)' 'font=times,bolditalic,18' 
+endscript 
+
+
+defscript toto 
+setaxesatt 'font=times,bold,18 fixedfontsize minorticks'
+# go
+zone 1 3
+prjkt h2_fa fa 100 
+prjkt h2_fa_d2m fa_d2m 200 
+prjkt h2hex hex 200 
+
+# Config A : Packed 2 meter dishes 100m x 100m   TPose ~ 6.25 10^4 sec  (6.25 = (5/2)^2)
+# Config B : Packed 5 meter dishes 100m x 100m   TPose ~ 10^4 sec
+# Config C : 5 meter dishes - config hexa, cercle R ~ 100 m   TPose ~ 10^4 sec
+# Config D : 5 meter dishes - config hexa, cercle R ~ 100 m  TPose ~ 50 10^4 sec
+
+# Nb de paires A , B , C 
+set npairs ( 2450 380 122 6000 )
+set kMpc ( k_Mpc_fa_d2m k_Mpc_fa k_Mpc_hex k_Mpc_hex ) 
+set maxprof ( ${k_Mpc_fa_d2m.vmax}  ${k_Mpc_fa.vmax} ${k_Mpc_hex.vmax} ${k_Mpc_hex.vmax} )  
+
+echo $maxprof
+set cols ( blue magenta red orange black )
+zone 1
+for i 0:4
+  echo "----- I= $i " 
+  ph = "$kMpc[i]"
+  maxph = $maxprof[i]
+  npair = $npairs[i]
+  set gropt "nsta same connectpoints line=solid,2 $cols[i]" 
+  if ( $i == 0 ) then 
+    set gropt "xylimits=0.005,1.,1,5000 logx logy nsta connectpoints line=solid,2 $cols[i] notit" 
+  endif
+  plot2d $ph x 1e4/val/$npair*$maxph val>5 "$gropt" 
+end
+settitle 'PNoise (mK^2) Mpc^2  z~1' ' ' 'font=times,bolditalic,18'
+setaxelabels 'k_Transverse Mpc^-1 ' ' PNoise (mK^2 Mpc^2)'  'font=times,bolditalic,18'
+
+addline 0.04 1000 0.4 20 'line=dashed,2 black' 
+set txt ( '(A) 100mx100m Packed 2m dishes, tint=6.2 10^4 s' '(B) 100mx100m Packed 5m dishes, tint=10^4 s' '(C) 400x5m dishes over 200m circle, , tint=10^4s' "(C') 400x5m dishes over 200m circle, , tint=50 10^4s" 'P(k) HI' ) 
+textdrawer txt cols 'font=times,bold,16 frame'
+
+endscript  
+ 
Index: /trunk/Cosmo/RadioBeam/pa.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/pa.cc	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/pa.cc	(revision 3489)
@@ -0,0 +1,75 @@
+char * names[10] = {"m","teta0","teta","wre","wim","wmod2","wx","wy","lwx","lwy"};
+NTuple lobes(10, names);
+int N = 7;
+int NO2 = N/2;
+int NTET=1000;
+// double dtet0 = Lambda/(NO2*DELX);
+
+/* #define LobeW1D LobeSquare 
+#define LobeW1D LobeGauss  
+double dtet = M_PI/30./NTET;
+double dtet0 = M_PI/60./NO2; 
+*/
+
+/*  */
+#define LobeW1D LobeDipoleDemiLambda 
+double dtet = M_PI/2./NTET;
+double dtet0 = M_PI/3./NO2; 
+DELX = Lambda/2.;
+/*  */
+
+cout << " DTeta= " << dtet << " DTeta0= " << dtet0 << endl;
+cout << " LobeW1D(0.) = " << LobeW1D(0.) << endl; 
+Timer tm("LOBES");
+
+Vector lobe1d(2*NTET+1),vteta(2*NTET+1);
+TVector< complex<r_8> >  lobeS(2*NTET+1);
+
+for(int i=-NTET; i<=NTET; i++) {
+  double teta = i*dtet;
+  lobe1d(i+NTET) = LobeW1D(teta);
+  vteta(i+NTET) = teta;
+  complex<r_8> ws(0., 0.);
+  for(int n=-NO2; n<=NO2; n++) {
+  	double phase = 2.*M_PI*DELX/Lambda*n*teta;	
+  	ws += complex< r_8 >( cos(phase), sin(phase) )*LobeW1D(teta);
+   }
+  lobeS(i+NTET) = ws; 
+}
+
+KeepObj(lobe1d);
+KeepObj(vteta);
+KeepObj(lobeS);
+
+
+double xnt[20];
+for(int m=-NO2; m<=NO2; m++) {
+  if (m%10 == 0) tm.Split();	
+  xnt[0] = m;  xnt[1] = m*dtet0;
+  for(int i=-NTET; i<=NTET; i++) {
+    double teta = i*dtet;
+    xnt[2] = teta;
+    complex<r_8> ws(0., 0.);
+    complex<r_8> wsck(0., 0.);
+    for(int n=-NO2; n<=NO2; n++) {
+      double phase = 2.*M_PI*DELX/Lambda*n*(teta-m*dtet0);	
+      ws += complex< r_8 >( cos(phase), sin(phase) )*LobeW1D(teta);
+    }
+    xnt[3] = ws.real();
+    xnt[4] = ws.imag();
+    //    xnt[5] = sqrt(ws.real()*ws.real()+ws.imag()*ws.imag());
+    xnt[5] = ws.real()*ws.real()+ws.imag()*ws.imag();
+    xnt[6] = xnt[5]*cos(teta); 
+    xnt[7] = xnt[5]*sin(teta); 
+    if (xnt[5] > 1.e-9) {
+      xnt[8] = 10.*log10(xnt[5])*cos(teta); 
+      xnt[9] = 10.*log10(xnt[5])*sin(teta); 
+    }
+    else xnt[8] = xnt[9] = -90.;
+
+    lobes.Fill(xnt);
+  }		
+}	
+
+cout << lobes;
+KeepObj(lobes);
Index: /trunk/Cosmo/RadioBeam/pa.h
===================================================================
--- /trunk/Cosmo/RadioBeam/pa.h	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/pa.h	(revision 3489)
@@ -0,0 +1,47 @@
+static double D = 4.;  // Taille du reflecteur (en m)
+static double Lambda = 0.4;  // Longueur d'onde en m 
+static double PiDoL = M_PI*D/Lambda;  // Pi * D / Lambda
+
+static double DELX = D; // Separation des reflecteurs
+
+// static double DELX = Lambda/2.; // Separation des reflecteurs
+
+
+double LobeSquare(double tet) 
+{
+if (fabs(tet) < 1.e-9) return 1.;
+double x = tet*PiDoL;
+return sin(x)/x;
+}
+
+double LobeGauss(double tet) 
+{
+double x = tet*PiDoL;
+ return exp(-x*x/10.);
+}
+
+
+double LobeDipoleDemiLambda(double tet) 
+{
+ double L = 0.5; 
+ double trad = acos(sin(tet));
+ double pil = M_PI*L;
+ double st = sin(trad), ct = cos(trad);
+ // L'antenne standing-wave
+ double v = 0.;
+ if(fabs(st)<1e-3) {
+   if(ct>0.) {  // theta ~= 0
+     v = pil*sin(pil)*trad/2.;  /* v *= v; */
+   } else {  // theta ~= 180
+     v = pil*sin(pil)*(trad-M_PI)/2.;  /* v *= v; */
+   }
+ } else {
+   v = (cos(pil*ct) - cos(pil))/st;  /* v *= v; */
+ }
+
+ //  cout << " teta= " << tet << " trad= " << trad << " v=" << v << endl;
+ return v;
+
+}
+
+
Index: /trunk/Cosmo/RadioBeam/pa.pic
===================================================================
--- /trunk/Cosmo/RadioBeam/pa.pic	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/pa.pic	(revision 3489)
@@ -0,0 +1,50 @@
+### Calcul et trace de lobes synthetises 
+### BAORadio --  R. Ansari  -   Avril 2008 
+### ===> Fichiers pa.pic + pa.cc + pa.h
+defscript go 
+  delobjs *
+  c++execfrf pa.cc pa.h 
+endscript
+
+go
+
+setaxesatt 'font=times,bolditalic,18 fixedfontsize minorticks'
+
+defscript tracelobe
+  set gropt "connectpoints xylimits=-90,90,50,280000 notit nsta"
+  if ( $# > 0 ) then
+    gropt = "$gropt $1" 
+  endif
+  plot2d lobes teta*180/M_PI wmod2 m==0 "$gropt" 
+  plot2d lobes teta*180/M_PI wmod2 m==50 'connectpoints same blue nsta'
+  plot2d lobes teta*180/M_PI wmod2 m==150 'connectpoints same orange nsta'
+  plot2d lobes teta*180/M_PI wmod2 m==200 'connectpoints same red nsta'
+  setaxelabels ' Angle Theta (degree) ' 'Synth. lobes' 'font=times,bold,18' 
+  settitle 'Reconstructed lobes (FFT-like synthesis)' ' '  'font=times,bold,18'
+  set txt ( 'Theta=0 deg' 'Theta=+12 deg' 'Theta=+36 deg' 'Theta=48 deg' )
+  set cols ( 'black' 'blue' 'orange' 'red' )
+  textdrawer txt cols 'font=times,bolditalic,18 frame' 
+endscript 
+
+
+defscript tracelobe2
+  plot2d lobes teta*180/M_PI wmod2 m==0 "black connectpoints notit nsta" 
+  plot2d lobes teta*180/M_PI wmod2 m==1 "blue same connectpoints notit nsta" 
+  plot2d lobes teta*180/M_PI wmod2 m==2 "magenta same connectpoints notit nsta" 
+  plot2d lobes teta*180/M_PI wmod2 m==3 "red same connectpoints notit nsta" 
+  setaxelabels ' Angle Theta (degree) ' 'Lobe' 'font=times,bold,18' 
+  settitle 'Reconstructed lobes - 7 Lambda/2 dipole sum' ' '  'font=times,bold,18'
+endscript 
+
+defscript tracelobe3
+  plot2d lobes wy wx m==0 'black connectpoints'
+  plot2d lobes wy wx m==2 'red same connectpoints'
+  plot2d lobes wy wx m==-2 'red same connectpoints'
+  plot2d lobes wy wx m==1 'orange same connectpoints'
+  plot2d lobes wy wx m==-1 'orange same connectpoints'
+  plot2d lobes wy wx m==3 'blue same connectpoints'
+  plot2d lobes wy wx m==-3 'blue same connectpoints'
+  plot2d lobes wy wx m==4 'violet same connectpoints'
+  plot2d lobes wy wx m==-4 'violet same connectpoints'
+endscript 
+
Index: /trunk/Cosmo/RadioBeam/vnpairs.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/vnpairs.cc	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/vnpairs.cc	(revision 3489)
@@ -0,0 +1,144 @@
+#define NP 10
+#define ND 20
+
+#define DIAMETRE 5 
+
+// Timer tm;
+
+Vector lesx(ND);
+
+/* int spos[NP] = {1,2,3,5,6,7,10,12,14,16};
+for(int i=0; i<NP; i++)  lesx(i) = -spos[NP-i-1];  
+for(int i=NP; i<2*NP; i++)  lesx(i) = spos[i-NP];  
+*/
+int spos[ND] = {-15,-13,-12,-10,-7,-6,-5,-3,-2,-1,1,2,3,5,6,7,10,12,14,15};
+for(int i=0; i<ND; i++)  lesx(i) = spos[i];  
+
+
+Matrix pos(2,ND*ND); 
+
+int k = 0;
+for(int i=0; i<ND; i++)  {
+  for(int j=0; j<ND; j++) {
+  	pos(0,k) = lesx(i); pos(1,k) = lesx(j);  k++;
+  }
+} 
+
+
+Histo h(-0.5, 31.5, 32);
+
+for(int i=0; i<ND-1; i++) 
+  for(int j=i+1; j<ND; j++) 
+    h.Add(fabs(lesx(i)-lesx(j)),1.);
+
+KeepObj(h);
+KeepObj(lesx);
+
+
+Histo2D h2(-0.5, 31.5, 32, -0.5, 31.5, 32);
+for(int i=0; i<ND*ND-1; i++) 
+  for(int j=i+1; j<ND*ND; j++) 
+    h2.Add(fabs(pos(0,i)-pos(0,j)),fabs(pos(1,i)-pos(1,j)), 1.);
+
+KeepObj(h2);
+
+
+Histo2D h2f(-0.5, 31.5, 96, -0.5, 31.5, 96);
+for(int i=0; i<ND*ND-1; i++) 
+  for(int j=i+1; j<ND*ND; j++) {
+    FillVis(h2f, fabs(pos(0,i)-pos(0,j)),fabs(pos(1,i)-pos(1,j)));
+  }
+KeepObj(h2f);
+
+cout << "h2f OK" << endl;
+
+//---------------------------------------------------------------------
+//----- Disposition compacte d'antennes circulaires -------------------
+//---------------------------------------------------------------------
+{
+double DIAM_FA = 5.;
+double L_FA = 100.;
+double DISTA_FA = DIAM_FA;
+
+int ND1_FA = L_FA/DISTA_FA;
+int ND_FA = ND1_FA*ND1_FA;
+
+Matrix pos_fa(2,ND_FA);
+int kk=0;
+for(int i=0; i<ND1_FA; i++) {
+  for(int j=0; j<ND1_FA; j++) {
+  	pos_fa(0, kk) = i*DISTA_FA;  pos_fa(1, kk) = j*DISTA_FA;  kk++;
+  }
+}
+KeepObj(pos_fa);
+cout << "pos_fa OK" << endl;
+
+Histo2D h2_fa(0., 2*ND1_FA*DIAM_FA, 3*ND1_FA, 0., 2*ND1_FA*DIAM_FA, 3*ND1_FA);
+Histo2D h2ns_fa(0., 2*ND1_FA*DIAM_FA, 4*ND1_FA, 0., 2*ND1_FA*DIAM_FA, 4*ND1_FA);
+
+for(int i=0; i<pos_fa.NCols()-1; i++) 
+  for(int j=i+1; j<pos_fa.NCols(); j++) {
+  	h2ns_fa.Add(fabs(pos_fa(0,i)-pos_fa(0,j)),fabs(pos_fa(1,i)-pos_fa(1,j)),1.);
+    FillVis(h2_fa, fabs(pos_fa(0,i)-pos_fa(0,j)),fabs(pos_fa(1,i)-pos_fa(1,j)), DIAM_FA);
+  }
+KeepObj(h2_fa);
+KeepObj(h2ns_fa);
+cout << "h2_fa OK" << endl;
+
+}
+//---------------------------------------------------------------------
+//------------------ Disposition hexagonale ---------------------------
+//---------------------------------------------------------------------
+{
+int nlayer = 4;
+int nhex = 1;
+for(int i=1; i<=nlayer; i++) {
+  nhex += i*6;	
+}
+int ndet = nhex*7;
+cout << " NLayer= " << nlayer << " NHex= " << nhex << " NDet= " << ndet << endl;
+ 
+Matrix posh(2,ndet);
+
+char * names[3] = {"xr", "yr", "sz"};
+NTuple recp(3, names);
+double xnt[5];
+xnt[0] = xnt[1] = xnt[2] = xnt[3] = xnt[4] = 0.;
+
+int kk = 0;
+for(int i=0; i<= nlayer; i++) {
+  double rlayer = i*4.5*DIAMETRE;
+//  if (i == 1) rlayer = i*4.*DIAMETRE;
+  Matrix ph = couche(i,rlayer) ; 
+  cout << " DEBUG - I=" << i << " / NLayer= " << nlayer << " ph.NCols= " << ph.NCols() << endl;
+  for(int j=0; j<ph.NCols(); j++) {
+  	Matrix hex = hexagon(ph(0,j), ph(1,j), DIAMETRE, ph(2,i));
+  	  for(int k=0; k<hex.NCols(); k++) {
+  	  	xnt[0] = posh(0, kk) = hex(0, k); 
+  		xnt[1] = posh(1, kk) = hex(1, k); 
+  		recp.Fill(xnt);
+	    kk++;
+  	  }  
+  }	
+}
+
+KeepObj(posh);
+KeepObj(recp);
+
+cout << "posh/recp OK" << endl;
+
+Histo2D h2hex(0., 40.*DIAMETRE, 120, 0., 40.*DIAMETRE, 120);
+Histo2D h2nshex(0., 20.*DIAMETRE, 120, 0., 20.*DIAMETRE, 120);
+
+for(int i=0; i<posh.NCols()-1; i++) 
+  for(int j=i+1; j<posh.NCols(); j++) {
+  	h2nshex.Add(fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j)), 1.);
+    FillVis(h2hex, fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j)), DIAMETRE);
+    FillVis(h2hex, fabs(posh(0,i)-posh(0,j))*0.75,fabs(posh(1,i)-posh(1,j)), DIAMETRE);
+    FillVis(h2hex, fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j))*0.75, DIAMETRE);
+  }
+KeepObj(h2hex);
+KeepObj(h2nshex);
+cout << "h2hex OK" << endl;
+}
+
Index: /trunk/Cosmo/RadioBeam/vnpairs.h
===================================================================
--- /trunk/Cosmo/RadioBeam/vnpairs.h	(revision 3489)
+++ /trunk/Cosmo/RadioBeam/vnpairs.h	(revision 3489)
@@ -0,0 +1,58 @@
+Matrix hexagon(double xc, double yc, double r=1., double theta=0.) 
+{
+	Matrix hex(2,7);
+	hex(0,0) = xc;  hex(1,0) = yc;
+	for(int i=1; i<=6; i++) {
+	  double tet = (i-1)*M_PI/3.+theta;
+	  hex(0,i) = r*cos(tet)+xc;
+	  hex(1,i) = r*sin(tet)+yc;
+	}
+	return hex;
+}
+
+Matrix couche(int num, double r=0.)
+{
+	if (num <= 0) {
+      Matrix pos(3,1);
+	  pos(0,0) = pos(1,0) = pos(2,0) = 0.;
+	  return pos;
+	}
+	else {
+	  if (r <= 0.0001)  r = 5.*num;  	
+	  int n = num*6;
+	  double tet0 = 0.;
+	  tet0 = M_PI/n;
+	  double dtet = M_PI*2./(double)n;
+	  Matrix pos(3, n);
+  	  for(int i=0; i<n; i++) {
+	    double tet = i*dtet+tet0;
+	    pos(0,i) = r*cos(tet);
+	    pos(1,i) = r*sin(tet);
+	    pos(2,i) = tet;
+	  }
+	  return pos;
+	}
+}
+
+int FillVis(Histo2D& h2, double x, double y, double d=1) 
+{
+  h2.Add(x, y, 1.);
+
+  double wa = 0.6;  
+  double wb = 0.3;  
+  for(int i=0; i<6; i++) {
+  	double dx1a = 0.35*cos(i*M_PI/3.)*d;
+  	double dy1a = 0.35*sin(i*M_PI/3.)*d;
+  	double dx1b = 2*dx1a; 
+  	double dy1b = 2*dy1a;
+  	for(int j=0; j<6; j++) {
+  	  double dx2a = 0.35*cos(j*M_PI/3.)*d;
+  	  double dy2a = 0.35*sin(j*M_PI/3.)*d;
+  	  double dx2b = 2*dx2a; 
+  	  double dy2b = 2*dy2a;
+  	  h2.Add((x+dx1a-dx2a), (y+dy1a-dy2a), wa); 
+  	  h2.Add((x+dx1b-dx2b), (y+dy1b-dy2b), wb);  	
+  	}
+  }
+  return 0;
+}
