Entry 3119

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define L 30 //lunghezza massima del loop
#define DELTA M_PI //pigreco math.h
#define CORSA 100000  //quanti cammini lunghi L voglio fare
#define CHUNK_SIZE 1000
#define SIGMA 0.02
 
inline void soluzione(const double *z,const double *r,double *v);
inline void croce(const double *x, const double *y, double *z);
inline double punto(const double *x,const double *y);
inline void normalizza(double *x);
void visualizza_rr(double rr[L][3]){
	int i,j;
	printf("\n\n\n");
	for(i=0;i<3;i++){
		printf("\n");
		for(j=0;j<L;j++)
			printf("%.4f\t",rr[j][i]);
	}
	printf("\n");
}

int main(){
	const gsl_rng_type *T;	 //dichiarazione tipo	
	gsl_rng *r;		//dichiarazione puntatore generatore
	gsl_rng_env_setup ();	//inizializzazione del generatore prendendo dalla bash le variabili: $ GSL_RNG_TYPE="..." GSL_RNG_SEED=... ./eseguibile   se fornite, altrimenti vedi default pag. 204 manuale gsl
			// ->	//in questo caso è ideale chiamare il programma con ad es GSL_RNG_SEED=date +%s  ./eseguibile
	T = gsl_rng_default;	//assegnazione tipo generatore
	r = gsl_rng_alloc (T);	//assegnazione generatore tipo T
	
	int chunk_size=CHUNK_SIZE;
	int i, j;
	double theta, phi, ct, st, cp, sp, tmp3, tmp4;
	/* scegliamo il primo vettore rr =asse y, il primo piano alfa =piano XY
	ad ogni passo k, alfa(k) è il piano passante per l'origine, che contiene il vettore rr(k)
	 ed il vettore v(k+1) perpendicolare ad rr(k) ed appartenete ad alfa(k-1)
	 z è l'asse di rotazione di theta; al primo passo se alfa=>XY rr=>asse y allora z=>asse z 
	 v giace sugli incroci fra alfa(k) e alfa(k-1) quindi è perpendicolare sia a z(k-1) che ad rr(k-1)
	 chiaramente alfa è determinato da z e viceversa quindi uso solo z
	*/
	double tmp1[3];
	double tmp2[3];
	double z[3]={0,0,1};
	double v[3];
	double rr[L][3], ss[L][3]; 	//standard C : nello stack per righe
	rr[0][0]=0;	//il primo vetore non viene mai toccato
	rr[0][1]=1;
	rr[0][2]=0;
	ss[0][0]=0;	//il primo vetore non viene mai toccato
	ss[0][1]=fabs(1.0+gsl_ran_gaussian(r, SIGMA));
	ss[0][2]=0;
	double m[3];
	double modulo;
	FILE *fp; 

    	fp=fopen("sequenza_dati.bin","wb");  //in questo file andremo ascrivere ogni sequenza lunga trenta vettori in formato binario

   	if(fp==NULL) return;
	
	setvbuf( fp , NULL , _IOFBF , 720 );  //questo dovrebbe forzare la scrittura su file ogno 720 byte che sono la dimensione della variabile rr
					      // cio dovrebbe essere utile usando open MP per evitare di incasinare l'ordine di scrittura sul file	


	
	for(i=1; i<= CORSA; i++){
	
		for(j=1; j<L; j++){
		
    			theta = -DELTA + 2 * DELTA * gsl_rng_uniform(r); //piatto fra -D e D
    			//printf("theta %.3f\t",theta);
    			ct = cos(theta);
    			st = sin(theta);
    			//printf("ct %.5f\t",ct);
    			//printf("st %.5f\t",st);
    			phi = 2 * M_PI * gsl_rng_uniform(r);     //piatto fra 0  e 2pi
    			cp=cos(phi);
    			sp=sin(phi);
    			
    			// 4 variabili temporanee per limitare i calcoli ripetuti
    			croce(z,rr[j-1],tmp1);
    			tmp3=punto(z,rr[j-1]);
    			
		 	//faccio le 2 rotazioni attorno ad un asse fissato
		 			 	
			        rr[j][0] = rr[j-1][0]*ct + tmp1[0]*st + tmp3 * (1-ct)*z[0];
				rr[j][1] = rr[j-1][1]*ct + tmp1[1]*st + tmp3 * (1-ct)*z[1];
			        rr[j][2] = rr[j-1][2]*ct + tmp1[2]*st + tmp3 * (1-ct)*z[2];
			 croce(rr[j-1],rr[j],tmp2);
			 tmp4= punto(rr[j],rr[j-1]);        
			        rr[j][0] = rr[j][0]*cp + tmp2[0]*sp + tmp4*(1-cp)*rr[j-1][0];
				rr[j][1] = rr[j][1]*cp + tmp2[1]*sp + tmp4*(1-cp)*rr[j-1][1];
				rr[j][2] = rr[j][2]*cp + tmp2[2]*sp + tmp4*(1-cp)*rr[j-1][2];
				
				
		    	//printf("tmp1 %.3f %.3f %.3f \t tmp2 %.3f %.3f %.3f \t tmp3 %.3f \ttmp4 %.3f \t",tmp1[0],tmp1[1],tmp1[2],tmp2[0],tmp2[1],tmp2[2],tmp3,tmp4);
		 	//printf("\n");
			normalizza(rr[j]); 	//questa ulteriore rinormalizzazione è fondamentale per controllare l'errore che si propaga
		    	
		    	//m[0]=gsl_ran_gaussian (r, SIGMA);
		    	//m[1]=gsl_ran_gaussian (r, SIGMA);
		    	//m[2]=gsl_ran_gaussian (r, SIGMA);
		    	//modulo=sqrt(punto(m, m));
		    	modulo=fabs(1+gsl_ran_gaussian (r, SIGMA));
		    	//printf("modulo %.5f\n",modulo);
				ss[j][0]=rr[j][0]*modulo;	
				ss[j][1]=rr[j][1]*modulo;
				ss[j][2]=rr[j][2]*modulo;

		    	
		    	//preparo i nuovi vettori
		    	soluzione(z,rr[j-1],v);	 // v giace sull' incrocio fra i due piani
        		normalizza(v);
        		croce(rr[j-1],v,z);
        	}
        	fwrite(ss,sizeof(double),90,fp); //90 è sizeof(rr)/sizeof(double)
        	//visualizza_rr(ss);

        } //fine sezione parallela
        fclose(fp);
   	return 0;
	
	
}

/////////////////////////////////////////////////////////////////////////////////
/*funzione che risolve con il metodo di kramer un sistema 3x3 NON singolare*/
/////////////////////////////////////////////////////////////////////////////////
inline void soluzione(const double *z, const double *r, double *v){
	double det_originale = z[0]*(r[1]-r[2]) - z[1]*(r[0]-r[2]) + z[2]*(r[0]-r[1]);
	v[0] = (z[1]*r[2] - r[1]*z[2]) / det_originale;
	v[1] = (r[0]*z[2] - z[0]*r[2]) / det_originale;
	v[2] = (z[0]*r[1] - r[0]*z[1]) / det_originale;

}

/////////////////////////////////////////////////////////////////////////////////
/*funzione che fa il prodotto a croce di 2 vettori*/
/////////////////////////////////////////////////////////////////////////////////
inline void croce(const double *x, const double *y, double *z){
	*z = x[1]*y[2]-y[1]*x[2];
	*(z+1) = x[0]*y[2]-y[0]*x[2];
	*(z+2) = x[0]*y[1]-y[0]*x[1];

}

/////////////////////////////////////////////////////////////////////////////////
/*funzione che restituisce il prodotto scalare*/
/////////////////////////////////////////////////////////////////////////////////
inline double punto(const double *x,const double *y){
	return x[0]*y[0] +x[1]*y[1] +x[2]*y[2];
}

/////////////////////////////////////////////////////////////////////////////////
/*funzione che normalizza un vettore di lunghezza 3 */
/////////////////////////////////////////////////////////////////////////////////
inline void normalizza(double *x){
	double norma=sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
	*x = *x / norma;
	*(x+1) = *(x+1) / norma;
	*(x+2) = *(x+2) / norma;
}