CALCUL DE SPIN

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

typedef float complex [2];

/* r = a + b */
void complex_add (complex a, complex b, complex *r)
{
	complex buffer = { a [0] + b [0],  a [1] + b [1]};

	*r [0] = buffer [0]; *r [1] = buffer [1];
}

/* r = a - b */
void complex_dif (complex a, complex b, complex *r)
{
	complex buffer = { a [0] - b [0],  a [1] - b [1]};

	*r [0] = buffer [0]; *r [1] = buffer [1];
}

/* r = a * b */
void complex_mul (complex a, complex b, complex *r)
{
	complex buffer = { a [0] * b [0] - a [1] * b [1],  a [0] * b [1] + a [1] * b [0]};

	*r [0] = buffer [0]; *r [1] = buffer [1];
}

/* r = sqrt (a) et a REEL !!! */
void complex_sqrt (float a, complex *r)
{
	if (a < 0)
	{
		*r [0] = 0.0f; *r [1] = (float)sqrt (a);
	}
	else
	{
		*r [0] = (float)sqrt (a); *r [1] = 0.0f;
	}
}

/* Structure pour caractériser les informations d'un spin */
typedef struct TSpin_
{
	float s;
	float *m;
	complex **mat_sx, **mat_sy, **mat_sz, **mat_ss;
} TSpin;


/* Pour libérer une strucure TSpin */
void detruire_spin (TSpin *data)
{
	free ((void*) (data)->m);
	free ((void*) (data)->mat_sx);
	free ((void*) (data)->mat_sy);
	free ((void*) (data)->mat_sz);
	free ((void*) (data)->mat_ss);
}

/* Pour créer une matrice de complex */
complex ** __get_matrice (int t)
{
	int i, j;
	complex **buffer;
	
	buffer = (complex**) malloc (sizeof(complex*)*t);
	
	for (i = 0; i < t; i ++)
	{
		buffer [i] = (complex*) malloc (sizeof(complex)*t);

		for (j = 0; j < t; j ++)
		{
			buffer [i] [j] [0] = 0.0f;
			buffer [i] [j] [1] = 0.0f;
		}
	}

	
	return buffer;		
}

/* Pour créer une structure TSpin */
void alloue_spin (TSpin *data, float s)
{
	int t = (int)2*s+1;
		
	data->s = s;
	data->m = (float*) malloc (sizeof(float)*t);
	
	data->mat_sx = __get_matrice (t);
	data->mat_sy = __get_matrice (t);
	data->mat_sz = __get_matrice (t);
	data->mat_ss = __get_matrice (t);
}

/* Pour afficher un nombre complex */
void __affiche_complex (complex a)
{
	if (a [0] == 0.0f && a [1] != 0.0f)
		printf ("%01f i", a [1]);
	else if (a [1] == 0.0f)
		printf ("%02f", a [0]);
	else
		printf ("%02f + %02f i", a [0], a [1]);
}

/* Pour afficher une matrice de complex */
void __affiche_matrice (complex **data, int t)
{
	int i, j;
	
	for (i = 0; i < t; i ++)
	{
		for (j = 0; j < t; j ++)
		{ __affiche_complex (data [i] [j]); printf ("\t");}
		
		printf ("\n");
	}		
}

/* Pour afficher une structure spin */
void affiche_spin (TSpin data)
{
	int i, t = (int)2*data.s+1;
	
	printf ("Spin %02f :\n\n", data.s);
	
	printf ("Base m = ");
	
	for (i = 0; i < t; i ++)
		printf ("|%02f>  ", data.m [i]);
	
	printf ("\n\nMatrice S^2 :\n\n");
	__affiche_matrice (data.mat_ss, t);
	
	printf ("\n\nMatrice Sx :\n\n");
	__affiche_matrice (data.mat_sx, t);
	
	printf ("\n\nMatrice Sy :\n\n");
	__affiche_matrice (data.mat_sy, t);
	
	printf ("\n\nMatrice Sz :\n\n");
	__affiche_matrice (data.mat_sz, t);
	
	printf ("\n\n");
}

/* Pour écrire un nombre complex dans un fichier */
void ecrire_complex (complex a, FILE *fichier)
{
	if (a [0] == 0.0f && a [1] != 0.0f)
		fprintf (fichier, "<font color=red><b>%01f i</b></font>", a [1]);
	else if (a [1] == 0.0f)
		if (a [0] == 0)
			fprintf (fichier, "%02f", a [0]);
		else
			fprintf (fichier, "<font color=red><b>%02f</b></font>", a [0]);
	else
		fprintf (fichier, "<font color=red><b>%02f + %02f i</b></font>", a [0], a [1]);
}

/* Pour écrire une matrice dans un fichier */
void ecrire_matrice (complex **data, int t, FILE *fichier)
{
	int i, j;
	
	fprintf (fichier, "<br><br><table border=1 width=90% align=center noWrap>");
	for (i = 0; i < t; i ++)
	{
		fprintf (fichier, "<tr>");
		for (j = 0; j < t; j ++)
		{ 
			fprintf (fichier, "<td noWrap>");
			ecrire_complex (data [i] [j], fichier); 
			fprintf (fichier, "</td>");
		}
		
		fprintf (fichier, "</tr>");
	}
	fprintf (fichier, "</table><br><br>");		
}


/* Pour écrire dans le fichier file_name une structure spin */
void sauver_spin (TSpin data, char *file_name)
{
	FILE *fichier;
	int i, t = (int)2*data.s+1;
	
	fichier = fopen (file_name, "w");
	
	fprintf (fichier, "<title>Etude de spin : generee par Decaux Thomas</title>");
	
	fprintf (fichier, "<br><br><h2>Etude Spin %02f :</h2>", data.s);
	
	fprintf (fichier, "Base m = ");
	
	for (i = 0; i < t; i ++)
		fprintf (fichier, "<b>|%02f> </b>", data.m [i]);
	
	fprintf (fichier, "<br><br>Matrice S^2 :\n\n");
	ecrire_matrice (data.mat_ss, t, fichier);
	
	fprintf (fichier, "\n\nMatrice Sx :\n\n");
	ecrire_matrice (data.mat_sx, t, fichier);
	
	fprintf (fichier, "\n\nMatrice Sy :\n\n");
	ecrire_matrice (data.mat_sy, t, fichier);
	
	fprintf (fichier, "\n\nMatrice Sz :\n\n");
	ecrire_matrice (data.mat_sz, t, fichier);
	
	fclose (fichier);
}

/* Pour calculer les informations d'un spin */
void calcul_spin (TSpin *spin)
{
	int i, j,  t = (int)2*(spin->s) + 1;
	complex **sp, **sm, demi = {0.5, 0}, ii = {2, 2};
	
	// Vecteurs de base
	spin->m [0] = -spin->s;
	
	for (i = 1; i < t; i ++)
		spin->m [i] = spin->m[i-1] + 1;
	
	// Calcul de S²
	for (i = 0; i < t; i ++)
		spin->mat_ss [i][i][0] = spin->s*(spin->s + 1);
	
	// Calcul de Sz
	for (i = 0; i < t; i ++)
		spin->mat_sz [i][i][0] = spin->m [i];
	
	// Calcul de S+
	sp = __get_matrice (t);
	for (i = 0; i < t-1; i ++)
		complex_sqrt (spin->s*(spin->s+1) - spin->m [i] * (1+spin->m [i]), &sp [i][i+1]);
	
	// Calcul de S-
	sm = __get_matrice (t);
	for (i = t-1; i > 0; i --)
		complex_sqrt (spin->s*(spin->s+1) - spin->m [i] * (1-spin->m [i]), &sp [i][i-1]);
	
	// Calcul de Sx = (S+ + S-) / 2
	for (i = 0; i < t; i ++)
		for (j = 0; j < t; j ++)
		{
			complex_add (sp [i][j], sm [i] [j], &spin->mat_sx [i] [j]);
			complex_mul (spin->mat_sx [i] [j], demi, &spin->mat_sx [i] [j]);
		}
		
	// Calcul de Sy = (S+ - S-) / 2*i
	for (i = 0; i < t; i ++)
		for (j = 0; j < t; j ++)
		{
			complex_dif (sp [i][j], sm [i] [j], &spin->mat_sy [i] [j]);
			complex_mul (spin->mat_sy [i] [j], demi, &spin->mat_sy [i] [j]);
			complex_mul (spin->mat_sy [i] [j], ii, &spin->mat_sy [i] [j]);
		}

	free ((void*) sp);
	free ((void*) sm);
}

int main (int argc, char *argv [])
{
	TSpin spin;
	float jaja;
		
	if (argc == 2)
		jaja = (float)atof (argv [1]);
	else
		jaja = 4;
	
	alloue_spin (&spin, jaja);
	
	calcul_spin (&spin);
	
	affiche_spin (spin);
	
	if (argc >= 3)
		sauver_spin (spin, argv [2]);
	sauver_spin (spin, "jaja.html");
	
	detruire_spin (&spin);

	getchar ();
	
	return 0;	
}