ah_resp.c

解吸SEED格式的源代码

/*===========================================================================*/
/* SEED reader     |            ah_resp                    |     sub program */
/*                 |                                       |                 */
/*===========================================================================*/
/*
    Name:       ah_resp
    Purpose:    Extract responses from Standard Exchange of Earthquake Data
            (SEED) format volume and write output data as L-DGO AH-file format.
            The routine fills response table into AH-file header.
    Usage:      usage is:
                ah_resp (&hed)

                    where &hed is pointer to ahhed


        AH format file header includes instrument response table
        in terms of poles and zeroes, and normalization factor (A0)
        and overall digital sensitivity (DS).

        following extra information is put into AH header:

        extra[0] = 1, for analog response only,
             = 2, analog and digital response,
        extra[1] = Normalization frequency for analog response,
        extra[2] = Frequency where sensitivity (DS) is measured,
        station.cali[i].pole.r = Number of numerator coefficients for FIR filter
                         i = number of poles(analog) + 1
        station.cali[j].pole.r = numerator coefficients are stored
        station.cali[j].pole.i = starting from j, in parallel

    Author:   Chris Laughbon 
            07/27/94

*/

/* -------------------------------------------------------------------- */
/* ----- Includes ----- */
#include "rdseed.h"
#include "ahhead.h"

/* ----- Defines ----- */

/* ----- Variables --- */
static struct type53sub *poles;
static struct type53sub *zeros;

/* ----- Prototypes ---- */ 

float calc_A0();
 
int cmp_floats();

struct type33 *find_type_33();
struct type34 *find_type_34();
struct type53 *find_type_53();
struct type58 *find_type_58_stage_0();

struct type48 *find_type_48();
struct type48 *find_type_48_stage_0();

struct type48 *get_type_48();
struct type48 *get_48();


struct type43 *find_type_43();
struct type43 *get_type_43();
struct type43 *get_43();



float get_A0();

void determine_gamma();
 
/* ------------------------------------------------------------------------- */
void fill_ah_resp(hed)
ahhed *hed;

{
	struct type33 *type_33;
	struct type34 *type_34;

	struct type53 *type_53;
	struct type43 *type_43;

	int i;

	int num_poles, num_zeros;

	float Sd, fn, fs, A0, calculated_A0;

	int gamma;

	poles = zeros = 0;

	A0 = calculated_A0 = 0;
	
	
	/* find out the response type - set the gamma, etc */
	type_34 = find_type_34(type34_head, current_channel->signal_units_code);

	if (type_34 != NULL)
		determine_gamma(type_34, &gamma);
	else
	{
		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");
		gamma = 0;

	}

	type_33 = find_type_33(type33_head, current_channel->instrument_code);

	/* transfer abbreviations */
	if (type_33 != NULL)
	{
		/* COMSIZE is defined in ahhead.h - currently 80 chars */

		/* instrument name */
		strncat(hed->record.rcomment, 
				type_33->abbreviation, COMSIZE);

	}
	else
	{
		sprintf(hed->station.stype, "N/A");
		sprintf(hed->record.rcomment, "Not Found");
	}

        type_53 = (struct type53 *)NULL;
        type_43 = (struct type43 *)NULL;
        
        /* one or the other will be there, if not.. . finished */
 
        type_53 = find_type_53(current_channel->response_head);
 
        type_43 = find_type_43(current_channel->response_head);
 
	if ((type_53 == NULL) && (type_43 == NULL))
	{
		fprintf(stderr, "AH output(): Unable to continue! Unable to calulate A0. No responses.\n");

		return;
	}

	A0 = get_A0(&poles, &zeros, 
			&num_poles, &num_zeros, 
			&Sd, gamma, 
			current_channel->response_head,
			sizeof(hed->station.cal)/sizeof(hed->station.cal[0]));


	if (A0 == -1) 	/* error */
	{
		if (poles != NULL)
			free(poles);

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

		poles = zeros = 0;

		return;

	}

	/* load up the poles and zeros */
	hed->station.cal[0].pole.r = num_poles;
	hed->station.cal[0].zero.r = num_zeros;

	for (i = 0; i < num_poles; i++)
	{
		hed->station.cal[i+1].pole.r = poles[i].real;

		hed->station.cal[i+1].pole.i = poles[i].imag;
	}

	for (i = 0; i < num_zeros; i++)
	{
		hed->station.cal[i+1].zero.r = zeros[i].real;

		hed->station.cal[i+1].zero.i = zeros[i].imag;
	}

	hed->station.DS = Sd;
	hed->station.A0 = A0;

	free(poles);
	free(zeros);

	return;
	
}

/* -------------------------------------------------------------- */
float get_A0(poles, zeros, num_ps, num_zs, 
		sd, gammas, res_head, max_pzs)
struct type53sub **poles;
struct type53sub **zeros;
int *num_ps, *num_zs;
float *sd;
int gammas;
struct response *res_head;
int max_pzs;

{
	int found_43;
	
	struct type58 *type_58;

	struct type48 *type_48;

	/* struct type60 *type_60; */

        struct type43 *type_43;
        struct type53 *type_53;

	int num_poles, num_zeros;
	float fn;
	float fs;
	float calculated_A0, A0;

	struct response *r = res_head;

	int i;
	int j;


	num_zeros = num_poles = 0;

	A0 = 1;

	type_53 = NULL;

	found_43 = FALSE;

	calculated_A0 = 0;

	/* make sure they are null so realloc will malloc 1st time only */

	*poles = NULL;
	*zeros = NULL;
 
	while (r != NULL) {
		if (r->type == 'P') {

			/* found polesnzeros blockette */

			type_53 = r->ptr.type53;

  			if (type_53->stage == 1) 
				fn = type_53->norm_freq;
 
    			A0 *= type_53->ao_norm;

			if (num_poles + type_53->number_poles >= max_pzs)
			{

				fprintf(stderr, "Warning, exceeded maximum number of poles. Clipping at stage =%d\n", type_53->stage);
                                break;
			}

			if (num_zeros + type_53->number_zeroes >= max_pzs) 
                        { 
 
                                fprintf(stderr, "Warning, exceeded maximum number of zeros. Clipping at stage =%d\nNetwork: %s, Station: %s, Channel: %s, Location:%s\n", 
					type_53->stage,
					current_station->network_code,
					current_station->station,
					current_channel->channel,
					current_channel->location);
                                break;
                        } 
 


			/* reallocate memory to hold the pNzs */
			if (type_53->number_poles > 0)
			if ((*poles = (struct type53sub *)
					realloc(*poles, 
						sizeof(struct type53sub) * 
							(type_53->number_poles + num_poles))) == NULL)
			{
				fprintf(stderr, "Unable to obtain memory for AH poles conversion!\nA0 set to zero for station %s; channel %s\n",
					current_station->station,
					current_channel->channel);

				return 0;
	
		
			}

			if ((type_53->stage == 1 ? type_53->number_zeroes + gammas : 
						type_53->number_zeroes) > 0)
			if ((*zeros = (struct type53sub *)
				realloc(*zeros, sizeof(struct type53sub) * 
						((type_53->stage == 1 ? type_53->number_zeroes + gammas : 
								type_53->number_zeroes) + num_zeros))) == NULL)
                	{
                       		fprintf(stderr, "Unable to obtain memory for AH zeros conversion!\nA0 set to zero for station %s; channel %s\n",
                                	current_station->station,
                                	current_channel->channel);
                        	return 0;
                	}

		
                	for (i = num_poles; i < type_53->number_poles + num_poles; i++)
                        	(*poles)[i] = type_53->pole[i - num_poles];

                	for (i = num_zeros; i < type_53->number_zeroes + num_zeros; i++)
                        	(*zeros)[i] = type_53->zero[i - num_zeros];

			/* add zeros */
			if (type_53->stage == 1)
			{
				for (i = 0; i < gammas; i++)
				{
                                	(*zeros)[num_zeros + type_53->number_zeroes + i].real = 0;
					(*zeros)[num_zeros + type_53->number_zeroes + i].imag = 0;
                        		(*zeros)[num_zeros + type_53->number_zeroes + i].real_error = 0;
                        		(*zeros)[num_zeros + type_53->number_zeroes + i].imag_error = 0;
				}

				num_zeros += gammas;
			}


    			/*
     		 	 * First, AH assumes the units of the poles and zeros are rad/sec,
     		 	 * so we convert Type B (Hz) to Type A (rad/sec) if necessary.
     		 	 *  
     		  	 * If Type==B then convert to type A format by:
     		 	 *
     		 	 * P(n) = 2*pi*P(n)      { n=1...Np }
     		 	 * Z(m) = 2*pi*Z(m)      { m=1...Nz }
    		 	 * A0   = A0 * (2*pi)**(Np-Nz)
     		 	 */ 

			if (*type_53->transfer == 'B')
    			{
 
				for (i = 0; i < type_53->number_poles; i++)
        			{
            				(*poles)[num_poles + i].real *= (2 * PI);
            				(*poles)[num_poles + i].imag *= (2 * PI);
        			}
 
        			for (i = 0; i < type_53->number_zeroes; i++)
        			{
            				(*zeros)[num_zeros + i].real *= (2 * PI);
            				(*zeros)[num_zeros + i].imag *= (2 * PI);
        			}
 
        			/* A0   = A0 * (2*pi)**(Np-Nz) */
        			A0 = A0 *
            				pow(2 * PI, 
					    (double)(type_53->number_poles - 
						type_53->number_zeroes));
 
    			}       /* if transfer function was analog - 'B' */

    			num_poles += type_53->number_poles;
    			num_zeros += type_53->number_zeroes;


		/* type 53s */
		} else if (r->type == 'R') {

			int stage;
			int i;

			/*for (i = 0, stage = 1; i < r->ptr.type60->number_stages; i++, stage++)*/
			for (stage = 1; stage <= r->ptr.type60->number_stages;stage++) {

				type_43 = get_type_43(r->ptr.type60->stage+(stage - 1));

				if (type_43 == NULL)
					continue;

				if (stage == 1)
    					fn = type_43->norm_freq;
 
    				A0 *= type_43->ao_norm;

				/* reallocate memory to hold the pNzs */
				if (type_43->number_poles > 0)
				if ((*poles = (struct type53sub *)
					realloc(*poles, 
						sizeof(struct type53sub) * 
							(type_43->number_poles + num_poles))) == NULL)
				{
					fprintf(stderr, "Unable to obtain memory for AH poles conversion!\nA0 set to zero for station %s; channel %s\n",
						current_station->station,
						current_channel->channel);

					return 0;
				}

			if ((stage == 1 ? type_43->number_zeroes + gammas : 
						type_43->number_zeroes) > 0)
				if ((*zeros = (struct type53sub *)
					realloc(*zeros, sizeof(struct type53sub) * 
						(stage == 1 ? type_43->number_zeroes + gammas :
							type_43->number_zeroes + num_zeros))) == NULL)
               			{
               				fprintf(stderr, "Unable to obtain memory for AH zeros conversion!\nA0 set to zero for station %s; channel %s\n",
                                		current_station->station,
                                		current_channel->channel);
                       			return 0;
               			}

		
                		for (i = num_poles; i < type_43->number_poles + num_poles; i++)
				{
					(*poles)[i].real = type_43->pole[i-num_poles].real;
                        		(*poles)[i].imag = type_43->pole[i-num_poles].imag;
                        		(*poles)[i].real_error = type_43->pole[i-num_poles].real_error;
                        		(*poles)[i].imag_error = type_43->pole[i-num_poles].imag_error;
				}

                		for (i = num_zeros; i < type_43->number_zeroes + num_zeros; i++)
				{
                        		(*zeros)[i].real = type_43->zero[i-num_zeros].real;
                        		(*zeros)[i].imag = type_43->zero[i-num_zeros].imag;
                        		(*zeros)[i].real_error = type_43->zero[i-num_zeros].real_error;
                        		(*zeros)[i].imag_error = type_43->zero[i-num_zeros].imag_error;
                		}


				if (stage == 1)
				{
                                	for (i = 0; i < gammas; i++)
                                	{
                                        	(*zeros)[num_zeros + type_43->number_zeroes + i].real = 0;
                                        	(*zeros)[num_zeros + type_43->number_zeroes + i].imag = 0;
                                        	(*zeros)[num_zeros + type_43->number_zeroes + i].real_error = 0;
                                        	(*zeros)[num_zeros + type_43->number_zeroes + i].imag_error = 0;
                                	}

                                	num_zeros += gammas;
				}
    				/*
     		 	 	 * First, AH assumes the units of the poles and zeros are rad/sec,
     		 	 	 * so we convert Type B (Hz) to Type A (rad/sec) if necessary.
     		 	 	 *  
     		  	 	 * If Type==B then convert to type A format by:
     		 	 	 *
     		 	 	 * P(n) = 2*pi*P(n)      { n=1...Np }
     		 	  	 * Z(m) = 2*pi*Z(m)      { m=1...Nz }
    		 	 	 * A0   = A0 * (2*pi)**(Np-Nz)
     		 	 	 */ 

				if (type_43->response_type == 'B')
    				{
 
					for (i = 0; i < type_43->number_poles; i++)
        				{
            					(*poles)[num_poles + i].real *= (2 * PI);
            					(*poles)[num_poles + i].imag *= (2 * PI);
        				}
 
        				for (i = 0; i < type_43->number_zeroes; i++)
        				{
            					(*zeros)[num_zeros + i].real *= (2 * PI);
            					(*zeros)[num_zeros + i].imag *= (2 * PI);
        				}
 
        				/* A0   = A0 * (2*pi)**(Np-Nz) */
        				A0 = A0 *
            					pow(2 * PI, 
						    (double)(type_43->number_poles - 
								type_43->number_zeroes));
 
    				}       /* if transfer function was analog - 'B' */

    				num_poles += type_43->number_poles;
    				num_zeros += type_43->number_zeroes;

				found_43 = TRUE;

			}

		}

		r = r->next;
	}

	if (type_53 == NULL && (!found_43))
	{
		fprintf(stderr, "Warning - couldn't find the Poles and Zeros!!\n");

		fprintf(stderr, 
			"For station: %s; channel: %s\n", 
			current_station->station, 
			current_channel->channel);
        
        	return -1;	/* flag error condition */
	}

 
    /*
     * Second, there is no place to specify the units of the response.
     * An AH file assumes that if an instrument response is deconvolved,
     * the seismogram will be displacement in meters.
     * Convert velocity or acceleration to displacement:
     *
     * Convert to displacement:
     *
     * if acceleration, gamma=2    \
     * elseif velocity,     gamma=1 \
     * elseif displacement, gamma=0  \___Done above
     * else  print error message     /
     * endif                        /
     *
     * Sd = Sd * (2*pi*fs)**gamma
     * Nz = Nz + gamma
     * set values of new zeros equal zero
     * A0 = A0 / (2*pi*fn)**gamma
     * Units = M - Displacement Meters
     */
                
    	type_58 = find_type_58_stage_0(res_head);
 
    	if (type_58 != NULL)
	{
		*sd = type_58->sensitivity;
    		fs = type_58->frequency;
	}
	else
	{
		type_48 = find_type_48_stage_0(res_head);
		
		if (type_48 != NULL)
		{
			*sd = type_48->sensitivity;
    			fs = type_48->frequency;

		}
		else
		{	

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

	} 
 
	*sd *= pow(2 * PI * fs, (double)gammas);
 
	A0 = A0 / pow((double)(2 * PI * fn), (double)gammas);

    /*
     * Third, there is no place in the AH header to specify either
     * the frequency of normalization or the frequency of the
     * digital sensitivity.  This is not a problem as long as these
     * two are the same.  If they are different then evaluate the