[C] Formula Bethe-Bloch
Salve ragazzi ,
Premetto che non scrivo programmi in C da 2 anni ormai , ora ripassando ho scritto tale codice :
Il programma gira , il problema è che stampa valori indipendenti dai parametri che io immetto...
Dove sta l'errore?
Poi un altra domanda.. come si creano i grafici? vorrei farne uno che ha in ascissa i valori di $bethe_i$ e in ordinata quelli di $\beta\gamma$ come in figura :
E' un argomento che non ho mai trattato...
Grazie mille per l'eventuale aiuto
Premetto che non scrivo programmi in C da 2 anni ormai , ora ripassando ho scritto tale codice :
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define D 0.307075
#define m 0.510998902
#define c 300000000
int i,e;
long double Z,A,u,p,M,z,I,G,bethe;
int main ()
{
printf("Inserisci il numero atomico Z del bersaglio :\n");
scanf("%Lf",&Z);
printf("Inserisci la densita' del bersaglio in [g]/[cm]^3 :\n");
scanf("%Lf",&p);
printf("Inserisci la massa M della particella in [MeV]/[c]^2 :\n");
scanf("%Lf",&M);
printf("Inserisci il numero di massa A del bersaglio :\n");
scanf("%Lf",&A);
printf("Inserisci il numero di atomico z della particella :\n");
scanf("%Lf",&z);
printf("Inserisci il potenziale medio di eccitazione I :\n");
scanf("%Lf",&I);
printf("Inserisci l'energia della particella in [MeV] :\n");
scanf("%d",&e);
for (i=e;i>=0;i--) {
G=1/pow((1-(2*i)/(m*c*c)),0.5);
bethe=D*(Z/A)*p*z*z*log((2*m*c*c*((G*G)-1))/(I*(1+G*(m/M))));
printf("L'energia dispersa e': %Lg [MeV][cm]^2[g]^-1 \n",bethe);
}
getchar();
return 0;
}
Il programma gira , il problema è che stampa valori indipendenti dai parametri che io immetto...
Dove sta l'errore?
Poi un altra domanda.. come si creano i grafici? vorrei farne uno che ha in ascissa i valori di $bethe_i$ e in ordinata quelli di $\beta\gamma$ come in figura :
E' un argomento che non ho mai trattato...
Grazie mille per l'eventuale aiuto
Risposte
Ti rispondo sui grafici perché non ho tempo di leggere il codice 
Stampa i valori \((x,y)\) come righe di un file di testo [es: "x1 y1 \n x2 y2 ...."] e poi usa gnuplot per fare i grafici.
Stampa i valori \((x,y)\) come righe di un file di testo [es: "x1 y1 \n x2 y2 ...."] e poi usa gnuplot per fare i grafici.
Non si può fare senza l'uso di programmi esterni?
Dubito che tu possa farlo senza una libreria esterna, che di fatto è come usare un programma esterno.
Probabilmente esiste qualche libreria per farlo, ma usare gnuplot è molto più rapido (volendo lo puoi anche chiamare direttamente dal programma C usando system).
Venendo al tuo programma non conosco la formula che ti serve ma verosimilmente hai sbagliato qualcosa nella formula. Comunque gran parte della tua formula rimane costante tra i vari cicli. È comunque meglio usare sqrt che pow (a dire il vero esistono algoritmi piuttosto buoni per calcolare l'inverso di una radice quadrata).
Siccome io non so che formula vuoi implementare ti traduco il tuo codice in formule matematiche.
Hai che \(\displaystyle G(i) = \frac{1}{\sqrt{1 - \frac{2i}{mc^2}}} = \frac{\sqrt{mc^2}}{\sqrt{mc^2 - 2i}} \). Ovviamente \(\displaystyle mc^2 \) e \(\displaystyle \sqrt{mc^2} \) sono costanti.
Per quanto riguarda \(\displaystyle \mathrm{Bethe} \), riordinando leggermente le operazioni, si ha che \(\displaystyle \mathrm{Bethe} = \frac{DZpz^2}{A}\ln\biggl[ \frac{2Mmc^2(G^2-1)}{I(M + Gm)}\biggr] = \frac{DZpz^2}{A}\Biggl\{ \ln\biggl( \frac{2Mmc^2}{I}\biggr) + \ln\biggl[ \frac{(G^2-1)}{(M + Gm)}\biggr] \Biggr\}\) dove ho separato parte che cambia da parte costante.
È la funzione che volevi implementare?
Comunque a me, mettendo numeri assolutamente a caso, da risultati diversi ogni volta.
Venendo al tuo programma non conosco la formula che ti serve ma verosimilmente hai sbagliato qualcosa nella formula. Comunque gran parte della tua formula rimane costante tra i vari cicli. È comunque meglio usare sqrt che pow (a dire il vero esistono algoritmi piuttosto buoni per calcolare l'inverso di una radice quadrata).
Siccome io non so che formula vuoi implementare ti traduco il tuo codice in formule matematiche.
Hai che \(\displaystyle G(i) = \frac{1}{\sqrt{1 - \frac{2i}{mc^2}}} = \frac{\sqrt{mc^2}}{\sqrt{mc^2 - 2i}} \). Ovviamente \(\displaystyle mc^2 \) e \(\displaystyle \sqrt{mc^2} \) sono costanti.
Per quanto riguarda \(\displaystyle \mathrm{Bethe} \), riordinando leggermente le operazioni, si ha che \(\displaystyle \mathrm{Bethe} = \frac{DZpz^2}{A}\ln\biggl[ \frac{2Mmc^2(G^2-1)}{I(M + Gm)}\biggr] = \frac{DZpz^2}{A}\Biggl\{ \ln\biggl( \frac{2Mmc^2}{I}\biggr) + \ln\biggl[ \frac{(G^2-1)}{(M + Gm)}\biggr] \Biggr\}\) dove ho separato parte che cambia da parte costante.
È la funzione che volevi implementare?
Comunque a me, mettendo numeri assolutamente a caso, da risultati diversi ogni volta.
Le librerie standard del C sono molto limitate perché è un linguaggio pensato per essere usato per programmazione di sistema. Le librerie sono quindi limitate ad operazioni che possono essere utili a quasi ogni programma. Fare dei grafici non fa certamente parte di queste operazioni di uso generale. Altri linguaggi, come Matlab, che sono stati pensati fin dall'inizio con scopi di calcolo scientifico ovviamente posseggono librerie per questo tipo. Per fare cose di questo tipo in C hai quindi bisogno di una libreria o un programma esterno oppure di implementare tale libreria da te (può essere un interessante esercizio ma direi che per il momento è oltre le tue capacità).
Grazie mille a tutti per l'interesse 
@Raptorista : mmm a questo punto allora meglio che mi sposti all'interno di ROOT cosi da fare tutto in un unico programma .
@vic85 : la formula è esatta , inserendo valori a caso ottengo i primi n-2 valori uguali e gli ultimi 2 "-inf" , forse ti può essere d'aiuto , se compilo da Code::Blocks ottengo ciò che ho detto in precedenza , da terminale invece mi dice che la funzione va in overflow e credo che il problema stia qui.. dato che maneggia numeri dell'ordine di $10^16$..
@apatriarca : Ciò che dici è corretto e mi stuzzica però non credo di esserne capace e inoltre non ho nemmeno il tempo materiale..
@Raptorista : mmm a questo punto allora meglio che mi sposti all'interno di ROOT cosi da fare tutto in un unico programma .
@vic85 : la formula è esatta , inserendo valori a caso ottengo i primi n-2 valori uguali e gli ultimi 2 "-inf" , forse ti può essere d'aiuto , se compilo da Code::Blocks ottengo ciò che ho detto in precedenza , da terminale invece mi dice che la funzione va in overflow e credo che il problema stia qui.. dato che maneggia numeri dell'ordine di $10^16$..
@apatriarca : Ciò che dici è corretto e mi stuzzica
però non credo di esserne capace e inoltre non ho nemmeno il tempo materiale..
Ah, siamo nel dominio della fisica allora? 
Per il problema di overflow, ti consiglio di normalizzare tutti i valori rispetto ad una qualche grandezza di riferimento in modo che siano tutti dello stesso ordine di grandezza, possibilmente l'unità.
Quando hai finito i calcoli, riporta tutto in unità fisiche. In questo modo eviti che la numerica ti rovini le feste.
Per il problema di overflow, ti consiglio di normalizzare tutti i valori rispetto ad una qualche grandezza di riferimento in modo che siano tutti dello stesso ordine di grandezza, possibilmente l'unità.
Quando hai finito i calcoli, riporta tutto in unità fisiche. In questo modo eviti che la numerica ti rovini le feste.
Sono iscritto a Fisica 
Mmm ci proverò
Mmm ci proverò
"MillesoliSamuele":
mmm a questo punto allora meglio che mi sposti all'interno di ROOT cosi da fare tutto in un unico programma .
Fidati.. non è quello che vuoi. Hai semplicemente bisogno di installare gnuplot in modo che sia nel PATH e poi richiamarlo usando system oppure da linea di comando (che poi è in pratica la stessa cosa).
Come ti ho detto, scrivere una libreria per fare grafici non è semplice. Oltretutto conviene comunque anche in quel caso affidarti a librerie esterne per questioni come l'esportazione o la visualizzazione delle immagini che creerai. Insomma.. lascia perdere e usa GNUPLOT o qualcosa di simile. Oppure lascia perdere il C e usa qualcosa come Matlab, Scilab, python+scipy, Julia..
In Python ho già in programma
Guardando la tua formula vedo un importante problema:
\[G^2 - 1\]
Questa formula VA EVITATA! A basse energie si ha infatti che \( G \approx 1/\sqrt{1 - 10^{-16}} \approx 1. \) Quindi quella sottrazione porta ad un valore molto vicino a zero che viene poi moltiplicato per un valore molto grande \( 2\,m\,c^2\). L'equazione andrebbe riscritta in modo da non avere quel tipo di cancellazione. Non saprei darti consigli adesso, dovrei studiarla un po' meglio. Forse potrebbe ad un certo punto convenire non usare \(G\) e trovare invece modi diversi di estrarre costanti dalla tua formula. Per il resto direi che dipende dai valori che inserisci.
\[G^2 - 1\]
Questa formula VA EVITATA! A basse energie si ha infatti che \( G \approx 1/\sqrt{1 - 10^{-16}} \approx 1. \) Quindi quella sottrazione porta ad un valore molto vicino a zero che viene poi moltiplicato per un valore molto grande \( 2\,m\,c^2\). L'equazione andrebbe riscritta in modo da non avere quel tipo di cancellazione. Non saprei darti consigli adesso, dovrei studiarla un po' meglio. Forse potrebbe ad un certo punto convenire non usare \(G\) e trovare invece modi diversi di estrarre costanti dalla tua formula. Per il resto direi che dipende dai valori che inserisci.
Il tuo codice ha vari problemi.
Prima di tutto i problemi numeri segnalati da apatriarca. Dopo di che ci sono problemi meno seri come il fatto di memorizzare i valori costanti nel #define in formato double o usare log al posto di logl.
Venendo al problema numerico si ha che:
\(\displaystyle 2mc^2(G^2 - 1) = 2mc^2\biggl( \frac{mc^2}{mc^2 - 2i} - 1 \biggr) = 2mc^2\frac{2i}{mc^2 - 2i} = 4G^2i \). Che è un numero molto più gestibile numericamente. Usando le proprietà del logaritmo io comunque scriverei il tutto così:
Prima di tutto i problemi numeri segnalati da apatriarca. Dopo di che ci sono problemi meno seri come il fatto di memorizzare i valori costanti nel #define in formato double o usare log al posto di logl.
Venendo al problema numerico si ha che:
\(\displaystyle 2mc^2(G^2 - 1) = 2mc^2\biggl( \frac{mc^2}{mc^2 - 2i} - 1 \biggr) = 2mc^2\frac{2i}{mc^2 - 2i} = 4G^2i \). Che è un numero molto più gestibile numericamente. Usando le proprietà del logaritmo io comunque scriverei il tutto così:
#include <stdio.h>
#include <math.h>
#define D 0.307075L
#define m 0.510998902L
#define c 300000000L
int main(void)
{
long double Z, A, p, M, z, I;
int e;
printf("Inserisci il numero atomico Z del bersaglio :\n");
scanf("%Lf", &Z);
printf("Inserisci la densita' del bersaglio in [g]/[cm]^3 :\n");
scanf("%Lf", &p);
printf("Inserisci la massa M della particella in [MeV]/[c]^2 :\n");
scanf("%Lf", &M);
printf("Inserisci il numero di massa A del bersaglio :\n");
scanf("%Lf", &A);
printf("Inserisci il numero di atomico z della particella :\n");
scanf("%Lf", &z);
printf("Inserisci il potenziale medio di eccitazione I :\n");
scanf("%Lf", &I);
printf("Inserisci l'energia della particella in [MeV] :\n");
scanf("%d", &e);
long double const mmc2 = m*c*c;
long double const k1 = D * (Z/A) * p * z * z;
long double const k2 = logl(4*M) - logl(I);
for (int i = e; i >= 0; i--)
{
long double const G2 = mmc2/(mmc2 - 2.0*i);
long double const bethe = k1 * (k2 + logl(G2) + logl(i) - logl(M + sqrt(G2)*m));
printf("L'energia dispersa e': %Lg [MeV][cm]^2[g]^-1 \n", bethe);
}
}
@Apatriarca ; si mi ero accorto di ciò , ma gia ho semplificato abbastanza la formula per via di termini uguali a G , non volevo ulteriormente approssimarla..
@vict85 grazie per il codice! Non ci sarei mai arrivato! Durante il corso ( durato poche ore ) ci è stato illustrato solo la base del C.. quindi non conoscevo i metodi introdotti da te , ne farò tesoro
@vict85 grazie per il codice! Non ci sarei mai arrivato! Durante il corso ( durato poche ore ) ci è stato illustrato solo la base del C.. quindi non conoscevo i metodi introdotti da te , ne farò tesoro
Ma non stavo consigliando di approssimare la formula, solo di evitare di trovarti in situazioni in cui numeri molto simili vengono sottratti tra di loro per poi essere moltiplicati per valori molto grandi. Scrivendo la formula in questo modo stai infatti dando molta importanza ad eventuali approssimazioni durante il calcolo (dovute al fatto che stai lavorando con valori approssimati). Esistono spesso molti modi diversi di scrivere una formula, ognuna con vantaggi diversi. Puoi per esempio cercare di fare meno calcoli, oppure leggere più leggibile l'espressione o anche cercare di ottenere un risultato il più preciso possibile. In questo caso il risultato potrebbe non essere molto preciso.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo