sffxx.cc 13.9 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
1
2
3
4
5
/*! \file sffxx.cc
 * \brief SFF library (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
6
 * $Id: sffxx.cc,v 1.9 2004-04-07 13:22:49 tforb Exp $
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
13
14
15
 * \author Thomas Forbriger
 * \date 21/12/2003
 * 
 * SFF library (implementation)
 * 
 * Copyright (c) 2003 by Thomas Forbriger (BFO Schiltach) 
 * 
 * REVISIONS and CHANGES 
 *  - 21/12/2003   V1.0   Thomas Forbriger
16
17
 *  - 23/02/2004   V1.1   changed DAST and STAT reading code
 *                        there are not necessarily any flag characters
18
19
20
 *  - 07/04/2004   V1.2  
 *                        - provide debug output
 *                        - correct reading of FREE block
thomas.forbriger's avatar
thomas.forbriger committed
21
22
23
24
 * 
 * ============================================================================
 */
#define TF_SFFXX_CC_VERSION \
25
  "TF_SFFXX_CC   V1.2"
thomas.forbriger's avatar
thomas.forbriger committed
26
#define TF_SFFXX_CC_CVSID \
27
  "$Id: sffxx.cc,v 1.9 2004-04-07 13:22:49 tforb Exp $"
thomas.forbriger's avatar
thomas.forbriger committed
28

29
#include<sstream>
thomas.forbriger's avatar
thomas.forbriger committed
30
#include <sffxx.h>
thomas.forbriger's avatar
thomas.forbriger committed
31
#include <gsexx.h>
thomas.forbriger's avatar
thomas.forbriger committed
32
33
34

namespace sff {

35
36
37
38
39
40
41
  namespace helper {
    //! Check GSE identifier at beginning of line.
    template<class C>
    bool IDmatch(const std::string& line)
    { return(line.substr(0,4)==std::string(C::LINEID)); }
  } // namespace helper

thomas.forbriger's avatar
thomas.forbriger committed
42
  //! Fortran library version (to ensure compatibility)
thomas.forbriger's avatar
thomas.forbriger committed
43
  const double STAT::libversion=1.10;
thomas.forbriger's avatar
thomas.forbriger committed
44
45
46
47
48
49
  const char* const STAT::LINEID="STAT";
  const char* const FREE::LINEID="FREE";
  const char* const SRCE::LINEID="SRCE";
  const char* const DAST::LINEID="DAST";
  const char* const INFO::LINEID="INFO";

thomas.forbriger's avatar
thomas.forbriger committed
50
51
  /*----------------------------------------------------------------------*/

thomas.forbriger's avatar
thomas.forbriger committed
52
53
54
55
56
57
  char coosysID(const Ecoosys& csid)
  {
    char retval;
    if (csid == CS_cartesian) retval='C';
    else if (csid == CS_spherical) retval='S';
    else throw
thomas.forbriger's avatar
thomas.forbriger committed
58
      GSE2::Terror("ERROR (sff::coosysID): library inconsistency!");
thomas.forbriger's avatar
thomas.forbriger committed
59
60
61
62
63
    return(retval);
  } // coosysID

  Ecoosys coosysID(const char& csid)
  {
thomas.forbriger's avatar
thomas.forbriger committed
64
65
66
    Ecoosys retval;
    if (csid=='C') { retval=CS_cartesian; }
    else if (csid=='S') { retval=CS_spherical; }
thomas.forbriger's avatar
thomas.forbriger committed
67
    else throw
thomas.forbriger's avatar
thomas.forbriger committed
68
69
      GSE2::Terror("ERROR (sff::coosysID): unknown coordinate system key!");
    return(retval);
thomas.forbriger's avatar
thomas.forbriger committed
70
71
  } // coosysID

thomas.forbriger's avatar
thomas.forbriger committed
72
73
74
75
76
77
78
/*======================================================================*/
// SFF structs
// -----------
//
// STAT
// ----

79
80
81
82
83
84
85
  /*! \struct STAT
   *
   * The STAT line is the first line in the file header.
   * It contains the version of the library that wrote the file,
   * a timestamp and flags indicating the presence of optional elements lika a
   * FREE block or an SRCE line.
   */
thomas.forbriger's avatar
thomas.forbriger committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  STAT::STAT(): hasfree(false), hassrce(false) 
  { setstamp(libtime::now()); }

  void STAT::setstamp(const libtime::TAbsoluteTime& date) const
  {
    char stamp[14];
    int yeardigits, century;
    century=int(date.year()/100);
    yeardigits=date.year()-100*century;
    sprintf(stamp, "%2.2i%2.2li%2.2li.%2.2li%2.2li%2.2li",
            yeardigits, date.month(), date.day(),
            date.hour(), date.minute(), date.second());
    timestamp=std::string(stamp);
  }

  std::string STAT::line() const
  {
    this->setstamp(libtime::now());
    char charline[40];
    std::string code("");
    if (this->hasfree) { code.append("F"); }
    if (this->hassrce) { code.append("S"); }
thomas.forbriger's avatar
thomas.forbriger committed
108
    sprintf(charline, "%-4s  %6.2f %-13s %-10s\n",
thomas.forbriger's avatar
thomas.forbriger committed
109
110
111
112
113
114
115
116
            STAT::LINEID, 
            STAT::libversion,
            timestamp.c_str(),
            code.c_str());
    std::string retval(charline);
    return(retval);
  } // STAT::line()

117
  void STAT::read(std::istream& fis)
118
  {
119
120
121
122
123
124
    std::string theline("");
    const int bufsize=81;
    char inputline[bufsize];
    fis.getline(inputline,bufsize);
    theline=inputline;
    std::istringstream is(inputline);
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
    std::string lineID;
    is >> lineID;
    if (!helper::IDmatch<STAT>(lineID)) throw
       GSE2::Terror("ERROR (STAT::read): missing STAT ID!");

    double inlibversion;
    is >> inlibversion;
    if (inlibversion>STAT::libversion) 
    { 
      throw 
        GSE2::Terror("ERROR (STAT::read): file library version too large!"); 
    }

    is >> timestamp;

    std::string code;
    is >> code;
    this->hasfree=(code.find('F')!=std::string::npos);
    this->hassrce=(code.find('S')!=std::string::npos);
  } // STAT::read

thomas.forbriger's avatar
thomas.forbriger committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
/*----------------------------------------------------------------------*/
// SRCE
// ----

  SRCE::SRCE(): 
    type("NSP"), date(libtime::now()),
    cs(CS_cartesian), cx(0.), cy(0.), cz(0.) { }

  std::string SRCE::line() const
  {
    char charline[95];
    int yeardigits, century;
    century=int(date.year()/100);
    yeardigits=date.year()-100*century;
thomas.forbriger's avatar
thomas.forbriger committed
160
    sprintf(charline, "%-4s %-20s %1c %15.6f%15.6f%15.6f "
thomas.forbriger's avatar
thomas.forbriger committed
161
162
163
164
165
166
167
168
169
170
            "%2.2i%2.2li%2.2li %2.2li%2.2li%2.2li.%3.3li\n",
            SRCE::LINEID, 
            type.c_str(), 
            coosysID(cs), cx, cy, cz,
            yeardigits, date.month(), date.day(),
            date.hour(), date.minute(), date.second(), date.milsec());
    std::string retval(charline);
    return(retval);
  } // SRCE::line()

171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
  void SRCE::read(std::istream& is)
  {
    std::string lineID;
    is >> lineID;
    if (!helper::IDmatch<SRCE>(lineID)) throw
       GSE2::Terror("ERROR (SRCE::read): missing SRCE ID!");

    char intype[21];
    is.get(intype, 21);
    type=intype;

    char cschar;
    is >> cschar;
    cs=coosysID(cschar);
    is >> cx;
    is >> cy;
    is >> cz;

    std::string datestring,timestring;
    is >> datestring;
    is >> timestring;
    std::string fulldate("");
193
    fulldate+=datestring.substr(0,2);
194
    fulldate+="/";
195
    fulldate+=datestring.substr(2,2);
196
    fulldate+="/";
197
    fulldate+=datestring.substr(4,2);
198
    fulldate+=" ";
199
    fulldate+=timestring.substr(0,2);
200
    fulldate+=":";
201
    fulldate+=timestring.substr(2,2);
202
    fulldate+=":";
203
    fulldate+=timestring.substr(4,6);
204
205
206
    date=libtime::TAbsoluteTime(fulldate);
  } // SRCE::read

thomas.forbriger's avatar
thomas.forbriger committed
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
/*----------------------------------------------------------------------*/
// DAST
// ----

  DAST::DAST(): 
    nchar(-1), ampfac(1.), hasfree(false), 
    hasinfo(false), last(false) { }

  std::string DAST::line() const
  {
    char charline[50];
    std::string code("");
    if (this->hasfree) { code.append("F"); }
    if (this->hasinfo) { code.append("I"); }
    if (!this->last) { code.append("D"); }
222
    // write -1 to nchar field in any case
thomas.forbriger's avatar
thomas.forbriger committed
223
    sprintf(charline, "%-4s  %10i %16.6E %-10s\n",
thomas.forbriger's avatar
thomas.forbriger committed
224
            DAST::LINEID, 
225
            -1, ampfac, code.c_str());
thomas.forbriger's avatar
thomas.forbriger committed
226
227
228
229
    std::string retval(charline);
    return(retval);
  } // DAST::line()

230
  void DAST::read(std::istream& fis, const bool& debug)
231
  {
232
233
234
235
236
237
    std::string theline("");
    const int bufsize=81;
    char inputline[bufsize];
    fis.getline(inputline,bufsize);
    theline=inputline;
    std::istringstream is(inputline);
238
239
240
241
242
243
244
245
    std::string lineID;
    is >> lineID;
    if (!helper::IDmatch<DAST>(lineID)) throw
       GSE2::Terror("ERROR (DAST::read): missing DAST ID!");
    is >> nchar;
    is >> ampfac;
    std::string code;
    is >> code;
246
247
248
249
250
    if (debug)
    {
      std::cerr << "DEBUG (DAST): read nchar=" << nchar
        << " ampfac=" << ampfac << " code=" << code << std::endl;
    }
251
252
253
254
255
    this->hasinfo=(code.find('I')!=std::string::npos);
    this->hasfree=(code.find('F')!=std::string::npos);
    this->last=(code.find('D')==std::string::npos);
  } // DAST::read

thomas.forbriger's avatar
thomas.forbriger committed
256
257
258
259
260
261
262
263
/*----------------------------------------------------------------------*/
// FREE
// ----

    FREE::FREE() { lines.clear(); }

    void FREE::write(std::ostream& os) const
    {
thomas.forbriger's avatar
thomas.forbriger committed
264
      os << FREE::LINEID << " " << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
265
266
267
      for(Tlines::const_iterator I=lines.begin();
          I != lines.end(); I++)
      { os << *I << std::endl; }
thomas.forbriger's avatar
thomas.forbriger committed
268
      os << FREE::LINEID << " " << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
269
270
    } // FREE::write

271
  void FREE::read(std::istream& is, const bool& debug)
272
273
  {
    std::string lineID;
274
    // getline(is,lineID);
275
    is >> lineID;
276
277
278
    if (debug) 
    { std::cerr << "DEBUG (FREE): lineID=" << lineID << std::endl; }
    if (!helper::IDmatch<FREE>(lineID.substr(0,4))) throw
279
       GSE2::Terror("ERROR (FREE::read): missing FREE ID!");
280
    lines.clear();
281
282
283
    is.ignore(10,'\n');
    getline(is,lineID);
    while (!helper::IDmatch<FREE>(lineID.substr(0,4)))
284
    {
285
286
      lines.push_back(lineID);
      getline(is,lineID);
287
288
289
    }
  } // FREE::read

thomas.forbriger's avatar
thomas.forbriger committed
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
/*----------------------------------------------------------------------*/
// WID2
// ----

  WID2::WID2():
    date(libtime::now()), station("NSP"), channel("NSP"), auxid("NSP"),
    nsamples(-1), dt(-1.), calib(-1.), calper(-1.), instype("NSP"),
    hang(-1.),vang(-1.) { }

  std::string WID2::line() const
  {
    GSE2::waveform::TWID2 wid2line;
    wid2line.Fyear=date.year();
    wid2line.Fmonth=date.month();
    wid2line.Fday=date.day();
    wid2line.Fhour=date.hour();
    wid2line.Fminute=date.minute();
    wid2line.Fseconds=double(date.second())+
      1.e-3*double(date.milsec()+1.e-3*double(date.micsec()));
    wid2line.Fstation=station;
    wid2line.Fchannel=channel;
    wid2line.Fauxid=auxid;
    wid2line.Fsamps=nsamples;
    wid2line.Fsamprate=1./dt;
    wid2line.Fcalib=calib;
    wid2line.Fcalper=calper;
    wid2line.Finstype=instype;
    wid2line.Fhang=hang;
    wid2line.Fvang=vang;
    return(wid2line.line());
  } // WID2::line()

322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
  void WID2::read(std::istream& is)
  {
    GSE2::waveform::TWID2 wid2line;
    wid2line.read(is);
    int second=int(wid2line.Fseconds);
    int milsec=int(1000*(wid2line.Fseconds-double(second)));
    date=libtime::TAbsoluteTime(wid2line.Fyear,
                                wid2line.Fmonth,
                                wid2line.Fday,
                                wid2line.Fhour,
                                wid2line.Fminute,
                                second, milsec);
    this->station=wid2line.Fstation;
    this->channel=wid2line.Fchannel;
    this->auxid=wid2line.Fauxid;
    this->nsamples=wid2line.Fsamps;
    this->dt=1./wid2line.Fsamprate;
    this->calib=wid2line.Fcalib;
    this->calper=wid2line.Fcalper;
    this->instype=wid2line.Finstype;
    this->hang=wid2line.Fhang;
    this->vang=wid2line.Fvang;
  } // WID2::read

thomas.forbriger's avatar
thomas.forbriger committed
346
347
348
349
350
351
352
353
354
355
/*----------------------------------------------------------------------*/
// INFO
// ----

  INFO::INFO():
    cs(CS_cartesian), cx(0.), cy(0.), cz(0.), nstacks(0) { }

  std::string INFO::line() const
  {
    char charline[60];
thomas.forbriger's avatar
thomas.forbriger committed
356
    sprintf(charline, "%-4s %1c %15.6f%15.6f%15.6f %4i\n",
thomas.forbriger's avatar
thomas.forbriger committed
357
358
359
360
361
362
            INFO::LINEID, 
            coosysID(cs), cx, cy, cz, nstacks);
    std::string retval(charline);
    return(retval);
  } // INFO::line()

363
364
365
366
  void INFO::read(std::istream& is)
  {
    std::string lineID;
    is >> lineID;
367
368
    if (!helper::IDmatch<INFO>(lineID)) throw
       GSE2::Terror("ERROR (INFO::read): missing INFO ID!");
369
370
371
372
373
374
375
376
377
    char cschar;
    is >> cschar;
    cs=coosysID(cschar);
    is >> cx;
    is >> cy;
    is >> cz;
    is >> nstacks;
  } // INFO::read

thomas.forbriger's avatar
thomas.forbriger committed
378
379
380
381
382
383
384
385
386
387
388
/*----------------------------------------------------------------------*/
// FileHeader
// ----------

  void FileHeader::write(std::ostream& os) const
  {
    os << Mstat.line();
    if (Mstat.hasfree) { Mfree.write(os); }
    if (Mstat.hassrce) { os << Msrce.line(); }
  }

389
390
391
392
393
394
395
  void FileHeader::read(std::istream& is)
  {
    Mstat.read(is);
    if (Mstat.hasfree) { Mfree.read(is); }
    if (Mstat.hassrce) { Msrce.read(is); }
  }

thomas.forbriger's avatar
thomas.forbriger committed
396
397
398
399
/*----------------------------------------------------------------------*/
// TraceHeader
// -----------

thomas.forbriger's avatar
thomas.forbriger committed
400
  void TraceHeader::writeheader(std::ostream& os) const
thomas.forbriger's avatar
thomas.forbriger committed
401
  {
402
    if (Mdebug) { std::cerr << "DEBUG: write DAST line" << std::endl; } 
thomas.forbriger's avatar
thomas.forbriger committed
403
    os << Mdast.line();
404
    if (Mdebug) { std::cerr << "DEBUG: write WID2 line" << std::endl; } 
thomas.forbriger's avatar
thomas.forbriger committed
405
    os << Mwid2.line();
thomas.forbriger's avatar
thomas.forbriger committed
406
407
408
409
  }

  void TraceHeader::writetrailer(std::ostream& os) const
  {
410
411
412
413
414
415
416
417
418
419
    if (Mdast.hasfree) 
    {
      if (Mdebug) { std::cerr << "DEBUG: write FREE block" << std::endl; } 
      Mfree.write(os);
    }
    if (Mdast.hasinfo) 
    {
      if (Mdebug) { std::cerr << "DEBUG: write INFO line" << std::endl; } 
      os << Minfo.line();
    }
thomas.forbriger's avatar
thomas.forbriger committed
420
421
  }

422
423
  void TraceHeader::readheader(std::istream& is) 
  {
424
425
426
    if (Mdebug) { std::cerr << "DEBUG: read DAST line" << std::endl; } 
    Mdast.read(is,Mdebug);
    if (Mdebug) { std::cerr << "DEBUG: read WID2 line" << std::endl; } 
427
428
429
430
431
    Mwid2.read(is);
  }

  void TraceHeader::readtrailer(std::istream& is) 
  {
432
433
434
435
436
437
438
439
440
441
442
443
    if (Mdebug) { std::cerr << "DEBUG: read trace trailer" << std::endl; } 
    if (Mdast.hasfree) 
    {
      if (Mdebug) { std::cerr << "DEBUG: read FREE block" << std::endl; } 
      Mfree.read(is, Mdebug); 
      if (Mdebug) { Mfree.write(std::cerr); }
    }
    if (Mdast.hasinfo) 
    {
      if (Mdebug) { std::cerr << "DEBUG: read INFO line" << std::endl; } 
      Minfo.read(is);
    }
444
445
  }

thomas.forbriger's avatar
thomas.forbriger committed
446
447
448
449
450
451
452
453
/*----------------------------------------------------------------------*/
// WaveformNormalizer
// ------------------

  const int WaveformNormalizer::limit=0x800000;

  WaveformNormalizer::WaveformNormalizer(const Enormmode& nm, 
                                         const double& maxval):
thomas.forbriger's avatar
thomas.forbriger committed
454
    Mmaxval(maxval), Mnorm(nm)
thomas.forbriger's avatar
thomas.forbriger committed
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
  {
    if (Mnorm == NM_one) 
    { 
      Mampfac=1.; 
      Mscale=false;
      if (Mmaxval > double(WaveformNormalizer::limit)) throw
        GSE2::Terror("ERROR (sff::WaveformNormalizer::scan): "
                     "dynamic range to large for non-normalizing mode");
    }
    else if (Mnorm == NM_ifneeded)
    { 
      Mampfac=1.;
      Mscale=false;
      if (Mmaxval > double(WaveformNormalizer::limit))
      {
        Mampfac=Mmaxval/double(WaveformNormalizer::limit); 
        Mscale=true;
      }
    }
    else if (Mnorm == NM_maxdyn)
    { 
      Mampfac=Mmaxval/double(WaveformNormalizer::limit); 
      Mscale=true;
    }
    else throw
      GSE2::Terror("ERROR (sff::WaveformNormalizer::scan): "
                   "library inconsistency!");
  }

484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
/*----------------------------------------------------------------------*/
// SkipWaveform
// ------------

    void SkipWaveform::read(std::istream& is) 
    {
      Mheader.readheader(is);;
      int nsamples=Mheader.wid2().nsamples;
      GSE2::waveform::TDAT2readCM6 freader(nsamples);
      int idata;
      for (int i=0; i<nsamples; i++)
      { idata=freader(is); }
      Mheader.readtrailer(is);
      Mvalid=true;
    } // SkipWaveform::read

thomas.forbriger's avatar
thomas.forbriger committed
500
501
502
} // namespace sff

/* ----- END OF sffxx.cc ----- */