T2rho (#2)

* Unit display according to SI specs

* Fix  bug in argument -t

* Printing XFe in rock composition tabel

* Update binary location

* Added T2Rho

* Update README.md

* Update travis badge URL

* Update CHANGELOG

.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

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
 

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

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.
 

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

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

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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************/
+#ifndef T2RHO_H_
+#define T2RHO_H_
+
+#include <QList>
+#include <QString>
+#include <QStringList>
+#include <QDateTime>
+#include <QRegExp>
+#include <math.h>
+#include <iostream>
+#include <stdlib.h>
+#include "PhysicalConstants.h"
+#include "PointClasses.h"
+#include "ERMs.h"
+
+class T2Rho {
+  QString PMethod;
+  QString file_in;
+  QString file_out;
+  EarthReferenceModel * ERM;
+  // Mineral properties
+  QList<double> m_rho;
+  QList<double> m_drhodX;
+  QList<double> m_K;
+  QList<double> m_dKdT;
+  QList<double> m_dKdP;
+  QList<double> m_dKdPdX;
+  QList<double> m_dKdX;
+  QList<double> m_alpha0;
+  QList<Point5D> data_out;
+  // Rock properties
+  double r_XFe;
+  QList<double> r_comp;
+  // Functions
+  bool SetPMethod(QString method);
+  void argsError(QString val, bool ok);
+  void setComp(QList<double> composition);
+  void setComp(int c);
+
+ public:
+  T2Rho();
+  void usage();
+  void info();
+  void readArgs(int &argc, char *argv[]);
+  void readFile();
+  void writeFile();
+};
+
+#endif  // T2RHO_H_

include/common/PointClasses.h

@@ -20,81 +20,85 @@
 #define POINTCLASSES_H_
 
 class Point2D {
+private:
   QList <double> 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 <double> 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 <double> 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 <double> 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 <double> 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 <double> 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 <double> 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_