selectsegy.c

地震记录微机格式读取

#include <malloc.h>
#include <sys/types.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include "stdafx.h"

//#define  NPOINT 1000
//#define	 NEWNPOINT  NPOINT/SAMP
//#define  SAMP 2

// the length of Reel Header (textual & binary header) in bytes
#define SEGY_LEN_TEXT 3200	// length of the Textual header (3200 byte) 
#define SEGY_LEN_BIN 400	// length of the Binary header (400 byte)
#define SEGY_LEN_REEL (SEGY_LEN_TEXT + SEGY_LEN_BIN) // length of the Reel header (3600 bytes)

#define SEGY_LEN_TRACE_HEAD 240	// length of the Trace header (240 bytes)

//void change_format4();
void sgy_swap_head(char* trHead);
//void sgy_swap_reel(char* reHead);
void sgy_swap_binHead(char* binHead);

int _tmain(int argc, char* argv[])
{
	union {
		char c[SEGY_LEN_TEXT];
		short s[1600];
		int i[800];
	}reel;
	union {
		char c[SEGY_LEN_BIN];
		short s[200];
		int i[100];
	}binHead;
	union {
		char c[SEGY_LEN_TRACE_HEAD];
		short s[120];
		int i[60];
	} hd;

int n,i,j,ntrace,ifile,SAMP;
long int length,NPOINT,NEWNPOINT;
float  sx,sy,rx,ry,data[10000],dis;
float  newdata[10000];
FILE   *fp,*fpo;
char  fnm[40],fnmo[40];
	printf("Input the name of the SEGY file :");
	scanf("%s",fnm);
	//sprintf(fnm,"2.sgy\0");
	printf("Input the SAMP=");
	scanf("%d",&SAMP);
	

	fp = fopen(fnm,"rb");
	fread(reel.c,1,SEGY_LEN_TEXT,fp);
	fread(binHead.c,1,SEGY_LEN_BIN,fp);	
	sgy_swap_binHead(binHead.c);
	NPOINT = binHead.s[10];
	fclose(fp);
	printf("\nin reel head : NPOINT=%d\n",NPOINT);

	fp = fopen(fnm,"rb");
	fseek(fp,0,SEEK_END);
	length = ftell(fp);
	fclose(fp);

	ntrace=(length-SEGY_LEN_REEL)/(NPOINT*4+SEGY_LEN_TRACE_HEAD);

	printf("%s %d\n",fnm,ntrace);

	fp = fopen(fnm,"rb");
	//fpo = fopen("hd.txt","wb");
	sprintf(fnmo,"CHHD_%s\0",fnm);
	fpo = fopen(fnmo,"wb");

	fread(reel.c,1,SEGY_LEN_TEXT,fp);
	fwrite(reel.c,1,SEGY_LEN_TEXT,fpo);
	fseek(fp,SEGY_LEN_TEXT, SEEK_SET);
	fread(binHead.c,1,SEGY_LEN_BIN,fp);
	
	sgy_swap_binHead(binHead.c);
	binHead.s[8]=binHead.s[8]*SAMP;
	binHead.s[10]=binHead.s[10]/SAMP;
	binHead.s[9]=binHead.s[9]*SAMP;
	binHead.s[11]=binHead.s[11]/SAMP;
	sgy_swap_binHead(binHead.c);
	fwrite(binHead.c,1,SEGY_LEN_BIN,fpo);

   for (i=0; i<ntrace; i++)
   {
     fread(hd.c,1,SEGY_LEN_TRACE_HEAD,fp);
	 fread(data,4,NPOINT,fp);

	 sgy_swap_head(hd.c);
	 //printf("hd.i[45](X coordinate)=%d, hd.i[46](Y coordinate)=%d\n",hd.i[45],hd.i[46]);
	 //fprintf(fpo,"hd.i[45](X coordinate)=%d, hd.i[46](Y coordinate)=%d\n",hd.i[45],hd.i[46]);
     hd.s[57]=hd.s[57]/SAMP;
     hd.s[58]=hd.s[58]*SAMP;
	 NEWNPOINT=NPOINT/SAMP;
	 if(i%10000==0) 
		 printf("%d of %d %d %d\n",i,ntrace,hd.s[57],hd.s[58]);

	 for (n=0; n<NEWNPOINT; n++)
	 {
		 newdata[n]=data[SAMP*n];
	 }
  
	sgy_swap_head(hd.c);
    fwrite(hd.c,1,SEGY_LEN_TRACE_HEAD,fpo);
	fwrite(newdata,4,NEWNPOINT,fpo);
		
  }
   fclose(fp);  fclose(fpo); 

}

  /*
 fp = fopen(fnm[0],"rb");        //
 fpo = fopen("RLine_1593","wb");     //

 fread(hd,4,900,fp);
 fwrite(hd,4,900,fpo);
 fclose(fp);

 for(ifile=0; ifile<8; ifile++)
 {
   fp = fopen(fnm[ifile],"rb");
   fseek(fp,0,SEEK_END);
   length = ftell(fp);
   fclose(fp);

   ntrace=(length-SEGY_LEN_TEXT)/(NPOINT*4+SEGY_LEN_TRACE_HEAD);
   nn = 0;

   fp = fopen(fnm[ifile],"rb");
   fread(hd,4,900,fp);
   for(i=0;i<ntrace;i++)
   {
     fread(hd,4,60,fp);
     fread(data,4,NPOINT,fp);
     change_format4(hd,60);

     if(hd[46]>1593000 && hd[46]<1593999)
     {
        change_format4(hd,60);
        fwrite(hd,4,60,fpo);
        fwrite(data,4,NPOINT,fpo);
        nn++;
     }
   }
   fclose(fp);
   printf("%s: total %d; output %d\n",fnm[ifile],ntrace,nn);
}
  fclose(fpo);
*/


//void change_format4(data,len)
//void change_format2(data,len)
//int *data;
//int len;
//{
//    int i;
//   register union {
//       int i ;
//       float f;
//       //unsigned char c[4];
//	   unsigned short s[2];
//   }temp;  
//   //register unsigned char ct;
//   register unsigned short ct;
//       
//
//     for ( i=0; i<len; i++ )
//	 {
//			temp.i = data[i];
//			ct=temp.s[0];
//			temp.s[0]=temp.s[1];
//			temp.s[1]=ct;
//			/*ct = temp.c[0];
//			temp.c[0] = temp.c[3];
//			temp.c[3] = ct;
//			ct = temp.c[1];
//			temp.c[1] = temp.c[2];
//			temp.c[2] = ct;*/
//			data[i] = temp.i;
//     }
//}

void swap_n_int (int dest[], const int src[], int n)
{
	register int i, iconv;

	for (i=0; i<n; i++)
	{
		iconv = src[i];
		iconv = (iconv<<24) | ((iconv>>24)&0xff) |
			((iconv&0xff00)<<8) | ((iconv&0xff0000)>>8);
		dest[i] = iconv;
	}
}

void swap_n_short(short data[], int n)
{
	/* in-place version */
	short *tempc = malloc(n * sizeof(short));
	if (! tempc) {
		fprintf(stderr, "alloc memory fail in  void swap_n_short(short*,int) !\n" );
		exit(1);
	}
	memcpy(tempc,data,n*sizeof(short));
	_swab ((char *)tempc, (char *)data, n*sizeof(short));
	if (tempc)
		free(tempc);
}

void sgy_swap_head(char* trHead)
/***********************************************************************
sgy_swap_head - Swap byte-order of a SEGY trace header
************************************************************************
Input::
trHead		input vector
************************************************************************
Note:
	the length of trHead is 240 bytes!
***********************************************************************/
{
	register union {
		short * s;
		int* i;
	} ptr;

	ptr.i = (int *)(trHead);//0

    swap_n_int(ptr.i,ptr.i,7);
	ptr.i += 7;

	//ptr.s=(short*)ptr.i;//28
	swap_n_short(ptr.s,4);
	ptr.s += 4;

	//ptr.i=(int*)ptr.s;//36
    swap_n_int(ptr.i,ptr.i,8);
	ptr.i += 8;

	//ptr.s=(short*)ptr.i;//68
	swap_n_short(ptr.s,2);
	ptr.s += 2;

	//ptr.i=(int*)ptr.s;//72
    swap_n_int(ptr.i,ptr.i,4);
	ptr.i += 4;

	//ptr.s=(short*)ptr.i;//88
	swap_n_short(ptr.s,46);
	ptr.s += 46;

	//--- changable ---
	//ptr.i=(int*)ptr.s;//180
    swap_n_int(ptr.i,ptr.i,5);
	ptr.i += 5;

	/*ptr.s=(short*)ptr.i;//200
	swap_n_short(ptr.s,2);
	ptr.s += 2;
	//204-239 ignore,swap as short anyway!
	swap_n_short(ptr.s,18);*/
    swap_n_int(ptr.i,ptr.i,10);
}

//void swapshort (short *si)
//{
///*  void swab(char *from, char *to, size_t nbytes); */
//	short tempc = *si;
//	_swab ((char *) &tempc, (char *) si, sizeof(short));
//}

/*void swapnshort (short *si,int n)
{
	short *tempc = malloc(n * sizeof(short));
	if (! tempc) {
		fprintf(stderr, "alloc memory fail in  void swapnshort(short*,int) !\n" );
		exit(1);
	}
	memcpy(tempc,si,n*sizeof(short));
	_swab ((char *)tempc,(char *)si,n*sizeof(short));
	if (tempc)
		free(tempc);
}

//void swaplong (int *li)
//{
// 
//	unsigned char temp;
//	unsigned char *ch = (unsigned char *)(li);
//	temp = ch[0]; ch[0] = ch[3]; ch[3] = temp;
//	temp = ch[1]; ch[1] = ch[2]; ch[2] = temp;
//}

void swapnlong (int *li, int n)
{
	register int i;
	unsigned char *ch;
	unsigned char temp;

	for (i=0; i<n; i++)
	{
		ch = (unsigned char *)(& (li[i]));
		temp = ch[0]; 
		ch[0] = ch[3]; 
		ch[3] = temp;

		temp = ch[1]; 
		ch[1] = ch[2]; 
		ch[2] = temp;
	}
}*/
void sgy_swap_binHead(char* binHead)
/***********************************************************************
sgy_swap_binHead - Swap byte-order of a SEGY binary header
************************************************************************
Input::
binHead		input vector
************************************************************************
Note:
	the length of binHead is 400 bytes!
***********************************************************************/
{
	register union {
		short * s;
		int* i;
	} ptr;

	ptr.i = (int*)(binHead);
	swap_n_int(ptr.i, ptr.i, 3);
	ptr.i += 3;

	swap_n_short(ptr.s, 24+170);
}


//void sgy_swap_reel(char* reHead)
///************************************************************************
//Note:
//	the length of reHead is 3600 bytes!
//***********************************************************************/
//{
//	register union {
//		short * s;
//		int* i;
//	} preel;
//
//	preel.i = (int *)(reHead);//0
//    swapnlong(preel.i,803);
//	preel.i += 803;
//
//	//ptr.s=(short*)ptr.i;//28
//	swapnshort(preel.s,194);
//	 
//}