ah_resp.c

解吸SEED格式的源代码

     * normalization at the frequency of the digital sensitivity.
     *
     *
     * if fn is not equal to fs then
     *  A0 = abs(prod{n=1...Np} [2*pi*i*fs - P(n)] /
                        prod{m=1..Nz} [2*pi*i*fs - Z(m)])
     *
     * endif
     * i = sqrt(-1)
     */

	/* default to the A0 as computed above */
	if ((num_zeros == 0) || (num_poles == 0))
		calculated_A0 = A0;
	else
    		calculated_A0 = calc_A0(num_poles, 
					num_zeros, 
					*poles,
					*zeros,
					fs);
 
	if (!cmp_floats(fn, fs))
        	A0 = calculated_A0;
	else
    	/* they are the same, perform consistancy check */
    	{
        	/* check to see if they differ by greater than .5% */
        	if (fabs((A0 - calculated_A0) / calculated_A0) > .005)
        	{
            		fprintf(stderr, "Warning, Normalization given for station: %s, channel %s is :%f.\nThis is inconsistent with the value calculated from poles and zeros: %f.\n\n", current_station->station, current_channel->channel, A0, calculated_A0);

            		/* use the calculated A0 */
            		A0 = calculated_A0;
        	}

	}

	*num_ps = num_poles;
	*num_zs = num_zeros;

	return A0;

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

void determine_gamma(t_34, add_zeros)
struct type34 *t_34;
int *add_zeros;

{
	char *p;

	*add_zeros = 0;

	if (t_34 != NULL)
	{

		if (t_34->description)
		{
			p = t_34->description;

			strupr(t_34->description);
			if (strstr(t_34->description, "VEL") != 0)
            		{
                		*add_zeros = 1;
				return;
            		}
			if (strstr(t_34->description, "ACCEL") != 0) 
            		{
                		*add_zeros = 2; 
				return;
            		} 

			if (strstr(t_34->description, "DISP") != 0) 
            		{
				/* no zeros - just go on */
				return;
			}
		}
		else
		{
			p = t_34->name;
			if (strcmp(t_34->name, "M") == 0)
				return;		/* displacement */
 
			if (strcmp(t_34->name, "M/S") == 0)
			{
				*add_zeros = 1;
				return;		/* Velocity */ 
			} 

			if (strcmp(t_34->name, "M/S**2") == 0)
			{
				*add_zeros = 2;
				return;		/* Acceleration */
			}

		}		
	}		/* if t_34 != NULL */
	
	/* if we got here - flag error ! */
	fprintf(stderr, "WARNING - unknown response type - we only know acceleration, velocity, and displacement.\nFound: %s\n", p);
            
   	fprintf(stderr, "For station: %s; channel: %s\n", current_station->station, current_channel->channel);
   
	fprintf(stderr, "Assuming a gamma of zero!\n");
 
	return;
	
}

/* --------------------------------------------------------------- */
struct type53 *find_type_53(r)
struct response *r;


{
	/*
	 * Type = Transfer function type (Blockette 53)
	 *  should be either A (Laplace transform analog response, in rad/sec)
	 *  or B (Analog response, in Hz)
	 */


	int low_stage = 999999;
	struct type53 *ptr = NULL;

    	while (r)
    	{
		if ((r->type == 'P') && 
			(*(r->ptr.type53->transfer) != 'D') &&
			(r->ptr.type53->stage != 0))
			if (r->ptr.type53->stage < low_stage)
			{
				low_stage = r->ptr.type53->stage;
				ptr = r->ptr.type53;
			}

		r = r->next;

    	}

    return ptr;
}

/* ----------------------------------------------------------------------- */
struct type43 *find_type_43(r)
struct response *r;


{
 
	struct type60 *type_60;
		
	struct type43 *type_43;

	int ix;
 
    	int low_stage = 999999;
    	struct type43 *ptr = NULL;

	while (r)
    	{
        	if (r->type == 'R')
		{
			type_60 = r->ptr.type60; 

			for (ix = 0; ix < type_60->number_stages; ix++)
			{		

				type_43 = get_type_43(type_60->stage + ix);

				if ((type_43 != NULL) &&
				   (type_60->stage[ix].value != 0) &&
				   (type_60->stage[ix].value < low_stage))
				{
					ptr = type_43;
					low_stage = type_60->stage[ix].value;
				}
					
			}
		}
               
        r = r->next;
    }
 
    return ptr;
}
 
/* ------------------------------------------------------------------------ */

struct type43 *get_type_43(stage)
struct type60sub1 *stage;

{
	
	struct type43 *type43_ptr;

	int i;

	for (i = 0; i < stage->number_responses; i++)
	{

		type43_ptr = get_43(stage->response[i].reference, type43_head);	
		
		if (type43_ptr)
			return type43_ptr;
	}

	return NULL;

}

/* ------------------------------------------------------------------------- */
struct type43 *get_43(code, t_43)
int code;
struct type43 *t_43;

{
	while (t_43)
	{
		if (t_43->response_code == code)
			break;

		t_43 = t_43->next;
	}

	return t_43;

}

/* ------------------------------------------------------------------------- */
struct type48 *find_type_48(r)
struct response *r;


{

    struct type60 *type_60;
       
    struct type48 *type_48;

    int ix;

    int low_stage = 999999;
    struct type48 *ptr = NULL;

    while (r)
    {
        if (r->type == 'R')
        {
            type_60 = r->ptr.type60;

            for (ix = 0; ix < type_60->number_stages; ix++)
            {

                type_48 = get_type_48(type_60->stage+ix);

                if (type_48 != NULL)
			if (type_60->stage[ix].value < low_stage)
			{
				ptr = type_48;
				low_stage = type_60->stage[ix].value;
			}
    
            }
        }
        r = r->next;
    }

    return ptr;
}

/* ------------------------------------------------------------------------ */
struct type48 *find_type_48_stage_0(r)
struct response *r;


{

    struct type60 *type_60;
       
    int ix;

    struct type48 *ptr = NULL;

    while (r)
    {
        if (r->type == 'R')
        {
            type_60 = r->ptr.type60;

            for (ix = 0; ix < type_60->number_stages; ix++)
            {

		if (type_60->stage[ix].value == 0)
		{
                	ptr = get_type_48(type_60->stage+ix);
		}
    
            }
        }

        r = r->next;
    }

    return ptr;
}

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


struct type48 *get_type_48(stage)
struct type60sub1 *stage;

{

    struct type48 *type48_ptr;

    int i;

    for (i = 0; i < stage->number_responses; i++)
    {

        type48_ptr = get_48(stage->response[i].reference, type48_head);

        if (type48_ptr)
            return type48_ptr;
    }
 
    return NULL;
 
}

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

struct type48 *get_48(code, t_48)
int code;
struct type48 *t_48;

{
    while (t_48)
    {
        if (t_48->response_code == code)
            break;

        t_48 = t_48->next;
    }

    return t_48;

}

/* ------------------------------------------------------------------------- */
struct type33 *find_type_33(t_33, code)
struct type33 *t_33;
int code;

{
	while (t_33)
	{
		if (t_33->code == code)
			return t_33;

		t_33 = t_33 ->next;
	}

	return NULL;

}
 
/* -------------------------------------------------------------- */
struct type34 *find_type_34(t_34, code)
struct type34 *t_34;
int code;
 
{
    while (t_34)
    {
        if (t_34->code == code)
            return t_34;
 
        t_34 = t_34 ->next;
    }
 
    return NULL;
 
}
 
/* -------------------------------------------------------------- */
/* routines to calculate the response A0 normalization from 
 * pole and zeros as supplied by Sandy Stromme.
 *
 */

/* 
struct complex {
        double x;
        double y;
};
*/ 
/* prototypes */
d_complex c_mult();
 
d_complex c_add();
 
d_complex c_sub();
 
d_complex c_div();
 
double c_abs();
/* -------------------------------------------------------------- */
/* entry routine */
float calc_A0(n_ps, n_zs, poles, zeros, ref_freq)
int n_ps;
int n_zs;
struct type53sub poles[];
struct type53sub zeros[]; 
float ref_freq;
 
{
	int i;

    	d_complex numer, denom, f0, hold;
    	double a0;

	 
    	f0.r = 0;
 
    	f0.i = 2 * PI * ref_freq;

	hold.i = zeros[0].imag;
	hold.r = zeros[0].real;
	denom = c_sub(f0, hold);

    	for (i = 1; i < n_zs; i++)
	{
		hold.i = zeros[i].imag;
    		hold.r = zeros[i].real;

       		denom = c_mult(denom, c_sub(f0, hold));

	}

	hold.i = poles[0].imag;
	hold.r = poles[0].real;

	numer = c_sub(f0, hold);
 
	for (i = 1; i < n_ps; i++)
	{
		hold.i = poles[i].imag; 
    		hold.r = poles[i].real;

        	numer = c_mult(numer, c_sub(f0, hold));

	}
   
    	a0 = c_abs(c_div(numer, denom));
 
    	return a0;
    
}
 
/* ----------------------------------------------------------- */
d_complex c_mult (a, b)
d_complex a, b;
{
 

        d_complex mult;
 
        mult.r = (a.r * b.r) - (a.i * b.i);
        mult.i = (a.r * b.i) + (a.i * b.r);
 
        return mult;
 
}
 
d_complex c_add (a, b)
d_complex a, b;
{
 
        d_complex add;
 
        add.r = a.r + b.r;
        add.i = a.i + b.i;
 
        return add;
 
}
 
d_complex c_sub (a, b)
d_complex a, b;
{
 
        d_complex sub;
 
        sub.r = a.r - b.r;
        sub.i = a.i - b.i;
 
        return sub;
 
}

 
 
d_complex c_div (a, b)
d_complex a, b;
 
{
 
        d_complex div;
 
        b.i = -b.i;
        div = c_mult(a, b);
        div.r /= b.r * b.r + b.i * b.i;
        div.i /= b.r * b.r + b.i * b.i;
 
        return div;
 
}
 
 
double c_abs(a)
d_complex a;
 
{
        return (pow((double)(a.r * a.r) + (a.i * a.i), 0.5));
}
/* --------------------------------------------------------------- */
struct type58 *find_type_58_stage_0(r)
struct response *r;
 
{
    while (r)
    {
        if (r->type == 'S')                 /* 'S' for type 58 */
            if (r->ptr.type58->stage == 0)
                return r->ptr.type58;
    
        r = r->next;
    }
 
    return NULL;
    
}
 
/* -------------------------------------------------------------- */
int cmp_floats(f1, f2)
float f1;
float f2;
 
{
    char f_1[30];
    char f_2[30];
 
    sprintf(f_1, "%6.6f", f1);
    sprintf(f_2, "%6.6f", f2);
 
    return (strcmp(f_1, f_2) == 0);
}