/*	NOTES:

		21-SEP-90 -- discovered that monitoring the D-norm value
		that it may jump down suddenly, this is due to the walking
		off of a reflection peak off the end, note that the output
		is shifted by length of the filter.  As this length gets
		longer, real reflection peaks can walk off the end.
		This occured for (2) at 42 filters.

		22-SEP-90 -- modified to accept JORDAN data.
		25-SEP-90 -- modified to accept multip[le scans/input
		26-SEP-90 -- modified plot_all to do specific scaling
		12-JUN-94 -- modified to work with C++
*/


#include<math.h>
#include<values.h>
#include <sys\stat.h>

#include "winlib.h"
#include "paramete.h"
#include "meddprm.h"
#include "general.h"
#include "inverse.h"
#include "files.h"

#define scale_min 0
#define scale_max 0

#define true 1
#define false 0


/*	#define L_Max 50   hey dumkopf in paramete.c L = 50 (max poss L) */
#define N 1
#define M 250
/*#define M 96*/

#define part_of_M 21

#define L_inc 1
#define L_begin 15
#define L_max 	16		/*minus one */


#define beginning_scan 1
#define ending_scan 616


#define plot_nothing 0     	// define plot level
#define plot_only_end_results 1
#define plot_each_iteration 2
#define plot_algorithm_tracking 4
#define plot_everything 5
#define debug1 6






	char input_filename[44];
	char output_filename[44];
	char ddf_name[44];

	int plot_query = plot_nothing;
//	int plot_query = plot_only_end_results;
//	int plot_query = plot_algorithm_tracking;
	int save_error_file = true;

	int data_source = realistic_simulation;
	unsigned long int scan_number;

	int L;

	double data[N][M];
	double dataX[N][M + L_Max];
	double Y_optimum[N][M + L_Max];
	double Y[N][M + L_Max];
	double f[N][M + L_Max];
	double C[N][L_Max];
	double R[L_Max], W[L_Max][L_Max], F[L_Max];
	double D_norm_array[L_Max];

	double data_max = 0.0;

int    num_files;		   /* No. of data files in data set	    */
int    first_file, last_file;	   /* Nos. of 1st & last files displayed    */
int    filenum; 		   /* No. of file currently active	    */
int    ndps;			   /* No. of data per scan (constant)	    */


struct				   /* Info describing whole data set	    */
{
    char  site[16];		   /* An alphanumeric site designator	    */
    char  contractor[21];	   /* Organization that performed survey    */
    char  crew[21];		   /* Names or initials of survey crew	    */
    char  survey_type[3];	   /* Code indicating type of geophys meas  */
    char  sensor_mfgr_model[35];   /* Instrument mfgr and model(s)	    */
    char  serial_nums[35];	   /* Serial numbers of instruments listed  */
    char  nom_freq[11]; 	   /* Frequency of sensor signal	    */
    char  transport_mode[21];	   /* Manual or vehicle transport for surv  */
    char  dist_units[8];	   /* Units used for distance measurement   */
    char  x_axis_incr[11];	   /* Along-trk dist between measurements   */
    char  trk_spacing[11];	   /* Distance between trks		    */
    char  z_axis_scale_factor[35]; /* Time-depth conversion factor (GPR)    */
    char  data_units[11];	   /* Units associated with recorded data   */
    char  data_conv_factor[35];    /* Conversion factor to get final units  */
    char  grid_orientation[35];    /* Dir of the more northerly primary axis*/
    char  data_set_num[21];	   /* File no. assigned during the survey   */
    char  record_cust[35];	   /* Name of person who keeps data files   */
    char  record_loc[151];	   /* Place where the data are stored	    */
    char  num_files[5]; 	   /* The number of data files in this set  */
    char  data_per_scan[5];	   /* Same as setup.num_conv_per_scan	    */
    char  num_scans_avgd[5];	   /* No. of scans averaged per output scan */
} hdr;


struct				   /* Info describing individual files	    */
{
    char  trk_num[5];		   /* Track or record number		    */
    char  trk_loc[12];		   /* Y coordinate of trk rel to grid axis  */
    char  trk_begin[8]; 	   /* X coordinate of starting point (+/-)  */
    char  trk_end[8];		   /* X coordinate of stopping point (+/-)  */
    char  trk_len[8];		   /* Track length in field units (ft or m) */
    char  dir[2];		   /* Direction the track was run (Norm/Rev)*/
    char  scans[5];		   /* Number of scans in the track (file)   */
    char  status[2];		   /* Data status code (Good/Bad/Proc/Cal)  */
    char  hr[3];		   /* Time and date of data file creation   */
    char  min[3];
    char  month[3];
    char  day[3];
    char  year[3];
    char  low_thresh[4];	   /* Lower threshold of ampl histogram     */
    char  hi_thresh[4]; 	   /* Upper threshold of ampl histogram     */
    char comment[50];		   /* Optional string added during editing  */
} list[5];		   /* Array of structs for MAXTRKS trks     */


	void MEDD_deconvolution(HWND hWnd);
int read_ddf(void);




void medd(HWND hWnd)
{
        #include "common.h"
	#include <FCNTL.H>



	int i,j;
	int l;
        double data_local_max, data_local_min;


	FILE *id_medded;
        int id_binary;

	if (chdir("c:\\jordan\\data") != 0)
	{
		sprintf(print_buffer,"Error, cannot change directory\nfrom medd");
		stop_screen(hWnd, print_buffer);

	}


//  Specify input filename


	input_filename[0] = 0;
    	strcpy( input_filename, "jord34.003" );
// 	strcpy( input_filename, "jnk03.dat" );

// Delete and re-create any output files of the same name

    	output_filename[0] = 0;
	strcpy( output_filename, "jord34.001" );

	ddf_name[0] = 0;
	strcpy( ddf_name, "jord34.ddf" );


//read_and_plot_data_file(hWnd,1, 10, M, output_filename);


	_fmode = O_BINARY;
	if (id_binary = creat( output_filename, S_IWRITE ) == -1)
	{
		sprintf(print_buffer,"Error, cannot delete binary file %s:\n %s\n%d\n", output_filename, sys_errlist[errno], id_binary);
		stop_screen(hWnd, print_buffer);
	}
////	close(id_binary);








// Begin reading in scans, processing then and saving the processed data
// to disk

/* if( read_ddf() == -1 )
    {
	sprintf(print_buffer,"Descriptor file not found");
	stop_screen(hWnd, print_buffer);
    }

	sprintf(print_buffer,"%s",hdr);
	stop_screen(hWnd, print_buffer);
	sprintf(print_buffer,"%s",list[0]);
	stop_screen(hWnd, print_buffer);
 */

for (scan_number = beginning_scan - 1; scan_number <= ending_scan - 1; scan_number ++)
{
//   data_source = simple_simulation;
//   data_source = real_simple_simulation;
//   data_source = realistic_simulation;
   data_source = jordan_data;
//   data_source = gray_simulation;



   if (scan_number == 1)
   {

//	printf("simple simulation from modified Cabrelli ex. 1 (1)\n");
//	printf("real simple simulation from     Cabrelli ex. 1 (2)\n");
//	printf("realistic simulation of GRP siginal            (3)\n");
//	printf("jordan data, excerpt from a real data file     (4)\n");
//	printf("gray simulation; digitized data: Gray ex. 1    (5)\n");
//	printf("\nenter data source: ");
//
//	scanf("%d", &data_source);


	if (data_source == 4)
	{
//		printf("\nenter starting scan number: ");
//		scanf("%d", &scan_number);
//

// May want unique scaling for this selection:
//
//		scale_min = 100;
//		scale_max = 150;

	}
   }                      

if (plot_query >= plot_everything)
{
sprintf(print_buffer,"Data Source Selected: %d\nscan #%d",data_source,scan_number);
stop_screen(hWnd, print_buffer);
}

	get_array_from_disk(hWnd, scan_number, M, &data[0][0], input_filename, binary_file);
//	get_data_file(hWnd,N,M,&data[0][0],data_source,scan_number);


	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M; j ++)
		{
			dataX[i][j] = data[i][j];

// debug			dataX[i][j] = exp((double)(j)/170.0 ) * sin((double)j/100.0 * 8.0 * pi * scan_number * 24);

/*			printf("dataX[%d][%d] = %f\n",i,j,dataX[i][j]);*/
/*			printf(" data[%d][%d] = %f\n",i,j,data[i][j]);*/
		}

	}


	if (machine == tower && plot_query >= plot_only_end_results)
        {
		sprintf(print_buffer, "MEDD Input, x(i)");
		plot_all(hWnd,-N,M,&data[0][0],print_buffer,
			scale_min,scale_max);


	}

	MEDD_deconvolution(hWnd);

		for (j = 0; j < M+L; j ++)
		{
			dataX[0][j] = absolute_value(dataX[0][j]);
		}
//  Normalize data so the entire range fits from +1 to -1
		normalize((int) M, &dataX[0][L/2]);

// And amplify the reflections and shift it down
		for (j = 0; j < M+L; j ++)
		{
			dataX[0][j] = dataX[0][j] * 2.0 - 1.0 ;
		}


//	save_array_to_disk(hWnd, scan_number, M+L, &dataX[0][0], output_filename, binary_file);
	save_array_to_disk(hWnd, scan_number, M, &dataX[0][L/2], output_filename, binary_file);



// which do you really want to do?
//???	if (machine == tower && plot_query >= plot_everything)
//        {
//		sprintf(print_buffer, "MEDD Output, Y(i) = f * x(i)");
//		plot_all(hWnd,-N,M+L,&dataX[0][0],print_buffer,
//			-0.5,0.5);
//        }


	if (machine == tower && plot_query >= plot_only_end_results)
        {
		sprintf(print_buffer, "MEDD Output, Y(i) = f * x(i)");
		plot_all(hWnd,-N,M,&dataX[0][L/2],print_buffer, scale_min,scale_max);
    	}
		
		data_local_max = maximum_of_array( (int) (M * .9), &dataX[0][9]);
		data_local_min = minimum_of_array( (int) (M * .9), &dataX[0][9]);
//		data_local_min = maximum(absolute_value(data_max_x),absolute_value(data_min_x));



//		sprintf(print_buffer, "Blow up of MEDD Output, Y(i) = f * x(i)");
//		plot_all(hWnd,-N,M+L,&dataX[0][0],print_buffer, data_local_min, data_local_max);


   }

	sprintf(print_buffer,"MEDD has Completed %d scans", scan_number );
	stop_screen(hWnd, print_buffer);


}


void MEDD_deconvolution(HWND hWnd)
{
	#include "common.h"

	int ickx,icky;

	int i,j,k,s,u;
	double expect_Y;
	double D_norm, D_norm_final;
	double medd_filter, medd_output;
	double optimum_D_norm, optimum_output;
	int optimum_L = L_begin;

	FILE *id_err;

	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L_Max - 1; j ++)
   		{
			Y_optimum[i][j] = 0;
		}

	}

	for (L = L_begin; L < L_max; L = L + L_inc)
	{

	D_norm_final = 0;

					/* calculate auto-corelation
					coefficients of input traces */


	for (k = 0; k < L; k ++)
	{

		for (i = 0; i < N; i ++)
		{
			C[i][k] = 0;
			for (j = 0; j < M; j ++)
			{
				if( j + k  >= 0 && j + k < M)
				{
		     
if (j+k < 0 || j+k > M + L_Max)
{
sprintf(print_buffer, "in C loop: i %d,j %d,j+k %d\n",i,j,j+k);
//sprintf("C: dataX[i][j + k] = dataX[%d][%d] = %f\n",i,j+k,dataX[i][j+k]);
stop_screen(hWnd, print_buffer);
}

				    C[i][k] =
				    C[i][k] + dataX[i][j] * dataX[i][j+k];
				}
			}

//sprintf(print_buffer,"C[%d][%d] = %f\n",i,k,C[i][k]);
//stop_screen(hWnd, print_buffer);

		}

	}


	for (k = 0; k < L; k ++)
	{
		R[k] = 0;
		for (i = 0; i < N; i ++)
		{
   			R[k] = R[k] + C[i][k];

		 }

//sprintf(print_buffer,"R[%d] = %f\n",k,R[k]);
//stop_screen(hWnd, print_buffer);

	}







					/* inversion of autocorrelation
					matrix R[I-J] using
					Levinson Recursion:
					     -1
					W = R   (L x L) - matrix.	*/



	invert_matrix(L,L-1,&R[0],&R[1],&R[1],&W[0][0]);


//					For debugging purposes
/*
for (ickx = 0; ickx < L; ickx++)
{
for(icky = 0; icky < L+1; icky++)
{
printf("W[%d][%d] = %f\n",ickx,icky,W[ickx][icky]);
}
}
*/
//					End of debugging


	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
		{

					/* calculation of [I] [J]
					filter, F	*/

			for (k = 0; k < L; k ++)
			{
				F[k] = 0;
				for (s = 0; s < L; s ++)
				{
					if ( j - s >= 0 && j - s < M)
					{

/* this assumes that j-s is terminative not cyclic*/


//		printf("dataX[i][j-s] = dataX[%d][%d] = %#g\n",
//		i,j-s,dataX[i][j-s]);
if (j-s < 0 || j-s > M + L_Max)
{
sprintf(print_buffer,"s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]:\n %d %d %d %d %f, %f, %f\n",
s,i,j-s,k,W[s][k],dataX[i][j-s],F[k]);
stop_screen(hWnd, print_buffer);
}

					F[k] = F[k] + dataX[i][j-s] * W[s][k];
					}




				}
//printf("\tFilter F[%d] = %#g\n",k,F[k]);
			}

/*
for (icky = 0; icky < M; icky ++)
printf("before s-k bussiness:	dataX[0][%d] = %#g\n",
icky,dataX[0][icky]);
*/


					/* calculation of the matrix Y
					corosponding to the output of the
					[I] [J] filter, F	*/

			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M+L-1; s ++)
				{
					Y[k] [s] = 0;
					for (u = 0; u < L; u ++)
					if (s - u + 1 >= 0 && s - u + 1 < M)
					{
/*printf("dataX[k][s-u+1] = dataX[%d][%d] = %f\n",k,s-u+1,dataX[k][s-u+1]);*/

						Y[k][s] =
					Y[k][s] + F[u] * dataX[k][s-u+1];
					}

//sprintf(print_buffer, "Y[%d][%d] = %#g\n",k,s,Y[k][s]);
//stop_screen(hWnd, print_buffer);
				}
			}
						/* Calculation of ||Y|| */

/*
for (ickx = 0; ickx < M+L-1; ickx ++)
printf("\t\t\tY[%d][%d] = %#g\n",0,ickx,Y[0][ickx]);
*/

			expect_Y = 0;
			for (k = 0; k < N; k ++)
			{
				for (s = 0; s < M + L - 1; s ++)
					expect_Y =
					expect_Y + Y[k][s] * Y[k][s];

//sprintf(print_buffer,"\t\t\tk = %d, j = %d, M + L = %d, ||Y|| = %#g\n",
//k,j,M+L,expect_Y);

			 }
			expect_Y = sqrt(expect_Y);



						/* Calculation of the [I][J]
						coefficients of the [I][J]
						output normalized	*/

			if(expect_Y != 0.0)
				D_norm = (absolute_value(Y[i][j]))/(expect_Y);
			else
				D_norm = 0;


			if (D_norm > D_norm_final)
			{
				D_norm_final = D_norm;
/*				medd_filter  = F[?];*/
				medd_output  = expect_Y;

			}


		}
	}

/*
// I don't know where this came from...???
// but it was commented out for as long as I can remember
 
	D_norm_array[L] = 0;
	for (i = 0; i < N; i ++)
	{
		for (j = 0; j < M + L - 1; j ++)
			D_norm_array[L] = D_norm_array[L] + Y[i][j] * Y[i][j];
	}
*/


	D_norm_array[L] = D_norm_final;



/* for (ickx = 0; ickx < M+L-1; ickx ++)
printf("Y[%d][%d] = %f\n",0,ickx,Y[0][ickx]);*/

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		if (j < L) f[0][j] = F[j];
	}

/*  Don't really need...
	if (machine == tower && plot_query >= plot_everything)
	{	sprintf(print_buffer, "MEDD filter values, f");
		plot_all(hWnd,-1,M,&f,print_buffer,scale_min,scale_max);
        }

		sprintf(print_buffer, "MEDD Input, x(i)");
//		plot_all(hWnd,-N,M,&data[0][0],print_buffer, scale_min,scale_max);
*/

	for (j = 0; j < M; j ++)
	{
		f[0][j] = 0;
		f[0][j] = Y[0][j];
//sprintf(print_buffer, "Y[0][j] = %#g\n", Y[0][j]);
//stop_screen(hWnd, print_buffer);

	}


	if (plot_query >=  plot_only_end_results)
        {
	  sprintf(print_buffer, "RESULTS for this Filter:\nD_norm: %f\nMEDD Output: %f\nCurrent number of Filters: %d\nError Factor: %#g\n",
	  D_norm_final, medd_output,L,D_norm_array[L]);
	  stop_screen(hWnd, print_buffer);

			sprintf(print_buffer, "MEDD output for this filter, f");
			plot_all(hWnd,-1,M+L,&f,print_buffer,scale_min, scale_max);


	}


///  Use this if unknown machine causes portability problems.
///  Also, can save a file of filter (L) information here.
//  
//	if (machine != unknown && plot_query >= plot_everything)
//	{
//		sprintf(message_text,"  Continue?");
//  		stop_screen(hWnd, message_text);
//
//	}
// ---------------------


		 if (L == L_begin) goto label1;

		 if (D_norm_array[L] > D_norm_array[optimum_L])
		 {
label1:		 optimum_L = L;
		 optimum_D_norm = D_norm_final;
		 optimum_output = medd_output;
			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + L - 1; j ++)
				{
				Y_optimum[i][j] = Y[i][j];
//sprintf(print_buffer, "Y_optimum[i][j] = %#g\n", Y_optimum[i][j]);
//stop_screen(hWnd, print_buffer);
				}
			}
		}
	}



// End of L loop -----------


	i = 0;

 if (machine == unknown || save_error_file == true)
		id_err = fopen("medderr.dat", "wt");

	for (L = L_begin; L < L_max; L = L + L_inc)
	{
		if (D_norm_array[L] != 0)
		{
			f[0][i] = log10(D_norm_array[L]);
                }
		else
		{
			sprintf(print_buffer, "Warning!  D_norm is zero.  Chech you data.\n");
			stop_screen(hWnd, print_buffer);
		}

		if (machine == unknown || save_error_file == true)
			fprintf(id_err, "%e\t%e\n", (double)(L), D_norm_array[L]);

		i++;
	}

	if (machine == unknown || save_error_file == true)
		fclose(id_err);


	if (machine == tower && plot_query >= plot_each_iteration)
	{
		sprintf(print_buffer, "Log of Error vs. filter number");
		plot_all(hWnd,-1,i,&f,print_buffer, scale_min, scale_max);
	}

if (plot_query >= plot_everything)
{
	sprintf(print_buffer, "FINAL RESULTS:\nD_norm: %f\nMEDD Output: %f\nBest Number of Filters: %d\nBest Error Factor: %#g\n",
	optimum_D_norm,optimum_output,optimum_L,D_norm_array[optimum_L]);
	stop_screen(hWnd, print_buffer);
}
///  Use this if unknown machine causes portability problems.
///  Also, can save a file of filter (L) information here.
//  
//	if (machine != unknown && plot_query >= plot_everything)
//	{
//		sprintf(message_text,"  Continue?");
//  		stop_screen(hWnd, message_text);
//	}
// ---------------------


			for (i = 0; i < N; i ++)
			{
				for (j = 0; j < M + optimum_L - 1; j ++)
				{
					dataX[i][j] = 0;
					dataX[i][j] = Y_optimum[i][j];

//sprintf(print_buffer, "dataX = %#g\n", dataX[i][j]);
//stop_screen(hWnd, print_buffer);

				}
			}

sprintf(message_text," End of Process.  Continue?");
//stop_screen(hWnd, message_text);

	return;
}





/*
------------------------------------------------------------------------------
---------------------------- FUNCTION READ_DDF -------------------------------

This function reads a data descriptor file (DDF) from disk if that file is
present.  The file should be present if the file descriptor information has
been previously edited.  If it is not present, the value -1 is returned to the
calling function.  The track-specific information read from the file is stored
in structure list[].

Note that the information relating to the first track, or data file, is placed
in list[1] rather than in list[0].

------------------------------------------------------------------------------
*/
int read_ddf(void)
{
    int status;
    int fh_ddf;                                   /* File handle for DDF    */

    fh_ddf = open( ddf_name,  O_BINARY );
    if( fh_ddf == -1 ) return -1;                 /* File does not exist    */

 /* ----------
    Read the DDF header.  This contains data that pertain to the entire data
    set.                                                                    */
// sizeof( hdr )
//  read( fh_ddf, &hdr, sizeof( hdr ) );
//    num_files = atoi( hdr.num_files );
//    ndps = atoi( hdr.data_per_scan );

 /* ----------
    Read the rest of the descriptor data.  Store the data in structure
    list[].                                                                 */

//    read( fh_ddf, &list[1], num_files * sizeof( list[1] ) );
    close( fh_ddf );
    return 0;
}