output_sac.c

解吸SEED格式的源代码

/*===========================================================================*/
/* SEED reader     |              output_sac               |    subprocedure */
/*===========================================================================*/
/*
	Name:	output_sac
	Purpose:build and write a seismogram file from the SEED header tables
		and data files; serve as an example for users who wish to write
		their own "output_data" procedure.  This procedure writes SAC
		files.
	Usage:	void output_sac();

	Input:	none (gets its data from globally-available tables and files)
	Output:	none (writes a seismogram file; seismogram files are named by
		beginning time, station, and component; for example,
		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.  Writes message
		to stderr describing the seismogram file being written.
	Externals:data_hdr - a table containing information about this seismogram
		current_station - a pointer to the SEED header tables for the
		station at which this seismogram was recorded
		current_channel - a pointer to the SEED header tables for the
		channel of the station at which this seismogram was recorded
	Warnings:unable to open seismogram file for writing
		failure to properly write the seismic data
	Errors:	none
	Called by:time_span_out
	Calls to:none
	Algorithm:build the seismogram file header by copying values from the
			structure "data_hdr" and the "current_station" and
			"current_channel" tables.  Write out the header and the data.
			This procedure may write either ASCII or binary data, depending
			upon the values of the variables "ascii" and "binary"; as
			implemented here, the programmer sets "ascii" to TRUE or FALSE
			and "binary" is automatically set to the negation of "ascii".
	Notes:	"output_data" is the name of the process called by
		"process_data" to write a seismogram file.  In the case shown
		here, the data are written in SAC format; the SAC format was
		chosen because many users have SAC available and will be able to
		use the output files directly.  If the user desires a different
		output format, this code and these notes should contain enough
		information for that user to construct his or her own 
		"output_data" procedure.

		The structure "data_hdr" contains the following information from
		the SEED fixed-length data header:

	structure element                     description
	====================================  ==========================
	struct input_data_hdr                 fixed data header
	{
		char station[5];                   station name
		char location[2];                  station location
		char channel[3];                   channel name
		char unused[2];                    not used
		struct input_time time;            time (see below)
		unsigned short int nsamples;       number samples
		short int sample_rate;             sample rate factor
		short int sample_rate_multiplier;  sample rate multiplier
		char activity_flags;               activity flags
		char io_flags;                     i/o flags
		char data_quality_flags;           data quality flags
		char number_blockettes;            # blockettes which follow
		long int number_time_corrections;  # .0001s time corrections
		unsigned short int bod;            beginning of data
		unsigned short int bofb;           beginning 1st blkt
	};
	struct input_time                     time from input stream
	{
		unsigned short int year;           year
		unsigned short int day;            Julian day
		char hour;                         hour
		char minute;                       minute
		char second;                       second
		char unused;
		unsigned short int fracsec;        10000ths of second
	};

	These header data are available to the programmer in the usual
	ways.  In the case shown here, a SAC header structure was
	defined in "sac.h" and various elements of the SAC header were
	filled in with data available from the SEED fixed-length data
	header.  See the section of code entitled "Build the SAC header"
	for examples.

	Other data are available from the tables pointed to by
	"current_station" and "current_channel".  These data are 
	described in Halbert et al. under the heading "Station Control
	Header Format".  The table structure is described in the leading
	comments to the main program ("rdseed.c") and in the 
	programmer's technical reference manual; the pointers
	"current_station" and "current_channel" are set by the procedure
	"process_data".  See the section of code within "Build the SAC 
	header" commented with "give values available from SEED tables 
	to SAC optional variables" for examples of the use of these
	station- and channel-specific items.

	Numeric data formats included here are correct for SAC.  The
	user may have to change them for output of data to be used by
	other seismic processing programs.
			
	Problems:	none known
	References:	Halbert, S. E., R. Buland, and C. R. Hutt (1988).  Standard for
			the Exchange of Earthquake Data (SEED), Version V2.0,
			February 25, 1988.  United States Geological Survey,
			Albuquerque Seismological Laboratory, Building 10002,
			Kirtland Air Force Base East, Albuquerque, New Mexico
			87115.  82 pp.
		O'Neill, D. (1987).  IRIS Interim Data Distribution Format
		(SAC ASCII), Version 1.0 (12 November 1987).  Incorporated
		Research Institutions for Seismology, 1616 North Fort Myer
		Drive, Suite 1440, Arlington, Virginia 22209.  11 pp.
		Tull, J. (1987).  SAC User's Manual, Version 10.2, October 7,
		1987.  Lawrence Livermore National Laboratory, L-205,
		Livermore, California 94550.  ??? pp.
	Language:	C, hopefully ANSI standard
	Author:		Dennis O'Neill
	Revisions:	07/15/88  Dennis O'Neill  Initial preliminary release 0.9
		11/21/88  Dennis O'Neill  Production release 1.0
		09/19/89  Dennis O'Neill  corrected output string length
		05/01/92  Allen Nance  corrected CMPINC Header variable for SEED Dip definition
		sach.cmpinc = current_channel->dip + 90.0;
		01/19/01  Doug Neuhauser - reintegrated SAC ascii and binary output into
		this single routine.  Made ascii output more efficient.  Reindented code.
*/

#include <stdio.h>
#include <sys/param.h>

#include "rdseed.h"	/* SEED tables and structures */
#include "sac.h"	/* SAC format structures */


extern int Output_PnZs;
extern int EventDirFlag;
extern int read_summary_flag;

extern struct type48 *find_type_48_stage_0(struct response *);

char *summary_file_get_event();
char *get_event_dir_name();
char *YMDtoYDDD();


struct type53 *find_type_53();
struct type43 *find_type_43();
struct type34 *find_type_34();


float get_A0();
 
char strlst();

int sac_reverse_seismic_data();

void output_sac(data_hdr, offset, asciiflag)
struct data_hdr *data_hdr;
int offset;
int asciiflag;
{
	FILE *outfile;			/* output file pointer */
	char outfile_name[100], channel[10];	/* output file name */
	struct sac sach;			/* SAC header structure */
	int i, j,k;				/* counter */
	char character[8+1];	/* for character transfer */
	int ascii;		/* switch for ASCII output */
	int binary;		/* switch for binary output */
	int reverse_flag;	/* Data reversal flag */
	double dip, azimuth;	/* temp holding for these variables */
	struct type48 *type_48;

	char wrkstr[1024], event_time[64];
	struct time ev_time;
	struct time trace_time;
	struct type34 *type_34;

	struct response *response; int flag;

	char orig_dir[MAXPATHLEN];

	getcwd(orig_dir, MAXPATHLEN);

	chdir(output_dir);

	/* select ASCII or binary output (binary chosen here) */
	ascii = asciiflag;		

	binary = !ascii;

	/* initialize SAC Header to all nulls */
	sach = sac_null;

	for (k=0;k<data_hdr->num_mux_chan;k++)
	{
		if (current_channel == NULL)
		{
			fprintf(stderr, "\tERROR - Unexpected End of Multiplexed Channel List\n");
			chdir(orig_dir);
			return;
		}

		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];
        	}

		dip = current_channel->dip;
		azimuth = current_channel->azimuth;

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


		flag = FALSE;
		for (response = current_channel->response_head; 
		     response != NULL; response = response->next)
		{
			reverse_flag &= 1;

			if (response->type == 'S') { 
				if ((response->ptr.type58)->stage == 0)
				{
					if (((response->ptr.type58)->sensitivity < 0.0) && 
					    (check_reverse & 2))
						reverse_flag |= 2;

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

					sach.scale = (float)type_48->sensitivity;
					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 & 1)
			switch (channel[2])
			{
			case 'Z':
				dip = -90.0;
				break;
			case 'N':
				azimuth = 0.0;
				break;
			case 'E':
				azimuth = 90.0;
				break;
			case 'A':
				dip = -60.0;
				azimuth = 0.0;
				break;
			case 'B':
				dip = 60.0;
				azimuth = 120.0;
				break;
			case 'C':
				dip = -60.0;
				azimuth = 240.0;
				break;
			default:
				break; 
			}

		if (reverse_flag & 2) sach.scale = -sach.scale;

/*                 +=======================================+                 */
/*=================|          Build the SAC header         |=================*/
/*                 +=======================================+                 */


		if (data_hdr->nsamples == 0)
		{
			fprintf(stderr,"\tWARNING (output_sac):  Output Data Contains No Samples\n");
			chdir(orig_dir);
			return;
		}

/* give values to SAC-required variables */
		/* give values to SAC-required variables */
		sach.delta = 1 / data_hdr->sample_rate;
		sach.b = 0.0;
		sach.e = ((double) data_hdr->nsamples - 1) / ((double) data_hdr->sample_rate);
		sach.npts = data_hdr->nsamples;
		sach.iftype = ITIME;			/* known a priori */
		sach.leven = TRUE;			/* known a priori */

		/* give values available from struct data_hdr to SAC optional variables */
		/* record start time */
		sach.nzyear = data_hdr->time.year;
		sach.nzjday = data_hdr->time.day;
		sach.nzhour = data_hdr->time.hour;
		sach.nzmin  = data_hdr->time.minute;
		sach.nzsec  = data_hdr->time.second;
		sach.nzmsec = data_hdr->time.fracsec / 10;	/* convert to msec */

		/* save trace time for later computation */
		trace_time.year    = data_hdr->time.year;
		trace_time.day     = data_hdr->time.day;
		trace_time.hour    = data_hdr->time.hour;
		trace_time.minute  = data_hdr->time.minute;
		trace_time.second  = data_hdr->time.second;
		trace_time.fracsec = data_hdr->time.fracsec;


		/* station and component names */
		sprintf (character, "%-8.8s", data_hdr->station);
		strncpy (sach.kstnm, character, 8);
		sprintf (character, "%-8.8s", channel);
		strncpy (sach.kcmpnm, character, 8);
		sprintf (character, "%-8.8s", data_hdr->network);
		strncpy(sach.knetwk, character, 8);
		sprintf (character, "%-8.8s", data_hdr->location);
		strncpy(sach.khole, character, 8);

		/* give values available from SEED tables to SAC optional variables */
		sach.stla = current_channel->latitude;
		sach.stlo = current_channel->longitude;
		sach.stel = current_channel->elevation;
		sach.stdp = current_channel->local_depth;

        	type_34 = find_type_34(type34_head,
				current_channel->signal_units_code);

		if (type_34 == NULL) {
			fprintf(stderr, 
			"Warning - couldn't find the abbrevation for the signal units code! Signal units code =%d\n", 
			current_channel->signal_units_code);
			fprintf(stderr, 
			"For station: %s; channel: %s\n\n", current_station->station, current_channel->channel);
			fprintf(stderr, "Idep set to unknown\n");
			sach.idep = 5; /* Unknown */
		} else {
			if (strcmp(type_34->name, "M") == 0) {
				sach.idep = 6; /* displacement */
			} else if (strcmp(type_34->name, "M/S") == 0) {
				sach.idep = 7; /* velocity */
			} else if (strcmp(type_34->name, "M/S**2") == 0) {
				sach.idep = 8; /* Acceleration */
			} else { 
				sach.idep = 5; /* Unknown */
			}

		}

		sach.cmpaz = azimuth;
		sach.cmpinc = dip + 90.0; /* SEED dip is 90 degrees off from SAC dip */

		if (read_summary_flag)
		{
			char *p;

			memset(wrkstr, 0, sizeof(wrkstr));

			/* recover current event from summary file */
			strncpy(wrkstr, summary_file_get_event(), 511);

			/* skip the source */

			/* source */
			p = strtok(wrkstr, ",");

			if (p == NULL)
			{
				fprintf(stderr, "Error - output_sac, unable to extract event info for sac header!\n");
				fprintf(stderr, "Event: %s\n", 
					summary_file_get_event());

				break;
			}
		
			/* date/time */
			p = strtok(NULL, ","); 

			if (p == NULL) 
			{        
				fprintf(stderr, "Error - output_sac, unable to extract event info for sac header!\n"); 
				fprintf(stderr, "Event: %s\n", 
					summary_file_get_event()); 
 
				break; 
			}

			/* save ev time for later computation */
			strcpy(event_time, p);

         
			/* latitude */
			p = strtok(NULL, ","); 
 
			if (p == NULL)  
			{         
				fprintf(stderr, "Error - output_sac, unable to extract event info for sac header!\n"); 
				fprintf(stderr, "Event: %s\n", 
                                        summary_file_get_event());
  
				break; 
			}
         
			sach.evla = atof(p); 
   		
			/* longitude */