segdread.c

su 的源代码库

 	n_sk =      bcd ((unsigned char *) &segd_general_header_1.sk, 0, 2);
 	n_ec =      bcd ((unsigned char *) &segd_general_header_1.ec, 0, 2);
 	n_ex =      bcd ((unsigned char *) &segd_general_header_1.ex, 0, 2);

 	if (verbose==2) info_gh1(&segd_general_header_1);

 	/* Additional general headers */

 	if((n_gh == 2) && (segd_general_header_1.m[0] == 0x13)) { /* Special case for Sercel SN358 */
        if( EXIT_FAILURE == get_gn_sn358(&segd_gen_head_sn358, tapeun) ) break;

 		if (verbose==2) info_gn_sn358(&segd_gen_head_sn358);
	        if(hdr1_i != 0) {
			nsamp_hdr358 = bcd(&(segd_gen_head_sn358.rec_length[0]),1,3)*100 /* 10ths sec to msec */ * 16 /* base scan per msec */ /hdr1_i + 1;
          if(verbose) warn("nsamp_hdr358=%d\n",nsamp_hdr358);
			if (ns != 0 && nsamp_hdr358 != 0 && ns != nsamp_hdr358)
				warn("General Header 1 nsamp %u not equal to Sercel Header nsamp %u\n",nsamp_hdr1, nsamp_hdr358);
			if (nsamp_hdr358 != 0) ns = nsamp_hdr358;
		}
 	}
 	else {
 		for (i = 0; i < n_gh; i++) {
 			if (i == 0) {

 				/* General header #2 */

				if ( EXIT_FAILURE == get_gh2(&segd_general_header_2, tapeun) ) break;
				if ((segd_general_header_2.rev[0] <= 1) && (segd_general_header_2.rev[0] != 0)) { /* looks like SEGD rev 1 */
 				if (tr.fldr == BCD_FFFF) tr.fldr = 65536 * segd_general_header_2.ef[0] +
					256 * segd_general_header_2.ef[1] + segd_general_header_2.ef[2];
 				if (n_cs == BCD_FF) n_cs = 256 * segd_general_header_2.en[0] + segd_general_header_2.en[1];
 				if (n_ec == BCD_FF) n_ec = 256 * segd_general_header_2.ecx[0] + segd_general_header_2.ecx[1];
 				if (n_ex == BCD_FF) n_ex = 256 * segd_general_header_2.eh[0] + segd_general_header_2.eh[1];
 				n_gt = segd_general_header_2.gt;
 				if (verbose==2) info_gh2(&segd_general_header_2);
			 if (hdr1_r == BCD_FFF && hdr1_i != 0)
				nsamp_hdr2 = segd_general_header_2.erl[2] +
			               256 * (segd_general_header_2.erl[1] +
				       256 * (segd_general_header_2.erl[0])) *
				       16 / hdr1_i + 1;
			 else nsamp_hdr2 = 0;
          if(verbose) warn("nsamp_hdr2=%d\n",nsamp_hdr2);
			 if(nsamp_hdr2 != 0 && ns != 0 && ns != nsamp_hdr2)
				warn("General Header 2 nsamp %u differs from previous General Header(s) nsamp %u\n", nsamp_hdr2, ns);
			 if(nsamp_hdr2 != 0) ns = nsamp_hdr2;
			} else {
 				if (verbose==2) info_gh2(&segd_general_header_2);
			}
 			}
 			else {

 				/* General header #n */

				if ( EXIT_FAILURE == get_ghn(&segd_general_header_n, tapeun) ) break;
 				if (verbose==2) info_ghn(&segd_general_header_n);
 			}
 	 	}
 	}

 	/* Verify the length of the first record */ 
        if(!buff) nread = sftell(tapeun) - startpos;
 	if (buff && nread != ((1 + n_gh + n_str * (n_cs + n_sk) + n_ec + n_ex) * 32))
 		err("Error with length of first record\n"
 			"\t... first record = %d bytes differs from ((1 + n_gh + n_str * (n_cs + n_sk) + n_ec + n_ex) * 32)\n"
 			"\t    with n_gh=%d, n_str=%d, n_cs=%d, n_sk=%d, n_ec=%d, n_ex=%d\n",
 			nread, n_gh, n_str, n_cs, n_sk, n_ec, n_ex);

 	/* Allocate space for array csd */
	if (csd == NULL)
 		if ((csd = (channel_set_header **) alloc2 ((size_t) n_cs, (size_t) n_str, (size_t) sizeof(channel_set_header))) == NULL)
 			err("error at csd allocation");

 	/* if gain allocate space for mmp array */
 	if (gain && (mmp == NULL)) 
 		if ((mmp = (float **) alloc2float ((size_t) n_cs, (size_t) n_str)) == NULL)
 			err("error at mmp allocation");

 	/* For each scan type */

 	n_chan = 0;
 	for (i_scan = 0; i_scan < n_str; i_scan ++) {

 		/* For each channel set */

 		for (i_cs = 0; i_cs < n_cs; i_cs ++) {
			segd_channel_set_header = &csd[i_scan][i_cs];
			if ( EXIT_FAILURE == get_csh(segd_channel_set_header, tapeun) ) break;
 			n_chan += bcd((unsigned char*) &(segd_channel_set_header->cs), 0, 4);
 			if (gain) {
 				mp = ((float) ((segd_channel_set_header->mp[1] & 0x7f) << 8 | segd_channel_set_header->mp[0])) / 1024.;
 				if (segd_channel_set_header->mp[1] >> 7) mp *= -1.;
 				mmp[i_scan][i_cs] = pow ((double) 2., (double) mp);

 				/* For the seismic traces */
 				if (verbose && segd_channel_set_header->c == 0x10)
 					warn("Multiplier value for channel set %d of scan type %d is : %7.3e", i_cs, i_scan, mmp[i_scan][i_cs]);
 			}
 			if (verbose==2) info_csh(segd_channel_set_header);
 		}

 		/* Sample skew header */
 		for (i_ss = 0; i_ss < n_sk; i_ss ++) {
			if ( EXIT_FAILURE ==  get_ssh(&segd_sample_skew, tapeun) ) break;
 			if (verbose==2) info_ssh(&segd_sample_skew);
 		}
 	}

 	/* Extended Header */
 	for (j = 0; j < n_ec; j ++) {
		if( EXIT_FAILURE == get_ech(&segd_extended_header, tapeun) ) break;
		/* Local decoding */
 		/* (void) sscanf(&segd_extended_header[11], "%2hd:%2hd:%2hd", &tr.hour, &tr.minute, &tr.sec); */
 		if (verbose==2) info_ech(&segd_extended_header);
 	}

 	/* External Header */
 	for (j = 0; j < n_ex; j ++) {
		if( EXIT_FAILURE == get_exh(&segd_external_header, tapeun) ) break;
 		if (verbose==2) info_exh(&segd_external_header);
 	}
 	if (verbose==2) warn ("there are %d channels", n_chan);

 	/*************************
 	* Read the n_chan traces *
 	*************************/

 	for (i_tr=0; i_tr<n_chan; i_tr++) {

 		if (buff) {
 			if (-1 == (nread = (int) read(tapefd, (void *) bloc1, (size_t) REC_L))){
 				if (verbose)
 					warn("tape read error on trace %d", itr);
 				if (++errcount > errmax)
 					err("exceeded maximum io errors");
 			} else { /* Reset counter on successful tape IO */
 				errcount = 0;
 			}
		(void) sfseek(tapeun,startpos,SEEK_SET);
		if (!nread) break;  /* middle exit loop instead of mile-long while */
		}

 			if (!(segd_general_header_1.y & 0x8000)) { /* multiplexed data */
				/* decode data regarding the format */
 				switch (segd_general_header_1.y) { 
 					case 0x0015:                           /* 20 bit binary multiplexed */
 						err("Format 0015 (20 bit binary multiplexed) not yet implemented"); break;
 					case 0x0022:                        /* 8 bit quaternary multiplexed */
 						err("Format 0022 (8 bit quaternary multiplexed) not yet implemented"); break;
 					case 0x0024:                       /* 16 bit quaternary multiplexed */
 						err("Format 0024 (16 bit quaternary multiplexed) not yet implemented"); break;
 					case 0x0036:           /* 24 bit 2's compliment integer multiplexed */
 						err("Format 0036 (24 bit 2's compliment integer multiplexed) not yet implemented"); break;
 					case 0x0038:           /* 32 bit 2's compliment integer multiplexed */
 						err("Format 0038 (32 bit 2's compliment integer multiplexed) not yet implemented"); break;
 					case 0x0042:                       /* 8 bit hexadecimal multiplexed */
 						err("Format 0042 (8 bit hexadecimal multiplexed) not yet implemented"); break;
 					case 0x0044:                      /* 16 bit hexadecimal multiplexed */
 						err("Format 0044 (16 bit hexadecimal multiplexed) not yet implemented"); break;
 					case 0x0048:                      /* 32 bit hexadecimal multiplexed */
 						err("Format 0048 (32 bit hexadecimal multiplexed) not yet implemented"); break;
 					case 0x0058:                             /* 32 bit IEEE multiplexed */
 						err("Format 0058 (32 bit IEEE multiplexed) not yet implemented"); break;
 					case 0x0200:                                             /* illegal */
 						err("Format 0200 illegal, do not use"); break;
 					case 0x0000:                                             /* illegal */
 						err("Format 0000 illegal, do not use"); break;
 					default:
 						err("Data format code: %04x not recognized",segd_general_header_1.y);
 				}
 			}
 			else { /* demultiplexed data */
				nsamp_cs = 0; nsamp_the = 0;
				if( EXIT_FAILURE == get_dth(&segd_dem_trace_header, tapeun) ) break;
 				if (verbose==2) info_dth(&segd_dem_trace_header);
 				scan_type = bcd ((unsigned char *) &segd_dem_trace_header.st, 0, 2) -1;
 				chan_set = bcd ((unsigned char *) &segd_dem_trace_header.cn, 0, 2) -1;
				nsamp_the = 0;
				nsamp_cs = 0;
				if (csd[scan_type][chan_set].te != 0 && hdr1_i != 0)
					nsamp_cs = 2*(csd[scan_type][chan_set].te - csd[scan_type][chan_set].tf)*(16<<bcd(&csd[scan_type][chan_set].sc_j,0,1))/hdr1_i + 1 ;
				else nsamp_cs = ((ns-1)<<bcd(&csd[scan_type][chan_set].sc_j,0,1)) + 1 ;
				if (nsamp_cs != 0 && ns_override == 0) ns = nsamp_cs;
          if(verbose) warn("nsamp_cs=%d\n",nsamp_cs);


				for (i=0; i < segd_dem_trace_header.the; i++) { /* read the trace header extension blocks */
					if ( EXIT_FAILURE == get_the(&segd_trace_header_ext, tapeun) ) break;
				warn ("segd_dem_trace_header.the = %d", segd_dem_trace_header.the);
				if (i == 0) { nsamp_the = segd_trace_header_ext.nbs[2] + 256*(segd_trace_header_ext.nbs[1]+256*(segd_trace_header_ext.nbs[0]));
          if(verbose) warn("nsamp_the=%d\n",nsamp_the); }
				if (nsamp_the != 0 && nsamp_cs != 0 && nsamp_cs != nsamp_the)
				warn("Demux trace header nsamp %u != Demux trace header nsamp %u\n", nsamp_cs, nsamp_the);
				if (nsamp_the != 0) ns = nsamp_the;
				}


  				/* set trace identification code from channel type */
 				switch (csd[scan_type][chan_set].c) {
 					case 0xc0: tr.trid = TDUMMY; break; /* Auxiliary data trailer */
 					case 0x90: tr.trid = TDUMMY; break; /* Signature, filtered */
 					case 0x80: tr.trid = TDUMMY; break; /* Signature, unfiltered */
 					case 0x70: tr.trid = TDUMMY; break; /* Other */
 					case 0x60: tr.trid = TDUMMY; break; /* External data */
 					case 0x50: tr.trid = TIMING; break; /* timing traces */
 					case 0x40: tr.trid = WBREAK; break; /* Water break traces */
 					case 0x30: tr.trid = UPHOLE; break; /* Uphole traces */
 					case 0x20: tr.trid = TBREAK; break; /* Time break traces*/
 					case 0x10: tr.trid = TREAL;  break; /* Real time traces*/
 					case 0x00: tr.trid = TDUMMY; break; /* Unused */
 					default: tr.trid = TDUMMY; warn("channel type %02x unknown", csd[scan_type][chan_set].c); break;
 				}

 				tr.tracf = bcd ((unsigned  char *) &(segd_dem_trace_header.tn[0]), 0, 4);
 				tr.delrt = csd[scan_type][chan_set].tf * 2;

				if(ns_override != 0) ns = ns_override;

				/* decode data regarding the format */
 				switch (segd_general_header_1.y) { 
 					case 0x8015:                        /* 20 bits binary demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the)*4)/10; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*5/2);
						ns = tr.ns;
 						F8015_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8022:                      /* 8 bit quaternary demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/1; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*1);
						ns = tr.ns;
 						F8022_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8024:                      /* 16 bit quaternary demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/2; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*2);
						ns = tr.ns;
 						F8024_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8036:          /* 24 bit 2's compliment integer demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/3; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*3);
						ns = tr.ns;
 						F8036_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8038:          /* 32 bit 2's compliment integer demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/4; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*4);
						ns = tr.ns;
 						F8038_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8042:                      /* 8 bit hexadecimal demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/1; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*1);
						ns = tr.ns;
 						F8042_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8044:                     /* 16 bit hexadecimal demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/2; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*2);
						ns = tr.ns;
 						F8044_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8048:                     /* 32 bit hexadecimal demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/4; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*4);
						ns = tr.ns;
 						F8048_to_float (tapeun, (float *) tr.data, ns); break;
 					case 0x8058:                            /* 32 bit IEEE demultiplexed */
						if(buff) tr.ns = ((nread-20-32*segd_dem_trace_header.the))/4; /* number of samples from block length */
						else tr.ns = ns;
 						if (tr.ns != ns) 
						   warn("Tape read length for trace %d inconsistent with trace length computed from header(s)\n",nread,ns*4);
						ns = tr.ns;
 						F8058_to_float (tapeun, (float *) tr.data, ns); break;
 					default:
 						err("Data format code: %04x not recognized",segd_general_header_1.y);
 				}

  				/* Apply gain if requested */
 				if (gain) {
 					for (i = 0; i<tr.ns; i++)
 						tr.data[i] *= mmp[scan_type][chan_set];
 				}
				if (ipt >= ptmin-1) {
 				/* Write the trace */
 				/* if aux = 0, skip auxiliary channels */
 				if (aux!=0 || csd[scan_type][chan_set].c == 0x10)
 				puttr(&tr);

 				/* Echo under verbose option */
 				if (verbose && ++itr % vblock == 0)
 					warn(" %d traces from tape", itr);

				}
 			}
 	}

 	/***************************
 	* Read the General Trailer *
 	***************************/

 	for (j = 0; j < n_gt; j ++) {
		if( EXIT_FAILURE == get_gt(&segd_general_trailer, tapeun) ) break;
 		if (verbose==2) info_gt(&segd_general_trailer);
 	}
 
 	/* EOF = end of shot */
 	if (buff) {
 		if (-1 == (nread = (int) read(tapefd, (void *) bloc1, (size_t) REC_L))){
 			if (verbose)
 				warn("tape read error on header block from shot %d", (ipt+1));
 			if (++errcount > errmax)
 				err("exceeded maximum io errors");
 		} else { /* Reset counter on successful tape IO */
 			errcount = 0;
 		}
	(void) sfseek(tapeun, startpos, SEEK_SET);
 	if (nread) warn("not at EOF as should be!");
	}
 	ipt++ ;
 }

 /* Clean up */
 ipt = ipt - ptmin + 1;
 if (verbose) warn ("%d shots (%d traces) from tape", ipt, itr);
 (void) sfclose(tapeun);
 if(buff) eclose(tapefd);

 if (verbose) warn("tape closed successfully");
 if(mmp != NULL) free2float (mmp);
 if(csd != NULL) free2 ((void **) csd);
 if(bloc1 != NULL) free1 (bloc1);

 return EXIT_SUCCESS;
}

/* bcd - convert bcd to int
 *
 * Credits:
 *      EOPG: Marc, Jdt
 *
 * Parameters:
 *    ptr    - address of first byte of the number
 *    begin  - 0 or 1, position of the first digit of the number
 *    n      - number of digits
 *
 */

static int bcd (unsigned char * ptr , int begin , int n)
{
 register int i;
 unsigned int val;

 val = 0;
 if (n == 0) return (val);

 for (i = 0; i<n; i++) {
 	val *= 10;
 	if (begin++ & 1) val += (*ptr++ & 15); 
 	else val += (*ptr >> 4) & 15;
 }
 return (val);
}

/* F0015_to_float - convert 20 bit binary multiplexed data into floating numbers