output_segy.c

解吸SEED格式的源代码

		/* find min and max values in trace */	
		segy_header.comment.p_stuff.max_samp = 0xF0000000;
		segy_header.comment.p_stuff.min_samp = 0x7FFFFFFF;

                j = k * (seis_buffer_length/data_hdr->num_mux_chan);

		for (i = 0; i < data_hdr->nsamples; i++)
		{
			if (seismic_data[offset+i+j] > segy_header.comment.p_stuff.max_samp)
				
                	segy_header.comment.p_stuff.max_samp = seismic_data[offset+i+j];

			if (seismic_data[offset+i+j] < segy_header.comment.p_stuff.min_samp) 
                                 
			segy_header.comment.p_stuff.min_samp = seismic_data[offset+i+j];
	
		}

		reverse_flag = 0;

		for (i=0;i<10;i++) 
			channel[i] = 0;

		strcpy(channel,current_channel->channel);

        	if (data_hdr->num_mux_chan > 1) 
		{
                	int i;

                	for (i=0; i<3; i++) {
                        	if (channel[i] == ' ' || 
				    channel[i] == '\0') 
				break;
                	}

                	if (i > 2) 
				i = 2;

                	channel[i] = data_hdr->mux_chan_name[k];
        	}


		reverse_flag = check_for_reversal(current_channel->dip,
						  current_channel->azimuth,
						  channel[2]);


		flag = FALSE;
        	for (response = current_channel->response_head;
                            response != NULL; response = response->next)
        	{
                	reverse_flag &= 1;
                
			if ((response->type != 'S') && response-> type != 'R')
				continue;

			if (response->type == 'S') {
               			if ((response->ptr.type58)->stage != 0)
					continue;

                        	if (((response->ptr.type58)->sensitivity < 0.0) &&
                                	(check_reverse & 2))
                                	reverse_flag |= 2;
 
                        	flag = TRUE;
			} else if (response->type == 'R') {
				type_48 = find_type_48_stage_0(response);
				if (type_48 == NULL) continue;

				if ((type_48->sensitivity < 0.0) && 
					    (check_reverse & 2))
						reverse_flag |= 2;

				flag = TRUE;
			}

                }
        
        	if (!flag && (check_reverse & 2))
                	fprintf(stderr, "Warning - Stage Zero Not Found on Gain Reversal check\n");
 
        	if (check_reverse & 1)
        	{
                	if (reverse_flag & 1)
                        	fprintf (stderr, "Warning... Azimuth/Dip Reversal found %s.%s, Data will be automatically inverted\n           Response Information will not be modified\n",
                        current_station->station, channel);
        	}
        	else
        	{
                	if (reverse_flag & 1)
                        	fprintf (stderr, "Warning... Azimuth/Dip Reversal found %s.%s, Data inversion was not selected\n", 
					current_station->station, channel);
                	reverse_flag &= 2;

        	}
 
        	if (check_reverse & 2)
        	{
                	if (reverse_flag & 2)
                        	fprintf (stderr, "Warning... Gain Reversal found %s.%s, Data will be automatically inverted\n           Response Information will not be modified\n",
                                current_station->station, channel);
        	}
        	else
        	{
                	if (reverse_flag & 2)
                        	fprintf (stderr, "Warning... Gain Reversal found %s.%s, Data inversion was not selected\n",
                                		current_station->station, channel);
                	reverse_flag &= 1;
        	}
       
		if (reverse_flag == 3)
        	{
                	reverse_flag = 0; 	/* Double inversion */
                	fprintf (stderr, "Warning... Gain and Dip/Azimuth Reversal found %s.%s, Data was not inverted\n",
                        	current_station->station, channel);
        	}
 
        if (reverse_flag & 2) 
		segy_header.comment.p_stuff.scale_factor = -segy_header.comment.p_stuff.scale_factor;
		
/*                 +=======================================+                 */
/*=================|  build name for and open output file  |=================*/
/*                 +=======================================+                 */

	/* 1988.023.15.34.08.2800.ANMO.SPZ is the filename for a seismogram 
  	 *	from year 1988, Julian day 23, 15:34:08.2800 UT, recorded 
	 *	at station ANMO from component SPZ. 
	 */

		/* if requested from WEED, place into subdirectory */
		if (EventDirFlag)
		{
			char dirname[MAXPATHLEN];

			strcpy(dirname, get_event_dir_name());

			if (!dir_exists(dirname))
                	{
                        	if (mkdir(dirname, 0777) == -1)
                        	{
                                	fprintf(stderr, "Unable to create event directory. Using current directory!\n");
                                
					perror("output_sac");
                                	strcpy(dirname, "./");
                        	}
 
                	}
		
			chdir(dirname);
 
		}

		sprintf (outfile_name, 
				"%04d.%03d.%02d.%02d.%02d.%03d.%s.%s.%s.%s.%c.SEGY",
			data_hdr->time.year,
			data_hdr->time.day,
			data_hdr->time.hour,
			data_hdr->time.minute,
			data_hdr->time.second,
			data_hdr->time.fracsec/10,
			type10.version >= 2.3 ? data_hdr->network : "",
			data_hdr->station,
			data_hdr->location,
			channel,
			input.type);

		if ((outfile = fopen (outfile_name, "w+")) == NULL)
		{
			fprintf (stderr, "\tWARNING (output_segy):  ");
			fprintf (stderr, "Output file %s is not available for writing.\n", outfile_name);

			perror("output_segy()");

			fprintf (stderr, "\tExecution continuing.\n");

			chdir(orig_dir);
		
			return;
		}

/*                 +=======================================+                 */
/*=================|       write a SEGY binary file        |=================*/
/*                 +=======================================+                 */

		/* describe the file being written */
		printf("Writing %s: %s: %s: %s, %5d samples (binary),",
			data_hdr->network,
			data_hdr->station,
			data_hdr->location,
			channel,
			data_hdr->nsamples);

		printf(" starting %04d,%03d %02d:%02d:%02d.%03d UT\n",
			data_hdr->time.year,
			data_hdr->time.day,
			data_hdr->time.hour,
			data_hdr->time.minute,
			data_hdr->time.second,
			data_hdr->time.fracsec/10);

		if (reverse_flag)
		{
			int tmpx;

			tmpx = segy_header.comment.p_stuff.max_samp;
			segy_header.comment.p_stuff.max_samp =
				segy_header.comment.p_stuff.min_samp;

        		segy_header.comment.p_stuff.min_samp = tmpx;
	
		}

		/* write the SAC header */
		if (fwrite (&segy_header, sizeof(struct segy), 1, outfile) != 1)
		{
			fprintf (stderr, "\tWARNING (output_segy_header):  ");
			fprintf (stderr, 
				"failed to properly write SEGY header to %s.\n",
					outfile_name);

			perror("output_segy()");

			fprintf(stderr, "\tExecution continuing.\n");

			fclose(outfile);

			chdir(orig_dir);

			return;

		}

		/* write the SEGY data */

		/* find index into start of data in buffer, if multiplexed
		 * it starts at locations other than zero
		 */
		j = k*(seis_buffer_length/data_hdr->num_mux_chan);

                if (reverse_flag)
                {
                        char wrkstr[200];

                        for (i=0;i<data_hdr->nsamples;i++)
                                seismic_data[offset+i+j] = -seismic_data[offset+i+j];
                        sprintf(wrkstr, "Data reversal initiated for SEGY file %s", outfile_name);
 
                        rdseed_alert_msg_out(wrkstr);
                        
                }
 

		for (i = 0; i < data_hdr->nsamples; i++)
		{
			this_sample = (int) seismic_data[offset+j+i];

			if (fwrite (&this_sample, sizeof(int), 1, outfile) != 1)
			{
				fprintf(stderr, "\tERROR - output_segy():failed to properly write SEGY data to %s.\n", outfile_name);

				fprintf(stderr, "Giving up on this file\n");

				perror("output_segy()");
			
				i = data_hdr->nsamples; /* drop out of loop */
			}

		}		/* for i < num_samples */

		fclose(outfile);

	} 		/* the big for loop - if multiplexed */


	chdir(orig_dir);
}

/* -------------------------------------------------------------------- */
struct type58 *find_type_58_stage_0();

struct type48 *find_type_48();


float get_overall_gain()

{
	struct type58 *t_58;
	struct type48 *t_48;

	t_58 = find_type_58_stage_0(current_channel->response_head);

	if (t_58 != NULL)
		return t_58->sensitivity;

	t_48 = find_type_48(current_channel->response_head);

	if (t_48 != NULL)
		return t_48->sensitivity;

	/*else , an error */	

	fprintf(stderr, "WARNING - couldn't find blockette (58/48)!!\n");
 
	fprintf(stderr, "For station: %s; channel: %s\n",
				current_station->station, 
				current_channel->channel);
 
	return -1;      /* error condition */

}

/* ------------------------------------------------------------------------- */