Changeset 4194


Ignore:
Timestamp:
04/27/11 12:20:43 (10 years ago)
Author:
scott
Message:

Ewspectra, compute_spectra & sniffspectra

Location:
trunk
Files:
15 added
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/ewdoc/WEB_DOC/modules.html

    r4091 r4194  
    343343    </tr> 
    344344    <tr> 
     345      <td align="center">compute_spectra </td> 
     346      <td>Posts a message for <a href="ovr/ewspectra_ovr.html">ewspectra</a> to begin processing</td> 
     347      <td style="text-align: center;"> <a href="ovr/compute_spectra_ovr.html">overview/setup</a></td> 
     348    </tr> 
     349    <tr> 
    345350      <td align="center">eqproc </td> 
    346351      <td>Initiates final event processing (head of mega-module).</td> 
     
    399404      </td> 
    400405      <td align="center"><a href="cmd/eqverify_assemble_cmd.html">commands</a> </td> 
     406    </tr> 
     407    <tr> 
     408      <td align="center">ewspectra </td> 
     409      <td>Takes data from one or more waveservers, computes and processes their spectra </td> 
     410      <td style="text-align: center;"><a href="ovr/ewspectra_ovr.html">overview</a> 
     411      </td> 
     412      <td align="center"><a href="cmd/ewspectra_cmd.html">commands</a> </td> 
    401413    </tr> 
    402414    <tr> 
     
    10461058      <td align="center">sniffring </td> 
    10471059      <td>Latches onto a user-defined transport ring, reads every 
    1048 nessage and prints logo to screen.</td> 
     1060message and prints logo to screen.</td> 
    10491061      <td><a href="ovr/sniffring_ovr.html">overview</a> </td> 
    10501062      <td align="center"> none </td> 
     1063    </tr> 
     1064    <tr> 
     1065      <td align="center">sniffspectra </td> 
     1066      <td>Latches onto a user-defined transport ring, reads every 
     1067spectra message (as from <a href="ovr/ewspectra_ovr.html">ewspectra</a>);  
     1068prints information to screen or to a file</td> 
     1069      <td style="text-align: center;"> <a href="ovr/sniffspectra_ovr.html">overview/setup</a></td> 
    10511070    </tr> 
    10521071    <tr> 
  • trunk/include/ew_spectra.h

    r4182 r4194  
     1/****************************************************************************** 
     2 * 
     3 *      File:                   ew_spectra.h 
     4 * 
     5 *      Function:               Header for ewspectra module 
     6 * 
     7 *      Author(s):              Scott Hunter, ISTI 
     8 * 
     9 *      Source:                 Started anew. 
     10 * 
     11 *      Notes:                   
     12 * 
     13 *      Change History: 
     14 *                      4/26/11 Started source 
     15 *       
     16 *****************************************************************************/ 
     17 
    118#ifndef SPECTRA_MSG_H 
    219#define SPECTRA_MSG_H 
    320 
    4 #define TRACE2_STA_LEN    7    /* SEED: 5 chars plus terminating NULL */ 
    5 #define TRACE2_NET_LEN    9    /* SEED: 2 chars plus terminating NULL */ 
    6 #define TRACE2_CHAN_LEN   4    /* SEED: 3 chars plus terminating NULL */ 
    7 #define TRACE2_LOC_LEN    3    /* SEED: 2 chars plus terminating NULL */ 
    8 #define LOC_NULL_STRING  "--"  /* NULL string for location code field */ 
     21#include <trace_buf.h> 
     22//#define       TRACE2_STA_LEN    7    /* SEED: 5 chars plus terminating NULL */ 
     23//#define       TRACE2_NET_LEN    9    /* SEED: 2 chars plus terminating NULL */ 
     24//#define       TRACE2_CHAN_LEN   4    /* SEED: 3 chars plus terminating NULL */ 
     25//#define       TRACE2_LOC_LEN    3    /* SEED: 2 chars plus terminating NULL */ 
     26//#define LOC_NULL_STRING  "--"  /* NULL string for location code field */ 
    927 
    1028typedef struct { 
     
    1331        char    chan[TRACE2_CHAN_LEN]; /* Component/channel code (NULL-terminated)*/ 
    1432        char    loc[TRACE2_LOC_LEN];   /* Location code (NULL-terminated) */ 
    15         char    units[2];              /* a for amplitude/phase or c for complex values ,  
     33        char    units[2];              /* a for amplitude/phase, c for complex values, p for peaks  
    1634                                          lower case is Intel byte order, upper case SPARC*/ 
    1735 
     
    2038                                          (seconds since midnight 1/1/1970) */ 
    2139        double  delta_frequency;       /* delta frequency of spectra (Hz) */ 
    22                 char    pad[16];        /* unused, but takes header to 64 bytes */ 
     40        int     numPeaks;              /* number of peaks (valid only if units is p or P */ 
     41                char    pad[12];        /* unused, but takes header to multiple of 8 bytes */ 
    2342} EW_SPECTRA_HEADER; 
    24 /* data follows header as 2 doubles per sample (amp phase or real imaginary) */ 
    25  
    26 #define MAX_EWS_PROCS   10      /* maximum number of spectra processing definitions */ 
     43/* In ring message, data follows header as 2 doubles per sample  
     44        (amp phase or real imaginary or freq amp) */ 
    2745 
    2846#endif 
  • trunk/include/ew_timeseries.h

    r4182 r4194  
     1/****************************************************************************** 
     2 * 
     3 *      File:                   ew_timeseries.h 
     4 * 
     5 *      Function:               Representation of & a library of functions that can handle  
     6 *                  a time-series for processing in Earthworm modules 
     7 * 
     8 *      Author(s):              Paul Friberg, ISTI 
     9 *                  Scott Hunter, ISTI 
     10 * 
     11 *      Source:                 Started anew. 
     12 * 
     13 *      Notes:                   
     14 * 
     15 *      Change History: 
     16 *                      4/26/11 Started source 
     17 *       
     18 *****************************************************************************/ 
    119 
    2  
    3 /* this structure is an attempt to make a library of functions that can 
    4         handle a time-series for processing in Earthworm modules (none exists) 
    5  
    6         February 8, 2011 
    7         Paul Friberg 
    8 */ 
     20#ifndef EW_TIME_SERIES_H 
     21#define EW_TIME_SERIES_H 
    922 
    1023#include <trace_buf.h> 
     
    1326typedef struct  { 
    1427        TRACE2X_HEADER trh2x; 
    15         int gaps;                /* # of gaps in the time-series */ 
    16         double gapTime;  /* amount of time in gaps */ 
    17         double gapValue; /* values used to fill gaps with */ 
     28        int gaps;                 /* # of gaps in the time-series */ 
     29        double gapTime;   /* amount of time in gaps */ 
     30        double gapValue;  /* values used to fill gaps with */ 
    1831        long   dataCount; /* # of entries in data */ 
    19         double *data;    /* a pointer to the data as an array of doubles, of trh2x.nsamp number of samples */ 
     32        int    dataType;  /* type of data in series */ 
     33        double *data;     /* a pointer to the data as an array of doubles, of trh2x.nsamp number of samples */ 
    2034} 
    2135EW_TIME_SERIES ; 
     36 
     37#define EWTS_TYPE_SIMPLE (0)    /* each entry is one double */ 
     38#define EWTS_TYPE_COMPLEX (1)   /* each entry is a complex value: 
     39                                                                        one double for the real part, 
     40                                                                        the next for the imaginary part */ 
     41#define EWTS_TYPE_AMPPHS (2)    /* each entry is a pair of values: 
     42                                                                        one double for the amplitude, 
     43                                                                        the next for the phase */ 
     44 
     45#define EWTS_MODE_FIRST (1)             /* process only 1st component of each entry */ 
     46#define EWTS_MODE_SECOND (2)    /* process only 2nd component of each entry */ 
     47#define EWTS_MODE_BOTH (3)              /* process both components of each entry */ 
    2248 
    2349/***************************************************************************** 
     
    2854int unary_calc_EWTS(  
    2955        EW_TIME_SERIES arg,                                     /* the EWTS to calculate from */ 
    30         void(*op)(const double, void* ),        /* the accumulator function */ 
     56        int mode,                                                       /* which component(s) to process */ 
     57        void(*op)(const double, const int, void* ),     /* the accumulator function */ 
    3158        void *oparg );                                          /* the "accumulator" argument for op */ 
    3259 
     
    3865int unary_modify_EWTS(  
    3966        EW_TIME_SERIES arg,                                             /* the EWTS to modify */ 
    40         double(*op)(const double, const double),/* the modify function */ 
     67        int mode,                                                               /* which component(s) to process */ 
     68        double(*op)(const double, const int, const double),/* the modify function */ 
    4169        const double oparg );                                   /* second argument to all calls to op */ 
    4270 
     
    4977        EW_TIME_SERIES input,                                   /* the EWTS to calculate from */ 
    5078        EW_TIME_SERIES output,                                  /* the EWTS to fill */ 
    51         double(*op)(const double, const double),/* function to compute fill values */ 
     79        int mode,                                                               /* which component(s) to process */ 
     80        double(*op)(const double, const int, const double),/* function to compute fill values */ 
    5281        const double oparg );                                   /* second arg to function */ 
    5382 
     
    6190        EW_TIME_SERIES arg1,                                            /* the EWTS to modify */ 
    6291        EW_TIME_SERIES arg2,                                            /* EWTS with second arguments */ 
    63         double(*op)(const double, const double) );      /* the modify function */ 
     92        int mode,                                                                       /* which component(s) to process */ 
     93        double(*op)(const double, const int, const double) );   /* the modify function */ 
    6494 
    6595/***************************************************************************** 
     
    73103        EW_TIME_SERIES arg2,                                            /* the 2nd EWTS to read */  
    74104        EW_TIME_SERIES output,                                          /* the EWTS to fill */ 
    75         double(*op)(const double, const double) );      /* function to compute fill values */ 
     105        int mode,                                                                       /* which component(s) to process */ 
     106        double(*op)(const double, const int, const double) );   /* function to compute fill values */ 
    76107 
    77108/* Taper types for taper_EWTS */ 
     
    88119        EW_TIME_SERIES ewts,    /* EWTS to taper */ 
    89120        int taper_type,                 /* type of taper to employ */ 
    90         double fraction );              /* fraction of EWTS to taper at each end */ 
     121        double fraction,                /* fraction of EWTS to taper at each end */ 
     122        int mode );                             /* which component(s) to process */ 
    91123 
    92124/***************************************************************************** 
     
    97129int subtract_from_EWTS(  
    98130        EW_TIME_SERIES ewts1,   /* EWTS to subtract from */ 
    99         EW_TIME_SERIES ewts2 ); /* EWTS of values to subtract */ 
     131        EW_TIME_SERIES ewts2,   /* EWTS of values to subtract */ 
     132        int mode );                             /* which component(s) to process */ 
    100133 
    101134/***************************************************************************** 
     
    113146*****************************************************************************/ 
    114147void print_EWTS( 
    115         EW_TIME_SERIES *ewts,   /* EWTS to print */ 
    116         char * filename);               /* name of file to print to */ 
     148        EW_TIME_SERIES *ewts,           /* EWTS to print */ 
     149        char *col1header,                       /* header for first column ("Time" if NULL) */ 
     150        /* map entry index/time to first column value time if NULL */ 
     151                double(*op)(int i, double t, void *arg),         
     152        void *entryMapArg,                      /* what to pass as 3rd arg to entryMap */ 
     153        FILE *fp);                                      /* open file to print to; stdout if NULL */ 
     154 
     155#endif 
  • trunk/include/ws2ts.h

    r4175 r4194  
     1/****************************************************************************** 
     2 * 
     3 *      File:                   ws2ts.h 
     4 * 
     5 *      Function:               Routines for getting data from waveserver(s) and 
     6 *                  converting it to timeseries data 
     7 * 
     8 *      Author(s):              Scott Hunter, ISTI 
     9 * 
     10 *      Source:                 Started anew. 
     11 * 
     12 *      Notes:                   
     13 * 
     14 *      Change History: 
     15 *                      4/26/11 Started source 
     16 *       
     17 *****************************************************************************/ 
     18 
    119#include <ew_timeseries.h> 
    220 
     
    2240 
    2341 
     42/************************************************************************* 
     43 *   init_ws()                                                           * 
     44 *      Reads config file that should contain at least one Waveserver    * 
     45 *          directive (IP and port and optional comment)                 * 
     46 *      callback called for each command not understood                  * 
     47 *      Build & return pointer to WShandle of waveserver info            * 
     48 *************************************************************************/ 
    2449WShandle * init_ws(char * config, char * program_name, int(*callback)(void)); 
    2550 
     51/************************************************************************* 
     52 *   getTraceBufsFromWS()                                                * 
     53 *      Fill TraceReq w/ tracebuf(s) from servers in wsh for SCNL        * 
     54 *      Return WS_ERR_NONE on success                                    * 
     55 *************************************************************************/ 
    2656int getTraceBufsFromWS(WShandle *wsh, char *S, char *C, char *N, char* L,  
    2757        TRACE_REQ * TraceReq); 
    2858 
     59/************************************************************************* 
     60 *   convertTraceBufReq2WEW()                                            * 
     61 *      Return a timeseries computed from  from tracebufs                * 
     62 *************************************************************************/ 
    2963EW_TIME_SERIES * convertTraceBufReq2WEW ( TRACE_REQ * tr, double gapFillValue ); 
    3064 
     65/************************************************************************* 
     66 *   ws2ts_purge()                                                       * 
     67 *      Free datastructures allocated by this library                    * 
     68 *      Any of the arguments that are NULL are ignored                   * 
     69 *************************************************************************/ 
    3170void ws2ts_purge( WShandle *wsh, TRACE_REQ *trh, EW_TIME_SERIES *ewtsh ); 
    3271 
  • trunk/params/earthworm.d

    r3298 r4194  
    9191 Module   MOD_ARC2TRIG       45 
    9292 Module   MOD_RCV_EW         48 
    93  Module   MOD_HELI1          49 
     93 Module   MOD_HELI1          49 
    9494 
    9595 
     
    127127 
    128128# Installation-specific message-type mappings (100-255): 
    129  Message  TYPE_SPECTRA       100 
    130  Message  TYPE_QUAKE2K       101 
    131  Message  TYPE_LINK          102 
    132  Message  TYPE_EVENT2K       103 
    133  Message  TYPE_PAGE          104 
    134  Message  TYPE_KILL          105 
    135  Message  TYPE_DSTDRINK      106 
    136  Message  TYPE_RESTART       107 
    137  Message  TYPE_REQSTATUS     108 
    138  Message  TYPE_STATUS        109 
    139  Message  TYPE_EQDELETE      110 
    140  Message  TYPE_EVENT_SCNL    111 
    141  Message  TYPE_RECONFIG      112 
    142  Message  TYPE_STOP          113  # stop a child. same as kill, except statmgr 
    143                                   # should not restart it 
     129 Message  TYPE_SPECTRA         100 
     130 Message  TYPE_QUAKE2K         101 
     131 Message  TYPE_LINK            102 
     132 Message  TYPE_EVENT2K         103 
     133 Message  TYPE_PAGE            104 
     134 Message  TYPE_KILL            105 
     135 Message  TYPE_DSTDRINK        106 
     136 Message  TYPE_RESTART         107 
     137 Message  TYPE_REQSTATUS       108 
     138 Message  TYPE_STATUS          109 
     139 Message  TYPE_EQDELETE        110 
     140 Message  TYPE_EVENT_SCNL      111 
     141 Message  TYPE_RECONFIG        112 
     142 Message  TYPE_STOP            113  # stop a child. same as kill, except statmgr 
     143                                    # should not restart it 
     144 Message  TYPE_COMPUTE_SPECTRA 114 
     145 Message  TYPE_SPECTRA_DATA    115 
     146 Message  TYPE_SPECTRA_PEAKS   116 
    144147 
    145148 
  • trunk/params/earthworm_global.d

    r4178 r4194  
    267267 Message  TYPE_TD_AMP           34  # time-domain reduced-rate amplitude summary 
    268268                                    #   produced by CISN RAD software & ada2ring 
     269 Message  TYPE_COMPUTE_SPECTRA  50  # produced by compute_spectra & read by ewspectra 
     270                                    #   to initiate reading waveserver(s) to  
     271                                    #   compute spectra 
     272 Message  TYPE_SPECTRA_DATA     51  # spectra results produced by ewspectra 
     273 Message  TYPE_SPECTRA_PEAKS    52  # peaks of spectra results produced by ewspectra 
     274 Message  TYPE_THRESHALARM      53  # alarm message produced by ewthresh 
    269275 
    270276 
  • trunk/src/libsrc/util/makefile.nt

    r4175 r4194  
    210210        complex_math.obj \ 
    211211        convertInstResponse.obj \ 
     212        ew_spectra_io.obj \ 
    212213        fft99.obj \ 
    213214        fft_prep.obj \ 
  • trunk/src/libsrc/util/makefile.sol

    r4175 r4194  
    159159        complex_math.c \ 
    160160        convertInstResponse.c \ 
     161        ew_spectra_io.c \ 
    161162        fft99.c \ 
    162163        fft_prep.c \ 
  • trunk/src/libsrc/util/makefile.ux

    r4175 r4194  
    164164        complex_math.c \ 
    165165        convertInstResponse.c \ 
     166        ew_spectra_io.obj \ 
    166167        fft_prep.c \ 
    167168        fft99.c \ 
  • trunk/src/libsrc/util/ws2ts.c

    r4182 r4194  
     1/****************************************************************************** 
     2 * 
     3 *      File:                   ws2ts.c 
     4 * 
     5 *      Function:               Routines for getting data from waveserver(s) and 
     6 *                  converting it to timeseries data 
     7 * 
     8 *      Author(s):              Scott Hunter, ISTI 
     9 * 
     10 *      Source:                 Started anew. 
     11 * 
     12 *      Notes:                   
     13 * 
     14 *      Change History: 
     15 *                      4/26/11 Started source 
     16 *       
     17 *****************************************************************************/ 
     18 
    119#include <kom.h> 
    220#include <stdio.h> 
     
    172190                        logit("e", "%s Could not get a connection to %s to get menu.\n", 
    173191                          whoami, server); 
     192                    final_ret = -1; 
    174193                    break; 
    175194                case WS_ERR_SOCKET: 
     
    215234                    logit("e", "%s Connection to %s returns error: %d\n", whoami, 
    216235                      server, ret); 
     236                    final_ret = -1; 
    217237            } 
    218238        } 
     
    228248} 
    229249 
    230  
    231  
    232 /************************************************************************* 
    233  *   init_ws (char * config, char * program_name)                        * 
    234  *      Reads config file that should contain at least one Waveserver    * 
    235  *          directive (IP and port and optional comment)                 * 
    236  *      Build & return pointer to WShandle of waveserver info            * 
    237  *************************************************************************/ 
    238250 
    239251WShandle * init_ws(char * config, char * program_name, int(*callback)(void))  
     
    299311} 
    300312 
    301  
    302 /************************************************************************* 
    303  *   getTraceBufsFromWS(WShandle *wsh, char *S, char *C, char *N, char* L, * 
    304  *      double start_time, double duration)                              * 
    305  *      Get tracebuf(s) from servers in wsh for SCNL                     * 
    306  *      Return WS_ERR_NONE on success                                    * 
    307  *************************************************************************/ 
    308313 
    309314int getTraceBufsFromWS(WShandle *wsh, char *S, char *C, char *N, char* L,  
     
    420425} 
    421426 
    422 /* Then to convert the tracebufs to a waveform: */ 
    423 /************************************************************************* 
    424  *   convertTraceBufReq2WEW ( TRACE_REQ * tr, double gapFillValue )      * 
    425  *      Return a waveform from tracebufs                                 * 
    426  *************************************************************************/ 
    427  
    428427EW_TIME_SERIES * convertTraceBufReq2WEW ( TRACE_REQ * tr,  double gapFillValue )  
    429428{ 
     
    436435        int          j; 
    437436        int          gap_count = 0; 
    438         long         nsamp, nfill; 
     437        long         nsamp, nfill, AH_nsamp; 
    439438        long         nfill_max = 0l; 
    440439        long         nsamp_this_scn = 0l; 
     
    483482        } 
    484483        AH_starttime = starttime; 
     484        AH_nsamp = 0; 
    485485        datatype = 'n'; 
    486486        if (wf->datatype[0] == 's' || wf->datatype[0] == 'i') 
     
    557557                } 
    558558                msg_p += eltSize * nsamp; 
     559                AH_nsamp += nsamp; 
    559560           
    560561                /* End-check based on length of snippet buffer */ 
     
    609610   
    610611  my_ewt = malloc( sizeof( *my_ewt ) ); 
     612  my_ewt->dataType = EWTS_TYPE_SIMPLE; 
    611613  my_ewt->data = mydata; 
    612614  my_ewt->gaps = gap_count; 
     
    614616  my_ewt->gapValue = fill; 
    615617  my_ewt->trh2x = *((TRACE2X_HEADER *)( tr->pBuf )); 
     618  my_ewt->trh2x.starttime = AH_starttime; 
     619  my_ewt->trh2x.endtime = endtime; 
     620  my_ewt->trh2x.nsamp = AH_nsamp; 
    616621  my_ewt->dataCount = dataLen; 
    617622  return my_ewt; 
     
    619624 
    620625 
    621 void ws2ts_purge( WShandle *wsh, TRACE_REQ *trh, EW_TIME_SERIES *ewtsh ) { 
     626void ws2ts_purge( WShandle *wsh, TRACE_REQ *trh, EW_TIME_SERIES *ewtsh )  
     627{ 
    622628        if ( wsh != NULL )  
    623629                free( wsh ); 
     
    634640} 
    635641 
    636 int unary_calc_EWTS( EW_TIME_SERIES arg, void(*op)(const double, void* ), void *oparg ) { 
     642int unary_calc_EWTS( EW_TIME_SERIES arg, int mode, 
     643                     void(*op)(const double, const int, void* ), void *oparg ) 
     644{ 
    637645        int i; 
    638         for ( i=0; i<arg.dataCount; i++ )  
    639                 op(arg.data[i], oparg); 
     646        if ( arg.dataType == EWTS_TYPE_SIMPLE ) 
     647                for ( i=0; i<arg.dataCount; i++ ) 
     648                        op(arg.data[i], 0, oparg); 
     649        else 
     650                for ( i=0; i<arg.dataCount; i++ ) { 
     651                        if ( mode & 1 ) 
     652                                op(arg.data[i*2], 0, oparg); 
     653                        if ( mode & 2 ) 
     654                                op(arg.data[i*2+1], 1, oparg); 
     655                } 
    640656        return 0; 
    641657} 
    642658 
    643 int unary_modify_EWTS( EW_TIME_SERIES arg, double(*op)(const double, const double), const double oparg ) { 
     659int unary_modify_EWTS( EW_TIME_SERIES arg, int mode, 
     660                       double(*op)(const double, const int, const double),  
     661                       const double oparg ) 
     662{ 
    644663        int i; 
    645         for ( i=0; i<arg.dataCount; i++ )  
    646                 arg.data[i] = op(arg.data[i], oparg); 
     664        if ( arg.dataType == EWTS_TYPE_SIMPLE ) 
     665                for ( i=0; i<arg.dataCount; i++ ) 
     666                        arg.data[i] = op(arg.data[i], 0, oparg); 
     667        else 
     668                for ( i=0; i<arg.dataCount; i++ ) { 
     669                        if ( mode & 1 ) 
     670                                arg.data[i*2] = op(arg.data[i*2], 0, oparg); 
     671                        if ( mode & 2 ) 
     672                                arg.data[i*2+1] = op(arg.data[i*2+1], 1, oparg); 
     673                } 
    647674        return 0; 
    648675} 
    649676 
    650 int unary_fill_EWTS( EW_TIME_SERIES input, EW_TIME_SERIES output, double(*op)(const double, const double), const double oparg ) { 
     677int unary_fill_EWTS( EW_TIME_SERIES input, EW_TIME_SERIES output, int mode, 
     678                                        double(*op)(const double, const int, const double),  
     679                                        const double oparg )  
     680{ 
    651681        int i; 
    652         if ( input.trh2x.nsamp != output.trh2x.nsamp ) 
     682        if ( input.trh2x.nsamp != output.trh2x.nsamp ) { 
     683                logit("e", "input & output timeseries' lengths differ\n"); 
    653684                return 1; 
    654         for ( i=0; i<input.dataCount; i++ ) 
    655                 output.data[i] = op(input.data[i], oparg); 
     685        } 
     686        if ( input.dataType != output.dataType ) { 
     687                logit("e", "input & output timeseries' dataTypes differ\n"); 
     688                return 1; 
     689        } 
     690        if ( input.dataType == EWTS_TYPE_SIMPLE ) 
     691                for ( i=0; i<input.dataCount; i++ ) 
     692                        output.data[i] = op(input.data[i], 0, oparg); 
     693        else 
     694                for ( i=0; i<input.dataCount; i++ ) { 
     695                        if ( mode & 1 )  
     696                                output.data[i*2] = op(input.data[i*2], 0, oparg); 
     697                        if ( mode & 2 )  
     698                                output.data[i*2+1] = op(input.data[i*2+1], 1, oparg); 
     699                } 
    656700        return 0; 
    657701} 
    658702 
    659 int binary_modify_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2, double(*op)(const double, const double) ) { 
     703int binary_modify_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2, int mode, 
     704                                        double(*op)(const double, const int, const double) )  
     705{ 
    660706        int i; 
    661707        if ( arg1.trh2x.nsamp != arg2.trh2x.nsamp ) 
    662708                return 1; 
    663         for ( i=0; i<arg1.dataCount; i++ ) 
    664                 arg1.data[i] = op(arg1.data[i], arg2.data[i]); 
     709        if ( arg1.trh2x.nsamp != arg2.trh2x.nsamp ) { 
     710                logit("e", "input timeseries' lengths differ\n"); 
     711                return 1; 
     712        } 
     713        if ( arg1.dataType != arg2.dataType ) { 
     714                logit("e", "input timeseries' dataTypes differ\n"); 
     715                return 1; 
     716        } 
     717        if ( arg1.dataType == EWTS_TYPE_SIMPLE ) 
     718                for ( i=0; i<arg1.dataCount; i++ ) 
     719                        arg1.data[i] = op(arg1.data[i], arg2.data[i], i%2); 
     720        else 
     721                for ( i=0; i<arg1.dataCount; i++ ) { 
     722                        if ( mode & 1 ) 
     723                                arg1.data[i*2] = op(arg1.data[i*2], arg2.data[i*2], 0); 
     724                        if ( mode & 2 ) 
     725                                arg1.data[i*2+1] = op(arg1.data[i*2+1], arg2.data[i*2+1], 1); 
     726                } 
    665727        return 0; 
    666728} 
    667729 
    668 int binary_fill_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2, EW_TIME_SERIES output, double(*op)(const double, const double) ) { 
     730int binary_fill_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2, EW_TIME_SERIES output,  
     731                                        int mode,  
     732                                        double(*op)(const double, const int, const double) )  
     733{ 
    669734        int i; 
    670         if ( arg1.dataCount != arg2.dataCount || arg1.dataCount != output.dataCount ) 
     735        if ( arg1.dataCount != arg2.dataCount ) { 
     736                logit("e", "input timeseries' lengths differ\n"); 
    671737                return 1; 
    672         for ( i=0; i<arg1.dataCount; i++ ) 
    673                 output.data[i] = op(arg1.data[i], arg2.data[i]); 
     738        } 
     739        if ( arg1.dataCount != output.dataCount ) { 
     740                logit("e", "input & output timeseries' lengths differ\n"); 
     741                return 1; 
     742        } 
     743        if ( arg1.dataType != arg2.dataType ) { 
     744                logit("e", "input timeseries' dataTypes differ\n"); 
     745                return 1; 
     746        } 
     747        if ( arg1.dataType != output.dataType ) { 
     748                logit("e", "input & output timeseries' dataTypes differ\n"); 
     749                return 1; 
     750        } 
     751        if ( arg1.dataType == EWTS_TYPE_SIMPLE ) 
     752                for ( i=0; i<arg1.dataCount; i++ ) 
     753                        output.data[i] = op(arg1.data[i], arg2.data[i], 0); 
     754        else  
     755                for ( i=0; i<arg1.dataCount; i++ ) { 
     756                        if ( mode & 1 ) 
     757                                output.data[i*2] = op(arg1.data[i*2], arg2.data[i*2], 0); 
     758                        if ( mode & 2 ) 
     759                                output.data[i*2+1] = op(arg1.data[i*2+1], arg2.data[i*2+1], 1); 
     760                } 
    674761        return 0; 
    675762} 
    676763 
    677 void taper_EWTS( EW_TIME_SERIES ewts, int taper_type, double fraction ) { 
     764void taper_EWTS( EW_TIME_SERIES ewts, int taper_type, double fraction, int mode )  
     765{ 
    678766        /* taper length in samples */ 
    679767        int nsamp = ewts.dataCount; 
     
    686774        switch (taper_type) { 
    687775        case BARTLETT: 
    688                 for (i= 0; i < (int) ffl; i++) { 
    689                         f= 1.0*((float) i)/ffl; 
    690                         fdata[i]*= f; 
    691                         fdata[nsamp-1-i]*= f; 
    692                 } 
     776                if ( ewts.dataType == EWTS_TYPE_SIMPLE )  
     777                        for (i= 0; i < (int) ffl; i++) { 
     778                                f= 1.0*((float) i)/ffl; 
     779                                fdata[i]*= f; 
     780                                fdata[nsamp-1-i]*= f; 
     781                        } 
     782                else 
     783                        for (i= 0; i < (int) ffl; i++) { 
     784                                f= 1.0*((float) i)/ffl; 
     785                                if ( mode & 1 ) {                                
     786                                        fdata[i*2]*= f; 
     787                                        fdata[(nsamp-1-i)*2]*= f; 
     788                                } 
     789                                if ( mode & 2 ) { 
     790                                        fdata[i*2+1]*= f; 
     791                                        fdata[(nsamp-1-i)*2+1]*= f; 
     792                                } 
     793                        }                        
    693794                break; 
    694795        case PARZAN: 
    695                 for (i= 0; i < (int) ffl; i++) { 
    696                         x= 1.0*((float) i)/(ffl+1.); 
    697                         if (x-0.5 <= 0) 
    698                                 f= 6.0*x*x-6.0*x*x*x; 
    699                         else 
    700                                 f= 1.0-2.0*pow((double) (1.0-x), (double) 3.0); 
    701                         fdata[i]*= f; 
    702                         fdata[nsamp-1-i]*= f; 
    703                 } 
     796                if ( ewts.dataType == EWTS_TYPE_SIMPLE )  
     797                        for (i= 0; i < (int) ffl; i++) { 
     798                                x= 1.0*((float) i)/(ffl+1.); 
     799                                if (x-0.5 <= 0) 
     800                                        f= 6.0*x*x-6.0*x*x*x; 
     801                                else 
     802                                        f= 1.0-2.0*pow((double) (1.0-x), (double) 3.0); 
     803                                fdata[i]*= f; 
     804                                fdata[nsamp-1-i]*= f; 
     805                        } 
     806                else  
     807                        for (i= 0; i < (int) ffl; i++) { 
     808                                x= 1.0*((float) i)/(ffl+1.); 
     809                                if (x-0.5 <= 0) 
     810                                        f= 6.0*x*x-6.0*x*x*x; 
     811                                else 
     812                                        f= 1.0-2.0*pow((double) (1.0-x), (double) 3.0); 
     813                                if ( mode & 1 ) {                                
     814                                        fdata[i*2]*= f; 
     815                                        fdata[(nsamp-1-i)*2]*= f; 
     816                                } 
     817                                if ( mode & 2 ) { 
     818                                        fdata[i*2+1]*= f; 
     819                                        fdata[(nsamp-1-i)*2+1]*= f; 
     820                                } 
     821                        } 
    704822                break; 
    705823        case HANNING: 
    706                 for (i= 0; i < (int) ffl; i++) { 
    707                         f= 0.5-0.5*cos(PI*((float) i)/(ffl+1.)); 
    708                         fdata[i]*= f; 
    709                         fdata[nsamp-1-i]*= f; 
    710                 } 
     824                if ( ewts.dataType == EWTS_TYPE_SIMPLE )  
     825                        for (i= 0; i < (int) ffl; i++) { 
     826                                f= 0.5-0.5*cos(PI*((float) i)/(ffl+1.)); 
     827                                fdata[i]*= f; 
     828                                fdata[nsamp-1-i]*= f; 
     829                        } 
     830                else 
     831                        for (i= 0; i < (int) ffl; i++) { 
     832                                f= 0.5-0.5*cos(PI*((float) i)/(ffl+1.)); 
     833                                if ( mode & 1 ) {                                
     834                                        fdata[i*2]*= f; 
     835                                        fdata[(nsamp-1-i)*2]*= f; 
     836                                } 
     837                                if ( mode & 2 ) { 
     838                                        fdata[i*2+1]*= f; 
     839                                        fdata[(nsamp-1-i)*2+1]*= f; 
     840                                } 
     841                        } 
    711842                break; 
    712843        case BMHARRIS: 
    713                 for (i= 0; i < (int) ffl; i++) { 
    714                         f= 0.35825- 
    715                         0.48829*cos(PI*((float) i)/(ffl+1.))+ 
    716                         0.14128*cos(2.0*PI*((float) i)/(ffl+1.))- 
    717                         0.01168*cos(3.0*PI*((float) i)/(ffl+1.)); 
    718                         fdata[i]*= f; 
    719                         fdata[nsamp-1-i]*= f; 
    720                 } 
     844                if ( ewts.dataType == EWTS_TYPE_SIMPLE )  
     845                        for (i= 0; i < (int) ffl; i++) { 
     846                                f= 0.35825- 
     847                                0.48829*cos(PI*((float) i)/(ffl+1.))+ 
     848                                0.14128*cos(2.0*PI*((float) i)/(ffl+1.))- 
     849                                0.01168*cos(3.0*PI*((float) i)/(ffl+1.)); 
     850                                if ( mode & 1 ) {                                
     851                                        fdata[i*2]*= f; 
     852                                        fdata[(nsamp-1-i)*2]*= f; 
     853                                } 
     854                                if ( mode & 2 ) { 
     855                                        fdata[i*2+1]*= f; 
     856                                        fdata[(nsamp-1-i)*2+1]*= f; 
     857                                } 
     858                        } 
     859                else 
     860                        for (i= 0; i < (int) ffl; i++) { 
     861                                f= 0.35825- 
     862                                0.48829*cos(PI*((float) i)/(ffl+1.))+ 
     863                                0.14128*cos(2.0*PI*((float) i)/(ffl+1.))- 
     864                                0.01168*cos(3.0*PI*((float) i)/(ffl+1.)); 
     865                                if ( mode & 1 ) {                                
     866                                        fdata[i*2]*= f; 
     867                                        fdata[(nsamp-1-i)*2]*= f; 
     868                                } 
     869                                if ( mode & 2 ) { 
     870                                        fdata[i*2+1]*= f; 
     871                                        fdata[(nsamp-1-i)*2+1]*= f; 
     872                                } 
     873                        } 
    721874                break; 
    722875        } 
    723876} 
    724877 
    725 double diff( const double a, const double b ) { 
     878double diff( const double a, int odd, const double b ) 
     879{ 
    726880        return a - b; 
    727881} 
    728882 
    729 int subtract_from_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2 ) { 
    730         return binary_modify_EWTS( arg1, arg2, diff ); 
    731 } 
    732  
    733 void addto( const double a, void * b ) { 
     883int subtract_from_EWTS( EW_TIME_SERIES arg1, EW_TIME_SERIES arg2, int mode )  
     884{ 
     885        return binary_modify_EWTS( arg1, arg2, mode, diff ); 
     886} 
     887 
     888void addto( const double a, int odd , void * b ) 
     889{ 
    734890        double * bptr = (double*)b; 
    735         *bptr += a; 
    736 } 
    737  
    738 int demean_EWTS( EW_TIME_SERIES arg ) { 
     891        bptr[odd] += a; 
     892} 
     893 
     894int demean_EWTS( EW_TIME_SERIES arg )  
     895{ 
    739896        int i; 
    740         double sum = 0, avg; 
    741         unary_calc_EWTS( arg, addto, &sum ); 
    742         avg = sum/arg.dataCount; 
    743         unary_modify_EWTS( arg, diff, avg ); 
     897        double sum[2] = {0,0}; 
     898        unary_calc_EWTS( arg, EWTS_MODE_BOTH, addto, sum ); 
     899        unary_modify_EWTS( arg, EWTS_MODE_FIRST, diff, sum[0]/arg.dataCount ); 
     900        if ( arg.dataType != EWTS_TYPE_SIMPLE ) 
     901                unary_modify_EWTS( arg, EWTS_MODE_SECOND, diff, sum[1]/arg.dataCount ); 
    744902        return 0; 
    745903} 
    746904 
    747 void print_EWTS(EW_TIME_SERIES *ts, char * filename) { 
    748         FILE *fp = stdout; 
     905static double identity( int i, double t, void *arg )  
     906{ 
     907        return t; 
     908} 
     909 
     910void print_EWTS(EW_TIME_SERIES *ts, char *col1header,  
     911        double(*op)(int i, double t, void *arg), 
     912        void *entryMapArg, 
     913        FILE *fp)  
     914{ 
    749915        int i; 
    750916        double step = 1.0/ts->trh2x.samprate; 
    751917        double t = ts->trh2x.starttime; 
    752         if ( filename != NULL )  
    753                 if (! (fp = fopen (filename, "w"))) { 
    754                         logit ("e", "Unable to open file %s\n", filename); 
    755                         return; 
    756                 } 
    757     for ( i=0; i<ts->dataCount; i++, t+=step ) 
    758         fprintf( fp, "%15.3lf\t%15.3lf\n", t, ts->data[i] ); 
    759     if ( filename != NULL ) 
    760             fclose( fp ); 
    761 } 
     918        if ( col1header == NULL ) 
     919                col1header = "Time"; 
     920        if ( op == NULL ) 
     921                op = identity; 
     922        if ( fp == NULL ) 
     923                fp == stdout; 
     924        if ( ts->dataType == EWTS_TYPE_SIMPLE ) { 
     925                fprintf( fp, "%15s\t%15s\n", col1header, "Value" ); 
     926            for ( i=0; i<ts->dataCount; i++, t+=step ) 
     927                fprintf( fp, "%15.3lf\t%15.3lf\n", op(i,t, entryMapArg), ts->data[i] ); 
     928    } else { 
     929        if ( ts->dataType == EWTS_TYPE_COMPLEX ) 
     930                        fprintf( fp, "%15s\t%15s\t%15s\n", col1header, "Real", "Imaginary" ); 
     931        else 
     932                        fprintf( fp, "%15s\t%15s\t%15s\n", col1header, "Amplitude", "Phase" ); 
     933            for ( i=0; i<ts->dataCount; i++, t+=step ) 
     934                fprintf( fp, "%15.3lf\t%15.3lf\t%15.3lf\n", op(i,t, entryMapArg), ts->data[i*2], ts->data[i*2+1] ); 
     935        } 
     936} 
  • trunk/src/seismic_processing/Makefile

    r4108 r4194  
    4949        carlsubtrig \ 
    5050        compress_UA \ 
     51        compute_spectra \ 
    5152        debias \ 
    5253        decimate \ 
     
    6566        ew2rsam \ 
    6667        ew2ssam \ 
     68        ewspectra \ 
    6769        fir \ 
    6870        localmag \ 
     
    7880        rayloc_ew \ 
    7981        raypicker \ 
     82        sniffspectra \ 
    8083        statrigfilter \ 
    8184        wftimefilter 
     
    8992        carlsubtrig \ 
    9093        compress_UA \ 
     94        compute_spectra \ 
    9195        debias \ 
    9296        decimate \ 
     
    104108        ew2rsam \ 
    105109        ew2ssam \ 
     110        ewspectra \ 
    106111        fir \ 
    107112        geqproc \ 
     
    116121        pkfilter \ 
    117122        raypicker \ 
     123        sniffspectra \ 
    118124        statrigfilter \ 
    119125        wftimefilter 
     
    132138        carlsubtrig \ 
    133139        compress_UA \ 
     140        compute_spectra \ 
    134141        debias \ 
    135142        decimate \ 
     
    148155        ew2rsam \ 
    149156        ew2ssam \ 
     157        ewspectra \ 
    150158        fir \ 
    151159        localmag \ 
     
    160168        rayloc_ew \ 
    161169        raypicker \ 
     170        sniffspectra \ 
    162171        statrigfilter \ 
    163172        wftimefilter 
  • trunk/src/seismic_processing/ewspectra/ewspectra.c

    r4184 r4194  
     1/****************************************************************************** 
     2 * 
     3 *      File:                   ewspectra.c 
     4 * 
     5 *      Function:               Earthworm module for computing spectra from waveserver data 
     6 * 
     7 *      Author(s):              Scott Hunter, ISTI 
     8 * 
     9 *      Source:                 Started anew. 
     10 * 
     11 *      Notes:                   
     12 * 
     13 *      Change History: 
     14 *                      4/26/11 Started source 
     15 *       
     16 *****************************************************************************/ 
     17 
    118#include <stdio.h> 
    219#include <stdlib.h> 
     
    825#include <kom.h> 
    926#include <ew_spectra.h> 
    10  
    11 /* 
    12 #include <ew_timeseries.h> 
    13 */ 
     27#include <transport.h> 
     28#include "iir.h" 
     29#include <ew_spectra_io.h> 
     30 
     31#define MAX_EWS_PROCS   10      /* maximum number of spectra processing definitions */ 
     32 
     33#define MSGSTK_OFF    0              /* MessageStacker has not been started      */ 
     34#define MSGSTK_ALIVE  1              /* MessageStacker alive and well            */ 
     35#define MSGSTK_ERR   -1              /* MessageStacker encountered error quit    */ 
     36volatile int MessageStackerStatus = MSGSTK_OFF; 
     37 
     38/* Error messages used by export 
     39 ***********************************/ 
     40#define  ERR_MISSMSG       0   /* message missed in transport ring        */ 
     41#define  ERR_TOOBIG        1   /* retreived msg too large for buffer      */ 
     42#define  ERR_NOTRACK       2   /* msg retreived; tracking limit exceeded  */ 
     43#define  ERR_BADMSG        3   /* message w/ bad timespan                 */ 
     44static char  errText[256];     /* string for log/error messages           */ 
     45 
     46#define backward 0 
    1447 
    1548char *progname; 
    1649 
    17 /*********************************************************************** 
    18  * ConvertTime () - given pStart return a double representing          * 
    19  *     number of seconds since 1970                                    * 
    20  ***********************************************************************/ 
    21 static int ConvertTime (char *pStart, double *start) 
    22 { 
    23     char                strDate[30]; 
    24     struct tm   tmDate; 
    25      
    26     if (pStart == NULL) { 
    27         logit ("e", "Invalid parameters passed in.\n"); 
    28         return EW_FAILURE; 
    29     } 
    30      
    31     strncpy (strDate, pStart, 14); 
    32     strDate[14] = '\0'; 
    33         tmDate.tm_sec   = atoi( strDate+12 ); 
    34         strDate[12] = '\0'; 
    35         tmDate.tm_min   = atoi( strDate+10 ); 
    36         strDate[10] = '\0'; 
    37         tmDate.tm_hour  = atoi( strDate+8 ); 
    38         strDate[ 8] = '\0'; 
    39         tmDate.tm_mday  = atoi( strDate+6 ); 
    40         strDate[ 6] = '\0'; 
    41         tmDate.tm_mon   = atoi( strDate+4 ); 
    42         strDate[ 4] = '\0'; 
    43         tmDate.tm_year  = atoi( strDate ) - 1900; 
    44     tmDate.tm_wday = tmDate.tm_yday = tmDate.tm_isdst = tmDate.tm_gmtoff = 0; 
    45     tmDate.tm_zone = "GMT"; 
    46  
    47     *start = (double)timegm( &tmDate ); 
    48      
    49     return EW_SUCCESS; 
    50 } 
     50/* Things to look up in the earthworm.h tables with getutil.c functions 
     51 **********************************************************************/ 
     52static long          InRingKey;     /* key of transport ring for input    */ 
     53static long          OutRingKey;    /* key of transport ring for output   */ 
     54static unsigned char InstId;        /* local installation id              */ 
     55static unsigned char MyModId;       /* Module Id for this program         */ 
     56 
     57static  SHM_INFO  InRegion;     /* shared memory region to use for input  */ 
     58static  SHM_INFO  OutRegion;    /* shared memory region to use for output */ 
     59 
     60MSG_LOGO  GetLogo;                      /* requesting module,type,instid */ 
     61pid_t MyPid;        /* Our own pid, sent with heartbeat for restart purposes */ 
     62         
     63time_t now;        /* current time, used for timing heartbeats */ 
     64time_t MyLastBeat;         /* time of last local (into Earthworm) hearbeat */ 
     65 
     66static int     HeartBeatInt;        /* seconds between heartbeats        */ 
     67static unsigned char TypeHeartBeat; 
     68static unsigned char TypeError; 
     69 
     70static char *Rawmsg = NULL;          /* "raw" retrieved msg for main thread      */ 
     71static char *MSrawmsg = NULL;        /* MessageStacker's "raw" retrieved message */ 
     72 
     73int static long    MaxMsgSize = 100;          /* max size for input msgs    */ 
     74 
     75WShandle *wsh; 
     76 
     77#define THREAD_STACK 8192 
     78static unsigned tidStacker;          /* Thread moving messages from transport */ 
     79 
     80int numPeaks = 0; 
     81double peakBand[2]; 
    5182 
    5283 
     
    179210int smooth;                                                     /* 0=no smoothing */ 
    180211double smooth_width;                            /* width (in secs) of smoothing window */ 
    181 char SCNL[MAX_EWS_PROCS][2][4][10];     /* SCNLs to process */ 
     212char SCNL[MAX_EWS_PROCS][3][4][10];     /* SCNLs to process */ 
    182213char binary[MAX_EWS_PROCS];                     /* 1=use 2 SCNLs for processing */ 
    183214char processing_mode[MAX_EWS_PROCS];/* what processing to do: 
     
    193224char outpath[100];                                      /* path for output file */ 
    194225FILE *fp;                                                       /* output file pointer */ 
     226char inRing[MAX_RING_STR];                      /* name of input ring */ 
     227char outRing[MAX_RING_STR];                     /* name of output ring */ 
     228int find_triggers; 
     229double highcutoff, lowcutoff;           /* freq for high and low cutoffs */ 
     230int highpoles, lowpoles;                        /* # poles for high & low cutoffs (0=none) */ 
     231 
     232char *taperName[5] = {"", "BARTLETT", "HANNING", "PARZAN", "BMHARRIS"}; 
     233char *Argv0 = "ewspectra"; 
     234 
     235static char    MyModName[MAX_MOD_STR];       /* speak as this module name/id      */ 
     236int spectraOut = 0; 
    195237 
    196238/*********************************************************************** 
     
    199241 ***********************************************************************/ 
    200242int parse_my_command() { 
    201         if (k_its("EWS_Scale")) { 
     243        if (k_its("Scale")) { 
    202244                scale = 1; 
    203         } else if (k_its("EWS_White")) { 
     245        } else if (k_its("White")) { 
    204246                white = 1; 
    205         } else if (k_its("EWS_Smooth")) { 
     247        } else if (k_its("Smooth")) { 
    206248                smooth = 1; 
    207249                smooth_width = k_val(); 
    208         } else if (k_its("EWS_Taper")) { 
     250        } else if (k_its("Taper")) { 
    209251                char *taperTypeName; 
    210252                taperTypeName = k_str(); 
    211                 if ( strcmp( taperTypeName, "BARTLETT" )==0 ) 
    212                         taperType = BARTLETT; 
    213                 else if ( strcmp( taperTypeName, "HANNING" )==0 ) 
    214                         taperType = HANNING; 
    215                 else if ( strcmp( taperTypeName, "PARZAN" )==0 ) 
    216                         taperType = PARZAN; 
    217                 else if ( strcmp( taperTypeName, "BMHARRIS" )==0 ) 
    218                         taperType = BMHARRIS; 
    219                 else { 
     253                for ( taperType = 4; taperType>0; taperType-- ) 
     254                        if ( strcmp( taperTypeName, taperName[taperType] )==0 ) 
     255                                break; 
     256                if ( taperType == 0 ) { 
    220257                        logit( "e", "Invalid type of taper: '%s'\n", taperTypeName ); 
    221258                        return 1; 
     
    223260                taperFrac = k_val(); 
    224261        } else if (k_its("MyModuleId")) { 
    225         } else if (k_its("DiffSpectraSCNLs")) { 
    226                 int i,j; 
     262                strcpy( MyModName, k_str() ); 
     263                if ( GetModId( MyModName, &MyModId ) != 0 ) { 
     264                  logit( "e", "%s: Invalid module name <%s>; exiting!\n", 
     265                                   Argv0, MyModName ); 
     266                  return 1; 
     267                } 
     268        } else if (k_its("DeconvolveSpectraSCNLs") || k_its("DiffSpectraSCNLs") 
     269                || k_its("PlainSpectraSCNL")) { 
     270                int i,j,step; 
    227271                if ( num_ews_procs >= MAX_EWS_PROCS ) { 
    228272                        logit( "e", "Too many processing requests; ignoring '%s'\n", k_com() ); 
    229273                        return 0; 
    230274                } 
    231                 for ( i=0; i<2; i++ ) 
    232                         for ( j=0; j<4; j++ ) 
    233                                 strcpy( SCNL[num_ews_procs][i][j], k_str() ); 
    234                 processing_mode[num_ews_procs] = 'd'; 
    235                 binary[num_ews_procs] = 1; 
    236                 num_ews_procs++; 
    237         } else if (k_its("PlainSpectraSCNL")) { 
    238                 int i,j; 
    239                 if ( num_ews_procs >= MAX_EWS_PROCS ) { 
    240                         logit( "e", "Too many processing requests; ignoring '%s'\n", k_com() ); 
    241                         return 0; 
    242                 } 
    243                 for ( j=0; j<4; j++ ) 
    244                         strcpy( SCNL[num_ews_procs][0][j], k_str() ); 
    245                 processing_mode[num_ews_procs] = 'p'; 
    246                 binary[num_ews_procs] = 0; 
     275                step = k_its("PlainSpectraSCNL") ? 2 : 1; 
     276                *SCNL[num_ews_procs][2][3] = 0; 
     277                for ( i=0; i<3; i+=step ) 
     278                        for ( j=0; j<4; j++ ) { 
     279                                char *str = k_str(); 
     280                                if ( str == NULL ) 
     281                                        if ( i < 2 ) { 
     282                                                logit( "e", "Incomplete source SCNL for %s\n", k_com() ); 
     283                                                return 0; 
     284                                        } else if ( j > 0 ) { 
     285                                                logit( "e", "Target SCNL must be fully specified for %s\n", k_com() ); 
     286                                                return 0; 
     287                                        } else 
     288                                                break; 
     289                                strcpy( SCNL[num_ews_procs][i][j], str ); 
     290                        } 
     291                processing_mode[num_ews_procs] = step==2 ? 'p' : k_its("DeconvolveSpectraSCNLs") ? 'd' : 's'; 
     292                binary[num_ews_procs] = (step == 1); 
    247293                num_ews_procs++; 
    248294        } else if (k_its("TimeSpan")) { 
     
    251297                if ( arg1 < 0 ) {  
    252298                        start = time(NULL) + arg1; 
    253                 } else if ( ConvertTime (arg1s, &start) == EW_FAILURE ) { 
     299                } else if ( EWSConvertTime (arg1s, &start) == EW_FAILURE ) { 
    254300                        return 1; 
    255301                } 
     
    257303        } else if (k_its("OutFile")) { 
    258304                strcpy( outpath, k_str() ); 
     305        } else if (k_its("InRing")) { 
     306                find_triggers = 1; 
     307                strcpy( inRing, k_str() ); 
     308                if ( ( InRingKey = GetKey(inRing) ) == -1 ) { 
     309                        logit( "e", "%s:  Invalid ring name <%s>; exiting!\n", 
     310                                         Argv0, inRing); 
     311                        return 1; 
     312                } 
     313        } else if (k_its("OutRing")) { 
     314                strcpy( outRing, k_str() ); 
     315            if ( ( OutRingKey = GetKey(outRing) ) == -1 ) { 
     316                        logit( "e", "%s:  Invalid ring name <%s>; exiting!\n", 
     317                                         Argv0, outRing); 
     318                        return 1; 
     319                } 
     320        } else if (k_its("HighCut")) { 
     321                highcutoff = k_val(); 
     322                if ( highcutoff <= 0 ) { 
     323                        logit( "e", "HighCut cutoff frequency must be > 0\n" ); 
     324                        return 1; 
     325                } 
     326                highpoles = k_int(); 
     327                if ( highpoles % 2 || highpoles < 2 || highpoles > 12 )  { 
     328                        logit( "e", "# poles for HighCut must be 2, 4, 6, 8, 10 or 12\n" ); 
     329                        return 1; 
     330                } 
     331        } else if (k_its("LowCut")) { 
     332                lowcutoff = k_val(); 
     333                if ( lowcutoff <= 0 ) { 
     334                        logit( "e", "LowCut cutoff frequency must be > 0\n" ); 
     335                        return 1; 
     336                } 
     337                lowpoles = k_int(); 
     338                if ( lowpoles % 2 || lowpoles < 2 || lowpoles > 12 )  { 
     339                        logit( "e", "# poles for LowCut must be 2, 4, 6, 8, 10 or 12\n" ); 
     340                        return 1; 
     341                } 
     342        } else if (k_its("ReportSpectra")) { 
     343                spectraOut = 1; 
     344        } else if (k_its("ReportPeaks")) { 
     345                numPeaks = k_int(); 
     346                peakBand[0] = k_val(); 
     347                peakBand[1] = k_val(); 
    259348        } else 
    260349                return 1; 
     
    263352 
    264353 
     354 
    265355/*********************************************************************** 
    266  * calc_spectra () - calculate the spectra for a timeseries 
     356 * deconvolve_spectra () - deconvolve spectra b into a 
    267357 *              return 0 if no problems, 1 otherwise 
    268358 ***********************************************************************/ 
    269 int calc_spectra( EW_TIME_SERIES * ewts ) { 
    270         char str[30]; 
    271         double amp, amp_sm, phase, nyquist; 
    272         double *spec, *specsm, mean, sum, timespace, timendata; 
    273         double amp_sum, phase_sum; 
    274         int i, j, ndata, smoother, err; 
    275         char *cp; 
    276         long ndata2; 
    277         double delta = 1/ewts->trh2x.samprate; 
    278  
    279         ndata = ewts->dataCount; 
    280  
    281         timendata = ndata; 
    282         timespace = delta; 
    283  
    284         /* get nearest power of 2 for listed data length */ 
    285         for (i= 1; (int) pow((double) 2., (double) i) < ewts->dataCount; i++); 
    286         ndata= (int) pow((double) 2.,(double) (i+1)); 
    287  
    288         /* allocate space for data and smoothing area  
    289         if ((spec= (double *) calloc((unsigned) (ndata+2), sizeof(double))) == NULL ) { */ 
    290         if ((spec= (double *) malloc((unsigned) (ndata+2)*sizeof(double))) == NULL ) { 
     359int deconvolve_spectra( EW_TIME_SERIES *spec_a, EW_TIME_SERIES *spec_b ) { 
     360        int ndata = spec_a->dataCount * 2; 
     361        int i; 
     362        double denom, result_0, result_1; 
     363        double *a = spec_a->data, *b = spec_b->data; 
     364        for ( i=0; i<=ndata; i+=2 ) { 
     365                denom = b[i]*b[i] + b[i+1]*b[i+1]; 
     366                if ( denom == 0.0 ) 
     367                        a[i] = a[i+1] = 0.0; 
     368                else { 
     369                        result_0 = a[i]*b[i] + a[i+1]*b[i+1]; 
     370                        result_1 = b[i]*a[i+1] - a[i]*b[i+1]; 
     371                        a[i]   = result_0 / denom; 
     372                        a[i+1] = result_1 / denom; 
     373                } 
     374                         
     375        } 
     376        return 0; 
     377} 
     378 
     379static double iScale( int i, double t, void* scale ) { 
     380        return i*(*(double*)scale); 
     381} 
     382 
     383/*********************************************************************** 
     384 * smooth_spectra () - smooth the spectra timeseries 
     385 *              return 0 if no problems, 1 otherwise 
     386 ***********************************************************************/ 
     387int smooth_spectra( EW_TIME_SERIES * ewspec, int mode ) { 
     388        double *specsm, even_sum, odd_sum; 
     389        int i, imod; 
     390        double delta = 1/ewspec->trh2x.samprate; 
     391        int smoother= smooth_width/delta; 
     392        int ndata = ewspec->dataCount * 2; 
     393        double *spec = ewspec->data; 
     394         
     395        /* make sure that smoother is even */ 
     396        smoother= (smoother/2)*2; 
     397        if ((specsm= (double *) malloc((smoother)*sizeof(double))) == NULL ) { 
    291398                logit("et","%s: memory allocation error\n", progname); 
    292399                return 1; 
    293400        } 
    294         /*for ( i=0; i<ewts->dataCount; i++ ) 
    295                 spec[i] = ewts->data[i];*/ 
    296         memcpy( spec, ewts->data, ewts->dataCount*sizeof(double) ); 
    297         bzero( spec+ewts->dataCount, sizeof(double)*(ndata-ewts->dataCount+2) ); 
    298  
    299         /* set header variables */ 
    300         nyquist= 1.0/(2.0*delta); 
    301         ndata2= (long) ((ndata/2.0)+1.0); 
    302         delta= nyquist/((double) (ndata2-1)); 
    303  
    304         if (smooth) { 
    305                 /*if ((specsm= (double *) calloc((unsigned) (ndata+2), sizeof(double))) == NULL ) { 
    306                 */ 
    307                 smoother= smooth_width/delta; 
    308                 /* make sure that smoother is even */ 
    309                 smoother= (smoother/2)*2; 
    310                 if ((specsm= (double *) malloc((smoother)*sizeof(double))) == NULL ) { 
    311                         logit("et","%s: memory allocation error\n", progname); 
    312                         free( spec ); 
    313                         return 1; 
    314                 } 
    315         } else 
    316                 specsm = NULL; 
    317  
    318  
    319         ah_fftr(spec, ndata); 
    320         for (i= 0, sum= 0.; i <= ndata; i+= 2) { 
     401        /* Create & fill smoothing window */ 
     402        memcpy( specsm, spec, smoother*2*sizeof(double) ); 
     403        even_sum = odd_sum = 0; 
     404        for ( i=0; i<smoother*2; i+=2) { 
     405                if ( mode & EWTS_MODE_FIRST ) 
     406                        even_sum   += spec[i]; 
     407                if ( mode & EWTS_MODE_SECOND ) 
     408                        odd_sum += spec[i+1]; 
     409        } 
     410                 
     411        for (i= smoother, imod=0; i <= ndata-smoother; i+= 2, imod = (imod+2)%(smoother*2) ) { 
     412                if ( mode & EWTS_MODE_FIRST ) { 
     413                        even_sum += spec[i+smoother];           /* Add next item to window sums */ 
     414                        spec[i]   = even_sum / (smoother+1);/* Store next smoothed result */ 
     415                        even_sum -= specsm[imod];                       /* Remove oldest value from sums... */ 
     416                        specsm[imod]   = spec[i+smoother];      /*... and replace it with newest in window */ 
     417                } 
     418                if ( mode & EWTS_MODE_SECOND ) { 
     419                        odd_sum  += spec[i+smoother+1];         /* Add next item to window sums */ 
     420                        spec[i+1] = odd_sum  / (smoother+1);/* Store next smoothed result */ 
     421                        odd_sum -= specsm[imod+1];              /* Remove oldest value from sums... */ 
     422                        specsm[imod+1] = spec[i+smoother+1];/*... and replace it with newest in window */ 
     423                } 
     424        } 
     425        free(specsm); 
     426        return 0; 
     427} 
     428 
     429/*********************************************************************** 
     430 * convert_spectra () - convert complex spectra to amp/phase representation 
     431 *              return 0 if no problems, 1 otherwise 
     432 ***********************************************************************/ 
     433int convert_spectra( EW_TIME_SERIES * ewspec ) { 
     434        int i, ndata = ewspec->dataCount * 2; 
     435        double *spec = ewspec->data; 
     436        double amp, phase; 
     437        for (i= 0; i <= ndata; i+= 2) { 
    321438                amp= (double) sqrt((double) (spec[i]*spec[i]+spec[i+1]*spec[i+1])); 
    322439                phase= (double) atan2(spec[i+1], spec[i]); 
    323440                spec[i]= amp; 
    324441                spec[i+1]= phase; 
    325  
    326                 /* old scaling factor (2.0/ndata); */ 
    327                 if (scale) 
    328                         spec[i]*= timespace; 
    329  
    330                 if (white) 
    331                         sum+= spec[i]; 
    332         } 
    333         if (white) { 
    334                 mean= sum/((ndata/2.)+1.); 
    335                 for (i= 0; i <= ndata; i+= 2) 
    336                         spec[i]/= mean; 
    337         } 
    338  
    339         if (smooth) { 
    340                 int imod = 0;   /* index into smoothing window */ 
    341                 /* Create & fill smoothing window */ 
    342                 memcpy( specsm, spec, smoother*2*sizeof(double) ); 
    343                 for ( i=0; i<smoother*2; i+=2) { 
    344                         amp_sum   += spec[i]; 
    345                         phase_sum += spec[i+1]; 
    346                 } 
    347                          
    348                 for (i= smoother; i <= ndata-smoother; i+= 2) { 
    349                         /* 
    350                         amp_sum   = 0.0; 
    351                         phase_sum = 0.0; 
    352                         for (j= i-smoother; j <= i+smoother; j+= 2) { 
    353                                         amp_sum+= spec[j]; 
    354                                         phase_sum+= spec[j+1]; 
     442        } 
     443        ewspec->dataType = EWTS_TYPE_AMPPHS; 
     444        return 0; 
     445} 
     446 
     447/* Similar thing for white */ 
     448 
     449/*********************************************************************** 
     450 * calc_spectra () - calculate the spectra for a timeseries, store in  
     451 *      provided timeseries 
     452 *              return 0 if no problems, 1 otherwise 
     453 ***********************************************************************/ 
     454int calc_spectra( EW_TIME_SERIES * ewts, EW_TIME_SERIES * ewspec )  
     455{ 
     456        double *spec; 
     457        int i, ndata; 
     458        /*double delta = 1/ewts->trh2x.samprate;*/ 
     459         
     460        *ewspec = *ewts; 
     461 
     462        ndata = ewts->dataCount; 
     463 
     464        /* get nearest power of 2 for listed data length */ 
     465        for (i= 1; (int) pow((double) 2., (double) i) < ewts->dataCount; i++); 
     466        ndata= (int) pow((double) 2.,(double) (i+1)); 
     467        ewspec->dataCount = ndata/2; 
     468        ewspec->dataType = EWTS_TYPE_COMPLEX; 
     469 
     470        /* allocate space for data and smoothing area */ 
     471        if ((spec= (double *) malloc((unsigned) (ndata+2)*sizeof(double))) == NULL ) { 
     472                logit("et","%s: memory allocation error\n", progname); 
     473                return 1; 
     474        } 
     475        ewspec->data = spec; 
     476        memcpy( spec, ewts->data, ewts->dataCount*sizeof(double) ); 
     477        bzero( spec+ewts->dataCount, sizeof(double)*(ndata-ewts->dataCount+2) ); 
     478 
     479        ah_fftr(spec, ndata); 
     480         
     481        return 0; 
     482} 
     483 
     484 
     485void findmax(EW_TIME_SERIES *ewspec, int top3[]) 
     486{ 
     487        double *arr = ewspec->data; 
     488        int n = ewspec->dataCount; 
     489        double prev = arr[0]; 
     490        int i, j, cap, locap; 
     491        int goup = 1; 
     492        double *top3val; 
     493 
     494        double delta = 1/ewspec->trh2x.samprate; 
     495        double nyquist= 1.0/(2.0*delta); 
     496        double ndata = ewspec->dataCount - 1; 
     497        double scale = nyquist/ndata; 
     498         
     499        printf( "Finding %d peaks\n", numPeaks ); 
     500        top3val = malloc( sizeof(double) * numPeaks ); 
     501        if ( top3val == NULL ) { 
     502                logit( "e", "Failed to allocate memory for peak-finding\n" ); 
     503                exit(1); 
     504        } 
     505        for ( i=0; i<numPeaks; i++ ) { 
     506                top3[i] = -1; 
     507                top3val[i] = -1; 
     508        } 
     509 
     510        locap = (peakBand[0] / scale)*2; 
     511        if ( locap < 0 ) 
     512                locap = 0; 
     513        if ( locap % 2 ) 
     514                locap--; 
     515        cap = (peakBand[1] / scale)*2; 
     516        if ( cap % 2 ) 
     517                cap++; 
     518        if ( cap > n ) 
     519                cap = n; 
     520        for( i=0; i <= cap; i+=2) 
     521        { 
     522                double l,r; 
     523                l = (i - 2 < 0)?arr[0]:arr[i-2]; 
     524                r = (i + 2 >= n)?arr[n-2]:arr[i+2]; 
     525                double cur = (l + 2*arr[i] + r)/4; 
     526                if(goup) { 
     527                        if(prev > cur && i > 0) { 
     528                                int newidx = i-2; 
     529                                double newval = arr[newidx]; 
     530                                if ( newidx >= locap && newval > top3val[numPeaks-1] ) { 
     531                                        for ( j=numPeaks-2; j>=0 && newval > top3val[j]; j-- ) { 
     532                                                top3val[j+1] = top3val[j]; 
     533                                                top3[j+1] = top3[j]; 
     534                                        } 
     535                                        top3[j+1] = newidx; 
     536                                        top3val[j+1] = newval; 
     537                                } 
     538                                goup = 0; 
    355539                        } 
    356                         amp_sum/= (double) (smoother+1); 
    357                         phase_sum/= (double) (smoother+1); 
    358                         specsm[i]= amp_sum; 
    359                         specsm[i+1]= phase_sum; 
    360                         */ 
    361                         /* Add next item to window sums */ 
    362                         amp_sum   += spec[i+smoother]; 
    363                         phase_sum += spec[i+smoother+1]; 
    364                         /* Store next smoothed result */ 
    365                         spec[i] = amp_sum / (smoother+1); 
    366                         spec[i+1] = phase_sum / (smoother+1); 
    367                         /* Remove oldest value from sums... */ 
    368                         imod = (imod + 2) % (smoother*2); 
    369                         amp_sum -= specsm[imod]; 
    370                         phase_sum -= specsm[imod+1]; 
    371                         /*... and replace it with newest in window */ 
    372                         specsm[imod] = spec[i+smoother]; 
    373                         specsm[imod+1] = spec[i+smoother+1]; 
    374                 } 
    375                 /* 
    376                 for (i= smoother; i <= ndata-smoother; i++) 
    377                                 spec[i]= specsm[i]; 
    378                 */ 
    379                 free(specsm); 
    380         } 
    381  
    382         if ( fp != NULL ) { 
    383                 fprintf( fp, "%15s\t%15s\t%15s\n", "Frequency","Amplitude","Phase" ); 
    384                 for ( i=0; i <= ndata; i+=2 ) { 
    385                         amp = spec[i]; 
    386                         phase = spec[i+1]; 
    387                         fprintf( fp, "%15.3lf\t%15.3lf\t%15.3lf\n", i*nyquist/ndata, amp, phase ); 
    388                 } 
    389         } 
    390         free(spec); 
    391  
    392         return 0; 
    393 }  
    394  
    395  
     540                } else { 
     541                        if(prev < cur) 
     542                                goup = 1; 
     543                } 
     544                prev = cur; 
     545        } 
     546        free( top3val ); 
     547} 
    396548 
    397549/*********************************************************************** 
     
    400552 *              in wsh 
    401553 ***********************************************************************/ 
    402 void process_timespan(  WShandle *wsh, double start, double stop ) { 
     554void process_timespan(  WShandle *wsh, double start, double stop )  
     555{ 
    403556        int ret, i, pn; 
    404         EW_TIME_SERIES *ewts = NULL, *ewts2 = NULL; 
     557        EW_TIME_SERIES *ewts = NULL, *ewts2 = NULL, ewspec, ewspec2; 
    405558        TRACE_REQ tr; 
     559        int *top3 = malloc( sizeof(int)*numPeaks ); 
    406560 
    407561        tr.pBuf = NULL; 
     
    432586                if ( ewts == NULL ) 
    433587                        continue; /* move on to next SCNL */ 
    434  
    435                 /* Print info about request & received data (if file output) */ 
    436                 if ( fp != NULL ) { 
    437                         fprintf( fp, "Requested start=%15.3lf\tend=%15.3lf\n", start, stop ); 
    438                         fprintf( fp, "nsamp=%15ld, start=%15.3lf, end=%15.3lf, samprate=%15.3lf\n",  
    439                                 ewts->dataCount, tr.actStarttime, tr.actEndtime,  
    440                                 ewts->trh2x.samprate ); 
    441                 } 
    442                  
     588                if ( 1 ) { 
     589                        char date[2][20]; 
     590                        EWSUnConvertTime( date[0], tr.actStarttime ); 
     591                        EWSUnConvertTime( date[1], tr.actEndtime); 
     592                        EWSUnConvertTime( date[0], ewts->trh2x.starttime ); 
     593                        EWSUnConvertTime( date[1], ewts->trh2x.endtime ); 
     594                } 
     595 
    443596                /* Handle the second timeseries, if present */ 
    444597                if ( binary[pn] ) { 
     
    455608                                continue; 
    456609                        } 
     610                        if ( 1 ) { 
     611                                char date[2][20]; 
     612                                EWSUnConvertTime( date[0], tr.actStarttime ); 
     613                                EWSUnConvertTime( date[1], tr.actEndtime); 
     614                                EWSUnConvertTime( date[0], ewts2->trh2x.starttime ); 
     615                                EWSUnConvertTime( date[1], ewts2->trh2x.endtime ); 
     616                        } 
    457617                } 
    458618                 
     
    461621                tr.pBuf = NULL; 
    462622                 
    463                 /* Complete second timeseries */ 
    464                 if ( binary[pn] ) { 
     623                /* Complete second timeseries for diff */ 
     624                if ( processing_mode[pn] == 's' ) { 
    465625                        demean_EWTS( *ewts ); 
     626                        if ( taperType ) 
     627                                taper_EWTS( *ewts, taperType, taperFrac, EWTS_MODE_BOTH ); 
    466628                        demean_EWTS( *ewts2 ); 
     629                        if ( taperType ) 
     630                                taper_EWTS( *ewts2, taperType, taperFrac, EWTS_MODE_BOTH ); 
    467631                         
    468                         if ( processing_mode[pn] == 'd' ) 
    469                                 subtract_from_EWTS( *ewts, *ewts2 ); 
     632                        subtract_from_EWTS( *ewts, *ewts2, EWTS_MODE_BOTH ); 
    470633                                 
    471634                        ws2ts_purge( NULL, NULL, ewts2 ); 
    472635                }                
    473636 
    474                 /* Demean, taper, & produce spectra */ 
     637                /* Demean, filter, taper, & produce spectra */ 
    475638                demean_EWTS( *ewts );    
     639                if ( highpoles | lowpoles ) { 
     640                        if ( 0 && taperType ) 
     641                                taper_EWTS( *ewts, taperType, taperFrac, EWTS_MODE_BOTH ); 
     642                        iir( ewts, lowcutoff, lowpoles, highcutoff, highpoles ); 
     643                        demean_EWTS( *ewts ); 
     644                } 
    476645                if ( taperType ) 
    477                         taper_EWTS( *ewts, taperType, taperFrac ); 
    478                 calc_spectra( ewts ); 
    479          
     646                        taper_EWTS( *ewts, taperType, taperFrac, EWTS_MODE_BOTH ); 
     647                calc_spectra( ewts, &ewspec );   
    480648                ws2ts_purge( NULL, NULL, ewts ); 
     649                 
     650                if ( processing_mode[pn] == 'd' ) { 
     651                        demean_EWTS( *ewts2 ); 
     652                        if ( highpoles | lowpoles ) { 
     653                                if ( 0 && taperType ) 
     654                                        taper_EWTS( *ewts2, taperType, taperFrac, EWTS_MODE_BOTH ); 
     655                                iir( ewts2, lowcutoff, lowpoles, highcutoff, highpoles ); 
     656                                demean_EWTS( *ewts2 ); 
     657                        } 
     658                        if ( taperType ) 
     659                                taper_EWTS( *ewts2, taperType, taperFrac, EWTS_MODE_BOTH ); 
     660                        calc_spectra( ewts2, &ewspec2 ); 
     661                        ws2ts_purge( NULL, NULL, ewts2 ); 
     662                         
     663                        deconvolve_spectra( &ewspec, &ewspec2 ); 
     664                        if ( ewspec2.data != NULL ) 
     665                                free( ewspec2.data ); 
     666                } 
     667                 
     668                convert_spectra( &ewspec ); 
     669                if ( smooth ) 
     670                        smooth_spectra( &ewspec, EWTS_MODE_FIRST | EWTS_MODE_SECOND ); 
     671                 
     672                double delta = 1/ewspec.trh2x.samprate; 
     673                double nyquist= 1.0/(2.0*delta); 
     674                double ndata = ewspec.dataCount - 1; 
     675                double scale = nyquist/ndata; 
     676                 
     677                if ( OutRingKey != -1 ) { 
     678                        memcpy( ewspec.trh2x.sta,  SCNL[pn][2][0], TRACE2_STA_LEN ); 
     679                        memcpy( ewspec.trh2x.net,  SCNL[pn][2][1], TRACE2_NET_LEN ); 
     680                        memcpy( ewspec.trh2x.chan, SCNL[pn][2][2], TRACE2_CHAN_LEN ); 
     681                        memcpy( ewspec.trh2x.loc,  SCNL[pn][2][3], TRACE2_LOC_LEN ); 
     682                } 
     683                 
     684                MSG_LOGO logo; 
     685                logo.instid = InstId; 
     686                logo.mod    = MyModId; 
     687                 
     688                char date[2][20]; 
     689                char buffer[2000], line[100], *specData = buffer, *peakData = buffer; 
     690                EWSUnConvertTime( date[0], start ); 
     691                EWSUnConvertTime( date[1], stop ); 
     692                sprintf( peakData, "Requested start=%15s\tend=%15s\n", date[0], date[1] ); 
     693                peakData += strlen(peakData); 
     694                EWSUnConvertTime( date[0], ewspec.trh2x.starttime ); 
     695                EWSUnConvertTime( date[1], ewspec.trh2x.endtime ); 
     696                sprintf( peakData, "Actual    start=%15s\tend=%15s\n",  
     697                        date[0], date[1] ); 
     698                peakData += strlen(peakData);    
     699                sprintf( peakData, "Sampling  nsamp=%15ld\trate=%15.3lf\n",  
     700                        ewspec.dataCount, ewspec.trh2x.samprate ); 
     701                peakData += strlen(peakData); 
     702                sprintf( peakData, "Source 1: %s.%s.%s.%s",  
     703                        SCNL[pn][0][0], SCNL[pn][0][1], SCNL[pn][0][2], SCNL[pn][0][3] ); 
     704                peakData += strlen(peakData); 
     705                if ( binary[pn] ) { 
     706                        sprintf( peakData, "\tSource 2: %s.%s.%s.%s\t%s\n",  
     707                                SCNL[pn][1][0], SCNL[pn][1][1], SCNL[pn][1][2], SCNL[pn][1][3], 
     708                                processing_mode[pn] == 's' ? "(difference)" : "(deconvolution)"); 
     709                        peakData += strlen(peakData); 
     710                } else  
     711                        *peakData++ = '\n'; 
     712                if ( taperType ) 
     713                        sprintf( peakData, "Tapered %5.2lf%% (%s)\n", taperFrac*100, taperName[taperType] ); 
     714                else 
     715                        sprintf( peakData, "Tapering disabled\n" ); 
     716                peakData += strlen(peakData); 
     717                if ( smooth )  
     718                        sprintf( peakData, "Smoothed (%5.2lf Hz)\n", smooth_width ); 
     719                else 
     720                        sprintf( peakData, "Smoothing disabled\n");                      
     721                peakData += strlen(peakData); 
     722                if ( lowpoles ) 
     723                        sprintf( peakData, "LowCut  (%d poles) %lf\n", lowpoles, lowcutoff ); 
     724                else 
     725                        sprintf( peakData, "LowCut disabled\n" ); 
     726                peakData += strlen(peakData); 
     727                if ( highpoles ) 
     728                        sprintf( peakData, "HighCut (%d poles) %lf\n", highpoles, highcutoff ); 
     729                else 
     730                        sprintf( peakData, "HighCut disabled\n" ); 
     731                peakData += strlen(peakData); 
     732                specData = peakData; 
     733 
     734                if ( numPeaks > 0 ) { 
     735                        findmax( &ewspec, top3 ); 
     736                         
     737                        sprintf( peakData, "Peak band: %15.3lf thru %15.3lf\n", peakBand[0], peakBand[1] ); 
     738                        peakData += strlen(peakData); 
     739                        for ( i=0; i<numPeaks; i++ ) { 
     740                                if ( top3[i] == -1 ) 
     741                                        sprintf( peakData, "Peak %d: Not in band\n", i+1 ); 
     742                                else 
     743                                        sprintf( peakData, "Peak %d: Freq = %lf, amp = %lf\n", i+1,  
     744                                                top3[i]*scale/2, ewspec.data[top3[i]] ); 
     745                                peakData += strlen( peakData ); 
     746                        } 
     747 
     748                        if ( OutRingKey != -1 ) { 
     749                                printf("Writing peaks to ring\n"); 
     750                                GetType( "TYPE_SPECTRA_PEAKS", &(logo.type) ); 
     751                                peakData[0] = 0; 
     752                                if( tport_putmsg( &OutRegion, &logo, peakData-buffer, buffer ) != PUT_OK ) { 
     753                                   logit("et", "%s:  Error sending writing peaks.\n", 
     754                                                  Argv0 ); 
     755                                } 
     756                        }  
     757                        if ( fp == NULL && OutRingKey == -1 ) 
     758                                logit( "w", "ewspectra: Peaks requested, but not output means specified\n" ); 
     759                } 
     760  
     761                if ( outpath[0] != 0 ) { 
     762                        fp = fopen( outpath, "a" ); 
     763                        if ( fp == NULL ) 
     764                                logit( "e", "Could not open output file: '%s'\n", outpath ); 
     765                        fwrite( buffer, peakData-buffer, 1, fp ); 
     766                } 
     767                         
     768                if ( spectraOut ) {                      
     769                        printf("ReportSpectra %lf %lf %lf %lf\n", delta, nyquist, ndata, scale); 
     770                        if ( fp != NULL ) { 
     771                                printf("Writing spectra to file\n"); 
     772                                write_EWTS_as_spectra_to_file( &ewspec, Argv0, fp ); 
     773                        }  
     774                        if ( OutRingKey != -1 ) { 
     775                                write_EWTS_as_spectra_to_ring( &ewspec, Argv0, &OutRegion, &logo ); 
     776                        } 
     777                        if ( fp == NULL && OutRingKey == -1 ) 
     778                                logit( "w", "ewspectra: Spectra output requested, but not output means specified\n" ); 
     779                } 
     780                 
     781                if ( fp != NULL ) { 
     782                        fclose( fp ); 
     783                        fp = NULL; 
     784                } 
     785                         
     786                if ( ewspec.data != NULL ) 
     787                        free( ewspec.data ); 
     788                 
    481789        } 
    482790         
     
    487795} 
    488796 
     797/*************************************************************************** 
     798 * ewspectra_status() builds a heartbeat or error message & puts it into      * 
     799 *                 shared memory.  Writes errors to log file & screen.     * 
     800 ***************************************************************************/ 
     801void ewspectra_status( unsigned char type, short ierr, char *note ) 
     802{ 
     803   MSG_LOGO    logo; 
     804   char        msg[256]; 
     805   long        size; 
     806   time_t      t; 
     807 
     808/* Build the message 
     809 *******************/ 
     810   logo.instid = InstId; 
     811   logo.mod    = MyModId; 
     812   logo.type   = type; 
     813 
     814   time( &t ); 
     815 
     816   if( type == TypeHeartBeat ) 
     817    sprintf( msg, "%ld %ld\n%c", (long) t, (long) MyPid, (char)0); 
     818   else if( type == TypeError ) 
     819   { 
     820    sprintf( msg, "%ld %hd %s\n%c", (long) t, ierr, note, 0); 
     821 
     822    logit( "et", "%s(%s): %s\n", Argv0, MyModName, note ); 
     823   } 
     824 
     825   size = strlen( msg );   /* don't include the null byte in the message */ 
     826 
     827/* Write the message to shared memory 
     828 ************************************/ 
     829   if( tport_putmsg( &InRegion, &logo, size, msg ) != PUT_OK ) 
     830   { 
     831        if( type == TypeHeartBeat ) { 
     832           logit("et","%s(%s):  Error sending heartbeat.\n", 
     833                  Argv0, MyModName ); 
     834    } 
     835    else if( type == TypeError ) { 
     836           logit("et", "%s(%s):  Error sending error:%d.\n", 
     837                  Argv0, MyModName, ierr ); 
     838    } 
     839   } 
     840 
     841   return; 
     842} 
     843 
     844/********************** Message Stacking Thread ******************* 
     845 *           Move messages from transport to memory queue         * 
     846 ******************************************************************/ 
     847thr_ret MessageStacker( void *dummy ) 
     848{ 
     849   long          recsize;   /* size of retrieved message             */ 
     850   MSG_LOGO      reclogo;       /* logo of retrieved message             */ 
     851   long      filteredSize;  /* size of message after user-filtering  */ 
     852   unsigned char filteredType;  /* type of message after filtering       */ 
     853   int       res; 
     854   int       ret; 
     855   char      msg[100]; 
     856   double    start, stop; 
     857 
     858   /* Tell the main thread we're ok 
     859   ********************************/ 
     860   MessageStackerStatus = MSGSTK_ALIVE; 
     861 
     862   /* Start main export service loop for current connection 
     863   ********************************************************/ 
     864   while( 1 ) 
     865   { 
     866      /* Get a message from transport ring 
     867      ************************************/ 
     868      res = tport_getmsg( &InRegion, &GetLogo, 1, 
     869                          &reclogo, &recsize, MSrawmsg, MaxMsgSize-1 ); 
     870 
     871      /* Wait if no messages for us 
     872       ****************************/ 
     873      if( res == GET_NONE ) {sleep_ew(100); continue;} 
     874       
     875      /* Check return code; report errors 
     876      ***********************************/ 
     877      if( res != GET_OK ) 
     878      { 
     879         if( res==GET_TOOBIG ) 
     880         { 
     881            sprintf( errText, "msg[%ld] i%d m%d t%d too long for target", 
     882                            recsize, (int) reclogo.instid, 
     883                (int) reclogo.mod, (int)reclogo.type ); 
     884            ewspectra_status( TypeError, ERR_TOOBIG, errText ); 
     885            continue; 
     886         } 
     887         else if( res==GET_MISS ) 
     888         { 
     889            sprintf( errText, "missed msg(s) i%d m%d t%d in %s",(int) reclogo.instid, 
     890                    (int) reclogo.mod, (int)reclogo.type, inRing ); 
     891            ewspectra_status( TypeError, ERR_MISSMSG, errText ); 
     892         } 
     893         else if( res==GET_NOTRACK ) 
     894         { 
     895            sprintf( errText, "no tracking for logo i%d m%d t%d in %s", 
     896                     (int) reclogo.instid, (int) reclogo.mod, (int)reclogo.type, 
     897                     inRing ); 
     898            ewspectra_status( TypeError, ERR_NOTRACK, errText ); 
     899         } 
     900      } 
     901 
     902      /* Process retrieved msg (res==GET_OK,GET_MISS,GET_NOTRACK) 
     903      ***********************************************************/ 
     904          strncpy( msg, MSrawmsg, recsize ); 
     905          msg[recsize] = 0; 
     906                start = atof( msg ); 
     907                if ( start < 0 ) {  
     908                        start += time(NULL); 
     909                } else if ( EWSConvertTime (msg, &start) == EW_FAILURE ) { 
     910            sprintf( errText, "COMPUTE_SPECTRA message w/ bad time" ); 
     911            ewspectra_status( TypeError, ERR_BADMSG, errText ); 
     912                } 
     913                stop = start + atoi( msg+14 ); 
     914           
     915                process_timespan( wsh, start, stop ); 
     916 
     917   } /* end of while */ 
     918 
     919   /* we're quitting 
     920   *****************/ 
     921error: 
     922   MessageStackerStatus = MSGSTK_ERR; /* file a complaint to the main thread */ 
     923   KillSelfThread(); /* main thread will restart us */ 
     924} 
     925 
    489926int main(int argc, char **argv) 
    490927{ 
    491         WShandle *wsh; 
     928/* Other variables: */ 
     929   int           res; 
     930   long          recsize;   /* size of retrieved message             */ 
     931   MSG_LOGO      reclogo;   /* logo of retrieved message             */ 
     932         
     933        /* Look up installations of interest 
     934        *********************************/ 
     935        if ( GetLocalInst( &InstId ) != 0 ) { 
     936          fprintf( stderr, 
     937                          "%s: error getting local installation id; exiting!\n", 
     938                           Argv0 ); 
     939          exit( -1 ); 
     940        } 
     941   if ( GetType( "TYPE_HEARTBEAT", &TypeHeartBeat ) != 0 ) { 
     942      fprintf( stderr, 
     943              "%s: Invalid message type <TYPE_HEARTBEAT>; exiting!\n", Argv0 ); 
     944      exit( -1 ); 
     945   } 
     946   if ( GetType( "TYPE_ERROR", &TypeError ) != 0 ) { 
     947      fprintf( stderr, 
     948              "%s: Invalid message type <TYPE_ERROR>; exiting!\n", Argv0 ); 
     949      exit( -1 ); 
     950   } 
    492951 
    493952        if ((progname= strrchr(*argv, (int) '/')) != (char *) 0) 
     
    504963    } 
    505964         
     965        /* Get our own Pid for restart purposes 
     966        ***************************************/ 
     967        MyPid = getpid(); 
     968        if(MyPid == -1) 
     969        { 
     970      logit("e", "%s: Cannot get pid; exiting!\n", Argv0); 
     971      return(0); 
     972        } 
     973 
    506974        /* Read-in & interpret config file */ 
    507975        scale= white= smooth= taperType= 0; 
    508976        num_ews_procs = 0; 
    509977        outpath[0] = 0; 
     978        find_triggers = 0; 
     979        lowpoles = highpoles = lowcutoff = highcutoff = 0; 
     980        InRingKey = OutRingKey = -1; 
     981        start = -1; 
    510982        wsh = init_ws( argv[1], progname, parse_my_command ); 
    511983         
     
    519991                exit(1); 
    520992        } 
    521          
    522         if ( outpath[0] != 0 ) { 
    523                 fp = fopen( outpath, "w" ); 
    524                 if ( fp == NULL ) 
    525                         logit( "e", "Could not open output file: '%s'\n", outpath ); 
    526         } 
    527          
    528         if ( start ) 
     993        if ( OutRingKey == -1 && outpath[0] == 0 ) { 
     994                logit( "e", "No output specified; exiting\n" ); 
     995                ws2ts_purge( wsh, NULL, NULL ); 
     996                exit(1); 
     997        } 
     998        if ( lowcutoff > 0 && highcutoff > 0 && lowcutoff >= highcutoff ) { 
     999                logit( "e", "Low Cutoff Frequency Must Be < High Cutoff Frequency.\n" ); 
     1000                ws2ts_purge( wsh, NULL, NULL ); 
     1001                exit(1); 
     1002        } 
     1003        if ( OutRingKey != -1 ) { 
     1004                int i; 
     1005                for ( i=num_ews_procs-1; i>=0; i-- ) 
     1006                        if ( *SCNL[i][2][3] == 0 ) 
     1007                                if ( outpath[0] != 0 ) 
     1008                                        logit( "w", "No SCNL for ring messages specified for %s.%s.%s.%s; will only write to file\n", 
     1009                                                SCNL[i][0][0], SCNL[i][0][1], SCNL[i][0][2], SCNL[i][0][3] ); 
     1010                                else { 
     1011                                        logit( "w", "No SCNL for ring messages specified for %s.%s.%s.%s; aborting this SCNL\n", 
     1012                                                SCNL[i][0][0], SCNL[i][0][1], SCNL[i][0][2], SCNL[i][0][3] ); 
     1013                                        num_ews_procs--; 
     1014                                        if ( i < num_ews_procs ) { 
     1015                                                int j,k; 
     1016                                                for ( j=0; j<3; j+=(binary[i]?1:2) ) 
     1017                                                        for ( k=0; j<4; k++ ) 
     1018                                                                strcpy( SCNL[i][j][k], SCNL[num_ews_procs][j][k] ); 
     1019                                                binary[i] = binary[num_ews_procs]; 
     1020                                                processing_mode[i] = processing_mode[num_ews_procs]; 
     1021                                        } 
     1022                                } 
     1023                if ( num_ews_procs == 0 ) { 
     1024                        logit( "e", "All SCNLs aborted; exiting\n" ); 
     1025                        ws2ts_purge( wsh, NULL, NULL ); 
     1026                        exit(1); 
     1027                } 
     1028                 
     1029        } 
     1030        if ( find_triggers ) { 
     1031                int geti = GetInst( "INST_WILDCARD", &(GetLogo.instid) ); 
     1032                int getm = GetModId( "MOD_WILDCARD", &(GetLogo.mod) ); 
     1033                int gett = GetType( "TYPE_COMPUTE_SPECTRA", &(GetLogo.type) ); 
     1034                int getMsg = 0; 
     1035                if ( ( MSrawmsg = (char *) malloc(MaxMsgSize) ) ==  NULL ) 
     1036                { 
     1037                  logit( "e", "%s: error allocating MSrawmsg; exiting!\n", 
     1038                                 Argv0 ); 
     1039                  getMsg = 1; 
     1040                } 
     1041                if ( geti || getm || gett || getMsg ) { 
     1042                        if ( geti ) 
     1043                                logit( "e", "%s: INST_WILDCARD unknown; exiting!\n", Argv0 ); 
     1044                        if ( getm ) 
     1045                                logit( "e", "%s: MOD_WILDCARD unknown; exiting!\n", Argv0 ); 
     1046                        if ( geti ) 
     1047                                logit( "e", "%s: TYPE_COMPUTE_SPECTRA unknown; exiting!\n", Argv0 ); 
     1048                        ws2ts_purge( wsh, NULL, NULL ); 
     1049                        exit(1); 
     1050                } 
     1051                tport_attach( &InRegion, InRingKey );            
     1052        } 
     1053        if ( OutRingKey != -1 ) { 
     1054                tport_attach( &OutRegion, OutRingKey ); 
     1055        } 
     1056         
     1057        if ( start != -1 ) 
    5291058                process_timespan( wsh, start, stop ); 
    530          
     1059 
    5311060        /* Here's where, if InputRing is defined, we'd listen for  
    5321061                COMPUTE_SPECTRA messages, and call process_timespan 
    5331062                for their timespans */ 
     1063        if ( find_triggers ) { 
     1064 
     1065           /* step over all messages from transport ring 
     1066           *********************************************/ 
     1067           /* As Lynn pointed out: if we're restarted by startstop after hanging, 
     1068                  we should throw away any of our messages in the transport ring. 
     1069                  Else we could end up re-sending a previously sent message, causing 
     1070                  time to go backwards... */ 
     1071           do 
     1072           { 
     1073                 res = tport_getmsg( &InRegion, &GetLogo, 1, 
     1074                                                         &reclogo, &recsize, MSrawmsg, MaxMsgSize-1 ); 
     1075           } while (res !=GET_NONE); 
     1076         
     1077           /* One heartbeat to announce ourselves to statmgr 
     1078           ************************************************/ 
     1079           ewspectra_status( TypeHeartBeat, 0, "" ); 
     1080           time(&MyLastBeat); 
     1081         
     1082         
     1083           /* Start the message stacking thread if it isn't already running. 
     1084                ****************************************************************/ 
     1085           if (MessageStackerStatus != MSGSTK_ALIVE ) 
     1086           { 
     1087                 if ( StartThread(  MessageStacker, (unsigned)THREAD_STACK, &tidStacker ) == -1 ) 
     1088                 { 
     1089                   logit( "e", 
     1090                                  "%s(%s): Error starting  MessageStacker thread; exiting!\n", 
     1091                          Argv0, MyModName ); 
     1092                   tport_detach( &InRegion ); 
     1093                   tport_detach( &OutRegion ); 
     1094                   return( -1 ); 
     1095                 } 
     1096                 MessageStackerStatus = MSGSTK_ALIVE; 
     1097           } 
     1098         
     1099         
     1100           /* Start main ringdup service loop 
     1101           **********************************/ 
     1102           while( tport_getflag( &InRegion ) != TERMINATE  && 
     1103                          tport_getflag( &InRegion ) != MyPid         ) 
     1104           { 
     1105                 /* Beat the heart into the transport ring 
     1106                  ****************************************/ 
     1107                  time(&now); 
     1108                  if (difftime(now,MyLastBeat) > (double)HeartBeatInt ) 
     1109                  { 
     1110                          ewspectra_status( TypeHeartBeat, 0, "" ); 
     1111                  time(&MyLastBeat); 
     1112                  } 
     1113         
     1114                  /* take a brief nap; added 970624:ldd 
     1115                   ************************************/ 
     1116                  sleep_ew(500); 
     1117           } /*end while of monitoring loop */ 
     1118                 
     1119        } 
    5341120                 
    5351121        /* Clean up after ourselves */ 
     1122        if ( InRingKey != -1 ) 
     1123                tport_detach( &InRegion ); 
     1124        if ( OutRingKey != -1 ) 
     1125                tport_detach( &OutRegion ); 
    5361126        ws2ts_purge( wsh, NULL, NULL ); 
    537         if ( fp != NULL ) 
    538                 fclose( fp ); 
    539 } 
     1127} 
  • trunk/src/seismic_processing/ewspectra/makefile.ux

    r4184 r4194  
    44LIBDIR = $(EW_HOME)/$(EW_VERSION)/lib 
    55 
    6 CFLAGS = $(GLOBALFLAGS) -I. 
     6CFLAGS = $(GLOBALFLAGS) -I. -g 
    77LDFLAGS = -lpthread -lc -lm 
    88 
    9 SRCS = ewspectra.c 
     9SRCS = ewspectra.c iir.c 
     10OBJS = $(SRCS:%.c=%.o) 
    1011 
    11 OBJS = $(SRCS:%.c=%.o) 
     12SRCSc = compute_spectra.c 
     13OBJSc = $(SRCSc:%.c=%.o) 
     14 
     15SRCSs = sniffspectra.c 
     16OBJSs = $(SRCSs:%.c=%.o) 
    1217 
    1318EW_LIBS = $(LIBDIR)/logit_mt.o $(LIBDIR)/kom.o $(LIBDIR)/threads_ew.o \ 
    1419          $(LIBDIR)/time_ew.o $(LIBDIR)/transport.o $(LIBDIR)/sleep_ew.o \ 
    1520          $(LIBDIR)/getutil.o $(LIBDIR)/sema_ew.o $(LIBDIR)/ws_clientII.o \ 
    16           $(LIBDIR)/socket_ew_common.o $(LIBDIR)/socket_ew.o \ 
    17           $(LIBDIR)/swap.o $(LIBDIR)/ws2ts.o 
     21          $(LIBDIR)/socket_ew_common.o $(LIBDIR)/socket_ew.o $(LIBDIR)/swap.o \ 
     22          $(LIBDIR)/ew_spectra_io.o 
    1823  
    1924.c.o: 
    2025        $(CC) $(CFLAGS) -c $< -o $@ 
    2126 
    22 all: ewspectra 
     27all: ewspectra compute_spectra sniffspectra 
    2328 
    24 ewspectra: $(OBJS) $(LIBDIR)/ws2ts.o  
    25         $(CC) $(GLOBALFLAGS) -o ewspectra $(OBJS) $(EW_LIBS) $(LDFLAGS) 
     29ewspectra: $(OBJS) 
     30        $(CC) $(GLOBALFLAGS) -o ewspectra $(OBJS) $(LIBDIR)/ws2ts.o $(EW_LIBS) $(LDFLAGS) 
    2631        cp ewspectra $(BINDIR) 
     32 
     33compute_spectra: $(OBJSc) 
     34        $(CC) $(GLOBALFLAGS) -o compute_spectra $(OBJSc) $(EW_LIBS) $(LDFLAGS) 
     35        cp compute_spectra $(BINDIR) 
     36 
     37sniffspectra: $(OBJSs) 
     38        $(CC) $(GLOBALFLAGS) -o sniffspectra $(OBJSs) $(EW_LIBS) $(LDFLAGS) 
     39        cp sniffspectra $(BINDIR) 
    2740 
    2841clean: 
    2942        rm *.o 
    30         rm ewspectra 
     43        rm ewspectra compute_spectra sniffspectra 
Note: See TracChangeset for help on using the changeset viewer.