Ave programmatori, sto scrivendo un programma che simuli il comportamento di un gas ideale e ne ricavi il coefficiente di diffusione (parte del programma è descritta nel testo "programmazione scientifica di M.Barone" a pg 319-329 nel caso qualcuno lo avesse).
Il programma deve soddisfare alcuni requisiti come il fatto che il reticolo contenente le particelle sia una mappa lineare e abbia condizioni periodiche al bordo inoltre devono essere stampate su file le posizioni iniziali di ogni particella quelle finali sul reticolo e quelle finali senza tener conto del reticolo(che poi verranno impiegate per ricavare il coefficiente di diffusione).
Il problema principale che ho riscontrato è quello di una inappropriata correlazione tra le posizioni iniziali delle particelle sul reticolo (come si può ben vedere dal grafico):
Qui ci sono le struct usate nelle funzioni:
typedef struct{
int x;
int y;
int z;
}coordinate;
//Definizione struct coordinate.
typedef struct{
coordinate p_0; /*posizione iniziale particella*/
coordinate p; /*posizione particella*/
coordinate tp; /*posizione reale particella*/
}particle;
//Definizione struct particle.
typedef struct{
int n; /*indicatore particella*/
int d; /*asse di movimento della particella*/
int m; /*verso particella*/
}scelta;
//Definizione struct scelta.
int l; /*lunghezza lato*/
int v; /*volume*/
int N; /*numero particelle*/
//Variabili globali.
Questa è la funzione che riempie il reticolo:
void riempimento(double dens, int *Reticolo, particle *P){
double random;
int r;
int n = 0;
int i;
for(i=0; i<v; i++){
Reticolo[i] = -1;
}
//Il reticolo è svuotato (-1 indica un volumetto vuoto).
while(n<N){
random = (double) rand()/RAND_MAX;
r = (int) (random*v); /*numero (pseoudo)casuale compreso tra 0 ed (v-1)*/
if(Reticolo[r] != -1){
Reticolo[r] = n;
P[n].p_0 = P[n].p = P[n].tp = trovacoord(r);
n++;
}
}
}
//Corpo funzione riempimento:riempie il reticolo e assegna le coordinate dell'array P[].
Qui c'è la funzione che converte l'indicatore del reticolo in un set di tre coordinate:
coordinate trovacoord(int i_v){
coordinate c;
c.z = i_v/(l*l); // c.z = parte intera di...
c.y = i_v/l - c.z*l; // c.y = parte intera di...
c.x = i_v - c.y*l - c.z*l*l;
return c;
}
//Corpo funzione findcord: converte l'indicatore del reticolo in un set di tre coordinate.
Normalmente cerco di rendere il codice "elegante", qui la situazione mi è decisamente sfuggita di mano, devo ammettere che quest'aggeggio oltre non fare il suo dovere è anche illeggibile, lo posto comunque nel caso qualcuno volesse farsi due risate
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
//Inclusione librerie.
#define D 3 /*numero dimensioni*/
//Define.
typedef struct{
int x;
int y;
int z;
}coordinate;
//Definizione struct coordinate.
typedef struct{
coordinate p_0; /*posizione iniziale particella*/
coordinate p; /*posizione particella*/
coordinate tp; /*posizione reale particella*/
}particle;
//Definizione struct particle.
typedef struct{
int n; /*indicatore particella*/
int d; /*asse di movimento della particella*/
int m; /*verso particella*/
}scelta;
//Definizione struct scelta.
int l; /*lunghezza lato*/
int v; /*volume*/
int N; /*numero particelle*/
//Variabili globali.
coordinate trovacord(int, int);
void riempimento(double, int*, particle*);
double s_quadratico(coordinate, coordinate);
scelta scelta_mossa ();
void mossa(particle*, int*);
//Dichiarazione funzioni.
/******************************************************/
int main (int argc, char* Argv[]){
if(argc != 2 && argc != 4 && argc != 6){
printf("\aInserire:\n lato(l)\n probabilita'riempimento (p)\n tempo\n passo t\n t finale.\n");
printf("Chiamare il programma con un parametro per maggiori informazioni\n");
exit(EXIT_FAILURE);
}
if(argc == 2){
printf("Questa applicazione .....");
exit(EXIT_FAILURE);
}
//Controllo numero argomenti inseriti.
srand(time(NULL));
//Inizializzazione seed per rand().
double dens; /*probabilita*/
double s_q_m = 0, c_diffusione; /*spostamento quadratico medio, coefficiente di diffusione*/
int stampa; /*se stampa e' pari a 1 viene generato un file contenente le posizioni di partenza e arrivo di ogni particella*/
register int n; /*indicatore particella*/
double t, delta_t, t_e, t_e_p, t_e_f; /*tempo, passo tempo, durata esperimento, passo durata esperimento, durata esperimento finale*/
register int i; /*variabile iterativa*/
int* Reticolo;
particle* P;
FILE *fp1, *fp2; /*file handlers*/
//Dichiarazione variabili.
do{
printf("stampare le posizioni iniziali e finali delle particelle? 0 = SI' 1 = NO\n ");
scanf("%i",&stampa);
}while(stampa != 0 && stampa != 1);
if(stampa == 0){
fp1 = fopen("Dati1.txt","w");
if(fp1 == NULL){
printf("\aERRORE APERTURA FILE!\n");
exit(EXIT_FAILURE);
}
fprintf(fp1,"# P_0 P TP\n");
fprintf(fp1,"# X Y Z X Y Z X Y Z\n");
}
fp2 = fopen("Dati2.txt","w");
if(fp2 == NULL){
printf("\aERRORE APERTURA FILE!\n");
exit(EXIT_FAILURE);
}
fprintf(fp2,"# t c\n");
//Apertura file di testo e scrittura intestazione.
l = strtod(Argv[1],NULL);
while (l == 0){
printf("Inserire l>0\n"); /*l NON puo' essere 0*/
scanf("%i",&l);
}
v =l*l*l; /*volume*/
dens = strtod(Argv[2],NULL);
N = floor(dens*v); /*numero particelle*/
delta_t = 1./N;
t_e = strtod(Argv[3],NULL);
if(argc == 6){
t_e_p = strtod(Argv[4],NULL);
t_e_f = strtod(Argv[5],NULL);
}
while(t_e_f<t_e && argc ==6){
printf("t_f deve essere maggiore di t!\n");
printf("Inserire t:\n");
scanf("%lf",&t_e);
printf("Inserire t_f:\n");
scanf("%lf",&t_e_f);
}
//Assegnazione e controllo argomenti inseriti.
Reticolo = (int*)malloc((v+1)*sizeof(int)); /*mappa lineare Reticolo[l][l][l]*/
P = (particle*)calloc(v+1,sizeof(particle));
//Allocazione dinamica di memoria per la mappa lineare Reticolo.
do{
riempimento(dens, Reticolo, P);
//Riempimento reticolo e assegnazione coordinate dell'array di particle.
t = 0.;
while(t<t_e){
mossa(P, Reticolo);
t += delta_t;
}
//Vengono eseguite tanti tentativi di muovere la particella quante sono le particelle.
t -= delta_t;
s_q_m = 0.;
for(n=0;n<N;n++){
s_q_m += s_quadratico(P[n].p_0,P[n].tp);
}
s_q_m = s_q_m/N;
c_diffusione = s_q_m/(2*D*t);
fprintf(fp2,"%12lf %3lf\n", t, c_diffusione);
//Calcolo coefficiente di diffusione.
if(stampa == 0){
fprintf(fp1, "# l = %i d = %lf t = %lf Numero particelle = %i\n",l, dens, t, N);
for(n=0;n<N;n++){
fprintf(fp1,"%6i %6i %6i",P[n].p_0.x, P[n].p_0.y, P[n].p_0.z);
fprintf(fp1," %6i %6i %6i",P[n].p.x, P[n].p.y, P[n].p.z);
fprintf(fp1," %6i %6i %6i\n",P[n].tp.x, P[n].tp.y, P[n].tp.z);
}
fprintf(fp1,"\n");
}
//Scrittura su file coordinate particelle se richiesto dall'utente.
if(argc == 6) t_e += t_e_p;
}while(t_e <= t_e_f && argc==6); /*in tal maniera se argc = 4 viene eseguito un solo ciclo*/
exit(EXIT_SUCCESS);
}
/******************************************************/
coordinate trovacoord(int i_v){
coordinate c;
c.z = i_v/(l*l); // c.z = parte intera di...
c.y = i_v/l - c.z*l; // c.y = parte intera di...
c.x = i_v - c.y*l - c.z*l*l;
return c;
}
//Corpo funzione findcord: converte l'indicatore del reticolo in un set di tre coordinate.
/****************************/
void riempimento(double dens, int *Reticolo, particle *P){
double random;
int r;
int n = 0;
int i;
for(i=0; i<v; i++){
Reticolo[i] = -1;
}
//Il reticolo è svuotato (-1 indica un volumetto vuoto).
while(n<N){
random = (double) rand()/RAND_MAX;
r = (int) (random*v); /*numero (pseoudo)casuale compreso tra 0 ed (v-1)*/
if(Reticolo[r] != -1){
Reticolo[r] = n;
P[n].p_0 = P[n].p = P[n].tp = trovacoord(r);
n++;
}
}
}
//Corpo funzione riempimento:riempie il reticolo e assegna le coordinate dell'array P[].
/****************************/
scelta scelta_mossa (){
double r;
scelta s;
r = (double)rand()/RAND_MAX;
s.n = (int) (r*N);
r = (double)rand()/(RAND_MAX+1.);
s.d = (int) (r*D+1);
r = (double)rand()/RAND_MAX;
if(r>0.5){
s.m = 1;
}else{
s.m = -1;
}
return s;
}
//Corpo funzione scelta_mossa restituisce le variabili per mossa(è chiamata internamente a move).
/****************************/
double s_quadratico(coordinate p_0, coordinate p_1){
double a,b,c; /*variabili di appoggio*/
double d; /*distanza*/
a = p_0.x - p_1.x;
b = p_0.y - p_1.y;
c = p_0.z - p_1.z;
d = a*a + b*b + c*c;
return d;
}
//Corpo funzione distanza: calcola la distanza tra due punti nel volume.
/****************************/
void mossa(particle *P, int *Reticolo){
int i_v0,i_v1; /*variabili di appoggio: indici posizione su Reticolo(vecchia e nuova)*/
scelta s = scelta_mossa(); /*chiamata funzione scelta mossa*/
particle P_nuovo = P[s.n]; /*variabile di appoggio (s.n è l'indice della particella scelta)*/
i_v0 = P_nuovo.p.x+P_nuovo.p.y*l+P_nuovo.p.z*l*l; /*P_nuovo NON e' ancora stato aggiornato*/
if(s.d == 1){
P_nuovo.tp.x = P_nuovo.p.x += s.m;
if(P_nuovo.p.x > l-1) P_nuovo.p.x = 0;
if(P_nuovo.p.x < 0) P_nuovo.p.x = l-1;
}
if(s.d == 2){
P_nuovo.tp.y = P_nuovo.p.y += s.m;
if(P_nuovo.p.y > l-1) P_nuovo.p.y = 0;
if(P_nuovo.p.y < 0) P_nuovo.p.y = l-1;
}
if(s.d == 3){
P_nuovo.tp.z = P_nuovo.p.z += s.m;
if(P_nuovo.p.z > l-1) P_nuovo.p.z = 0;
if(P_nuovo.p.z < 0) P_nuovo.p.z = l-1;
}
//Aggiornamento P_nuovo.
i_v1 = P_nuovo.p.x+P_nuovo.p.y*l+P_nuovo.p.z*l*l; /*indice Reticolo corrispondente alle coordinate di P_nuovo*/
if(Reticolo[i_v1] == -1){ /*se la nuova posizione sul reticolo è libera*/
P[s.n] = P_nuovo; /*la particella s.n cambia coordinate*/
Reticolo [i_v1] = s.n; /*il volume del reticolo corrispondente alle nuove cooridnate di P[s.n] e' occupato da s.n*/
Reticolo[i_v0] = -1; /*il volume del reticolo corrispondente alle vecchie coordinate di P[s.n] e' svuotato*/
}
}
//Corpo funzione move:tenta di muovere la particella n in direzione d con verso m,
//in caso di successo aggiorna il reticolo e le coordinate della particella.
/******************************************FINE PROGRAMMA**************************************/
Allegati:
Ottenuto con gnuplot,era tanto carino prima di cambiare taglia :(