polymodel.cc 12.7 KB
Newer Older
thomas.forbriger's avatar
thomas.forbriger committed
1
2
3
4
5
/*! \file polymodel.cc
 * \brief polynomial subsurface model (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
6
 * $Id: polymodel.cc,v 1.7 2003-01-03 19:38:05 forbrig Exp $
thomas.forbriger's avatar
thomas.forbriger committed
7
8
9
10
11
12
13
14
15
 * \author Thomas Forbriger
 * \date 30/12/2002
 * 
 * polynomial subsurface model (implementation)
 * 
 * Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt) 
 * 
 * REVISIONS and CHANGES 
 *  - 30/12/2002   V1.0   Thomas Forbriger
thomas.forbriger's avatar
thomas.forbriger committed
16
17
18
 *  - 03/01/2003   V1.1   
 *                        - version with means compiles
 *                        - reinvented value function
thomas.forbriger's avatar
thomas.forbriger committed
19
20
21
22
 * 
 * ============================================================================
 */
#define TF_POLYMODEL_CC_VERSION \
thomas.forbriger's avatar
thomas.forbriger committed
23
  "TF_POLYMODEL_CC   V1.1"
thomas.forbriger's avatar
thomas.forbriger committed
24
#define TF_POLYMODEL_CC_CVSID \
25
  "$Id: polymodel.cc,v 1.7 2003-01-03 19:38:05 forbrig Exp $"
thomas.forbriger's avatar
thomas.forbriger committed
26
27
28

#include <gremlin1/polymodel.h>
#include <tfxx/error.h>
thomas.forbriger's avatar
thomas.forbriger committed
29
#include <aff/shaper.h>
thomas.forbriger's avatar
thomas.forbriger committed
30
#include <aff/slice.h>
thomas.forbriger's avatar
thomas.forbriger committed
31
32
#include <stdio.h>
#include <string>
33
#include <fstream>
thomas.forbriger's avatar
thomas.forbriger committed
34
35
36

namespace gremlin1 {

thomas.forbriger's avatar
thomas.forbriger committed
37
  typedef PolynomialModelFile::Tvalue Tvalue;
thomas.forbriger's avatar
thomas.forbriger committed
38
  typedef PolynomialModelFile::Tarray Tarray;
thomas.forbriger's avatar
thomas.forbriger committed
39
40
41

  //! define identifier
  const char PolynomialModelFile::modversion2[];
thomas.forbriger's avatar
thomas.forbriger committed
42
  //! define static
thomas.forbriger's avatar
thomas.forbriger committed
43
  const int PolynomialModelFile::norder;
thomas.forbriger's avatar
thomas.forbriger committed
44
  //! define static
thomas.forbriger's avatar
thomas.forbriger committed
45
  const int PolynomialModelFile::nparameters;
thomas.forbriger's avatar
thomas.forbriger committed
46

thomas.forbriger's avatar
thomas.forbriger committed
47
48
  //! default constructor
  PolynomialModelFile::PolynomialModelFile():
thomas.forbriger's avatar
thomas.forbriger committed
49
50
51
    Mpara(norder,1,nparameters), Mdepth(aff::Shaper(0,1)),
    Mfollow(1,nparameters), Mnpol(1,nparameters),
    Mmeans(), Mmeans_are_valid(false)
thomas.forbriger's avatar
thomas.forbriger committed
52
53
54
55
56
    { 
      Mpara=0.;
      Mdepth=0.;
      Mfollow=false;
      Mnpol=1;
thomas.forbriger's avatar
thomas.forbriger committed
57
      clean_unused_and_check(); 
thomas.forbriger's avatar
thomas.forbriger committed
58
59
60
    }

  /*----------------------------------------------------------------------*/
61
62
63
64
65
66
67
68
69

  //! read from file
  PolynomialModelFile::PolynomialModelFile(const char* filename)
  {
    std::ifstream ifs(filename);
    read_from_stream(ifs);
  }

  /*----------------------------------------------------------------------*/
thomas.forbriger's avatar
thomas.forbriger committed
70
71
72
73
 
  //! check consistency
  void PolynomialModelFile::check_consistency() const
  {
thomas.forbriger's avatar
thomas.forbriger committed
74
75
    const char message[]="PolynomialModelFile: illegal index range";
     try {
thomas.forbriger's avatar
thomas.forbriger committed
76
      // check Mpara
thomas.forbriger's avatar
thomas.forbriger committed
77
      TFXX_assert(((Mpara.f(0)==1)&&(Mpara.l(0)==norder)), message);
78
      TFXX_assert(((Mpara.f(1)==1)&&(Mpara.l(1)>=Mpara.f(1))), message);
thomas.forbriger's avatar
thomas.forbriger committed
79
      TFXX_assert(((Mpara.f(2)==1)&&(Mpara.l(2)==nparameters)), message);
thomas.forbriger's avatar
thomas.forbriger committed
80
81
82
      // check Mdepth
      TFXX_assert(((Mdepth.f(0)==0)&&(Mdepth.l(0)>Mdepth.f(0))), message);
      // check Mfollow
83
      TFXX_assert(((Mfollow.f(0)==1)&&(Mfollow.l(0)>=Mfollow.f(0))), message);
thomas.forbriger's avatar
thomas.forbriger committed
84
      TFXX_assert(((Mfollow.f(1)==1)&&(Mfollow.l(1)==nparameters)), message);
thomas.forbriger's avatar
thomas.forbriger committed
85
      // check Mnpol
86
      TFXX_assert(((Mnpol.f(0)==1)&&(Mnpol.l(0)>=Mnpol.f(0))), message);
thomas.forbriger's avatar
thomas.forbriger committed
87
      TFXX_assert(((Mnpol.f(1)==1)&&(Mnpol.l(1)==nparameters)), message);
thomas.forbriger's avatar
thomas.forbriger committed
88
89
90
91
      // check for consistent number of sections
      TFXX_assert(((Mpara.l(1)==Mdepth.l(0))&&(Mfollow.l(0)==Mdepth.l(0))&&
                   (Mnpol.l(0)==Mdepth.l(0))), message);
    }
thomas.forbriger's avatar
thomas.forbriger committed
92
    catch ( ... ) {
thomas.forbriger's avatar
thomas.forbriger committed
93
94
95
96
97
98
99
100
101
      TFXX_abort("PolynomialModelFile: internal error!");
    }
  }

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

  //! return bottom depth of section
  Tvalue PolynomialModelFile::bottom(const int& i) const
  {
thomas.forbriger's avatar
thomas.forbriger committed
102
    TFXX_assert((!(ishalfspace(i))), "illegal section index!");
thomas.forbriger's avatar
thomas.forbriger committed
103
104
    return(Mdepth(valid_section_index(i))); 
  }
thomas.forbriger's avatar
thomas.forbriger committed
105
106
107
108
109
110
111
112
113
114

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

  //! return bottom depth of section
  Tvalue PolynomialModelFile::top(const int& i) const
  {
    TFXX_assert((!(ishalfspace(i))), "illegal section index!");
    return(Mdepth(valid_section_index(i)-1)); 
  }

thomas.forbriger's avatar
thomas.forbriger committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
  /*----------------------------------------------------------------------*/

  //! return checked section index
  int PolynomialModelFile::valid_section_index(const int& i) const
  {
    TFXX_assert(((i>0)&&(i<=Mdepth.l(0))), "illegal section index!");
    return(i);
  }

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

  //! return checked parameter index
  int PolynomialModelFile::valid_parameter_index(const int& i) const
  {
thomas.forbriger's avatar
thomas.forbriger committed
129
    TFXX_assert(((i>0)&&(i<=nparameters)), "illegal parameter index!");
thomas.forbriger's avatar
thomas.forbriger committed
130
131
132
133
134
135
136
137
    return(i);
  }

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

  //! return checked polynomial index
  int PolynomialModelFile::valid_polynomial_index(const int& i) const
  {
thomas.forbriger's avatar
thomas.forbriger committed
138
    TFXX_assert(((i>0)&&(i<=norder)), "illegal polynomial index!");
thomas.forbriger's avatar
thomas.forbriger committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    return(i);
  }

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

  //! find section index
  int PolynomialModelFile::secindex(const Tvalue& depth) const
  {
    static int retval=1;
    if (!((Mdepth(retval-1)<=depth)&&(depth<Mdepth(retval))))
    {
      if (depth>Mdepth(retval))
      {
        while ((depth>Mdepth(retval))&&(retval<Mdepth.l(0))) { retval++; }
      } else {
        while ((depth<Mdepth(retval-1))&&(retval>1)) { retval--; }
      }
    }
    return(retval);
  }

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

  //! calculate parameter value at given depth
  Tvalue PolynomialModelFile::polyarg(const Tvalue& depth,
                                      const int& isec) const
  {
    int i=valid_section_index(isec);
thomas.forbriger's avatar
thomas.forbriger committed
167
    TFXX_assert((!(ishalfspace(i))), "isec is in halfpsace"); 
168
    return(depth-0.5*(Mdepth(i)+Mdepth(i-1)));
thomas.forbriger's avatar
thomas.forbriger committed
169
170
171
172
  }

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

thomas.forbriger's avatar
thomas.forbriger committed
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
  //! set unused polynomial coefficients to zero and check values
  void PolynomialModelFile::clean_unused_and_check()
  {
    // check consistency
    check_consistency();
    // set depth of upper halfspace
    Mdepth(Mdepth.f(0))=0.;
    // cycle all sections
    for (int isec=Mpara.f(1); isec<=Mpara.l(1); ++isec)
    {
      // cycle all parameters
      for (int ipar=Mpara.f(2); ipar<=Mpara.l(2); ++ipar)
      {
        if (ishalfspace(isec))
        {
          Mnpol(isec,ipar)=1;
          Mfollow(isec,ipar)=true;
        }
        if (isec==Mpara.f(1))
        {
          Mfollow(isec,ipar)=false;
        }
        switch (Mnpol(isec,ipar)) {
          // intentionally fall through
          case 1: Mpara(2, isec,ipar)=0.;
          case 2: Mpara(3, isec,ipar)=0.;
          case 3: break;
          default: TFXX_abort("illegal polynomial order!");
        }
      }
      // check depth
      TFXX_assert(((Mdepth(isec)>Mdepth(isec-1)) || ishalfspace(isec)),
                  "depth must increase!");
    }
    make_follow();
  }

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

thomas.forbriger's avatar
thomas.forbriger committed
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
  //! calculate parameter value at given depth
  Tvalue PolynomialModelFile::value(const int& ipar,
                                    const Tvalue& depth,
                                    const int& isec) const
  {
    int is=valid_section_index(isec);
    int ip=valid_parameter_index(ipar);
    Tvalue retval=Mpara(1,is,ip);
    if (!ishalfspace(is))
    {
      Tvalue x=polyarg(depth,is);
      retval=Mpara(1,is,ip)+x*(Mpara(2,is,ip)+x*Mpara(3,is,ip));
    }
    return(retval);
  }

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

thomas.forbriger's avatar
thomas.forbriger committed
230
231
232
  //! concatenate follow-sections
  void PolynomialModelFile::make_follow()
  {
thomas.forbriger's avatar
thomas.forbriger committed
233
234
235
236
237
238
239
240
241
242
243
    // cycle all sections
    for (int isec=Mpara.f(1); isec<=Mpara.l(1); ++isec)
    {
      // cycle all parameters
      for (int ipar=Mpara.f(2); ipar<=Mpara.l(2); ++ipar)
      {
        if (Mfollow(isec,ipar))
        {
          TFXX_assert((isec!=Mpara.f(1)),
                      "first section has nothing to follow!");
          // calculate value at bottom of upper section
244
          Tvalue vup=value(ipar, Mdepth(isec-1), (isec-1));
thomas.forbriger's avatar
thomas.forbriger committed
245
246
247
248
249
250
251
          // calculate value at top of this section
          Tvalue vlo=value(ipar, Mdepth(isec-1), isec);
          // correct mean of this section
          Mpara(1,isec,ipar)+=(vup-vlo);
        }
      }
    }
thomas.forbriger's avatar
thomas.forbriger committed
252
    Mmeans_are_valid=false;
thomas.forbriger's avatar
thomas.forbriger committed
253
254
255
256
257
258
259
260
261
262
263
  }

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

  //! output
  void PolynomialModelFile::write_to_stream(std::ostream& os) const
  {
    check_consistency();
    os << TF_POLYMODEL_CC_CVSID << std::endl << std::endl;
    os << modversion2 << std::endl;
    os.width(5);
264
265
266
    os << nsections();
    os.width(15);
    os << " " << "<-- number of sections" << std::endl << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
267
268
269
270
271
272
273
274
275
    os << "    Vp             Vs             density        Qp             Qs"
      << std::endl;
    const int bufsize=15;
    char buffer[bufsize];
    for (int isec=Mpara.f(1); isec<Mpara.l(1); ++isec)
    {
      os << std::endl;
      std::snprintf(buffer, bufsize, "%12.4f", Mdepth(isec));
      os.width(12);
276
277
278
279
280
      os << buffer;
      os.width(8);
      os << " " << "<-- bottom of section";
      os.width(5);
      os << isec << " / polynomial expansion:" << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
281
      for (int ipar=1; ipar<=nparameters; ++ipar)
thomas.forbriger's avatar
thomas.forbriger committed
282
283
284
285
286
287
288
      {
        os.width(12);
        int npol=Mnpol(isec,ipar);
        if (Mfollow(isec,ipar)) npol *= -1;
        os << npol << "   ";
      }
      os << std::endl << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
289
      for (int ipol=1; ipol<=norder; ++ipol)
thomas.forbriger's avatar
thomas.forbriger committed
290
      {
thomas.forbriger's avatar
thomas.forbriger committed
291
        for (int ipar=1; ipar<=nparameters; ++ipar)
thomas.forbriger's avatar
thomas.forbriger committed
292
293
294
        {
          if (ipar<4)
          {
295
            std::snprintf(buffer, bufsize, "%12.6f", Mpara(ipol,isec,ipar));
thomas.forbriger's avatar
thomas.forbriger committed
296
297
298
299
300
301
          } else {
            std::snprintf(buffer, bufsize, "%12.4f", Mpara(ipol,isec,ipar));
          }
          os.width(12);
          os << buffer << "   ";
        }
302
        os << std::endl;
thomas.forbriger's avatar
thomas.forbriger committed
303
304
305
306
307
308
309
310
311
312
      }
    }
  }

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

  //! input 
  void PolynomialModelFile::read_from_stream(std::istream& is)
  {
    // skip comment
313
314
315
    const int maxskip=500;
    is.ignore(maxskip,'\n');
    is.ignore(maxskip,'\n');
thomas.forbriger's avatar
thomas.forbriger committed
316
    std::string identifier;
317
318
    getline(is, identifier);
    TFXX_assert((identifier.find(modversion2)==0),
thomas.forbriger's avatar
thomas.forbriger committed
319
320
321
322
                "missing or wrong identifier!");
    // number of sections
    int nsec;
    is >> nsec;
323
    is.ignore(maxskip,'\n');
thomas.forbriger's avatar
thomas.forbriger committed
324
    // prepare arrays
325
    try {
thomas.forbriger's avatar
thomas.forbriger committed
326
      Mpara=Tarray(norder,nsec+1,nparameters);
327
      Mdepth=Tarray(aff::Shaper(0,nsec+1));
thomas.forbriger's avatar
thomas.forbriger committed
328
329
      Mfollow=Tbarray(nsec+1,nparameters);
      Mnpol=Tiarray(nsec+1,nparameters);
330
331
332
333
334
335
    }
    catch ( ... ) {
      std::cerr << "PolynomialModelFile: could not create arrays for "
        << nsec << " sections!" << std::endl;
      TFXX_abort("aborting...");
    }
thomas.forbriger's avatar
thomas.forbriger committed
336
337
338
339
340
    Mdepth=0.;
    Mpara=0.;
    Mfollow=false;
    Mnpol=1;
    // skip headline
341
342
    is.ignore(maxskip,'\n');
    is.ignore(maxskip,'\n');
thomas.forbriger's avatar
thomas.forbriger committed
343
344
345
    // read
    for (int isec=1; isec<=nsec; ++isec)
    {
346
      is.ignore(maxskip,'\n');
thomas.forbriger's avatar
thomas.forbriger committed
347
      is >> Mdepth(isec);
348
      is.ignore(maxskip,'\n');
thomas.forbriger's avatar
thomas.forbriger committed
349
      for (int ipar=1; ipar<=nparameters; ++ipar)
thomas.forbriger's avatar
thomas.forbriger committed
350
351
352
353
354
355
356
357
      {
        is >> Mnpol(isec, ipar);
        if (Mnpol(isec,ipar)<0)
        {
          Mnpol(isec,ipar) *= -1;
          Mfollow(isec,ipar)=true;
        }
      }
thomas.forbriger's avatar
thomas.forbriger committed
358
      for (int ipol=1; ipol<=norder; ++ipol)
thomas.forbriger's avatar
thomas.forbriger committed
359
      {
thomas.forbriger's avatar
thomas.forbriger committed
360
        for (int ipar=1; ipar<=nparameters; ++ipar)
thomas.forbriger's avatar
thomas.forbriger committed
361
362
363
364
365
366
367
        {
          is >> Mpara(ipol, isec, ipar);
        }
      }
    }
    // finalize
    clean_unused_and_check();
thomas.forbriger's avatar
thomas.forbriger committed
368
    Mmeans_are_valid=false;
thomas.forbriger's avatar
thomas.forbriger committed
369
  }
thomas.forbriger's avatar
thomas.forbriger committed
370
371
372
373
374
375
376
377
378

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

  //! calculate mean values for all sections
  void PolynomialModelFile::calculate_means() const
  {
    if (!Mmeans_are_valid)
    {
      Mmeans=Tarray(norder,nsections(),nparameters);
thomas.forbriger's avatar
thomas.forbriger committed
379
      for (int isec=Mmeans.f(1); isec<=Mmeans.l(1); ++isec)
thomas.forbriger's avatar
thomas.forbriger committed
380
      {
381
        for (int ipar=Mmeans.f(2); ipar<=Mmeans.l(2); ++ipar)
thomas.forbriger's avatar
thomas.forbriger committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
        {
          Mmeans(3,isec,ipar)=2.*Mpara(3,isec,ipar);
          Mmeans(2,isec,ipar)=Mpara(2,isec,ipar);
          Mmeans(1,isec,ipar)=Mpara(1,isec,ipar)+
            Mpara(3,isec,ipar)*thickness(isec)*thickness(isec)/3.;
        }
      }
      Mmeans_are_valid=true;
    }
  }

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

  //! calculate mean values for depth range
  /*
  TCarray PolynomialModelFile::means(const Tvalue& top
                                     const Tvalue& bottom) const
  {
  }
  */

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

  //! calculate derivative values at given depth
  Tarray PolynomialModelFile::values(const Tvalue& depth,
                                     const int& isec) const
  {
    int is=valid_section_index(isec);
    Tarray result(norder,nparameters);
    Tarray secpar(aff::slice(Mpara)()(is));
    result.copyin(secpar);
    if (!ishalfspace(is))
    {
      Tvalue x=polyarg(depth,is);
      for (int ipar=result.f(1); ipar<=result.l(1); ++ipar)
      {
        result(3,ipar)=2.*secpar(3,ipar);
        result(2,ipar)=secpar(2,ipar)+2.*secpar(3,ipar)*x;
        result(1,ipar)=secpar(1,ipar)+x*(secpar(2,ipar)+x*secpar(3,ipar));
      }
    }
    // NOTICE: returns reference (but result should have no reference to data
    // members)
    return(result);
  }

428
} // namespace gremlin1
thomas.forbriger's avatar
thomas.forbriger committed
429
430

/* ----- END OF polymodel.cc ----- */