/****************************************************************************/
/****    binthedata.c                                                    ****/
/****                                                                    ****/
/****    Lecture des donnees mars et "binnage" au format mgsab           ****/
/****    format (temps), lat, lon, alt, Br, Btheta, Bphi                 ****/
/****************************************************************************/

/***********************/
/****    includes    ***/
/***********************/

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

/**********************/
/****    Define    ****/
/**********************/

#define M_PI   3.14159265358979323846
#define PI_180 (M_PI / 180.0)
#define RM 3393.5
#define AMP (int)(alt_max-alt_min)
/***********************/
/****    Formats    ****/
/***********************/

#define READ_bin "%lf %lf %lf %lf %lf %lf %lf %d %lf %lf %lf\n" 
#define READs "%lf %lf %lf %lf %lf %lf\n"
#define READl "%lf %lf %lf %lf %lf %lf %lf\n"
#define READld "%lf %lf %lf %lf %lf %lf %lf %1s\n"
#define LECTl1 &dday, &lon, &lat, &alt, &br, &bt, &bp
#define LECTl2 &dday, &lat, &lon, &alt, &br, &bt, &bp
#define LECTld &dday, &lat, &lon, &alt, &br, &bt, &bp, &dl
#define READb "%lf %lf %lf %lf %lf %lf %d %lf %lf %lf\n" 
#define LECTb &lon, &lat, &alt, &br, &bt, &bp, &d, &dd, &dd, &dd

#define WRITE_bin "%13.9lf %7.3lf %7.3lf %7.3lf %12.5lf %12.5lf %12.5lf %3d %10.3lf %10.3lf %10.3lf\n"

#define WRITEs "%13.9lf %8.3lf %8.3lf %8.3lf %10.3lf %10.3lf %10.3lf\n"
#define WRITEl "%13.9lf %8.3lf %8.3lf %8.3lf %10.3lf %10.3lf %10.3lf\n"
#define ECRITl dday, lat, lon, alt, br, bt, bp
#define WRITEld "%13.9lf %8.3lf %8.3lf %8.3lf %10.3lf %10.3lf %10.3lf %1s\n"
#define ECRITld dday, lat, lon, alt, br, bt, bp, dl

#define NMAX 50

/****************************/
/****    Main Program    ****/
/****************************/

int main(argc, argv)
int argc;
char *argv[];
      {
	FILE   *fin, *fout;
	FILE   *fout_00, *fout_01, *fout_02, *fout_03, *fout_04;
	FILE   *fout_05, *fout_06, *fout_07, *fout_08, *fout_09;
	FILE   *fout_10, *fout_11, *fout_12, *fout_13, *fout_14;
	FILE   *fout_15, *fout_16, *fout_17, *fout_18, *fout_19;
	FILE   *fout_20, *fout_21, *fout_22, *fout_23, *fout_24;
	FILE   *fout_25, *fout_26, *fout_27, *fout_28, *fout_29;
	FILE   *fout_30, *fout_31, *fout_32, *fout_33, *fout_34;
	FILE   *fout_35, *fout_36, *fout_37, *fout_38, *fout_39;
	FILE   *fout_40, *fout_41, *fout_42, *fout_43, *fout_44;
	FILE   *fout_45, *fout_46, *fout_47, *fout_48, *fout_49;
	char   str[200];
	double alt_max, alt_min;
	double lat_max, lat_min;
	double lon_max, lon_min;
	int    choix;
	int    data_in, data_out;
	int    filled_bin=0;
	int    filled_bin1=0;
	int    filled_bin2=0;
	double inc_lon, inc_lat, inc_alt;
	double amp_lon, amp_lat, amp_alt;
	int    npt_lon, npt_lat, npt_alt;
	int    l, L, A;
	double dday, lat, lon, alt, br, bt, bp;
	int    index_alt;
	int    npt_tot;
	double **BR, **BT, **BP;
	double **alt_var;
	double *BR_max, *BT_max, *BP_max;
	double *BR_min, *BT_min, *BP_min;
	double *BR_moy, *BT_moy, *BP_moy;
	double *BR_med, *BT_med, *BP_med;
	double *BR_std, *BT_std, *BP_std;
	double std_br, std_bt, std_bp;
	double sum_br, sum_bt, sum_bp;
	double dif_br, dif_bt, dif_bp;
	int    i, j, k, index;
	int    *NPT;
	char dl[1];
	double temp;
	int d;
	double dd;
	dday = 0.0;

	/* Introduction */
    
	printf("############################################################\n");
	printf("binthedata.exe\n\n");

	/* Test the number of parameters */
	
	if (argc != 6 && argc != 10 && argc != 13)
	  {
	    printf("\tYou must call 6, 10 or 13 parameters!\n");
	    printf("\tbinthedata.exe <choix> <fin> <fout> <alt_min> <alt_max>\n");
	    printf("\tand eventually <lat_min> <lat_max> <lon_min> <<lon_max>\n");
	    printf("\tand eventually <inc_lat> <inc_lon> <inc_alt>\n");
	    printf("\t<choix> = 1 (median) or 2 (average) or 3 (minimum) or 4 (maximum)\n");
	    printf("\t<choix> = 5 : one data per orbit per bin.\n");
	    exit(-1);
	  }
	
	choix = atoi(argv[1]);
	if (choix != 1 && choix != 2 && choix != 3 && choix != 4 && choix != 5)
	  {
	    printf("\tBad <choix> value. Must be 1 <-> 5. Please check!\n");
	    exit(-1);
	  }

	fin = fopen(argv[2], "r");
	if (fin == NULL)
	  {
	    printf("Reading error on input file <%s>\n. Please check!\n", argv[2]);
	    exit(-1);
	  }

	fout = fopen(argv[3], "w");
	if (fin == NULL)
	  {
	    printf("Writing error on output file <%s>\n. Please check!\n", argv[3]);
	    fclose(fin);
	    exit(-1);
	  }

	/* Minimum and Maximum altitude */
	
	alt_min = atof(argv[4]);
	alt_max = atof(argv[5]);
	
	if (argc == 10 || argc == 13)
	  {
	    lat_min = atof(argv[6]);
	    lat_max = atof(argv[7]);
	    lon_min = atof(argv[8]);
	    lon_max = atof(argv[9]);
	  }
	else
	  {
	    lat_min = -90.;
	    lat_max = 90.;
	    lon_min = 0.;
	    lon_max = 360.;
	  }
	
	if (argc == 13)
	  {
	    inc_lat = atof(argv[10]);
	    inc_lon = atof(argv[11]);
	    inc_alt = atof(argv[12]);
	  }
	else
	  {
	    inc_lat = 1.;
	    inc_lon = 1.;
	    inc_alt = 10.;
	  }
	
	/****    Verification    ****/
	
	if (alt_min < 0.)
	  {
	    do
	      {
		printf("\t<alt_min> lower than 0! Please enter a new value (0.)\n");
		fgets(str, 150, stdin);
		if (str[0] == '\n')	alt_min = 0.;
		else 
		  {
		    sscanf(str, "%lf", &alt_min);
		  }
	      }
	    while (alt_max < 0.);
	  }
	if (alt_max > 1000.)
	  {
	    do
	      {
		printf("\t<alt_max> greater than 1000! Please enter a new value (1000.)\n");
		fgets(str, 150, stdin);
		if (str[0] == '\n')	alt_max = 1000.;
		else 
		  {
		    sscanf(str, "%lf", &alt_max);
		  }
	      }
	    while (alt_max > 1000.);
	  }
	
	if (alt_min > alt_max)
	  {
	    printf("\t<alt_min> greater than <alt_max>.\n");
	    printf("\tPlease check. Exit!\n");
	    fclose(fin);
	    exit(-1);
	  }
	
	if (lat_min > lat_max)
	  {
	    printf("\tError in latitude range. Please check!\n");
	    exit(-1);
	  }

	/* parametres */

	/* Initialisation */
	
	amp_lat = lat_max - lat_min;
	if (lon_max < lon_min)
	    amp_lon = 360. + lon_max - lon_min;
	else
	    amp_lon = lon_max - lon_min;
	amp_alt = alt_max - alt_min;

	npt_lat = (int)(amp_lat / (double)inc_lat);
	npt_lon = (int)(amp_lon / (double)inc_lon);
	npt_alt = (int)(amp_alt / (double)inc_alt);

	fprintf(stdout, "Nombre de points = %d (lat) x %d (lon) x %d (alt)\n",
		npt_lat, npt_lon, npt_alt);

	l = 0      ; L = 0      ; A = 0      ;
	lat = lat_min + (double)l*inc_lat;
	fprintf(stdout, "\tL'intervalle lat est %7.3lf/", lat);
	l = npt_lat; L = 0      ; A = 0      ;
	lat = lat_min + (double)l*inc_lat;
	fprintf(stdout, "%7.3lf\n", lat); 
	l = 0      ; L = 0      ; A = 0      ;
	lon = lon_min + (double)L*inc_lon;
	fprintf(stdout, "\tL'intervalle lon est %7.3lf/", lon);
	l = 0      ; L = npt_lon; A = 0      ;
	lon = lon_min + (double)L*inc_lon;
	fprintf(stdout, "%7.3lf\n", lon); 
	l = 0      ; L = 0      ; A = 0      ;
	alt = alt_min + (double)A*inc_alt;
	fprintf(stdout, "\tL'intervalle alt est %7.3lf/", alt);
	l = 0      ; L = 0      ; A = npt_alt;
	alt = alt_min + (double)A*inc_alt;
	fprintf(stdout, "%7.3lf\n", alt); 
	fflush(stdout);

	/****    Lecture et ecriture des donnees triees par altitude    ****/

	fprintf(stdout, "\n\tOuverture des fichiers"); 
	fflush(stdout);
	fout_00 = fopen("temp_00.alt", "w");
	fout_01 = fopen("temp_01.alt", "w");
	fout_02 = fopen("temp_02.alt", "w");
	fout_03 = fopen("temp_03.alt", "w");
	fout_04 = fopen("temp_04.alt", "w");
	fout_05 = fopen("temp_05.alt", "w");
	fout_06 = fopen("temp_06.alt", "w");
	fout_07 = fopen("temp_07.alt", "w");
	fout_08 = fopen("temp_08.alt", "w");
	fout_09 = fopen("temp_09.alt", "w");
	fout_10 = fopen("temp_10.alt", "w");
	fout_11 = fopen("temp_11.alt", "w");
	fout_12 = fopen("temp_12.alt", "w");
	fout_13 = fopen("temp_13.alt", "w");
	fout_14 = fopen("temp_14.alt", "w");
	fout_15 = fopen("temp_15.alt", "w");
	fout_16 = fopen("temp_16.alt", "w");
	fout_17 = fopen("temp_17.alt", "w");
	fout_18 = fopen("temp_18.alt", "w");
	fout_19 = fopen("temp_19.alt", "w");
	fout_20 = fopen("temp_20.alt", "w");
	fout_21 = fopen("temp_21.alt", "w");
	fout_22 = fopen("temp_22.alt", "w");
	fout_23 = fopen("temp_23.alt", "w");
	fout_24 = fopen("temp_24.alt", "w");
	fout_25 = fopen("temp_25.alt", "w");
	fout_26 = fopen("temp_26.alt", "w");
	fout_27 = fopen("temp_27.alt", "w");
	fout_28 = fopen("temp_28.alt", "w");
	fout_29 = fopen("temp_29.alt", "w");
	fout_30 = fopen("temp_30.alt", "w");
	fout_31 = fopen("temp_31.alt", "w");
	fout_32 = fopen("temp_32.alt", "w");
	fout_33 = fopen("temp_33.alt", "w");
	fout_34 = fopen("temp_34.alt", "w");
	fout_35 = fopen("temp_35.alt", "w");
	fout_36 = fopen("temp_36.alt", "w");
	fout_37 = fopen("temp_37.alt", "w");
	fout_38 = fopen("temp_38.alt", "w");
	fout_39 = fopen("temp_39.alt", "w");
	fout_40 = fopen("temp_40.alt", "w");
	fout_41 = fopen("temp_41.alt", "w");
	fout_42 = fopen("temp_42.alt", "w");
	fout_43 = fopen("temp_43.alt", "w");
	fout_44 = fopen("temp_44.alt", "w");
	fout_45 = fopen("temp_45.alt", "w");
	fout_46 = fopen("temp_46.alt", "w");
	fout_47 = fopen("temp_47.alt", "w");
	fout_48 = fopen("temp_48.alt", "w");
	fout_49 = fopen("temp_49.alt", "w");
		
	data_in = data_out = 0;

	fprintf(stdout,"\tLecture et tri par altitude...\n");
	fflush(stdout);
	while(!feof(fin))
	  {
/* 	    fscanf(fin, READl, LECTl1); */
/* 	    fscanf(fin, READl, LECTl2); */
	    fscanf(fin, READld, LECTld);
	    data_in ++;
	    if (lon < 0.)lon += 360.;
	    else if (lon >= 360.)lon -= 360.;
	    if (lat >= lat_min && lat < lat_max && alt >= alt_min && alt < alt_max)
	      {
		if ((lon_min < lon_max && lon >= lon_min && lon < lon_max) ||
		    (lon_min > lon_max && (lon >= lon_min || lon < lon_max)))
		  {
		    index_alt = (int)((alt - alt_min)/inc_alt);
		    data_out ++;
		    switch(index_alt)
		      {
		      case 0:
			{
			  fprintf(fout_00, WRITEl, ECRITl);
			  break;
			}
		      case 1:
			{
			  fprintf(fout_01, WRITEl, ECRITl);
			  break;
			}
		      case 2:
			{
			  fprintf(fout_02, WRITEl, ECRITl);
			  break;
			}
		      case 3:
			{
			  fprintf(fout_03, WRITEl, ECRITl);
			  break;
			}
		      case 4:
			{
			  fprintf(fout_04, WRITEl, ECRITl);
			  break;
			}
		      case 5:
			{
			  fprintf(fout_05, WRITEl, ECRITl);
			  break;
			}
		      case 6:
			{
			  fprintf(fout_06, WRITEl, ECRITl);
			  break;
			}
		      case 7:
			{
			  fprintf(fout_07, WRITEl, ECRITl);
			  break;
			}
		      case 8:
			{
			  fprintf(fout_08, WRITEl, ECRITl);
			  break;
			}
		      case 9:
			{
			  fprintf(fout_09, WRITEl, ECRITl);
			  break;
			}
		      case 10:
			{
			  fprintf(fout_10, WRITEl, ECRITl);
			  break;
			}
		      case 11:
			{
			  fprintf(fout_11, WRITEl, ECRITl);
			  break;
			}
		      case 12:
			{
			  fprintf(fout_12, WRITEl, ECRITl);
			  break;
			}
		      case 13:
			{
			  fprintf(fout_13, WRITEl, ECRITl);
			  break;
			}
		      case 14:
			{
			  fprintf(fout_14, WRITEl, ECRITl);
			  break;
			}
		      case 15:
			{
			  fprintf(fout_15, WRITEl, ECRITl);
			  break;
			}
		      case 16:
			{
			  fprintf(fout_16, WRITEl, ECRITl);
			  break;
			}
		      case 17:
			{
			  fprintf(fout_17, WRITEl, ECRITl);
			  break;
			}
		      case 18:
			{
			  fprintf(fout_18, WRITEl, ECRITl);
			  break;
			}
		      case 19:
			{
			  fprintf(fout_19, WRITEl, ECRITl);
			  break;
			}
		      case 20:
			{
			  fprintf(fout_20, WRITEl, ECRITl);
			  break;
			}
		      case 21:
			{
			  fprintf(fout_21, WRITEl, ECRITl);
			  break;
			}
		      case 22:
			{
			  fprintf(fout_22, WRITEl, ECRITl);
			  break;
			}
		      case 23:
			{
			  fprintf(fout_23, WRITEl, ECRITl);
			  break;
			}
		      case 24:
			{
			  fprintf(fout_24, WRITEl, ECRITl);
			  break;
			}
		      case 25:
			{
			  fprintf(fout_25, WRITEl, ECRITl);
			  break;
			}
		      case 26:
			{
			  fprintf(fout_26, WRITEl, ECRITl);
			  break;
			}
		      case 27:
			{
			  fprintf(fout_27, WRITEl, ECRITl);
			  break;
			}
		      case 28:
			{
			  fprintf(fout_28, WRITEl, ECRITl);
			  break;
			}
		      case 29:
			{
			  fprintf(fout_29, WRITEl, ECRITl);
			  break;
			}
		      case 30:
			{
			  fprintf(fout_30, WRITEl, ECRITl);
			  break;
			}
		      case 31:
			{
			  fprintf(fout_31, WRITEl, ECRITl);
			  break;
			}
		      case 32:
			{
			  fprintf(fout_32, WRITEl, ECRITl);
			  break;
			}
		      case 33:
			{
			  fprintf(fout_33, WRITEl, ECRITl);
			  break;
			}
		      case 34:
			{
			  fprintf(fout_34, WRITEl, ECRITl);
			  break;
			}
		      case 35:
			{
			  fprintf(fout_35, WRITEl, ECRITl);
			  break;
			}
		      case 36:
			{
			  fprintf(fout_36, WRITEl, ECRITl);
			  break;
			}
		      case 37:
			{
			  fprintf(fout_37, WRITEl, ECRITl);
			  break;
			}
		      case 38:
			{
			  fprintf(fout_38, WRITEl, ECRITl);
			  break;
			}
		      case 39:
			{
			  fprintf(fout_39, WRITEl, ECRITl);
			  break;
			}
		      case 40:
			{
			  fprintf(fout_40, WRITEl, ECRITl);
			  break;
			}
		      case 41:
			{
			  fprintf(fout_41, WRITEl, ECRITl);
			  break;
			}
		      case 42:
			{
			  fprintf(fout_42, WRITEl, ECRITl);
			  break;
			}
		      case 43:
			{
			  fprintf(fout_43, WRITEl, ECRITl);
			  break;
			}
		      case 44:
			{
			  fprintf(fout_44, WRITEl, ECRITl);
			  break;
			}
		      case 45:
			{
			  fprintf(fout_45, WRITEl, ECRITl);
			  break;
			}
		      case 46:
			{
			  fprintf(fout_46, WRITEl, ECRITl);
			  break;
			}
		      case 47:
			{
			  fprintf(fout_47, WRITEl, ECRITl);
			  break;
			}
		      case 48:
			{
			  fprintf(fout_48, WRITEl, ECRITl);
			  break;
			}
		      case 49:
			{
			  fprintf(fout_49, WRITEl, ECRITl);
			  break;
			}
		      }
		  }
	      }
	  }
	
	fclose(fout_00);
	fclose(fout_01);
	fclose(fout_02);
	fclose(fout_03);
	fclose(fout_04);
	fclose(fout_05);
	fclose(fout_06);
	fclose(fout_07);
	fclose(fout_08);
	fclose(fout_09);
	fclose(fout_10);
	fclose(fout_11);
	fclose(fout_12);
	fclose(fout_13);
	fclose(fout_14);
	fclose(fout_15);
	fclose(fout_16);
	fclose(fout_17);
	fclose(fout_18);
	fclose(fout_19);
	fclose(fout_20);
	fclose(fout_21);
	fclose(fout_22);
	fclose(fout_23);
	fclose(fout_24);
	fclose(fout_25);
	fclose(fout_26);
	fclose(fout_27);
	fclose(fout_28);
	fclose(fout_29);
	fclose(fout_30);
	fclose(fout_31);
	fclose(fout_32);
	fclose(fout_33);
	fclose(fout_34);
	fclose(fout_35);
	fclose(fout_36);
	fclose(fout_37);
	fclose(fout_38);
	fclose(fout_39);
	fclose(fout_40);
	fclose(fout_41);
	fclose(fout_42);
	fclose(fout_43);
	fclose(fout_44);
	fclose(fout_45);
	fclose(fout_46);
	fclose(fout_47);
	fclose(fout_48);
	fclose(fout_49);

	fclose(fin);

	/****    Lecture fichier par fichier    ****/

	fprintf(stdout,"\n\tLecture OK. Travail par bande d'altitude\n");
	fflush(stdout);
	
	npt_tot = npt_lat * npt_lon;

	BR = (double **)malloc(npt_tot*sizeof(double*));
	BT = (double **)malloc(npt_tot*sizeof(double*));
	BP = (double **)malloc(npt_tot*sizeof(double*));
	alt_var = (double **)malloc(npt_tot*sizeof(double*));
	for (i = 0; i < npt_tot; i++)
	  {
	    BR[i] = (double *)malloc(NMAX*sizeof(double));
	    BT[i] = (double *)malloc(NMAX*sizeof(double));
	    BP[i] = (double *)malloc(NMAX*sizeof(double));
	    alt_var[i] = (double *)malloc(NMAX*sizeof(double));
	  }

	if (choix == 1)
	  {
	    BR_med = (double *)malloc(npt_tot*sizeof(double));
	    BT_med = (double *)malloc(npt_tot*sizeof(double));
	    BP_med = (double *)malloc(npt_tot*sizeof(double));
	  }
	else if (choix == 2)
	  {
	    BR_moy = (double *)malloc(npt_tot*sizeof(double));
	    BT_moy = (double *)malloc(npt_tot*sizeof(double));
	    BP_moy = (double *)malloc(npt_tot*sizeof(double));
	  }
	else if (choix == 3)
	  {
	    BR_min = (double *)malloc(npt_tot*sizeof(double));
	    BT_min = (double *)malloc(npt_tot*sizeof(double));
	    BP_min = (double *)malloc(npt_tot*sizeof(double));
	  }
	else if (choix == 4)
	  {
	    BR_max = (double *)malloc(npt_tot*sizeof(double));
	    BT_max = (double *)malloc(npt_tot*sizeof(double));
	    BP_max = (double *)malloc(npt_tot*sizeof(double));
	  }
	BR_std = (double *)malloc(npt_tot*sizeof(double));
	BT_std = (double *)malloc(npt_tot*sizeof(double));
	BP_std = (double *)malloc(npt_tot*sizeof(double));
	
	NPT = (int *)malloc(npt_tot*sizeof(int*));

	BR[0][0] = 0.;
	for (A = 0; A < 50 ; A ++)
	  {
	    fprintf(stdout, "\t\tAltitude = %lf\n", (double)(A * inc_alt) + alt_min);
	    sprintf(str, "temp_%02d.alt", A);
	    fin = fopen(str, "r");

	    for (i = 0 ; i < npt_tot ; i ++)
	      {
		NPT[i] = 0;
	      }
	    while (!feof(fin))
	      {
		d = fscanf(fin, READl, LECTl2);
		if (d == 7)
		  {
		    l = (lat - lat_min) / inc_lat;
		    if (lon_max > lon_min)
		      {
			L = (int)((lon - lon_min ) / inc_lon);
		      }
		    else
		      {
			if (lon >= lon_min)
			  {
			    L = (int)((lon - lon_min ) /inc_lon);
			  }
			else
			  {
			    L = (int)((lon + 360. - lon_min ) /inc_lon);
			  }
		      }


		    index = l * npt_lon + L;
		    for (i = 0 ; i < NPT[index]; i++)
		      {
			if (br < BR[index][i])
			  {
			    temp = BR[index][i];
/* 			    j = i; */
/* 			    fprintf(stdout," %d\n",j);fflush(stdout); */
			    BR[index][i] = br;
			    br = temp;
			    for (j = i+1 ; j < NPT[index]; j++)
			      {
				temp = BR[index][j];
				BR[index][j] = br;
				br = temp;
			      }
			    break;
			  }
		      }
		    BR[index][NPT[index]] = br;
		    
		    for (i = 0 ; i < NPT[index]; i++)
		      {
			if (bt < BT[index][i])
			  {
			    temp = BT[index][i];
			    BT[index][i] = bt;
			    bt = temp;
			    for (j = i+1 ; j < NPT[index]; j++)
			      {
				temp = BT[index][j];
				BT[index][j] = bt;
				bt = temp;
			      }
			    break;
			  }
		      }
		    BT[index][NPT[index]] = bt;
		    
		    for (i = 0 ; i < NPT[index]; i++)
		      {
			if (bp < BP[index][i])
			  {
			    temp = BP[index][i];
			    BP[index][i] = bp;
			    bp = temp;
			    for (j = i+1 ; j < NPT[index]; j++)
			      {
				temp = BP[index][j];
				BP[index][j] = bp;
				bp = temp;
			      }
			    break;
			  }
		      }
		    BP[index][NPT[index]] = bp;

		    alt_var[index][NPT[index]] = alt;

		    NPT[index] += 1;

		    
/* 		fprintf(stdout, "Case %d, %d points lus\n", index, NPT[index]);fflush(stdout); */
		  }
	      }
	    fclose(fin);

	    if (choix == 1)
	      {
		for (l = 0 ; l < npt_lat; l++)
		  {
		    for (L = 0 ; L < npt_lon ; L++)
		      {
			index = l * npt_lon + L;
			if (NPT[index] > 0)
			  {
/* 			    for (k = 0 ; k < NPT[index]; k++) */
/* 			      { */
/* 				fprintf(stdout, "%3d %8.3lf %8.3lf %8.3lf\n", k,BR[index][k],BT[index][k],BP[index][k]);fflush(stdout); */
/* 			      } */
			    if (NPT[index] == 1)
			      {
				BR_med[index]=BR[index][NPT[index]-1];
				BT_med[index]=BT[index][NPT[index]-1];
				BP_med[index]=BP[index][NPT[index]-1];
			      }
			    else if ( NPT[index]%2 == 0)
			      {
				BR_med[index]=(BR[index][NPT[index]/2-1]+BR[index][NPT[index]/2])/2.;
				BT_med[index]=(BT[index][NPT[index]/2-1]+BT[index][NPT[index]/2])/2.;
				BP_med[index]=(BP[index][NPT[index]/2-1]+BP[index][NPT[index]/2])/2.;
			      }
			    else if ( NPT[index]%2 != 0)
			      {
				BR_med[index]=BR[index][(NPT[index]-1)/2];
				BT_med[index]=BT[index][(NPT[index]-1)/2];
				BP_med[index]=BP[index][(NPT[index]-1)/2];
			      }

			    std_br = std_bt = std_bp = 0;

			    for (k = 0; k < NPT[index] ; k ++)
			      {
				dif_br = (BR[index][k] - BR_med[index]);
				std_br  += dif_br*dif_br;
				dif_bt = (BT[index][k] - BT_med[index]);
				std_bt  += dif_bt*dif_bt;
				dif_bp = (BP[index][k] - BP_med[index]);
				std_bp  += dif_bp*dif_bp;
			      }
			    if (NPT[index] > 2)
			      {
				BR_std[index] = sqrt(std_br/(NPT[index]-1));
				BT_std[index] = sqrt(std_bt/(NPT[index]-1));
				BP_std[index] = sqrt(std_bp/(NPT[index]-1));
			      }
			    else 
			      {
				BR_std[index] = 3.;
				BT_std[index] = 3.;
				BP_std[index] = 3.;
			      }
			    lat = lat_min + (double)l*inc_lat + inc_lat/2.;
			    lon = lon_min + (double)L*inc_lon + inc_lon/2.;
			    alt = alt_min + (double)A*inc_alt + inc_alt/2.;
			    
			    fprintf(fout, WRITE_bin, dday, lon, lat, alt, 
				    BR_med[index], BT_med[index], BP_med[index],
				    NPT[index],
				    BR_std[index], BT_std[index], BP_std[index]);

			  }
		      }
		  }
	      }

	    else if (choix == 2)
	      {
		for (l = 0 ; l < npt_lat; l++)
		  {
		    for (L = 0 ; L < npt_lon ; L++)
		      {
			index = l * npt_lon + L;
			if (NPT[index] > 0)
			  {
			    sum_br = sum_bt = sum_bp = 0.;
			    for (k = 0; k < NPT[index] ; k ++)
			      {
				sum_br += BR[index][k];
				sum_bt += BT[index][k];
				sum_bp += BP[index][k];
			      }
			    BR_moy[index] = sum_br / (double)NPT[index];
			    BT_moy[index] = sum_bt / (double)NPT[index];
			    BP_moy[index] = sum_bp / (double)NPT[index];

			    std_br = std_bt = std_bp = 0;

			    for (k = 0; k < NPT[index] ; k ++)
			      {
				dif_br = (BR[index][i] - BR_moy[index]);
				std_br  += dif_br*dif_br;
				dif_bt = (BT[index][i] - BT_moy[index]);
				std_bt  += dif_bt*dif_bt;
				dif_bp = (BP[index][i] - BP_moy[index]);
				std_bp  += dif_bp*dif_bp;
			      }
			    if (NPT[index] > 2)
			      {
				BR_std[index] = sqrt(std_br/(NPT[index]-1));
				BT_std[index] = sqrt(std_bt/(NPT[index]-1));
				BP_std[index] = sqrt(std_bp/(NPT[index]-1));
			      }
			    else 
			      {
				BR_std[index] = 3.;
				BT_std[index] = 3.;
				BP_std[index] = 3.;
			      }
			    lat = lat_min + (double)l*inc_lat + inc_lat/2.;
			    lon = lon_min + (double)L*inc_lon + inc_lon/2.;
			    if (lon >= 360.)
			      lon -= 360.;
			    alt = alt_min + (double)A*inc_alt + inc_alt/2.;

			    fprintf(fout, WRITE_bin, dday, lon, lat, alt, 
				    BR_moy[index], BT_moy[index], BP_moy[index],
				    NPT[index],
				    BR_std[index], BT_std[index], BP_std[index]);
			  }
		      }
		  }
	      }


	    else if (choix == 3)
	      {
		for (l = 0 ; l < npt_lat; l++)
		  {
		    for (L = 0 ; L < npt_lon ; L++)
		      {
			index = l * npt_lon + L;
			if (NPT[index] > 0)
			  {
			    BR_min[index] = BR[index][0];
			    BT_min[index] = BT[index][0];
			    BP_min[index] = BP[index][0];

			    std_br = std_bt = std_bp = 0;

			    for (k = 0; k < NPT[index] ; k ++)
			      {
				dif_br = (BR[index][i] - BR_min[index]);
				std_br  += dif_br*dif_br;
				dif_bt = (BT[index][i] - BT_min[index]);
				std_bt  += dif_bt*dif_bt;
				dif_bp = (BP[index][i] - BP_min[index]);
				std_bp  += dif_bp*dif_bp;
			      }
			    if (NPT[index] > 2)
			      {
				BR_std[index] = sqrt(std_br/(NPT[index]-1));
				BT_std[index] = sqrt(std_bt/(NPT[index]-1));
				BP_std[index] = sqrt(std_bp/(NPT[index]-1));
			      }
			    else 
			      {
				BR_std[index] = 3.;
				BT_std[index] = 3.;
				BP_std[index] = 3.;
			      }
			    lat = lat_min + (double)l*inc_lat + inc_lat/2.;
			    lon = lon_min + (double)L*inc_lon + inc_lon/2.;
			    if (lon >= 360.)
			      lon -= 360.;
			    alt = alt_min + (double)A*inc_alt + inc_alt/2.;

			    fprintf(fout, WRITE_bin, dday, lon, lat, alt, 
				    BR_min[index], BT_min[index], BP_min[index],
				    NPT[index],
				    BR_std[index], BT_std[index], BP_std[index]);
			  }
		      }
		  }
	      }

	    else if (choix == 4)
	      {
		for (l = 0 ; l < npt_lat; l++)
		  {
		    for (L = 0 ; L < npt_lon ; L++)
		      {
			index = l * npt_lon + L;
			if (NPT[index] > 0)
			  {
			    BR_max[index] = BR[index][NPT[index]-1];
			    BT_max[index] = BT[index][NPT[index]-1];
			    BP_max[index] = BP[index][NPT[index]-1];

			    std_br = std_bt = std_bp = 0;

			    for (k = 0; k < NPT[index] ; k ++)
			      {
				dif_br = (BR[index][i] - BR_max[index]);
				std_br  += dif_br*dif_br;
				dif_bt = (BT[index][i] - BT_max[index]);
				std_bt  += dif_bt*dif_bt;
				dif_bp = (BP[index][i] - BP_max[index]);
				std_bp  += dif_bp*dif_bp;
			      }
			    if (NPT[index] > 2)
			      {
				BR_std[index] = sqrt(std_br/(NPT[index]-1));
				BT_std[index] = sqrt(std_bt/(NPT[index]-1));
				BP_std[index] = sqrt(std_bp/(NPT[index]-1));
			      }
			    else 
			      {
				BR_std[index] = 3.;
				BT_std[index] = 3.;
				BP_std[index] = 3.;
			      }
			    lat = lat_min + (double)l*inc_lat + inc_lat/2.;
			    lon = lon_min + (double)L*inc_lon + inc_lon/2.;
			    if (lon >= 360.)
			      lon -= 360.;
			    alt = alt_min + (double)A*inc_alt + inc_alt/2.;

			    fprintf(fout, WRITE_bin, dday, lon, lat, alt, 
				    BR_max[index], BT_max[index], BP_max[index],
				    NPT[index],
				    BR_std[index], BT_std[index], BP_std[index]);
			  }
		      }
		  }
	      }


	    for (l = 0; l < npt_lat; l++)
	      {
		for (L = 0 ; L < npt_lon ; L++)
		  {
		    index = l * npt_lon + L;
		    if (NPT[index] > 0)
		      {
			filled_bin ++;
			if (NPT[index] > 2)
			  {
			    filled_bin2 ++;
			  }
			if (NPT[index] > 1)
			  {
			    filled_bin1 ++;
			  }
		      }
		  }
	      }

	    for (l = 0; l < npt_lat; l++)
	      {
		for (L = 0 ; L < npt_lon ; L++)
		  {
		    index = l * npt_lon + L;
		    if (NPT[index] > 1)
		      {
			lat = lat_min + (double)l*inc_lat + inc_lat/2.;
			lon = lon_min + (double)L*inc_lon + inc_lon/2.;
			if (lon >= 360.)
			  lon -= 360.;
			fprintf(fout, "%7.3lf %7.3lf", lon, lat);
			for (i=0;i<NPT[index]; i++)
			  {
			    fprintf(fout, " %7.3lf", alt_var[index][i]);
			  }
			fprintf(fout, "\n");
		      }
		  }
	      }
	  }

	fclose(fout);

	fprintf(stdout, "\tNombre de donnees lues      : %10d\n", data_in);
	fprintf(stdout, "\tNombre de donnees utilisees : %10d\n", data_out);
	fprintf(stdout, "\tNombre de bin remplis       : %10d\n", filled_bin);
	fprintf(stdout, "\tNombre de bin remplis (N>1) : %10d\n", filled_bin1);
	fprintf(stdout, "\tNombre de bin remplis (N>2) : %10d\n", filled_bin2);
      }