/*
	Programme pour calculer la forme d'un miroir
	d'apres les mesures a l'appareil de Foucault.
	
	Ce prog demande la focale du miroir,
	et les tirages mesures a l'appareil de Foucault
	pour chaque zone de l'ecran de Couderc.
	Il trouve automatiquement la constante a retrancher aux mesures
	pour ajuster le foyer au cercle de moindre aberration.
	
	La forme du miroir apparait a l'affichage lors de l'execution
	du programme, mais avec peu de details.
	
	Le prog ecrit en clair la coupe du front d'onde dans un fichier.
	Ces valeurs peuvent servir pour afficher plus finement la forme
	du miroir par une courbe en les copiant/collant dans un grapheur.
	
	Le fichier disque contenant les resultats du calcul est cree
	dans le dossier local lors de l'execution du programme.
	Son contenu est remis a jour s'il existe deja.
	Le nom du fichier cree est : "coupe.txt"
Historique:	
	9 1989  version init ecrite avec Think C sur Mac par L.Koechlin
	9 2002  LK reprise C++ Code Warrior Mac et PC et sortie fichier
*/

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

// ---------- Constantes a ajuster a votre situation -----------
const float DIAMETRE = 0.208;  // toutes les valeurs sont en metres

const float R_ZONE1 =   37.3e-3; // rayon externe de chaque zone
const float R_ZONE2 =   64.5e-3;
const float R_ZONE3 =   83.5e-3;
const float R_ZONE4 =  104.0e-3;
// -------------------------------------------------------------



class miroir {
	float diam;    // diametre du miroir
	float foc;     // focale du miroir
	int nb_zones;  // nb de zones decoupees dans ecran devant miroir
	float lambda ; // longueur d'onde utilisee pour le calcul
	
	float *hm;         // rayon (hauteur) du milieu d'une zone
	float *bord_zone;  // rayons bords zones sur decoupes ecran Couder
	float *hm2_sur_r;  // loi "de la parabole" : p = hm^2 / r
	float *mesures;    // valeurs mesurees au Foucault
	float *ecarts;     // differences : mesures - hm^2 / r
	float *aberr_trans;// aberrations transverses calculees au foyer

public:
	int   nb_pt_h; // nb pts calcules sur onde le long du r
	float resol_h;    // taille en metres de l'element de resolution
	float *onde;      // coupe radiale du front d'onde calcule
	float *onde_pv;   // coupe radiale fr d'onde ajust pour Pic-Vallee
	float onde_min, onde_max;  // defauts du front d'onde
// ------------------
	miroir (void);
	void calcul_aberr_trans (void);
	void calcul_onde (void);
	void enreg_onde (void);
	void imprime_rustic_onde (void);
};


/* -------------------------- main -------------------------- */
main ()
{
	miroir *le_miroir = new miroir ();
	le_miroir->calcul_aberr_trans ();
	le_miroir->calcul_onde ();
	le_miroir->enreg_onde ();
	le_miroir->imprime_rustic_onde ();
	return 0;
}



/* --------------------- methodes miroir --------------------- */
miroir::miroir (void)
{
	nb_pt_h = 76;      // nb de points echantillonnes sur rayon miroir
	lambda = 0.55e-6;  // longueur d'onde moyenne pour le visible
	diam = DIAMETRE;   // du miroir
	printf ("\nProgramme ajustŽ pour un miroir de diametre %3.0f mm\n",
	diam*1000);
	
	printf ("Entrer la focale de votre miroir en millimetres\n");
	printf ("(par exemple 1220)\n");
	scanf ("%f", &foc);    // lire la focale en mm
	foc /= 1000.;          // la mettre en metres
	getc (stdin);
	
	nb_zones = 4;  // changer cette valeur si ecran a plus de 4 zones
	bord_zone = new float [nb_zones];
	bord_zone [0] = R_ZONE1;
	bord_zone [1] = R_ZONE2;
	bord_zone [2] = R_ZONE3;
	bord_zone [3] = R_ZONE4;
	// ajouter ici eventuellement les rayons des zones supplementaires

	hm =        new float [nb_zones];
	hm2_sur_r = new float [nb_zones];
	mesures =   new float [nb_zones];
	ecarts =    new float [nb_zones];
	aberr_trans = new float [nb_zones];

	resol_h = diam/2 / (float)nb_pt_h; // le pas de l'integration
	onde = new float [nb_pt_h +1];     // contiendra la coupe radiale
	onde_pv = new float [nb_pt_h +1];  // meme coupe mais refocalisee
}

// ----------- ajustage foyer et aberrations transverses --------------
/*
Cette premiere partie trouve le tirage du meilleur foyer
et calcule les aberrations transverses.
Elle donne ensuite les defauts approximatifs du front d'onde.
On peut tracer a la main le dessin a partir
des quatre valeurs donnes dans la colone "defauts approx"
*/
void miroir::calcul_aberr_trans (void)
{
	float x_minaber; //abscis crcle moindr aber: cst a retr.aux msures
	int i, j, k;
	float r = 2*foc;		// rayon de courbure du miroir

	printf ("Entrer les mesures en mm pour les zones 1, 2, 3, 4\n");
	printf ("(par exemple  0.27  0.65  2.38  3.65)\n");
	for (i = 0 ; i < nb_zones ; i++)
	{
		scanf ("%f", &(mesures [i]));  // lire mesures entrees en mm
		mesures [i] /= 1000.;          // mettre les mesures en metres
		getc (stdin);
	}
	// calcul du rayon moyen des zones a partir de leur rayon externe
	hm [0] = bord_zone [0] /1.414;
	for (i = 1 ; i < nb_zones ; i++)
			hm [i] = (bord_zone [i] + bord_zone [i-1]) /2;

	for (i = 0 ; i < nb_zones ; i++)
	{
		hm2_sur_r [i] = hm [i] * hm [i] / r;      // calcul hm^2/r
		ecarts [i] = mesures [i] - hm2_sur_r [i]; // calcul differnces
	}	
	/*
Trouver le cercle de moindre aberration :
sa position est donnee par la seule intersection de 2 rayons
qui soit a l'exterieur de tous les autres rayons.
des rayons sont "opposes" s'ils viennent de cotes opposes du miroir.
	 */
	// D'abord trouver l' intersec de chaque rayon avec les "opposes"
	int minaber;  // Boolean
	float x, y;
	for (i = 0 ; i < nb_zones ; i++)
	{
		for (j = i+1 ; j < nb_zones ; j++)
		{
			x = (ecarts [i] * hm [i] + ecarts [j] * hm [j])
                                        / (hm [i] + hm [j]);
			y = fabs (hm [i] * (x - ecarts [i]));
			
			// chercher si elle est a l'exterieur des autres rayons
			minaber = true;
			for (k = 0 ; k < nb_zones ; k++)
			{
				if ((k != i) && (k != j))
				{
					if (y < fabs (hm [k] * (x - ecarts [k])))
					{
						minaber = false; // Il y a un rayon +loin
						break;
					}
				}
			}
			if (minaber)  { x_minaber = x; goto fin; }
		}
	}
fin : ;

	// Calcul des aberrations transverses
	for (i = 0 ; i< nb_zones ; i++)
		aberr_trans[i] = -1 * (ecarts[i] - x_minaber) * hm[i] / (2*r);	
	/*
Calcul approx des defauts lineaires sur l'onde. 
Le calcul ci-dessous n'est pas une vraie integ, il est approximatif.
Les valeurs affichees ici correspondent a celles de bulletin de controle
de Texereau.
	*/
	float *onde_zigzag = new float [nb_zones]; // estimation de l'onde
	onde_zigzag [0] = bord_zone [0] * aberr_trans [0] /foc;
	for (i = 1 ; i< nb_zones ; i++)
		onde_zigzag [i] =
				onde_zigzag [i-1]
				+ ((bord_zone [i]
				- bord_zone [i-1])  * aberr_trans [i] /foc);
	// lister les resultats
	printf ("\nPour lambda =%4.0f nm, rayon de diffraction =  %3.2e\n",
	 lambda * 1e9, 1.22 * lambda * foc / diam);
	printf ("                                                  m\n");
	printf ("zone bord  milieu  loi  mesure  ecart aberr-L   aber-T\
   pente  onde approx\n");
	printf ("      mm     mm     mm    mm     mm      mm       m   \
                 m\n");
	printf ("---  ----  ------  ---  ------  ----- ------- --------\
   ------  ----------\n");
	for (i = 0 ; i< nb_zones ; i++)
printf("%1d %6.1f %6.1f %6.2f %6.2f %6.2f %6.2f %10.2e %10.2e %9.2e\n",
			i+1, bord_zone [i]*1000,
			hm [i]*1000,
			hm2_sur_r [i]*1000,
			mesures [i]*1000,
			ecarts [i]*1000,
			(ecarts [i] - x_minaber)*1000,
			aberr_trans [i],
			aberr_trans [i]/foc,
			onde_zigzag [i]);
	printf("\n");
}
	
// ------- Integration plus fine des defauts de front d'onde --------
/*
A partir des aberrations transverses calculees ci_dessus, cette
routine calcule la coupe radiale du front d'onde.
Dans le tableau onde [], le front d'onde est calcule pour un foyer
ajuste sur le cercle de moindre aberration.
Dans le tableau onde_pv [], le front d'onde est calcule pour un foyer
ajuste a une surface de reference en contact avec le front d'onde
a la fois au centre et au bord du miroir. Ceci permet de calculer
le defaut max "pic-vallee" du front d'onde.
*/
void miroir::calcul_onde (void)
{
	int i;		
	
	float aber_zn_inf = 0;	             // aber tr. nulle au centre
	float aber_zn_sup = aberr_trans [0]; // init aber tr. de zone 1

	float hm_zone_inf = 0;               // centre du miroir	
	float hm_zone_sup = hm[0];           // rayon moyen zone 1
	
	
	onde [0] = 0; // cale le defaut sur l'onde (piston) a 0 au centre
	int zone = 0; // initialement on est dans la zone 1
	
	// interp (extrapoler quand necessaire) les aberr transverses
	for (i = 0 ; i<nb_pt_h ; i++)
	{
		float h = i*resol_h;    // hauteur d'un rayon
		if ((h > hm [zone]) && (zone < nb_zones-1))
		{
			zone++;   // on entre dans la zone suivante
			aber_zn_inf = aber_zn_sup;
			aber_zn_sup = aberr_trans [zone];
			
			hm_zone_inf = hm_zone_sup;
			hm_zone_sup = hm [zone];
		}
		// aberration interpolee
		float aber_inter = aber_zn_inf
                   + (aber_zn_sup - aber_zn_inf)
                   * (h - hm_zone_inf) / (hm_zone_sup - hm_zone_inf);
		
		// integ les pentes deduites des aber, donne defaut sur l'onde
		onde [i+1] = onde [i] + aber_inter * resol_h / foc;
	}
	/*
La coupe de front d'onde qui vient d'etre calculee dans le tableau
onde [] est celle correspondant au foyer de moindre aberation.
C'est celle qui a les ecarts de pente les plus faibles.
Ce n'est pas forcement celle qu a les plus petits defauts "pic-vallee"
c'est a dire le paraboloide de reference en contact avec le front d'onde
a la fois au bord et au centre du miroir.
Pour calculer les defauts "Pic-Vallee" du front d'onde il faut
prendre le paraboloide ci-dessus comme surface de reference.
	*/
	// soustraction de la parabole de ref pour foyer "pic-Vallee"
	float ref_au_bord = onde[nb_pt_h] /nb_pt_h/nb_pt_h;
	for (i = 0 ; i<= nb_pt_h ; i++)
		onde_pv [i] = onde [i] - i*i* ref_au_bord;
		
	// calcul du "lambda sur x pic-vallee"
	onde_min = onde_max = 0;
	for (i = 0 ; i<= nb_pt_h ; i++)
	{
		if      (onde_pv [i] > onde_max) onde_max = onde_pv [i];
		else if (onde_pv [i] < onde_min) onde_min = onde_pv [i];
	}
	printf ("Defaut max :   %4.0f nanometres.\
        Miroir a lambda sur %4.1f Pic-Vallee\n",
		(onde_max - onde_min)*1e9, lambda/(onde_max - onde_min) );
		
	// calcul du "lambda sur x   R M S"
	float defauts_rms = 0;
	for (i = 0 ; i<= nb_pt_h ; i++)
		defauts_rms += onde_pv[i] * onde_pv[i];
	defauts_rms /= nb_pt_h;  // moyenne des carres des defauts
	defauts_rms = sqrt (defauts_rms);  // Root Mean Square
	printf ("Defaut R M S : %4.0f nanometres.\
        Miroir a lambda sur %4.1f R M S\n",
		defauts_rms*1e9, lambda/defauts_rms);
	}

// ------------ enregistrement sur disque --------------
/*
Enregistre sur disque (en clair) un tableau de deux colones contenant
1ere colone: coupe du front d'onde pour le foyer de moindre aberration
2eme colone: coupe du front d'onde pour un foyer ajuste donnant les plus 
petits defauts "pic-vallee".
Le nom du fichier disque est defini dans le code ci-dessous.
*/
void miroir::enreg_onde (void)
{
	int i;		
	FILE	*fp = fopen ("coupe.txt","w");
	
	fprintf (fp, "abscis \tcoupe MA\tcoupe PV \n");
	for (i = 0 ; i<= nb_pt_h ; i++)
		fprintf (fp, "%5.3e\t%5.3e\t%5.3e\n", i*resol_h,
				onde [i], onde_pv [i]);
	fclose (fp);
}

// ------------ affichage rustique sur console --------------
/*
Dessine avec des lettres sur la fenetre console la forme du front d'onde
L'amplitude est normalisee verticalement a 15 lignes et le rayon miroir
a 76 caracteres
*/
void miroir::imprime_rustic_onde (void)
{
	int i, j;
	int ligne;                   // hauteur du front d'onde en lignes
	const int max_lign = 15;     // normalisation onde en lignes impr
	char *onde_impr[max_lign +1];// dessin onde imprim
	
	for (i = 0 ; i<= max_lign ; i++)  // alloc tableau 2D pour dessin
		onde_impr[i] = new char [nb_pt_h +1];
	
	for (i = 0 ; i<= max_lign ; i++)
	for (j = 0 ; j<= nb_pt_h ; j++) onde_impr[i][j] = ' ';

	// dessin de l'axe horizontal dans le tableau
	int ligne_axe = -onde_min * max_lign / (onde_max - onde_min);
	if (ligne_axe <0) ligne_axe = 0;                // securite
	if (ligne_axe > max_lign) ligne_axe = max_lign; // securite
	for (j = 0 ; j<= nb_pt_h ; j++) onde_impr [ligne_axe][j] = '_';

	// dessin de l'onde dans le tableau
	for (j = 0 ; j<= nb_pt_h ; j++)
	{
		ligne = (onde_pv[j] - onde_min)
		              * max_lign / (onde_max - onde_min);
		if (ligne <0) ligne = 0;                // securite
		if (ligne > max_lign) ligne = max_lign; // securite
		onde_impr [ligne][j] = '.';
	}

	// affichage du tableau
	for (i = max_lign ; i>=0 ; i--)
	{
		printf ("\n");
		for (j = 0 ; j<= nb_pt_h ; j++) 
			printf ("%c", onde_impr[i][j]);
	}
	printf ("\n\nCoupe de l'onde reflechie par le miroir ");
	printf ("du centre (a gauche) au bord (a droite)\n");
	printf ("%d lignes en hauteur correspondent a %3.0f nanometres\n",
		max_lign, (onde_max - onde_min)* 1e9);
}