#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define DT 0.001

double deltaE(double fE_0);
double deltaY(double fL);
void step(double *fdeltaE, double *fdeltaY, double fE_0, double fL);

int main(){
    printf("il programma simula una particella che entra in un quadrato di lato l deciso dall-utente con energia iniziale decisa dall'utente e poi ad ogni passo toglie un energia delta E rispetto a quella data e faccia variare la posizione su y positivamente o negativamente \n \n ");
    double t=0, y=0, E_0, L, deltaE=0, deltaY=0,y_maxx=0, temp;
    int seed, j;
    FILE *fp;
    seed=time(0);
    srand48(seed);
    printf("immettere valore positivo di E_0 minore di 10 MeV \n");
    scanf("%lf", &E_0);
    while(E_0<0||E_0>10){
        printf("errore! reimmettere valore positivo di E_0 minore di 10 MeV \n");
        scanf("%lf", &E_0);
    }
    temp=E_0;
    printf("immettere valore positivo di L minore di 17 cm \n");
    scanf("%lf", &L);
    while(L<0||L>17){
        printf("errore! reimmettere valore positivo di L minore di 17 cm \n");
        scanf("%lf", &L);
    }
    fp=fopen("risultati.txt", "w");
    temp=E_0;
    for(j=0;j<1000;j++){
        while((y<L)&&(E_0>=0.5)){
            step(&deltaE, &deltaY, E_0, L);
            y+=deltaY;
            E_0-=deltaE;
            t+=DT;
            if(y>y_maxx){y_maxx=y;}
        }
        fprintf(fp, "alla %d* propagazione erano rispettivamete y.max= %lf tempo finale= %lf ed energia finale= %lf \n", j+1, y_maxx, t, E_0);
        y=0; E_0=temp; t=0;
    }
    fclose(fp);
}

double deltaE(double fE_0){
    double f, deltaE;
    f=0.03+0.7*(lrand48()/(double)RAND_MAX);
    deltaE=f*fE_0;
    return deltaE;
}

double deltaY(double fL){
    double v, deltaL;
    deltaL=0.05*fL+(0.03*fL)*((lrand48()/(double)RAND_MAX));
    v=(lrand48()/(double)RAND_MAX);
    if(v>0.5){
        deltaL=-deltaL;
    }
    return deltaL;
}

void step(double *fdeltaE, double *fdeltaY, double fE_0, double fL){
    *fdeltaE=deltaE(fE_0);
    *fdeltaY=deltaY(fL);
}