#include <cmath>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

using namespace std;

/*--------------------------------------------------------------------------------------------
--------------------------------------- Presentation : ---------------------------------------
Ce programme permet de simuler une perturbation dans un milieu 2D d'incice optique inhomogene.
Le milieu considere est de l'eau avec un gradient d'indice dependant uniquement de la hauteur.
La perturbation sera represente par une onde module par une fonction gaussienne se deplacant.
----------------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------*/

/*---------------------------- Liste des variables globales : ------------------------------*/
/*------------------------------- Caracteristiques eau : -----------------------------------*/
    double xmax = 0.6;                                  //- longueur de la cuve
    double zmax = 0.1;                                  //- hauteur maximale de la cuve
    double n_0 = 1.34;                                  //- indice optique en z = 0
    double n_zmax = 1.32;                               //- indice optique en zmmax
/*--------------------------- Caracteristiques Perturbation : ------------------------------*/
    /* Position initiale perturbation (aleatoire) : */
    double xp;
    double zp;
    /* Longueur d'onde perturbation : */
    double lambda;
    /* Modulation gaussienne : */
    double amplitude_i;                                 //- amplitude initiale
    double sigma_i;                                     //- ecart type initial
    //double sigma_i_x;                                 //- ecart type initial en x
    //double sigma_i_z;                                 //- ecart type initial en z
    /* Variation temporelle modulation : */
    double v;                                           //- vitesse deplacement de la modulation
    double alpha;                                       //- vitesse d'evolution de l'ecart type
    //double alpha_x;                                   //- vitesse d'evolution de l'ecart type en x
    //double alpha_z;                                   //- vitesse d'evolution de l'ecart type en z
/*------------------------------------------------------------------------------------------*/

double alea(double a, double b){return (((rand()/((double)RAND_MAX+1))*(b-a))+a);}  /* Renvoie un double compris entre a et b avec une probabilite uniforme. */
double r(double x, double z){return sqrt(x*x+z*z);}

double n_z(double z){                               /* Renvoie uniquement la partie stationnaire du gradient d'indice optique (dependant de la position */
    return  n_0+(n_zmax-n_0)*z;                     /* Modele lineaire */
}

double sigma(double t){return sigma_i+alpha*sqrt(t);}                       /* Etalement de la gaussienne */
double amplitude(double t){return amplitude_i/(sigma(t)*sqrt(2*M_PI));}     /* Ampltude totale de la gaussienne */
double R(double t){return v*t;}                                             /* Deplacement de la gaussienne */

double n(double x, double z, double t){                                     /* Renvoie l'indice optique à la position (x,z) à l'instant t */
    return n_z(z)+cos(2.*M_PI/lambda*sqrt(fabs(pow(x-xp,2)+pow(z-zp,2)-pow(R(t),2))))*amplitude(t)*exp(-(fabs(pow(x-xp,2)+pow(z-zp,2)-pow(R(t),2)))/(2*pow(sigma(t),2)));
    /* Somme du gradient et de la perturbation : un cos de longueur d'onde lambda, module par une fonction gaussienne s'etalant au cours du temps, tout deux se deplacant dans la cuve */
}

int main(){
    for(int i=0;i<75;i++){double a=alea(0,1);}                  /* Permet d'aller plus loin dans la graine generatrice de nombre aleatoire */
    /*-------------------------- Initialisation Perturbation : -----------------------------*/
    xp = alea(0.1,0.4);                                         /* - position aleatoire de la perturbation */
    zp = alea(0.025,0.075);

    lambda = alea(0.0001,0.005);                                /* - longueur d'onde de la perturbation */

    sigma_i = alea(0.0005,0.005);                               /* - ecart type de la modulation de la perturbation */
    //double sigma_x_i = alea(0.001,0.01);
    //double sigma_z_i = alea(0.001,0.01);
    amplitude_i = n_z(zp)*alea(0.00001,0.001);                  /* - amplitude initiale de la perturbation aleatoire entre 0.001% et 0.1% de l'indice optique à l'origine de la perturbation */

    v = alea(0.01,0.05);                                        /* vitesse de deplacement de la perturbation, entre 1 cm/s et 5 cm/s */
    alpha = alea(0.001,0.1);                                    /* vitesse d etalement de la gaussienne */
    //alpha_x = alea(0.01,1);                                   /* vitesse d etalement de la gaussienne */
    //alpha_z = alea(0.01,alpha_x);                             /* vitesse d etalement de la gaussienne */

    cout << "xp = " << xp << endl;
    cout << "zp = " << zp << endl;
    cout << "lambda = " << lambda << endl;
    cout << "sigma_i = " << sigma_i << endl;
    cout << "amplitude_i = " << amplitude_i << endl;
    cout << "v = " << v << endl;
    cout << "alpha = " << alpha << endl << endl;

    int j=0; // variable pour les sorties de fichiers
    ofstream fichier;

    /*------------------------------------ Traitement : --------------------------------------*/

    for(double t=0; t<=1; t+=0.1){

        /* Creation du nom de fichier : */
        string deb("indice."), fin(".res");
        char buf[21];
        sprintf(buf, "%i", j);
        string filename(deb+buf+fin);
        /* ouverture du fichier j : */
        fichier.open(filename.c_str());

        for(double x=0;x<=xmax;x+=0.0005){
            for(double z=0;z<=zmax;z+=0.0005){
                fichier << x << " " << z << " " << n(x,z,t) << endl;
            }
        }
        fichier.close();
        cout << t << endl;
        cout << sigma(t) << endl;
        cout << amplitude(t) << endl;
        j++;
    }
}

reset

set size ratio 0.1666

set xrange[0:59.9]
set yrange[0:9.9]

set xlabel "position x (cm)"
set ylabel "position z (cm)"
set title "Evolution du gradient d'indice optique dans la cuve - t=1s"

set pm3d map
splot 'indice.10.res' u ($1*100):($2*100):3 with image

set term png
set output "Grad.10.png"
replot