#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

/* Les longueurs sont données en m, la résistance en ohm, la capacité en farad, la tension en volt,
 l'intensité en amperes et les masses en kg, les vitesses en m/s et les energies en joules*/

//DEBUT DU MAIN
int main()

{
    //////////////////////////DECLARATIONS//////////////////////////
    
    double L, v; //Position et vitesse
    double a, aa; //Accélérations et accélérations au temps précédent
    double F, t; //projection des forces et temps
    int npas;
    
    /////////////////////////CONSTANTES ET PARAMETRES//////////////////////////
    
    int i;
    
    double I; //Intensité dans les plaques, Intensité dans la bille
    double u0 = 4*3.14*pow(10,-7); //Perméabilité magnétique du vide
    
    double R = 1*pow(10,-3); //Rayon de la bille
    double m = 2700*4*3.141592654*R*R*R/3;//masse de la bille
    double l = 0.15; //Position initiale de la bille
    
    double h = 1*pow(10,-2); //Hauteur des plaques
    double Lmax = 0.20; //Longueur des plaques
    double N = 4; //Nombre de câbles dans une plaque

    double Ra = 0.0008;  //Resistance
    double C = 1500*pow(10,-6), U = 250.0; //Capacité et tension d'un condensateur
    double Nc = 15; //Nombre de condensateurs
    double dt = C*Ra*0.001; //Pas de temps
    
    double Coef, log1, log2, b1;
    
    //CONDITIONS INITIALES
	L = l;
    v = 0.0;
    a = 0.0;
    aa = 0.0;
    t = 0.0;
    npas = 0;
    
	//BOUCLE PRINCIPALE
    
    while (L <= Lmax) {
        
        ////////CALCUL DES NOUVELLES POSITIONS ET MISE A ZERO DES FORCES
        L = L + v*dt + (2/3)*a*dt*dt - (1/6)*aa*dt*dt;
        F=0.0;
		
        ////////CALCUL DES NOUVELLES FORCES
        //Calcul de l'intensité au temps t dans les plaques
        I = (U/Ra)*exp(-t/(Ra*Nc*C))/N; 
        
        
        Coef = (u0*I*I)/(2*3.14*N);
        
        for(i=1;i<=N;i++){
            b1 = - R - (2*i - 1)*(h/2);
            log1 = log(  fabs( (sqrt((R-b1)*(R-b1) + L*L) - L)/(R-b1) ) );          
            log2 = log(  fabs( (sqrt((R+b1)*(R+b1) + L*L) - L)/(R+b1) ) );
            F += Coef*(log1-log2);
        }
        
        ////////CALCUL DES NOUVELLES VITESSES
        v = v + (1/(3*m))*F*dt + (5/6)*a*dt - (1/6)*aa*dt;
        
        //ACTUALISATION DES ACCELERATIONS
        aa = a;
        a= F/m;
        
        
        //Calcul du temps ecoule 
        t = t + dt;
        npas ++;    
        
    } //FIN BOUCLE PRINCIPALE  
    
    printf("V finale: %lf m/s, temps de parcours: %lf ms\n",v,t*1000);
    printf("%d",npas);
	
}//FIN DU MAIN