#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#define length 130
#define NShot 1
#define Channel 40
#define Sample_Point 2000


static float data[NShot][Channel][Sample_Point];

void ricker(int Lr,int L,float dalt,float fp,float *w)
{
	double temp;
	double pi;
	int Ll;
	int k,i;
	Ll = L - Lr - 1;
	pi = 4.0 * atan(1.0);
	for(i = -Lr; i <= Ll; i++)
	{
		k = i + Lr;
		temp = pi * fp * (dalt * i);
		temp *= temp;
		w[k] = (float) ( (1.0 - 2.0 * temp) *  exp(-temp) );
		//printf("%5d %10.6f\n",k,w[k]);

	}
	return;
}

int making_data ( FILE *fp,int num_trace, int sample, int num_shot )
{
	int i,j,k,m;
	
	static float w[130];

	

	if((fp = fopen("data.txt","wb")) == NULL)
	{
		printf("----cannot open  data file --\n");
		return 0; 
	}

		
	ricker(30,100,0.001,30.0,w);

	for ( i=0;i<num_shot;i++ ) {
		for ( j=0;j<num_trace;j++) {
			for ( k=0;k<sample;k++) {
				 data[i][j][k]=0;
			}
			fwrite(&data[j+j*i], sizeof(float), sample, fp);
		}
	}

	fclose(fp);
}


int out_segy(FILE *fp,
 			float* trace,
			 char* outsgyfile, 
			 int xline_s, int xline_e, 
			 int n_sample, int T_rate, int nshot)
{

	FILE *outsgy;
	char cha[80];
	int istart, i, iy, iadd;
	short trhdS[200], idhdS[200] ;
	long  trhdL[200], idhdL[200], ibyte;
	float trhdR[200];

	for(i = 0; i < 200; i++)
	{
		trhdL[i] = 0;
		trhdS[i] = 0;
		trhdR[i] = 0;
	}


	if((outsgy = fopen(outsgyfile,"wb")) == NULL)
	{
		printf("----cannot open segy datafile in output segy--\n");
		return 0; 
	}

	//對前面的3200個字節(jié)進行寫入
	fseek(outsgy, 0L, 0);
	for(i = 0; i < 80; i++) cha[i] = ' ';
	cha[0] = 'C';
	cha[5] = 'k';
	for(i = 0; i < 40; i++)
	{
		fseek(outsgy, i * 80L, 0);
		fwrite(cha,1, 80, outsgy);
	}
	fseek(outsgy, 3200L, 0);


	//開始操作400個字節(jié)
	idhdL[0] = 10;
	idhdL[1] = 99;
	idhdL[2] = 99;

	fwrite(idhdL, 4, 3,outsgy);


	fseek(outsgy,3212L,0);

	for(i = 0;i < 194; i++)
	{
		idhdS[i] = 0;
	}
	/* short idhds */
	idhdS[0]  = (xline_e - xline_s + 1)*nshot;                /*3213-3214 number 0f data trace per record   */
	idhdS[1] =  0;                 /* 16 number of auxuliary traces  */
	idhdS[2] = T_rate*1000;        /* 18 sample rate in microsecond  */
	idhdS[3] = T_rate*1000;        /* 20 sample rate in microsecond  */
	idhdS[4] = n_sample;            /*22 sample number per trace     */
	idhdS[5] = n_sample;            /*24  sample number per trace    */
	idhdS[6] = 1;         /* 26 data format 1:float point(4)         */
	idhdS[7] = 1;         /* 28 cdp fold          */
	idhdS[8] = 1;        /* 30 trace sort 1: no sort    */
	idhdS[9] = 1;        /* 32 vertically stacked code 1: no sum  */
	idhdS[10] = 3;       /* 34 sweep frequence at start   */
	idhdS[11] = 60;      /* 36 sweep frequence at end     */
	idhdS[12] = 300;     /* 38 sweep length(in ms)        */
	idhdS[13] = 1;       /* 40 sweep type code 1:linear  */
	idhdS[14] = 1;      /* 42 trace number per sweep channal */
	idhdS[15] = 20;     /* 44 sweep trace taper length in ms at start  */
	idhdS[16] = 20;     /* 46 sweep trace taper length in ms at end    */
	idhdS[17] = 1;      /* 48 taper type 1:linear  */
	idhdS[18] = 1;      /* 50 correlated data trace 1:no 2:yes    */
	idhdS[19] = 1;      /* 52 binary gain recorved 1:y 2:n   */
	idhdS[20] = 1;      /* 54 amplitude recovery method 1:none  */
	idhdS[21] = 1;      /* 56 measurement system 1: metre 2:feet  */
	idhdS[22] = 1;      /* 58 input signal polarity   */
	idhdS[23] = 1;      /* 60 vibratory polarity code  */

	fwrite(idhdS, 2, 194, outsgy);
	fseek(outsgy, 3600L, 0);

/*-----------------  Here write segy traces data  */

	iadd = 0;
	rewind(fp);
	for (int ishot = 0; ishot < nshot; ishot++)
	{
	for(iy = xline_s; iy <= xline_e; iy++)
	{
		istart =  (iy - xline_s) + (xline_e - xline_s + 1)*ishot;
		ibyte = istart * ( n_sample * sizeof(float) + 240) + 3600;
		iadd = iadd+1;

		trhdL[0] = iadd;	  // TRACENO  
		trhdL[1] = iadd;	  // trace number in the reel  
		trhdL[2] = 1;	  //  ffid  
		trhdL[3] = iadd;		  // chan  
		trhdL[4] = iadd;	  // source number  
		trhdL[5] =  1 + istart;		// cdp number 


		fseek(outsgy, ibyte, 0);
		fwrite(trhdL, 4, 7, outsgy);

		trhdS[0] = 1;     // 道識別碼 

		fseek(outsgy, ibyte+28, 0);
		fwrite(trhdS, 2, 4, outsgy);

		for(i = 0; i < 68; i++)trhdS[i] = 0;
		//對應這一道的樣點數(shù)和采樣間隔,對應115和117
		trhdS[5] = n_sample;
		trhdS[6] = T_rate*1000;

		fseek(outsgy, ibyte+104L, 0);
		fwrite(trhdS, 2, 68, outsgy);

		fseek(outsgy, ibyte+240L, 0);
		fread(trace, sizeof(float), n_sample, fp);
		fwrite(trace, sizeof(float), n_sample, outsgy);

	}
	rewind(fp);
	}

	fclose(outsgy);

	return 1;

}
main()
{
	int i;

	making_data( Channel, 2000,1, Nshot );

	out_segy(out,data,segy.txt, 1, Channel, 2, 0.001, Nshot);
	


}