Index: /trunk/Cosmo/RadioBeam/README
===================================================================
--- /trunk/Cosmo/RadioBeam/README	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/README	(revision 3160)
@@ -0,0 +1,22 @@
+===============================================================================
+  Projet BAORadio / HSHS
+  Simulation program for radio lobe reconstruction 
+ 
+  2006-2007  R. Ansari (UPS/LAL-IN2P3)    C. Magneville (DAPNIA/SPP)
+              (C) IN2P3/CNRS              (C)  DAPNIA/CEA
+ ===============================================================================
+
+Note :  Uses SOPHYA class library  (http://www.sophya.org )
+
+List of files and classes :
+---------------------------
+
+- treccyl.cc : main program (for tests)
+- multicyl.h .cc : Class MultiCylindres , 3D (freq,AngX,AngY) signal simulation 
+  and reconstruction. Uses class MultiBeamCyl
+- mbeamcyl.h .cc : Class MultiBeamCyl 2D (freq, AngX) signal simulation and reconstruction
+- brsrc.h .cc : Class BRSourceGen for managing list of sources 
+  (frequencies and angular position)
+- makefile 
+
+- recbeam.pic : script with example commands for display with piapp
Index: /trunk/Cosmo/RadioBeam/brsrc.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/brsrc.cc	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/brsrc.cc	(revision 3160)
@@ -0,0 +1,138 @@
+#include "brsrc.h"
+#include "srandgen.h"
+#include "fioarr.h"
+#include "vector3d.h"
+
+//=================================================
+BRSourceGen::BRSourceGen()
+{
+  int ns = 3*9*5;
+  freq.SetSize(ns);
+  amp.SetSize(ns);
+  angX.SetSize(ns);
+  angY.SetSize(ns);
+  phase.SetSize(ns);
+
+  amp = 1.;
+  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);
+
+  double fr[3] = {0.1, 0.25, 0.4};   // attention 0 <= f < 0.5
+  int k = 0;
+  for(int kf=0; kf<3; kf++) {
+    for(int ja=-4; ja<=4; ja++) {
+      for(int ia=-2; ia<=2; ia++) {
+	freq(k) = fr[kf];
+	// Pour lambda=50 cm, angle M_PI/12. -> delta phi = 2 pi sur ~ 1.9 m
+	angX(k) = M_PI/12.*(double)ja; 
+	// Pour lambda=50 cm, angle M_PI/5000. -> delta phi = 2 pi sur ~ 800 m
+	angY(k) = M_PI/3000.*(double)ia;
+	k++;
+      }
+    }
+  }
+}
+
+BRSourceGen::BRSourceGen(int ns, double maxangXrad, double maxangYrad, 
+			 double minfreq, double minamp, double maxamp)
+{
+  freq.SetSize(ns);
+  amp.SetSize(ns);
+  angX.SetSize(ns);
+  angY.SetSize(ns);
+  phase.SetSize(ns);
+
+  if (minfreq>=0.5) {
+    cout << "BRSourceGen::BRSourceGen()/Warning minfreq=" << minfreq << " >=0.5 --> =0." << endl;
+    minfreq = 0.;
+  }
+  for(int i=0; i<ns; i++) {
+    freq(i) = minfreq+drand01()*(0.5-minfreq);
+    amp(i) = minamp+drand01()*(maxamp-minamp);
+    angX(i) = drandpm1()*maxangXrad;
+    angY(i) = drandpm1()*maxangYrad;
+  }
+  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);
+}
+
+BRSourceGen::BRSourceGen(vector<double> f, int nsf, double maxangXrad, double maxangYrad, 
+			 double minfreq, double minamp, double maxamp)
+{
+  int ns = nsf*f.size();
+
+  freq.SetSize(ns);
+  amp.SetSize(ns);
+  angX.SetSize(ns);
+  angY.SetSize(ns);
+  phase.SetSize(ns);
+  if (minfreq>=0.5) {
+    cout << "BRSourceGen::BRSourceGen()/Warning minfreq=" << minfreq << " >=0.5 --> =0." << endl;
+    minfreq = 0.;
+  }
+
+  int j = 0;
+  for(int k=0; k<f.size(); k++) {
+    if (f[k] >= 0.5) {
+      cout << "BRSourceGen::BRSourceGen() f[k=" << k << "]=" << f[k] << " >=0.5 ignored " << endl;
+      continue;
+    }
+    for(int i=0; i<nsf; i++) {
+      freq(j) = f[k];
+      amp(j) = minamp+drand01()*(maxamp-minamp);
+      angX(j) = drandpm1()*maxangXrad;
+      angY(j) = drandpm1()*maxangYrad;
+      j++;
+    }
+  }
+  phase = RandomSequence(RandomSequence::Flat, 0., M_PI);
+}
+
+BRSourceGen::BRSourceGen(string ppfname)
+{
+  ReadPPF(ppfname);
+}
+
+void BRSourceGen::WritePPF(string ppfname)
+{
+  POutPersist po(ppfname);
+  po << freq << amp << angX << angY << phase;
+}
+
+void BRSourceGen::ReadPPF(string ppfname)
+{
+  PInPersist pi(ppfname);
+  pi >> freq >> amp >> angX >> angY >> phase;
+}
+
+void BRSourceGen::Print(ostream& os)
+{
+  os << "---- BRSourceGen::Print() - NSources= " << freq.Size() << endl; 
+  for(int k=0; k<freq.Size(); k++) {
+    os << k << " - f=" << freq(k) << " amp=" << amp(k) << " angX=" << angX(k) 
+       << " rad =" << Angle(angX(k)).ToDegree() << " deg angY="
+       << " rad =" << Angle(angX(k)).ToArcMin() << " arcmin" 
+       << " phase=" << phase(k) << endl; 
+  }
+}
+
+
+NTuple BRSourceGen::Convert2Table(double freq0)
+{
+  char* nomsntsrc[9] = {"num","freq","freqred","amp","angrad","angdeg","phase",
+			"angyrad","angymin"};
+  NTuple nt(9, nomsntsrc);
+  double xnt[15];
+  for(int kk=0; kk<freq.Size(); kk++) {
+    xnt[0] = kk;
+    xnt[1] = freq(kk)+freq0;
+    xnt[2] = freq(kk);
+    xnt[3] = amp(kk);
+    xnt[4] = angX(kk);
+    xnt[5] = Angle(angX(kk)).ToDegree();
+    xnt[6] = phase(kk);
+    xnt[7] = angY(kk);
+    xnt[8] = Angle(angY(kk)).ToArcMin();
+    nt.Fill(xnt);
+  }
+
+  return nt;
+}
Index: /trunk/Cosmo/RadioBeam/brsrc.h
===================================================================
--- /trunk/Cosmo/RadioBeam/brsrc.h	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/brsrc.h	(revision 3160)
@@ -0,0 +1,62 @@
+#ifndef BRSOURCE_H
+#define BRSOURCE_H
+
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include "tvector.h"
+#include "ntuple.h"
+
+/* 
+   Projet BAORadio / HSHS
+   Programme de simulation pour reconstruction de lobe radio.
+   Reconstruction de lobe a partir de plusieurs cylindres
+
+   R. Ansari - LAL      Sep 2006
+
+*/
+
+//--------------------------------
+// Classe auxiliaire pour generer une liste de sources
+class BRSourceGen {
+public:
+  // Pour fabriquer une liste de ns=3*9*5=135 sources, repartis en 3 frequences (0.1,0.25,0.4)
+  // regulierement repartis en angle X : 9 angX en +/- pas de Pi/12
+  //                           angle Y : 5 angY en +/- pas de Pi/3000
+  BRSourceGen(); // fabrique un pattern fixe de 20 sources
+
+  // Generation d'une liste de source aleatoirement reparties en freq,angX,angY
+  // @ lambda ~ 0.5 m et D = 50 m , limite de diffraction ~ pi/300
+  BRSourceGen(int ns, double maxangXrad=M_PI/3., double maxangYrad=M_PI/500.,
+	      double minfreq=0.05, double minamp=0.3, double maxamp=1.5);
+
+  // Generation de nsf sources pour chaque frequence f[k], reparties aleatoirement en angX,Y
+  BRSourceGen(vector<double> f, int nsf, 
+	      double maxangXrad=M_PI/3., double maxangYrad=M_PI/500.,
+	      double minfreq=0.05, double minamp=0.3, double maxamp=1.5);
+
+  // Pour relire les sources depuis un fichier
+  BRSourceGen(string ppfname);
+
+  inline int NbSources() { return freq.Size(); }
+
+  void WritePPF(string ppfname);
+  void ReadPPF(string ppfname);
+
+  void Print(ostream& os);
+
+  // Retourne la liste des sources sous forme de NTuple
+  // freq0= frequence de base qui s'ajoute aux freq[k]
+  NTuple Convert2Table(double freq0=0.);   
+
+  //-------- Variables / objets membres 
+ Vector freq; // frequence des sources
+ Vector amp;  // amplitude des sources
+ Vector angX;  // direction des sources en radian, selon l'axe X=N-S : kx=sin(ang)
+ Vector angY; // direction des sources selon Y=E-O : ky=sin(angY)
+
+ Vector phase;  // Phase des sources (choisi aleatoirement pour le moment) 
+};
+
+
+
+#endif
Index: /trunk/Cosmo/RadioBeam/makefile
===================================================================
--- /trunk/Cosmo/RadioBeam/makefile	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/makefile	(revision 3160)
@@ -0,0 +1,31 @@
+include $(SOPHYABASE)/include/sophyamake.inc
+
+all : treccyl
+
+clean :
+	rm Objs/*
+depend :
+	mkdir Objs
+
+treccyl : Objs/treccyl
+	echo 'makefile : treccyl made'
+
+######
+
+Objs/treccyl : Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o  Objs/brsrc.o
+	$(CXXLINK) -o Objs/treccyl \
+	Objs/treccyl.o Objs/multicyl.o Objs/mbeamcyl.o Objs/brsrc.o \
+	$(SOPHYAEXTSLBLIST)
+
+Objs/treccyl.o : treccyl.cc  multicyl.h mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/treccyl.o  treccyl.cc
+
+Objs/brsrc.o : brsrc.cc  brsrc.h 
+	$(CXXCOMPILE) -o Objs/brsrc.o  brsrc.cc
+
+Objs/mbeamcyl.o : mbeamcyl.cc  mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/mbeamcyl.o  mbeamcyl.cc
+
+Objs/multicyl.o : multicyl.cc  multicyl.h mbeamcyl.h brsrc.h 
+	$(CXXCOMPILE) -o Objs/multicyl.o  multicyl.cc
+
Index: /trunk/Cosmo/RadioBeam/mbeamcyl.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/mbeamcyl.cc	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/mbeamcyl.cc	(revision 3160)
@@ -0,0 +1,197 @@
+#include "mbeamcyl.h"
+#include "fftpserver.h"
+#include "vector3d.h"
+#include "matharr.h"
+#include "srandgen.h"
+
+//=================================================
+
+MultiBeamCyl::MultiBeamCyl(int nr, int ns, double posy)
+  : texact(ns) , tjitt(ns) , toffset(nr) ,
+    signal(ns), sigjitt(ns) , gain(nr)
+{
+  NR = nr;
+  NS = ns;
+  posY = posy;
+
+  SetPrintLevel(0);
+
+  SetBaseFreqDa(2., 0.25);
+  SetNoiseSigma(0.);
+  SetTimeJitter(0.);
+  SetTimeOffsetSigma(0.);
+  SetGains(1., 0., 0);
+  adfg = false; src = NULL;
+  SetSources(new BRSourceGen, true);
+}
+
+MultiBeamCyl::~MultiBeamCyl()
+{
+  if (adfg && src) delete src;
+}
+
+void MultiBeamCyl::SetSources(BRSourceGen* brs, bool ad)
+{
+  if (brs == NULL) return;
+  if (adfg && src) delete src;
+  src = brs;  adfg=ad;
+  if (PrtLev > 1) 
+    cout << " MultiBeamCyl::SetSources(brs=" << brs <<" ,bool) NbSrc=" 
+	 << src->NbSources() << endl;
+
+}
+
+void MultiBeamCyl::SetGains(double g, double sigg, int nzerogain)
+{
+  if (sigg < 1.e-6)  gain = g;
+  else gain = RandomSequence(RandomSequence::Gaussian, g, sigg);
+  int k;
+  for (k=0; k<NR; k++) if (gain(k) < 0) gain(k) = 0.;
+  for(k=0; k<nzerogain; k++) {
+    int zg = random()%NR;
+    if ((zg >=0) && (zg < NR))  gain(zg) = 0.; 
+  }
+  if (PrtLev > 1) 
+    cout << " MultiBeamCyl::SetGains(g=" << g <<" ,sigg=" << sigg 
+	 << " ,nzg=" << nzerogain << " )" << endl;
+}
+
+void MultiBeamCyl::ComputeTimeVectors()
+{
+  texact = RegularSequence(0., 1.);
+  toffset = RandomSequence(RandomSequence::Gaussian, 0., toffsig);
+  NewTJitVector();
+}
+
+void MultiBeamCyl::NewTJitVector(int num)
+{
+  if (timejitter > 1.e-19) {
+    tjitt = RandomSequence(RandomSequence::Gaussian, 0., timejitter);
+    tjitt += texact;  
+  }
+  else tjitt = texact;  
+  if (num >= 0) tjitt += toffset(num);
+  
+}
+
+int MultiBeamCyl::ComputeSignalVector(int num, bool fgsignojit)
+{
+  int nok = 0;
+  signal = 0.;
+  sigjitt = 0.;
+  for(int is=0; is<src->freq.Size(); is++) {
+    double fr = src->freq(is);
+    if ((fr < 0.) || (fr > 0.5)) continue;
+    nok++;
+    // Pour le dephasage entre recepteurs, on doit utiliser la frequence vraie,
+    // pas celle apres shift (freq-reduite)
+    // lambda = c T = c/freq avec c = 1, dephasage = 2*pi*num*Da*sin(ang)/lambda
+    double dephasage = num*Da*sin(src->angX(is)) + posY*sin(src->angX(is)) ;
+    dephasage *= (2*M_PI*(fr+freq0));
+    // On ajoute alors la phase propre de chaque source
+    dephasage += src->phase(is);
+    double amprep = src->amp(is)*AngResponse(src->angX(is), src->angY(is));
+    for(int k=0; k<NS; k++) {
+      sigjitt(k) += amprep*sin(2.*M_PI*fr*tjitt(k)+dephasage);
+      if (fgsignojit) 
+	signal(k) += amprep*sin(2.*M_PI*fr*texact(k)+dephasage);
+    }
+  }
+  // Application du gain du detecteur 
+  r_4 ga = gain(num);
+  if (fabs(ga-1.) > 1.e-9) {
+    signal *= ga;
+    sigjitt *= ga;
+  }
+  // Ajout de bruit (ampli ...) 
+  if (signoise > 1.e-18) {
+    for(int k=0; k<NS; k++) {
+      sigjitt(k) += GauRnd(0., signoise);
+      if (fgsignojit) signal(k) += GauRnd(0., signoise);
+    }
+  }
+
+  FFTPackServer ffts;
+
+  ffts.FFTForward(sigjitt, f_sigjit);
+  if (fgsignojit) ffts.FFTForward(signal, f_sig);
+
+  return nok;
+}
+
+
+/*  --- a supprimer ?
+inline float myZmodule(complex<r_4>& z) 
+{
+  return (float)sqrt((double)(z.real()*z.real()+z.imag()*z.imag()));
+}
+----- */
+
+void MultiBeamCyl::ReconstructSourcePlane(bool fgzerocentre)
+{
+  ComputeTimeVectors();
+  int noksrc = ComputeSignalVector(0, false);
+  vector<TVector< complex<r_4> > >  vvfc;
+  cout << "MultiBeamCyl::ReconstructSourcePlane()/Info NR=" << NR << " NFreq=" 
+       << f_sigjit.Size()-1 << " NOkSrc=" << noksrc << endl;
+  if (PrtLev > 0) {
+    cout << " ... posY= " << posY << " MeanGain=" << Mean(gain) << " MeanAmpSrc=" 
+	 << Mean(src->amp) << endl;
+    cout << " ...SigNoise=" << signoise << " TimeJitter=" << timejitter 
+	 << " TOffSig=" << toffsig << " Da=" << Da << " Freq0=" << freq0 << endl;
+  }
+  // On ne s'occupe pas de la composante continue
+  for(int jf=1; jf<f_sigjit.Size(); jf++) {
+    TVector<complex<r_4> > cf(NR);
+    vvfc.push_back(cf);
+  }
+  cout << "ReconstructSourcePlane()/Info: computing s(t) for each receptor ..." << endl;
+  int pmod = NR/10;
+  for(int ir=0; ir<NR; ir++) {
+    if (timejitter > 1.e-19) NewTJitVector(ir);
+    noksrc = ComputeSignalVector(ir, false);
+    for(int jf=1; jf<f_sigjit.Size(); jf++)
+      vvfc[jf-1](ir) = f_sigjit(jf);
+    if ( (PrtLev>0) && (ir%pmod == 0) )
+      cout << " OK s(t) for ir=" << ir << " / NR=" << NR << " NOkSrc=" << noksrc << endl;
+  }
+
+  cout << "ReconstructSourcePlane()/Info: computing s(ang) for each freq by FFT" << endl;
+  cmplx_srcplane.SetSize(f_sigjit.Size()-1, NR);
+  // rec_srcplane.SetSize(f_sigjit.Size()-1, NR);
+  FFTPackServer ffts;
+  TVector<complex<r_4> > fcf;
+  pmod = vvfc.size()/10;
+  for(int jf=0; jf<vvfc.size(); jf++) {
+    ffts.FFTForward(vvfc[jf], fcf);
+    if (fcf.Size() != NR) {
+      cout << "ReconstructSourcePlane()/BUG jf=" << jf << " fcf.Size() != NR "
+	   << fcf.Size() << " != " << NR << endl;
+      continue;
+    }
+    if (fgzerocentre) {  // On veut avoir la direction angle=0 au milieu de la matrice
+      int milieu = (NR-1)/2;
+      if (NR%2 == 0) milieu++;
+      int decal = NR - milieu - 1;
+      for (int ir=0; ir<=milieu; ir++) 
+	cmplx_srcplane(jf, ir+decal) = fcf(ir);
+      //	rec_srcplane(jf, ir+decal) = myZmodule(fcf(ir));
+
+      for (int ir=milieu+1; ir<NR; ir++) 
+	cmplx_srcplane(jf, decal-(NR-ir)) = fcf(ir);
+      //	rec_srcplane(jf, decal-(NR-ir)) = myZmodule(fcf(ir));
+      
+    }
+    else for (int ir=0; ir<NR; ir++) 
+      cmplx_srcplane(jf, ir) = fcf(ir);
+    //      rec_srcplane(jf, ir) = myZmodule(fcf(ir));
+
+    if ( (PrtLev > 0) && (jf%pmod == 0)) 
+      cout << " OK rec_srcplane(jf, ir) for jf=" << jf << endl;
+  }
+
+  
+  cout << "ReconstructSourcePlane()/Info: rec_srcplane computed OK" << endl;
+}
+
+
Index: /trunk/Cosmo/RadioBeam/mbeamcyl.h
===================================================================
--- /trunk/Cosmo/RadioBeam/mbeamcyl.h	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/mbeamcyl.h	(revision 3160)
@@ -0,0 +1,114 @@
+#ifndef MBEAMCYL_H
+#define MBEAMCYL_H
+
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include "tvector.h"
+#include "brsrc.h"
+
+/* 
+   Classe de simulation simple pour estimer les effets de bruit 
+   (ampli, ...) et decalages, bruit de phase d'horloge pour la 
+   reconstruction de la position des sources (manip BAO-radio)
+   R. Ansari - LAL      Sept 2006
+
+   Note : On travaille ds un systeme sans dimension
+   unite de temps = 1 = temps d'echantillonnage (f_echant = 1)
+   vitesse de propagation (c) = 1
+   Donc, lambda=1  pour f = 1
+
+   On specifie les frequences des sources (BRSourceGen) en frequence reduite
+   dans la premiere zone de Nyquist (f = 0 ... 0.5) 
+   Neanmoins, on a besoin de la frequence vrai des sources --> methode SetBaseFreqDa
+*/
+
+
+//-----------------------------------------------------------------------
+//--- Classe de reconstruction de plan source correspondant a un cylindre
+class MultiBeamCyl {
+public :
+  // nr = nb de recepteurs (antennes) le long du cylindre
+  // ns = nb d'echantillons en temps de chaque paquet
+  // posy = position spatiale en Y pour chaque cylindre (unite ou lambda=1 pour f=1)
+ MultiBeamCyl(int nr=1024, int ns=4096, double posy=0.);
+ ~MultiBeamCyl();
+
+ // Niveau de print de debug
+ inline void SetPrintLevel(int prl=0) { PrtLev=prl; }
+
+ // Reponse angulaire pour un recepteur (cylindre+antenne)
+ // Modifier cette fonction pour tenir compte de la response angulaire
+ inline double AngResponse(double angx, double angy) { return 1.; }
+
+ inline double getCylinderYPos() { return posY; }
+ // Specification de la frequence de base f0 et espacement des recepteurs
+ // freq_vrai = freq_BRSourceGen + f0
+ // pour f_vrai = 2, lambda = 0.5 et lambda/2 = 0.25
+ inline void SetBaseFreqDa(double f0=2., double a=0.25) 
+   {  freq0 = f0;   Da = a; }
+
+ // frequences reduites (entre 0 ... 0.5) , angle en radian, amp ~ 1
+ void SetSources(BRSourceGen* brs, bool ad=true);
+ inline void SetSources(BRSourceGen& brs)
+   { SetSources(&brs, false); }
+
+ // Definition du sigma du bruit gaussien sur les echantillons
+ inline void SetNoiseSigma(double sig=0.) {  signoise = sig; }
+ // Definition du sigma du jitter d'horloge (typ 0.01)
+ inline void SetTimeJitter(double tjit=0.) { timejitter = tjit; } 
+ // Definition du sigma des offsets d'horloge entre recepteurs (typ 0.02)
+ inline void SetTimeOffsetSigma(double tsig=0.) { toffsig = tsig; }
+
+ // Definition du gain et sigma de fluctuations de gain 
+ // nzerogain: nb de recepteurs (choisis au hasard) avec gain mis a zero
+ void SetGains(double g=1., double sigg=0., int nzerogain=0);
+
+ // Calcule les vecteurs de temps (0 1 2 ...ns-1)
+ void ComputeTimeVectors();
+ // Recalcule un nouveau vecteur temps - avec jitter
+ void NewTJitVector(int num=-1);
+
+ // Calcul du vecteur signal (s(t)) -  et les vecteurs TF 
+ // pour le recepteur numero numdet
+ // si fgsignojit==true, on calcule aussi pour Temps_Sans_Jitter
+ int ComputeSignalVector(int numdet=0, bool fgsignojit=true);
+
+ // Calcul de la matrice reconstruction du plan source (freq vs angle)
+ // fgzerocentre -> angle zero au mileiu de la matrice
+ void ReconstructSourcePlane(bool fgzerocentre=true);
+
+ inline  TMatrix< complex<r_4> > & getRecSrcPlane() { return cmplx_srcplane; }
+
+ //-------------- Variables membres ------------ 
+ int NR, NS;  // nb recepteurs, nb d'echantillons ds chaque paquet
+ double posY;    // Position selon Y (E-O) de la ligne de recepteurs
+
+ int PrtLev;   // Niveau de print de debug
+
+ double Da; // distance entre recepteurs
+ double freq0;  // frequence de base (freqvrai=f+freq0)
+ 
+ Vector texact;  // les temps exact
+ Vector tjitt;   // temps avec jitter
+ Vector toffset;  // Offset en temps entre recepteurs 
+ double timejitter; // le sigma du jitter en temps des coups d'horloge
+ double toffsig;    // sigma des offsets en temps
+ 
+ BRSourceGen* src;  // Les sources
+ bool adfg; // if true, delete src 
+
+ double signoise; // sigma du bruit additif (bruit ampli ... )
+ 
+ TVector<r_4> gain;  // Gain de chaque recepteur 
+
+ TVector<r_4> signal; // signal sans jitter en temps 
+ TVector<r_4> sigjitt; // signal avec jitter en temps
+
+ TVector< complex<r_4> > f_sig, f_sigjit;  // TF des vecteurs signal , sigjitt
+
+ TMatrix< complex<r_4> > cmplx_srcplane; // composantes complexe du plan-source reconstruit
+ // TMatrix<r_4> rec_srcplane;  // Matrice plan source colonnes->angX, lignes->freq 
+ };
+
+
+#endif
Index: /trunk/Cosmo/RadioBeam/multicyl.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/multicyl.cc	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/multicyl.cc	(revision 3160)
@@ -0,0 +1,154 @@
+#include "multicyl.h"
+#include "fftpserver.h"
+#include "vector3d.h"
+#include "matharr.h"
+#include "srandgen.h"
+#include "ctimer.h"
+#include "resusage.h"
+
+
+//=================================================
+
+MultiCylinders::MultiCylinders(int nr, int ns)
+{
+  NR = nr;
+  NS = ns;
+
+  SetPrintLevel(0);
+
+  SetBaseFreqDa(2., 0.25);
+  SetNoiseSigma(0.);
+  SetTimeJitter(0.);
+  SetTimeOffsetSigma(0.);
+  SetGains(1., 0., 0);
+  adfg = false; src = NULL;
+  SetSources(new BRSourceGen, true);
+}
+
+MultiCylinders::~MultiCylinders()
+{
+  if (adfg && src) delete src;
+  for(int k=0; k<mCyl.size(); k++) delete mCyl[k];
+}
+
+MultiBeamCyl& MultiCylinders::GetCylinder(int k)
+{
+  if ((k < 0) || (k >= mCyl.size()))
+    throw RangeCheckError("MultiCylinders::GetCylinder(k) k<0 OR k>=NCyl");
+  return (*mCyl[k]);
+}
+
+void MultiCylinders::SetSources(BRSourceGen* brs, bool ad)
+{
+  if (brs == NULL) return;
+  if (adfg && src) delete src;
+  src = brs;  adfg=ad;
+}
+
+
+void MultiCylinders::ConfigureCylinders()
+{
+  cout << " MultiCylinders::ConfigureCylinders() with NCyl= " << mCyl.size() << endl;
+  for(int k=0; k<mCyl.size(); k++) {
+    mCyl[k]->SetPrintLevel(PrtLev);
+    mCyl[k]->SetBaseFreqDa(freq0, Da);
+    mCyl[k]->SetNoiseSigma(signoise);
+    mCyl[k]->SetTimeJitter(timejitter);
+    mCyl[k]->SetTimeOffsetSigma(toffsig);
+    mCyl[k]->SetGains(gain, siggain, ngainzero);
+    mCyl[k]->SetSources(src, false);
+  }
+}
+
+
+void MultiCylinders::ReconstructCylinderPlaneS(bool fgzerocentre)
+{
+  Timer tm("RecCylPlaneS");
+  ResourceUsage resu;
+  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - NCyl= " << mCyl.size() 
+       << " MemSize=" << resu.getMemorySize() << " kb" << endl;       
+  ConfigureCylinders();
+
+  char buff[128];
+  for(int k=0; k<mCyl.size(); k++) {
+    cout << "---- Cyl[" << k << "] Calling MultiBeamCyl.ReconstructSourcePlane() ..." << endl;
+    mCyl[k]->ReconstructSourcePlane(fgzerocentre);
+    sprintf(buff,"Cyl[%d].RecSrcPlane()",k);
+    tm.Split(buff);
+  }
+  resu.Update();
+  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - Done " 
+       << " MemSize=" << resu.getMemorySize() << " kb" << endl;
+}
+
+void MultiCylinders::ReconstructSourceBox(int halfny, double stepangy)
+{
+  Timer tm("RecSrcBox");
+
+  ReconstructCylinderPlaneS(true);
+  TMatrix< complex<r_4> > & mtx = GetCylinder(0).getRecSrcPlane(); 
+  // boite 3D X:angX, Y:angY , Z: freq
+  sa_size_t sz[5] = {0,0,0,0,0};
+  sz[0] = mtx.NCols();
+  sz[1] = halfny*2+1;
+  sz[2] = mtx.NRows();
+  TArray< complex<r_4> > & box = getRecSrcBox();
+  box.ReSize(3, sz);
+  //  box = complex< r_4 >( 0.f, 0.f );
+  cout << " MultiCylinders::ReconstructSourceBox(" << halfny << "," << stepangy
+       << "=" << Angle(stepangy).ToArcMin() << " ) srcbox:" << endl;
+  box.Show(cout);
+
+
+  int pmod = mtx.NRows()/10;
+  
+  double fstep = 1./(double)NS;  // pas en frequence, attention, on a vire la composante continu
+  for(int kf=0; kf<mtx.NRows(); kf++) { // Loop over frequencies
+    double frq = (double)(kf+1.)*fstep + freq0;  // + 1 car f=0 a ete vire
+    
+    for(int lc=0; lc<mCyl.size(); lc++) {  // Loop over cylinders 
+      MultiBeamCyl& cyl = GetCylinder(lc);
+
+      TMatrix< complex<r_4> > & recp = cyl.getRecSrcPlane(); 
+
+      double facl = - 2*M_PI*frq*cyl.getCylinderYPos();  // attention au signe -
+      double dphi;
+      complex< r_4 > phasefactor;
+      int jyy = 0;
+      for(int jy=-halfny; jy<=halfny; jy++) {  // Loop over Y angular steps
+	double dphi = facl * sin( (double)jy*stepangy );
+	phasefactor = complex< r_4 > ((r_4) cos(dphi) , (r_4) sin(dphi));
+	for(int ix=0; ix<sz[0]; ix++) {  // Loop over AngX directions
+	  // sur recp : index ligne -> frequence  ,    index colonne -> angX ,
+	  box(ix, jyy, kf) += recp(kf, ix)*phasefactor;
+	}  //  
+	jyy++;
+      } // End of Loop over Y angles
+    } // End of loop over cylinders
+    tm.Nop();
+    if ( (PrtLev>0) && (kf%pmod == 0) )   
+      cout << " OK box(angx,angy, freq=kf) done for kf=" << kf 
+	   << " / Max_kf=" << mtx.NRows() << endl;
+  } // End of loop over over frequencies
+
+  cout << " MultiCylinders::ReconstructSourceBox() done " << endl; 
+} 
+
+inline float myZmodule(complex<r_4>& z) 
+{
+  return (float)sqrt((double)(z.real()*z.real()+z.imag()*z.imag()));
+}
+
+TMatrix< r_4 >  MultiCylinders::getRecXYSlice(int kfmin, int kfmax)
+{
+  TArray< complex<r_4> > & box = getRecSrcBox();  
+  if ((kfmin < 0) || (kfmin >=  box.SizeZ()) || (kfmax < kfmin) || (kfmax >=  box.SizeZ()) )
+    throw RangeCheckError("MultiCylinders::getRecXYSlice(kfmin, kfmax)");
+  TMatrix< r_4> rmtx(box.SizeY(), box.SizeX());
+  for(int kf=kfmin; kf<=kfmax; kf++) {
+    for(int jy=0; jy<box.SizeY(); jy++) 
+      for(int ix=0; ix<box.SizeX(); ix++) 
+	rmtx(jy, ix) += myZmodule(box(ix, jy, kf));
+  }
+  return(rmtx);
+}
Index: /trunk/Cosmo/RadioBeam/multicyl.h
===================================================================
--- /trunk/Cosmo/RadioBeam/multicyl.h	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/multicyl.h	(revision 3160)
@@ -0,0 +1,107 @@
+#ifndef MULTICYL_H
+#define MULTICYL_H
+
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include "mbeamcyl.h"
+
+/* 
+   Projet BAORadio / HSHS
+   Programme de simulation pour reconstruction de lobe radio.
+   Reconstruction de lobe a partir de plusieurs cylindres
+
+   R. Ansari - LAL      Jan 2007
+
+   Note : On travaille ds un systeme sans dimension
+   unite de temps = 1 = temps d'echantillonnage (f_echant = 1)
+   vitesse de propagation (c) = 1
+*/
+
+
+//-----------------------------------------------------------------------
+//--- Classe de reconstruction de plan source correspondant a un cylindre
+class MultiCylinders {
+public :
+  // nr = nb de recepteurs
+  // ns = nb d'echantillons en temps de chaque paquet
+ MultiCylinders(int nr=256, int ns=1024);
+ ~MultiCylinders();
+
+ // Niveau de print de debug
+ inline void SetPrintLevel(int prl=0) { PrtLev=prl; }
+
+ // Ajout d'un cylindre, en position posY
+ inline int AddCylinder(double posy)
+   {  mCyl.push_back( new MultiBeamCyl(NR, NS, posy) ); return mCyl.size(); }
+
+ MultiBeamCyl & GetCylinder(int i);
+
+ // Specification de la frequence de base f0 et espacement des recepteurs
+ inline void SetBaseFreqDa(double f0=2., double a=0.25) 
+   {  freq0 = f0;   Da = a;  }
+
+ // frequences reduites (entre 0 ... 0.5) , angle en radian, amp ~ 1
+ void SetSources(BRSourceGen* brs, bool ad=true);
+ inline void SetSources(BRSourceGen& brs)
+   { SetSources(&brs, false); }
+
+ // Definition du sigma du bruit gaussien sur les echantillons
+ inline void SetNoiseSigma(double sig=0.) {  signoise = sig; }
+ // Definition du sigma du jitter d'horloge (typ 0.01)
+ inline void SetTimeJitter(double tjit=0.) { timejitter = tjit; } 
+ // Definition du sigma des offsets d'horloge entre recepteurs (typ 0.02)
+ inline void SetTimeOffsetSigma(double tsig=0.) { toffsig = tsig; }
+
+ // Definition du gain et sigma de fluctuations de gain 
+ // nzerogain: nb de recepteurs (choisis au hasard) avec gain mis a zero
+ inline void SetGains(double g=1., double sigg=0., int nzerogain=0)
+   {  gain=g;  siggain=sigg;  ngainzero = nzerogain; }
+
+
+ // Fait la configuration des cylindres, et reconstruit le plan source
+ // pour chaque cylindre
+ void ReconstructCylinderPlaneS(bool fgzerocentre=true);
+
+ // Reconstruction de la boite 3D des sources 
+ // - halfny : nb de bin en angY = 2*halfny+1 
+ // - stepangy : pas en angle Y (radian) 1 arcmin= 0.00029 rad (~Pi/11000)
+ // Cette methode appelle ReconstructCylinderPlaneS(true);
+ void ReconstructSourceBox(int halfny=15, double stepangy=M_PI/10000);
+
+
+ // Pour recuperer une tranche reconstruite angX, angY 
+ // avec la moyenne des modules entre kfmin <= k (selon z) <= kfmax 
+ TMatrix< r_4 > getRecXYSlice(int kfmin, int kfmax);
+
+ // Acces a la boite 3D de sources reconstruite
+ inline TArray< complex<r_4> > & getRecSrcBox() { return cmplx_srcbox; }
+
+ // Configure chaque cylindre a partir des parametres de la classe 
+ // est appele automatiquement par ReconstructCylinderPlaneS
+ void ConfigureCylinders();
+
+ //-------------- Variables membres ------------ 
+ int NR, NS;  // nb recepteurs, nb d'echantillons ds chaque paquet
+
+ int PrtLev;   // Niveau de print de debug
+
+ double Da; // distance entre recepteurs
+ double freq0;  // frequence de base (freqvrai=f+freq0)
+ 
+ double timejitter; // le sigma du jitter en temps des coups d'horloge
+ double toffsig;    // sigma des offsets en temps
+ 
+ BRSourceGen* src;  // Les sources
+ bool adfg; // if true, delete src 
+
+ double signoise; // sigma du bruit additif (bruit ampli ... )
+ double gain, siggain;
+ int ngainzero;
+
+ vector<MultiBeamCyl *> mCyl;  // Les differents cylindres 
+
+ TArray< complex<r_4> > cmplx_srcbox; // boite3D des sources (angX,Y,freq) reconstruit
+ };
+
+
+#endif
Index: /trunk/Cosmo/RadioBeam/recbeam.pic
===================================================================
--- /trunk/Cosmo/RadioBeam/recbeam.pic	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/recbeam.pic	(revision 3160)
@@ -0,0 +1,50 @@
+### Script piapp pour faire tourner le programme de simulation 
+### de reconstruction de beam par la technique de double FFT
+###  Reza , Sep 2006
+#  Le code (2 classes) se trouve dans mbeamcyl.h mbeamcyl.cc
+#  et une classe auxiliaire brsrc.h brsrc.cc 
+#  Classe  MultiBeamCyl et la classe utilitaire BRSourceGen
+# --------------------
+# 0/ Verification, modif des parametres
+#  Modifier si besoin le code dans trecbeam.cc qui instancie
+#  les classes (MultiBeamCyl en particulier) avec le parametres
+#  necessaires (nb de recepteurs, nb d'echantillons ...) 
+#  Avant de faire la simulation et la reconstruction de beam
+
+# 1/ Faire tourner le programme - reconstruction plan X-freq (en dehors de piapp)
+Objs/treccyl X onecyl.ppf
+
+
+# 2/ Affichage resultat
+#  Afficher le plan source reconstruit (les frequences le 
+# long de l'axe Y (les lignes de la matrice), et les directions
+# des sources (angles) suivant l'axe Ox (colonnes)
+# Le zero des angles se trouve au milieu
+delobjs *
+openppf onecyl.ppf
+disp recsrcplane 'lut=sqrt,0.,0.1 invy grey32 win'
+
+# 2.b/ le ntuple ntsrc contient la liste des sources de depart
+print ntsrc
+# Attention, ce qui est reconstruit n'est pas directement l'angle, mais
+# sin(angle)*frequence - pour une matrice 512lignesx256colonnes
+plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'
+addtitle 'PlanSource (Y=freq,X=AngNS)' ' ' 'font=times,bold,18 red'
+w2eps recbeam.eps
+
+
+# 3/ Faire tourner le programme - reconstruction cube X-Y-freq (en dehors de piapp)
+Objs/treccyl XY fourcyl.ppf
+
+# 4/ Affichage resultats
+delobjs *
+openppf fourcyl.ppf
+listobjs
+
+disp recsrcplane0 'lut=sqrt,0.,0.1 invy grey32 win'
+plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'
+disp recsrcplane3 'lut=sqrt,0.,0.1 invy grey32 win'
+plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'
+
+
+
Index: /trunk/Cosmo/RadioBeam/treccyl.cc
===================================================================
--- /trunk/Cosmo/RadioBeam/treccyl.cc	(revision 3160)
+++ /trunk/Cosmo/RadioBeam/treccyl.cc	(revision 3160)
@@ -0,0 +1,227 @@
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <math.h>
+#include <iostream>
+#include <typeinfo>
+
+#include "tvector.h"
+#include "srandgen.h"
+#include "fioarr.h"
+#include "sopemtx.h"
+#include "pexceptions.h"
+#include "matharr.h"
+
+#include "sambainit.h"
+
+// #include "tarrinit.h"
+
+#include "timing.h"
+
+#include "multicyl.h"
+#include "mbeamcyl.h"
+
+
+// Declaration des fonctions de ce fichier
+static int test1cyl(string& ppfname);
+static int testmulticyl(string& ppfname);
+
+//-----------------------------------------------------------
+// -------------- Parametres de simulation  -----------------
+//-----------------------------------------------------------
+static int MR = 256;  // Nombre de recepteur
+static int NE = 1024;  // Nombre d'echantillon en temps;
+static double freq0 = 2.;  // frequence de base
+static double da = 0.25;     // pas des antennes le long du cylindre
+static double snoise = 1.0;  // sigma du bruit 
+static double tjit = 0.05;   // sigma du jitter en temps
+static double tos = 0.02;    // sigma des offsets en temps
+static double gmean = 1.;    // gain moyen
+static double gsig = 0.;     // sigma des gains
+static int nantgz = 0;       // nb d'antennes morts (-> gain=0)
+static int prtlevel = 0;     // niveau de print 
+//-----------------------------------------------------------
+
+
+/* --------------------------------------------------------
+  Le main programme de test des classes de reconstruction
+  multilobe radio - R. Ansari , Sep06 -- 2007
+  --------------------------------------------------------- */
+   
+int main(int narg, char* arg[])
+{
+
+SophyaInit();
+InitTim();   // Initializing the CPU timer
+
+string ppfname = "treccyl.ppf";
+int act = 1;
+if (narg < 3) {
+  cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
+       << "  act= X ou XY  , ppfname=  treccyl.ppf par defaut" << endl;
+  return 1;
+}
+if (strcmp(arg[1],"XY") == 0) act = 2;
+if (narg > 2)  ppfname = arg[2];
+if (narg > 3)  prtlevel = atoi(arg[3]);
+
+int rc = 0;
+try {
+  if (act == 2) rc = testmulticyl(ppfname);
+  else rc = test1cyl(ppfname);
+  cout << "treccy/Info - FIN  " << endl;
+}
+  catch (PThrowable& exc) {
+    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
+    rc = 77;
+  }  
+  catch (std::exception& sex) {
+    cerr << "\n treccyl.cc std::exception :" 
+         << (string)typeid(sex).name() << "\n msg= " 
+         << sex.what() << endl;
+  }
+  catch (...) {
+    cerr << " treccyl.cc catched unknown (...) exception  " << endl; 
+    rc = 78; 
+  } 
+
+ return rc;
+}
+
+
+//--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
+int test1cyl(string& ppfname)
+{
+
+  // BRSourceGen sg;
+  int nsrc = 60; 
+  BRSourceGen sg(nsrc);
+  //  sg.WritePPF(string("brsrc1.ppf"));
+
+  cout << "=== test1cyl: BRSourceGen NbSrc= " << sg.NbSources()
+       << " NbRecep=" << MR << " NSamples=" << NE << endl;
+
+  // BRSourceGen  sg(string("brsrc1.ppf"));
+  if (prtlevel > 1)  sg.Print(cout);
+  
+  
+  MultiBeamCyl  mb(MR, NE);
+  mb.SetPrintLevel(prtlevel);
+  mb.SetBaseFreqDa(freq0, da);
+  mb.SetNoiseSigma(snoise);
+  mb.SetTimeJitter(tjit);
+  mb.SetTimeOffsetSigma(tos);
+  mb.SetGains(gmean, gsig, nantgz);
+
+  mb.SetSources(sg);
+  
+  mb.ComputeTimeVectors();
+  mb.ComputeSignalVector(0, true);
+  cout << "treccy/test1cyl:  signal vectors OK " << endl;
+  PrtTim("test1cyl:[1] ");
+
+  POutPersist po(ppfname);
+  po << PPFNameTag("signal") << mb.signal;
+  po << PPFNameTag("sigjitt") << mb.sigjitt;
+  po << PPFNameTag("f_sig") << mb.f_sig;
+  po << PPFNameTag("f_sigjit") << mb.f_sigjit;
+
+  NTuple ntsrc = sg.Convert2Table(freq0);
+  po << PPFNameTag("ntsrc") << ntsrc;
+
+  cout << "treccy/test1cyl: - sig/f_sig,ntsrc to OutPPF OK  " << endl;
+  PrtTim("test1cyl[2] ");
+  
+  mb.ReconstructSourcePlane(true);
+  {
+  TMatrix<r_4> srcplane = module(mb.getRecSrcPlane() );
+  po << PPFNameTag("recsrcplane") << srcplane;
+  }
+  PrtTim("test1cyl[3] ");
+
+  return 0;
+
+}
+
+
+//--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
+int testmulticyl(string& ppfname)
+{
+
+  // BRSourceGen sg;
+  int nsf = 10;
+  vector<double> frq;
+  frq.push_back(0.1);
+  frq.push_back(0.27);
+  frq.push_back(0.38);
+  
+  
+  cout << "treccy/testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf << ")" << endl;
+  BRSourceGen sg(frq, nsf, M_PI/3, M_PI/1000);
+  // sg.WritePPF(string("brsrcm.ppf"));
+  // BRSourceGen  sg(string("brsrcm.ppf"));
+  cout << "=== testmulticyl: BRSourceGen NbSrc= " << sg.NbSources() 
+       << " NbRecep=" << MR << " NSamples=" << NE << endl;
+  if (prtlevel > 1)  sg.Print(cout);
+  
+  
+  MultiCylinders  mcyl (MR, NE);
+  mcyl.SetPrintLevel(prtlevel);
+  mcyl.SetBaseFreqDa(freq0, da);
+  mcyl.SetNoiseSigma(snoise);
+  mcyl.SetTimeJitter(tjit);
+  mcyl.SetTimeOffsetSigma(tos);
+  mcyl.SetGains(gmean, gsig, nantgz);
+
+  mcyl.AddCylinder(0.);
+  mcyl.AddCylinder(100.);
+  mcyl.AddCylinder(300.);
+  mcyl.AddCylinder(700.);
+  mcyl.AddCylinder(1500.);
+
+  mcyl.SetSources(sg);
+  
+  PrtTim("testmulticyl[1] ");
+
+  //  mcyl.ReconstructCylinderPlaneS(true);
+  mcyl.ReconstructSourceBox();
+
+  POutPersist po(ppfname);
+
+  NTuple ntsrc = sg.Convert2Table(freq0);
+  po << PPFNameTag("ntsrc") << ntsrc;
+
+  {
+  TMatrix<r_4> srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane());
+  po << PPFNameTag("recsrcplane0") << srcplane0;
+  }
+
+  {
+  TMatrix<r_4> srcplane3 = module(mcyl.GetCylinder(3).getRecSrcPlane());
+  po << PPFNameTag("recsrcplane3") << srcplane3;
+  }
+
+  PrtTim("testmulticyl[2] ");
+
+  int kfmin, kfmax;
+  po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();
+  kfmin = mcyl.getRecSrcBox().SizeZ()*frq[0] - 3;  kfmax = kfmin+6;
+  {
+  TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
+  po << PPFNameTag("recXYf0") << slice0;
+  }
+  kfmin = mcyl.getRecSrcBox().SizeZ()*frq[1] - 3;  kfmax = kfmin+6;
+  {
+  TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
+  po << PPFNameTag("recXYf1") << slice1;
+  }
+  kfmin = mcyl.getRecSrcBox().SizeZ()*frq[2] - 3;  kfmax = kfmin+6;
+  {
+  TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
+  po << PPFNameTag("recXYf2") << slice2;
+  }
+
+  PrtTim("testmulticyl[3] ");
+
+  return 0;
+
+}
