output_sac.c

解吸SEED格式的源代码

	chdir(output_dir);

	while (current_station)
	{
		/* check to see if this station was requested by the user */
		if (!chk_station(current_station->station))
		{
			current_station = current_station->next;
			continue;
		}

		if ((type10.version >= 2.3) && 
			!chk_network(current_station->network_code))
		{
			current_station = current_station->next;
            		continue;
		}

		/* check to see if this station lies within the requested time
		 * spans. Convert station time into time structure
		 */
		timecvt(&start, current_station->start);
		timecvt(&end, current_station->end);

		/* for chk_time to work, we need to convert zeros in the
		 * end time struct into really big numbers.
	 	 */

		/* probably should be done in timecvt() */
		cnvt_end_time(&end);
		
		if (!chk_time(start, end))
		{
			current_station = current_station->next;
           		continue;
		}
		
		/* now we have a good station, do the same for its channels */

		current_channel = current_station->type52_head;

		while (current_channel)
		{
			/* check to see if this station was requested by the user */ 
			if (!chk_channel(current_channel->channel)) 
			{   
				current_channel = current_channel->next;
				continue; 
			}    

			/* ignore LOG channels */
			if (strcasecmp(current_channel->channel, "LOG") == 0)
			{
				current_channel = current_channel->next;
                                continue;
			}
 
			/* check to see if this channel lies within the requested time   
			 * spans. Convert channel time into time structure   
			 */ 
			timecvt(&start, current_channel->start);            
			timecvt(&end, current_channel->end);    

			cnvt_end_time(&end);
 
			if (!chk_time(start, end))   
			{
				current_channel = current_channel->next; 
				continue;
			}
        
			/* now we have a good station/channel; write out poles */     
			write_out_PnZs(start);

			current_channel = current_channel->next;

		}		/* while channels */

		current_station = current_station->next;	
	
	}				/* while stations */

	current_station = t_50;	/* there maybe downwind ramifications, so
				 * leave as is 
				 */
	current_channel = t_52;

	chdir(orig_dir);

}
/* ----------------------------------------------------------------------- */
char determine_orient_code();


int write_out_PnZs(struct time start)

{
	
	struct type34 *type_34;

	struct type53sub *poles;
	struct type53sub *zeros;

	float Sd, fn, fs, A0, calculated_A0;

	int num_zeros;
	int num_poles;
	int add_zeros;
	int finished;


	char transfer_type;	/* from SEED manual A laplace, B- analog, etc */	
	FILE *fptr;

	/* file name : SAC_PZS_Network_Station_Channel */
	char fname[200];

	struct response *response;

	char orient;

	poles = zeros = 0; 

	finished = FALSE;

	/* if multiplexed data, check for orientation for this channel */
	if (current_channel->subchannel > 0)
	{
		current_channel->channel[2] = 
			determine_orient_code(current_channel); 

		current_channel->channel[3] = 0;
	}

	sprintf(fname, "SAC_PZs_%s_%s_%s_%s_%04d.%03d.%02d.%02d.%02d.%04d",
                                current_station->network_code ?
                                                current_station->network_code :
                                                "NA",
                                current_station->station,
                                current_channel->channel,
                                current_channel->location,
				start.year,
				start.day,
				start.hour,
				start.minute,
				start.second, 
				start.fracsec);

	fptr = fopen(fname, "a");

	if (fptr == NULL)
	{
		fprintf(stderr, "Unable to open poles and zeros file. Going on.\n\n");

		perror("output_polesNzeros()");

		return -1;

	}

        /* from the type 34 we determine the ground motion type,
         * displacement, velocity, or acceleration
         */
        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, "Setting the number of zeros to add to 0\n");
                add_zeros = 0;
        }
        else
                determine_gamma(type_34, &add_zeros);  /* see ah_resp.c */

	
	A0 = get_A0(&poles, &zeros,
                        &num_poles, &num_zeros,
                        &Sd, add_zeros, 
			current_channel->response_head,
			99999);

	if (A0 == -1) 
	{
		/* error condition */
	       	fclose(fptr);

		if (poles != NULL)
			free(poles);

		if (zeros != NULL)
			free(zeros);


	        return 1;

	}


	print_PnZs(poles,
       			zeros,
              		num_poles,
               		num_zeros,
               		fptr);
			
       	if (fprintf(fptr, "CONSTANT %e\n", Sd * A0) == -1)
       	{
		fprintf(stderr, "\tWARNING (output_polesNzeros():  ");
               	fprintf(stderr, "\tUnable to write to file\n");

               	perror("output_resp()");
 
               	fprintf(stderr, "\tExecution continuing.\n");
 
               	fclose(fptr);

               	return;
 
       	}

	if (poles != NULL) 
		free(poles);

	if (zeros != NULL) 
		free(zeros); 

	fclose(fptr);


	return 1;

} 				/* writeout_PnZs() */

/* ------------------------------------------------------------------------- */
int print_PnZs(ps, zs, n_ps, n_zs, fptr)
struct type53sub *ps;	/* array of poles */
struct type53sub *zs;	/* array of zeros */
int n_ps;		/* number of poles */
int n_zs;		/* number of zeros */
FILE *fptr;

{
	int i;

	/* Print number of zeros */
	fprintf(fptr, "ZEROS %d\n", n_zs);

	for (i = 0; i < n_zs; i++)
		if ((zs[i].real == 0) && (zs[i].imag == 0))
			continue;
		else
                        if (fprintf(fptr, "%4.4f  %4.4f\n",
                                        zs[i].real,
                                        zs[i].imag) == -1)
                        {
                                fprintf(stderr, "\tWARNING (print_PnZs():  ");
                                fprintf(stderr, "\tUnable to write to file\n");
 
                                perror("print_PnZs()");
 
                                fprintf (stderr, "\tExecution continuing.\n");
                        }
 

	/* Print number of poles */ 
        fprintf(fptr, "POLES %d\n", n_ps);    

        for (i = 0; i < n_ps; i++)   
		if ((ps[i].real == 0) && (ps[i].imag == 0)) 
                        continue;
                else 
			if (fprintf(fptr, "%4.4f  %4.4f\n",
						ps[i].real,
						ps[i].imag) == -1)
			{
				fprintf(stderr, "\tWARNING (print_PnZs():  ");
				fprintf(stderr, "\tUnable to write to file\n");

				perror("printPnZs()");

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

			}

	return;

}
 
/* -------------------------------------------------------------- */
/* -------------------------------------------------------------- */

void cnvt_end_time(e)
struct time *e;

{

		if (e->year == 0)
			e->year = 99999;	/* that outta do it */
		if (e->day == 0)
			e->day = 9999;
		if (e->hour == 0)
			e->hour = 24;
		if (e->minute == 0)
			e->minute = 60;
		if (e->second == 0)
			e->second = 60;
		if (e->fracsec == 0)
			e->fracsec = 99999;
}


char *get_event_dir_name()

{

	static char dirname[100];

	char this_event[128], this_date_time[23], wrkstr[40];
	char *p;

	memset(this_event, 0, sizeof(this_event));

	/* we only need the event line up to the data and time */
	strncpy(this_event, summary_file_get_event(), 100);

	/* skip the source */
	p = strtok(this_event, ",");

	/* copy entire date/time field */
	p = strtok(NULL, ",");

	if (p == NULL)
	{
		p = "1900.00.00.00.00.00.0";
	}
		
	strcpy(this_date_time, p);
	
	/* date and time */
       	p = strtok(this_date_time, "/");

	if (p == NULL)
		p = "";

	/* use atoi() to strip spaces */
	sprintf(dirname, "Event_%d.", atoi(p));

	/* month */
	p = strtok(NULL, "/");

	if (p == NULL)
		p = "";

	sprintf(wrkstr, "%02d", atoi(p));
	strcat(dirname, wrkstr);

	/* day */
	p = strtok(NULL, " ");

	if (p == NULL)
		p = "";
        
	sprintf(wrkstr, ".%02d", atoi(p)); 
	strcat(dirname, wrkstr);
 
	/* hour:minute:second.frac */
	p = strtok(NULL, ":");

	if (p == NULL)
		p = "";

	sprintf(wrkstr, ".%02d", atoi(p)); 
	strcat(dirname, wrkstr);

	/* minute */
	p = strtok(NULL, ":");
 
	if (p == NULL)
		p = "";
	
	sprintf(wrkstr, ".%02d", atoi(p)); 
	strcat(dirname, wrkstr);
 
	/* seconds */
	p = strtok(NULL, ".");
 
	if (p == NULL)
		p = "";

	sprintf(wrkstr, ".%02d", atoi(p)); 
	strcat(dirname, wrkstr);

	/* fractional part */
	p = strtok(NULL, ",");
 
	if (p == NULL)
		p = "";
 
	sprintf(wrkstr, ".%02d", atoi(p));
	strcat(dirname, wrkstr);
 

	return dirname;

}		/* get_event_dir_name() */

/* -------------------------------------------------------------------- */
int check_for_reversal(dip, azimuth, c)
float dip; 
float azimuth;
char c;

{

	int reverse_flag = 0;

	switch (c)
	{
		case 'Z':
			if (dip==90.0 && azimuth==0.0)
				reverse_flag = 1;
			break;
		case 'N':
			if (dip==0.0 && azimuth==180.0)
				reverse_flag = 1;
			break;
		case 'E':
			if (dip==0.0 && azimuth==270.0)
				reverse_flag = 1;
			break;
		case 'A':
			if (dip==60.0 && azimuth==180.0)
				reverse_flag = 1;
			break;
		case 'B':
			if (dip==60.0 && azimuth==300.0)
				reverse_flag = 1;
			break;
		case 'C':
			if (dip==60.0 && azimuth==60.0)
				reverse_flag = 1;
			break;
		case 'H' : case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0': case 'O':
			break;
		default:
			fprintf(stderr, "Warning... Invalid Channel Direction character on %s,%s\n", current_station->station, current_channel->channel);
			break; 

	}		/* switch */

	return reverse_flag;

}

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