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

// ----
// ---- Pattern de diffraction Est-Ouest pour HSHS
// ----
// ...Figure de DIFFRACTION
//  Le detecteur de chacune des antennes reste toujours dans le plan focal
//    sur l'axe optique E-W.
//  La puissance mesuree par l'antenne est donc celle mesuree dans le plan focal
//    sur l'axe optique.
//  Pour un rayon arrivant avec l'incidence theta, la valeur mesuree sur l'axe
//    optique est celle de la figure de diffraction d'un rayon arrivant
//    paralellement a l'axe optique (theta=0) mais mesuree en -theta et MULTIPLIEE
//    par cos(theta) puisque la surface collectrice effective est plus faible.
//  Le programme renvoie la figure de diffraction "diff" dont le maximum est en "t=0"
//    NON-MULTIPLIEE par cos(theta)
// ...Figure d'INTERFERENCE
//  Si tous les detecteurs de chaque antenne sont lus et additionnes en phase,
//    la figure d'interference est donnee par les differences de marche
//    dues a la variation de l'angle d'incidence des rayons.
//    Cette figure d'interference est celle donnee dans le programme ci-apres par "intf".
//  Mais l'experience mesure les amplitudes elementaires sur chaque antenne et les additionne
//    en ajoutant ("electroniquement") une phase "phi" a chaque antenne.
//    La figure d'interference est donc decalee telle que son maximum d'ordre
//    zero (rayons arrivant en phase) soit centre en "phi".
//    Dans le programme le dephasage "phi" est exprime par l'angle d'incidence theta0
//    des rayons qui vont arriver en phase.
//    La difference de marche appliquee a une antenne
//    a la position X sera delta=X*sin(theta0)
// ...Normalisation
//  Les courbes "intf" et "diff" sont chacunes telles que leurs valeurs maximales
//    sur tout le domaine en angle sont egales a 1
// ...Resultat
//  Le resultat final pour obtenir la contribution a la puissance mesuree
//    par l'antenne en fonction de la direction du rayon incident est donc:
//           ----  Puissance = intf * diff * cos(theta)  ----
//

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

// Lambda = longueur d'onde (metres)
double Lambda = 0.21;

// Tmax = angle maximum de scan a partir du zenith (degres)
//        (si <0   alors   -Tmax* diffraction_first_zero)
// Nang nombre de points entre deux maxi d'interference
//double Tmax = 90.;
double Tmax = -5.;
double Nang = 100;

// Theta0 = angle d'incidence (degre) des rayons qui arrivent en phase
//       (traduit le dephasage de phase electronique applique aux antennes)
double Theta0 = 0.;

// Dcyl = largeur du cylindre en metres
// Ncyl = nombre cylindres
// Xcyl = distance entre les cylindres en metres
double Dcyl = 10.;
const int Ncyl = 10;
double Xcyl[Ncyl];
Xcyl[0] = 0.;
 for(int i=1;i<Ncyl;i++) Xcyl[i] = Xcyl[i-1] + Dcyl;
 //for(int i=1;i<Ncyl;i++) Xcyl[i] = Xcyl[i-1] + Dcyl + (i-1)*Dcyl/2.;
 //for(int i=1;i<Ncyl;i++) Xcyl[i] = Xcyl[i-1] + Dcyl + (i-1)*Dcyl/Ncyl;
 // for(int i=1;i<Ncyl;i++) Xcyl[i] = Xcyl[i-1] -Dcyl/2. + Dcyl + (i-1)*Dcyl/Ncyl;

//---------------------------------------------
cout<<"Lambda="<<Lambda<<" m"<<endl;
cout<<"Dcyl="<<Dcyl<<" m"<<endl;
cout<<"Ncyl="<<Ncyl<<endl;
for(int i=0;i<Ncyl;i++) {
  cout<<"   Xcyl["<<i<<"] = "<<Xcyl[i]<<" m";
  if(i>0) cout<<"  d = "<<Xcyl[i]-Xcyl[i-1];
  cout<<endl;
}

//---------------------------------------------
if(Dcyl<=0. || Ncyl<=0) return -1;
double dxmin = 1.e99, dxmax = -1.e99;
for(int i=0;i<Ncyl-1;i++) for(int j=i+1;j<Ncyl;j++) {
  double d = fabs(Xcyl[j]-Xcyl[i]);
  if(d<dxmin) dxmin=d;
  if(d>dxmax) dxmax=d;
}
cout<<"Minimum distance between cylinders = "<<dxmin<<endl;
if(dxmin<Dcyl) cout<<"WARNIG: cylinder width > minimum cylinder distance"<<endl;
cout<<"Maximum distance between cylinders = "<<dxmax<<endl;
if(dxmin<=0. || dxmax<=0.) return -2;

double samin = Lambda/dxmin;
double samax = Lambda/dxmax;
cout<<"Interference secondary maxi for max distance "<<samax<<" in sin() units -> "
    <<asin(samax)/torad*60.<<"\'"<<endl;
cout<<"Interference principal maxi for min distance "<<samin<<" in sin() units -> "
    <<asin(samin)/torad*60.<<"\'"<<endl;

double s0diff = Lambda/Dcyl;
cout<<"Diffraction first zero "<<s0diff<<" in sin() units -> "
    <<asin(s0diff)/torad<<" deg"<<endl;

double s_theta0 = sin(Theta0*torad);
cout<<"Dephasage electronique exprime en angle "<<Theta0<<" deg -> sin="<<s_theta0<<endl;

if(Tmax<=0.) Tmax = -Tmax *asin(s0diff)/torad;
if(Tmax>90.) Tmax = 90.;
double sainc = samax/Nang;
cout<<"Scan 0 to "<<Tmax<<" deg from zenith"<<endl;
cout<<"Using "<<Nang<<" points between interference maxi"<<endl;
cout<<"Sinus increment="<<sainc<<" -> "<<asin(sainc)/torad*60.<<"\'"<<endl;
if(sainc<=0.) return -1;

Tmax *= torad;
if(Tmax>M_PI_2) Tmax=M_PI_2;
int nest = int(2.*sin(Tmax)/sainc)+1;
cout<<"Estimated number of entries "<<nest<<endl;

//---------------------------------------------
const char *namev[3] = {"t","intf","diff"};
float xnt[3];
NTuple nt(3,namev);

int lpmod = (nest>10)? nest/10: 1;
int nent = 0;
// On parcourt en sin(angle)
for(double sa=-sin(Tmax);sa<sin(Tmax) && sa<=1.;sa+=sainc) {
  double a = asin(sa);
  // Interference (add dephasage electronique)
  double sumc = 0., sums = 0.;
  for(int n=0;n<Ncyl;n++) {
    double d = 2.*M_PI*(sa-s_theta0)*(Xcyl[n]-Xcyl[0])/Lambda;
    sumc += cos(d);
    sums += sin(d);
  }
  double intf = (sumc*sumc+sums*sums)/(Ncyl*Ncyl);
  // Diffraction (le detecteur est toujour sur l'axe optique)
  double diffr = M_PI*sa*Dcyl/Lambda;
  diffr = SinXsX_Sqr(diffr);
  //
  if(nent%lpmod==0)
    cout<<nent<<" sa="<<sa<<" a="<<a/torad<<" intf="<<intf<<" diff="<<diffr<<endl;
  // Fill NTuple
  xnt[0] = a/torad;
  xnt[1] = intf;
  xnt[2] = diffr;
  nt.Fill(xnt);
  nent++;
}
cout<<"Nent="<<nent<<endl;

POutPersist pos("cmvhshsew.ppf");
pos << PPFNameTag("nt") << nt;

return 0;
}

/*
del nt
openppf cmvhshsew.ppf

set norm cos(t*M_PI/180.)

set t t
set t sin(t*M_PI/180.)

n/plot nt.intf%$t ! ! "nsta connectpoints"
n/plot nt.diff%$t ! ! "nsta connectpoints same red"
n/plot nt.$norm%$t ! ! "nsta connectpoints same orange"
n/plot nt.diff*intf*$norm%$t ! ! "nsta connectpoints same blue"

# Plot en dB
n/plot nt.10.*log10(diff*intf*$norm)%$t diff*intf>0. ! "nsta connectpoints"
addline -90 0 90 0 "red"

# Dephasage electronique d'angle t (deg)
set t0 0.6
n/plot nt.intf%sin(t*M_PI/180.)-sin($t0*M_PI/180.) ! ! "nsta connectpoints"
n/plot nt.diff%sin(t*M_PI/180.) ! ! "nsta connectpoints same red"

 */