source: Sophya/trunk/Cosmo/RadioBeam/pa.cc

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