diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..d368b34 --- /dev/null +++ b/.clang-format @@ -0,0 +1,108 @@ +--- +Language: Cpp +# BasedOnStyle: LLVM +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Right +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: false +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: false +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 80 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + - Regex: '.*' + Priority: 1 +IncludeIsMainRegex: '(Test)?$' +IndentCaseLabels: false +IndentWidth: 2 +IndentWrappedFunctionNames: false +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Right +ReflowComments: true +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 8 +UseTab: Never +... + diff --git a/.gitignore b/.gitignore index 8eca4ee..d3ff22b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,17 +1,14 @@ bin/* -V2RhoT/*.o -V2RhoT/Makefile -V2RhoT/V2RhoT -V2RhoT/test.out -V2RhoT/Makefile -V2T/*.o -V2T/Makefile -V2T/V2T -V2T/docs/*synctex.gz -V2T/docs/*.log -V2T/docs/*.aux -V2T/docs/*.pdf -V2T/Makefile +/V2T +/V2RhoT +*.a +*.o +Makefile +src/*/Makefile +docs/*synctex.gz +docs/*.log +docs/*.aux +docs/*.pdf .DS_Store .qmake.stash *.app diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..6bd35d3 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,22 @@ +language: cpp +env: +- QT_BASE=4 +- QT_BASE=5 +install: +- sudo apt-get update -qq +- sudo apt-get install -qq qt${QT_BASE}-default qt${QT_BASE}-qmake +script: +- qmake +- make +after_success: +- tar -zcvf bin/bin-VeloDT-prerelease-ubuntu-qt${QT_BASE}-x86_64.tar.gz bin/V2RhoT + bin/V2T +deploy: + provider: releases + api_key: + secure: ORqoxcPCHo19J3npN9/lFvwYMOuz7qifGg8G75OKIj3p3LlFdidmAHvWD1hd0Q3+zPXzNxG/CiAuwALjWFQdLIcB0vsKeeOlUII8/eGhDkaPjyW/VOdQT5XxqSevvtRdTIqJO2WAxO80TdZ6YZoz7Q9lkxSnGkpBH569pIIC2OkgcC62HGXT3jGlJ8I82fVl32uKTh2954kRRT5eXv24rYuGhp+jarc5n2io02HLOHLdz0qg3OcemFYJfbGJ7uPjQUh3lFJd6VfLX1cdkRYBFedsWST3y7i2+SnVUpIY+mr/bEGgf2xrE8OLXw4ABA8PhFloqkiHiOrkZ6Xc4Az+cmhZjOLAT04yL1HBpXrpUb8bpJVIfH5eda8ORwey/RST8JF/eqFM7pu+AXJpG9Fs6lEm/GczuTxMOcF3cTyxM3/Q96VZq3ZEVe8REuLqqyyy5J7IMifWZPH3hkWHjTYM2NsWroG45duFW+w3ts9QrLPOOwoXYgPMQ9CJJoo5lOrlfM4rOBujbXk+tn6QoyRFhne2mxIrHE4pcyfIwDz5bsi7ZxpfFa7sWaAro/DY+OGCuneNI1DtFJHGG0QzoCNVZXBTUUgeqQQf3Sp3XOF/O++yJ804+QS+fIMaAGB/fAo6l4d13zCpWqqQBjf3a5CcBve0aovSQHntid7xfW5tBNs= + file: bin/bin-VeloDT-prerelease-ubuntu-qt${QT_BASE}-x86_64.tar.gz + on: + tags: true + repo: cmeessen/VeloDT + skip_cleanup: 'true' diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 index 0000000..e00442b --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,22 @@ +# Changelog +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). + +## [1.1.0] - 2020-06-16 + +### Added + +- CHANGELOG +- Build binaries for Ubuntu using QT4 and QT5 and release them when a tag is + added + +### Changed + +- Repository structure changed, using include and src folders + +### Fixed + +- Fixed an bug when using argument `-t` +- Fixed an typo in `readFile()` that made it impossible to use topography or + crustal thickness maps diff --git a/Example/Makefile b/Example/Makefile index cbe4a97..be61071 100644 --- a/Example/Makefile +++ b/Example/Makefile @@ -1,5 +1,5 @@ -V2RhoT = ../V2RhoT/V2RhoT -V2T = ../V2T/V2T +V2RhoT = ../bin/V2RhoT +V2T = ../bin/V2T # Plot properties FIGURE = results.ps ROI_T = 450/1450/-200/-50 @@ -9,7 +9,7 @@ MAPWIDTH = 6 # Targets -clean: plot +clean: rm Vs.dat Vp.dat V2RhoT_Vs.dat V2RhoT_Vp.dat V2T_Vs.dat gmt.* plot: convert diff --git a/README.md b/README.md index eb5e53d..43844e2 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # VeloDT -[![DOI](https://zenodo.org/badge/121383426.svg)](https://zenodo.org/badge/latestdoi/121383426) +[![DOI](https://zenodo.org/badge/121383426.svg)](https://zenodo.org/badge/latestdoi/121383426) +[![Build Status](https://travis-ci.org/cmeessen/VeloDT.svg?branch=restructure-repo)](https://travis-ci.org/cmeessen/VeloDT) *VeloDT* contains two C++ implementations of *Vp*- and/or *Vs*-conversions to temperature and/or density valid in the upper mantle. @@ -13,14 +14,17 @@ different approaches. ## Getting started Either clone the github repository with -``` + +```bash git clone https://github.com/cmeessen/VeloDT.git ``` -or download the [latest release](https://github.com/cmeessen/VeloDT/releases/latest) as source code or binary and extract it somewhere. If you downloaded the source code, you will have to compile the code. Therefore, navigate into the folders `V2RhoT` or `V2T` which both contain a README.md with instructions on compilation and usage -- [**V2RhoT**](./V2RhoT/) closely follows [Goes et al. (2000)](https://doi.org/10.1029/1999JB900300) to convert vp or +or download the [latest release](https://github.com/cmeessen/VeloDT/releases/latest) as source code or binary and extract it somewhere. If you downloaded the source code, you will have to compile the code (see below). + +- [**V2RhoT**](./V2RhoT.md) closely follows [Goes et al. (2000)](https://doi.org/10.1029/1999JB900300) to convert vp or vs to temperature and density assuming a mineral assemblage -- [**V2T**](./V2T/) implements the method by [Priestley and McKenzie + +- [**V2T**](./V2T.md) implements the method by [Priestley and McKenzie (2006)](https://doi.org/10.1016/j.epsl.2006.01.008) and converts *Vs* to temperature Please refer to the original publications for validity and pitfalls of the methods. @@ -41,13 +45,15 @@ Requirements for plotting the example: ### Compiling -Both tools are compiled in the same manner. Open the directory `V2RhoT` or `V2T` in a terminal and execute in the following order +Navigate to the base folder `./VeloDT` and execute -``` +```bash qmake make ``` +This will build the `V2RhoT` and `V2T` executables and put them in the `./VeloDT` folder. + ## License VeloDT is published under the terms of [**GNU General Public License v3.0**](./LICENSE). diff --git a/V2RhoT/README.md b/V2RhoT.md similarity index 100% rename from V2RhoT/README.md rename to V2RhoT.md diff --git a/V2T/README.md b/V2T.md similarity index 100% rename from V2T/README.md rename to V2T.md diff --git a/VeloDT.pro b/VeloDT.pro new file mode 100644 index 0000000..dc724bd --- /dev/null +++ b/VeloDT.pro @@ -0,0 +1,6 @@ +TEMPLATE = subdirs +CONFIG -= app_bundle +CONFIG += ordered +SUBDIRS += src/common src/V2RhoT src/V2T +V2RhoT.depends = common +V2T.depends = common diff --git a/common/PointClasses.h b/common/PointClasses.h deleted file mode 100644 index de20f7f..0000000 --- a/common/PointClasses.h +++ /dev/null @@ -1,104 +0,0 @@ -/******************************************************************************* -* Copyright (C) 2017 by Christian Meeßen * -* * -* This file is part of VeloDT. * -* * -* VeloDT is free software: you can redistribute it and/or modify * -* it under the terms of the GNU General Public License as published by * -* the Free Software Foundation version 3 of the License. * -* * -* VeloDT is distributed in the hope that it will be useful, but * -* WITHOUT ANY WARRANTY; without even the implied warranty of * -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * -* General Public License for more details. * -* * -* You should have received a copy of the GNU General Public License * -* along with VeloDT. If not, see . * -*******************************************************************************/ -#include -#ifndef POINTCLASSES_H_ -#define POINTCLASSES_H_ - -class Point2D { -private: - QList Vals; -public: - Point2D(); - Point2D(double x, double y); - // Assigning properties - void setX(double x){Vals[0]=x;} - void setY(double y){Vals[1]=y;} - // Retrieving properties - double x(){return Vals[0];} - double y(){return Vals[1];} - // Operators - double &operator[](int idx); -}; - -class Point3D { - private: - QList Vals; - public: - Point3D(); - Point3D(double x, double y, double z); - // Assigning properties - void setX(double x){Vals[0]=x;} - void setY(double y){Vals[1]=y;} - void setZ(double z){Vals[2]=z;} - // Retrieving properties - double x(){return Vals[0];} - double y(){return Vals[1];} - double z(){return Vals[2];} - // Operators - double &operator[](int idx); - Point3D &operator=(Point3D p); -}; - -class Point4D { - private: - QList Vals; - public: - Point4D(); - Point4D(double x, double y, double z, double v); - // Assigning properties - void setX(double x){Vals[0]=x;} - void setY(double y){Vals[1]=y;} - void setZ(double z){Vals[2]=z;} - void setV(double v){Vals[3]=v;} - // Retrieving properties - double x(){return Vals[0];} - double y(){return Vals[1];} - double z(){return Vals[2];} - double v(){return Vals[3];} - double *p(int i){return &Vals[i];} - // Operators - double &operator[](int idx); - Point4D &operator=(Point4D p); -}; - - -class Point5D { - private: - QList Vals; - public: - Point5D(); - Point5D(double x, double y, double z, double v, double prop); - // Assigning properties - void setX(double x){Vals[0]=x;} - void setY(double y){Vals[1]=y;} - void setZ(double z){Vals[2]=z;} - void setV(double v){Vals[3]=v;} - void setProp(double prop){Vals[4]=prop;} - // Retrieving properties - double x(){return Vals[0];} - double y(){return Vals[1];} - double z(){return Vals[2];} - double v(){return Vals[3];} - double prop(){return Vals[4];} - - double &operator[](int idx); - Point5D &operator=(Point5D p); - double *p(int i){return &Vals[i];} -}; - -#endif // POINTCLASSES_H_ diff --git a/V2T/docs/constants.png b/docs/constants.png similarity index 100% rename from V2T/docs/constants.png rename to docs/constants.png diff --git a/V2T/docs/eqn1.png b/docs/eqn1.png similarity index 100% rename from V2T/docs/eqn1.png rename to docs/eqn1.png diff --git a/V2T/docs/eqn2.png b/docs/eqn2.png similarity index 100% rename from V2T/docs/eqn2.png rename to docs/eqn2.png diff --git a/V2T/docs/eqn3.png b/docs/eqn3.png similarity index 100% rename from V2T/docs/eqn3.png rename to docs/eqn3.png diff --git a/V2T/docs/eqn4.png b/docs/eqn4.png similarity index 100% rename from V2T/docs/eqn4.png rename to docs/eqn4.png diff --git a/V2T/docs/eqn5.png b/docs/eqn5.png similarity index 100% rename from V2T/docs/eqn5.png rename to docs/eqn5.png diff --git a/V2T/docs/eqn6.png b/docs/eqn6.png similarity index 100% rename from V2T/docs/eqn6.png rename to docs/eqn6.png diff --git a/V2T/docs/eqn7.png b/docs/eqn7.png similarity index 100% rename from V2T/docs/eqn7.png rename to docs/eqn7.png diff --git a/V2T/docs/eqn8.png b/docs/eqn8.png similarity index 100% rename from V2T/docs/eqn8.png rename to docs/eqn8.png diff --git a/V2T/docs/eqn9.png b/docs/eqn9.png similarity index 100% rename from V2T/docs/eqn9.png rename to docs/eqn9.png diff --git a/V2T/docs/eqns.tex b/docs/eqns.tex similarity index 100% rename from V2T/docs/eqns.tex rename to docs/eqns.tex diff --git a/V2RhoT/MineraldRhodT.h b/include/V2RhoT/MineraldRhodT.h similarity index 96% rename from V2RhoT/MineraldRhodT.h rename to include/V2RhoT/MineraldRhodT.h index e15f2b6..85a3868 100644 --- a/V2RhoT/MineraldRhodT.h +++ b/include/V2RhoT/MineraldRhodT.h @@ -32,7 +32,6 @@ This class hosts the mineral property dRho/dT It must be calculated numerically to avoid errors far away from the reference temperature. The values will be stored in arrays. **/ -private: // Variables double Tmin, Tmax; QList > vals; @@ -43,10 +42,10 @@ temperature. The values will be stored in arrays. double Gnt); void fill(); // Calculate dRho/dT tables -public: + public: MineraldRhodT(); bool set_AlphaMode(int mode); - double dRhodT(double T, int mineral); // Returns dRhodT at T + double dRhodT(double T, int mineral); // Returns dRhodT at T void exportTable(); }; diff --git a/V2RhoT/Rock.h b/include/V2RhoT/Rock.h similarity index 87% rename from V2RhoT/Rock.h rename to include/V2RhoT/Rock.h index 6677dce..020a5ff 100644 --- a/V2RhoT/Rock.h +++ b/include/V2RhoT/Rock.h @@ -28,12 +28,11 @@ #include "PhysicalConstants.h" class Rock { -private: // Variables bool verbose; int AlphaMode; QString MineralPropertyDB; - MineraldRhodT * dRhodT; // Table that stores dRho/dT(T) + MineraldRhodT * dRhodT; // Table that stores dRho/dT(T) QList Composition; QList minerals_rho; QList minerals_rhoXFe; @@ -64,7 +63,7 @@ class Rock { // Anelasticity parameters QString Qmode; - double c_a, c_A, c_H, c_V, c_omega; + double c_a, c_A, c_H, c_V, c_omega, c_frequency; // P/T independent parameters for calculation of d/dT for s-waves double rock_dmudT; @@ -99,22 +98,22 @@ class Rock { void Vsyn_PT(QString VelType); void dVdTsyn_PT(QString VelType); -public: + public: Rock(); ~Rock(); void printProperties(); void printComposition(); void writedRdT(); - void setVerbose(bool val){verbose=val;} + void setVerbose(bool val) {verbose = val;} // Defining properties bool set_AlphaMode(int mode); bool set_MineralPropertyDB(QString db); bool set_XFe(double val); - void set_T0(double val){c_T0=val;} - void set_P0(double val){c_P0=val;} + void set_T0(double val) {c_T0 = val;} + void set_P0(double val) {c_P0 = val;} void set_omega(QString VelType); - void set_omega(double f){c_omega=2.*M_PI*f;} + void set_omega(double f); void set_Comp_init(double Ol, double Opx, double Cpx, double Sp, double Gnt); void set_Comp(double Ol, double Opx, double Cpx, double Sp, double Gnt); @@ -142,17 +141,18 @@ class Rock { bool setQ(int mode); // Obtaining values - int get_AlphaMode(){return AlphaMode;} + int get_AlphaMode() {return AlphaMode;} QString get_AlphaModeStr(); - QString get_MineralPropertyDB(){return MineralPropertyDB;} - double get_Vsyn_PT(){return rock_Vsyn_PT;} - double get_dVdTsyn_PT(){return rock_dVdTsyn_PT;} - double getComposition(int idx){return Composition[idx];} - double getRho(){return rock_rho_PT;} - double getOmega(){return c_omega;} - QString getQ(){return Qmode;} - double getXFe(){return rock_XFe;} - bool CustomComposition(){return UseCustomComposition;} + QString get_MineralPropertyDB() {return MineralPropertyDB;} + double get_Vsyn_PT() {return rock_Vsyn_PT;} + double get_dVdTsyn_PT() {return rock_dVdTsyn_PT;} + double getComposition(int idx) {return Composition[idx];} + double getRho() {return rock_rho_PT;} + double getOmega() {return c_omega;} + double get_frequency() {return c_frequency;} + QString getQ() {return Qmode;} + double getXFe() {return rock_XFe;} + bool CustomComposition() {return UseCustomComposition;} }; #endif // ROCK_H_ diff --git a/V2RhoT/V2RhoT.h b/include/V2RhoT/V2RhoT.h similarity index 82% rename from V2RhoT/V2RhoT.h rename to include/V2RhoT/V2RhoT.h index ff49244..39b3978 100644 --- a/V2RhoT/V2RhoT.h +++ b/include/V2RhoT/V2RhoT.h @@ -19,6 +19,7 @@ #define V2RHOT_H_ #include +#include #include #include #include //exit @@ -53,7 +54,7 @@ class V2RhoT { bool verbose; // True for debugging bool petrel; // Output Petrel points with attributes Rock * MantleRock; // The object that hosts the rock properties - EarthReferenceModel * ERM; // Calculates pressure from an ERM + EarthReferenceModel * ERM; // Calculates pressure from an ERM // Input data properties - 1: data_Vs, 2: t_crust, 3:z_topo double x_min1, x_max1, y_min1, y_max1, z_min1, z_max1, x_min2, x_max2, @@ -65,25 +66,24 @@ class V2RhoT { QList data_V; QList data_T; - public: - V2RhoT(); - ~V2RhoT(); - QString FileIn(){return File_In;} - QString FileOut(){return File_Out;} - bool readFile(QString InName, QString InType); - bool saveFile(QString OutName); - void readArgs(int &argc, char *argv[]); - void usage(); - bool Iterate(); - void Info(); + bool SetPMethod(QString method); + double pressure(double x, double y, double z); + double pressure_crust(double x, double y, double z); + double pressure_simple(double z); + void argsError(QString val, bool ok); + void help(); - private: - bool SetPMethod(QString method); - double pressure(double x, double y, double z); - double pressure_crust(double x, double y, double z); - double pressure_simple(double z); - void argsError(QString val, bool ok); - void help(); + public: + V2RhoT(); + ~V2RhoT(); + QString FileIn() {return File_In;} + QString FileOut() {return File_Out;} + bool readFile(QString InName, QString InType); + bool saveFile(QString OutName); + void readArgs(int &argc, char *argv[]); + void usage(); + bool Iterate(); + void Info(); }; #endif //V2RHOT_H_ diff --git a/V2T/V2T.h b/include/V2T/V2T.h similarity index 79% rename from V2T/V2T.h rename to include/V2T/V2T.h index df64508..cd3d725 100644 --- a/V2T/V2T.h +++ b/include/V2T/V2T.h @@ -76,32 +76,31 @@ class V2T { QList ERMz; // List for depth values of the Earth reference model QList ERMrho; // List for density values of the ERM - public: - V2T(); - bool newton(); - bool readFile(QString InName, QString InType); - bool saveFile(QString OutName); - void readArgs(int &argc, char *argv[]); - bool test_data(); - void Info(); - void usage(); + void usage_extended(); + void argsError(QString val, bool ok); + bool SetPMethod(QString method); + double ftheta(double VsS, double P, double T); + double dfdtheta(double P, double T); + double pressure(double x, double y, double z); + double pressure_simple(double z); + double pressure_crust(double x, double y, double z); + void WRITE_P(QString method); - QString FileIn(){return File_In;} - QString FileZTopo(){return File_z_topo;} - QString FileTCrust(){return File_t_crust;} - QString FileOut(){return File_Out;} - bool UseCrust(){return use_t_crust;} + public: + V2T(); + bool newton(); + bool readFile(QString InName, QString InType); + bool saveFile(QString OutName); + void readArgs(int &argc, char *argv[]); + bool test_data(); + void Info(); + void usage(); - private: - void usage_extended(); - void argsError(QString val, bool ok); - bool SetPMethod(QString method); - double ftheta(double VsS, double P, double T); - double dfdtheta(double P, double T); - double pressure(double x, double y, double z); - double pressure_simple(double z); - double pressure_crust(double x, double y, double z); - void WRITE_P(QString method); + QString FileIn() {return File_In;} + QString FileZTopo() {return File_z_topo;} + QString FileTCrust() {return File_t_crust;} + QString FileOut() {return File_Out;} + bool UseCrust() {return use_t_crust;} }; #endif // V2T_H_ diff --git a/common/ANSIICodes.h b/include/common/ANSIICodes.h similarity index 100% rename from common/ANSIICodes.h rename to include/common/ANSIICodes.h diff --git a/common/ERMs.h b/include/common/ERMs.h similarity index 86% rename from common/ERMs.h rename to include/common/ERMs.h index b4491bc..b34f382 100644 --- a/common/ERMs.h +++ b/include/common/ERMs.h @@ -33,17 +33,16 @@ class EarthReferenceModel { QList ERMz; QList ERMrho; QString ERMtype; - private: - bool INIT_AK135(); - bool INIT_PREM(); - public: - EarthReferenceModel(); - EarthReferenceModel(QString type); - bool set(QString type); - double pressure(double z); - QString type(){return ERMtype;} - bool writeP(); - bool writeP(double dz); + bool INIT_AK135(); + bool INIT_PREM(); + public: + EarthReferenceModel(); + EarthReferenceModel(QString type); + bool set(QString type); + double pressure(double z); + QString type() {return ERMtype;} + bool writeP(); + bool writeP(double dz); }; #endif // ERMS_H_ diff --git a/common/PhysicalConstants.h b/include/common/PhysicalConstants.h similarity index 100% rename from common/PhysicalConstants.h rename to include/common/PhysicalConstants.h diff --git a/include/common/PointClasses.h b/include/common/PointClasses.h new file mode 100644 index 0000000..bed7e58 --- /dev/null +++ b/include/common/PointClasses.h @@ -0,0 +1,100 @@ +/******************************************************************************* +* Copyright (C) 2017 by Christian Meeßen * +* * +* This file is part of VeloDT. * +* * +* VeloDT is free software: you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation version 3 of the License. * +* * +* VeloDT is distributed in the hope that it will be useful, but * +* WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * +* General Public License for more details. * +* * +* You should have received a copy of the GNU General Public License * +* along with VeloDT. If not, see . * +*******************************************************************************/ +#include +#ifndef POINTCLASSES_H_ +#define POINTCLASSES_H_ + +class Point2D { + QList Vals; + public: + Point2D(); + Point2D(double x, double y); + // Assigning properties + void setX(double x) {Vals[0] = x;} + void setY(double y) {Vals[1] = y;} + // Retrieving properties + double x() {return Vals[0];} + double y() {return Vals[1];} + // Operators + double &operator[](int idx); +}; + +class Point3D { + QList Vals; + public: + Point3D(); + Point3D(double x, double y, double z); + // Assigning properties + void setX(double x) {Vals[0] = x;} + void setY(double y) {Vals[1] = y;} + void setZ(double z) {Vals[2] = z;} + // Retrieving properties + double x() {return Vals[0];} + double y() {return Vals[1];} + double z() {return Vals[2];} + // Operators + double &operator[](int idx); + Point3D &operator=(Point3D p); +}; + +class Point4D { + QList Vals; + public: + Point4D(); + Point4D(double x, double y, double z, double v); + // Assigning properties + void setX(double x) {Vals[0] = x;} + void setY(double y) {Vals[1] = y;} + void setZ(double z) {Vals[2] = z;} + void setV(double v) {Vals[3] = v;} + // Retrieving properties + double x() {return Vals[0];} + double y() {return Vals[1];} + double z() {return Vals[2];} + double v() {return Vals[3];} + double *p(int i) {return &Vals[i];} + // Operators + double &operator[](int idx); + Point4D &operator=(Point4D p); +}; + + +class Point5D { + QList Vals; + public: + Point5D(); + Point5D(double x, double y, double z, double v, double prop); + // Assigning properties + void setX(double x) {Vals[0] = x;} + void setY(double y) {Vals[1] = y;} + void setZ(double z) {Vals[2] = z;} + void setV(double v) {Vals[3] = v;} + void setProp(double prop) {Vals[4] = prop;} + // Retrieving properties + double x() {return Vals[0];} + double y() {return Vals[1];} + double z() {return Vals[2];} + double v() {return Vals[3];} + double prop() {return Vals[4];} + + double &operator[](int idx); + Point5D &operator=(Point5D p); + double *p(int i){return &Vals[i];} +}; + +#endif // POINTCLASSES_H_ diff --git a/V2RhoT/MineraldRhodT.cpp b/src/V2RhoT/MineraldRhodT.cpp similarity index 83% rename from V2RhoT/MineraldRhodT.cpp rename to src/V2RhoT/MineraldRhodT.cpp index cc72ff3..36552ca 100644 --- a/V2RhoT/MineraldRhodT.cpp +++ b/src/V2RhoT/MineraldRhodT.cpp @@ -17,7 +17,8 @@ *******************************************************************************/ #include "MineraldRhodT.h" -using namespace std; +using std::cout; +using std::endl; MineraldRhodT::MineraldRhodT() { AlphaMode = 0; @@ -32,7 +33,7 @@ MineraldRhodT::MineraldRhodT() { } bool MineraldRhodT::set_AlphaMode(int mode) { - switch(mode) { + switch (mode) { case 0: // Alpha = const. AlphaMode = mode; @@ -53,25 +54,25 @@ bool MineraldRhodT::set_AlphaMode(int mode) { void MineraldRhodT::set_alpha(int idx, double Ol, double Opx, double Cpx, double Sp, double Gnt) { - if(idx==0) { + if (idx == 0) { alpha0.append(Ol); alpha0.append(Opx); alpha0.append(Cpx); alpha0.append(Sp); alpha0.append(Gnt); - } else if(idx==1) { + } else if (idx == 1) { alpha1.append(Ol); alpha1.append(Opx); alpha1.append(Cpx); alpha1.append(Sp); alpha1.append(Gnt); - } else if(idx==2) { + } else if (idx == 2) { alpha2.append(Ol); alpha2.append(Opx); alpha2.append(Cpx); alpha2.append(Sp); alpha2.append(Gnt); - } else if(idx==3) { + } else if (idx == 3) { alpha3.append(Ol); alpha3.append(Opx); alpha3.append(Cpx); @@ -106,29 +107,29 @@ void MineraldRhodT::fill() { // Fill values for every mineral after Goes et al (2000) // Initial values are for T=0degC=273K - Rho_i0.append(3222.); //Ol - Rho_i0.append(3198.); //Opx - Rho_i0.append(3280.); //Cpx - Rho_i0.append(3578.); //Sp - Rho_i0.append(3565.); //Gnt + Rho_i0.append(3222.0); //Ol + Rho_i0.append(3198.0); //Opx + Rho_i0.append(3280.0); //Cpx + Rho_i0.append(3578.0); //Sp + Rho_i0.append(3565.0); //Gnt - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { Rho_i1.append(0.); alpha_i.append(0.); dRhodT_i.append(0.); vals_T.append(0.); } - vals_T.append(0.); // vals_T needs one additional item for the temperature + vals_T.append(0.); // vals_T needs one additional item for the temperature // Calculate - for(int i=0; i(i); vals_T[0] = T_i; - for(int j=0; j<5; j++) { - if(AlphaMode==0) { + for (int j=0; j < 5; j++) { + if (AlphaMode == 0) { alpha_i[j] = alpha0[j]; - } else if(AlphaMode==1) { + } else if (AlphaMode == 1) { alpha_i[j] = alpha0[j] + alpha1[j]*T_i + alpha2[j]/T_i + alpha3[j]/T_i/T_i; } @@ -159,14 +160,14 @@ double MineraldRhodT::dRhodT(double T, int mineral) { cout << PRINT_WARNING "Iteration failed. Temperature out of bounds T=" << T << " K\n"; return 999; - } else if ( fmod(T,1) == 0 && T < 2273.) { + } else if (fmod(T, 1) == 0 && T < 2273.0) { // If temperature is already in table - return vals[(int) T - 273][mineral+1]; + return vals[static_cast(T) - 273][mineral+1]; } else { - T_low = (int) (T - fmod(T,1) - 273); + T_low = static_cast(T - fmod(T, 1) - 273); dRdT_low = vals[T_low][mineral+1]; dRdT_high = vals[T_low+1][mineral+1]; - return (T - (double) T_low)*(dRdT_high - dRdT_low) + dRdT_low; + return (T - static_cast(T_low))*(dRdT_high - dRdT_low) + dRdT_low; } } @@ -177,7 +178,7 @@ void MineraldRhodT::exportTable() { cout << "Writing dRho/dT to " << OutName.toUtf8().data() << endl; QFile tmp(OutName); - if(!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { + if (!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() << endl; exit(1); @@ -189,18 +190,18 @@ void MineraldRhodT::exportTable() { fout.setRealNumberNotation(QTextStream::FixedNotation); fout << T_header.toUtf8().data() << endl; - for(int i=0; i object) { cout << left << setw(width) << setfill(' ') << title.toUtf8().data(); cout << left << setw(10) << setfill(' ') << unit.toUtf8().data(); - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { cout << right << setw(10) << setfill(' ') << object[i]; } cout << endl; @@ -386,23 +396,37 @@ void Rock::printline(int width, QString title, QString unit, void Rock::printComposition() { QList header; - int width = 15; + QList info_output; + const unsigned int width_col_1 = 12; + const unsigned int width_col_2 = 5; + const unsigned int width_col_c = 8; + const unsigned int width_tot = width_col_1 + width_col_2 + 6*width_col_c; + + info_output << Composition; + info_output << rock_XFe; + header << ("Ol"); header << ("Opx"); header << ("Cpx"); header << ("Sp"); header << ("Gnt"); + header << ("XFe"); - cout << endl; - cout << left << setw(width) << setfill(' ') << "Property"; - cout << left << setw(10) << setfill(' ') << "Unit"; - for(int i=0; i<5; i++) { - cout << right << setw(10) << setfill(' ') << header[i].toUtf8().data(); + cout << endl + << left << setw(width_col_1) << setfill(' ') << "Property" + << left << setw(width_col_2) << setfill(' ') << "Unit"; + for (int i=0; i < 6; i++) { + cout << right << setw(width_col_c) << setfill(' ') + << header[i].toUtf8().data(); } - cout << endl; - cout << setw(75) << setfill('-') << "" << endl; - printline(width, "Composition", "rel.", Composition); - cout << endl; + cout << endl + << setw(width_tot) << setfill('-') << "" << endl + << left << setw(width_col_1) << setfill(' ') << "Composition" + << left << setw(width_col_2) << setfill(' ') << "rel."; + for (int i = 0; i < 6; i++) { + cout << right << setw(width_col_c) << setfill(' ') << info_output[i]; + } + cout << endl << endl; } void Rock::writedRdT() { @@ -412,7 +436,7 @@ void Rock::writedRdT() { void Rock::printProperties() { QList header; - int width=15; + const int width = 15; header << ("Ol"); header << ("Opx"); header << ("Cpx"); @@ -425,7 +449,7 @@ void Rock::printProperties() { << endl << left << setw(width) << setfill(' ') << "Property" << left << setw(10) << setfill(' ') << "Unit"; - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { cout << right << setw(10) << setfill(' ') << header[i].toUtf8().data(); } cout << endl; @@ -452,23 +476,22 @@ void Rock::printProperties() { void Rock::alpha() { // Calculates alpha(rock_T) - if(verbose){ + if (verbose) { cout << endl << "Calculate alpha(T)" << endl << "> AlphaMode " << AlphaMode << endl; } - switch(AlphaMode) { + switch (AlphaMode) { case 0: - if(verbose){cout << "> alpha=const." << endl;} - for(int i=0; i<5; i++) { + if (verbose) cout << "> alpha=const." << endl; + for (int i=0; i < 5; i++) minerals_alpha_T[i] = minerals_alpha0[i]; - } break; case 1: - if(verbose){cout << "> alpha(T)" << endl;} - for(int i=0; i<5; i++) { + if (verbose) cout << "> alpha(T)" << endl; + for (int i=0; i < 5; i++) minerals_alpha_T[i] = minerals_alpha0[i] + minerals_alpha1[i]*rock_T - + minerals_alpha2[i]/rock_T + minerals_alpha3[i]/rock_T/rock_T; - } + + minerals_alpha2[i]/rock_T + + minerals_alpha3[i]/rock_T/rock_T; break; } } @@ -480,14 +503,14 @@ bool Rock::set_XFe(double val) { } void Rock::rho_PT() { - //Sets VRH average rock density + // Sets VRH average rock density QList rho_minerals_PT; double rho_avrg_PT; // Calculate rho(P,T) for each mineral rho_avrg_PT = 0.0; - if(verbose){cout << endl << "Calculate rho_PT" << endl;} - for(int i=0; i<5; i++) { - if(verbose){ + if (verbose) cout << endl << "Calculate rho_PT" << endl; + for (int i=0; i < 5; i++) { + if (verbose) { cout << "Mineral index i " << i << endl << "> minerals_rhoXFe[i] " << minerals_rhoXFe[i] << endl << "> minerals_alpha_T[i] " << minerals_alpha_T[i] << endl @@ -500,29 +523,28 @@ void Rock::rho_PT() { << "> rho_minerals_PT[i] " << rho_minerals_PT[i] << endl; } rho_minerals_PT.append(minerals_rhoXFe[i]*(1. - minerals_alpha_T[i]*(rock_T - - c_T0) + (rock_P - c_P0)/minerals_K_PT[i])); + - c_T0) + (rock_P - c_P0)/minerals_K_PT[i])); rho_avrg_PT = rho_avrg_PT + Composition[i]*rho_minerals_PT[i]; } rock_rho_PT = rho_avrg_PT; } void Rock::mu_PT() { - //Returns VRH averaged rock property + // Returns VRH averaged rock property double mu_voigt, mu_reuss; - if(verbose){cout << endl << "Calculate mu(P,T)\n";} + if (verbose) cout << endl << "Calculate mu(P,T)\n"; // Calculate mu(P,T) for each mineral - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { minerals_mu_PT[i] = minerals_mu[i] + (rock_T - c_T0)*minerals_dmudT[i] + (rock_P - c_P0)*minerals_dmudP[i] + rock_XFe*minerals_dmudX[i]; - if(verbose){ + if (verbose) cout << "> mu[" << i << "]: " << minerals_mu_PT[i] << endl; - } } // Calculate mu_voigt and mu_reuss mu_voigt = 0.0; mu_reuss = 0.0; - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { mu_voigt = mu_voigt + Composition[i]*minerals_mu_PT[i]; mu_reuss = mu_reuss + Composition[i]/minerals_mu_PT[i]; } @@ -535,23 +557,22 @@ void Rock::K_PT() { // Returns VRH averaged rock property double K_voigt, K_reuss; - if(verbose){cout << endl << "Calculate K(P,T)\n";} + if (verbose) cout << endl << "Calculate K(P,T)\n"; // Calculate K(P,T) for each mineral - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { minerals_K_PT[i] = minerals_K[i] + (rock_T - c_T0)*minerals_dKdT[i] + (rock_P - c_P0)*(minerals_dKdP[i] + rock_XFe*minerals_dKdPdX[i]) + rock_XFe*minerals_dKdX[i]; - if(verbose){ + if (verbose) cout << "K[" << i << "]: " << minerals_K_PT[i] << endl; - } } // Calculate K_voigt and K_reuss K_voigt = 0.0; K_reuss = 0.0; - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { K_voigt = K_voigt + Composition[i]*minerals_K_PT[i]; K_reuss = K_reuss + Composition[i]/minerals_K_PT[i]; } @@ -563,30 +584,27 @@ void Rock::K_PT() { void Rock::Qmu_T() { // Calculate Q_mu, Eqn A6 double val, Estar; - if(verbose){cout << endl << "Calculate Qmu_T" << endl;} + if (verbose) cout << endl << "Calculate Qmu_T" << endl; Estar = c_H + rock_P*c_V; val = c_A*pow(c_omega, c_a)*exp(c_a*Estar/c_R/rock_T); - if(verbose){cout << "> Qmu(T) = " << val << endl;} + if (verbose) cout << "> Qmu(T) = " << val << endl; rock_Qmu_T = val; } void Rock::QP_T() { // Calculate Q_P, Eqn A7 - double val, Estar, L; - if(verbose){cout << endl << "Calculate QP_T" << endl;} - Estar = c_H + rock_P*c_V; - val = c_A*pow(c_omega,c_a)*exp(c_a*Estar/c_R/rock_T); + double L; + if (verbose) cout << endl << "Calculate QP_T" << endl; L = 4.*rock_mu_PT/(3.*rock_K_PT + 4.*rock_mu_PT); rock_QP_T = rock_Qmu_T/L; - if(verbose){cout << "> QP(T) = " << rock_QP_T << endl;} + if (verbose) cout << "> QP(T) = " << rock_QP_T << endl; } void Rock::dmudT() { double anh_mu_reuss = 0; - if(verbose){cout << endl << "Calculate dmudT" << endl;} - for(int i=0; i<5; i++) { + if (verbose) cout << endl << "Calculate dmudT" << endl; + for (int i=0; i < 5; i++) anh_mu_reuss = anh_mu_reuss + Composition[i]/minerals_mu_PT[i]; - } anh_mu_reuss = 1./anh_mu_reuss; rock_dmudT = anh_sum1 + anh_sum2/pow(anh_mu_reuss, 2); } @@ -595,16 +613,17 @@ void Rock::dMdT() { // Calculate d/dT (P,T) for P-waves, Eqn. A5b double M_reuss = 0.0; anh_sum2 = 0.0; - if(verbose){cout << "Calculate d/dT\n";} - for(int i=0; i<5; i++) { - M_reuss = M_reuss + Composition[i]/(minerals_K_PT[i] + + 4./3*minerals_mu_PT[i]); + if (verbose) cout << "Calculate d/dT" << endl; + for (int i=0; i < 5; i++) { + M_reuss = M_reuss + Composition[i]/(minerals_K_PT[i] + + 4.0/3.0*minerals_mu_PT[i]); anh_sum2 = anh_sum2 + (Composition[i]/(minerals_K_PT[i] + 4./3*minerals_mu_PT[i])*(minerals_dKdT[i] + 4./3*minerals_dmudT[i])); } M_reuss = 1./M_reuss; - if(verbose) { + if (verbose) { cout << "M_reuss: " << M_reuss << endl << "anh_sum1: " << anh_sum1 << endl << "anh_sum2: " << anh_sum2 << endl; @@ -617,17 +636,17 @@ void Rock::Vsyn_PT(QString VelType) { // Calculate synthethic velocity, Eqn A8 double Vanh, Vanel; - if(verbose){cout << endl << "Calculate Vsyn_PT:\n";} - if(VelType=="S") { + if (verbose) cout << endl << "Calculate Vsyn_PT:\n"; + if (VelType == "S") { Vanh = sqrt(rock_mu_PT/rock_rho_PT); Vanel = 1. - 2./rock_Qmu_T/tan(M_PI*c_a/2.); } else { - if(verbose){cout << "VelType: P\n";} + if (verbose) cout << "VelType: P\n"; Vanh = sqrt((rock_K_PT + 4./3*rock_mu_PT)/rock_rho_PT); Vanel = 1. - 2./rock_QP_T/tan(M_PI*c_a/2.); } - if(verbose) { + if (verbose) { cout << "rock_K_PT: " << rock_K_PT << endl << "rock_mu_PT: " << rock_mu_PT << endl << "rock_rho_PT: " << rock_rho_PT << endl @@ -644,8 +663,8 @@ void Rock::drhodT_T() { QList mineral_drhodT_T; result = 0.; - for(int i=0; i<5; i++) { - mineral_drhodT_T.append(dRhodT->dRhodT(rock_T,i)); + for (int i=0; i < 5; i++) { + mineral_drhodT_T.append(dRhodT->dRhodT(rock_T, i)); // Calculate average drho/dT for rock result = result + Composition[i]*mineral_drhodT_T[i]; } @@ -655,16 +674,17 @@ void Rock::drhodT_T() { void Rock::dVdTsyn_PT(QString VelType) { // Calculate dV/dT_syn, Eqn. A9 (anel) and A4 (anh) double dVdTanel, dVdTanh; - if(VelType=="S") { + if (VelType == "S") { dVdTanel = c_A*c_H/rock_Qmu_T/2./c_R/rock_T/rock_T/tan(M_PI*c_a/2.); dVdTanh = (rock_dmudT - pow(rock_Vsyn_PT, 2)*rock_drhodT_T) /(2.*rock_rho_PT*rock_Vsyn_PT); } else { dVdTanel = c_A*c_H/rock_QP_T/2./c_R/rock_T/rock_T/tan(M_PI*c_a/2.); - dVdTanh = (rock_dMdT_PT - pow(rock_Vsyn_PT,2)*rock_drhodT_T)/(2.*rock_rho_PT*rock_Vsyn_PT); + dVdTanh = (rock_dMdT_PT - pow(rock_Vsyn_PT, 2)*rock_drhodT_T)/ + (2.*rock_rho_PT*rock_Vsyn_PT); } - if(verbose) { + if (verbose) { cout << "Calculate dV/dT_syn\n" << "dVdTanel: " << dVdTanel << endl << "dVdTanh: " << dVdTanh << endl; @@ -684,17 +704,20 @@ bool Rock::calc_prop_PT(double Pressure, double Temperature, QString VelType) { mu_PT(); Qmu_T(); - if(VelType=="P"){ QP_T();} + if (VelType == "P") + QP_T(); Vsyn_PT(VelType); - if(VelType=="S"){ dmudT();} - else { dMdT();} + if (VelType == "S") + dmudT(); + else + dMdT(); drhodT_T(); dVdTsyn_PT(VelType); - if(verbose) { + if (verbose) { cout << endl << "Calculated properties:" << endl << "rock_T: " << rock_T << endl @@ -716,16 +739,16 @@ bool Rock::calc_prop_PT(double Pressure, double Temperature, QString VelType) { bool Rock::calc_prop(QString VelType) { // Calculates P/T independent parameters for Eqn A5b - if(verbose){ + if (verbose) { cout << "----------------------------------------------" << endl << "Calculating P/T independent parameters" << endl << "----------------------------------------------" << endl; } // Mineral density including iron content - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) { minerals_rhoXFe[i] = minerals_rho[i] + minerals_drhodX[i]*rock_XFe; - if(verbose){ + if (verbose) { cout << "minerals_rho[i] " << minerals_rho[i] << endl << "minerals_drhodX[i] " << minerals_drhodX[i] << endl << "rock_XFe " << rock_XFe << endl @@ -734,16 +757,15 @@ bool Rock::calc_prop(QString VelType) { } anh_sum1 = 0.; - if(VelType=="S") { - for(int i=0; i<5; i++) { + if (VelType == "S") { + for (int i=0; i < 5; i++) anh_sum1 = anh_sum1 + Composition[i]*minerals_dmudT[i]; - } } else { - for(int i=0; i<5; i++) { + for (int i=0; i < 5; i++) anh_sum1 = anh_sum1 + Composition[i]*(minerals_dKdT[i] - + 4./3.*minerals_dmudT[i]); - } + + 4.0/3.0*minerals_dmudT[i]); } - if(verbose){cout << "----------------------------------------------" << endl;} + if (verbose) + cout << "----------------------------------------------" << endl; return true; } diff --git a/V2RhoT/V2RhoT.cpp b/src/V2RhoT/V2RhoT.cpp similarity index 82% rename from V2RhoT/V2RhoT.cpp rename to src/V2RhoT/V2RhoT.cpp index 6e42a28..68da753 100644 --- a/V2RhoT/V2RhoT.cpp +++ b/src/V2RhoT/V2RhoT.cpp @@ -20,7 +20,8 @@ #define _USE_MATH_DEFINES #define EOL "\n" -using namespace std; +using std::cout; +using std::endl; const QString compilationTime = QString("%1 %2").arg(__DATE__).arg(__TIME__); @@ -70,7 +71,7 @@ void V2RhoT::Info() { << "Input\n" << "-----\n" << "Vs : " << File_In.toUtf8().data() << "\n"; - if(use_t_crust){ + if (use_t_crust) { cout << "Topography : " << File_z_topo.toUtf8().data() << "\n" << "Crustal thickness : " << File_t_crust.toUtf8().data() << "\n"; } @@ -79,7 +80,7 @@ void V2RhoT::Info() { << "------\n" << "Temperatures : " << File_Out.toUtf8().data() << "\n" << endl; - if(use_t_crust) { + if (use_t_crust) { cout << "Densities\n" << "---------\n" << "Crust : "<< rho_crust << " kg/m3\n" @@ -89,15 +90,15 @@ void V2RhoT::Info() { } cout << "Anelasticity parameters\n" << "-----------------------\n" - << "Q : " << MantleRock->getQ().toUtf8().data() << endl + << "Q : " << MantleRock->getQ().toUtf8().data() << endl << endl; - if(MantleRock->CustomComposition() && !verbose) { + if (MantleRock->CustomComposition() && !verbose) { cout << "Using custom rock composition\n" << "-----------------------------\n"; MantleRock->printComposition(); } - if(verbose){ + if (verbose) { cout << "Mantle rock composition" << endl << "-----------------------" << endl; MantleRock->printComposition(); @@ -106,7 +107,7 @@ void V2RhoT::Info() { cout << "Other\n" << "-----\n" << "Velocity type : " << VelType.toUtf8().data() << endl - << "Omega : " << MantleRock->getOmega()/2./M_PI << " Hz\n" + << "Frequency : " << MantleRock->get_frequency() << " Hz\n" << "Verbose : " << (verbose ? "true" : "false") << endl << "z-scaling factor : " << scaleZ << endl << "V-scaling factor : " << scaleVs << endl @@ -183,7 +184,7 @@ void V2RhoT::help() { } void V2RhoT::argsError(QString val, bool ok) { - if(!ok) { + if (!ok) { cout << PRINT_ERROR "Invalid value in " << val.toUtf8().data() << endl; exit(1); } @@ -206,133 +207,134 @@ void V2RhoT::readArgs(int &argc, char *argv[]) { definedPMethod = false; QStringList arg; - for(int i=0; i1) && (arg[1]=="-h" || arg[1]=="-help") ) { + if ((argc > 1) && (arg[1] == "-h" || arg[1] == "-help")) { help(); - } else if( (argc>1) && (arg[1]=="-prop")) { + } else if ((argc > 1) && (arg[1] == "-prop")) { MantleRock->printProperties(); - } else if( (argc>1) && (arg[1]=="-writedRdT")) { + } else if ((argc > 1) && (arg[1] == "-writedRdT")) { MantleRock->writedRdT(); - } else if( argc>1 && arg[1]=="--WRITEP_AK135") { + } else if (argc > 1 && arg[1] == "--WRITEP_AK135") { ERM->set("AK135"); ok = ERM->writeP(); argsError(arg[1], ok); exit(0); - } else if( argc>1 && arg[1]=="--BUILD") { + } else if (argc > 1 && arg[1] == "--BUILD") { cout << endl << "Build date: " << compilationTime.toUtf8().data() << endl << endl; exit(0); - }else if(argc<4) { + } else if (argc < 4) { cout << endl << PRINT_ERROR "Not enough arguments!"; usage(); } else { File_In = arg[1].toUtf8().data(); File_Out = arg[2].toUtf8().data(); - for(int i=3; iset_AlphaMode(1); - } else if(arg[i]=="-compc") { - for(int j=1; j<6; j++) { + } else if (arg[i] == "-compc") { + for (int j=1; j < 6; j++) { comp_input.append(arg[i+j].toDouble(&ok)); argsError(arg[i+j], ok); } MantleRock->set_Comp(comp_input[0], comp_input[1], comp_input[2], comp_input[3], comp_input[4]); i+=5; - } else if(arg[i]=="-compp") { + } else if (arg[i] == "-compp") { MantleRock->set_Comp(arg[i+1].toDouble(&ok)); argsError(arg[i+1], ok); i++; - } else if(arg[i]=="-ERM") { + } else if (arg[i] == "-ERM") { ok = SetPMethod(arg[i+1]); argsError(arg[i], ok); definedPMethod = true; i++; - } else if(arg[i]=="-f") { + } else if (arg[i] == "-f") { UseCustomOmega = true; CustomFreq = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-fdamp") { + } else if (arg[i] == "-fdamp") { c_Fdamp = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-minDB"){ + } else if (arg[i] == "-minDB") { int DBarg = arg[i+1].toInt(&ok); argsError(arg[i], ok); - if(DBarg==1){ + if (DBarg == 1) { ok = MantleRock->set_MineralPropertyDB("Cammarano"); - } else if(DBarg==2){ + } else if (DBarg == 2) { ok = MantleRock->set_MineralPropertyDB("Goes"); } argsError(arg[i], ok); i++; - } else if(arg[i]=="-petrel") { - petrel=true; - } else if(arg[i]=="-Q") { + } else if (arg[i] == "-petrel") { + petrel = true; + } else if (arg[i] == "-Q") { int Qarg = arg[i+1].toInt(&ok); argsError(arg[i], ok); ok = MantleRock->setQ(Qarg); argsError(arg[i], ok); i++; - } else if(arg[i]=="-rc") { + } else if (arg[i] == "-rc") { rho_crust = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-rm") { + } else if (arg[i] == "-rm") { rho_mantle = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-ra") { + } else if (arg[i] == "-ra") { rho_avrg = arg[i+1].toDouble(&ok); SetPMethod("simple"); argsError(arg[i], ok); i++; - } else if(arg[i]=="-scaleZ") { + } else if (arg[i] == "-scaleZ") { scaleZ = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-scaleV") { + } else if (arg[i] == "-scaleV") { scaleVs = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-scatter") { + } else if (arg[i] == "-scatter") { ArbitraryPoints = true; - } else if(arg[i]=="-Tstart") { + } else if (arg[i] == "-Tstart") { T_start = arg[i+1].toDouble(&ok); argsError(arg[i], ok); i++; - } else if(arg[i]=="-t_crust") { + } else if (arg[i] == "-t_crust") { File_t_crust = arg[i+1]; okCrust = true; SetPMethod("crust"); definedPMethod = true; i++; - } else if(arg[i]=="-t") { + } else if (arg[i] == "-t") { threshold = arg[i+1].toDouble(&ok); argsError(arg[i], ok); - } else if(arg[i]=="-v") { + i++; + } else if (arg[i] == "-v") { verbose = true; MantleRock->setVerbose(verbose); - } else if(arg[i]=="-xfe") { + } else if (arg[i] == "-xfe") { double ArgXFe = arg[i+1].toDouble(&ok); argsError(arg[i], ok); ok = MantleRock->set_XFe(ArgXFe); argsError(arg[i], ok); i++; - } else if(arg[i]=="-z_topo") { + } else if (arg[i] == "-z_topo") { File_z_topo = arg[i+1]; okTopo = true; SetPMethod("crust"); @@ -345,45 +347,35 @@ void V2RhoT::readArgs(int &argc, char *argv[]) { } } - // Some logical checks - //if(ArbitraryPoints) - //{ - // printf("ArbitraryPoints: TRUE\n"); - //} - //if(use_t_crust) - //{ - // printf("use_t_crust: TRUE\n"); - //} - if(!definedPMethod) { + // Some logic checks + if (!definedPMethod) { // If not defined by user initiate the default method ERM->set(PMethod); } - if(ArbitraryPoints && use_t_crust) { + if (ArbitraryPoints && use_t_crust) { cout << PRINT_ERROR "Defined crustal thickness/topography using scattered" << " data. Exit.\n"; exit(1); } - if(okTopo ^ okCrust) { + if (okTopo ^ okCrust) { cout << PRINT_ERROR "Not enough arguments:" << endl; - if(okTopo) { + if (okTopo) cout << "No crustal thickness file defined!" << endl; - } else { + else cout << "No topographic elevation file defined!" << endl; - } exit(1); } - if(use_t_crust && (okTopo || okCrust)) { + if (use_t_crust && (okTopo || okCrust)) { cout << PRINT_ERROR "Pressure calculation method set to 'crust' but "; - if(okTopo) { + if (okTopo) cout << "crustal thickness not defined!" << endl; - } else { + else cout << "topography not defined!" << endl; - } exit(1); } // Define Omega - if(UseCustomOmega) { + if (UseCustomOmega) { MantleRock->set_omega(CustomFreq); } else { // If no custom frequency given define omega according to Goes et al. (2000) @@ -394,10 +386,10 @@ void V2RhoT::readArgs(int &argc, char *argv[]) { bool V2RhoT::readFile(QString InName, QString InType) { // InType "topo", "crust" or "vox" - double x,y,z,val; + double x, y, z, val; bool okx, oky, okz, okval, okGrid; - okGrid = false; // Used to check if grid size was extracted from input voxel + okGrid = false; // Used to check if grid size was extracted from input voxel cout << "Reading file: " << InName.toUtf8().data() << endl; QFile file(InName); @@ -418,17 +410,17 @@ bool V2RhoT::readFile(QString InName, QString InType) { y_max3 = x_max1; int n = 0; - if(file.open(QIODevice::ReadOnly)) { - QTextStream stream( &file ); + if (file.open(QIODevice::ReadOnly)) { + QTextStream stream(&file); - while(!stream.atEnd()) { + while (!stream.atEnd()) { n++; QString t = stream.readLine().simplified(); - if(!t.startsWith("#") && (InType == "topo" || InType == "crust")) { + if (!t.startsWith("#") && (InType == "topo" || InType == "crust")) { /// Topography or crustal thickness QStringList vals = t.split(" "); - if(vals.count() != 6) { + if (vals.count() != 6) { file.close(); cout << PRINT_ERROR "In header of " << InName.toUtf8().data() << " - grid size." << endl; @@ -437,9 +429,9 @@ bool V2RhoT::readFile(QString InName, QString InType) { x = vals[0].toDouble(&okx); y = vals[1].toDouble(&oky); z = vals[2].toDouble(&okz); - if(okx && oky && okz) { - Point3D p3D(x,y,z); - if(InType == "crust") { + if (okx && oky && okz) { + Point3D p3D(x, y, z); + if (InType == "crust") { t_crust.append(p3D); } else if (InType == "topo") { z_topo.append(p3D); @@ -449,47 +441,47 @@ bool V2RhoT::readFile(QString InName, QString InType) { cout << PRINT_ERROR "In value conversion line " << n << endl; return false; } - //Determine min/max - if(InType == "crust") { - if(x < x_min2) { + // Determine min/max + if (InType == "crust") { + if (x < x_min2) { x_min2 = x; - } else if(x > x_max2) { + } else if (x > x_max2) { x_max2 = x; } - if(y < y_min2) { + if (y < y_min2) { y_min2 = y; - } else if(y > y_max2) { + } else if (y > y_max2) { y_max2 = y; } } else if (InType =="topo") { - if(x < x_min3) { + if (x < x_min3) { x_min3 = x; - } else if(x > x_max3) { + } else if (x > x_max3) { x_max3 = x; } - if(y < y_min3) { + if (y < y_min3) { y_min3 = y; - } else if(y > y_max3) { + } else if (y > y_max3) { y_max3 = y; } } } - } else if(t.startsWith("# Grid_size:") && InType == "vox") { + } else if (t.startsWith("# Grid_size:") && InType == "vox") { // velocity grid - t.remove(0,12); + t.remove(0, 12); QStringList gridSize = t.split("x", QString::SkipEmptyParts); nX = gridSize[0].toInt(&okx); nY = gridSize[1].toInt(&oky); nZ = gridSize[2].toInt(&okz); - if(!okx || !oky || !okz) { + if (!okx || !oky || !okz) { cout << PRINT_ERROR "In header of " << InName.toUtf8().data() << " - grid size."; exit(1); } okGrid = true; - } else if(!t.startsWith("#") && InType == "vox") { + } else if (!t.startsWith("#") && InType == "vox") { QStringList vals = t.split(" "); - if(vals.count() != 4) { + if (vals.count() != 4) { file.close(); cout << PRINT_ERROR "In line " << n << ": unkown amount of columns." << endl; @@ -499,15 +491,15 @@ bool V2RhoT::readFile(QString InName, QString InType) { y = vals[1].toDouble(&oky); z = scaleZ*vals[2].toDouble(&okz); val = scaleVs*vals[3].toDouble(&okval); - if(val < 50) { + if (val < 50) { cout << endl << endl << PRINT_WARNING "Imported velocity is < 50 m/s! Maybe imported " << "velocities are in km/s?\n" << "To convert to m/s use option -scaleV 1000\n" << endl; exit(1); - } else if(okx && oky && okz && okval) { - Point5D p5D(x,y,z,val,0.0); + } else if (okx && oky && okz && okval) { + Point5D p5D(x, y, z, val, 0.0); data_V.append(p5D); } else { file.close(); @@ -515,26 +507,26 @@ bool V2RhoT::readFile(QString InName, QString InType) { exit(1); } // Get minima and maxima - if(x < x_min1) { + if (x < x_min1) { x_min1 = x; - } else if(x > x_max1) { + } else if (x > x_max1) { x_max1 = x; } - if(y < y_min1) { + if (y < y_min1) { y_min1 = y; - } else if(y > y_max1) { + } else if (y > y_max1) { y_max1 = y; } - if(z < z_min1) { + if (z < z_min1) { z_min1 = z; - } else if(z > z_max1) { + } else if (z > z_max1) { z_max1 = z; } } } } - if(!ArbitraryPoints && !okGrid) { + if (!ArbitraryPoints && !okGrid) { cout << PRINT_WARNING "Could not find grid information. Set output to " << "scattered data.\n"; ArbitraryPoints = true; @@ -554,11 +546,10 @@ bool V2RhoT::saveFile(QString OutName) { QDateTime currentDateTime = QDateTime::currentDateTime(); QString timestamp = currentDateTime.toString(); - if(use_t_crust) { + if (use_t_crust) usecrust = "yes"; - } else { + else usecrust = "no"; - } Info_header = QString("# Transformation settings:\n"); Info_header += QString("# Date created: %1\n").arg(timestamp); @@ -571,23 +562,23 @@ bool V2RhoT::saveFile(QString OutName) { Info_header += QString("# Gnt - %1\n").arg(MantleRock->getComposition(4), 5, 'f', 2); Info_header += QString("# Iron content XFe: %1\n").arg(MantleRock->getXFe(),3,'f',2); Info_header += QString("# Pressure calculation method: %1\n").arg(PMethod); - if(use_t_crust) { + if (use_t_crust) { Info_header += QString("# Use crustal thickness for pressure calculation: %1\n").arg(usecrust); Info_header += QString("# Density crust/mantle/average: %1%2%3\n") - .arg(rho_crust,7,'f',1).arg(rho_mantle,7,'f',1).arg(rho_avrg,7,'f',1); + .arg(rho_crust, 7, 'f', 1).arg(rho_mantle, 7, 'f', 1).arg(rho_avrg, 7, 'f', 1); } Info_header += QString("# Alpha: %1\n").arg(MantleRock->get_AlphaModeStr()); - if(use_t_crust) { + if (use_t_crust) { Info_header += QString("# Topography: %1\n").arg(File_z_topo); Info_header += QString("# Crustal thickness: %1\n").arg(File_t_crust); } - Info_header += QString("# Wave frequency (Omega) / Hz: %1\n").arg(MantleRock->getOmega()/2/M_PI); + Info_header += QString("# Wave frequency / Hz: %1\n").arg(MantleRock->get_frequency()); Info_header += QString("# Dampening factor: %1\n").arg(c_Fdamp); Info_header += QString("# Iteration starting temperature / K: %1\n").arg(T_start); Info_header += QString("# Anelasticity parameters: %1\n").arg(MantleRock->getQ()); Info_header += QString("# Average iteration steps: %1\n").arg(count_avrg,0,'f',1); - if(petrel) { + if (petrel) { T_header = QString("# Petrel Points with attributes\n"); T_header += QString("# Unit in X and Y direction: m\n"); T_header += QString("# Unit in depth: m\n"); @@ -607,21 +598,21 @@ bool V2RhoT::saveFile(QString OutName) { T_header += QString("# Columns:\n"); T_header += QString("# 1 - X\n"); T_header += QString("# 2 - Y\n"); - T_header += QString("# 3 - Z [m]\n"); - T_header += QString("# 4 - V_%1 [m/s]\n").arg(VelType); - T_header += QString("# 5 - T [degC]\n"); - T_header += QString("# 6 - Rho [kg/m3]"); + T_header += QString("# 3 - Z / m\n"); + T_header += QString("# 4 - V_%1 / m/s\n").arg(VelType); + T_header += QString("# 5 - T / degC\n"); + T_header += QString("# 6 - Rho / kg/m3"); } else { T_header = QString("# Type: GMS GridPoints\n"); T_header += QString("# Version: 2\n"); T_header += QString("# Description:\n"); T_header += QString("# Format: free\n"); - T_header += QString("# Field: 1 X [m]\n"); - T_header += QString("# Field: 2 Y [m]\n"); - T_header += QString("# Field: 3 Z [m]\n"); - T_header += QString("# Field: 4 V_%1 [m/s]\n").arg(VelType); - T_header += QString("# Field: 5 T [degC]\n"); - T_header += QString("# Field: 6 Rho [kg/m3]\n"); + T_header += QString("# Field: 1 X / m\n"); + T_header += QString("# Field: 2 Y / m\n"); + T_header += QString("# Field: 3 Z / m\n"); + T_header += QString("# Field: 4 V_%1 / m/s\n").arg(VelType); + T_header += QString("# Field: 5 T / degC\n"); + T_header += QString("# Field: 6 Rho / kg/m3\n"); T_header += QString("# Projection: Local Rectangular\n"); T_header += QString("# Information from grid:\n"); T_header += QString("# Grid_size: %1 x %2 x %3\n").arg(nX).arg(nY).arg(nZ); @@ -635,7 +626,7 @@ bool V2RhoT::saveFile(QString OutName) { cout << "Writing temperature file " << OutName.toUtf8().data() << endl; QFile tmp(OutName); - if(!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { + if (!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() << endl; exit(1); @@ -646,7 +637,7 @@ bool V2RhoT::saveFile(QString OutName) { fout.setFieldAlignment(QTextStream::AlignRight); fout.setRealNumberNotation(QTextStream::FixedNotation); fout << T_header.toUtf8().data() << endl; - for(int i=0; ipressure(z); - } else if(use_t_crust) { + } else if (use_t_crust) { return pressure_crust(x, y, z); } else { return pressure_simple(z); } } -double V2RhoT::pressure_simple(double z){ +double V2RhoT::pressure_simple(double z) { // Simple lithostatic pressure with average density return rho_avrg * c_g * abs(z); } @@ -717,12 +708,14 @@ double V2RhoT::pressure_crust(double x, double y, double z) { double P_crust, t_mantle, P_mantle; // Get index in data_t_crust and data_z_topo - int i=0; - while( (t_crust[i].x()!=x) && (t_crust[i].y()!=y)) {i++;} + int i = 0; + while ((t_crust[i].x() != x) && (t_crust[i].y() != y)) + i++; i_crust = i; i = 0; - while( (z_topo[i].x()!=x) && (z_topo[i].y()!=y)) {i++;} - i_topo=i; + while ((z_topo[i].x() != x) && (z_topo[i].y() != y)) + i++; + i_topo = i; P_crust = rho_crust*c_g*t_crust[i_crust].z(); // Pressure crust t_mantle = z_topo[i_topo].z() - t_crust[i_crust].z() - z; // Mantle thickness @@ -752,20 +745,21 @@ bool V2RhoT::Iterate() { << "T_start: " << T_start << " K\n"; n_V = data_V.length(); - for(int i=0; icalc_prop(VelType); // Calculate all P/T independent rock properties + // Calculate all P/T independent rock properties + MantleRock->calc_prop(VelType); T_n = T_start; // Temperature at step n T_n1 = 0.; // Temperature at step n+1 counter = 0; deltaT = threshold + 1; - while(deltaT > threshold) { + while (deltaT > threshold) { // Calculate rock properties - if(verbose){ + if (verbose) { cout << endl << endl << "Point " << i+1 << ", Step " << counter << endl << "X " << x << endl @@ -781,7 +775,7 @@ bool V2RhoT::Iterate() { T_n1 = T_n + c_Fdamp*(V - Vsyn)/dVdTsyn; deltaT = abs(T_n - T_n1); counter = counter + 1; - if(counter>10000) { + if (counter > 10000) { cout << "Too many iterations at point " << i << endl << "X(" << x << ") Y(" << y <<") Z(" << z << ") V(" << V << ")\n" << "Set T=-1\n"; @@ -791,7 +785,7 @@ bool V2RhoT::Iterate() { T_n = T_n1; } - if(verbose) { + if (verbose) { cout << "Iteration finished, deltaT = " << deltaT << ", T = " << T_n1 << endl; } @@ -802,10 +796,11 @@ bool V2RhoT::Iterate() { p5D.setZ(z); p5D.setProp(MantleRock->getRho()); data_T.append(p5D); - count_total.append(double(counter)); + count_total.append(static_cast(counter)); - progress = (int)round( (double(i)/double(n_V)) * 100.); - if(progress % 5 == 0) { + progress = static_cast( + round(static_cast(i)/static_cast(n_V))*100.0); + if (progress % 5 == 0) { printf("\rProgress: %i ", progress); } } @@ -813,7 +808,7 @@ bool V2RhoT::Iterate() { // Calculate average counts count_avrg = 0; - for(int i=0; i 0) { - VelTemp.readArgs(argc,argv); + if (argc > 0) { + VelTemp.readArgs(argc, argv); VelTemp.Info(); - VelTemp.readFile(VelTemp.FileIn(),"vox"); + VelTemp.readFile(VelTemp.FileIn(), "vox"); VelTemp.Iterate(); VelTemp.saveFile(VelTemp.FileOut()); } else { diff --git a/V2RhoT/V2RhoT.pro b/src/V2RhoT/V2RhoT.pro similarity index 76% rename from V2RhoT/V2RhoT.pro rename to src/V2RhoT/V2RhoT.pro index 629566f..696d564 100644 --- a/V2RhoT/V2RhoT.pro +++ b/src/V2RhoT/V2RhoT.pro @@ -16,24 +16,17 @@ # along with VeloDT. If not, see . # ################################################################################ TEMPLATE = app +DESTDIR = ../../bin TARGET = V2RhoT CONFIG -= app_bundle -DEPENDPATH += . -DEFINES += -INCLUDEPATH += . ../ ../common -LIBS += -# Input -HEADERS += V2RhoT.h \ - Rock.h \ - MineraldRhodT.h \ - ../common/ANSIICodes.h \ - ../common/PointClasses.h \ - ../common/PhysicalConstants.h \ - ../common/ERMs.h -SOURCES += V2RhoT.cpp \ - Rock.cpp \ - MineraldRhodT.cpp \ - ../common/PointClasses.cpp \ - ../common/ERMs.cpp -RESOURCES += +INCLUDEPATH += ../../include/common ../../include/V2RhoT + +LIBS += -L../common -lcommon + +SOURCES += V2RhoT.cpp Rock.cpp MineraldRhodT.cpp + +HEADERS += ../../include/V2RhoT/V2RhoT.h \ + ../../include/V2RhoT/Rock.h \ + ../../include/V2RhoT/MineraldRhodT.h + diff --git a/V2T/V2T.cpp b/src/V2T/V2T.cpp similarity index 83% rename from V2T/V2T.cpp rename to src/V2T/V2T.cpp index 27a174f..31dae70 100644 --- a/V2T/V2T.cpp +++ b/src/V2T/V2T.cpp @@ -68,7 +68,7 @@ void V2T::Info() { << "Output\n" << "------\n" << "Temperatures " << File_Out.toUtf8().data() << endl << endl; - if (use_t_crust || (PMethod=="simple")) { + if (use_t_crust || (PMethod == "simple")) { cout << "Densities\n" << "---------\n"; } @@ -76,7 +76,7 @@ void V2T::Info() { cout << "Crust " << rho_crust << "kg/m3\n" << "Mantle " << rho_mantle << "kg/m3\n"; } - if (PMethod=="simple") { + if (PMethod == "simple") { cout << "Average " << rho_avrg << " kg/m3\n\n"; } cout << "Other\n" @@ -184,7 +184,7 @@ void V2T::usage_extended() { bool V2T::test_data() { /// Check minima and maxima of data_Vs, t_crust and z_topo - if( use_t_crust && ( + if (use_t_crust && ( (x_min1 != x_min2) || (x_min2 != x_min3) || (x_max1 != x_max2) || (x_max2 != x_max3) || (y_min1 != y_min2) || (y_min2 != y_min3) || @@ -197,7 +197,7 @@ bool V2T::test_data() { } void V2T::argsError(QString val, bool ok) { - if(!ok) { + if (!ok) { cout << PRINT_ERROR "Invalid value in " << val.toUtf8().data() << endl; exit(1); } @@ -216,88 +216,84 @@ void V2T::readArgs(int &argc, char *argv[]) { okTopo = false; QStringList arg; - for(i = 0; i < argc; i++) { + for (i = 0; i < argc; i++) arg << argv[i]; - } - if( (argc>1) && (arg[1]=="-h" || arg[1]=="-help") ) { + if ((argc > 1) && (arg[1] == "-h" || arg[1] == "-help")) { usage(); - } - else if( (argc>1) && (arg[1]=="--help") ) { + } else if ((argc > 1) && (arg[1] == "--help")) { usage_extended(); - } - else if(argc<3) { + } else if (argc < 3) { cout << endl << endl << PRINT_ERROR "Not enough arguments!\n\n"; usage(); - } - else { - File_In=arg[1].toUtf8().data(); - File_Out=arg[2].toUtf8().data(); + } else { + File_In = arg[1].toUtf8().data(); + File_Out = arg[2].toUtf8().data(); - for(i = 3; i < argc; i++) { - if(arg[i]=="-ERM") { - ok=SetPMethod(arg[i+1]); + for (i=3; i < argc; i++) { + if (arg[i] == "-ERM") { + ok = SetPMethod(arg[i+1]); argsError(arg[i], ok); i++; - } else if(arg[i]=="--WRITE_P") { + } else if (arg[i] == "--WRITE_P") { WRITE_P(arg[i+1]); - } else if(arg[i]=="-outVs"){ + } else if (arg[i] == "-outVs") { outVs = true; - } else if(arg[i]=="-rc") { + } else if (arg[i] == "-rc") { rho_crust = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); + argsError(arg[i], ok); i++; - } else if(arg[i]=="-rm") { + } else if (arg[i] == "-rm") { rho_mantle = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); + argsError(arg[i], ok); i++; - } else if(arg[i]=="-ra") { + } else if (arg[i] == "-ra") { rho_avrg = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); + argsError(arg[i], ok); i++; - } else if(arg[i]=="-scaleZ") { + } else if (arg[i] == "-scaleZ") { scaleZ = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); + argsError(arg[i], ok); i++; - } else if(arg[i]=="-scaleVs") { + } else if (arg[i] == "-scaleVs") { scaleVs = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); + argsError(arg[i], ok); i++; - } else if(arg[i]=="-t_crust") { + } else if (arg[i] == "-t_crust") { File_t_crust = arg[i+1]; SetPMethod("crust"); okCrust = true; - } else if(arg[i]=="-z_topo") { + } else if (arg[i] == "-z_topo") { File_z_topo = arg[i+1]; SetPMethod("crust"); okTopo = true; - } else if(arg[i]=="-scatter") { + } else if (arg[i] == "-scatter") { ArbitraryPoints = true; - } else if(arg[i]=="-t") { + } else if (arg[i] == "-t") { threshold = arg[i+1].toDouble(&ok); - argsError(arg[i],ok); - } else if(arg[i]=="-v") { - verbose=true; + argsError(arg[i], ok); + } else if (arg[i] == "-v") { + verbose = true; } } } // Some logical checks - if(ArbitraryPoints && use_t_crust) { + if (ArbitraryPoints && use_t_crust) { cout << PRINT_ERROR "Defined crustal thickness/topography using scattered " "data. Exit.\n"; exit(1); } - if(okTopo ^ okCrust) { + if (okTopo ^ okCrust) { cout << PRINT_ERROR "Not enough arguments:\n"; - if(okTopo) { + if (okTopo) { cout << "No crustal thickness file defined!\n"; } else { cout << "No topographic elevation file defined!\n"; } exit(1); } - if( (use_t_crust == false) && (ERM->type() == "Undefined") ) { + if ((use_t_crust == false) && (ERM->type() == "Undefined")) { ERM->set("AK135"); } } @@ -308,7 +304,7 @@ bool V2T::readFile(QString InName, QString InType) { double vs_min = -1; bool okx, oky, okz, okval, okGrid; - okGrid = false; // Used to check if grid size was extracted from input voxel + okGrid = false; // Used to check if grid size was extracted from input voxel cout << "Reading file: " << InName.toUtf8().data() << endl; QFile file(InName); @@ -317,29 +313,28 @@ bool V2T::readFile(QString InName, QString InType) { x_max1 = x_max2 = x_max3 = y_max1 = y_max2 = y_max3 = z_max1 = -x_min1; int n = 0; - if(file.open(QIODevice::ReadOnly)) { + if (file.open(QIODevice::ReadOnly)) { QTextStream stream( &file ); - while(!stream.atEnd()) { + while (!stream.atEnd()) { n++; QString t = stream.readLine().simplified(); /// Topography or crustal thickness - - if(!t.startsWith("#") && (InType == "topo" || InType == "crust")) { + if (!t.startsWith("#") && (InType == "topo" || InType == "crust")) { QStringList vals = t.split(" "); - if(vals.count() != 3) { + if (vals.count() != 3) { file.close(); cout << PRINT_ERROR "In header of " << InName.toUtf8().data() << " - grid size."; exit(1); } else { - x=vals[0].toDouble(&okx); - y=vals[1].toDouble(&oky); - z=vals[2].toDouble(&okz); - if(okx && oky && okz) { - Point3D p3D(x,y,z); - if(InType == "crust") { + x = vals[0].toDouble(&okx); + y = vals[1].toDouble(&oky); + z = vals[2].toDouble(&okz); + if (okx && oky && okz) { + Point3D p3D(x, y, z); + if (InType == "crust") { t_crust.append(p3D); } else if (InType == "topo") { z_topo.append(p3D); @@ -350,50 +345,50 @@ bool V2T::readFile(QString InName, QString InType) { return false; } - //Determine min/max - if(InType == "crust") { - if(x < x_min2) { + // Determine min/max + if (InType == "crust") { + if (x < x_min2) { x_min2 = x; - } else if(x > x_max2) { + } else if (x > x_max2) { x_max2 = x; } - if(y < y_min2) { + if (y < y_min2) { y_min2 = y; - } else if(y > y_max2) { + } else if (y > y_max2) { y_max2 = y; } } else if (InType =="topo") { - if(x < x_min3) { + if (x < x_min3) { x_min3 = x; - } else if(x > x_max3) { + } else if (x > x_max3) { x_max3 = x; } - if(y < y_min3) { + if (y < y_min3) { y_min3 = y; - } else if(y > y_max3) { + } else if (y > y_max3) { y_max3 = y; } } } - } - /// S-wave velocity grid - else if(t.startsWith("# Grid_size:") && InType == "vox") { - t.remove(0,12); + } else if (t.startsWith("# Grid_size:") && InType == "vox") { + // S-wave velocity grid + t.remove(0, 12); QStringList gridSize = t.split("x", QString::SkipEmptyParts); nX = gridSize[0].toInt(&okx); nY = gridSize[1].toInt(&oky); nZ = gridSize[2].toInt(&okz); - if(!okx || !oky || !okz) { + if (!okx || !oky || !okz) { cout << PRINT_ERROR "In header of " << InName.toUtf8().data() << " - grid size." << endl; exit(1); } okGrid = true; - } else if(!t.startsWith("#") && InType == "vox") { + } else if (!t.startsWith("#") && InType == "vox") { QStringList vals = t.split(" "); - if(vals.count() != 4) { + if (vals.count() != 4) { file.close(); - cout << PRINT_ERROR "In line " << n << ": unkown amount of columns.\n"; + cout << PRINT_ERROR "In line " << n << ": unkown amount of columns." + << endl; exit(1); } else { x = vals[0].toDouble(&okx); @@ -401,8 +396,8 @@ bool V2T::readFile(QString InName, QString InType) { z = scaleZ*vals[2].toDouble(&okz); val = scaleVs*vals[3].toDouble(&okval); - if(okx && oky && okz && okval) { - Point4D p4D(x,y,z,val); + if (okx && oky && okz && okval) { + Point4D p4D(x, y, z, val); data_Vs.append(p4D); } else { file.close(); @@ -412,36 +407,38 @@ bool V2T::readFile(QString InName, QString InType) { // Get minima and maxima // X - if(x < x_min1) { + if (x < x_min1) { x_min1 = x; - } else if(x > x_max1) { + } else if (x > x_max1) { x_max1 = x; } // Y - if(y < y_min1) { + if (y < y_min1) { y_min1 = y; - } else if(y > y_max1) { + } else if (y > y_max1) { y_max1 = y; } // Z - if(z < z_min1) { + if (z < z_min1) { z_min1 = z; - } else if(z > z_max1) { + } else if (z > z_max1) { z_max1 = z; } // Get minimum vs to check for km/s or m/s - if(vs_min == -1){ vs_min = val; } - if(val < vs_min) { vs_min = val; } + if (vs_min == -1) + vs_min = val; + if (val < vs_min) + vs_min = val; } } } - if(!ArbitraryPoints && !okGrid) { + if (!ArbitraryPoints && !okGrid) { cout << PRINT_WARNING "Could not find grid information. Set output to " "scattered data.\n"; ArbitraryPoints = true; } file.close(); - if(vs_min > 10) { + if (vs_min > 10) { cout << PRINT_WARNING "Minimum Vs is " << vs_min << " which is unusually " "high. Vs must be in km/s. Use -scaleVs to correct.\n"; } @@ -465,14 +462,14 @@ bool V2T::saveFile(QString OutName) { T_info += QString("# z-factor: %1\n").arg(scaleZ, 0, 'f'); T_info += QString("# Vs-factor: %1\n").arg(scaleVs, 0, 'f'); T_info += QString("# Pressure calculation method: %1\n").arg(PMethod); - if(PMethod == "simple") { + if (PMethod == "simple") { T_info += QString("# Average density: %1 kg/m3\n").arg(rho_avrg, 0, 'f'); } else if (PMethod == "crust") { T_info += QString("# Density crust: %1 kg/m3\n").arg(rho_crust, 0, 'f'); T_info += QString("# Density mantle: %1 kg/m3").arg(rho_mantle, 0, 'f'); } // GMS compatible file - if(!ArbitraryPoints) { + if (!ArbitraryPoints) { T_header = QString("# Type: GMS GridPoints\n" "# Version: 2\n" "# Description:\n" @@ -481,19 +478,19 @@ bool V2T::saveFile(QString OutName) { "# Field: 2 Y / m\n" "# Field: 3 Z / m\n" "# Field: 4 T / degC\n"); - if(outVs) { + if (outVs) { T_header += QString("# Field: 4 VsObs / degC\n" "# Field: 5 VsCalc / degC\n"); } T_header += QString("# Projection: Local Rectangular\n"); T_header += QString("# Information from grid:\n"); T_header += QString("# Grid_size: %1 x %2 x %3\n").arg(nX).arg(nY).arg(nZ); - T_header += QString("# Grid_X_range: %1 to %2\n").arg(x_min1,0,'f') - .arg(x_max1,0,'f'); - T_header += QString("# Grid_Y_range: %1 to %2\n").arg(y_min1,0,'f') - .arg(y_max1,0,'f'); - T_header += QString("# Grid_Z_range: %1 to %2\n").arg(z_min1,0,'f') - .arg(z_max1,0,'f'); + T_header += QString("# Grid_X_range: %1 to %2\n").arg(x_min1, 0, 'f') + .arg(x_max1, 0, 'f'); + T_header += QString("# Grid_Y_range: %1 to %2\n").arg(y_min1, 0, 'f') + .arg(y_max1, 0, 'f'); + T_header += QString("# Grid_Z_range: %1 to %2\n").arg(z_min1, 0, 'f') + .arg(z_max1, 0, 'f'); T_header += T_info; T_header += QString("# End:"); } else { @@ -503,7 +500,7 @@ bool V2T::saveFile(QString OutName) { "# 2 - Y\n" "# 3 - Z / m\n" "# 4 - T / degC\n"); - if(outVs) { + if (outVs) { T_header += QString("# 5 - VsObs / degC\n" "# 6 - VsCalc / degC\n"); } @@ -513,8 +510,9 @@ bool V2T::saveFile(QString OutName) { cout << "Writing temperature file " << OutName.toUtf8().data() << endl; QFile tmp(OutName); - if(!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { - cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() << endl; + if (!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { + cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() + << endl; exit(1); } @@ -523,13 +521,13 @@ bool V2T::saveFile(QString OutName) { out.setFieldAlignment(QTextStream::AlignRight); out.setRealNumberNotation(QTextStream::FixedNotation); out << T_header.toUtf8().data(); - for(int i=0; ipressure(z); } else { return pressure_simple(z); @@ -608,15 +606,17 @@ double V2T::pressure_crust(double x, double y, double z) { double P_crust, t_mantle, P_mantle; // Get index in data_t_crust and data_z_topo - i=0; - while( (t_crust[i].x()!=x) && (t_crust[i].y()!=y)){i++;} + i = 0; + while ((t_crust[i].x() != x) && (t_crust[i].y() != y)) + i++; i_crust = i; i = 0; - while( (z_topo[i].x()!=x) && (z_topo[i].y()!=y)){i++;} + while ((z_topo[i].x() != x) && (z_topo[i].y() != y)) + i++; i_topo=i; - P_crust = rho_crust*c_g*t_crust[i_crust].z(); //P from crust - t_mantle = z_topo[i_topo].z() - t_crust[i_crust].z() - z; //Thickness mantle + P_crust = rho_crust*c_g*t_crust[i_crust].z(); // P from crust + t_mantle = z_topo[i_topo].z() - t_crust[i_crust].z() - z; // Thickness mantle if (t_mantle < 0) { cout << "Error: Mantle thickness < 0\nExit\n"; @@ -668,12 +668,12 @@ bool V2T::newton() { count_fail = 0; // Counts how often solution could not be found j = 0; - for(int i=0; i threshold) { - numerator = ftheta(VsS,P,theta_i1); - denominator = dfdtheta(P,theta_i1); + while (delta_theta > threshold) { + numerator = ftheta(VsS, P, theta_i1); + denominator = dfdtheta(P, theta_i1); - if(denominator == 0) { + if (denominator == 0) { // If denominator = 0, add 1 to counter and continue with next value count_zero += 1; theta_i1 = -1; @@ -707,7 +708,7 @@ bool V2T::newton() { j++; - if(j > 10000) { + if (j > 10000) { printf("\r" PRINT_WARNING "Too many iterations!\n"); printf("Coordinates: %f, %f, %f, Vs: %f\n",x,y,z,Vs); count_fail += 1; @@ -715,7 +716,7 @@ bool V2T::newton() { break; } } - if(verbose){ + if (verbose){ cout << "Estimated temperature " << theta_i1 << endl << "delta_theta " << delta_theta << endl << "Iteration steps " << j << endl << endl; @@ -724,7 +725,7 @@ bool V2T::newton() { theta_i1 = (VsS-c_c)/c_m; } - if(verbose) { + if (verbose) { cout << endl << "Point #" << i << endl << "Depth z / m " << z << endl @@ -734,8 +735,9 @@ bool V2T::newton() { << "Calculated temperature theta_i / degC " << theta_i1 << endl << endl; } else { - progress = (int)round( (double(i)/double(n)) * 100); - if(progress % 5 == 0) { + progress = static_cast( + round(static_cast(i)/static_cast(n))*100.0); + if (progress % 5 == 0) { printf("\rProgress: %i ", progress); } } @@ -748,7 +750,7 @@ bool V2T::newton() { data_T.append(p4D); // Calculate synthetic velocity from temperature - if(outVs) { + if (outVs) { double VsSCalc = c_m*theta_i1 + c_c + c_A*exp(-(c_E*1000 + P*c_Va)/c_R/(theta_i1 + c_K)); double VsCalc = VsSCalc*(1 + c_bV*(fabs(z)/1000 - 50)); @@ -776,15 +778,16 @@ void V2T::WRITE_P(QString method) { zi = zmin; while (zi <= zmax) { TEST_Z.append(zi); - TEST_P.append(pressure(0,0,zi)); + TEST_P.append(pressure(0, 0, zi)); zi+=dz; } cout << "Writing pressures to " << OutName.toUtf8().data() << endl; QFile tmp(OutName); - if(!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { - cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() << endl; + if (!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) { + cout << PRINT_ERROR "Could not open file " << OutName.toUtf8().data() + << endl; exit(1); } @@ -796,7 +799,7 @@ void V2T::WRITE_P(QString method) { out.setFieldAlignment(QTextStream::AlignRight); out.setRealNumberNotation(QTextStream::FixedNotation); out << T_header.toUtf8().data() << endl; - for(i = 0; i < TEST_P.length(); i++) { + for (i=0; i < TEST_P.length(); i++) { out << TEST_Z[i]; out << "\t"; out << TEST_P[i]; @@ -812,13 +815,13 @@ void V2T::WRITE_P(QString method) { //############################################################################## int main(int argc, char *argv[]) { V2T VelTemp; - if(argc > 0) { - VelTemp.readArgs(argc,argv); + if (argc > 0) { + VelTemp.readArgs(argc, argv); VelTemp.Info(); - VelTemp.readFile(VelTemp.FileIn(),"vox"); - if(VelTemp.UseCrust()) { - VelTemp.readFile(VelTemp.FileTCrust(),"crust"); - VelTemp.readFile(VelTemp.FileZTopo(),"topo"); + VelTemp.readFile(VelTemp.FileIn(), "vox"); + if (VelTemp.UseCrust()) { + VelTemp.readFile(VelTemp.FileTCrust(), "crust"); + VelTemp.readFile(VelTemp.FileZTopo(), "topo"); } VelTemp.test_data(); VelTemp.newton(); diff --git a/V2T/V2T.pro b/src/V2T/V2T.pro similarity index 84% rename from V2T/V2T.pro rename to src/V2T/V2T.pro index 7b54c4d..d6a76df 100644 --- a/V2T/V2T.pro +++ b/src/V2T/V2T.pro @@ -16,17 +16,15 @@ # along with VeloDT. If not, see . # ################################################################################ TEMPLATE = app +DESTDIR = ../../bin TARGET = V2T CONFIG -= app_bundle -DESTDIR = ./ -INCLUDEPATH += ../ ../common +INCLUDEPATH += ../../include/common ../../include/V2T -HEADERS += V2T.h \ - ../common/PointClasses.h \ - ../common/ANSIICodes.h \ - ../common/ERMs.h +LIBS += -L../common -lcommon + +SOURCES += V2T.cpp + +HEADERS += V2T.h -SOURCES += V2T.cpp \ - ../common/PointClasses.cpp \ - ../common/ERMs.cpp diff --git a/common/ERMs.cpp b/src/common/ERMs.cpp similarity index 93% rename from common/ERMs.cpp rename to src/common/ERMs.cpp index 03b8da7..fa8e2cc 100644 --- a/common/ERMs.cpp +++ b/src/common/ERMs.cpp @@ -17,7 +17,8 @@ *******************************************************************************/ #include "ERMs.h" -using namespace std; +using std::cout; +using std::endl; const int AK135_LEN = 27; const int PREM_LEN = 107; @@ -59,7 +60,7 @@ EarthReferenceModel::EarthReferenceModel(QString type) { bool EarthReferenceModel::INIT_AK135() { // Initialise AK135 from array - for(int i = 0; i < AK135_LEN; i++){ + for (int i = 0; i < AK135_LEN; i++) { ERMz.append(ARR_AK135z[i]); ERMrho.append(ARR_AK135rho[i]); } @@ -69,7 +70,7 @@ bool EarthReferenceModel::INIT_AK135() { bool EarthReferenceModel::INIT_PREM() { // Inititalise PREM from array - for(int i = 0; i < PREM_LEN; i++){ + for (int i = 0; i < PREM_LEN; i++) { ERMz.append(ARR_PREMz[i]); ERMrho.append(ARR_PREMrho[i]); } @@ -88,7 +89,7 @@ double EarthReferenceModel::pressure(double z) { bool DepthCalculated = false; Pcalc = 0.; z_abs = fabs(z); - while ( !DepthCalculated ) { + while (!DepthCalculated) { z1 = ERMz[i]*1000; z2 = ERMz[i+1]*1000; Rho1 = ERMrho[i]*1000; @@ -108,9 +109,9 @@ double EarthReferenceModel::pressure(double z) { bool EarthReferenceModel::set(QString type) { // Define a reference model - if(type=="AK135") { + if (type == "AK135") { INIT_AK135(); - } else if(type=="PREM") { + } else if (type == "PREM") { INIT_PREM(); } else { std::cout << "Unknown reference model " << type.toUtf8().data() << endl; @@ -130,13 +131,13 @@ bool EarthReferenceModel::writeP(double dz) { double nzfloat = (zmax - zmin)/dz + 1; cout << "Calculating pressures" << endl << "ERM type is " << ERMtype.toUtf8().data() << endl; - if(fmod(nzfloat,1) == 0.) { + if (fmod(nzfloat, 1) == 0.) { // Calculate pressures - int nz = (int) nzfloat; + int nz = static_cast(nzfloat); QVector depths; QVector pressures; float z; - for(int i=0; i. * *******************************************************************************/ #include "PointClasses.h" + /******************************************************************************* Point2D *******************************************************************************/ @@ -49,13 +50,13 @@ Point3D::Point3D(double x, double y, double z) { } Point3D &Point3D::operator=(const Point3D p) { - Vals[0]=p.Vals[0]; - Vals[1]=p.Vals[1]; - Vals[2]=p.Vals[2]; + Vals[0] = p.Vals[0]; + Vals[1] = p.Vals[1]; + Vals[2] = p.Vals[2]; return *this; } -double &Point3D::operator [](int idx) { +double &Point3D::operator[](int idx) { return Vals[idx]; } @@ -67,24 +68,24 @@ Point4D::Point4D() { Vals.append(0.0); Vals.append(0.0); Vals.append(0.0); -}; +} Point4D::Point4D(double x, double y, double z, double v) { Vals.append(x); // at least 4 values, representing x,y,z,v Vals.append(y); Vals.append(z); Vals.append(v); -}; +} -double &Point4D::operator [](int idx) { +double &Point4D::operator[](int idx) { return Vals[idx]; } Point4D &Point4D::operator=(const Point4D p) { - Vals[0]=p.Vals[0]; - Vals[1]=p.Vals[1]; - Vals[2]=p.Vals[2]; - Vals[3]=p.Vals[3]; + Vals[0] = p.Vals[0]; + Vals[1] = p.Vals[1]; + Vals[2] = p.Vals[2]; + Vals[3] = p.Vals[3]; return *this; } @@ -107,14 +108,14 @@ Point5D::Point5D(double x, double y, double z, double v, double prop) { Vals.append(prop); } -double &Point5D::operator [](int idx) { +double &Point5D::operator[](int idx) { return Vals[idx]; } Point5D &Point5D::operator=(const Point5D p) { - Vals[0]=p.Vals[0]; - Vals[1]=p.Vals[1]; - Vals[2]=p.Vals[2]; - Vals[3]=p.Vals[3]; + Vals[0] = p.Vals[0]; + Vals[1] = p.Vals[1]; + Vals[2] = p.Vals[2]; + Vals[3] = p.Vals[3]; return *this; } diff --git a/src/common/common.pro b/src/common/common.pro new file mode 100644 index 0000000..435dbc3 --- /dev/null +++ b/src/common/common.pro @@ -0,0 +1,9 @@ +INCLUDEPATH += ../../include/common +WARNINGS += -Wall +TEMPLATE = lib +CONFIG += staticlib +SOURCES += ERMS.cpp PointClasses.cpp +HEADERS += ../../include/common/ERMS.h \ + ../../include/common/PointClasses.h \ + ../../include/common/ANSIICodes.h \ + ../../include/common/PhysicalConstants.h