#include<iostream.h>
#include<stdlib.h>
#include<fstream.h>
#include<math.h>
#include "Rayon.h"

const double PI=3.14159265359;

Rayon::Rayon(double d, double s, int n, int ech){
dist=d;   //distance de mesure (m)
seuil=s;   //seuil pour le calcul du lcr (dB)
nb_ondes=n;   //nombre d'ondes incidentes considérées
echelle=ech;   //echelle pour le calcul des distributions
E0=10E-6;   //amplitude du champ (V)
beta=0;   //angle entre l'onde incidente et le sol (radian)
gamma=0;   //angle entre le vecteur vitesse du mobile et l'axe des x (radian)
vit=10;   // vitesse du mobile (m/s)
Fe=1000;   //fréquence d'échantillonnage (Hz)
freq=900E6;   //fréquence de la communication (Hz)
c=3E8;   //vitesse de la lumière (m/s)
lambda=c/freq;   //longueur d'onde de la communication (m)
nb_echant=(int)(dist/vit*Fe);   //nombre d'échantillons
lcr=0;   //taux de franchissement d'un seuil (Hz)
tau=0;   //durée moyenne d'évanouissement (s)
alpha=new double[nb_ondes];   //angle azimutal d'arrivée (radian)
phi=new double[nb_ondes];   //phase de l'onde (radian)
I=new double[nb_echant];   //partie en phase (V)
Q=new double[nb_echant];   //partie en quadrature (V)
r=new double[nb_echant];   //enveloppe du signal (V)
loi_I=new double[echelle];   //distribution de I(t)
loi_Q=new double[echelle];   //distribution de Q(t)
loi_r=new double[echelle];   //distribution de r(t)
r_dB=new double[nb_echant];	  //enveloppe du signal (dB)
}

Rayon::~Rayon(){
delete alpha;
delete phi;
delete I;
delete Q;
delete r;
delete loi_I;
delete loi_Q;
delete loi_r;
cout<<"Fin de la simulation\n";
}

void Rayon::calcul_IQ(){
int tps=(int)(nb_echant/nb_ondes);
for(int i=0;i<nb_echant;i++){
I[i]=0;
Q[i]=0;
for(int j=0;j<nb_ondes;j++){
if(i%tps==0){
alpha[j]=((double)rand()/RAND_MAX)*2*PI;
phi[j]=((double)rand()/RAND_MAX)*2*PI;
}
I[i]+=E0*cos(2*PI*vit/lambda*cos(gamma-alpha[j])*cos(beta)*i/Fe+phi[j]);
Q[i]+=E0*sin(2*PI*vit/lambda*cos(gamma-alpha[j])*cos(beta)*i/Fe+phi[j]);
}
}
ecrire("I.xls", I, nb_echant);
ecrire("Q.xls", Q, nb_echant);
}

void Rayon::calcul_r(){
for(int i=0;i<nb_echant;i++){
r[i]=sqrt(I[i]*I[i]+Q[i]*Q[i]);
r_dB[i]=10*log(r[i]);
}
ecrire("r.xls", r, nb_echant);
ecrire("r_dB.xls", r_dB, nb_echant);
}

void Rayon::calcul_loi_IQ(){
double minI, maxI, minQ, maxQ, increment_I, increment_Q;
minI=I[0];
maxI=I[0];
minQ=Q[0];
maxQ=Q[0];
for(int i=1;i<nb_echant;i++){
if(I[i]<minI) minI=I[i];
if(I[i]>maxI) maxI=I[i];
if(Q[i]<minQ) minQ=Q[i];
if(Q[i]>maxQ) maxQ=Q[i];
}
increment_I=(maxI-minI)/echelle;
increment_Q=(maxQ-minQ)/echelle;
for(int k=0;k<echelle;k++){
loi_I[k]=0;
loi_Q[k]=0;
for(int j=0;j<nb_echant;j++){
if(minI+k*increment_I<I[j] && I[j]<minI+(k+1)*increment_I) loi_I[k]++;
if(minQ+k*increment_Q<Q[j] && Q[j]<minQ+(k+1)*increment_Q) loi_Q[k]++;
}
}
ecrire("loi_I.xls", loi_I, echelle);
ecrire("loi_Q.xls", loi_Q, echelle);
}

void Rayon::calcul_loi_r(){
double minr, maxr, increment_r;
minr=r[0];
maxr=r[0];
for(int i=1;i<nb_echant;i++){
if(r[i]<minr) minr=r[i];
if(r[i]>maxr) maxr=r[i];
}
increment_r=(maxr-minr)/echelle;
for(int k=0;k<echelle;k++){
loi_r[k]=0;
for(int j=0;j<nb_echant;j++){
if(minr+k*increment_r<r[j] && r[j]<minr+(k+1)*increment_r) loi_r[k]++;
}
}
ecrire("loi_r.xls", loi_r, echelle);
}

double Rayon::moyenne(double *t, int taille, int puissance){
double m=0;
for(int i=0;i<taille;i++) m+=pow(*(t+i), puissance);
m/=taille;
return m;
}

double Rayon::variance(double *t, int taille){
double v=0;
v=moyenne(t, taille, 2)-moyenne(t, taille)*moyenne(t, taille);
return v;
}

void Rayon::calcul_lcr(){
double lcr_form, Fm, rho;
lcr_form=0;
for(int i=0;i<nb_echant-1;i++) if(10*log(r[i])<seuil && 10*log(r[i+1])>seuil) lcr++;
lcr/=(dist/vit);
Fm=vit/lambda;
rho=fabs(seuil)/(sqrt(variance(loi_r, echelle))*sqrt(2));
lcr_form=sqrt(2*PI)*rho*Fm*exp(-rho*rho);
cout<<"lcr par simulation : "<<lcr<<"\n";
cout<<"lcr par calcul : "<<lcr_form<<"\n";
}

void Rayon::calcul_tau(){
double tau_form, Fm, rho;
tau_form=0;
for(int i=0;i<nb_echant;i++) if(10*log(r[i])<seuil) tau++;
tau/=Fe;
tau/=(lcr*dist/vit);
Fm=vit/lambda;
rho=fabs(seuil)/(sqrt(variance(loi_r, echelle))*sqrt(2));
tau_form=(exp(rho*rho)-1)/(Fm*rho*sqrt(2*PI));
cout<<"tau par simulation : "<<tau<<"\n";
cout<<"tau par calcul : "<<tau_form<<"\n";
}

void Rayon::ecrire(char *nom, double *t, int taille){
ofstream sortie;
sortie.open(nom);
for(int i=0;i<taille;i++) sortie<<*(t+i)<<"\n";
sortie<<moyenne(t, taille)<<"\n";
sortie<<sqrt(variance(t, taille))<<"\n";
sortie.close();
}

void main(){
Rayon ray(200,-107,20,100);
ray.calcul_IQ();
ray.calcul_loi_IQ();
ray.calcul_r();
ray.calcul_loi_r();
ray.calcul_lcr();
ray.calcul_tau();
}