diff --git a/.travis.yml b/.travis.yml
index 6bd35d3..d96fbda 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -9,13 +9,13 @@ script:
- qmake
- make
after_success:
-- tar -zcvf bin/bin-VeloDT-prerelease-ubuntu-qt${QT_BASE}-x86_64.tar.gz bin/V2RhoT
+- tar -zcvf bin/bin-VeloDT-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
+ file: bin/bin-VeloDT-ubuntu-qt${QT_BASE}-x86_64.tar.gz
on:
tags: true
repo: cmeessen/VeloDT
diff --git a/CHANGELOG b/CHANGELOG
index e00442b..49918c8 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -3,7 +3,19 @@ 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
+## [v1.2.0] - 2020-06-16
+
+### Added
+
+- T2Rho applies the equations of
+ [Goes et al. (2000)](https://doi.org/10.1029/1999JB900300) to compute
+ densities from a given mineral composition, temperature and pressure
+
+### Fixed
+
+- File name of binaries is not "prerelease" anymore
+
+## [v1.1.0] - 2020-06-16
### Added
diff --git a/Example/Makefile b/Example/Makefile
index be61071..735c928 100644
--- a/Example/Makefile
+++ b/Example/Makefile
@@ -1,5 +1,6 @@
V2RhoT = ../bin/V2RhoT
V2T = ../bin/V2T
+T2Rho = ../bin/T2Rho
# Plot properties
FIGURE = results.ps
ROI_T = 450/1450/-200/-50
@@ -31,13 +32,16 @@ plot: convert
# GGV Rho Vs
gmt psxy V2RhoT_Vs.dat -i5,2+s0.001 -J -R$(ROI_R) -W1p,red -X7 -BwSne+t'Density' -Bx+l'@~r@~ / kg/m@+3@+' -Baf -O -P -K >> $(FIGURE)
# GGV Rho Vp
- gmt psxy V2RhoT_Vp.dat -i5,2+s0.001 -J -R -W1p,blue -O -P >> $(FIGURE)
+ gmt psxy V2RhoT_Vp.dat -i5,2+s0.001 -J -R -W1p,blue -O -P -K >> $(FIGURE)
+ # Rho from PMK
+ gmt psxy T2Rho_V2T_Vs.dat -i5,2+s0.001 -J -R -W1p,red,. -O -P >> $(FIGURE)
gmt psconvert -A -Tg -Z $(FIGURE) -E100 -Qt -Qg
convert: createfiles Vs.dat Vp.dat
$(V2RhoT) Vs.dat V2RhoT_Vs.dat -type S -ERM PREM -scaleZ -1000 -scaleV 1000 -compp 0 -xfe 0.086
$(V2RhoT) Vp.dat V2RhoT_Vp.dat -type P -ERM PREM -scaleZ -1000 -scaleV 1000 -compp 0 -xfe 0.086
$(V2T) Vs.dat V2T_Vs.dat -ERM PREM -scaleZ -1000
+ $(V2RhoT) V2T_Vs.dat T2Rho_V2T_Vs.dat -ERM PREM
createfiles:
awk '{if(NR>3){print($$1,$$2,$$3,$$4)}}' PREM.dat > Vp.dat
diff --git a/Example/results.png b/Example/results.png
index b1acad5..b570c23 100644
Binary files a/Example/results.png and b/Example/results.png differ
diff --git a/README.md b/README.md
index 43844e2..507d953 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
# VeloDT
[![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)
+[![Build Status](https://travis-ci.org/cmeessen/VeloDT.svg?branch=master)](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.
@@ -23,9 +23,9 @@ or download the [latest release](https://github.com/cmeessen/VeloDT/releases/lat
- [**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.md) implements the method by [Priestley and McKenzie
(2006)](https://doi.org/10.1016/j.epsl.2006.01.008) and converts *Vs* to temperature
+- [**T2Rho**](./T2Rho.md) applies the equations of [Goes et al. (2000)](https://doi.org/10.1029/1999JB900300) to compute densities from a given mineral composition, temperature and pressure
Please refer to the original publications for validity and pitfalls of the methods.
diff --git a/T2Rho.md b/T2Rho.md
new file mode 100644
index 0000000..760b731
--- /dev/null
+++ b/T2Rho.md
@@ -0,0 +1,58 @@
+# T2Rho
+
+This tool helps to directly compute densities from the output of [`V2T`](./V2T.md), using the equations given in [Goes et al. (2000)](https://doi.org/10.1029/1999JB900300).
+
+## Compilation
+
+Before compilation make sure to have Qt and qmake installed. Then go to this folder and run
+
+```
+qmake
+make
+```
+
+This generates the `T2Rho` executable.
+
+## Using T2Rho
+
+Executing `T2Rho -h` displays the available console commands together with two examples.
+
+```
+usage: T2Rho File_In File_Out [options]
+
+T2Rho is designed to directly process output files from V2T in
+order to obtain densities from temperatures assuming a specific
+mantle composition.
+
+ Required input parameters:
+ --------------------------
+ File_In Path and name of grid file containing x y z Vs
+ File_Out Output file name and path
+
+ Option Value Default Description
+ ------ ----- ------- -----------
+ -h This information
+ -ERM string AK135 P calculation method AK135 or PREM
+ -compc vals Custom rock composition
+ -compc Ol Opx Cpx Sp Gnt
+ -compp val 0 Use predefined rock compositions:
+ 0 - Garnet Lherzolite after (Jordan,
+ 1979; Goes et al, 2000)
+ 1 - On-cratonic (Shapiro and Ritzwoller, 2004)
+ 2 - Off-cratonic (Shapiro and Ritzwoller, 2004)
+ 3 - Oceanic (Shapiro and Ritzwoller, 2004)
+ -xfe val 0.1 Define iron content of the rock in mole fraction
+```
+
+### Mandatory arguments
+
+`T2Rho` requires the input file `File_In`, containing x y z and vs, and the name of the output file `File_Out`. Input units for z is masl, for vs km/s.
+
+### Pressure calculation
+
+Standard calculation of **pressure** uses the earth reference model AK135.
+- options are `AK135`, `PREM` or `simple`
+- `-ERM PREM` activates pressure calculation with PREM
+- `-ERM simple` uses the average density defined with `-ra`
+- an experimental feature is the pressure calculation using topography and crustal thickness. This is activated by using `-t_crust FILENAME` and `-z_topo FILENAME`, which both require EarthVision formatted grids containing crustal thickness and topographic elevation. The pressure is then calculated assuming constant density for the crust (`-rc 2890`) and mantle (`-rm 3300`)
+- `-ra` defines an average density which is then used to calculate the pressure
diff --git a/VeloDT.pro b/VeloDT.pro
index dc724bd..d63a62d 100644
--- a/VeloDT.pro
+++ b/VeloDT.pro
@@ -1,6 +1,6 @@
TEMPLATE = subdirs
CONFIG -= app_bundle
CONFIG += ordered
-SUBDIRS += src/common src/V2RhoT src/V2T
+SUBDIRS += src/common src/V2RhoT src/V2T src/T2Rho
V2RhoT.depends = common
V2T.depends = common
diff --git a/include/T2Rho/T2Rho.h b/include/T2Rho/T2Rho.h
new file mode 100644
index 0000000..18c657e
--- /dev/null
+++ b/include/T2Rho/T2Rho.h
@@ -0,0 +1,66 @@
+/*******************************************************************************
+* Copyright (C) 2018 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 . *
+*******************************************************************************/
+#ifndef T2RHO_H_
+#define T2RHO_H_
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "PhysicalConstants.h"
+#include "PointClasses.h"
+#include "ERMs.h"
+
+class T2Rho {
+ QString PMethod;
+ QString file_in;
+ QString file_out;
+ EarthReferenceModel * ERM;
+ // Mineral properties
+ QList m_rho;
+ QList m_drhodX;
+ QList m_K;
+ QList m_dKdT;
+ QList m_dKdP;
+ QList m_dKdPdX;
+ QList m_dKdX;
+ QList m_alpha0;
+ QList data_out;
+ // Rock properties
+ double r_XFe;
+ QList r_comp;
+ // Functions
+ bool SetPMethod(QString method);
+ void argsError(QString val, bool ok);
+ void setComp(QList composition);
+ void setComp(int c);
+
+ public:
+ T2Rho();
+ void usage();
+ void info();
+ void readArgs(int &argc, char *argv[]);
+ void readFile();
+ void writeFile();
+};
+
+#endif // T2RHO_H_
diff --git a/include/common/PointClasses.h b/include/common/PointClasses.h
index bed7e58..de20f7f 100644
--- a/include/common/PointClasses.h
+++ b/include/common/PointClasses.h
@@ -20,81 +20,85 @@
#define POINTCLASSES_H_
class Point2D {
+private:
QList Vals;
- public:
+public:
Point2D();
Point2D(double x, double y);
// Assigning properties
- void setX(double x) {Vals[0] = x;}
- void setY(double y) {Vals[1] = y;}
+ 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];}
+ 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);
+ 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 {
- 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);
+ 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 {
- 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];}
+ 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];}
+ double &operator[](int idx);
+ Point5D &operator=(Point5D p);
+ double *p(int i){return &Vals[i];}
};
#endif // POINTCLASSES_H_
diff --git a/src/T2Rho/T2Rho.cpp b/src/T2Rho/T2Rho.cpp
new file mode 100644
index 0000000..cfc0ac0
--- /dev/null
+++ b/src/T2Rho/T2Rho.cpp
@@ -0,0 +1,367 @@
+/*******************************************************************************
+* Copyright (C) 2018 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 "T2Rho.h"
+
+using std::cout;
+using std::endl;
+
+T2Rho::T2Rho() {
+ PMethod = "AK135";
+ ERM = new EarthReferenceModel(PMethod);
+
+ // Mineral properties Tab. A1 in Goes et al (2000)
+ m_rho.append(3222);
+ m_rho.append(3198);
+ m_rho.append(3280);
+ m_rho.append(3578);
+ m_rho.append(3565);
+
+ m_drhodX.append(1182);
+ m_drhodX.append(804);
+ m_drhodX.append(377);
+ m_drhodX.append(702);
+ m_drhodX.append(758);
+
+ m_K.append(1.29E+11);
+ m_K.append(1.11E+11);
+ m_K.append(1.05E+11);
+ m_K.append(1.98E+11);
+ m_K.append(1.73E+11);
+
+ m_dKdT.append(-1.60E+07);
+ m_dKdT.append(-1.20E+07);
+ m_dKdT.append(-1.30E+07);
+ m_dKdT.append(-2.80E+07);
+ m_dKdT.append(-2.10E+07);
+
+ m_dKdP.append(4.2);
+ m_dKdP.append(6);
+ m_dKdP.append(6.2);
+ m_dKdP.append(5.7);
+ m_dKdP.append(4.9);
+
+ m_dKdPdX.append(0);
+ m_dKdPdX.append(0);
+ m_dKdPdX.append(-1.9);
+ m_dKdPdX.append(0);
+ m_dKdPdX.append(0);
+
+ m_dKdX.append(0);
+ m_dKdX.append(-10E9);
+ m_dKdX.append(13E9);
+ m_dKdX.append(12E9);
+ m_dKdX.append(7E9);
+
+ m_alpha0.append(2.01E-5);
+ m_alpha0.append(3.871E-5);
+ m_alpha0.append(3.206E-5);
+ m_alpha0.append(6.969E-5);
+ m_alpha0.append(9.91E-6);
+
+ // Initiate rock properties
+ r_XFe = 0.1;
+ r_comp.append(0);
+ r_comp.append(0);
+ r_comp.append(0);
+ r_comp.append(0);
+ r_comp.append(0);
+ setComp(0);
+}
+
+void T2Rho::usage() {
+ cout << endl
+ << "usage: T2Rho File_In File_Out [options]" << endl
+ << endl
+ << "T2Rho is designed to directly process output files from V2T in" << endl
+ << "order to obtain densities from temperatures assuming a specific" << endl
+ << "mantle composition." << endl
+ << endl
+ << " Required input parameters:" << endl
+ << " --------------------------" << endl
+ << " File_In Path and name of grid file containing x y z Vs" << endl
+ << " File_Out Output file name and path" << endl
+ << endl
+ << " Option Value Default Description" << endl
+ << " ------ ----- ------- -----------" << endl
+ << " -h This information" << endl
+ << " -ERM string AK135 P calculation method AK135 or PREM" << endl
+ << " -compc vals Custom rock composition" << endl
+ << " -compc Ol Opx Cpx Sp Gnt" << endl
+ << " -compp val 0 Use predefined rock compositions:" << endl
+ << " 0 - Garnet Lherzolite after (Jordan," << endl
+ << " 1979; Goes et al, 2000)" << endl
+ << " 1 - On-cratonic (Shapiro and Ritzwoller, 2004)" << endl
+ << " 2 - Off-cratonic (Shapiro and Ritzwoller, 2004)" << endl
+ << " 3 - Oceanic (Shapiro and Ritzwoller, 2004)" << endl
+ << " -xfe val 0.1 Define iron content of the rock in mole fraction" << endl
+ << endl;
+ exit(0);
+}
+
+bool T2Rho::SetPMethod(QString method) {
+ if (ERM->set(method)) {
+ PMethod = method;
+ return true;
+ } else {
+ return false;
+ }
+}
+
+void T2Rho::setComp(QList composition) {
+ // Check length
+ if (!(composition.length() == 5)) {
+ cout << PRINT_ERROR "Number of mineral phases must be equal to 5." << endl;
+ }
+ // Check sum
+ double compsum = 0;
+ for (int i=0; i < 5; i++) {
+ compsum += composition[i];
+ }
+ if (!(compsum == 1)) {
+ cout << PRINT_ERROR "Sum of the phases is unequal to 1." << endl;
+ exit(1);
+ } else {
+ for (int i=0; i < 5; i++) {
+ r_comp[i] = composition[i];
+ }
+ }
+}
+
+void T2Rho::setComp(int c) {
+ switch (c) {
+ case 0:
+ // Default composition after Goes et al (2000)
+ r_comp[0] = 0.67;
+ r_comp[1] = 0.225;
+ r_comp[2] = 0.045;
+ r_comp[3] = 0;
+ r_comp[4] = 0.06;
+ break;
+ case 1:
+ // On-cratonic Shapiro and Ritzwoller (2004)
+ r_comp[0] = 0.83;
+ r_comp[1] = 0.15;
+ r_comp[2] = 0;
+ r_comp[3] = 0;
+ r_comp[4] = 0.02;
+ break;
+ case 2:
+ // Off-cratonic Shapiro and Ritzwoller (2004)
+ r_comp[0] = 0.68;
+ r_comp[1] = 0.18;
+ r_comp[2] = 0.11;
+ r_comp[3] = 0;
+ r_comp[4] = 0.03;
+ break;
+ case 3:
+ // Oceanic Shapiro and Ritzwoller (2004)
+ r_comp[0] = 0.75;
+ r_comp[1] = 0.21;
+ r_comp[2] = 0.035;
+ r_comp[3] = 0.005;
+ r_comp[4] = 0.0;
+ break;
+ default:
+ cout << PRINT_ERROR "Unknown composition " << c << "\n. Exit.\n";
+ exit(1);
+ }
+}
+
+void T2Rho::argsError(QString val, bool ok) {
+ if (!ok) {
+ cout << PRINT_ERROR "Invalid value in " << val.toUtf8().data() << endl;
+ exit(1);
+ }
+}
+
+void T2Rho::readArgs(int &argc, char *argv[]) {
+ int i;
+ bool ok;
+ QStringList arg;
+
+ for (i=0; i < argc; i++)
+ arg << argv[i];
+
+ if ( (argc > 1) && (arg[1] == "-h" || arg[1] == "-help") ) {
+ usage();
+ } 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();
+
+ for (i=3; i < argc; i++) {
+ if (arg[i] == "-ERM") {
+ ok = SetPMethod(arg[i+1]);
+ argsError(arg[i], ok);
+ i++;
+ } else if (arg[i] == "-h") {
+ usage();
+ } else if (arg[i] == "-xfe") {
+ r_XFe = arg[i+1].toDouble(&ok);
+ argsError(arg[i], ok);
+ } else if (arg[i] == "-compc") {
+ QList compc;
+ for (int j=1; j < 6; j++) {
+ compc.append(arg[i+j].toDouble(&ok));
+ argsError(arg[i+j], ok);
+ }
+ setComp(compc);
+ } else if (arg[i] == "-compp") {
+ int compp = arg[i+1].toInt(&ok);
+ argsError(arg[i], ok);
+ setComp(compp);
+ }
+ }
+ }
+}
+
+void T2Rho::info() {
+ cout << "***********************************************************" << endl
+ << endl
+ << " Temperature to Density Conversion" << endl
+ << endl
+ << "***********************************************************" << endl
+ << endl
+ << "Settings" << endl
+ << "--------" << endl
+ << "Rock composition (ol,opx,cpx,gnt,sp):" << endl;
+ for (int i=0; i < 5; i++)
+ cout << r_comp[i] << " ";
+ cout << endl
+ << "Iron content XFe: " << r_XFe << endl
+ << "Earth reference model: " << ERM->type().toUtf8().data() << endl
+ << endl;
+}
+
+void T2Rho::readFile() {
+ /**
+ * Basically everything happens in this function. This function reads the
+ * input file and directly computes the density. The result is written into
+ * the list data_out.
+ **/
+ double x, y, z, T, P;
+ bool okx, oky, okz, okT;
+ double T0 = 293.5; // Reference temperature
+ double P0 = 0; // Reference pressure
+ QRegExp separator("(\\s)");
+
+ cout << "Reading file: " << file_in.toUtf8().data() << endl;
+ QFile file(file_in);
+
+ if (file.open(QIODevice::ReadOnly)) {
+ QTextStream stream(&file);
+
+ while (!stream.atEnd()) {
+ QString t = stream.readLine().simplified();
+
+ if (!t.startsWith("#")) {
+ QStringList vals = t.split(separator);
+ if (vals.count() != 4) {
+ file.close();
+ cout << PRINT_ERROR "In header of " << file_in.toUtf8().data()
+ << ". Counting " << vals.count() << " entries instead of 4."
+ << endl;
+ exit(1);
+ } else {
+ x = vals[0].toDouble(&okx);
+ y = vals[1].toDouble(&oky);
+ z = vals[2].toDouble(&okz);
+ T = vals[3].toDouble(&okT);
+ // Directly compute density
+ P = ERM->pressure(z+1);
+ // Compute phase properties and density
+ double rho_avrg = 0;
+ for (int i=0; i < 5; i++) {
+ double dT = T - T0;
+ double dP = P - P0;
+ double m_K_PX = m_dKdP[i] + r_XFe*m_dKdPdX[i];
+ double K_T = m_K[i] + dT*m_dKdT[i] + dP*m_K_PX;
+ double rho_PTX = m_rho[i]*(1 - m_alpha0[i]*dT + dP/K_T)
+ + m_drhodX[i]*r_XFe;
+ rho_avrg += r_comp[i]*rho_PTX;
+ }
+ // Add result to output data list
+ data_out.append(Point5D(x, y, z, T, rho_avrg));
+ }
+ }
+ }
+ }
+}
+
+void T2Rho::writeFile() {
+ QString header;
+ // Get timestamp
+ QDateTime CurrentTime = QDateTime::currentDateTime();
+ QString timefmt = "yyyy-MM-dd hh:mm:ss";
+ QString timestamp = CurrentTime.toString(timefmt);
+ // General information string
+ header = QString("# Created: %1\n").arg(timestamp);
+ header += QString("# Input: %1\n").arg(file_in);
+ header += QString("# Pressure calculation method: %1\n").arg(PMethod);
+ header += QString("# Mantle composition:\n");
+ header += QString("# Ol - %1\n").arg(r_comp[0], 5, 'f', 2);
+ header += QString("# Opx - %1\n").arg(r_comp[1], 5, 'f', 2);
+ header += QString("# Cpx - %1\n").arg(r_comp[2], 5, 'f', 2);
+ header += QString("# Sp - %1\n").arg(r_comp[3], 5, 'f', 2);
+ header += QString("# Gnt - %1\n").arg(r_comp[4], 5, 'f', 2);
+ header += QString("# Iron content XFe: %1\n").arg(r_XFe, 3, 'f', 2);
+ header += QString("# Point data\n"
+ "# Columns:\n"
+ "# 1 - X\n"
+ "# 2 - Y\n"
+ "# 3 - Z / m\n"
+ "# 4 - T / degC\n"
+ "# 5 - Rho / kg/m3\n");
+
+ cout << "Writing output file " << file_out.toUtf8().data() << endl;
+
+ QFile tmp(file_out);
+ if (!tmp.open(QIODevice::WriteOnly | QIODevice::Text)) {
+ cout << PRINT_ERROR "Could not open file " << file_out.toUtf8().data()
+ << endl;
+ exit(1);
+ }
+
+ QTextStream out(&tmp);
+ out.setRealNumberPrecision(2);
+ out.setFieldAlignment(QTextStream::AlignRight);
+ out.setRealNumberNotation(QTextStream::FixedNotation);
+ out << header.toUtf8().data();
+ for (int i=0; i < data_out.length(); i++) {
+ out << data_out[i][0] << "\t"
+ << data_out[i][1] << "\t"
+ << data_out[i][2] << "\t"
+ << data_out[i][3] << "\t"
+ << data_out[i][4] << "\n";
+ }
+ tmp.close();
+}
+
+int main(int argc, char *argv[]) {
+ T2Rho converter;
+ if (argc > 0) {
+ converter.readArgs(argc, argv);
+ converter.info();
+ converter.readFile();
+ converter.writeFile();
+ } else {
+ converter.usage();
+ }
+}
\ No newline at end of file
diff --git a/src/T2Rho/T2Rho.pro b/src/T2Rho/T2Rho.pro
new file mode 100644
index 0000000..7b160d8
--- /dev/null
+++ b/src/T2Rho/T2Rho.pro
@@ -0,0 +1,29 @@
+################################################################################
+# Copyright (C) 2018 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 . #
+################################################################################
+TEMPLATE = app
+DESTDIR = ../../bin
+TARGET = T2Rho
+CONFIG -= app_bundle
+
+INCLUDEPATH += ../../include/common ../../include/T2Rho
+
+LIBS += -L../common -lcommon
+
+SOURCES += T2Rho.cpp
+
+HEADERS += ../../include/T2Rho/T2Rho.h