Code:
// -- Program --------------------------------------------------------------------------------------
//
// getvel
// Creates an "x z v" file for plotting the sound speed values using GMT.
//
// -- Usage ----------------------------------------------------------------------------------------
//
// ./getvel Args
//
// Mandatory Arguments:
//
// Optional Arguments:
//
// -- History --------------------------------------------------------------------------------------
//
// V01 - DEC 1999 - Luca D'Auria
// Initial release.
// V02 - SEP 2008 - Luca Elia
// Version fixed for g++ 4.3.0.
// V03 - MAY 2009 - Christopher Dimech
// Adopted C++ coding standards based on the IBM Standard C++
// Library Reference Guide and the Ellementel Telecommunication
// Systems Laboratory C++ Style Guide.
// V04 - NOV 2009 - Luca D'Auria
// Code changes after visit to RISSC-Lab.
// V05 - JUL 2011 - Christopher Dimech
// Adopted UML 2.3 for documenting class structure.
//
// -------------------------------------------------------------------------------------------------
#include "Parsing.hh"
#include "GenFunc.hh"
#include "VelMod.hh"
#include "common.hh"
////////////////////////////////////////////////////////////////////////////////////////////////////
const int DefKx = 20;
const int DefKz = 20;
const int DefKI = 20;
const char *DOC[] = {
"",
"-- Program --------------------------------------------------------------------------------------",
"",
" GetVel",
" Creation of a \"x z v\" file",
"",
"-- Usage ----------------------------------------------------------------------------------------",
"",
" getvel model [arguments]",
"",
" model",
" Velocity model file (required)",
"",
" vdx",
" Interval in x direction in velocity file",
"",
" kx",
" Interval calculated with vdx = Lx / kx (Default is 20).",
"",
" vdz",
" Interval in z direction in velocity file or",
"",
" kz",
" Interval calculated with vdz = Lz / kz (Default is 20).",
"",
" vdi",
" Interval along interface",
"",
" ki",
" Interval calculated with vdi = Li / ki (Default is 20).",
"",
" vel or velp",
" File with P wave velocity in 'x z v' format",
"",
" vels",
" File with S wave velocity in 'x z v' format",
"",
" lay",
" (on/off). If on output files are in \"x z v l\" format, where l",
" is the active layer.",
"",
" intf",
" Output file for interfaces in 'x z' format, separated by '>'",
"",
" see HTML help for more details",
"",
"-------------------------------------------------------------------------------------------------",
"",
NULL
};
////////////////////////////////////////////////////////////////////////////////////////////////////
int main(
int argc,
char* argv[]) {
if (argc < 2) { // Print documentation
FILE* pipe;
const char** ptr = DOC;
pipe = popen("more", "w");
while (*ptr) { (void)fprintf(pipe, "%s\n", *ptr++); }
pclose(pipe);
exit(1);
}
// -- Reading command line arguments -------------------------------------------------------------
bool lay = false;
int Kx;
int Kz;
int KI;
Real Lx;
Real Lz;
Real Vdx;
Real Vdz;
Real Idx;
Vect2 Xi;
Vect2 Xf;
String S;
Velmod vm;
Parsing Pc(argc, argv);
// -----------------------------------------------------------------------------------------------
if (Pc.GetString("CMOD", S)) {
ifstream model(S);
Parsing Pm(model);
vm.Read(Pm);
Pm.GetVect2("XI", Xi);
Pm.GetVect2("XF", Xf);
Lx = Xf.X - Xi.X;
Lz = Xf.Z - Xi.Z;
} else {
error("Bad Vdx parameter");
}
// -----------------------------------------------------------------------------------------------
if ( !Pc.GetReal("VDX", Vdx) ) {
if ( !Pc.GetInt("KX", Kx) ) { Kx = DefKx; }
if (Kx < 0) { error("Bad Kx parameter"); }
Vdx = Lx / Kx;
} else {
Kx = (int)rint(Lx / Vdx);
Vdx = Lx / Kx;
}
if (Vdx < 0) { error("Bad Vdx parameter"); }
// -----------------------------------------------------------------------------------------------
if ( !Pc.GetReal("VDZ", Vdz) ) {
if ( !Pc.GetInt("KZ", Kz) ) { Kz = DefKz; }
if (Kz < 0) { error("Bad Kz parameter"); }
Vdz = Lz / Kz;
} else {
Kz = (int)rint(Lz / Vdz);
Vdz = Lz / Kz;
}
if (Vdz < 0) { error("Bad Vdz parameter"); }
// -----------------------------------------------------------------------------------------------
if ( !Pc.GetReal("IDX", Idx) ) {
if ( !Pc.GetInt("KI", KI) ) { KI = DefKI; }
if (KI < 0) { error("Bad KI parameter"); }
Idx = Lx / KI;
} else {
KI = (int)rint(Lx / Idx);
Idx = Lx / KI;
}
if (Idx < 0) { error("Bad Idx parameter"); }
// -----------------------------------------------------------------------------------------------
bool bvp = true;
bool bvs = true;
bool bintf = true;
ofstream velp;
ofstream vels;
ofstream intf;
if (Pc.GetString("VEL", S)) {
velp.open(S);
if (velp.bad()) { error("Bad VEL file"); }
} else if (Pc.GetString("VELP", S)) {
velp.open(S);
if (velp.bad()) { error("Bad VELP file"); }
} else {
bvp = false;
}
// -----------------------------------------------------------------------------------------------
if (Pc.GetString("VELS", S)) {
vels.open(S);
if ( vels.bad() ) { error("Bad VELS file"); }
} else {
bvs = false;
}
// -----------------------------------------------------------------------------------------------
if (Pc.GetString("INTF", S)) {
intf.open(S);
if (intf.bad()) { error("Bad INTF file"); }
} else {
bintf = false;
}
// -----------------------------------------------------------------------------------------------
if (Pc.GetString("LAY", S)) {
S.ToUpper();
if (S == String("ON")) { lay = true;
} else if (S != String("OFF")) { error("Bad LAY parameter"); }
}
// -----------------------------------------------------------------------------------------------
int i, j;
Real v, x, z;
if (bvp) {
for (i = 0; i <= Kx; i++) {
for (j = 0; j <= Kz; j++) {
x = Xi.X + (i * Vdx);
z = Xi.Z + (j * Vdz);
v = vm.VP(x, z);
velp << x << ' ' << z << ' ' << v;
if (lay) { velp << ' ' << vm.GetLayer(x, z) << endl;
} else { velp << endl; }
}
}
velp.close();
}
// -----------------------------------------------------------------------------------------------
if (bvs) {
for (i = 0; i <= Kx; i++) {
for (j = 0; j <= Kz; j++) {
x = Xi.X + (i * Vdx);
z = Xi.Z + (j * Vdz);
v = vm.VS(x, z);
vels << x << ' ' << z << ' ' << v;
if (lay) { vels << ' ' << vm.GetLayer(x, z) << endl;
} else { vels << endl; }
}
}
vels.close();
}
// -----------------------------------------------------------------------------------------------
if (bintf) {
int NLay = vm.GetNLayers(), Lay;
intf << '>' << endl;
for (Lay = 1; Lay <= NLay; Lay++) {
for (i = 0; i <= KI; i++) {
x = Xi.X + i * Idx;
z = vm.Z(x, i);
intf << x << ' ' << z << endl;
}
intf << '>' << endl;
}
intf.close();
}
return (0);
}
////////////////////////////////////////////////////////////////////////////////////////////////////