1 #ifdef HAVE_CONFIG_H
2 #include <config.h>
3 #endif
4
5 /*===================================================================
6 Name: evresp_ Version 3.0
7 Purpose:
8 FORTRAN callable interface to the evresp routine (below)
9 Reference:
10 SEED. Standard for the Exchange of Earthquake Data
11 Reference Manual
12 SEED Format Version 2.3 or later
13 ??? 1995
14 Author: Thomas J. McSweeney
15
16 Usage (from FORTRAN):
17
18 nmatch = evresp(sta,cha,net,datime,units,file,freq,nfreqs,resp,rtype,
19 1 verbose, start_stage, stop_stage)
20
21 Notes:
22 C users should call 'evresp' directly, rather than using this interface.
23 This interface includes extra arguments that are required by the FORTRAN
24 compiler (the length of each string in the argument list, in the order
25 that they appear in the argument list), which C programmers will probably
26 not want to include in their call)
27
28 whereas the C function returns a linked list of responses (one for each
29 response that matched the user's request), this routine returns the
30 response for one (1) station-channel-network for one (1) effective time.
31 If more than one match is found for a given station-channel-network-time,
32 an error condition is raised (and a value of -1 is returned to the calling
33 routine to indicate failure). Likewise, a value of 1 is returned if no
34 match is found for the given station-channel-network-time. If a unique
35 match is found, a value of 0 is returned to the calling routine
36
37 *=================================================================*/
38 /*
39 8/28/2001 -- [ET] Moved several variable definitions from 'evresp.h'
40 into this file.
41 10/19/2005 -- [ET] Added 'evresp_itp()' function with support for
42 List blockette interpolation; made 'evresp()'
43 call 'evresp_itp()' function with default
44 values for List blockette interpolation parameters.
45 2/13/2006 -- [ET] Moved 'use_delay()' function from 'evalresp.c'
46 to 'evresp.c'; modified to close input file
47 when a single response file is specified.
48 3/28/2006 -- [ET] Added "free(freqs_orig)" at end of 'evresp_itp()'
49 function; added "free_matched_files(output_files)"
50 in 'evresp_itp()' function.
51 8/21/2006 -- [IGD] Version 3.2.36: Added support for TESLA units.
52 8/23/2006 -- [ET] Improved detection of pressure ("PA") input units
53 for in/out units check.
54 10/16/2006 -- [ET] Modified to free array allocated in 'evresp_itp()'
55 function.
56 02/27/2007 -- [IGD] Added return (#ifdef LIB_MODE) if the input file is not
57 found
58 */
59
60 #include "./evresp.h"
61 #include <stdlib.h>
62 #include <string.h>
63
64 /* define a global flag to use if using "default" units */
65 int def_units_flag;
66
67 /* define a pointer to a channel structure to use in determining the input and
68 output units if using "default" units and for use in error output*/
69 struct channel *GblChanPtr;
70 float unitScaleFact = 1.0;
71
72 /* define global variables for use in printing error messages */
73 char *curr_file;
74 int curr_seq_no;
75
76 /* and set a global variable to contain the environment for the setjmp/longjmp
77 combination for error handling */
78 jmp_buf jump_buffer;
79
80 char myLabel[20];
81
evresp_(char * sta,char * cha,char * net,char * locid,char * datime,char * units,char * file,float * freqs,int * nfreqs_in,float * resp,char * rtype,char * verbose,int * start_stage,int * stop_stage,int * stdio_flag,int lsta,int lcha,int lnet,int llocid,int ldatime,int lunits,int lfile,int lrtype,int lverbose,int useTotalSensitivityFlag)82 int evresp_(char *sta, char *cha, char *net, char *locid, char *datime,
83 char *units, char *file, float *freqs, int *nfreqs_in, float *resp,
84 char *rtype, char *verbose, int *start_stage, int *stop_stage,
85 int *stdio_flag, int lsta, int lcha, int lnet, int llocid,
86 int ldatime, int lunits, int lfile, int lrtype, int lverbose, int useTotalSensitivityFlag)
87 {
88 struct response *first = (struct response *)NULL;
89 double *dfreqs;
90 int i,j, nfreqs, start, stop, flag;
91
92 /* add null characters to end of input string arguments (remove trailing
93 spaces first */
94
95 add_null(sta, lsta-1, 'a');
96 add_null(cha, lcha-1, 'a');
97 add_null(net, lnet-1, 'a');
98 add_null(locid, llocid-1, 'a');
99 add_null(datime, ldatime-1, 'a');
100 add_null(units, lunits-1, 'a');
101 add_null(file, lfile-1, 'a');
102 add_null(rtype, lrtype-1, 'a');
103 add_null(verbose, lverbose-1, 'a');
104
105 nfreqs = *nfreqs_in;
106 start = *start_stage;
107 stop = *stop_stage;
108 flag = *stdio_flag;
109
110 dfreqs = alloc_double(nfreqs);
111 for(i = 0; i < nfreqs; i++)
112 dfreqs[i] = freqs[i];
113
114 /* then call evresp */
115
116 first = evresp(sta, cha, net, locid, datime, units, file, dfreqs, nfreqs,
117 rtype, verbose, start, stop, flag, useTotalSensitivityFlag);
118
119 /* free up the frequency vector */
120
121 free(dfreqs);
122
123 /* check the output. If no response found, return 1, else if more than one response
124 found, return -1 */
125
126 if(first == (struct response *)NULL) {
127 return(1);
128 }
129 else if(first->next != (struct response *)NULL) {
130 free_response(first);
131 return(-1);
132 }
133
134 /* if only one response found, convert from complex output vector into multiplexed
135 real output for FORTRAN (real1, imag1, real2, imag2, ..., realN, imagN) */
136
137 for(i = 0, j = 0; i < nfreqs; i++) {
138 resp[j++] = (float) first->rvec[i].real;
139 resp[j++] = (float) first->rvec[i].imag;
140 }
141
142 /* free up dynamically allocated space */
143
144 free_response(first);
145
146 /* and return to FORTRAN program */
147
148 return(0);
149
150 }
151
152 /* IGD 03/01/05 Small function to set and return
153 * a static flag to use or not use the estimated delay in
154 * response computation
155 * Input: NEGATIVE means that we want to query the value of the flag
156 * TRUE or FALSE means that we want to set corresponding values
157 * The reason we want to use this global variable is because we don't
158 * want to change the number of arguments in evresp() function which
159 * is used in users programs
160 */
161 /* 2/6/2006 -- [ET] Moved from 'evalresp.c' to 'evresp.c' */
use_delay(int flag)162 int use_delay(int flag)
163 {
164 /* WE USE THOSE WEIRD magic numbers here because
165 * there is a chance that use_delay_flag is not
166 * defined: in user program which uses evresp()
167 * when use_delay() is not used before evresp().
168 */
169 int magic_use_delay = 35443647;
170 int magic_dont_use_delay = -90934324;
171 static int use_delay_flag = FALSE;
172 if (TRUE == flag)
173 use_delay_flag = magic_use_delay;
174 if (FALSE == flag)
175 use_delay_flag = magic_dont_use_delay;
176
177 if (use_delay_flag == magic_use_delay)
178 return TRUE;
179 return FALSE;
180 }
181
182
183
184
185 /*===================================================================
186 Name: evresp Version 3.0
187 Purpose:
188 Extract channel response parameters from either ASCII
189 files produced by rdseed -r ("response" file) or
190 rdseed -d ("sta-cha" files) and calculate the complex
191 response.
192 Reference:
193 SEED. Standard for the Exchange of Earthquake Data
194 Reference Manual
195 SEED Format Version 2.3 or later
196 ??? 1995
197 Author: Thomas J. McSweeney
198 Modofications: Ilya Dricker (i.dricker@isti.com) IGD for versions of evalresp 3.2.17
199 Notes:
200 ???. Version 3.0
201 - modified to parse "new" rdseed RESP file output (includes a new
202 field that contains the blockette and field numbers for each of
203 the items in the RESP file)
204 - is a very substantial change over the previous releases of
205 evresp. The code has been completely rewritten from the original
206 form authored by Jean-Francios Fels to support several new features.
207 among them are:
208 (a) a "new" RESP file format that contains the blockette and
209 field numbers as prefixes to each line. This allows for
210 quick determination of whether or not the program is
211 parsing the correct information without relying on searching
212 for non-standardized character strings in the RESP file
213 (b) support for the blockette [61] responses
214 (c) support for the response-reference style responses (i.e.
215 a blockette [60] followed by a series of blockette [41] or
216 blockette [43] through blockette [48] responses)
217 - the code has been rewritten so that the calculations are all confined
218 to this function and the functions that it calls. All the user
219 has to do us supply the appropriate control parameters to this function
220 - the parsing has been entirely reworked so that each blockette style is
221 parsed in a seperate. This should make the code easier to maintain and
222 allow for changes in the output from RDSEED (either in number of fields
223 on a line or in which fields are output from a given blockette)
224 - the code has been converted to ANSI standard C, rather than K&R style C
225
226 Thomas J. McSweeney: tjm@iris.washington.edu
227
228 *=================================================================*/
229
230 double Pi = 3.141592653589793;
231 double twoPi = 6.283185307179586;
232
233 /* IGD 08/21/06 Added Tesla */
234 char SEEDUNITS[][UNITS_STR_LEN] = {"Undef Units", "Displacement", "Velocity",
235 "Acceleration", "Counts", "Volts", "", "Pascals", "Tesla"};
236
237 char FirstLine[MAXLINELEN];
238 int FirstField;
239
240 /* This version of the function includes
241 the 'listinterp...' parameters */
242
evresp_itp(char * stalst,char * chalst,char * net_code,char * locidlst,char * date_time,char * units,char * file,double * freqs,int nfreqs,char * rtype,char * verbose,int start_stage,int stop_stage,int stdio_flag,int listinterp_out_flag,int listinterp_in_flag,double listinterp_tension,int useTotalSensitivityFlag)243 struct response *evresp_itp(char *stalst, char *chalst, char *net_code,
244 char *locidlst, char *date_time, char *units,
245 char *file, double *freqs, int nfreqs,
246 char *rtype, char *verbose, int start_stage,
247 int stop_stage, int stdio_flag,
248 int listinterp_out_flag, int listinterp_in_flag,
249 double listinterp_tension, int useTotalSensitivityFlag)
250 {
251 struct channel this_channel;
252 struct scn *scn;
253 struct string_array *sta_list, *chan_list;
254 struct string_array *locid_list;
255 int i, j, k, count = 0, which_matched, mode, new_file;
256 volatile int test = 1;
257 int err_type;
258 char out_name[MAXLINELEN], locid[LOCIDLEN+1];
259 char *locid_ptr, *end_locid_ptr;
260 size_t locid_len;
261 struct matched_files *flst_head = (struct matched_files *)NULL;
262 struct matched_files *flst_ptr = NULL, *output_files = NULL;
263 struct file_list *lst_ptr = NULL, *tmp_ptr = NULL, *out_file = NULL, *tmp_file = NULL;
264 struct response *resp = NULL, *next_ptr = NULL;
265 struct response *prev_ptr = (struct response *)NULL;
266 struct response *first_resp = (struct response *)NULL;
267 struct complex *output = NULL;
268 struct scn_list *scns = NULL;
269 FILE *fptr = NULL;
270 double *freqs_orig = NULL; /* for saving the original frequencies */
271 int nfreqs_orig;
272
273 /* Let's save the original frequencies requested by a user since they can be overwritten */
274 /* if we process blockette 55 IGD for version 3.2.17 of evalresp*/
275
276 nfreqs_orig = nfreqs;
277 freqs_orig = (double *) malloc(sizeof(double) * nfreqs_orig);
278 memcpy (freqs_orig, freqs, sizeof(double) * nfreqs_orig);
279
280 /* set 'GblChanPtr' to point to 'this_channel' */
281
282 GblChanPtr = &this_channel;
283
284 /* clear out the FirstLine buffer */
285
286 memset(FirstLine, 0, sizeof(FirstLine));
287
288
289 /* if the verbose flag is set, then print some diagnostic output (other than
290 errors) */
291
292 if(verbose && !strcmp(verbose,"-v")) {
293 fprintf(stderr, "<< EVALRESP RESPONSE OUTPUT V%s >>\n", REVNUM);
294 fflush(stderr);
295 }
296
297 /* first, determine the values of Pi and twoPi for use in evaluating
298 the instrument responses later */
299
300 Pi = acos(-1.0);
301 twoPi = 2.0 * Pi;
302
303 /* set the values of first_units and last_units to null strings */
304
305 strncpy(this_channel.staname,"",STALEN);
306 strncpy(this_channel.network,"",NETLEN);
307 strncpy(this_channel.locid,"",LOCIDLEN);
308 strncpy(this_channel.chaname,"",CHALEN);
309 strncpy(this_channel.beg_t,"",DATIMLEN);
310 strncpy(this_channel.end_t,"",DATIMLEN);
311 strncpy(this_channel.first_units,"",MAXLINELEN);
312 strncpy(this_channel.last_units,"",MAXLINELEN);
313
314 /* and initialize the linked list of pointers to filters */
315
316 this_channel.first_stage = (struct stage *)NULL;
317
318 /* parse the "stalst" string to form a list of stations */
319
320 for(i = 0; i < (int)strlen(stalst); i++) {
321 if(stalst[i] == ',')
322 stalst[i] = ' ';
323 }
324 sta_list = ev_parse_line(stalst);
325
326 /* remove any blank spaces from the beginning and end of the string */
327
328 locid_ptr = locidlst;
329 strncpy(locid, "", LOCIDLEN);
330 while(*locid_ptr && *locid_ptr == ' ') locid_ptr++;
331 end_locid_ptr = locid_ptr + strlen(locid_ptr) - 1;
332 while(end_locid_ptr > locid_ptr && *end_locid_ptr == ' ') end_locid_ptr--;
333 locid_len = end_locid_ptr - locid_ptr + 1;
334 if (locid_len > LOCIDLEN)
335 locid_len = LOCIDLEN;
336 strncpy(locid, locid_ptr, locid_len);
337 locid[LOCIDLEN] = '\0';
338
339 /* parse the "locidlst" string to form a list of channels */
340
341 locid_list = parse_delim_line(locid,",");
342
343 /* parse the "chalst" string to form a list of channels */
344
345 for(i = 0; i < (int)strlen(chalst); i++) {
346 if(chalst[i] == ',')
347 chalst[i] = ' ';
348 }
349 chan_list = ev_parse_line(chalst);
350
351 /* then form a set of network-station-locid-channel tuples to search for */
352
353 scns = alloc_scn_list(chan_list->nstrings*sta_list->nstrings*locid_list->nstrings);
354 for(i = 0; i < sta_list->nstrings; i++) {
355 for(j = 0; j < locid_list->nstrings; j++) {
356 for(k = 0; k < chan_list->nstrings; k++, count++) {
357 scn = scns->scn_vec[count];
358 strncpy(scn->station, sta_list->strings[i], STALEN);
359 if(strlen(locid_list->strings[j]) == strspn(locid_list->strings[j], " "))
360 memset(scn->locid, 0, LOCIDLEN);
361 else
362 strncpy(scn->locid, locid_list->strings[j], LOCIDLEN);
363 strncpy(scn->channel, chan_list->strings[k], CHALEN);
364 strncpy(scn->network, net_code, NETLEN);
365 }
366 }
367 }
368 #ifdef LOG_LABEL
369 sprintf(myLabel, "[%s.%s.%s.%s]", scn->network, scn->station, scn->locid, scn->channel);
370 #else
371 myLabel[0] = '\0';
372 #endif
373 /* if input is from stdin, set fptr to stdin, else find whatever matching
374 files there are */
375
376 if(stdio_flag) {
377 fptr = stdin;
378 mode = 0;
379 }
380 else {
381 flst_head = find_files(file, scns, &mode);
382 flst_ptr = flst_head;
383 }
384
385 /* find the responses for each of the station channel pairs as they
386 occur in the file */
387
388 if(!mode && !stdio_flag) {
389 curr_file = file;
390 if((fptr = fopen(file,"r")) == (FILE *)NULL) {
391 #ifdef LIB_MODE
392 fprintf(stderr, "%s failed to open file %s\n", myLabel, file);
393 return NULL;
394 #else
395 error_exit(OPEN_FILE_ERROR,"failed to open file %s", file);
396 #endif
397 }
398 }
399
400 /* allocate space for the first response */
401
402 resp = alloc_response(nfreqs);
403
404 for(i = 0; i < scns->nscn && (mode || test); i++) {
405
406 /* allocate space for 'matched_files' pointer used to determine if a
407 file has already been read */
408
409 if(!stdio_flag)
410 output_files = alloc_matched_files();
411
412 /* then check the mode to determine if are parsing one file or a list
413 of files (note: if input is from stdin, is one file) */
414
415 if(!mode) {
416 which_matched = 0;
417 while(test && which_matched >= 0) {
418 new_file = 0;
419 if(!(err_type = setjmp(jump_buffer))) {
420 which_matched = find_resp(fptr, scns, date_time, &this_channel);
421 #ifdef LIB_MODE /* IGD 25-Sep-2007 Looks like we do not need this: function returns anyway */
422 // if(which_matched < 0) {
423 // if(!stdio_flag) /* if not input from console then */
424 // fclose(fptr); /* close input file
425 // return NULL;
426 // }
427 #endif
428
429 /* found a station-channel-network that matched. First construct
430 an output filename and compare to other output files. If this
431 filename doesn't match any of them, (or if it is the first
432 file found) parse the channel's response information.
433 Otherwise skip it (since a match has already been found) */
434
435 sprintf(out_name,"%s.%s.%s.%s",this_channel.network,
436 this_channel.staname,this_channel.locid,
437 this_channel.chaname);
438 #ifdef LOG_LABEL
439 sprintf(myLabel, "[%s]", out_name);
440 #else
441 myLabel[0] = '\0';
442 #endif
443 if(!stdio_flag) {
444 tmp_file = output_files->first_list;
445 for(k = 0; k < output_files->nfiles; k++) {
446 out_file = tmp_file;
447 if(!strcmp(out_file->name,out_name))
448 break;
449 tmp_file = out_file->next_file;
450 }
451 }
452 if((stdio_flag && !new_file) || !output_files->nfiles) {
453 if(!stdio_flag) {
454 output_files->nfiles++;
455 out_file = alloc_file_list();
456 output_files->first_list = out_file;
457 out_file->name = alloc_char(strlen(out_name)+1);
458 strcpy(out_file->name,out_name);
459 }
460 new_file = 1;
461 }
462 else if((stdio_flag && !new_file) || k == output_files->nfiles) {
463 if(!stdio_flag) {
464 output_files->nfiles++;
465 out_file->next_file = alloc_file_list();
466 tmp_file = out_file->next_file;
467 out_file = tmp_file;
468 out_file->name = alloc_char(strlen(out_name)+1);
469 strcpy(out_file->name,out_name);
470 }
471 new_file = 1;
472 }
473 else
474 new_file = 0;
475
476 if(new_file && which_matched >= 0) {
477
478 /* fill in station-channel-net information for the response */
479
480 strncpy(resp->station,this_channel.staname,STALEN);
481 strncpy(resp->locid,this_channel.locid,LOCIDLEN);
482 strncpy(resp->channel,this_channel.chaname,CHALEN);
483 strncpy(resp->network,this_channel.network,NETLEN);
484 output = resp->rvec;
485
486 /* found a station channel pair that matched a response, so parse
487 the response into a channel/filter list */
488
489 test = parse_channel(fptr, &this_channel);
490
491 if(listinterp_in_flag &&
492 this_channel.first_stage->first_blkt->type == LIST)
493 { /* flag set for interpolation and stage type is "List" */
494 interpolate_list_blockette(
495 &(this_channel.first_stage->first_blkt->blkt_info.list.freq),
496 &(this_channel.first_stage->first_blkt->blkt_info.list.amp),
497 &(this_channel.first_stage->first_blkt->blkt_info.list.phase),
498 &(this_channel.first_stage->first_blkt->blkt_info.list.nresp),
499 freqs,nfreqs,listinterp_tension);
500 }
501
502 /* check the filter sequence that was just read */
503 check_channel(&this_channel);
504
505 /* If we process blockette 55, we should recompute resp->rvec */
506 /* because the number of output responses is generally different from */
507 /* what is the user requested */
508 /*if we don't use blockette 55, we should set the frequencies to the original */
509 /* user defined position if we did mess up with frequencies in -possible - blockette 55*/
510 /* containing previous file. Modifications by I.Dricker IGD*/
511
512
513 free(resp->rvec);
514 /* 'freqs' array is passed in and should not be freed -- 10/18/2005 -- [ET] */
515 /* free(freqs); */
516
517 if (this_channel.first_stage->first_blkt != NULL && this_channel.first_stage->first_blkt->type == LIST)
518 {
519 /*to prevent segmentation in case of bogus input files */
520 nfreqs = this_channel.first_stage->first_blkt->blkt_info.list.nresp;
521 freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
522 memcpy (freqs, this_channel.first_stage->first_blkt->blkt_info.list.freq, sizeof(double) * nfreqs); /*cp*/
523 resp->rvec = alloc_complex(nfreqs);
524 output=resp->rvec;
525 resp->nfreqs = nfreqs;
526 resp->freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
527 memcpy (resp->freqs, this_channel.first_stage->first_blkt->blkt_info.list.freq, sizeof(double) * nfreqs); /*cp*/
528 }
529 else
530 {
531 nfreqs = nfreqs_orig;
532 freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
533 memcpy (freqs, freqs_orig, sizeof(double) * nfreqs); /*cp*/
534 resp->rvec = alloc_complex(nfreqs);
535 output=resp->rvec;
536 resp->nfreqs = nfreqs;
537 resp->freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
538 memcpy (resp->freqs, freqs_orig, sizeof(double) * nfreqs); /*cp*/
539 }
540
541
542 /* normalize the response of the filter sequence */
543
544 norm_resp(&this_channel, start_stage, stop_stage, 0);
545
546 /* calculate the response at the requested frequencies */
547
548 calc_resp(&this_channel, freqs, nfreqs, output, units,
549 start_stage, stop_stage, useTotalSensitivityFlag);
550
551 /* diagnostic output, if the user requested it */
552
553 if(verbose && !strcmp(verbose,"-v")) {
554 print_chan(&this_channel, start_stage, stop_stage,
555 stdio_flag, listinterp_out_flag, listinterp_in_flag, useTotalSensitivityFlag);
556 }
557
558 free(freqs); /* free array that was allocated above */
559
560 /* and, finally, free the memory associated with this channel/filter
561 list and continue searching for the next match */
562
563 free_channel(&this_channel);
564 if(first_resp == (struct response *)NULL) {
565 first_resp = resp;
566 }
567 next_ptr = alloc_response(nfreqs);
568 resp->next = next_ptr;
569 prev_ptr = resp;
570 resp = next_ptr;
571 }
572 else {
573 strncpy(FirstLine,"",MAXLINELEN);
574 test = next_resp(fptr);
575 }
576 }
577 else {
578 if(new_file)
579 output_files->nfiles--;
580 free_channel(&this_channel);
581 /* catch errors that cause parsing to fail midstream */
582 if(err_type == PARSE_ERROR || err_type == UNRECOG_FILTYPE ||
583 err_type == UNDEF_SEPSTR || err_type == IMPROP_DATA_TYPE ||
584 err_type == RE_COMP_FAILED || err_type == UNRECOG_UNITS) {
585 strncpy(FirstLine,"",MAXLINELEN);
586 test = next_resp(fptr);
587 }
588 else if(err_type == UNDEF_PREFIX) {
589 test = 0;
590 }
591 }
592 }
593 if(!stdio_flag)
594 free_matched_files(output_files); /* added 3/28/2006 -- [ET] */
595
596 /* allocated one too many responses */
597
598 free_response(resp);
599 if(prev_ptr != (struct response *)NULL)
600 prev_ptr->next = (struct response *)NULL;
601 break;
602
603 }
604 else if(mode) {
605 lst_ptr = flst_ptr->first_list;
606 scn = scns->scn_vec[i];
607 next_scn:
608 for(j = 0; j < flst_ptr->nfiles; j++) {
609 if(!stdio_flag) {
610 fptr = fopen(lst_ptr->name,"r");
611 }
612 if(fptr != (FILE *)NULL) {
613 curr_file = lst_ptr->name;
614 look_again:
615 new_file = 0;
616 if(!(err_type = setjmp(jump_buffer))) {
617 which_matched = get_resp(fptr, scn, date_time, &this_channel);
618 #ifdef LIB_MODE /* IGD 25-Sep-2007 Looks like we do not need this: function returns anyway */
619 // if(which_matched < 1) {
620 // if(!stdio_flag) /* if not input from console then */
621 // fclose(fptr); /* close input file */
622 // return NULL;
623 // }
624 #endif
625 if(which_matched >= 0) {
626
627 /* found a station-channel-network that matched. First construct
628 an output filename and compare to other output files. If this
629 filename doesn't match any of them, (or if it is the first
630 file found) parse the channel's response information.
631 Otherwise skip it (since a match has already been found) */
632
633 sprintf(out_name,"%s.%s.%s.%s",this_channel.network,
634 this_channel.staname,this_channel.locid,
635 this_channel.chaname);
636 #ifdef LOG_LABEL
637 sprintf (myLabel, "[%s]", out_name);
638 #else
639 myLabel[0] = '\0';
640 #endif
641 tmp_file = output_files->first_list;
642 for(k = 0; k < output_files->nfiles; k++) {
643 out_file = tmp_file;
644 if(!strcmp(out_file->name,out_name))
645 break;
646 tmp_file = out_file->next_file;
647 }
648 if(!output_files->nfiles) {
649 output_files->nfiles++;
650 out_file = alloc_file_list();
651 output_files->first_list = out_file;
652 out_file->name = alloc_char(strlen(out_name)+1);
653 strcpy(out_file->name,out_name);
654 new_file = 1;
655 }
656 else if(k == output_files->nfiles) {
657 output_files->nfiles++;
658 out_file->next_file = alloc_file_list();
659 tmp_file = out_file->next_file;
660 out_file = tmp_file;
661 out_file->name = alloc_char(strlen(out_name)+1);
662 strcpy(out_file->name,out_name);
663 new_file = 1;
664 }
665 else
666 new_file = 0;
667
668 if(new_file) {
669
670 /* fill in station-channel-net information for the response */
671
672 strncpy(resp->station,this_channel.staname,STALEN);
673 strncpy(resp->locid,this_channel.locid,LOCIDLEN);
674 strncpy(resp->channel,this_channel.chaname,CHALEN);
675 strncpy(resp->network,this_channel.network,NETLEN);
676 output = resp->rvec;
677
678 /* parse the response into a channel/filter list */
679
680 test = parse_channel(fptr, &this_channel);
681
682 /* IGD 01/04/01 Add code preventing a user from defining output units as DIS and ACC if
683 the input units are PRESSURE after */
684 if (strncmp (this_channel.first_units, "PA -",4) == 0) {
685 if (strcmp(units, "VEL") != 0) {
686 if(strcmp(units, "DEF") != 0) {
687 fprintf(stderr, "%s WARNING: OUTPUT %s does not make sense if INPUT is PRESSURE\n",
688 myLabel, units);
689 strcpy (units, "VEL");
690 fprintf(stderr, "%s OUTPUT units are reset and interpreted as PRESSURE\n", myLabel);
691 }
692 }
693 }
694 /* IGD 08/21/06 Add code preventing a user from defining output units as DIS and ACC if
695 the input units are TESLA */
696 if (strncmp (this_channel.first_units, "T -", 3) == 0) {
697 if (strcmp(units, "VEL") != 0) {
698 if(strcmp(units, "DEF") != 0) {
699 fprintf(stderr, "%s WARNING: OUTPUT %s does not make sense if INPUT is MAGNETIC FLUX\n",
700 myLabel, units);
701 strcpy (units, "VEL");
702 fprintf(stderr, "%s OUTPUT units are reset and interpreted as TESLA\n", myLabel);
703 }
704 }
705 }
706
707 if(listinterp_in_flag &&
708 this_channel.first_stage->first_blkt->type == LIST) {
709 /* flag set for interpolation and stage type is "List" */
710 interpolate_list_blockette(
711 &(this_channel.first_stage->first_blkt->blkt_info.list.freq),
712 &(this_channel.first_stage->first_blkt->blkt_info.list.amp),
713 &(this_channel.first_stage->first_blkt->blkt_info.list.phase),
714 &(this_channel.first_stage->first_blkt->blkt_info.list.nresp),
715 freqs,nfreqs,listinterp_tension);
716 }
717
718 /* check the filter sequence that was just read */
719 check_channel(&this_channel);
720
721 /* If we process blockette 55, we should recompute resp->rvec */
722 /* because the number of output responses is generally different from */
723 /* what is the user requested */
724 /*if we don't use blockette 55, we should set the frequencies to the original */
725 /* user defined position if we did mess up with frequencies in -possible - blockette 55*/
726 /* containing previous file. Modifications by I.Dricker / IGD */
727
728 free(resp->rvec);
729 /* 'freqs' array is passed in and should not be freed -- 10/18/2005 -- [ET] */
730 /* free(freqs); */
731 if (this_channel.first_stage->first_blkt != NULL && this_channel.first_stage->first_blkt->type == LIST) {
732 /* This is to prevent segmentation if the response input is bogus responses */
733 nfreqs = this_channel.first_stage->first_blkt->blkt_info.list.nresp;
734 freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
735 memcpy (freqs, this_channel.first_stage->first_blkt->blkt_info.list.freq, sizeof(double) * nfreqs); /*cp*/
736 resp->rvec = alloc_complex(nfreqs);
737 output=resp->rvec;
738 resp->nfreqs = nfreqs;
739 resp->freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
740 memcpy (resp->freqs, this_channel.first_stage->first_blkt->blkt_info.list.freq, sizeof(double) * nfreqs); /*cp*/
741 }
742 else {
743 nfreqs = nfreqs_orig;
744 freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
745 memcpy (freqs, freqs_orig, sizeof(double) * nfreqs); /*cp*/
746 resp->rvec = alloc_complex(nfreqs);
747 output=resp->rvec;
748 resp->nfreqs = nfreqs;
749 resp->freqs = (double *) malloc(sizeof(double) * nfreqs); /* malloc a new vector */
750 memcpy (resp->freqs, freqs_orig, sizeof(double) * nfreqs); /*cp*/
751 }
752
753 /* normalize the response of the filter sequence */
754
755 norm_resp(&this_channel, start_stage, stop_stage, 0);
756
757 /* calculate the response at the requested frequencies */
758
759 calc_resp(&this_channel, freqs, nfreqs, output, units,
760 start_stage, stop_stage, useTotalSensitivityFlag);
761
762 /* diagnostic output, if the user requested it */
763
764 if(verbose && !strcmp(verbose,"-v")) {
765 print_chan(&this_channel, start_stage, stop_stage,
766 stdio_flag, listinterp_out_flag, listinterp_in_flag, useTotalSensitivityFlag);
767 }
768
769 free(freqs); /* free array that was allocated above */
770
771 /* and, finally, free the memory associated with this
772 channel/filter list and continue searching for the
773 next match */
774
775 free_channel(&this_channel);
776 if(first_resp == (struct response *)NULL) {
777 first_resp = resp;
778 }
779 next_ptr = alloc_response(nfreqs);
780 resp->next = next_ptr;
781 prev_ptr = resp;
782 resp = next_ptr;
783 }
784 FirstField = 0;
785 strncpy(FirstLine,"",MAXLINELEN);
786 if(!stdio_flag) {
787 fclose(fptr);
788 }
789 }
790 else {
791 strncpy(FirstLine,"",MAXLINELEN);
792 test = next_resp(fptr);
793 if (!test) {
794 if(!stdio_flag) {
795 fclose(fptr);
796 }
797 }
798 }
799
800 /* if not the last file in the list, move on to the next one */
801
802 if(lst_ptr->next_file != (struct file_list *)NULL) {
803 tmp_ptr = lst_ptr->next_file;
804 lst_ptr = tmp_ptr;
805 }
806 }
807 else {
808 if (new_file)
809 output_files->nfiles--;
810 /* catch errors that cause parsing to fail midstream */
811 if(err_type == PARSE_ERROR || err_type == UNRECOG_FILTYPE ||
812 err_type == UNDEF_SEPSTR || err_type == IMPROP_DATA_TYPE ||
813 err_type == RE_COMP_FAILED || err_type == UNRECOG_UNITS) {
814 strncpy(FirstLine,"",MAXLINELEN);
815 test = next_resp(fptr);
816 }
817 else if(err_type == UNDEF_PREFIX) {
818 test = 0;
819 }
820 free_channel(&this_channel);
821 if (!test) {
822 FirstField = 0;
823 strncpy(FirstLine,"",MAXLINELEN);
824 if(!stdio_flag) {
825 fclose(fptr);
826 }
827 }
828 else
829 goto look_again;
830
831 /* if not the last file in the list, move on to the next one */
832
833 if(lst_ptr->next_file != (struct file_list *)NULL) {
834 tmp_ptr = lst_ptr->next_file;
835 lst_ptr = tmp_ptr;
836 }
837
838 }
839 }
840 }
841 /* if not the last station-channel-network in the list,
842 move on to the next one */
843 if(i < (scns->nscn-1)) {
844 flst_ptr = flst_ptr->ptr_next;
845 lst_ptr = flst_ptr->first_list;
846 i++;
847 scn = scns->scn_vec[i];
848 goto next_scn;
849 }
850 if(!stdio_flag)
851 free_matched_files(output_files);
852
853 /* allocated one too many responses */
854
855 free_response(resp);
856 if(prev_ptr != (struct response *)NULL)
857 prev_ptr->next = (struct response *)NULL;
858
859 } /* end else if mode */
860
861 } /* end for loop */
862
863 /* added file close if single input file -- 2/13/2006 -- [ET]: */
864 if(!mode && !stdio_flag) /* if single file was opened then */
865 fclose(fptr); /* close input file */
866
867 /* and print a list of WARNINGS about the station-channel pairs that were not
868 found in the input RESP files */
869
870 for(i = 0; i < scns->nscn; i++) {
871 scn = scns->scn_vec[i];
872 if(!scn->found) {
873 fprintf(stderr,"%s WARNING: no response found for NET=%s,STA=%s,LOCID=%s,CHAN=%s,DATE=%s\n",
874 myLabel, scn->network, scn->station, scn->locid, scn->channel, date_time);
875 fflush(stderr);
876 }
877 }
878 free_scn_list(scns);
879 if(flst_head != (struct matched_files *)NULL)
880 free_matched_files(flst_head);
881 free_string_array(chan_list);
882 free_string_array(locid_list);
883 free_string_array(sta_list);
884
885 free(freqs_orig); /* added 3/28/2006 -- [ET] */
886
887 return(first_resp);
888
889 }
890
891 /* This version of the function does not include
892 the 'listinterp...' parameters */
893
evresp(char * stalst,char * chalst,char * net_code,char * locidlst,char * date_time,char * units,char * file,double * freqs,int nfreqs,char * rtype,char * verbose,int start_stage,int stop_stage,int stdio_flag,int useTotalSensitivityFlag)894 struct response *evresp(char *stalst, char *chalst, char *net_code,
895 char *locidlst, char *date_time, char *units,
896 char *file, double *freqs, int nfreqs,
897 char *rtype, char *verbose, int start_stage,
898 int stop_stage, int stdio_flag, int useTotalSensitivityFlag)
899 {
900 return evresp_itp(stalst,chalst,net_code,locidlst,date_time,units,file,
901 freqs,nfreqs,rtype,verbose,start_stage,stop_stage,
902 stdio_flag,0,0,0.0, 0);
903 }
904
905