CoolProp/include/IdealCurves.h

#include "AbstractState.h"
#include "crossplatform_shared_ptr.h"
#include "Solvers.h"
#include "CoolPropTools.h"
#include <string>

namespace CoolProp{

	class CurveTracer : public FuncWrapper1D
	{
	public:
		AbstractState *AS;
        double p0, T0, lnT, lnp, rho_guess;
		std::vector<double> T, p;
		enum OBJECTIVE_TYPE {OBJECTIVE_INVALID = 0, OBJECTIVE_CIRCLE, OBJECTIVE_T };
		OBJECTIVE_TYPE obj;
		CurveTracer(AbstractState *AS, double p0, double T0) : AS(AS), p0(p0), T0(T0), lnT(_HUGE), lnp(_HUGE), rho_guess(_HUGE), obj(OBJECTIVE_INVALID)
		{
			this->p.push_back(p0);			
		};
		void init(){
			// Solve for Temperature for first point
			this->obj = OBJECTIVE_T;
			this->rho_guess = -1;
			this->T.push_back(Secant(this, T0, 0.001*T0, 1e-10, 100));
		}

		virtual double objective(void) = 0;

		virtual double starting_direction(){
			return M_PI / 2.0;
		}

		double call(double t){
			if (this->obj == OBJECTIVE_CIRCLE){
				double T2, P2;
				this->TPcoords(t, lnT, lnp, T2, P2);
				this->AS->update(PT_INPUTS, P2, T2);
			}
			else{
				if (this->rho_guess < 0)
					this->AS->update(PT_INPUTS, this->p[this->p.size()-1], t);
				else{
					GuessesStructure guesses;
					guesses.rhomolar = this->rho_guess;
					this->AS->update_with_guesses(PT_INPUTS, this->p[this->p.size() - 1], t, guesses);
				}
			}
			double r = this->objective();
			return r;
		}

		void TPcoords(double t, double lnT, double lnp, double &T, double &p){
			double rlnT = 0.1, rlnp = 0.1;
			T = exp(lnT + rlnT*cos(t));
			p = exp(lnp + rlnp*sin(t));
		}

		void trace(std::vector<double> &T, std::vector<double> &p)
		{
			double t = this->starting_direction();
			for (int i = 0; i < 1000; ++i){
				try{
					this->lnT = log(this->T[this->T.size() - 1]);
					this->lnp = log(this->p[this->p.size() - 1]);
					this->obj = OBJECTIVE_CIRCLE;
					t = Brent(this, t - M_PI / 2.0, t + M_PI / 2.0, DBL_EPSILON, 1e-10, 100);
					double T2, P2;
					this->TPcoords(t, this->lnT, this->lnp, T2, P2);
					this->T.push_back(T2);
					this->p.push_back(P2);
					if (this->T[this->T.size() - 1] < this->AS->keyed_output(iT_triple) || this->p[this->p.size() - 1] > 1000 * this->AS->keyed_output(iP_critical)){
						break;
					}
				}
				catch (std::exception &){
					break;
				}
			}
			T = this->T;
			p = this->p;
		}
	};

	class IdealCurveTracer : public CurveTracer{
	public:
		IdealCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
		/// Z = 1
		double objective(void){ return this->AS->keyed_output(iZ) - 1; };
	};

	class BoyleCurveTracer : public CurveTracer{
	public:
		BoyleCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
		/// dZ/dv|T = 0
		double objective(void){
			double r = (this->AS->p() - this->AS->rhomolar()*this->AS->first_partial_deriv(iP, iDmolar, iT)) / (this->AS->gas_constant()*this->AS->T());
			return r;
		};
	};
	class JouleInversionCurveTracer : public CurveTracer{
	public:
		JouleInversionCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
		/// dZ/dT|v = 0
		double objective(void){
			double r = (this->AS->gas_constant()*this->AS->T() * 1 / this->AS->rhomolar()*this->AS->first_partial_deriv(iP, iT, iDmolar) - this->AS->p()*this->AS->gas_constant() / this->AS->rhomolar()) / POW2(this->AS->gas_constant()*this->AS->T());
			return r;
		};
	};
	class JouleThomsonCurveTracer : public CurveTracer{
	public:
		JouleThomsonCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
		/// dZ/dT|p = 0
		double objective(void){
			double dvdT__constp = -this->AS->first_partial_deriv(iDmolar, iT, iP) / POW2(this->AS->rhomolar());
			double r = this->AS->p() / (this->AS->gas_constant()*POW2(this->AS->T()))*(this->AS->T()*dvdT__constp - 1 / this->AS->rhomolar());
			return r;
		};
	};

} /* namespace CoolProp */