output_ah.c

解吸SEED格式的源代码

/*===========================================================================*/
/* SEED reader     |              output_ah                |    subprocedure */
/*===========================================================================*/
/*
	Name:		output_ah
	Purpose:	build a seismogram file from the SEED header tables
			and data files and writes output seismogram files.
			This procedure writes output in L-DGO AH-files format
			in binary, which
			include response table (A0, DS, Poles & Zeroes of 
			analog response and coefficients of numerator and 
			denominator of digital FIR filter used.
	Usage:		void output_ah ();
			output_ah (data_hdr, offset);
	Input:		pointer to current data and offset in data buffer and 
			gets its data from globally-available tables and files
	Output:		none (writes a seismogram file; seismogram files are named by
			station, component, beginning time, and a tag; for example,
			ANMO.SPZ.1988.01.23.15.34.08.AH is the filename for a
			seismogram from year 1988, January 23, 15:34:08 UT,
			recorded at station ANMO from component SPZ.  Writes message
			to stderr describing the seismogram file being written.
	Called by:	time_span_out();
	Calls to:	fill_ah_resp(); ahio_routine();
	Author:		Won-Young Kim, Lamont-Doherty Geol. Obs
	Revisions:	03/31/91  Base on Dennis O'Neill's (09/19/89) output_data.
			05/06/92  Allen Nance - added "stype" variable in header
*/

#include "rdseed.h"		/* SEED tables and structures */
#include <stdio.h>
#include <stdlib.h>
#include <sys/param.h>
#include <rpc/rpc.h>

#include <math.h>
#include "ahhead.h"		/* LDGO AH format structures */

char *ddd2mmdd();

static int days_in_month[]= {31,28,31,30,31,30,31,31,30,31,30,31,31};

extern int read_summary_flag;
extern int EventDirFlag;
extern struct type48 *find_type_48_stage_0(struct response *);

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

void output_ah (data_hdr, offset)
struct data_hdr *data_hdr;
int offset;
{
	register struct ah_time *t;		/* AH stuff begins */
	ahhed 	hed;
	XDR	xdr_out;
	void 	month_day();
	int 	ier;
	long 	jday;
	void	fill_ah_resp();			/* AH stuff ends */
	char	*get_type33 (), *p;
	FILE 	*outfile;			/* output file pointer */
	char 	outfile_name[100], channel[10];	/* output file name */
	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 for dip and azimuth */
	struct type48 *type_48;

	char wrkstr[512];
        int y, ddd;

	char dirname[MAXPATHLEN];

	char orig_dir[MAXPATHLEN];

	getcwd(orig_dir, MAXPATHLEN);

	chdir(output_dir);

	if (EventDirFlag)
	{
		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_ah");
				strcpy(dirname, "./");
			}

		}

		chdir(dirname);

	}
		
	/* select ASCII or binary output (binary chosen here) */

	ascii = FALSE;
	binary = !ascii;

	for (k=0;k<data_hdr->num_mux_chan;k++)
	{
		if (data_hdr->nsamples == 0)
		{
			fprintf(stderr,"Output Data Contains No Samples\n");
			chdir(orig_dir);
			return;
		}		

		if (current_channel == NULL)
		{
			fprintf(stderr, "\tERROR - Unexpected End of Multiplexed Channel List\n");
			chdir(orig_dir); 

			return;
		}

		reverse_flag = 0;

		memset(channel,0,sizeof(channel));
                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]);
	
		{
			struct response *response; int flag;

			flag = FALSE;
			for (response = current_channel->response_head; response != NULL; response = response->next)
			{
				if (response->type == 'S') {
					if ((response->ptr.type58)->stage == 0) {
						if ((response->ptr.type58)->sensitivity < 0.0)
							reverse_flag |= 2;
						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) 
							reverse_flag |= 2;
	
						flag = TRUE;
					}
				}
			}
			if (!flag && (check_reverse & 2)) fprintf(stderr, "Warning - Stage Zero Not Found on Gain Reversal check\n");
		}

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

		/* get null ah header */

		get_null_head(&hed);

		/* Fill AH staion instrument response	*/
		/* if not channels == AT, SOH, LOG */

		if ((strcmp(current_channel->channel, "AT")  != 0) &&
		    (strcmp(current_channel->channel, "LOG") != 0) && 
            	    (strcmp(current_channel->channel, "SHO") != 0))
		{
			fill_ah_resp (&hed);	

		}
		if (check_reverse & 1)
		{
			if (reverse_flag & 1)
				fprintf (stderr, "Warning... Azimuth/Dip Reversal found %s.%s, Data will be automatically inverted\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           Header Scale Negated, but Response File will not be modified\n\n", current_station->station, channel);

			

			hed.station.DS = -hed.station.DS;
			}

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

		/* Change on both Dip/Azimuth and Sensitivity reversals */

		if (reverse_flag & 2) 
			hed.station.DS = -hed.station.DS;

		if (reverse_flag & 1)
        		switch (strlst(channel))
        		{
                		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;
        		}
 

		/* place the azimuth and the dip into record 
		 * comment. Units is always DISP (m)
		 */
			
		sprintf(hed.record.rcomment, 
				"Comp azm=%3.1f,inc=%3.1f; Disp (m);",
				azimuth, dip);

	/* AH record information */

		hed.record.delta=1/data_hdr->sample_rate;
		hed.record.ndata= data_hdr->nsamples;
		
		hed.record.abstime.yr=data_hdr->time.year;
		jday= data_hdr->time.day;/* given Julian day get Month and Day */
		month_day(&(hed.record.abstime), (int)jday);
		hed.record.abstime.hr= (short) data_hdr->time.hour;
		hed.record.abstime.mn= (short) data_hdr->time.minute;
		hed.record.abstime.sec= (float) data_hdr->time.second+(float)data_hdr->time.fracsec*0.0001;
		hed.record.type=FLOAT;
		if(reverse_flag) 			
			strncat(hed.record.rcomment, 
				"Signal Reversal found & Data inverted;",
				 COMSIZE - strlen(hed.record.rcomment));

		/* AH station information	*/

		sprintf (character, "%-5.5s\0", data_hdr->station);
		strncpy (hed.station.code, character, 6);
		sprintf (character, "%-5.5s\0", channel);
		strncpy (hed.station.chan, character, 6);
		p = get_type33(current_station->owner_code);
		if (p != NULL) strncpy (hed.station.stype, p, 8); else hed.station.stype[0] = 0;
		hed.station.slat=current_channel->latitude;
		hed.station.slon=current_channel->longitude;
		hed.station.elev=current_channel->elevation;

		/* AH event information	*/
		if (read_summary_flag)
		{