Commit c5b2f2de authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

compiles...

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1247
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 012721c6
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id: Makefile,v 1.1 2002-12-30 13:03:40 forbrig Exp $
# $Id: Makefile,v 1.2 2002-12-30 16:04:32 forbrig Exp $
#
# Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
#
......@@ -209,9 +209,7 @@ DOXYWWWPATH=$(WWWBASEDIR)/libgremlin1xx
doxyclean: ;/bin/rm -rfv docfull/* docbrief docfull
DOXYSRC=$(README) $(HEADERS) $(SRC) \
tests/f77procs.P tests/f77procs.f \
tests/f77common.inc tests/f77common_com.P
DOXYSRC=$(README) $(HEADERS) $(SRC)
doc%/html/index.html: doxy%.cfg $(DOXYSRC)
mkdir -vp $(DOXYWWWPATH)
......
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: polymodel.cc,v 1.1 2002-12-30 13:03:41 forbrig Exp $
* $Id: polymodel.cc,v 1.2 2002-12-30 16:04:32 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -17,15 +17,23 @@
* ============================================================================
*/
#define TF_POLYMODEL_CC_VERSION \
"TF_POLYMODEL_CC V1.0 "
"TF_POLYMODEL_CC V1.0"
#define TF_POLYMODEL_CC_CVSID \
"$Id: polymodel.cc,v 1.1 2002-12-30 13:03:41 forbrig Exp $"
"$Id: polymodel.cc,v 1.2 2002-12-30 16:04:32 forbrig Exp $"
#include <gremlin1/polymodel.h>
#include <tfxx/error.h>
#include <aff/shaper.h>
#include <stdio.h>
#include <string>
namespace gremlin1 {
typedef PolynomialModelFile::Tvalue Tvalue;
//! define identifier
const char PolynomialModelFile::modversion2[];
//! default constructor
PolynomialModelFile::PolynomialModelFile():
Mpara(3,1,5), Mdepth(aff::Shaper(0,1)), Mfollow(1,5), Mnpol(1,5)
......@@ -34,7 +42,7 @@ namespace gremlin1 {
Mdepth=0.;
Mfollow=false;
Mnpol=1;
check_consistency();
clean_unused_and_check();
}
/*----------------------------------------------------------------------*/
......@@ -42,8 +50,8 @@ namespace gremlin1 {
//! check consistency
void PolynomialModelFile::check_consistency() const
{
const char[] message="PolynomialModelFile: illegal index range";
try {
const char message[]="PolynomialModelFile: illegal index range";
try {
// check Mpara
TFXX_assert(((Mpara.f(0)==1)&&(Mpara.l(0)==3)), message);
TFXX_assert(((Mpara.f(1)==0)&&(Mpara.l(1)>Mpara.f(1))), message);
......@@ -60,8 +68,7 @@ namespace gremlin1 {
TFXX_assert(((Mpara.l(1)==Mdepth.l(0))&&(Mfollow.l(0)==Mdepth.l(0))&&
(Mnpol.l(0)==Mdepth.l(0))), message);
}
catch (tfxx::Exception)
{
catch ( ... ) {
TFXX_abort("PolynomialModelFile: internal error!");
}
}
......@@ -71,7 +78,7 @@ namespace gremlin1 {
//! return bottom depth of section
Tvalue PolynomialModelFile::bottom(const int& i) const
{
TFXX_assert((!(ishalfspace())), "illegal section index!");
TFXX_assert((!(ishalfspace(i))), "illegal section index!");
return(Mdepth(valid_section_index(i)));
}
/*----------------------------------------------------------------------*/
......@@ -126,32 +133,195 @@ namespace gremlin1 {
const int& isec) const
{
int i=valid_section_index(isec);
TFXX_assert((!(ishalfspace(isec))), "isec is in halfpsace");
return(depth-0.5(Mdepth(isec)-Msec(isec-1)));
TFXX_assert((!(ishalfspace(i))), "isec is in halfpsace");
return(depth-0.5*(Mdepth(i)-Mdepth(i-1)));
}
/*----------------------------------------------------------------------*/
//! calculate parameter value at given depth
Tvalue PolynomialModelFile::value(const Epara& ipar,
const Tvalue& depth) const
Tvalue PolynomialModelFile::value(const int& ipar,
const Tvalue& depth,
const int& isec) const
{
int isec=secindex(depth);
int is=valid_section_index(isec);
int ip=valid_parameter_index(ipar);
Tvalue retval=Mpara(1,isec,ip);
if (!ishalfspace(isec))
Tvalue retval=Mpara(1,is,ip);
if (!ishalfspace(is))
{
Tvalue x=polyarg(depth,isec);
retval=Mpara(1,isec,ip)+x*(Mpara(2,isec,ip)+x*Mpara(3,isec,ip));
Tvalue x=polyarg(depth,is);
retval=Mpara(1,is,ip)+x*(Mpara(2,is,ip)+x*Mpara(3,is,ip));
}
return(retval);
}
/*----------------------------------------------------------------------*/
//! 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();
}
/*----------------------------------------------------------------------*/
//! concatenate follow-sections
void PolynomialModelFile::make_follow()
{
// 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
Tvalue vup=value(ipar, Mdepth(isec-1), isec-1);
// 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);
}
}
}
}
/*----------------------------------------------------------------------*/
//! 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);
os << nsections()
<< " <-- number of sections" << std::endl << std::endl;
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);
os << buffer << " <-- bottom of section " << isec
<< " polynomial expansion:" << std::endl;
for (int ipar=1; ipar<6; ++ipar)
{
os.width(12);
int npol=Mnpol(isec,ipar);
if (Mfollow(isec,ipar)) npol *= -1;
os << npol << " ";
}
os << std::endl << std::endl;
for (int ipol=1; ipol<4; ++ipol)
{
for (int ipar=1; ipar<6; ++ipar)
{
if (ipar<4)
{
std::snprintf(buffer, bufsize, "%12.4g", Mpara(ipol,isec,ipar));
} else {
std::snprintf(buffer, bufsize, "%12.4f", Mpara(ipol,isec,ipar));
}
os.width(12);
os << buffer << " ";
}
os << std::endl << std::endl;
}
}
}
/*----------------------------------------------------------------------*/
//! input
void PolynomialModelFile::read_from_stream(std::istream& is)
{
// skip comment
is.ignore(500,'\n');
is.ignore(500,'\n');
std::string identifier;
is >> identifier;
TFXX_assert((identifier == std::string(modversion2)),
"missing or wrong identifier!");
is.ignore(500,'\n');
// number of sections
int nsec;
is >> nsec;
is.ignore(500,'\n');
// prepare arrays
Mpara=Tarray(3,nsec+1,5);
Mdepth=Tarray(aff::Shaper(0,nsec+1));
Mfollow=Tbarray(nsec+1,5);
Mnpol=Tiarray(nsec+1,5);
Mdepth=0.;
Mpara=0.;
Mfollow=false;
Mnpol=1;
// skip headline
is.ignore(500,'\n');
is.ignore(500,'\n');
// read
for (int isec=1; isec<=nsec; ++isec)
{
is.ignore(500,'\n');
is >> Mdepth(isec);
is.ignore(500,'\n');
for (int ipar=1; ipar<6; ++ipar)
{
is >> Mnpol(isec, ipar);
if (Mnpol(isec,ipar)<0)
{
Mnpol(isec,ipar) *= -1;
Mfollow(isec,ipar)=true;
}
}
for (int ipol=1; ipol<4; ++ipol)
{
for (int ipar=1; ipar<6; ++ipar)
{
is >> Mpara(ipol, isec, ipar);
}
}
}
// finalize
clean_unused_and_check();
}
}
......
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: polymodel.h,v 1.1 2002-12-30 13:03:41 forbrig Exp $
* $Id: polymodel.h,v 1.2 2002-12-30 16:04:32 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -23,7 +23,7 @@
#define TF_POLYMODEL_H_VERSION \
"TF_POLYMODEL_H V1.0 "
#define TF_POLYMODEL_H_CVSID \
"$Id: polymodel.h,v 1.1 2002-12-30 13:03:41 forbrig Exp $"
"$Id: polymodel.h,v 1.2 2002-12-30 16:04:32 forbrig Exp $"
#include<iostream>
#include<aff/array.h>
......@@ -46,6 +46,10 @@ namespace gremlin1 {
typedef aff::Array<bool> Tbarray;
//! int array
typedef aff::Array<int> Tiarray;
//! model file identifier
static const char modversion2[]="ModVersion 2";
//! physical model parameters
enum Epara {
//! P-velocity
......@@ -69,7 +73,7 @@ namespace gremlin1 {
void write_to_stream(std::ostream& os) const;
//! return number of sections
int nsections() const { return(Mpara,size(1)-2); }
int nsections() const { return(Mpara.size(1)-2); }
//! lower interface of section \p i
Tvalue bottom(const int& i) const;
//! is this section index the half-space
......@@ -86,19 +90,29 @@ namespace gremlin1 {
int secindex(const Tvalue& depth) const;
//! return parameter value at given depth
Tvalue value(const Epara& ipar,
const Tvalue& depth) const;
const Tvalue& depth) const
{ return(value(ipar, depth, secindex(depth))); }
//! return parameter value at given depth calculated from
//! polynomial for given section \p isec
Tvalue value(const int& ipar,
const Tvalue& depth,
const int& isec) const;
//! return agrument of polynomial for given depth and section
Tvalue polyarg(const Tvalue& depth, const int& isec) const;
private:
//! check consistency of stored data
void check_consistency() const;
//! check for valid section index
int valid_section_index(const int& i) const;
//! check for valid polynomial index
int valid_polynomial_index(const int& i) const;
//! check for valid parameter index
int valid_parameter_index(const int& i) const;
//! clean all data (should be called after modifications)
//! calls check_consistency and make_follow
void clean_unused_and_check();
//! check consistency of stored data
void check_consistency() const;
//! make parameters continuous, where follow flag is set
void make_follow();
......@@ -136,12 +150,12 @@ namespace gremlin1 {
//! output operator
std::ostream& operator<<(std::ostream& os,
const PolynomialModelFile& pmf)
{ pmf.write_to_stream(os); }
{ pmf.write_to_stream(os); return(os); }
//! output operator
std::istream& operator>>(std::istream& is,
PolynomialModelFile& pmf)
{ pmf.read_from_stream(is); }
std::istream& operator>>(std::istream& is,
PolynomialModelFile& pmf)
{ pmf.read_from_stream(is); return(is); }
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment