// 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"

 */