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

nice and working

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: 1271
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 17be7596
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: mocox.cc,v 1.3 2002-12-31 17:20:52 forbrig Exp $
* $Id: mocox.cc,v 1.4 2003-01-03 22:21:50 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -17,12 +17,15 @@
*
* ============================================================================
*/
#define MOCOX_SHORT_VERSION \
"MOCOX V1.1"
#define MOCOX_VERSION \
"MOCOX V1.1 convert gremlin1 model (modversion2)"
MOCOX_SHORT_VERSION " convert gremlin1 model (modversion2)"
#define MOCOX_CVSID \
"$Id: mocox.cc,v 1.3 2002-12-31 17:20:52 forbrig Exp $"
"$Id: mocox.cc,v 1.4 2003-01-03 22:21:50 forbrig Exp $"
#include <iostream>
#include <stdio.h>
#include <cmath>
#include <fstream>
#include <gremlin1/polymodel.h>
......@@ -32,11 +35,14 @@
#include <aff/slice.h>
#include <aff/subarray.h>
#include <aff/iterator.h>
#include <aff/dump.h>
using std::cout;
using std::cerr;
using std::endl;
const char flnode_id[]="FLNODE";
typedef gremlin1::PolynomialModelFile::Tarray Tarray;
typedef gremlin1::PolynomialModelFile::Tiarray Tiarray;
typedef gremlin1::PolynomialModelFile::TCarray TCarray;
......@@ -45,16 +51,129 @@ typedef gremlin1::PolynomialModelFile::TCarray TCarray;
// functions
// =========
// find nodes for flnode file
Tarray defnodes(const PolynomialModelFile& modelfile,
const bool& verbose)
// create a complete set of nodes
void dump_nodes(const TCarray& nodes)
{
cout << endl
<< "Converted model:" << endl;
const int fieldwidth=12;
char buffer[fieldwidth];
// captions
cout << " ";
cout.width(fieldwidth); cout << "depth (m)";
cout.width(fieldwidth); cout << "v_P (km/s)";
cout.width(fieldwidth); cout << "v_S (km/s)";
cout.width(fieldwidth); cout << "rho (g/cm^3)";
cout.width(fieldwidth); cout << "Q_P";
cout.width(fieldwidth); cout << "Q_S";
cout << endl;
// values
for (int i=nodes.f(0); i<=nodes.l(0); ++i)
{
cout << " ";
for (int j=nodes.f(1); j<=nodes.l(1); ++j)
{
std::snprintf(buffer, fieldwidth, "%10.4f", nodes(i,j));
cout.width(fieldwidth);
cout << buffer;
}
cout << endl;
}
} // dump_nodes
/*----------------------------------------------------------------------*/
// create a complete set of nodes
Tarray convert_flnode(const TCarray& nodes, const bool& verbose)
{
Tarray flnodes(nodes.size(0)+1,nodes.size(1));
int nf=nodes.f(0);
int nl=nodes.l(0);
flnodes=0.;
Tarray setflnodes=subarray(flnodes)(nf,nl);
TCarray indepth =slice(nodes)()(1);
TCarray invp =slice(nodes)()(2);
TCarray invs =slice(nodes)()(3);
TCarray inrho =slice(nodes)()(4);
TCarray inqp =slice(nodes)()(5);
TCarray inqs =slice(nodes)()(6);
Tarray fldepth =slice(setflnodes)()(1);
Tarray flrho =slice(setflnodes)()(2);
Tarray flvp =slice(setflnodes)()(3);
Tarray flvs =slice(setflnodes)()(4);
Tarray flqk =slice(setflnodes)()(5);
Tarray flqm =slice(setflnodes)()(6);
aff::deepcopy(indepth,fldepth);
aff::deepcopy(invp ,flvp );
aff::deepcopy(invs ,flvs );
aff::deepcopy(inrho ,flrho );
aff::deepcopy(inqp ,flqk );
aff::deepcopy(inqs ,flqm );
Tarray fllast=slice(flnodes)(flnodes.l(0));
Tarray setfllast=slice(setflnodes)(setflnodes.l(0));
aff::deepcopy(setfllast,fllast);
return(flnodes);
} // convert_flnode
/*----------------------------------------------------------------------*/
// create a complete set of nodes
void write_flnode(const std::string& filename,
const TCarray& nodes,
const bool& verbose,
const std::string& sourcefile="unknown")
{
Tarray flnodes=convert_flnode(nodes, verbose);
std::ofstream ofs(filename.c_str());
ofs << flnode_id;
ofs << "converted by " << MOCOX_SHORT_VERSION
<< " from " << sourcefile << endl;
const double reference_frequency=0.;
const double earth_radius=6371.;
const int damping_mode=1; // 1 = no dispersion
const int anisotropy_mode=0; // 0 = isotropic model
const int fieldwidth=12;
char buffer[fieldwidth];
std::snprintf(buffer, fieldwidth, "%10.4f", reference_frequency);
ofs << buffer << " ";
ofs << damping_mode << " ";
ofs << anisotropy_mode << " ";
std::snprintf(buffer, fieldwidth, "%10.4f", earth_radius);
ofs << buffer << endl;
ofs << flnodes.size(0) << endl;
for (int i=flnodes.f(0); i<=flnodes.l(0); ++i)
{
ofs << " ";
for (int j=flnodes.f(1); j<=flnodes.l(1); ++j)
{
std::snprintf(buffer, fieldwidth, "%10.4f", flnodes(i,j));
ofs.width(fieldwidth);
ofs << buffer;
}
ofs << endl;
}
} // write_flnode
/*----------------------------------------------------------------------*/
// create a complete set of nodes
Tarray create_nodes(const gremlin1::PolynomialModelFile& modelfile,
const bool& verbose)
{
if (verbose) cout << " define nodes:" << endl;
// create local alias
int nsections=modelfile.nsections();
TCarray para=modelfile.parameters();
TCarray depth=modelfile.depth();
TCarray means=modelfile.sectionmeans();
/* scale factor
* May be passed as function argument in future versions.
......@@ -80,31 +199,92 @@ Tarray defnodes(const PolynomialModelFile& modelfile,
double scalefactor=0.01;
if (verbose) cout << " find appropriate number of nodes per section"
<< endl;
<< endl << " scale-factor: " << scalefactor << endl;
// array to store first and last node index
Tiarray ifirst(nsections);
Tiarray ilast(nsections);
// array to store number of node layers per section
Tiarray nlay(nsections);
nlay=1;
int totalnlay=0;
int totalnodes=0;
// find number of layers per section
double bottom_of_upper=0.;
for (int i=1,i<=nsections; ++i)
for (int i=1; i<=nsections; ++i)
{
if (verbose) cout << " sec #" << i <<":";
double d=depth(i)-bottom_of_upper;
if (verbose) cout << " d=" << d << "m";
ifirst(i)=totalnodes+1;
double d=modelfile.thickness(i);
if (verbose) {
cout << " d=";
cout.width(8); cout.precision(4);
cout << d << "m";
cout << " z=";
cout.width(8); cout.precision(4);
cout << modelfile.bottom(i) << "m";
}
// all second order coefficients in this section
TCarray c(aff::slice(para)(3)(i));
for (int j=c.f(0); j<=c.l(0); ++)
TCarray c(aff::slice(means)(3)(i));
for (int j=c.f(0); j<=c.l(0); ++j)
{
int n=1+int(d*sqrt(2.*abs(c(j))/(3.*scalefactor)));
int n=1+static_cast<int>(d*sqrt(std::abs(c(j))/(3.*scalefactor)));
nlay(i)= nlay(i)>n ? nlay(i) : n;
}
if (verbose) cout << " nlay=" << nlay(i);
totalnlay += nlay(i);
bottom_of_upper=depth(i);
if (verbose) {
cout << " nlay=";
cout.width(3);
cout << nlay(i);
}
// we need a node per layer and an extra node at each interface
totalnodes += nlay(i)+1;
ilast(i)=totalnodes;
if (verbose) {
cout << " first=";
cout.width(3);
cout << ifirst(i);
cout << " last=";
cout.width(3);
cout << ilast(i);
}
if (verbose) cout << endl;
}
} // defnodes
// fill result with depth and values
// result array:
Tarray result(totalnodes, gremlin1::PolynomialModelFile::nparameters+1);
if (verbose) cout << " define node depth and fill in values" << endl;
for (int i=1; i<=nsections; ++i)
{
double top=modelfile.top(i);
double dz=modelfile.thickness(i)/nlay(i);
for (int j=0; j<=nlay(i); ++j)
{
if (verbose) cout << " sec #" << i <<":";
// index of node
int inode=ifirst(i)+j;
// depth of node
double depth=j*dz+top;
result(inode,1)=depth;
if (verbose) {
cout << " node=";
cout.width(3);
cout << inode;
cout << " depth=";
cout.width(8);
cout.precision(4);
cout << depth << "m";
}
for (int k=result.f(1); k<result.l(1); ++k)
{
result(inode, 1+k)=modelfile.value(k,depth,i);
}
if (verbose) cout << endl;
}
}
return(result);
} // create_nodes
/*======================================================================*/
// main program
......@@ -118,6 +298,7 @@ int main(int iargc, char* argv[])
{
MOCOX_VERSION "\n"
"usage: mocox file [-w file] [-v] [-dump] [-flnode file]" "\n"
" [-nodes]" "\n"
" or: mocox --help|-h" "\n"
};
......@@ -132,6 +313,7 @@ int main(int iargc, char* argv[])
"-dump dump model to terminal" "\n"
"-w file write model to gremlin1 format \"file\"" "\n"
"-flnode file convert to FLNODE model \"file\"" "\n"
"-nodes convert to nodes file and dump" "\n"
"\n"
MOCOX_CVSID
};
......@@ -150,6 +332,8 @@ int main(int iargc, char* argv[])
{"dump",arg_no,"-"},
// 4: convert to FLNODE model
{"flnode",arg_yes,"-"},
// 5: dump input model to terminal
{"nodes",arg_no,"-"},
{NULL}
};
......@@ -180,6 +364,7 @@ int main(int iargc, char* argv[])
bool opt_dump=cmdline.optset(3);
bool opt_flnode=cmdline.optset(4);
std::string arg_flnode=cmdline.string_arg(4);
bool opt_nodes=cmdline.optset(5);
TFXX_assert(cmdline.extra(), "ERROR (mocox): missing file name!");
std::string filename=cmdline.next();
......@@ -212,11 +397,23 @@ int main(int iargc, char* argv[])
/*======================================================================*/
// convert to flnode is requested
if (opt_flnode)
if (opt_flnode||opt_nodes)
{
if (opt_verbose) cout << endl << "convert to FLNODE model:" << endl;
if (opt_verbose) cout << endl << "convert to nodes model:" << endl;
// first find out the number of nodes we beed
Tarray nodes=defnodes(modelfile, opt_verbose);
Tarray nodes=create_nodes(modelfile, opt_verbose);
// dump nodes if requested
if (opt_nodes)
{
dump_nodes(nodes);
}
// write flnode file if requested
if (opt_flnode)
{
write_flnode(arg_flnode,nodes,opt_verbose,filename);
}
} // if (opt_flnode)
......
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