#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<\tc\med\medlib\paramete.c>
#include<\tc\recipes\utilitys\nr.h>
#include<\tc\recipes\utilitys\nrutil.h>
#include<\tc\recipes\convlv.c>

/*
#include<\tc\med\medlib\inverse.c>
#include<\tc\med\medlib\printa.c>
*/

/*	#define L_Max 50   hey dumkopf in paramete.c L = 50 (max poss L) */






#define N 1
#define M 128




#define simple_simulation 1
#define real_simple_simulation 2
#define realistic_simulation 3
#define jordan_data 4


/*
#define L_inc 1
#define L_begin 2
#define L_max 21		/*minus one */


#define L_inc 4
#define L_begin 5
#define L_max 6			/*minus one */


	FILE *id_trace;
	FILE *id_reflect;
	FILE *id_convolved;



	int L;

	double data[N][M];
	double dataX[N][M+L_Max];

	double data_max = 0.0;



main()
{
#include <\tc\include\dos.h>



	int i,j;
	int gra = 1;



	get_data(M,&dataX[0][0]);

	if (gra == 1)
	{
		for (i = 0; i < N; i ++)
			for (j = 0; j < M; j ++)

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

	}


	if(gra == 1)
	{
	plot_all(0,M,&data[0][0],
	"Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
	}

	get_old_data(M,&dataX[0][0]);

	if (gra == 1)
	{
		for (i = 0; i < N; i ++)
			for (j = 0; j < M; j ++)

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

	}


	if(gra == 1)
	{
	plot_all(0,M,&data[0][0],
	"Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
	}

	get_realistic_data(M,&dataX[0][0]);

	if (gra == 1)
	{
		for (i = 0; i < N; i ++)
			for (j = 0; j < M; j ++)

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

	}


	if(gra == 1)
	{
	plot_all(0,M,&data[0][0],
	"Source wavelet convolved with reflection signal, x(i) = w(i) * q(i)");
	}

/*	if (gra == 1)
	{
		for (j = 0; j < M; j ++)
			dataX[0][j] = dataX[0][j] * .7;
	}
*/


	exit(0);
}

get_data(int lngth, double *in_data)
{

	int wave_source = 0;			/* 0 => Cabrelli ex.1 */
						/* 1 => Gray ex. 1 */

	int i, j,k;
	double source_wavelet[M], spike_trace[M];

	double real_i, real_j, real_k, sum;

	double alpha,phi,time;
	double decay_rate;


	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise[M];

	double init_amp_0 =  1.0;
	double init_amp_1 =   .3880;
	double init_amp_2 =   .8274;
	double init_amp_3 =  -.6939;
	double init_amp_4 =   .9687;


	reflect_1 = lngth/2 * .1059;
	reflect_2 = lngth/2 * .1595;
	reflect_3 = lngth/2 * .3050;
	reflect_4 = lngth/2 * .6041;

	alpha = -18.778;
	alpha = -25;
	phi   = -.1960;
	decay_rate = 1.46;
	decay_rate = 5.5;
	decay_rate = 7.5;


	id_trace     = fopen("med2_tr1.dat","wt");
	id_reflect   = fopen("med2_rf1.dat","wt");
	id_convolved = fopen("med2_cv1.dat","wt");

	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/(lngth/2);
		printf("+");

		srand(j + 42);
		k = random(1000) - 500;
		real_k = k;
		noise[j] = (real_k/1000) * .02;

		if (j > lngth/2 + 1) noise[j] = 0;

					/* zero padding */

		source_wavelet[j] = 0;
		spike_trace[j] = 0;

		if(wave_source == 1) goto Gray1;

	/* create source wavelet (Cabrelli ex. 1) */

		source_wavelet[j] =
		(sin(alpha * time + phi) * exp(-decay_rate * time) *
		init_amp_0);

		if (j > lngth/2 + 1) source_wavelet[j] = 0;

					/* zero padding */

		fprintf(id_trace, "%le\t%e\n",
			(double)(j), source_wavelet[j]);


	/* create spike trace (as ex. 1) */

		if(j == reflect_1)
		spike_trace[j] = init_amp_1;

		if(j == reflect_2)
		spike_trace[j] = init_amp_2;

		if(j == reflect_3)
		spike_trace[j] = init_amp_3;

		if(j == reflect_4)
		spike_trace[j] = init_amp_4;

		if(wave_source == 0) goto Cab1;
		/* create wavelet as per Gray (ex. 1) */

Gray1:		source_wavelet[ 0] = 0.0;
		source_wavelet[ 1] = .033175;
		source_wavelet[ 2] = .076922;
		source_wavelet[ 3] = .165922;
		source_wavelet[ 4] = .52652;
		source_wavelet[ 5] = .80902;
		source_wavelet[ 6] = -.0039796;
		source_wavelet[ 7] = -.30637;
		source_wavelet[ 8] = -.30637;
		source_wavelet[ 9] = -.15650;
		source_wavelet[10] = .035809;
		source_wavelet[11] = .47878;
		source_wavelet[12] = -.17639;
		source_wavelet[13] = -.15385;
		source_wavelet[14] = -.013263;
		source_wavelet[15] = .31300;
		source_wavelet[16] = -.027852;
		source_wavelet[17] = -.14854;
		source_wavelet[18] = .014589;
		source_wavelet[19] = .13793;
		source_wavelet[20] = -.026526;
		source_wavelet[21] = -.12600;
		source_wavelet[22] = -.022547;
		source_wavelet[23] = .007957;


	/* create spike trace (as ex. 1) */

		spike_trace[ 1] = -1.00;
		spike_trace[ 8] =  -.582;
		spike_trace[22] =  -.900;
		spike_trace[28] =   .440;
		spike_trace[40] =  -.420;
		spike_trace[44] = -1.00;


Cab1:		fprintf(id_reflect, "%le\t%e\n",
			(double)(j), spike_trace[j]);


	/* create convolution of source and spike */

		sum = 0;
		for (i = 0; i <= j; i ++)
		{

			if (j-i >= 0 && j-i < lngth)
			sum =
			sum + source_wavelet[i] * (spike_trace[j-i]);
		}

		*(in_data + j) = sum + noise[j];

		fprintf(id_convolved, "%le\t%e\n",
			(double)(j), sum);

	}

/*	for (j = 0; j < lngth; j ++)
		*(in_data + j) = *(in_data + j) - sum;*/

	fclose( id_trace     );
	fclose(	id_reflect   );
	fclose(	id_convolved );

	plot_all(0,lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_all(0,lngth,&spike_trace,"Reflection signals, q(i)");
}

get_digitized_data(int lngth, double *in_data)
{

	int wave_source = 1;			/* 0 => Cabrelli ex.1 */
						/* 1 => Gray ex. 1 */

	int i, j,k;
	double source_wavelet[M], spike_trace[M];

	double real_i, real_j, real_k, sum;

	double alpha,phi,time;
	double decay_rate;


	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise[M];

	double init_amp_0 =  1.0;
	double init_amp_1 =   .3880;
	double init_amp_2 =   .8274;
	double init_amp_3 =  -.6939;
	double init_amp_4 =   .9687;


	reflect_1 = lngth/2 * .1059;
	reflect_2 = lngth/2 * .1595;
	reflect_3 = lngth/2 * .3050;
	reflect_4 = lngth/2 * .6041;

	alpha = -18.778;
	alpha = -25;
	phi   = -.1960;
	decay_rate = 1.46;
	decay_rate = 5.5;
	decay_rate = 7.5;


	id_trace     = fopen("med2_tr5.dat","wt");
	id_reflect   = fopen("med2_rf5.dat","wt");
	id_convolved = fopen("med2_cv5.dat","wt");

	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/(double)(lngth/2);
		printf("+");

		srand(j + 42);
		k = random(1000) - 500;
		real_k = k;
		noise[j] = (real_k/1000) * .02;

		if (j > lngth/2 + 1) noise[j] = 0;

					/* zero padding */

		source_wavelet[j] = 0;
		spike_trace[j] = 0;

		if(wave_source == 1) goto Gray1;

	/* create source wavelet (Cabrelli ex. 1) */

		source_wavelet[j] =
		(sin(alpha * time + phi) * exp(-decay_rate * time) *
		init_amp_0);

		if (j > lngth/2 + 1) source_wavelet[j] = 0;

					/* zero padding */

		fprintf(id_trace, "%le\t%e\n",
			(double)(j), source_wavelet[j]);


	/* create spike trace (as ex. 1) */

		if(j == reflect_1)
		spike_trace[j] = init_amp_1;

		if(j == reflect_2)
		spike_trace[j] = init_amp_2;

		if(j == reflect_3)
		spike_trace[j] = init_amp_3;

		if(j == reflect_4)
		spike_trace[j] = init_amp_4;

		if(wave_source == 0) goto Cab1;
		/* create wavelet as per Gray (ex. 1) */

Gray1:		source_wavelet[ 0] = 0.0;
		source_wavelet[ 1] = .033175;
		source_wavelet[ 2] = .076922;
		source_wavelet[ 3] = .165922;
		source_wavelet[ 4] = .52652;
		source_wavelet[ 5] = .80902;
		source_wavelet[ 6] = -.0039796;
		source_wavelet[ 7] = -.30637;
		source_wavelet[ 8] = -.30637;
		source_wavelet[ 9] = -.15650;
		source_wavelet[10] = .035809;
		source_wavelet[11] = .47878;
		source_wavelet[12] = -.17639;
		source_wavelet[13] = -.15385;
		source_wavelet[14] = -.013263;
		source_wavelet[15] = .31300;
		source_wavelet[16] = -.027852;
		source_wavelet[17] = -.14854;
		source_wavelet[18] = .014589;
		source_wavelet[19] = .13793;
		source_wavelet[20] = -.026526;
		source_wavelet[21] = -.12600;
		source_wavelet[22] = -.022547;
		source_wavelet[23] = .007957;


	/* create spike trace (as ex. 1) */

		spike_trace[ 1] = -1.00;
		spike_trace[ 8] =  -.582;
		spike_trace[22] =  -.900;
		spike_trace[28] =   .440;
		spike_trace[40] =  -.420;
		spike_trace[44] = -1.00;


Cab1:		fprintf(id_reflect, "%le\t%e\n",
			(double)(j), spike_trace[j]);


	/* create convolution of source and spike */

		sum = 0;
		for (i = 0; i <= j; i ++)
		{

			if (j-i >= 0 && j-i < lngth)
			sum =
			sum + source_wavelet[i] * (spike_trace[j-i]);
		}

		*(in_data + j) = sum + noise[j];


		fprintf(id_convolved, "%le\t%e\n",
			(double)(j), sum);

	}

/*	for (j = 0; j < lngth; j ++)
		*(in_data + j) = *(in_data + j) - sum;*/

	fclose( id_trace     );
	fclose(	id_reflect   );
	fclose(	id_convolved );

	plot_all(0,lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_all(0,lngth,&spike_trace,"Reflection signals, q(i)");
}


get_realistic_data(int lngth, double *in_data)
{

	int wave_source = 0;			/* 0 => Cabrelli ex.1 */
						/* 1 => Gray ex. 1 */

	int i, j,k;
	double source_wavelet[M], spike_trace[M];

	double real_i, real_j, real_k, sum;

	double alpha,phi,time;
	double decay_rate;


	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise[M];

	double init_amp_0 =  1.0;
	double init_amp_1 =   .990;
	double init_amp_2 =   .800;
	double init_amp_3 =   .0500;
	double init_amp_4 =   .500;
	double alpha2 = -15;
	double alpha3 = -35;
	double alpha4 = -55;
	double alpha5 = -20;
	double init_amp2_0 =  .2;


	reflect_1 = lngth/2 * .050;
	reflect_2 = lngth/2 * .400;
	reflect_3 = lngth/2 * .500;
	reflect_4 = lngth/2 * .700;

	alpha = -250;
	phi   = -.2;
	decay_rate = 15.0;

	id_trace     = fopen("med2_tr3.dat","wt");
	id_reflect   = fopen("med2_rf3.dat","wt");
	id_convolved = fopen("med2_cv3.dat","wt");


	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/(lngth/2);
		printf("+");

		srand(j + 42);
		k = random(1000) - 500;
		real_k = k;
		noise[j] = (real_k/1000) * .02;

		if (j > lngth/2 + 1) noise[j] = 0;

					/* zero padding */


		source_wavelet[j] = 0;
		spike_trace[j] = 0;

		if(wave_source == 1) goto Gray1;

	/* create source wavelet (Cabrelli ex. 1) */

		source_wavelet[j] =
		(sin(alpha * time + phi) * exp(-decay_rate * time) *
		init_amp_0);

		if (j > lngth/2 + 1) source_wavelet[j] = 0;

					/* zero padding */

		fprintf(id_trace, "%le\t%e\n",
			(double)(j), source_wavelet[j]);


	/* create spike trace (as ex. 1) */

		if(j == reflect_1)
		spike_trace[j] = init_amp_1;

		if(j == reflect_2)
		spike_trace[j] = init_amp_2;

		if(j == reflect_3)
		spike_trace[j] = init_amp_3;

		if(j == reflect_4)
		spike_trace[j] = init_amp_4;

		if(wave_source == 0) goto Cab1;
		/* create wavelet as per Gray (ex. 1) */

Gray1:		source_wavelet[ 0] = 0.0;
		source_wavelet[ 1] = .033175;
		source_wavelet[ 2] = .076922;
		source_wavelet[ 3] = .165922;
		source_wavelet[ 4] = .52652;
		source_wavelet[ 5] = .80902;
		source_wavelet[ 6] = -.0039796;
		source_wavelet[ 7] = -.30637;
		source_wavelet[ 8] = -.30637;
		source_wavelet[ 9] = -.15650;
		source_wavelet[10] = .035809;
		source_wavelet[11] = .47878;
		source_wavelet[12] = -.17639;
		source_wavelet[13] = -.15385;
		source_wavelet[14] = -.013263;
		source_wavelet[15] = .31300;
		source_wavelet[16] = -.027852;
		source_wavelet[17] = -.14854;
		source_wavelet[18] = .014589;
		source_wavelet[19] = .13793;
		source_wavelet[20] = -.026526;
		source_wavelet[21] = -.12600;
		source_wavelet[22] = -.022547;
		source_wavelet[23] = .007957;


	/* create spike trace (as ex. 1) */

		spike_trace[ 1] = -1.00;
		spike_trace[ 8] =  -.582;
		spike_trace[22] =  -.900;
		spike_trace[28] =   .440;
		spike_trace[40] =  -.420;
		spike_trace[44] = -1.00;

	/* create convolution of source and spike */

Cab1:		fprintf(id_reflect, "%le\t%e\n",
			(double)(j), spike_trace[j]);

		sum = 0;
		for (i = 0; i <= j; i ++)
		{

			if (j-i >= 0 && j-i < lngth)
			sum =
			sum + source_wavelet[i] * (spike_trace[j-i]);
		}

		real_j = j;
		time = real_j/(lngth/2);
		if (time > .65)
		{
			sum =
			((sin(alpha2 * time + phi) + sin(alpha3 * time)
			+ sin(alpha4 * time) + sin(alpha5 * time))/4*
			init_amp2_0);
		}

		*(in_data + j) = sum + noise[j];


		fprintf(id_convolved, "%le\t%e\n",
			(double)(j), *(in_data + j));
	}

/*	for (j = 0; j < lngth; j ++)
		*(in_data + j) = *(in_data + j) - sum;*/


	fclose( id_trace     );
	fclose(	id_reflect   );
	fclose(	id_convolved );

	plot_all(0,lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_all(0,lngth,&spike_trace,"Reflection signals, q(i)");
}


get_old_data(int lngth, double *in_data)
{


	int i, j;
	float source_wavelet[M], spike_trace[M], out_data[M];

	double real_i, real_j, sum;

	double alpha,phi,time;
	double decay_rate;

	int reflect_1, reflect_2, reflect_3,reflect_4;

	double noise = 0;

	double init_amp_0 =  1.0;
	double init_amp_1 =   .3880;
	double init_amp_2 =   .8274;
	double init_amp_3 =  -.6939;
	double init_amp_4 =   .9687;


	reflect_1 = lngth/2 * .1059;
	reflect_2 = lngth/2 * .1595;
	reflect_3 = lngth/2 * .3050;
	reflect_4 = lngth/2 * .6041;

	alpha = -18.778;
	phi   = -.1960;
	decay_rate = 1.46;



	id_trace     = fopen("med2_tr2.dat","wt");
	id_reflect   = fopen("med2_rf2.dat","wt");
	id_convolved = fopen("med2_cv2.dat","wt");

	for (j = 0; j < lngth; j ++)
	{
		real_j = j;
		time = real_j/(lngth/2);


	/* create source wavelet (Cabrelli ex. 1) */
		source_wavelet[j] =
		(sin(alpha * time + phi) * exp(-decay_rate * time) *
								init_amp_0);

		if (j > lngth/2 + 1) source_wavelet[j] = 0;

					/* zero padding */


		fprintf(id_trace, "%le\t%e\n",
			(double)(j), source_wavelet[j]);

	/* create spike trace (as ex. 1) */

		spike_trace[j] = .0000000;

		if(j == reflect_1)
		spike_trace[j] = init_amp_1;

		if(j == reflect_2)
		spike_trace[j] = init_amp_2;

		if(j == reflect_3)
		spike_trace[j] = init_amp_3;

		if(j == reflect_4)
		spike_trace[j] = init_amp_4;


		fprintf(id_reflect, "%le\t%e\n",
			(double)(j), spike_trace[j]);
	}

	/* create convolution of source and spike */

	convlv( source_wavelet, lngth, spike_trace, lngth/2, 1, out_data);

	for (j = 0; j  < lngth; j ++)
	{
		*(in_data + j) = (double)(out_data[j]);

		fprintf(id_convolved, "%le\t%e\n",
			(double)(j), *(in_data + j) );

	}

/*	for (j = 0; j < lngth; j ++)
		*(in_data + j) = *(in_data + j) - sum;*/


	fclose( id_trace     );
	fclose(	id_reflect   );
	fclose(	id_convolved );

	plot_all(0,lngth,&source_wavelet,"Source wavelet, w(i)");
	plot_all(0,lngth,&spike_trace,"Reflection signals, q(i)");
}