// 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
// > cmvhshsns -n -d 0.105 -g 2,0.105,0. -i 2,0.41,0. -t 50,180,0 -f 1420 1410 1430
// > cmvhshsns -n -D -1 -d 10 -g 1 -i 10,10,0. -t 50,180,0 -f 1420 1410 1430
// - Pittsburgh Sun 24 Nov
//   > cmvhshsns -n -d 0.147 -D 0 -g 1 -i 8,0.147,61.06 -t 360,90,0 -f 1420 
#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
  <<"         si \"-D -1\" largeur du cylindre (etude E-W)"<<endl
  <<"  -D : +1 lobe du dipole approximation de Hertz"<<endl
  <<"       -1 le lobe du dipole est remplace par le sinc() du cylindre (etude E-W)"<<endl
  <<"        0 lobe du dipole a 2 brins de longueur \"L\""<<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 de temps"<<endl
  <<"..distances L,D: >0 en m , <0 en unites de longeur d'onde"<<endl
  <<"..angles Theta: en deg"<<endl
  <<"                en picosecondes de temps si option \"-p\""<<endl
  <<endl;
}

double thetafromdt(double &theta,double dt,double lambda,double dconsec);
double dtfromtheta(double theta,double lambda,double dconsec);
double princpiclarge(double lambda,int N, double d,double sth=0.);

//----------------------------------------------------------------
int main(int narg,char *arg[])
{
 const double torad = M_PI/180.;

 // --- longueur d'onde en 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
 int dipolType = 0; // 1=dipole de Hertz, -1=sinc(), sinon dipole a 2 brins

 // --- 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:d:g:i:t:")) != -1) {
  switch (c) {
  case 'f' :
    argfreq = true;
    break;
  case 'p' :
    thetaps = true;
    break;
  case 'n' :
    normone = true;
    break;
  case 'D' :
    dipolType = atoi(optarg);
    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 : type="<<dipolType<<", longueur (ou largeur) 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  (type="<<dipolType<<")"<<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);

   //... largeur du pic principale (distance entre les 2 zeros de part et d'autre)
   if(N_g>1) {
     double t = princpiclarge(lambda,N_g,dg,sthg);
     if(t>0.) cout<<"groupe: largeur du lobe principal: "<<t<<" rad   "<<t/torad<<" deg"<<endl;
   }
   if(N_i>1) {
     double t = princpiclarge(lambda,N_i,di,sthi);
     if(t>0.) cout<<"interf: largeur du lobe principal: "<<t<<" rad   "<<t/torad<<" deg"<<endl;
   }

   //... distance approx entre 2 zeros (hors pic principale)
   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;
   double vhmax = -1.e20, thmax=0.;
   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;
     if(dipolType<0)        ant = LobeSinc(ld/lambda,t);
       else if(dipolType>0) ant = AntDipole(ld/lambda,ta);
         else               ant = 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;
     if(fabs(t)<M_PI/2. && intf>vhmax) {vhmax=intf; thmax=t;}

     xnt[0] = t/torad;
     xnt[1] = ant;
     xnt[2] = intfg;
     xnt[3] = intfi;
     xnt[4] = intf;
     nt.Fill(xnt);
   }
   cout<<"Found maxi ("<<vhmax<<") for t = "<<thmax<<" rad = "<<thmax/torad<<" deg"<<endl;

   //...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("dipolType") = dipolType;
   dvl("Ng") = N_g; dvl("Dg") = dg; dvl("Thg") = thg/torad; dvl("Tg") = tg;
   dvl("Ni") = N_i; dvl("Di") = di; dvl("Thi") = thi/torad; dvl("Ti") = ti;
   dvl("Tmax") = Tmax; dvl("Tcent") = Tcent;
   dvl("vhmax") = vhmax; dvl("thmax") = thmax/torad;
   sprintf(str,"dvl_%d",il);
   pos << PPFNameTag(str) << dvl;
   }

   //... remplissage des zeros et des maximas principaux
   for(int i=0;i<2;i++) {  // i==0 groupes , i==1 interferences
     int N = N_g;
     double d = dg, sth = sthg, norme = norme_g;
     sprintf(str,"ntzg_%d",il);
     if(i==1) {
       N = N_i;
       d = di; sth = sthi; norme = norme_i;
       sprintf(str,"ntzi_%d",il);
     }
     if(N<=1) continue;
     const int nnt=2; double xnt[nnt];
     const char *namev[nnt] = {"t","ztyp"};
     NTuple nt(nnt,namev);
     double no = (normone) ? 1.: norme;
     for(int is=-1;is<=1;is+=2) {  // ki>=0 ou k<0
       int k0 = (is==1) ? 0: -1;
       for(int k=k0;;k+=is) {
         xnt[0] =  k*lambda/(N*d)+sth;
         xnt[1] = (k%N==0) ? no: 0.;  // zero ou maxi?
         if(fabs(xnt[0])>1.) break;
         xnt[0]=asin(xnt[0])/torad; nt.Fill(xnt);
       }
     }
     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);
}

double princpiclarge(double lambda,int N, double d,double sth)
// Input:
//   lambda : longueur d'onde
//   N : nombre de dipole dans le regroupement
//   d : distance entre deux dipoles
//   sth : sin(theta) equivalent pour dephasage "electronique"
// Return: si >0 distance angulaire en radians
//         si -1 zero haut (droite) non-existant
//         si -2 zero bas (gauche) non-existant
{
 double st1 = lambda/(N*d)+sth;
 if(fabs(st1)>1.) return -1.;
 double st2 = -lambda/(N*d)+sth;
 if(fabs(st2)>1.) return -2.;
 return fabs(asin(st1)-asin(st2));
}

/******************************
delobjs *
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

c++exec cout<<dvl_${l}<<endl;

n/plot nt_$l.$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"

 */