Index: trunk/Cosmo/SimLSS/Makefile
===================================================================
--- trunk/Cosmo/SimLSS/Makefile	(revision 3622)
+++ trunk/Cosmo/SimLSS/Makefile	(revision 3632)
@@ -21,5 +21,5 @@
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec $(EXE)cmvdefsurv \
        $(EXE)cmvtintfun $(EXE)cmvconcherr $(EXE)cmvtinterp $(EXE)cmvtstagn \
-       $(EXE)cmvschdist
+       $(EXE)cmvschdist $(EXE)cmvhshsns
  
 PROGSOBJ = \
@@ -28,5 +28,6 @@
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
           $(OBJ)cmvtintfun.o $(OBJ)cmvtinterp.o \
-          $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o $(OBJ)cmvtstagn.o $(OBJ)cmvschdist.o
+          $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o $(OBJ)cmvtstagn.o $(OBJ)cmvschdist.o \
+          $(OBJ)cmvhshsns.o
  
 LIBROBJ = \
@@ -272,2 +273,10 @@
 
 ##############################################################################
+cmvhshsns: $(EXE)cmvhshsns
+	echo $@ " done"
+$(EXE)cmvhshsns: $(OBJ)cmvhshsns.o $(LIBR)
+	$(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvhshsns.o $(MYLIB) 
+$(OBJ)cmvhshsns.o: cmvhshsns.cc
+	$(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvhshsns.cc
+
+##############################################################################
Index: trunk/Cosmo/SimLSS/cmvhshsns.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvhshsns.cc	(revision 3632)
+++ trunk/Cosmo/SimLSS/cmvhshsns.cc	(revision 3632)
@@ -0,0 +1,361 @@
+// lobes pour HSHS N-S
+// > cmvhshsns -n -d 0.105 -g 4,0.105,0. -i 2,0.41,0. -t 50,180,0 -f 1420 1410 1430
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "ntuple.h"
+#include "constcosmo.h"
+#include "geneutils.h"
+
+//----------------------------------------------------------------
+void usage(void);
+void usage(void)
+{
+cout<<"cmvhshsns [...] val1 val2 ..."<<endl
+  <<"  -f : val1... sont des frequences en MHz"<<endl
+  <<"       (defaut: longueurs d'onde en m"<<endl
+  <<"  -n : hauteur des lobes normalises a 1"<<endl
+  <<"  -d L : longueur totale dipole"<<endl
+  <<"         L==0. alors approximation du  dipole"<<endl
+  <<"  -g N_g,D_g,Theta_g : regroupement de dipoles"<<endl
+  <<"                       N_g<=1 pas de regroupement"<<endl
+  <<"  -i N_i,D_i,Theta_i : interferences entre regroupements de dipoles"<<endl
+  <<"                       N_i<=1 pas d'interference"<<endl
+  <<"  -t Nang,Tmax,Tcent : nombre de pts entre 2 zeros consecutifs (def=25)"<<endl
+  <<"                       angle maxi (deg, def=180), angle central (deg, def=0)"<<endl
+  <<"  -p : interprete Theta_{g,i} en picosecondes"<<endl
+  <<"..distances L,D: >0 en m , <0 en unites de longeur d'onde"<<endl
+  <<"..angles Theta: en deg"<<endl
+  <<"                en secondes si option \"-p\""<<endl
+  <<endl;
+}
+
+double thetafromdt(double &theta,double dt,double lambda,double dconsec);
+double dtfromtheta(double theta,double lambda,double dconsec);
+
+//----------------------------------------------------------------
+int main(int narg,char *arg[])
+{
+ const double torad = M_PI/180.;
+
+ // --- longueur d'onde em m
+ vector<double> Lambda, Nu;
+ bool argfreq = false;
+ bool thetaps = false;
+ bool normone = false;
+
+ // --- dipole de longeur totale L (2 brins de L/2)
+ double L_d = -0.5;   // >0 en m, <0 en unite de lambda
+
+ // --- groupes: regroupement des dipoles
+ // N est le nombre de dipoles regroupes
+ // D est la distance entre deux dipoles consecutifs
+ // Theta est le dephasage entre 2 dipoles consecutifs (dephasage electronique)
+ //  i.e. c'est l'angle d'arrivee des rayons qui sont combines en phase
+ int N_g = 1;  // nombre de dipoles regroupes
+ double D_g = -0.5;    // >0 en m, <0 en unite de lambda
+ double Theta_g = 0.;  // en deg
+
+ // --- interference des groupes
+ // N est le nombre de groupes sur la ligne focale
+ // D est la distance entre deux groupes consecutifs
+ // Theta est le dephasage entre 2 groupes consecutifs (dephasage electronique)
+ int N_i = 100;  // nombre de groupes
+ double D_i = (N_g>1) ? N_g*D_g:  -0.5;    // >0 en m, <0 en en unite de lambda
+ double Theta_i = 0.;  // en deg
+
+ // --- Tmax = angle maximum de scan a partir du zenith (degres)
+ double Tcent = 0., Tmax = 180.;
+
+ // --- Nang nombre de points entre deux zeros de la figure d'interfrence
+ int Nang = -1;
+
+ // Decodage des arguments
+ char c;
+ while((c = getopt(narg,arg,"hpnfd:g:i:t:")) != -1) {
+  switch (c) {
+  case 'f' :
+    argfreq = true;
+    break;
+  case 'p' :
+    thetaps = true;
+    break;
+  case 'n' :
+    normone = true;
+    break;
+  case 'd' :
+    sscanf(optarg,"%lf",&L_d);
+    break;
+  case 'i' :
+    sscanf(optarg,"%d,%lf,%lf",&N_i,&D_i,&Theta_i);
+    break;
+  case 'g' :
+    sscanf(optarg,"%d,%lf,%lf",&N_g,&D_g,&Theta_g);
+    break;
+  case 't' :
+    sscanf(optarg,"%d,%lf,%lf",&Nang,&Tmax,&Tcent);
+    break;
+  case 'h' :
+  default :
+    usage();
+    return -1;
+    break;
+  }
+ }
+
+ if(optind>=narg) {usage(); return -2;}
+ for(int i=optind;i<narg;i++) {
+   double v = atof(arg[i]);
+   if(v<=0.) continue;
+   if(argfreq) {
+     v *= 1.e6;
+     Nu.push_back(v);
+     Lambda.push_back(SpeedOfLight_Cst*1.e3/v);
+   } else {
+     Lambda.push_back(v);
+     Nu.push_back(SpeedOfLight_Cst*1.e3/v);
+   }
+ }
+ cout<<"Nombre de longueurs d'onde a traiter "<<Lambda.size()<<" (unites: m et Hz)"<<endl;
+ if(Lambda.size()==0) return -3;
+ for(unsigned short i=0;i<Lambda.size();i++)printf(" %.3f m ,  %.3f MHz\n",Lambda[i],Nu[i]/1.e6);
+
+ cout<<"Dipole : longueur totale L="<<L_d<<endl;
+ if(L_d==0.) return -4;
+
+ if(N_g<=0) N_g = 1;
+ cout<<"Regroupements N="<<N_g<<" D="<<D_g<<" Theta="<<Theta_g<<endl;
+ if(N_g>1 && D_g==0.) return -4;
+
+ if(N_i<=0) N_i = 1;
+ cout<<"Interferences N="<<N_i<<" D="<<D_i<<" Theta="<<Theta_i<<endl;
+ if(N_i>1 && D_i==0.) return -4;
+
+ if(Nang<=0) Nang = 25;
+ while(Tcent<-180.) Tcent += 360.;
+ while(Tcent>180.)  Tcent -= 360.;
+ Tmax = fabs(Tmax); if(Tmax>180.) Tmax=180.;
+ cout<<"Display: Nang="<<Nang<<" Tmax="<<Tmax<<" Tcent="<<Tcent<<endl;
+
+ double norme_g=1., norme_i=1.;
+ if(normone) {norme_g=N_g*N_g; norme_i=N_i*N_i;}
+ cout<<"Normalisation: norme_g="<<norme_g<<"  norme_i="<<norme_i<<endl;
+
+ char str[32];
+ POutPersist pos("cmvhshsns.ppf");
+
+ //------- Boucle sur les longeurs d'ondes
+ for(unsigned short il=0;il<Lambda.size();il++) {
+
+   //... mise en forme des parametres pour la longueur d'onde
+   double lambda = Lambda[il], nu = Nu[il];
+   cout<<"\n\n>>> Lambda = "<<lambda<<" m ,  nu = "<<nu/1.e6<<" MHz"<<endl;
+   double ld = (L_d<0.) ? -L_d*lambda : L_d;
+   cout<<"dipole: ld="<<ld<<" m"<<endl;
+   double dg = (D_g<0.) ? -D_g*lambda : D_g;
+   cout<<"groupe: ("<<N_g<<"), dg="<<dg<<" m -> "<<dg*N_g<<" m"<<endl;
+   double di = (D_i<0.) ? -D_i*lambda : D_i;
+   cout<<"interf: ("<<N_i<<"), di="<<di<<" m -> "<<di*N_i<<" m"<<endl;
+   double thg, thi, tg, ti;
+   if(thetaps) {
+     tg = Theta_g;
+     double rcg = thetafromdt(thg,tg,lambda,dg);
+     ti = Theta_i;
+     double rci = thetafromdt(thi,ti,lambda,di);
+     cout<<"dephasage electronique: groupe "<<tg<<" sec -> "<<thg/torad<<" deg pour rc="<<rcg<<endl
+         <<"                        interf "<<ti<<" sec -> "<<thi/torad<<" deg pour rc="<<rci<<endl;
+     if(rcg>1 || rci>1) {cout<<"!!!! Lambda NON-TRAITE"<<endl; continue;}
+   } else {
+     thg = Theta_g*torad; tg  = dtfromtheta(thg,lambda,dg);
+     thi = Theta_i*torad; ti  = dtfromtheta(thi,lambda,di);
+     cout<<"dephasage electronique: groupe "<<thg/torad<<" deg -> "<<tg<<" sec"<<endl
+         <<"                        interf "<<thi/torad<<" deg -> "<<ti<<" sec"<<endl;
+   }
+   double sthg = sin(thg);
+   double sthi = sin(thi);
+
+   //... distance approx entre 2 zeros
+   double dzero = M_PI/2.;
+   if(N_g>1) {
+     double z = lambda/(N_g*dg);
+     cout<<"groupe: distance entre 2 zeros: d(sin(t))="<<z<<" (approx "<<z/torad<<" deg)"<<endl;
+     if(z<dzero) dzero = z;
+   }
+   if(N_i>1) {
+     double z = lambda/(N_i*di);
+     cout<<"interf: distance entre 2 zeros: d(sin(t))="<<z<<" (approx "<<z/torad<<" deg)"<<endl;
+     if(z<dzero) dzero = z;
+   }
+   cout<<"distance approx entre 2 zeros "<<dzero<<"  rad = "<<dzero/torad<<" deg"<<endl;
+
+   //... remplissage des angles
+   {
+   const int nnt=5; float xnt[nnt];
+   const char *namev[nnt] = {"t","ant","intfg","intfi","intf"};
+   NTuple nt(nnt,namev);
+
+   long npt = long(Tmax*torad/dzero*Nang +0.5);
+   double dt = Tmax*torad/npt;
+   cout<<"nombre de points dans la boucle "<<2*npt+1<<" , dt="<<dt/torad<<endl;
+   for(int i=-npt;i<=npt;i++) {
+     double t = Tcent*torad + i*dt;
+     double st = sin(t);
+     double ta = acos(st); // angle par rapport au fil de l'antenne
+     double deltag = M_PI*dg/lambda*(st-sthg);
+     double deltai = M_PI*di/lambda*(st-sthi);
+     double ant = (ld==0.) ? AntDipole(ld/lambda,ta): AntCentFed(ld/lambda,ta);
+     double intfg = (N_g==1) ? 1.: SinNXsX_Sqr(deltag,N_g)/norme_g;
+     double intfi = (N_i==1) ? 1.: SinNXsX_Sqr(deltai,N_i)/norme_i;
+     double intf = ant*intfg*intfi;
+
+     xnt[0] = t/torad;
+     xnt[1] = ant;
+     xnt[2] = intfg;
+     xnt[3] = intfi;
+     xnt[4] = intf;
+     nt.Fill(xnt);
+   }
+
+   //...ecriture ppf
+   sprintf(str,"nt_%d",il);
+   cout<<"writing "<<str<<" into ppf file"<<endl;
+   pos << PPFNameTag(str) << nt;
+   DVList dvl;
+   dvl("Lambda") = lambda; dvl("Nu") = nu;
+   dvl("Ld") = ld;
+   dvl("Ng") = N_g; dvl("Dg") = dg; dvl("Thg") = thg; dvl("Tg") = tg;
+   dvl("Ni") = N_i; dvl("Di") = di; dvl("Thi") = thi; dvl("Ti") = ti;
+   dvl("Tmax") = Tmax; dvl("Tcent") = Tcent;
+   sprintf(str,"dvl_%d",il);
+   pos << PPFNameTag(str) << dvl;
+   }
+
+   //... remplissage des zeros et des maximas principaux
+   {
+   const int nnt=2; double xnt[nnt];
+   const char *namev[nnt] = {"t","ztyp"};
+   if(N_g>1) {
+     NTuple nt(nnt,namev);
+     double no = (normone) ? 1.: norme_g;
+     for(int is=-1;is<=1;is+=2) {
+       int k0 = (is==1) ? 0: -1;
+       for(int k=k0;;k+=is) {
+         xnt[0] =  k*lambda/(N_g*dg)+sthg;
+         xnt[1] = (k%N_g==0) ? no: 0.;
+         if(fabs(xnt[0])>1.) break;
+         xnt[0]=asin(xnt[0])/torad; nt.Fill(xnt);
+       }
+     }
+     sprintf(str,"ntzg_%d",il);
+     pos << PPFNameTag(str) << nt;
+   }
+   if(N_i>1) {
+     NTuple nt(nnt,namev);
+     double no = (normone) ? 1.: norme_i;
+     for(int is=-1;is<=1;is+=2) {
+       int k0 = (is==1) ? 0: -1;
+       for(int k=k0;;k+=is) {
+         xnt[0] =  k*lambda/(N_i*di)+sthi;
+         xnt[1] = (k%N_i==0) ? no: 0.;
+         if(fabs(xnt[0])>1.) break;
+         xnt[0]=asin(xnt[0])/torad; nt.Fill(xnt);
+       }
+     }
+     sprintf(str,"ntzi_%d",il);
+     pos << PPFNameTag(str) << nt;
+   }
+   }
+
+ }
+
+ return 0;
+}
+
+//-----------------
+double thetafromdt(double &theta,double dt,double lambda,double dconsec)
+// Input:
+//   dt : decalage en seconde
+//   lambda : longueur d'onde en m
+//   dconsec : distance entre 2 dipoles consecutifs en m
+// Output et Return code "rc":
+//   on cherche theta tel que: "C*dt = dconsec*sin(theta)"
+//   ...si on peut trouver theta on renvoie
+//        theta (en rad) angle equivalent
+//        rc = 1
+//   ...si on ne peut pas trouver theta on renvoie
+//        theta (en rad) angle equivalent
+//        rc > 1
+//      pour la premiere solution de "C*dt = rc*dconsec*sin(theta)"
+//   ...mauvais parametre on renvoie
+//        rc=0
+{
+  if(dconsec<=0. || lambda<=0.) return 0.;
+
+  // conversion de dt en longueur
+  dt *= SpeedOfLight_Cst*1.e3;
+
+  // translation dans une longueur d'onde
+  dt = (dt/lambda - trunc(dt/lambda))*lambda;
+
+  // calcul de sin(theta)
+  double st = dt/dconsec;
+
+  // calcul du nombre de dconsec pour avoir la premiere solution
+  double rc = ceil(fabs(st));
+  if(rc==0) rc=1;
+  
+  // equivalent en angle possible
+  theta = asin(st/rc);
+  return rc;
+ }
+
+double dtfromtheta(double theta,double lambda,double dconsec)
+// Input:
+//   theta : angle en rad
+//   lambda : longueur d'onde en m
+//   dconsec : distance entre 2 dipoles consecutifs en m
+// Return: dt = dconsec/C * sin(theta) remis dans la periode
+{
+  double dt = dconsec*sin(theta);
+
+  // translation dans une longueur d'onde
+  dt = (dt/lambda - trunc(dt/lambda))*lambda;
+
+  return dt/(SpeedOfLight_Cst*1.e3);
+}
+
+/******************************
+openppf cmvhshsns.ppf
+
+set l 0
+
+set t t
+set t sin(t*M_PI/180.)
+
+set cut 1
+set cut -90<t&&t<90
+
+n/plot $nt.$t%_nl
+
+n/plot nt_$l.ant%$t $cut ! "nsta cpts green"
+n/plot nt_$l.intfg%$t $cut ! "nsta cpts same red"
+n/plot nt_$l.intfi%$t $cut ! "nsta cpts same blue"
+n/plot nt_$l.intf%$t $cut ! "nsta cpts same black"
+
+n/plot ntzg_$l.ztyp%$t $cut ! "nsta same marker=circle,9 red"
+n/plot ntzi_$l.ztyp%$t $cut ! "nsta same marker=star,9 blue"
+
+
+# compare frequences
+n/plot nt_0.intf%$t $cut ! "nsta cpts black"
+n/plot nt_1.intf%$t $cut ! "nsta cpts same blue"
+n/plot nt_2.intf%$t $cut ! "nsta cpts same red"
+n/plot nt_3.intf%$t $cut ! "nsta cpts same orange"
+
+ */
Index: trunk/Cosmo/SimLSS/geneutils.cc
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.cc	(revision 3622)
+++ trunk/Cosmo/SimLSS/geneutils.cc	(revision 3632)
@@ -361,4 +361,78 @@
 
 //-------------------------------------------------------------------
+/*=========================================================
+Antenne de longueur totale d (lue au milieu cad 2 brins de d/2)
+see Jackson p401  (antenne selon Oz, centre z=0 en O)
+k = w/C = 2 Pi Nu / C = 2 Pi / Lambda
+- densite de courant: J = I sin(kd/2 - k|z|)*delta(x)*delta(y)
+                        (selon Oz) et pour |z|<d/2
+                        I current peak value if kd>=Pi
+- current at gap (center): I0 = I sin(kd/2)
+#
+Z=sqrt(Mu/Epsilon) is the impedance of the vacum (376.7 Ohms)
+  cad pour le vecteur de pointing: |S| = |E0|^2/Z
+- Time average radiated power by solid angle unit:
+ t is theta the polar angle, t=0 along antenna, Pi/2 perpendicular
+ dP/dOmega = Z*I^2/(8Pi^2) * (cos(kd/2*cos(t)) - cos(kd/2))^2 / sin(t)^2
+           if we define  L = d/lambda ,  kd/2= Pi * L
+ dP/dOmega = Z*I^2/(8Pi^2) * (cos(Pi*L*cos(t)) - cos(Pi*L))^2 / sin(t)^2
+#
+- pour t -> 0
+          dP/dOmega = Z*I^2/(8Pi^2) * ( 1/4 * (Pi*L*sin(Pi*L)*t)^2 )
+- pour t -> Pi
+          dP/dOmega = Z*I^2/(8Pi^2) * ( 1/4 * (Pi*L*sin(Pi*L)*(t-Pi))^2 )
+#
+- Polarisation in plane containing antenna and the direction of propagation
+#
+- pour kd << 1  (L->0)  -> approximation du dipole:
+          dP/dOmega = Z*I^2/(8Pi^2) * ( (Pi*L)^4/4 * sin(t)^2 )
+- pour kd=Pi  (L=1/2)  -> half-wave length:
+          dP/dOmega = Z*I^2/(8Pi^2) * cos(Pi/2*cos(t))^2/sin(t)^2
+- pour kd=2*Pi  (L=1)  -> full-wave length:
+          dP/dOmega = Z*I^2/(8Pi^2) * 4*cos(Pi/2*cos(t))^4/sin(t)^2
+#
+- Le programme retourne Val tel que: dP/dOmega = Z*I^2/(8Pi^2) * Val
+    L<0 (arg[0]<0) alors Val est normalise tel que le maximum = 1
+#
+La carte antenne fournie par Peterson:
+         d=8cm, distance entre 2 centres consecutifs 10cm
+#
+# ---- AntCentFed et AntDipole
+# Input:
+#   t = angle entre la direction de la radiation et le fil de l'antenne (radian)
+#   L = d/Lambda  (d=longueur totale de l'antenne)
+# Output:
+#   (8Pi^2) dP/dOmega / (Z*I^2)
+========================================================= */
+
+double AntCentFed(double L, double trad)
+{
+ 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;
+ }
+
+ return v;
+}
+
+double AntDipole(double L, double trad)
+{
+ double pil = M_PI*L;
+ double st = sin(trad);
+ // L'approximation du dipole
+ double vd = pil*pil*st/2.;  vd *= vd;
+ return vd;
+}
+
+//-------------------------------------------------------------------
 
 static unsigned short IntegrateFunc_GlOrder = 0;
Index: trunk/Cosmo/SimLSS/geneutils.h
===================================================================
--- trunk/Cosmo/SimLSS/geneutils.h	(revision 3622)
+++ trunk/Cosmo/SimLSS/geneutils.h	(revision 3632)
@@ -94,4 +94,7 @@
 double SinNXsX_Sqr(double x,unsigned long N,bool app=false);
 
+double AntCentFed(double L, double trad);
+double AntDipole(double L, double trad);
+
 //----------------------------------------------------
 double IntegrateFunc(GenericFunc& func,double xmin,double xmax
