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

using namespace std;

double c =2.99792458e8;     //vitesse lumière
double lambda = 632.8e-3;   //longueur d'onde du laser

int Npas;           // Discrétisation de l'espace
double Tmin, Tmax;

int conv(double y, double ymax){    // converti une position vers l'indice de discrétisation de l'espace
    return int(y*Npas/ymax);
}

void gradient(int Npas, double grad[]){   //initialisation du tableau de gradient avec une fonction
    int i;
    for(i=0;i<Npas;i++){
        grad[i]=Tmin+(Tmax-Tmin)*double(i)/double(Npas); //fonction que l'on peut changer
    }
}

double n(double y, double ymax, double grad[]){   // renvoie l'indice optique correspondant à la position entré
    return 1+(1e-6 *(288.02+1.478/lambda/lambda+0.0316/pow(lambda,4))/(1+0.003716*grad[conv(y,ymax)]));
}

void deriv(double ymax, double dt, double grad[],double y[], double dy[]){ //Equation du mouvement d'un " grain de laser "
    dy[0]=y[2];                                                 // x
    dy[1]=y[3];                                                 // z
    dy[2]=0;                                                    // vx, calculer après de manière analytique à partir de vz et c²/n²
    dy[3]=pow(y[2],2)/n(y[1],ymax,grad)*(n(y[1]+dy[1]*dt,ymax,grad)-n(y[1],ymax,grad))/(1e-100+dy[1]*dt);           // vz
}

void rk4(double grad[],double ymax,int m, double x, double y[] ,double dx, void deriv(double,double, double[],double[], double[])){
    // Algorythme de runge kutta d'ordre 4
    int i ;
    double ddx ;
    double d1[m], d2[m], d3[m], d4[m], yp[m];
    ddx = dx/2.;
    deriv(ymax,ddx,grad,y,d1) ;
    for( i = 0; i< m; i++){ yp[i] = y[i] + d1[i]*ddx ; }
    yp[2]=sqrt(c*c/n(yp[1],ymax,grad)/n(yp[1],ymax,grad)-yp[3]*yp[3]);  // calcul de la valeur de vx en fonction de la valeur de vy
    deriv(ymax,ddx,grad,yp,d2) ;
    for( i = 0; i < m; i++){ yp[i] = y[i] + d2[i]*ddx ; }
    yp[2]=sqrt(c*c/n(yp[1],ymax,grad)/n(yp[1],ymax,grad)-yp[3]*yp[3]);
    deriv(ymax,ddx,grad,yp,d3) ;
    for( i = 0; i< m; i++){ yp[i] = y[i] + d3[i]*dx ;}
    yp[2]=sqrt(c*c/n(yp[1],ymax,grad)/n(yp[1],ymax,grad)-yp[3]*yp[3]);
    deriv(ymax,ddx,grad,yp,d4) ;
    for( i = 0; i < m ; i++)
    { y[i] = y[i] + dx*( d1[i] + 2*d2[i] + 2*d3[i] + d4[i] )/6 ; }
    y[2]=sqrt(c*c/n(y[1],ymax,grad)/n(y[1],ymax,grad)-y[3]*y[3]);
}

int main(){
    double xmax = 0.6;
    double zmax = 0.1;

    //initialisation:

    cout << "Tmin (Celcius) (z=0) ?" << endl;
    cin >> Tmin;
    cout << "Tmax (Celcius) (z=zmax)?" << endl;
    cin >> Tmax;
    cout << "Proportion de la cuve : 10x10x60" << endl;
    Npas=100000;
    double theta;
    cout <<"Angle incident (degre) " << endl;
    cin >> theta;
    theta *= M_PI/180.;      //conversion en rad

    double dz = zmax/Npas;

    // Evolution rayon :

    double grad[Npas];  // n(z)
    gradient(Npas,grad);
    double y[4];
    y[0]=0; //x
    y[1]=zmax/2; //z
    y[2]=cos(theta)*c/n(y[1],zmax,grad);   //vx
    y[3]=sin(theta)*c/n(y[1],zmax,grad);   //vz
    double dt =1e-14;


    ofstream fichier("Air_1_nograd.res");     //on enregistre les différents paramètres dans un fichier
    for(double t=0;t<=2.2e-9;t+=dt)
    {

        if(y[0]>xmax){break;}
        if(y[1]>0.1-1.5*dz){y[3]=-y[3];}
        if(y[1]<1.5*dz){y[3]=-y[3];}
        fichier << t ;
        cout << t << endl ;       //affiche le temps (juste pour montrer que le programme tourne !
        for(int i=0;i<4;i++)    fichier << " " << y[i]  ;   //enregistre chaque paramètres
        fichier << endl;
        rk4(grad,zmax,4,t,y,dt,deriv);      //évolution des coordonnées
    }
    fichier.close();

}