--- LuaCoolProp core module. -- -- Copyright (C) 2026 Christophe Jorssen. -- Distributed under LPPL 1.3c or later. -- -- This file implements the Lua side of the `luacoolprop` package. The module is -- intentionally independent of any particular TeX format: it can be loaded from -- plain LuaTeX, LuaLaTeX, ConTeXt, or `texlua`. Its responsibilities are: -- -- * loading the CoolProp shared library through LuaTeX FFI; -- * exposing a small, idiomatic Lua wrapper around CoolProp's C ABI; -- * computing thermodynamic curves for registered diagram types; -- * emitting PGFPlots code for families of isolines and process paths; -- * validating thermodynamic states and reporting actionable diagnostics. -- -- Units at the Lua boundary follow CoolProp conventions unless explicitly -- stated otherwise: pressure in Pa, temperature in K, mass enthalpy in J/kg, -- mass entropy in J/(kg K), density in kg/m^3, and specific volume in m^3/kg. -- The PGFPlots layer normally rescales enthalpy to kJ/kg and pressure to bar. -- -- @module luacoolprop -- @author Christophe Jorssen -- @license LPPL-1.3c-or-later for LuaCoolProp code; CoolProp is a separate -- MIT-licensed dependency. -- @field _NAME Stable module name. -- @field _VERSION LuaCoolProp package version; this is independent of the -- loaded CoolProp version. local M = { _NAME = "luacoolprop", _VERSION = "1.1.0", } local ok_ffi, ffi_or_err = pcall(require, "ffi") if not ok_ffi then error("luacoolprop requires the LuaTeX/texlua FFI module; require('ffi') failed: " .. tostring(ffi_or_err), 0) end local ffi = ffi_or_err -- C declarations generated from CoolPropLib.h by scripts/update-coolprop.sh. -- They correspond to CoolProp's exported C ABI, not to the C++ API. -- BEGIN GENERATED COOLPROP FFI DECLARATIONS ffi.cdef[[ double Props1SI(const char* FluidName, const char* Output); void Props1SImulti(const char* Outputs, char* backend, const char* FluidNames, const double* fractions, const long length_fractions, double* result, long* resdim1); double PropsSI(const char* Output, const char* Name1, double Prop1, const char* Name2, double Prop2, const char* FluidName); void PropsSImulti(const char* Outputs, const char* Name1, double* Prop1, const long size_Prop1, const char* Name2, double* Prop2, const long size_Prop2, char* backend, const char* FluidNames, const double* fractions, const long length_fractions, double* result, long* resdim1, long* resdim2); long PhaseSI(const char* Name1, double Prop1, const char* Name2, double Prop2, const char* FluidName, char* phase, int n); long get_global_param_string(const char* param, char* Output, int n); long get_parameter_information_string(const char* param, char* Output, int n); long get_fluid_param_string(const char* fluid, const char* param, char* Output, int n); long get_fluid_param_string_len(const char* fluid, const char* param); void set_config_string(const char* key, const char* val); void set_config_double(const char* key, const double val); void set_config_bool(const char* key, const _Bool val); void set_departure_functions(const char* string_data, long* errcode, char* message_buffer, const long buffer_length); int set_reference_stateS(const char* Ref, const char* reference_state); int set_reference_stateD(const char* Ref, double T, double rhomolar, double hmolar0, double smolar0); void propssi_(const char* Output, const char* Name1, const double* Prop1, const char* Name2, const double* Prop2, const char* FluidName, double* output); double F2K(double T_F); double K2F(double T_K); long get_param_index(const char* param); long get_input_pair_index(const char* pair); long redirect_stdout(const char* file); int get_debug_level(); void set_debug_level(int level); double saturation_ancillary(const char* fluid_name, const char* output, int Q, const char* input, double value); double HAPropsSI(const char* Output, const char* Name1, double Prop1, const char* Name2, double Prop2, const char* Name3, double Prop3); double cair_sat(double T); void hapropssi_(const char* Output, const char* Name1, const double* Prop1, const char* Name2, const double* Prop2, const char* Name3, const double* Prop3, double* output); double HAProps(const char* Output, const char* Name1, double Prop1, const char* Name2, double Prop2, const char* Name3, double Prop3); void haprops_(const char* Output, const char* Name1, const double* Prop1, const char* Name2, const double* Prop2, const char* Name3, const double* Prop3, double* output); long AbstractState_factory(const char* backend, const char* fluids, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_fluid_names(const long handle, char* fluids, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_free(const long handle, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_set_fractions(const long handle, const double* fractions, const long N, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_get_mole_fractions(const long handle, double* fractions, const long maxN, long* N, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_get_mole_fractions_satState(const long handle, const char* saturated_state, double* fractions, const long maxN, long* N, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_get_fugacity(const long handle, const long i, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_get_fugacity_coefficient(const long handle, const long i, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_update(const long handle, const long input_pair, const double value1, const double value2, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_specify_phase(const long handle, const char* phase, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_unspecify_phase(const long handle, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_keyed_output(const long handle, const long param, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_first_saturation_deriv(const long handle, const long Of, const long Wrt, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_first_partial_deriv(const long handle, const long Of, const long Wrt, const long Constant, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_second_two_phase_deriv(const long handle, const long Of1, const long Wrt1, const long Constant1, const long Wrt2, const long Constant2, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_second_partial_deriv(const long handle, const long Of1, const long Wrt1, const long Constant1, const long Wrt2, const long Constant2, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_first_two_phase_deriv_splined(const long handle, const long Of, const long Wrt, const long Constant, const double x_end, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_first_two_phase_deriv(const long handle, const long Of, const long Wrt, const long Constant, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_update_and_common_out(const long handle, const long input_pair, const double* value1, const double* value2, const long length, double* T, double* p, double* rhomolar, double* hmolar, double* smolar, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_update_and_1_out(const long handle, const long input_pair, const double* value1, const double* value2, const long length, const long output, double* out, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_update_and_5_out(const long handle, const long input_pair, const double* value1, const double* value2, const long length, long* outputs, double* out1, double* out2, double* out3, double* out4, double* out5, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_set_binary_interaction_double(const long handle, const long i, const long j, const char* parameter, const double value, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_set_cubic_alpha_C(const long handle, const long i, const char* parameter, const double c1, const double c2, const double c3, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_set_fluid_parameter_double(const long handle, const long i, const char* parameter, const double value, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_build_phase_envelope(const long handle, const char* level, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_get_phase_envelope_data(const long handle, const long length, double* T, double* p, double* rhomolar_vap, double* rhomolar_liq, double* x, double* y, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_get_phase_envelope_data_checkedMemory(const long handle, const long length, const long maxComponents, double* T, double* p, double* rhomolar_vap, double* rhomolar_liq, double* x, double* y, long* actual_length, long* actual_components, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_build_spinodal(const long handle, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_get_spinodal_data(const long handle, const long length, double* tau, double* delta, double* M1, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_all_critical_points(const long handle, const long length, double* T, double* p, double* rhomolar, long* stable, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_keyed_output_satState(const long handle, const char* saturated_state, const long param, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_backend_name(const long handle, char* backend, long* errcode, char* message_buffer, const long buffer_length); void AbstractState_fluid_param_string(const long handle, const char* param, char* return_buffer, const long return_buffer_length, long* errcode, char* message_buffer, const long buffer_length); int AbstractState_phase(const long handle, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_saturated_liquid_keyed_output(const long handle, const long param, long* errcode, char* message_buffer, const long buffer_length); double AbstractState_saturated_vapor_keyed_output(const long handle, const long param, long* errcode, char* message_buffer, const long buffer_length); void add_fluids_as_JSON(const char* backend, const char* fluidstring, long* errcode, char* message_buffer, const long buffer_length); int C_is_valid_fluid_string(const char* fluidName); int C_extract_backend(const char* fluid_string, char* backend, const long backend_length, char* fluid, const long fluid_length); double PropsS(const char* Output, const char* Name1, double Prop1, const char* Name2, double Prop2, const char* Ref); double Props(const char* Output, const char Name1, double Prop1, const char Name2, double Prop2, const char* Ref); double Props1(const char* FluidName, const char* Output); ]] -- END GENERATED COOLPROP FFI DECLARATIONS local function getenv(name) local v = os.getenv(name) if v == nil or v == "" then return nil end return v end local function file_exists(path) if not path or path == "" then return false end local f = io.open(path, "rb") if f then f:close(); return true end return false end local function basename_is_path(s) return s and (s:find("/", 1, true) or s:find("\\", 1, true)) end local function tex_find_file(name) -- ConTeXt deliberately prevents direct kpse access; use its resolver when -- available. LuaLaTeX, plain LuaTeX, and texlua expose kpse.find_file. if type(resolvers) == "table" and type(resolvers.findfile) == "function" then local ok, found = pcall(resolvers.findfile, name) if ok and found and found ~= "" then return found end elseif type(kpse) == "table" and type(kpse.find_file) == "function" then local ok, found = pcall(kpse.find_file, name) if ok and found and found ~= "" then return found end end return nil end local function add_candidate(list, seen, candidate) if type(candidate) == "string" then candidate = candidate:match("^%s*(.-)%s*$") end if candidate and candidate ~= "" and not seen[candidate] then seen[candidate] = true list[#list + 1] = candidate end end local function library_candidates(user_path) local candidates, seen = {}, {} add_candidate(candidates, seen, user_path) add_candidate(candidates, seen, getenv("LUACOOLPROP_LIB")) add_candidate(candidates, seen, getenv("COOLPROP_LIB")) -- Local-directory candidates are useful for project-local and texlua-based workflows. add_candidate(candidates, seen, "./libCoolProp.so") add_candidate(candidates, seen, "./libCoolProp.dylib") add_candidate(candidates, seen, "./CoolProp.dll") add_candidate(candidates, seen, "./coolprop.dll") -- Format-aware lookup for TeX trees. This is intentionally conservative: -- the resolver only knows platform-dependent shared-library extensions when -- they are actually installed in a searchable tree. for _, name in ipairs({"libCoolProp.so", "libCoolProp.dylib", "CoolProp.dll", "coolprop.dll"}) do add_candidate(candidates, seen, tex_find_file(name)) end -- Names delegated to the platform dynamic loader. ffi.load("CoolProp") -- maps to the usual platform-specific prefix/suffix rules. add_candidate(candidates, seen, "CoolProp") add_candidate(candidates, seen, "libCoolProp") return candidates end local lib = nil local loaded_from = nil local load_errors = {} --- Load the CoolProp shared library. -- The search order is: explicit path, `LUACOOLPROP_LIB`, `COOLPROP_LIB`, the -- current directory, the active TeX resolver, and finally the platform dynamic -- loader. -- @tparam[opt] string path Explicit library path. -- @treturn cdata CoolProp FFI library handle. -- @treturn string Location or loader name that succeeded. function M.load_library(path) if lib ~= nil then return lib, loaded_from end for _, candidate in ipairs(library_candidates(path)) do local should_try = true if basename_is_path(candidate) and not file_exists(candidate) then should_try = false load_errors[#load_errors + 1] = candidate .. " : file not found" end if should_try then local ok, handle = pcall(ffi.load, candidate) if ok then lib = handle loaded_from = candidate M.C = lib return lib, loaded_from else load_errors[#load_errors + 1] = candidate .. " : " .. tostring(handle) end end end error("luacoolprop: unable to load CoolProp shared library. " .. "Pass library={/absolute/path/to/the/library}, or set LUACOOLPROP_LIB " .. "to libCoolProp.so, libCoolProp.dylib, or CoolProp.dll.\n" .. table.concat(load_errors, "\n"), 0) end --- Return the library path or loader name used by `load_library`. -- @treturn string|nil Loaded library path, or nil if the library has not yet been loaded. function M.loaded_library() return loaded_from end local function C() return M.load_library() end local function cstr(s, argname) if type(s) ~= "string" then error((argname or "argument") .. " must be a string", 3) end return s end local function number(x, argname) if type(x) ~= "number" then error((argname or "argument") .. " must be a number", 3) end return x end local DEFAULT_BUFFER_LENGTH = 4096 local function buffer(n) n = n or DEFAULT_BUFFER_LENGTH return ffi.new("char[?]", n), n end local function as_lua_string(buf) return ffi.string(buf) end local function err_context(n) local err = ffi.new("long[1]", 0) local msg, len = buffer(n or DEFAULT_BUFFER_LENGTH) return err, msg, len end local function check_err(err, msg, context) local code = tonumber(err[0]) if code ~= 0 then local text = as_lua_string(msg) if text == "" then text = "CoolProp error code " .. tostring(code) end error("luacoolprop: " .. (context or "CoolProp call") .. " failed: " .. text, 3) end end --- Return a CoolProp global string parameter. -- @tparam string param CoolProp global parameter name, for example `version`. -- @tparam[opt=4096] number n Buffer length. -- @treturn string Parameter value. function M.global_param_string(param, n) local out, len = buffer(n or DEFAULT_BUFFER_LENGTH) local ret = C().get_global_param_string(cstr(param, "param"), out, len) if tonumber(ret) < 0 then error("luacoolprop: get_global_param_string failed for " .. param, 2) end return as_lua_string(out) end --- Return the CoolProp version string. -- @treturn string CoolProp version. function M.version() return M.global_param_string("version") end --- Return the Git revision embedded in the loaded CoolProp library. -- This value identifies the external runtime dependency; it is not the -- LuaCoolProp revision. -- @treturn string CoolProp Git revision, as reported by CoolProp. function M.gitrevision() return M.global_param_string("gitrevision") end --- Return descriptive information for a CoolProp parameter. -- Parameter names are CoolProp identifiers such as `Hmass`, `P`, or `T`. -- The exact information strings are defined by the loaded CoolProp version. -- @tparam string param CoolProp parameter identifier. -- @tparam[opt=4096] number n Output-buffer length in bytes. -- @treturn string Parameter information supplied by CoolProp. function M.parameter_information_string(param, n) local out, len = buffer(n or DEFAULT_BUFFER_LENGTH) local ret = C().get_parameter_information_string(cstr(param, "param"), out, len) if tonumber(ret) < 0 then error("luacoolprop: get_parameter_information_string failed for " .. param, 2) end return as_lua_string(out) end --- Return a string-valued property of a CoolProp fluid. -- Common parameters include `aliases`, `CAS`, and `formula`; available names -- are controlled by the loaded CoolProp version. -- @tparam string fluid CoolProp fluid name. -- @tparam string param String parameter name. -- @tparam[opt=4096] number n Output-buffer length in bytes. -- @treturn string Fluid metadata supplied by CoolProp. function M.fluid_param_string(fluid, param, n) local out, len = buffer(n or DEFAULT_BUFFER_LENGTH) local ret = C().get_fluid_param_string(cstr(fluid, "fluid"), cstr(param, "param"), out, len) if tonumber(ret) < 0 then error("luacoolprop: get_fluid_param_string failed for " .. fluid .. ", " .. param, 2) end return as_lua_string(out) end --- Return the buffer size requested for a string-valued fluid property. -- Use this only when the default 4096-byte buffer accepted by -- `fluid_param_string` is insufficient. -- @tparam string fluid CoolProp fluid name. -- @tparam string param String parameter name. -- @treturn number Required length in bytes as reported by CoolProp. function M.fluid_param_string_len(fluid, param) return tonumber(C().get_fluid_param_string_len(cstr(fluid, "fluid"), cstr(param, "param"))) end --- Evaluate a thermodynamic property with CoolProp `PropsSI`. -- @tparam string output Requested property key. -- @tparam string name1 First input property key. -- @tparam number prop1 First input value in SI units. -- @tparam string name2 Second input property key. -- @tparam number prop2 Second input value in SI units. -- @tparam string fluid CoolProp fluid name. -- @treturn number Requested property in SI units. function M.propsSI(output, name1, prop1, name2, prop2, fluid) return tonumber(C().PropsSI( cstr(output, "output"), cstr(name1, "name1"), number(prop1, "prop1"), cstr(name2, "name2"), number(prop2, "prop2"), cstr(fluid, "fluid") )) end --- Alternate spelling for `propsSI` using CoolProp's capitalization. -- @function PropsSI M.PropsSI = M.propsSI --- Evaluate a one-argument CoolProp fluid property with `Props1SI`. -- @tparam string fluid CoolProp fluid name. -- @tparam string output Requested property key. -- @treturn number Requested property in SI units. function M.props1SI(fluid, output) return tonumber(C().Props1SI(cstr(fluid, "fluid"), cstr(output, "output"))) end --- Alternate spelling for `props1SI` using CoolProp's capitalization. -- @function Props1SI M.Props1SI = M.props1SI --- Evaluate a humid-air property with CoolProp `HAPropsSI`. -- All numeric values and the result use CoolProp's SI-unit convention. -- @tparam string output Requested humid-air property key. -- @tparam string name1 First input property key. -- @tparam number prop1 First input value in SI units. -- @tparam string name2 Second input property key. -- @tparam number prop2 Second input value in SI units. -- @tparam string name3 Third input property key. -- @tparam number prop3 Third input value in SI units. -- @treturn number Requested property in SI units. function M.hapropsSI(output, name1, prop1, name2, prop2, name3, prop3) return tonumber(C().HAPropsSI( cstr(output, "output"), cstr(name1, "name1"), number(prop1, "prop1"), cstr(name2, "name2"), number(prop2, "prop2"), cstr(name3, "name3"), number(prop3, "prop3") )) end --- Alternate spelling for `hapropsSI` using CoolProp's capitalization. -- @function HAPropsSI M.HAPropsSI = M.hapropsSI --- Return the phase name for a state defined by two independent properties. -- The input names, values, and returned phase vocabulary are defined by -- CoolProp. Numeric inputs must be expressed in SI units. -- @tparam string name1 First input property key. -- @tparam number prop1 First input value in SI units. -- @tparam string name2 Second input property key. -- @tparam number prop2 Second input value in SI units. -- @tparam string fluid CoolProp fluid name. -- @tparam[opt=4096] number n Output-buffer length in bytes. -- @treturn string CoolProp phase name. function M.phaseSI(name1, prop1, name2, prop2, fluid, n) local out, len = buffer(n or DEFAULT_BUFFER_LENGTH) local ret = C().PhaseSI( cstr(name1, "name1"), number(prop1, "prop1"), cstr(name2, "name2"), number(prop2, "prop2"), cstr(fluid, "fluid"), out, len) if tonumber(ret) < 0 then error("luacoolprop: PhaseSI failed", 2) end return as_lua_string(out) end --- Alternate spelling for `phaseSI` using CoolProp's capitalization. -- @function PhaseSI M.PhaseSI = M.phaseSI --- Test whether CoolProp recognizes a fluid specification. -- @tparam string fluid CoolProp fluid or backend-qualified fluid string. -- @treturn boolean True when the loaded library accepts the string. function M.is_valid_fluid_string(fluid) return C().C_is_valid_fluid_string(cstr(fluid, "fluid")) ~= 0 end --- Resolve a CoolProp output-parameter identifier to its numeric index. -- Numeric indices are version-specific and must not be persisted. -- @tparam string param CoolProp parameter identifier. -- @treturn number Numeric parameter index, or CoolProp's negative sentinel. function M.param_index(param) return tonumber(C().get_param_index(cstr(param, "param"))) end --- Resolve a CoolProp input-pair identifier to its numeric index. -- Numeric indices are version-specific and must not be persisted. -- @tparam string pair Input-pair identifier, for example `PT_INPUTS`. -- @treturn number Numeric input-pair index, or CoolProp's negative sentinel. function M.input_pair_index(pair) return tonumber(C().get_input_pair_index(cstr(pair, "pair"))) end --- Set a string-valued global CoolProp configuration option. -- This mutates process-wide CoolProp state and therefore affects subsequent -- calls made through every LuaCoolProp object in the same Lua process. -- @tparam string key CoolProp configuration key. -- @tparam string value New string value. function M.set_config_string(key, value) C().set_config_string(cstr(key, "key"), cstr(value, "value")) end --- Set a numeric global CoolProp configuration option. -- This mutates process-wide CoolProp state. -- @tparam string key CoolProp configuration key. -- @tparam number value New numeric value. function M.set_config_double(key, value) C().set_config_double(cstr(key, "key"), number(value, "value")) end --- Set a boolean global CoolProp configuration option. -- Lua truthiness is used: only `false` and `nil` map to C false. This mutates -- process-wide CoolProp state. -- @tparam string key CoolProp configuration key. -- @tparam boolean value New boolean value. function M.set_config_bool(key, value) C().set_config_bool(cstr(key, "key"), value and true or false) end --- Set CoolProp's process-wide diagnostic level. -- @tparam number level Numeric debug level accepted by CoolProp. function M.set_debug_level(level) C().set_debug_level(assert(tonumber(level), "level must be numeric")) end --- Return CoolProp's current process-wide diagnostic level. -- @treturn number Current numeric debug level. function M.get_debug_level() return tonumber(C().get_debug_level()) end --- Evaluate a saturation ancillary correlation. -- This is a direct binding to CoolProp's C API. Inputs and output use the SI -- convention associated with the selected CoolProp property identifiers. -- @tparam string fluid CoolProp fluid name. -- @tparam string output Requested output property. -- @tparam number Q Saturated branch selector, normally 0 or 1. -- @tparam string input Input property identifier. -- @tparam number value Input value in SI units. -- @treturn number Ancillary result in SI units. function M.saturation_ancillary(fluid, output, Q, input, value) return tonumber(C().saturation_ancillary( cstr(fluid, "fluid"), cstr(output, "output"), assert(tonumber(Q), "Q must be numeric"), cstr(input, "input"), number(value, "value") )) end --- Managed CoolProp AbstractState method table. -- Instances own exactly one native handle and are constructed with -- `M.AbstractState`. -- @type State local State = {} State.__index = State --- Release this state's CoolProp handle. -- The operation is idempotent. Prefer an explicit call when the state is no -- longer needed; the FFI finalizer is only a safety net whose execution time is -- deliberately unspecified by Lua. -- @treturn nil function State:free() if self.handle ~= nil and self.handle >= 0 then local err, msg, len = err_context() C().AbstractState_free(self.handle, err, msg, len) self.handle = nil if self._owner ~= nil then self._owner.handle = -1 end check_err(err, msg, "AbstractState_free") end end --- Return the fluid specification owned by this state. -- @treturn string CoolProp fluid or mixture string. function State:fluid_names() local err, msg, len = err_context() local out = buffer(DEFAULT_BUFFER_LENGTH) C().AbstractState_fluid_names(self.handle, out, err, msg, len) check_err(err, msg, "AbstractState_fluid_names") return as_lua_string(out) end --- Update this state from a CoolProp input pair. -- String pair identifiers are resolved for each call; callers performing a hot -- loop may resolve the numeric index once with `input_pair_index`. -- @tparam string|number input_pair CoolProp input-pair name or numeric index. -- @tparam number value1 First input in the pair's SI units. -- @tparam number value2 Second input in the pair's SI units. -- @treturn State This state, to allow method chaining. function State:update(input_pair, value1, value2) local pair = type(input_pair) == "string" and M.input_pair_index(input_pair) or tonumber(input_pair) local err, msg, len = err_context() C().AbstractState_update(self.handle, pair, number(value1, "value1"), number(value2, "value2"), err, msg, len) check_err(err, msg, "AbstractState_update") return self end --- Read one property from the current state. -- @tparam string|number param CoolProp parameter name or numeric index. -- @treturn number Property value in CoolProp SI units. function State:keyed_output(param) local p = type(param) == "string" and M.param_index(param) or tonumber(param) local err, msg, len = err_context() local out = C().AbstractState_keyed_output(self.handle, p, err, msg, len) check_err(err, msg, "AbstractState_keyed_output") return tonumber(out) end --- Return the numeric CoolProp phase index for the current state. -- Use `phaseSI` when a human-readable phase string is required. -- @treturn number CoolProp phase index. function State:phase() local err, msg, len = err_context() local out = C().AbstractState_phase(self.handle, err, msg, len) check_err(err, msg, "AbstractState_phase") return tonumber(out) end local state_gc_ct = ffi.metatype("struct { long handle; }", { __gc = function(obj) if obj.handle ~= nil and obj.handle >= 0 and lib ~= nil then local err, msg, len = err_context() pcall(function() lib.AbstractState_free(obj.handle, err, msg, len) end) obj.handle = -1 end end }) --- Construct a managed wrapper around CoolProp `AbstractState`. -- The returned object releases its handle through `:free()`; when supported by -- the LuaTeX engine, a finalizer is also attached as a safety net. Callers -- should still use deterministic cleanup and must not access the object after -- `:free()`. -- @tparam string backend CoolProp backend, for example `HEOS`. -- @tparam string fluids Fluid or mixture specification. -- @treturn State Managed AbstractState wrapper. function M.AbstractState(backend, fluids) local err, msg, len = err_context() local handle = C().AbstractState_factory(cstr(backend, "backend"), cstr(fluids, "fluids"), err, msg, len) check_err(err, msg, "AbstractState_factory") if tonumber(handle) < 0 then error("luacoolprop: AbstractState_factory returned an invalid handle", 2) end -- Keep a cdata owner for finalization. CoolProp handles are integer -- indices and may legitimately start at 0. A negative value is used only -- on the Lua side as the "already freed" sentinel. local owner = state_gc_ct(handle) local self = setmetatable({ handle = tonumber(handle), _owner = owner }, State) return self end --- Expose the managed AbstractState method table for introspection. -- Construct instances with `AbstractState`; do not call methods with a -- hand-built table. -- @table State M.State = State -- -------------------------------------------------------------------------- -- Thermodynamic diagram helpers -- -------------------------------------------------------------------------- M.diagram = M.diagram or {} M.tex = M.tex or {} M.diagram.types = M.diagram.types or {} --- Register a thermodynamic diagram implementation. -- Diagram codes are uppercase coordinate pairs such as `PH`, `PV`, and `TS`. -- Implementations expose format-independent renderers and a `families` table. -- @tparam string diagram_type Canonical uppercase diagram code. -- @tparam table implementation Diagram implementation table. function M.diagram.register_type(diagram_type, implementation) if type(diagram_type) ~= "string" or not diagram_type:match("^[A-Z][A-Z0-9]*$") then error("luacoolprop: diagram type must be an uppercase identifier such as PH, PV, or TS", 2) end if type(implementation) ~= "table" then error("luacoolprop: diagram implementation for " .. diagram_type .. " must be a table", 2) end if M.diagram.types[diagram_type] ~= nil then error("luacoolprop: diagram type " .. diagram_type .. " is already registered", 2) end implementation.type = diagram_type M.diagram.types[diagram_type] = implementation return implementation end --- Return a registered thermodynamic diagram implementation. -- @tparam string diagram_type Canonical uppercase diagram code. -- @treturn table Diagram implementation table. function M.diagram.get_type(diagram_type) local implementation = M.diagram.types[diagram_type] if implementation == nil then local available = {} for name in pairs(M.diagram.types) do available[#available + 1] = name end table.sort(available) local suffix = #available > 0 and "; available types: " .. table.concat(available, ", ") or "" error("luacoolprop: unknown diagram type '" .. tostring(diagram_type) .. "'" .. suffix, 2) end return implementation end local function finite_number(x) return type(x) == "number" and x == x and x ~= math.huge and x ~= -math.huge end -- Keep related parser and consistency helpers in one private namespace. Apart -- from making their role explicit, this avoids Lua 5.3's limit of 200 locals -- per function for this intentionally self-contained runtime module. M._strict = {} -- Parse only conventional decimal notation. Lua's tonumber also accepts -- implementation-specific spellings such as hexadecimal numbers; TeX-facing -- keys deliberately do not. A full stop is the only decimal separator, while -- either `e` or `E` introduces an optional signed decimal exponent. function M._strict.decimal_number(value) if type(value) == "number" then if finite_number(value) then return value end return nil end if type(value) ~= "string" then return nil end local text = value:gsub("^%s+", ""):gsub("%s+$", "") if text:find(",", 1, true) then return nil end local valid = text:match("^[%+%-]?%d+%.?%d*[eE][%+%-]?%d+$") or text:match("^[%+%-]?%.%d+[eE][%+%-]?%d+$") or text:match("^[%+%-]?%d+%.?%d*$") or text:match("^[%+%-]?%.%d+$") if not valid then return nil end local number_value = tonumber(text) if not finite_number(number_value) then return nil end return number_value end function M._strict.require_decimal(value, label) local parsed = M._strict.decimal_number(value) if parsed == nil then error("luacoolprop: " .. tostring(label or "value") .. " must be a full finite decimal number using '.'; scientific notation is accepted", 3) end return parsed end function M._strict.decimal_number_prefix(text) local patterns = { "^([%+%-]?%d+%.?%d*[eE][%+%-]?%d+)(.*)$", "^([%+%-]?%.%d+[eE][%+%-]?%d+)(.*)$", "^([%+%-]?%d+%.?%d*)(.*)$", "^([%+%-]?%.%d+)(.*)$", } for _, pattern in ipairs(patterns) do local numeric_text, suffix = text:match(pattern) if numeric_text then local numeric_value = M._strict.decimal_number(numeric_text) if numeric_value then return numeric_value, suffix, numeric_text end end end return nil end local function to_number(x, default) if x == nil or x == "" or x == "auto" then return default end local n = M._strict.decimal_number(x) if n == nil then error("luacoolprop: invalid numeric value '" .. tostring(x) .. "'; use a full finite decimal number with '.' as the decimal " .. "separator (scientific notation such as 1.2E5 is accepted)", 3) end return n end local function to_string(x, default) if x == nil or x == "" then return default end return tostring(x) end local function to_bool(x, default) if x == nil or x == "" or x == "auto" then return default end if type(x) == "boolean" then return x end x = tostring(x):lower() if x == "true" or x == "yes" or x == "on" or x == "1" then return true end if x == "false" or x == "no" or x == "off" or x == "0" then return false end error("luacoolprop: invalid boolean value '" .. tostring(x) .. "'; expected true or false (aliases: yes/no, on/off, 1/0)", 3) end function M._strict.choice_value(value, default, label, accepted) local selected = to_string(value, default):lower() for _, candidate in ipairs(accepted) do if selected == candidate then return selected end end error("luacoolprop: invalid " .. label .. " '" .. tostring(value) .. "'; expected one of: " .. table.concat(accepted, ", "), 3) end function M._strict.domain_options(source, target, context) target.domain_policy = source.domain_policy or source.discontinuity_policy target.domain_context = context return target end -- Reference-state changes affect every later CoolProp call for a fluid. The -- diagram layer therefore treats them as initialization, records the selected -- convention, and locks it as soon as the first diagram calculation starts. M._reference_states = {} M._reference_locks = {} --- Select a standard CoolProp enthalpy/entropy reference state. -- Call this during initialization, before creating any diagram for `fluid`. -- Repeating the current choice is harmless; changing a locked choice is an -- error because it would make one document contain incompatible coordinates. -- @tparam string fluid Pure-fluid CoolProp identifier. -- @tparam[opt=DEF] string reference One of `DEF`, `IIR`, `ASHRAE`, or `NBP`. -- @tparam[opt] string library Explicit shared-library path. -- @treturn string Normalized selected reference-state name. function M.set_reference_state(fluid, reference, library) M.load_library(library) fluid = cstr(fluid, "fluid") local base = fluid:match("^[^:]+::(.+)$") or fluid local selected = M._strict.choice_value(reference, "DEF", "reference state", {"def", "iir", "ashrae", "nbp"}):upper() local current = M._reference_states[base] if M._reference_locks[base] and current ~= selected then error("luacoolprop: reference state for '" .. base .. "' is locked to " .. tostring(current) .. "; select " .. selected .. " during initialization before the first diagram calculation", 2) end if current ~= selected then local result = tonumber(C().set_reference_stateS(cstr(base, "fluid"), cstr(selected, "reference state"))) if result ~= 1 then local ok, detail = pcall(M.global_param_string, "errstring") error("luacoolprop: CoolProp rejected reference state " .. selected .. " for '" .. base .. "'" .. (ok and detail ~= "" and (": " .. detail) or ""), 2) end M._reference_states[base] = selected end return selected end --- Return the reference state selected through LuaCoolProp for a fluid. -- `nil` means no LuaCoolProp diagram has initialized this fluid yet. -- @tparam string fluid CoolProp fluid identifier. -- @treturn string|nil `DEF`, `IIR`, `ASHRAE`, `NBP`, or nil. function M.reference_state(fluid) fluid = cstr(fluid, "fluid") local base = fluid:match("^[^:]+::(.+)$") or fluid return M._reference_states[base] end --- Reject mixtures and initialize the reference state for a diagram call. -- Low-level wrappers remain mixture-capable; this guard applies only to the -- diagram layer, whose saturation and critical-point algorithms require one -- CoolProp fluid reported as pure, with one quality dome. -- @tparam[opt={}] table opts Diagram options containing `fluid`, `library`, -- and optionally `reference_state`. -- @treturn string Validated fluid identifier. function M.diagram.prepare_fluid(opts) opts = opts or {} M.load_library(opts.library) local fluid = to_string(opts.fluid, "R134a") local base = fluid:match("^[^:]+::(.+)$") or fluid local syntactic_mixture = base:find("&", 1, true) or base:lower():match("%.mix$") or base:find("[", 1, true) local pure = false if not syntactic_mixture then local ok, value = pcall(M.fluid_param_string, base, "pure") pure = ok and tostring(value):lower() == "true" end if not pure then error("luacoolprop: diagrams require a pure fluid; mixture or unsupported " .. "fluid string '" .. fluid .. "' was rejected. Low-level CoolProp " .. "wrappers remain available for deliberate mixture calculations", 2) end local selected = opts.reference_state or M._reference_states[base] or "DEF" M.set_reference_state(base, selected) M._reference_locks[base] = true return fluid end M.diagram.axis_properties = { PH = {enthalpy = true, pressure = true}, PV = {specific_volume = true, pressure = true}, TS = {entropy = true, temperature = true}, HS = {entropy = true, enthalpy = true}, PT = {temperature = true, pressure = true}, } M.diagram.axis_units = { pressure = {option = "pressure_axis_unit", scale = "p_scale", values = {pa = 1, kpa = 1e-3, mpa = 1e-6, bar = 1e-5}}, enthalpy = {option = "enthalpy_axis_unit", scale = "h_scale", values = {jkg = 1, kjkg = 1e-3}}, entropy = {option = "entropy_axis_unit", scale = "s_scale", values = {jkgk = 1, kjkgk = 1e-3}}, specific_volume = {option = "specific_volume_axis_unit", scale = "v_scale", values = {m3kg = 1, lkg = 1000}}, temperature = {option = "temperature_axis_unit", scale = "temperature_scale", values = {kelvin = 1, k = 1, celsius = 1, c = 1, degc = 1}}, } M.diagram.family_switches = { phase_envelope = {"phase_envelope", "phase_envelopes", "saturation"}, quality = {"quality", "isoquality"}, isotherm = {"isotherm", "isotherms"}, isentrope = {"isentrope", "isentropes"}, isochore = {"isochore", "isochores"}, isenthalp = {"isenthalp", "isenthalps"}, isobar = {"isobar", "isobars"}, } --- Normalize shared options and enforce a diagram's declared capabilities. -- Every registered type uses this same boundary; no diagram is the template -- or fallback implementation for another one. function M.diagram.normalize_options(diagram_type, opts, aliases) local out = {} for key, value in pairs(opts or {}) do out[key] = value end for canonical, internal in pairs(aliases or {}) do if out[canonical] ~= nil then out[internal] = out[canonical] end end local implementation = M.diagram.get_type(diagram_type) for family, switches in pairs(M.diagram.family_switches) do if not (implementation.families or {})[family] then for _, switch in ipairs(switches) do if out[switch] ~= nil and to_bool(out[switch], false) then error("luacoolprop: diagram type " .. diagram_type .. " does not provide the '" .. family .. "' background family; use a process path if that conserved " .. "property is required", 2) end end end end local coordinates = M.diagram.axis_properties[diagram_type] or {} for property, specification in pairs(M.diagram.axis_units) do local raw = out[specification.option] if raw ~= nil and raw ~= "" and raw ~= "auto" then if not coordinates[property] then error("luacoolprop: " .. specification.option:gsub("_", " ") .. " is not an axis option of diagram type " .. diagram_type, 2) end local unit = tostring(raw):lower():gsub("[%s/%^%-%(%)]", "") local factor = specification.values[unit] if not factor then local accepted = {} for name in pairs(specification.values) do accepted[#accepted + 1] = name end table.sort(accepted) error("luacoolprop: invalid " .. specification.option:gsub("_", " ") .. " '" .. tostring(raw) .. "'; expected one of: " .. table.concat(accepted, ", "), 2) end out[specification.option] = unit out[specification.scale] = factor if property == "temperature" then out.temperature_offset = (unit == "celsius" or unit == "c" or unit == "degc") and 273.15 or 0 end end end return out end local function log10(x) return math.log(x) / math.log(10) end local function trim_trailing_zeros(s) s = s:gsub("(%..-)0+$", "%1") s = s:gsub("%.$", "") if s == "-0" then s = "0" end return s end local function fmt_num(x, digits) digits = digits or 6 if not finite_number(x) then return "nan" end local ax = math.abs(x) local s if ax ~= 0 and (ax < 1e-4 or ax >= 1e6) then s = string.format("%." .. tostring(math.max(1, digits)) .. "e", x) s = s:gsub("0+e", "e") else s = string.format("%." .. tostring(math.max(0, digits)) .. "f", x) s = trim_trailing_zeros(s) end return s end local function fmt_q(q) if math.abs(q) < 5e-13 then return "0" end if math.abs(q - 1) < 5e-13 then return "1" end return trim_trailing_zeros(string.format("%.3f", q)) end local function fmt_t_value(v) local av = math.abs(v) if av >= 100 or math.abs(v - math.floor(v + 0.5)) < 2e-8 then return trim_trailing_zeros(string.format("%.0f", v)) end return trim_trailing_zeros(string.format("%.1f", v)) end local function fmt_s_value(v) local av = math.abs(v) if av >= 10 then return trim_trailing_zeros(string.format("%.1f", v)) elseif av >= 1 then return trim_trailing_zeros(string.format("%.2f", v)) end return trim_trailing_zeros(string.format("%.3f", v)) end local function fmt_v_value(v) local av = math.abs(v) if av ~= 0 and (av < 1e-3 or av >= 1) then local s = string.format("%.3g", v) return s:gsub("e%+", "e") elseif av < 0.01 then return trim_trailing_zeros(string.format("%.4f", v)) elseif av < 0.1 then return trim_trailing_zeros(string.format("%.3f", v)) end return trim_trailing_zeros(string.format("%.2f", v)) end local function append_coord(out, point, digits) if point._break_before and #out > 0 then out[#out + 1] = "(nan,nan)" end out[#out + 1] = "(" .. fmt_num(point.x, digits) .. "," .. fmt_num(point.y, digits) .. ")" end local function parse_number_list(s) if s == nil or s == "" or s == "auto" then return nil end local values = {} for token in tostring(s):gmatch("[^,%s;]+") do local n = M._strict.decimal_number(token) if not n then error("luacoolprop: invalid numeric list element '" .. token .. "'; use '.' as the decimal separator and separate values with " .. "commas, spaces, or semicolons", 3) end values[#values + 1] = n end if #values == 0 then return nil end return values end local function make_linear_values(vmin, vmax, vstep) local values = {} if vstep <= 0 then error("luacoolprop: step must be positive", 3) end if vmax < vmin then error("luacoolprop: grid maximum must be greater than or equal to its minimum", 3) end local n = math.floor((vmax - vmin) / vstep + 0.5) for i = 0, n do local v = vmin + i * vstep if v >= vmin - 1e-10 and v <= vmax + 1e-10 then values[#values + 1] = v end end if #values == 0 or math.abs(values[#values] - vmax) > 1e-9 then values[#values + 1] = vmax end return values end local function make_count_values(vmin, vmax, count) count = count or 1 if not finite_number(count) or count < 1 or count ~= math.floor(count) then error("luacoolprop: grid count must be a positive integer", 3) end if vmax < vmin then error("luacoolprop: grid maximum must be greater than or equal to its minimum", 3) end if count == 1 then return {0.5 * (vmin + vmax)} end local values = {} for i = 0, count - 1 do values[#values + 1] = vmin + (vmax - vmin) * i / (count - 1) end return values end local function make_log_values(vmin, vmax, count) count = count or 1 if not finite_number(count) or count < 1 or count ~= math.floor(count) then error("luacoolprop: logarithmic grid count must be a positive integer", 3) end if vmax < vmin then error("luacoolprop: logarithmic grid maximum must be greater than or equal to its minimum", 3) end if vmin <= 0 or vmax <= 0 then error("luacoolprop: logarithmic value grids require strictly positive bounds", 3) end if count == 1 then return {math.sqrt(vmin * vmax)} end local lmin, lmax = log10(vmin), log10(vmax) local values = {} for i = 0, count - 1 do values[#values + 1] = 10 ^ (lmin + (lmax - lmin) * i / (count - 1)) end return values end local function normalize_sorted_unique(values, snap) table.sort(values) local out = {} for _, v in ipairs(values) do if snap then v = snap(v) end if #out == 0 or math.abs(v - out[#out]) > 1e-10 * math.max(1, math.abs(v)) then out[#out + 1] = v end end return out end local function make_grid_values(spec) local mode = M._strict.choice_value(spec.mode, "linear", "grid mode", { "linear", "log", "logarithmic", "list", "explicit", "preset", "predefined", "auto", }) local explicit = parse_number_list(spec.values) if explicit or mode == "list" or mode == "explicit" then if not explicit then error("luacoolprop: list grid mode requires at least one explicit numeric value", 3) end return normalize_sorted_unique(explicit, spec.snap) end if spec.preset_values and (mode == "preset" or mode == "predefined" or mode == "auto") then return normalize_sorted_unique({table.unpack(spec.preset_values)}, spec.snap) end local vmin, vmax = assert(spec.min, "missing grid minimum"), assert(spec.max, "missing grid maximum") local values if mode == "log" or mode == "logarithmic" then values = make_log_values(vmin, vmax, spec.count or 10) elseif spec.step and spec.step ~= "auto" then values = make_linear_values(vmin, vmax, spec.step) else values = make_count_values(vmin, vmax, spec.count or 10) end return normalize_sorted_unique(values, spec.snap) end local function make_q_values(opts) opts = opts or {} local qmin = to_number(opts.q_min, 0) local qmax = to_number(opts.q_max, 1) local qstep = to_number(opts.q_step, 0.1) local qcount = to_number(opts.q_count, nil) local qmode = to_string(opts.quality_mode or opts.q_mode, "linear") local qvalues = opts.quality_values or opts.q_values local qpreset = M._strict.choice_value(opts.quality_preset or opts.q_preset, "default", "quality preset", {"default", "dense", "sparse", "boundary", "boundaries"}) local preset_values if qpreset == "default" or qpreset == "dense" then preset_values = make_linear_values(0, 1, 0.1) elseif qpreset == "sparse" then preset_values = {0, 0.25, 0.5, 0.75, 1} elseif qpreset == "boundary" or qpreset == "boundaries" then preset_values = {0, 1} end local values = make_grid_values({ mode = qmode, values = qvalues, min = qmin, max = qmax, step = qstep, count = qcount, preset_values = preset_values, snap = function(q) if q < 0 and q > -1e-12 then q = 0 end if q > 1 and q < 1 + 1e-12 then q = 1 end if math.abs(q) < 5e-13 then q = 0 end if math.abs(q - 1) < 5e-13 then q = 1 end return q end, }) for _, q in ipairs(values) do if q < 0 or q > 1 then error("luacoolprop: quality grid values must lie between 0 and 1", 3) end end return values end local function join_style(...) local out = {} for i = 1, select("#", ...) do local s = select(i, ...) if s ~= nil then s = tostring(s) if s ~= "" and s ~= "none" then out[#out + 1] = s end end end return table.concat(out, ",") end local function slug_number(x, digits) digits = digits or 4 if math.abs(x) < 5e-12 then x = 0 end local s local ax = math.abs(x) if ax ~= 0 and (ax < 1e-3 or ax >= 1e5) then s = string.format("%." .. digits .. "e", x) s = s:gsub("0+e", "e") else s = trim_trailing_zeros(string.format("%." .. digits .. "f", x)) end s = s:gsub("%-", "m") :gsub("%+", "") :gsub("%.", "p") :gsub("e", "e") return s end local function curve_id(prefix, kind, value, suffix, digits) return table.concat({prefix or "lcp-ph", kind, slug_number(value, digits or 4) .. (suffix or "")}, "-") end local function short_curve_id(id) if type(id) ~= "string" then return id end return (id:gsub("^lcp%-ph%-", "")) end local function parse_curve_overrides(opts) if opts._curve_overrides_parsed then return opts._curve_overrides_parsed end local parsed = {} local s = opts.curve_overrides or opts.curve_override or "" for line in tostring(s):gmatch("[^\r\n]+") do local id, key, value = line:match("^%s*([^|]+)|([^|]+)|(.-)%s*$") if id and key then id = id:gsub("^%s+", ""):gsub("%s+$", "") key = key:gsub("^%s+", ""):gsub("%s+$", ""):gsub("[%s%-]+", "_") parsed[id] = parsed[id] or {} parsed[id][key] = value or "" local sid = short_curve_id(id) parsed[sid] = parsed[sid] or parsed[id] end end opts._curve_overrides_parsed = parsed return parsed end local function curve_override(opts, id) local ov = parse_curve_overrides(opts or {}) return ov[id] or ov[short_curve_id(id)] or {} end local function override_bool(ov, key, default) if ov and ov[key] ~= nil then return to_bool(ov[key], default) end return default end local function override_number(ov, key, default) if ov and ov[key] ~= nil then return to_number(ov[key], default) end return default end local function override_string(ov, key, default) if ov and ov[key] ~= nil and ov[key] ~= "" then return ov[key] end return default end local function plot_name_style(id) return "name path=" .. id end local function safe_propsSI(output, name1, prop1, name2, prop2, fluid) local ok, value = pcall(M.propsSI, output, name1, prop1, name2, prop2, fluid) if ok and finite_number(value) then return value end return nil, value end local function safe_props1SI(fluid, output, default) local ok, value = pcall(M.props1SI, fluid, output) if ok and finite_number(value) then return value end return default end --- Return basic constants exposed by CoolProp for a fluid string. -- This metadata helper is deliberately outside the diagram pure-fluid guard. -- @tparam string fluid CoolProp fluid name. -- @tparam[opt] string library Explicit shared-library path. -- @treturn table Constants including `pcrit`, `ptriple`, `Tcrit`, `Ttriple`, -- `Tmin`, `Tmax`, and `rhocrit` when available. function M.fluid_constants(fluid, library) M.load_library(library) fluid = fluid or "R134a" return { fluid = fluid, pcrit = safe_props1SI(fluid, "pcrit"), ptriple = safe_props1SI(fluid, "ptriple"), Tcrit = safe_props1SI(fluid, "Tcrit"), Ttriple = safe_props1SI(fluid, "Ttriple", safe_props1SI(fluid, "Tmin")), Tmin = safe_props1SI(fluid, "Tmin"), Tmax = safe_props1SI(fluid, "Tmax"), rhocrit = safe_props1SI(fluid, "rhomass_critical", safe_props1SI(fluid, "rhocrit")), } end local function resolve_pressure_limits(opts, fluid) M.diagram.prepare_fluid(opts) local c = M.fluid_constants(fluid, opts.library) local pcrit = c.pcrit or 1e6 local ptriple = c.ptriple or pcrit * 1e-4 local pmax_factor = to_number(opts.pressure_max_factor, 1.10) if pmax_factor <= 0 then error("luacoolprop: pressure_max_factor must be strictly positive", 3) end local pmin = to_number(opts.pressure_min, math.max(ptriple * 1.05, pcrit * 2e-3)) local pmax = to_number(opts.pressure_max, pcrit * pmax_factor) if pmin <= 0 then error("luacoolprop: pressure_min must be strictly positive", 3) end if pmax <= pmin then error("luacoolprop: pressure_max must be greater than pressure_min", 3) end return pmin, pmax, c end local function ph_point_from_logp_quality(fluid, q, t, h_scale, p_scale) local p_si = 10 ^ t local h_si = safe_propsSI("H", "P", p_si, "Q", q, fluid) if not h_si then return nil end local x = h_si * h_scale local y = p_si * p_scale if y <= 0 then return nil end return {x = x, y = y, p = p_si, h = h_si, t = t} end -- Return a single, quality-independent representation of the critical state. -- P,Q is singular at that state, so evaluating q=0 and q=1 independently -- cannot produce a reliable common endpoint. T,rho is well-defined for a -- pure fluid and gives both saturation branches exactly the same PH point. local function ph_critical_point(fluid, constants, h_scale, p_scale, fallback_pressure) if not constants or not finite_number(constants.pcrit) or not finite_number(constants.Tcrit) or not finite_number(constants.rhocrit) or constants.pcrit <= 0 or constants.Tcrit <= 0 or constants.rhocrit <= 0 then return nil end local h_si = safe_propsSI("H", "T", constants.Tcrit, "Dmass", constants.rhocrit, fluid) local approximated = false if not h_si and fallback_pressure then -- A few third-party CoolProp backends do not accept the exact critical -- (T,rho) pair. The two saturation enthalpies converge to the same -- value, so their mean at the closest common valid pressure is a stable, -- quality-independent fallback. The plotted pressure remains exactly pc. local h_liquid = safe_propsSI("H", "P", fallback_pressure, "Q", 0, fluid) local h_vapour = safe_propsSI("H", "P", fallback_pressure, "Q", 1, fluid) if h_liquid and h_vapour then h_si = 0.5 * (h_liquid + h_vapour) approximated = true end end if not h_si then return nil end return { x = h_si * h_scale, y = constants.pcrit * p_scale, p = constants.pcrit, h = h_si, T = constants.Tcrit, rho = constants.rhocrit, v = 1 / constants.rhocrit, t = log10(constants.pcrit), critical = true, critical_limit = true, quality_defined = false, critical_approximated = approximated, } end -- Find the closest pressure below pc at which both limiting saturation states -- can be evaluated. Testing both q=0 and q=1 makes the selected pressure -- independent of the particular quality curve being generated. The ordered -- fallback sequence accommodates different CoolProp backends without using a -- fluid-specific constant. local quality_critical_epsilons = { 1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, } local function quality_subcritical_pressure(fluid, pcrit) for _, epsilon in ipairs(quality_critical_epsilons) do local pressure = pcrit * (1 - epsilon) local h_liquid = safe_propsSI("H", "P", pressure, "Q", 0, fluid) local h_vapour = safe_propsSI("H", "P", pressure, "Q", 1, fluid) if h_liquid and h_vapour then return pressure end end return nil end local function ph_point_from_logp_temperature(fluid, T, t, h_scale, p_scale) local p_si = 10 ^ t local h_si = safe_propsSI("H", "T", T, "P", p_si, fluid) if not h_si then return nil end local x = h_si * h_scale local y = p_si * p_scale if y <= 0 then return nil end return {x = x, y = y, p = p_si, h = h_si, t = t, T = T} end local function ph_point_from_logp_entropy(fluid, s, t, h_scale, p_scale) local p_si = 10 ^ t local h_si = safe_propsSI("H", "P", p_si, "S", s, fluid) if not h_si then return nil end local x = h_si * h_scale local y = p_si * p_scale if y <= 0 then return nil end return {x = x, y = y, p = p_si, h = h_si, t = t, s = s} end local function ph_point_from_logp_specific_volume(fluid, v, t, h_scale, p_scale) if not v or v <= 0 then return nil end local p_si = 10 ^ t local rho_si = 1 / v local h_si = safe_propsSI("H", "P", p_si, "Dmass", rho_si, fluid) if not h_si then return nil end local x = h_si * h_scale local y = p_si * p_scale if y <= 0 then return nil end return {x = x, y = y, p = p_si, h = h_si, t = t, v = v, rho = rho_si} end local function pv_point_from_logp_pair(fluid, input, value, t, v_scale, p_scale) local p_si = 10 ^ t local rho_si = safe_propsSI("Dmass", "P", p_si, input, value, fluid) if not rho_si or rho_si <= 0 then return nil end local v_si = 1 / rho_si local point = { x = v_si * v_scale, y = p_si * p_scale, p = p_si, v = v_si, rho = rho_si, t = t, } if input == "Q" then point.q = value elseif input == "T" then point.T = value elseif input == "S" then point.s = value elseif input == "H" then point.h = value end return point end local function pv_saturation_point_from_T(fluid, T, q, v_scale, p_scale) local p_si = safe_propsSI("P", "T", T, "Q", q, fluid) local rho_si = safe_propsSI("Dmass", "T", T, "Q", q, fluid) if not p_si or not rho_si or p_si <= 0 or rho_si <= 0 then return nil end local v_si = 1 / rho_si return { x = v_si * v_scale, y = p_si * p_scale, p = p_si, v = v_si, rho = rho_si, T = T, q = q, t = log10(p_si), } end local function pv_critical_point(constants, v_scale, p_scale) if not constants or not finite_number(constants.pcrit) or not finite_number(constants.Tcrit) or not finite_number(constants.rhocrit) or constants.pcrit <= 0 or constants.Tcrit <= 0 or constants.rhocrit <= 0 then return nil end local v_si = 1 / constants.rhocrit return { x = v_si * v_scale, y = constants.pcrit * p_scale, p = constants.pcrit, v = v_si, rho = constants.rhocrit, T = constants.Tcrit, t = log10(constants.pcrit), critical = true, critical_limit = true, quality_defined = false, } end local function ph_saturation_point_from_T(fluid, T, q, h_scale, p_scale) local p_si = safe_propsSI("P", "T", T, "Q", q, fluid) local h_si = safe_propsSI("H", "T", T, "Q", q, fluid) if not p_si or not h_si or p_si <= 0 then return nil end return {x = h_si * h_scale, y = p_si * p_scale, p = p_si, h = h_si, T = T, q = q, t = log10(p_si)} end local function midpoint_error(a, b, m, opts) local xlin local dx if opts.log_x then if a.x <= 0 or b.x <= 0 or m.x <= 0 then return math.huge end xlin = 0.5 * (log10(a.x) + log10(b.x)) dx = (log10(m.x) - xlin) * (opts.log_x_weight or 30) else xlin = 0.5 * (a.x + b.x) dx = m.x - xlin end local dy if opts.log_y == false then local ylin = 0.5 * (a.y + b.y) dy = (m.y - ylin) * (opts.y_weight or 1) else local ylin = 0.5 * (log10(a.y) + log10(b.y)) dy = (log10(m.y) - ylin) * opts.log_weight end return math.sqrt(dx*dx + dy*dy) end local function adaptive_segment(events, point_func, ta, tb, pa, pb, opts, depth) local tm = 0.5 * (ta + tb) local pm = point_func(tm) events.evaluations = events.evaluations + 1 if not pm then events.omitted = events.omitted + 1 end local refine = depth < opts.max_depth and (not pa or not pb or not pm or midpoint_error(pa, pb, pm, opts) > opts.tolerance) if refine then adaptive_segment(events, point_func, ta, tm, pa, pm, opts, depth + 1) adaptive_segment(events, point_func, tm, tb, pm, pb, opts, depth + 1) return end -- A false event is an explicit topological break. In particular, never use -- an interpolated fake endpoint when CoolProp rejects the midpoint: doing so -- could reconnect two physically disconnected valid domains. events[#events + 1] = pa or false if not (pa and pm and pb) then events[#events + 1] = pm or false end events[#events + 1] = pb or false end local function adaptive_logp_curve(t0, t1, point_func, opts) local intervals = math.max(1, math.floor(opts.initial_intervals or 18)) local events = {evaluations = 0, omitted = 0} for i = 0, intervals - 1 do local ta = t0 + (t1 - t0) * i / intervals local tb = t0 + (t1 - t0) * (i + 1) / intervals local pa, pb = point_func(ta), point_func(tb) events.evaluations = events.evaluations + 2 if not pa then events.omitted = events.omitted + 1 end if not pb then events.omitted = events.omitted + 1 end if i == 0 then events.starts_invalid = not pa end if i == intervals - 1 then events.ends_invalid = not pb end adaptive_segment(events, point_func, ta, tb, pa, pb, opts, 0) end local segments, current = {}, nil for _, event in ipairs(events) do if not event then current = nil else if not current then current = {} segments[#segments + 1] = current end local previous = current[#current] if not previous or math.abs((event.t or 0) - (previous.t or 0)) > 1e-14 * math.max(1, math.abs(event.t or 0)) then current[#current + 1] = event end end end local points, retained = {}, {} for _, segment in ipairs(segments) do if #segment >= 2 then if #retained > 0 then segment[1]._break_before = true end retained[#retained + 1] = segment for _, point in ipairs(segment) do points[#points + 1] = point end end end points.segments = retained points.domain_gap_count = math.max(0, #retained - 1) points.omitted_sample_count = events.omitted points.partially_omitted = events.omitted > 0 and #points > 0 points.starts_invalid = events.starts_invalid points.ends_invalid = events.ends_invalid local policy = M._strict.choice_value(opts.domain_policy or opts.discontinuity_policy, "ignore", "domain policy", {"ignore", "silent", "warning", "warn", "error"}) if events.omitted > 0 and policy ~= "ignore" and policy ~= "silent" then local message = "luacoolprop: " .. tostring(opts.domain_context or "curve") .. " omitted " .. tostring(events.omitted) .. " invalid CoolProp samples and retained " .. tostring(#retained) .. " disconnected valid segment(s)" if policy == "error" then error(message, 3) end if type(texio) == "table" and type(texio.write_nl) == "function" then texio.write_nl("term and log", "warning: " .. message) else io.stderr:write("warning: ", message, "\n") end end return points end -- Refine the high-pressure part of a quality curve in the natural critical -- coordinate epsilon = (pc-p)/pc. A fixed epsilon grid does not honour the -- adaptive tolerance after its points have been merged into the main curve. -- Running the same midpoint test in log10(epsilon) resolves the singular -- curvature without a fluid-specific pressure interval. The pass begins no -- lower than pc/2; ordinary log-pressure sampling remains responsible for the -- smooth, low-pressure part of the curve. local function refine_quality_critical_tail(points, point_func, pmin, pmax, pcrit, sampling) if not finite_number(pmin) or not finite_number(pmax) or not finite_number(pcrit) or pcrit <= 0 then return points end local epsilon_near = 1 - pmax / pcrit local epsilon_far = math.min(0.5, 1 - pmin / pcrit) if epsilon_near <= 0 or epsilon_far <= epsilon_near * (1 + 1e-12) then return points end local tail = adaptive_logp_curve(log10(epsilon_far), log10(epsilon_near), function(logepsilon) local pressure = pcrit * (1 - 10 ^ logepsilon) local roundoff = 2e-14 * pcrit if math.abs(pressure - pmin) <= roundoff then pressure = pmin end if math.abs(pressure - pmax) <= roundoff then pressure = pmax end if pressure < pmin or pressure > pmax then return nil end return point_func(log10(pressure)) end, sampling) for _, point in ipairs(tail) do points[#points + 1] = point end points.omitted_sample_count = (points.omitted_sample_count or 0) + (tail.omitted_sample_count or 0) points.partially_omitted = (points.omitted_sample_count or 0) > 0 and #points > 0 table.sort(points, function(a, b) return a.p < b.p end) local write = 0 for _, point in ipairs(points) do if write == 0 or math.abs(point.p - points[write].p) > 1e-13 * math.max(point.p, points[write].p) then write = write + 1 points[write] = point else points[write] = point end end for index = #points, write + 1, -1 do points[index] = nil end return points end --- Rebuild connected-component metadata after a specialised sampler mutates -- the flat point array (for example by appending an exact critical point). -- @tparam table points Flat curve points using `_break_before` markers. -- @treturn table The same table, with fresh `segments` and gap metadata. function M.diagram.refresh_segments(points) local segments, current = {}, nil for _, point in ipairs(points or {}) do if point._break_before or not current then current = {} segments[#segments + 1] = current end current[#current + 1] = point end points.segments = segments points.domain_gap_count = math.max(0, #segments - 1) return points end --- Append sampled or exact points while preserving domain boundaries. -- This is used by compound curves whose physical branches are evaluated -- independently. An invalid endpoint on either side forces a topological -- break before the next retained component. -- @tparam table target Destination point array. -- @tparam table source Sampled point array or exact point list. -- @treturn table `target`. function M.diagram.append_curve_points(target, source) source = source or {} if #source > 0 then if #target > 0 and (target._ends_invalid or source.starts_invalid) then source[1]._break_before = true end for _, point in ipairs(source) do target[#target + 1] = point end end target.omitted_sample_count = (target.omitted_sample_count or 0) + (source.omitted_sample_count or 0) target._ends_invalid = source.ends_invalid or (#source == 0 and (source.omitted_sample_count or 0) > 0) target.partially_omitted = (target.omitted_sample_count or 0) > 0 and #target > 0 return target end --- Generate a PH quality curve. -- @tparam table opts Curve options. Required fields include `fluid`, -- `pressure_min`, and `pressure_max`; optional scales are `h_scale` and `p_scale`. -- @tparam number q Vapor quality in [0,1]. -- @treturn table Array of points `{x=..., y=..., p=..., h=...}`. When the -- requested interval reaches `pcrit`, the final record is a common limiting -- endpoint with `critical=true`, `critical_limit=true`, and -- `quality_defined=false`. It does not assign a vapor quality at the critical -- state. The high-pressure branch is adaptively refined in logarithmic -- reduced pressure. function M.diagram.ph_isoquality_curve(opts, q) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local requested_pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local requested_pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") q = number(q, "q") if requested_pmin <= 0 or requested_pmax <= requested_pmin then error("luacoolprop: invalid pressure interval for PH quality curve", 2) end if not finite_number(q) or q < 0 or q > 1 then error("luacoolprop: q must be a finite vapor quality between 0 and 1", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local pmin = requested_pmin if constants.ptriple and pmin < constants.ptriple then pmin = constants.ptriple * (1 + 1e-10) end local pmax = requested_pmax local critical_point if constants.pcrit and requested_pmax >= constants.pcrit then pmax = quality_subcritical_pressure(fluid, constants.pcrit) or constants.pcrit * 0.999 critical_point = ph_critical_point(fluid, constants, opts.h_scale or 1e-3, opts.p_scale or 1e-5, pmax) elseif constants.pcrit then pmax = math.min(pmax, constants.pcrit * (1 - 1e-10)) end if not pmax or pmax <= pmin then return {} end local t0, t1 = log10(pmin), log10(pmax) local point_func = function(t) return ph_point_from_logp_quality(fluid, q, t, opts.h_scale or 1e-3, opts.p_scale or 1e-5) end local sampling = { initial_intervals = to_number(opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.25), log_weight = to_number(opts.log_weight, 30), } M._strict.domain_options(opts, sampling, "PH quality curve") local points = adaptive_logp_curve(t0, t1, point_func, sampling) if constants.pcrit then refine_quality_critical_tail(points, point_func, pmin, pmax, constants.pcrit, sampling) end if critical_point and #points > 0 then -- The final subcritical sample captures the steep local curvature. The -- critical limit then closes every quality curve at one bit-identical -- endpoint. Quality itself is undefined at that endpoint. points[#points + 1] = critical_point end return M.diagram.refresh_segments(points) end --- Generate a PH isotherm. -- In the two-phase region the isotherm is an isobar; the generator therefore -- emits only the two saturated endpoints for that plateau. -- @tparam table opts Curve options. -- @tparam number T Temperature in K. -- @treturn table Array of plot points. function M.diagram.ph_isotherm_curve(opts, T) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") T = number(T, "T") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for PH isotherm curve", 2) end if not finite_number(T) or T <= 0 then error("luacoolprop: T must be a finite absolute temperature greater than zero", 2) end local c = opts.constants or M.fluid_constants(fluid, opts.library) local h_scale = opts.h_scale or 1e-3 local p_scale = opts.p_scale or 1e-5 local points = {} local eps = to_number(opts.saturation_pressure_epsilon, 1e-5) local sampling = { initial_intervals = to_number(opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.35), log_weight = to_number(opts.log_weight, 30), } M._strict.domain_options(opts, sampling, "PH isotherm") local has_plateau = c.Tcrit and c.Ttriple and T > c.Ttriple and T < c.Tcrit local psat, sat_l, sat_v if has_plateau then sat_l = ph_saturation_point_from_T(fluid, T, 0, h_scale, p_scale) sat_v = ph_saturation_point_from_T(fluid, T, 1, h_scale, p_scale) if sat_l and sat_v then psat = sat_l.p else has_plateau = false end end local function append_points(src) M.diagram.append_curve_points(points, src) end if has_plateau and psat and psat > pmin and psat < pmax then -- Superheated-vapour branch up to just below saturation. local pv1 = pmin local pv2 = psat * (1 - eps) if pv2 > pv1 then append_points(adaptive_logp_curve(log10(pv1), log10(pv2), function(t) return ph_point_from_logp_temperature(fluid, T, t, h_scale, p_scale) end, sampling)) end -- Two-phase part: exactly the two saturation endpoints. No pointless -- sampling is done inside the horizontal isobaric plateau. append_points({sat_v, sat_l}) -- Compressed/subcooled-liquid branch above saturation. local pl1 = psat * (1 + eps) local pl2 = pmax if pl2 > pl1 then append_points(adaptive_logp_curve(log10(pl1), log10(pl2), function(t) return ph_point_from_logp_temperature(fluid, T, t, h_scale, p_scale) end, sampling)) end else append_points(adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return ph_point_from_logp_temperature(fluid, T, t, h_scale, p_scale) end, sampling)) end points._ends_invalid = nil return M.diagram.refresh_segments(points) end --- Generate a PH isentrope. -- @tparam table opts Curve options. -- @tparam number s Specific entropy in J/(kg K). -- @treturn table Array of plot points. function M.diagram.ph_isentrope_curve(opts, s) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") s = number(s, "s") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for PH isentrope curve", 2) end if not finite_number(s) then error("luacoolprop: s must be a finite specific entropy", 2) end local point_func = function(t) return ph_point_from_logp_entropy(fluid, s, t, opts.h_scale or 1e-3, opts.p_scale or 1e-5) end local gen = { initial_intervals = to_number(opts.isentrope_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isentrope_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isentrope_tolerance or opts.tolerance, 0.30), log_weight = to_number(opts.log_weight, 30), } M._strict.domain_options(opts, gen, "PH isentrope") return adaptive_logp_curve(log10(pmin), log10(pmax), point_func, gen) end --- Generate a PH isochore. -- CoolProp is called with density, `D=1/v`, because specific volume is not a -- native input key in the high-level API. -- @tparam table opts Curve options. -- @tparam number v Specific volume in m^3/kg. -- @treturn table Array of plot points. function M.diagram.ph_isochore_curve(opts, v) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") v = number(v, "v") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for PH isochore curve", 2) end if not finite_number(v) or v <= 0 then error("luacoolprop: v must be a finite specific volume greater than zero", 2) end local point_func = function(t) return ph_point_from_logp_specific_volume(fluid, v, t, opts.h_scale or 1e-3, opts.p_scale or 1e-5) end local gen = { initial_intervals = to_number(opts.isochore_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isochore_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isochore_tolerance or opts.tolerance, 0.28), log_weight = to_number(opts.log_weight, 30), } M._strict.domain_options(opts, gen, "PH isochore") return adaptive_logp_curve(log10(pmin), log10(pmax), point_func, gen) end --- Generate a PV quality curve. -- Both plot coordinates are positive and are intended for logarithmic axes. -- The endpoint at the critical pressure is computed once from the canonical -- critical density, so all qualities meet at exactly the same plotted point. -- @tparam table opts Curve options, including `fluid`, `pressure_min`, -- `pressure_max`, `v_scale`, and `p_scale`. -- @tparam number q Vapour quality in [0,1]. -- @treturn table Points `{x=v*v_scale,y=p*p_scale,p=...,v=...,rho=...}`; -- the high-pressure branch is adaptively refined in logarithmic reduced -- pressure. function M.diagram.pv_quality_curve(opts, q) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local requested_pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local requested_pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") q = number(q, "q") if requested_pmin <= 0 or requested_pmax <= requested_pmin then error("luacoolprop: invalid pressure interval for PV quality curve", 2) end if not finite_number(q) or q < 0 or q > 1 then error("luacoolprop: q must be a finite vapour quality between 0 and 1", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) local pmin = requested_pmin if constants.ptriple and pmin < constants.ptriple then pmin = constants.ptriple * (1 + 1e-10) end local pmax = requested_pmax local critical_point if constants.pcrit and requested_pmax >= constants.pcrit then pmax = quality_subcritical_pressure(fluid, constants.pcrit) or constants.pcrit * 0.999 critical_point = pv_critical_point(constants, v_scale, p_scale) elseif constants.pcrit then pmax = math.min(pmax, constants.pcrit * (1 - 1e-10)) end if pmax <= pmin then return {} end local point_func = function(t) return pv_point_from_logp_pair(fluid, "Q", q, t, v_scale, p_scale) end local sampling = { initial_intervals = to_number( opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x = true, log_x_weight = to_number(opts.log_x_weight, 30), } M._strict.domain_options(opts, sampling, "PV quality curve") local points = adaptive_logp_curve(log10(pmin), log10(pmax), point_func, sampling) if constants.pcrit then refine_quality_critical_tail(points, point_func, pmin, pmax, constants.pcrit, sampling) end if critical_point then if #points > 0 then points[#points + 1] = critical_point end end return M.diagram.refresh_segments(points) end --- Generate a PV isotherm. -- Below the critical temperature the two saturated endpoints are joined by -- the physical constant-pressure two-phase segment. -- @tparam table opts Curve options. -- @tparam number T Absolute temperature in K. -- @treturn table Array of PV plot points. function M.diagram.pv_isotherm_curve(opts, T) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") T = number(T, "T") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for PV isotherm curve", 2) end if not finite_number(T) or T <= 0 then error("luacoolprop: T must be a finite absolute temperature greater than zero", 2) end local c = opts.constants or M.fluid_constants(fluid, opts.library) local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) local sampling = { initial_intervals = to_number( opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x = true, log_x_weight = to_number(opts.log_x_weight, 30), } M._strict.domain_options(opts, sampling, "PV isotherm") local points = {} local function append_points(source) M.diagram.append_curve_points(points, source) end local epsilon = to_number(opts.saturation_pressure_epsilon, 1e-5) local has_plateau = c.Tcrit and c.Ttriple and T > c.Ttriple and T < c.Tcrit local psat, sat_l, sat_v if has_plateau then sat_l = pv_saturation_point_from_T(fluid, T, 0, v_scale, p_scale) sat_v = pv_saturation_point_from_T(fluid, T, 1, v_scale, p_scale) if sat_l and sat_v then psat = sat_l.p else has_plateau = false end end if has_plateau and psat and psat > pmin and psat < pmax then local vapour_max = psat * (1 - epsilon) if vapour_max > pmin then append_points(adaptive_logp_curve(log10(pmin), log10(vapour_max), function(t) return pv_point_from_logp_pair(fluid, "T", T, t, v_scale, p_scale) end, sampling)) end append_points({sat_v, sat_l}) local liquid_min = psat * (1 + epsilon) if pmax > liquid_min then append_points(adaptive_logp_curve(log10(liquid_min), log10(pmax), function(t) return pv_point_from_logp_pair(fluid, "T", T, t, v_scale, p_scale) end, sampling)) end else append_points(adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return pv_point_from_logp_pair(fluid, "T", T, t, v_scale, p_scale) end, sampling)) end points._ends_invalid = nil return M.diagram.refresh_segments(points) end --- Generate a PV isentrope. -- @tparam table opts Curve options. -- @tparam number s Specific entropy in J/(kg K). -- @treturn table Array of PV plot points. function M.diagram.pv_isentrope_curve(opts, s) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") s = number(s, "s") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for PV isentrope curve", 2) end if not finite_number(s) then error("luacoolprop: s must be a finite specific entropy", 2) end local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) local point_func = function(t) return pv_point_from_logp_pair(fluid, "S", s, t, v_scale, p_scale) end return adaptive_logp_curve(log10(pmin), log10(pmax), point_func, M._strict.domain_options(opts, { initial_intervals = to_number( opts.isentrope_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isentrope_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isentrope_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x = true, log_x_weight = to_number(opts.log_x_weight, 30), }, "PV isentrope")) end -- Generate a PV isenthalp for process paths. It is intentionally private: -- only quality, isotherm, and isentrope are public background families. local function pv_isenthalp_curve(opts, h) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") if pmin <= 0 or pmax <= pmin or not finite_number(h) then return {} end local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) return adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return pv_point_from_logp_pair(fluid, "H", h, t, v_scale, p_scale) end, M._strict.domain_options(opts, { initial_intervals = to_number(opts.initial_intervals, 14), max_depth = to_number(opts.max_depth, 7), tolerance = to_number(opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x = true, log_x_weight = to_number(opts.log_x_weight, 30), }, "PV isenthalpic process")) end -- Build a complete TS point from a CoolProp input pair. The plotted -- coordinates are deliberately kept separate from the raw SI properties: -- x=s*s_scale and y=T*temperature_scale. Keeping p, h, s, T, v and q on -- every valid point makes the curve generators useful to texlua clients as -- well as to the TeX renderer. local function ts_point_from_pair(fluid, name1, value1, name2, value2, parameter, s_scale, temperature_scale, temperature_offset) local entropy = safe_propsSI("S", name1, value1, name2, value2, fluid) local temperature = safe_propsSI("T", name1, value1, name2, value2, fluid) if not entropy or not temperature or temperature <= 0 then return nil end local pressure = safe_propsSI("P", name1, value1, name2, value2, fluid) local enthalpy = safe_propsSI("H", name1, value1, name2, value2, fluid) local density = safe_propsSI("Dmass", name1, value1, name2, value2, fluid) local quality = safe_propsSI("Q", name1, value1, name2, value2, fluid) if quality and (quality < 0 or quality > 1) then quality = nil end return { x = entropy * s_scale, y = (temperature - (temperature_offset or 0)) * temperature_scale, p = pressure, h = enthalpy, s = entropy, T = temperature, rho = density, v = density and density > 0 and 1 / density or nil, q = quality, t = parameter, } end local function ts_point_from_logp_pair(fluid, input, value, t, s_scale, temperature_scale, temperature_offset) local pressure = 10 ^ t local point = ts_point_from_pair(fluid, "P", pressure, input, value, t, s_scale, temperature_scale, temperature_offset) if point then point.p = pressure end if point and input == "Q" then point.q = value end if point and input == "H" then point.h = value end return point end -- P,Q is singular at the critical state. As for PH and PV, use the unique -- (Tcrit,rhocrit) state so all quality curves terminate at one bit-identical -- TS coordinate. The mean of the two closest saturation entropies is a -- backend-independent fallback when the exact pair is rejected. local function ts_critical_point(fluid, constants, s_scale, temperature_scale, temperature_offset, fallback_pressure) if not constants or not finite_number(constants.pcrit) or not finite_number(constants.Tcrit) or not finite_number(constants.rhocrit) or constants.pcrit <= 0 or constants.Tcrit <= 0 or constants.rhocrit <= 0 then return nil end local entropy = safe_propsSI("S", "T", constants.Tcrit, "Dmass", constants.rhocrit, fluid) local approximated = false if not entropy and fallback_pressure then local liquid = safe_propsSI("S", "P", fallback_pressure, "Q", 0, fluid) local vapour = safe_propsSI("S", "P", fallback_pressure, "Q", 1, fluid) if liquid and vapour then entropy = 0.5 * (liquid + vapour) approximated = true end end if not entropy then return nil end return { x = entropy * s_scale, y = (constants.Tcrit - (temperature_offset or 0)) * temperature_scale, p = constants.pcrit, s = entropy, T = constants.Tcrit, rho = constants.rhocrit, v = 1 / constants.rhocrit, t = log10(constants.pcrit), critical = true, critical_limit = true, quality_defined = false, critical_approximated = approximated, } end --- Generate a TS constant-quality curve. -- Entropy is plotted horizontally and absolute temperature vertically. Raw -- thermodynamic values remain in SI units on each returned point. -- @tparam table opts Options including pressure bounds, `s_scale`, and -- `temperature_scale`. -- @tparam number quality Vapour quality in [0,1]. -- @treturn table TS points `{x,y,p,h,s,T,v,rho,q}`. Sampling close to the -- critical state is adaptive in logarithmic reduced pressure and uses the TS -- display-coordinate error metric. function M.diagram.ts_quality_curve(opts, quality) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local requested_pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local requested_pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") quality = number(quality, "quality") if requested_pmin <= 0 or requested_pmax <= requested_pmin then error("luacoolprop: invalid pressure interval for TS quality curve", 2) end if not finite_number(quality) or quality < 0 or quality > 1 then error("luacoolprop: quality must be a finite number between 0 and 1", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local s_scale = to_number(opts.s_scale, 1e-3) local temperature_scale = to_number(opts.temperature_scale or opts.t_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local pmin = requested_pmin if constants.ptriple and pmin < constants.ptriple then pmin = constants.ptriple * (1 + 1e-10) end local pmax = requested_pmax local critical if constants.pcrit and requested_pmax >= constants.pcrit then pmax = quality_subcritical_pressure(fluid, constants.pcrit) or constants.pcrit * 0.999 critical = ts_critical_point(fluid, constants, s_scale, temperature_scale, temperature_offset, pmax) elseif constants.pcrit then pmax = math.min(pmax, constants.pcrit * (1 - 1e-10)) end if pmax <= pmin then return {} end local sampling = { initial_intervals = to_number( opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 8), log_y = false, y_weight = to_number(opts.temperature_weight, 0.01), } M._strict.domain_options(opts, sampling, "TS quality curve") local point_func = function(t) return ts_point_from_logp_pair(fluid, "Q", quality, t, s_scale, temperature_scale, temperature_offset) end local points = adaptive_logp_curve(log10(pmin), log10(pmax), point_func, sampling) if constants.pcrit then refine_quality_critical_tail(points, point_func, pmin, pmax, constants.pcrit, sampling) end if critical then if #points > 0 then points[#points + 1] = critical end end return M.diagram.refresh_segments(points) end --- Generate a TS isenthalp. -- @tparam table opts Options including pressure bounds and coordinate scales. -- @tparam number enthalpy Specific enthalpy in J/kg. -- @treturn table TS points. function M.diagram.ts_isenthalp_curve(opts, enthalpy) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") enthalpy = number(enthalpy, "enthalpy") if pmin <= 0 or pmax <= pmin then error("luacoolprop: invalid pressure interval for TS isenthalp", 2) end if not finite_number(enthalpy) then error("luacoolprop: enthalpy must be a finite specific enthalpy", 2) end local s_scale = to_number(opts.s_scale, 1e-3) local temperature_scale = to_number(opts.temperature_scale or opts.t_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) return adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return ts_point_from_logp_pair(fluid, "H", enthalpy, t, s_scale, temperature_scale, temperature_offset) end, M._strict.domain_options(opts, { initial_intervals = to_number( opts.isenthalp_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isenthalp_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isenthalp_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 8), log_y = false, y_weight = to_number(opts.temperature_weight, 0.01), }, "TS isenthalp")) end -- Complete one HS point from any CoolProp input pair. Plot coordinates are -- scaled at this boundary while the thermodynamic record retains SI values. local function hs_point_from_pair(fluid, name1, value1, name2, value2, parameter, s_scale, h_scale) local entropy = safe_propsSI("S", name1, value1, name2, value2, fluid) local enthalpy = safe_propsSI("H", name1, value1, name2, value2, fluid) if not entropy or not enthalpy then return nil end local pressure = safe_propsSI("P", name1, value1, name2, value2, fluid) local temperature = safe_propsSI("T", name1, value1, name2, value2, fluid) local density = safe_propsSI("Dmass", name1, value1, name2, value2, fluid) local quality = safe_propsSI("Q", name1, value1, name2, value2, fluid) if quality and (quality < 0 or quality > 1) then quality = nil end return { x = entropy * s_scale, y = enthalpy * h_scale, p = pressure, h = enthalpy, s = entropy, T = temperature, rho = density, v = density and density > 0 and 1 / density or nil, q = quality, t = parameter, } end local function hs_point_from_logp_pair(fluid, input, value, t, s_scale, h_scale) local pressure = 10 ^ t local point = hs_point_from_pair(fluid, "P", pressure, input, value, t, s_scale, h_scale) if point then point.p = pressure end if point and input == "Q" then point.q = value end if point and input == "T" then point.T = value end return point end local function hs_critical_point(fluid, constants, s_scale, h_scale, fallback_pressure) if not constants or not finite_number(constants.pcrit) or not finite_number(constants.Tcrit) or not finite_number(constants.rhocrit) or constants.pcrit <= 0 or constants.Tcrit <= 0 or constants.rhocrit <= 0 then return nil end local entropy = safe_propsSI("S", "T", constants.Tcrit, "Dmass", constants.rhocrit, fluid) local enthalpy = safe_propsSI("H", "T", constants.Tcrit, "Dmass", constants.rhocrit, fluid) local approximated = false if (not entropy or not enthalpy) and fallback_pressure then local sl = safe_propsSI("S", "P", fallback_pressure, "Q", 0, fluid) local sv = safe_propsSI("S", "P", fallback_pressure, "Q", 1, fluid) local hl = safe_propsSI("H", "P", fallback_pressure, "Q", 0, fluid) local hv = safe_propsSI("H", "P", fallback_pressure, "Q", 1, fluid) if sl and sv and hl and hv then entropy = entropy or 0.5 * (sl + sv) enthalpy = enthalpy or 0.5 * (hl + hv) approximated = true end end if not entropy or not enthalpy then return nil end return { x = entropy * s_scale, y = enthalpy * h_scale, p = constants.pcrit, h = enthalpy, s = entropy, T = constants.Tcrit, rho = constants.rhocrit, v = 1 / constants.rhocrit, t = log10(constants.pcrit), critical = true, critical_limit = true, quality_defined = false, critical_approximated = approximated, } end --- Generate an HS constant-quality curve. -- @tparam table opts Pressure bounds and coordinate scales. -- @tparam number quality Vapour quality in [0,1]. -- @treturn table HS points with SI thermodynamic fields. Sampling close to -- the critical state is adaptive in logarithmic reduced pressure and uses the -- HS display-coordinate error metric. function M.diagram.hs_quality_curve(opts, quality) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local requested_pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local requested_pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") quality = number(quality, "quality") if requested_pmin <= 0 or requested_pmax <= requested_pmin then error("luacoolprop: invalid pressure interval for HS quality curve", 2) end if not finite_number(quality) or quality < 0 or quality > 1 then error("luacoolprop: quality must be a finite number between 0 and 1", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local s_scale = to_number(opts.s_scale, 1e-3) local h_scale = to_number(opts.h_scale, 1e-3) local pmin = requested_pmin if constants.ptriple and pmin < constants.ptriple then pmin = constants.ptriple * (1 + 1e-10) end local pmax = requested_pmax local critical if constants.pcrit and requested_pmax >= constants.pcrit then pmax = quality_subcritical_pressure(fluid, constants.pcrit) or constants.pcrit * 0.999 critical = hs_critical_point(fluid, constants, s_scale, h_scale, pmax) elseif constants.pcrit then pmax = math.min(pmax, constants.pcrit * (1 - 1e-10)) end if pmax <= pmin then return {} end local sampling = { initial_intervals = to_number( opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.enthalpy_weight, 0.01), log_y = false, y_weight = to_number(opts.enthalpy_weight, 0.01), } M._strict.domain_options(opts, sampling, "HS quality curve") local point_func = function(t) return hs_point_from_logp_pair(fluid, "Q", quality, t, s_scale, h_scale) end local points = adaptive_logp_curve(log10(pmin), log10(pmax), point_func, sampling) if constants.pcrit then refine_quality_critical_tail(points, point_func, pmin, pmax, constants.pcrit, sampling) end if critical then if #points > 0 then points[#points + 1] = critical end end return M.diagram.refresh_segments(points) end --- Generate an HS isochore. -- @tparam table opts Pressure bounds and coordinate scales. -- @tparam number volume Mass-specific volume in m^3/kg. -- @treturn table HS points. function M.diagram.hs_isochore_curve(opts, volume) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") volume = number(volume, "volume") if pmin <= 0 or pmax <= pmin or volume <= 0 then error("luacoolprop: invalid HS isochore arguments", 2) end local s_scale = to_number(opts.s_scale, 1e-3) local h_scale = to_number(opts.h_scale, 1e-3) return adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return hs_point_from_logp_pair(fluid, "Dmass", 1 / volume, t, s_scale, h_scale) end, M._strict.domain_options(opts, { initial_intervals = to_number( opts.isochore_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isochore_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isochore_tolerance or opts.tolerance, 0.10), log_weight = 1, log_y = false, y_weight = to_number(opts.enthalpy_weight, 0.01), }, "HS isochore")) end --- Generate an HS isotherm, including its exact two-phase segment. -- @tparam table opts Pressure bounds and coordinate scales. -- @tparam number temperature Absolute temperature in K. -- @treturn table HS points. function M.diagram.hs_isotherm_curve(opts, temperature) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") temperature = number(temperature, "temperature") if pmin <= 0 or pmax <= pmin or temperature <= 0 then error("luacoolprop: invalid HS isotherm arguments", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local s_scale = to_number(opts.s_scale, 1e-3) local h_scale = to_number(opts.h_scale, 1e-3) local sampling = { initial_intervals = to_number( opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.10), log_weight = 1, log_y = false, y_weight = to_number(opts.enthalpy_weight, 0.01), } M._strict.domain_options(opts, sampling, "HS isotherm") local points = {} local function append(source) M.diagram.append_curve_points(points, source) end local has_plateau = constants.Ttriple and constants.Tcrit and temperature > constants.Ttriple and temperature < constants.Tcrit local liquid, vapour, psat if has_plateau then liquid = hs_point_from_pair(fluid, "T", temperature, "Q", 0, 0, s_scale, h_scale) vapour = hs_point_from_pair(fluid, "T", temperature, "Q", 1, 0, s_scale, h_scale) psat = liquid and liquid.p or nil if not liquid or not vapour or not psat then has_plateau = false end end local epsilon = to_number(opts.saturation_pressure_epsilon, 1e-5) if has_plateau and psat > pmin and psat < pmax then local vapour_max = psat * (1 - epsilon) if vapour_max > pmin then append(adaptive_logp_curve(log10(pmin), log10(vapour_max), function(t) return hs_point_from_logp_pair(fluid, "T", temperature, t, s_scale, h_scale) end, sampling)) end vapour.t, liquid.t = log10(psat), log10(psat) append({vapour, liquid}) local liquid_min = psat * (1 + epsilon) if pmax > liquid_min then append(adaptive_logp_curve(log10(liquid_min), log10(pmax), function(t) return hs_point_from_logp_pair(fluid, "T", temperature, t, s_scale, h_scale) end, sampling)) end else append(adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return hs_point_from_logp_pair(fluid, "T", temperature, t, s_scale, h_scale) end, sampling)) end points._ends_invalid = nil return M.diagram.refresh_segments(points) end --- Generate an HS isobar over a temperature interval. -- Entropy is used as the actual sampling parameter, which crosses the -- two-phase region continuously and includes its exact constant-pressure -- mixture segment. -- @tparam table opts Options including temperature bounds and scales. -- @tparam number pressure Constant pressure in Pa. -- @treturn table HS points. function M.diagram.hs_isobar_curve(opts, pressure) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" pressure = number(pressure, "pressure") if pressure <= 0 then error("luacoolprop: HS isobar pressure must be positive", 2) end local constants = opts.constants or M.fluid_constants(fluid, opts.library) local temperature_unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) local temperature_offset = (temperature_unit == "celsius" or temperature_unit == "c" or temperature_unit == "degc") and 273.15 or 0 local Tmin_user = to_number(opts.isobar_temperature_min, nil) local Tmax_user = to_number(opts.isobar_temperature_max, nil) local Tmin = Tmin_user and (Tmin_user + temperature_offset) or math.max((constants.Tmin or constants.Ttriple or 200) + 1, (constants.Ttriple or 200) + 1) local Tmax = Tmax_user and (Tmax_user + temperature_offset) or math.min((constants.Tmax or 600) - 1, (constants.Tcrit or 400) + 150) if Tmax <= Tmin then if Tmin_user or Tmax_user then error("luacoolprop: isobar_temperature_max must be greater than isobar_temperature_min", 2) end return {} end local s1 = safe_propsSI("S", "P", pressure, "T", Tmin, fluid) local s2 = safe_propsSI("S", "P", pressure, "T", Tmax, fluid) if not s1 or not s2 or s1 == s2 then return {} end local s_scale = to_number(opts.s_scale, 1e-3) local h_scale = to_number(opts.h_scale, 1e-3) return adaptive_logp_curve(s1, s2, function(entropy) return hs_point_from_pair(fluid, "P", pressure, "S", entropy, entropy, s_scale, h_scale) end, M._strict.domain_options(opts, { initial_intervals = to_number( opts.isobar_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isobar_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isobar_tolerance or opts.tolerance, 0.10), log_weight = 1, log_y = false, y_weight = to_number(opts.enthalpy_weight, 0.01), }, "HS isobar")) end -- Complete one PT point from pressure and a second independent property. -- Temperature and pressure are the plotted coordinates; all other fields are -- retained in SI units when CoolProp can resolve them. PT itself is not used -- as an input pair on the saturation line because that pair is singular for a -- pure fluid there. function M.diagram._pt_point_from_logp_pair(fluid, input, value, t, temperature_scale, temperature_offset, p_scale) local pressure = 10 ^ t local temperature = safe_propsSI("T", "P", pressure, input, value, fluid) if not temperature or temperature <= 0 then return nil end local point = { x = (temperature - (temperature_offset or 0)) * temperature_scale, y = pressure * p_scale, p = pressure, T = temperature, t = t, } point.h = safe_propsSI("H", "P", pressure, input, value, fluid) point.s = safe_propsSI("S", "P", pressure, input, value, fluid) point.rho = safe_propsSI("Dmass", "P", pressure, input, value, fluid) point.v = point.rho and point.rho > 0 and 1 / point.rho or nil local quality = safe_propsSI("Q", "P", pressure, input, value, fluid) if quality and quality >= 0 and quality <= 1 then point.q = quality end if input == "S" then point.s = value elseif input == "H" then point.h = value elseif input == "Dmass" then point.rho = value point.v = value > 0 and 1 / value or nil elseif input == "Q" then point.q = value end return point end -- Evaluate one point of the pure-fluid liquid--vapour coexistence curve. -- Both limiting qualities have the same PT coordinates. Evaluating both is -- nevertheless useful: it prevents a backend-specific failure of one branch -- from silently producing an incomplete phase envelope. function M.diagram._pt_saturation_point_from_logp(fluid, t, temperature_scale, temperature_offset, p_scale) local pressure = 10 ^ t local liquid_temperature = safe_propsSI("T", "P", pressure, "Q", 0, fluid) local vapour_temperature = safe_propsSI("T", "P", pressure, "Q", 1, fluid) if not liquid_temperature or not vapour_temperature then return nil end local temperature = 0.5 * (liquid_temperature + vapour_temperature) return { x = (temperature - (temperature_offset or 0)) * temperature_scale, y = pressure * p_scale, p = pressure, T = temperature, t = t, phase_equilibrium = true, quality_defined = false, } end function M.diagram._pt_endpoint(constants, which, temperature_scale, temperature_offset, p_scale) local pressure = which == "triple" and constants.ptriple or constants.pcrit local temperature = which == "triple" and constants.Ttriple or constants.Tcrit if not finite_number(pressure) or not finite_number(temperature) or pressure <= 0 or temperature <= 0 then return nil end return { x = (temperature - (temperature_offset or 0)) * temperature_scale, y = pressure * p_scale, p = pressure, T = temperature, t = log10(pressure), phase_equilibrium = true, quality_defined = false, triple = which == "triple" or nil, critical = which == "critical" or nil, critical_limit = which == "critical" or nil, } end function M.diagram._pt_subcritical_pressure(fluid, critical_pressure) for _, epsilon in ipairs(quality_critical_epsilons) do local pressure = critical_pressure * (1 - epsilon) local liquid = safe_propsSI("T", "P", pressure, "Q", 0, fluid) local vapour = safe_propsSI("T", "P", pressure, "Q", 1, fluid) if liquid and vapour then return pressure end end return nil end -- Find the closest usable saturation pressure above the triple point. Some -- CoolProp backends expose a minute mismatch between the tabulated Ttriple -- constant and the saturation solver exactly at its lower validity boundary. -- Starting the numerical branch only once its temperature reaches Ttriple -- keeps the exact endpoint and the sampled coexistence locus monotone without -- embedding a fluid-specific pressure offset. function M.diagram._pt_supertriple_pressure(fluid, constants) local ptriple, Ttriple = constants.ptriple, constants.Ttriple local pcrit = constants.pcrit local tolerance = math.max(1e-9, math.abs(Ttriple) * 1e-12) local epsilon_lower = 0 for _, epsilon_upper in ipairs({ 1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 0.5, }) do local pressure = ptriple * (1 + epsilon_upper) if pressure >= pcrit then break end local liquid = safe_propsSI("T", "P", pressure, "Q", 0, fluid) local vapour = safe_propsSI("T", "P", pressure, "Q", 1, fluid) if liquid and vapour and 0.5 * (liquid + vapour) >= Ttriple - tolerance then -- Bisect in the dimensionless relative-pressure offset. The resulting -- gap is the smallest one justified by the backend, not a plotting -- heuristic or a fluid table. for _ = 1, 32 do local epsilon_mid = 0.5 * (epsilon_lower + epsilon_upper) local midpoint_pressure = ptriple * (1 + epsilon_mid) local midpoint_liquid = safe_propsSI( "T", "P", midpoint_pressure, "Q", 0, fluid) local midpoint_vapour = safe_propsSI( "T", "P", midpoint_pressure, "Q", 1, fluid) if midpoint_liquid and midpoint_vapour and 0.5 * (midpoint_liquid + midpoint_vapour) >= Ttriple - tolerance then epsilon_upper = epsilon_mid else epsilon_lower = epsilon_mid end end return ptriple * (1 + epsilon_upper) end epsilon_lower = epsilon_upper end return nil end function M.diagram._resolve_pt_pressure_limits(opts, fluid) M.diagram.prepare_fluid(opts) local constants = M.fluid_constants(fluid, opts.library) local ptriple = constants.ptriple local pcrit = constants.pcrit if not finite_number(ptriple) or not finite_number(pcrit) or ptriple <= 0 or pcrit <= ptriple then error("luacoolprop: CoolProp did not provide a valid triple-to-critical " .. "pressure interval for pure fluid '" .. tostring(fluid) .. "'", 3) end local pmin = to_number(opts.pressure_min, ptriple) local factor = to_number(opts.pressure_max_factor, 1.10) if factor <= 0 then error("luacoolprop: pressure_max_factor must be strictly positive", 3) end local pmax = to_number(opts.pressure_max, pcrit * factor) if pmin <= 0 or pmax <= pmin then error("luacoolprop: pressure_max must be greater than a strictly positive pressure_min", 3) end return pmin, pmax, constants end --- Generate the pure-fluid liquid--vapour equilibrium curve in a PT diagram. -- The curve is sampled between the triple and critical points, clipped only -- by explicit pressure bounds. Its endpoints come directly from CoolProp's -- fluid constants and are therefore independent of vapour-quality branch. -- @tparam table opts PT options and coordinate scales. -- @treturn table PT points with SI `T` and `p` fields. function M.diagram.pt_phase_envelope_curve(opts) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local requested_min = M._strict.require_decimal(opts.pressure_min, "pressure_min") local requested_max = M._strict.require_decimal(opts.pressure_max, "pressure_max") local constants = opts.constants or M.fluid_constants(fluid, opts.library) local pmin = math.max(requested_min, constants.ptriple) local pmax = math.min(requested_max, constants.pcrit) if pmax <= pmin then return {} end local temperature_scale = to_number(opts.temperature_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local p_scale = to_number(opts.p_scale, 1e-5) local subcritical = pmax local append_critical = requested_max >= constants.pcrit if append_critical then subcritical = M.diagram._pt_subcritical_pressure(fluid, constants.pcrit) or constants.pcrit * (1 - 1e-7) end local sampling = M._strict.domain_options(opts, { initial_intervals = to_number( opts.phase_envelope_initial_intervals or opts.initial_intervals, 20), max_depth = to_number(opts.phase_envelope_max_depth or opts.max_depth, 9), tolerance = to_number( opts.phase_envelope_tolerance or opts.tolerance, 0.08), log_weight = to_number(opts.log_weight, 30), log_y = true, }, "PT liquid-vapour phase envelope") local point_function = function(t) return M.diagram._pt_saturation_point_from_logp(fluid, t, temperature_scale, temperature_offset, p_scale) end local sampling_min = math.max(pmin, constants.ptriple * (1 + 1e-10)) if requested_min <= constants.ptriple then sampling_min = math.max(sampling_min, M.diagram._pt_supertriple_pressure(fluid, constants) or sampling_min) end local points = {} if subcritical > sampling_min then points = adaptive_logp_curve(log10(sampling_min), log10(subcritical), point_function, sampling) else points.segments = {} points.domain_gap_count = 0 points.omitted_sample_count = 0 points.partially_omitted = false end if constants.pcrit and subcritical > sampling_min then refine_quality_critical_tail(points, point_function, pmin, subcritical, constants.pcrit, sampling) end if requested_min <= constants.ptriple then local triple = M.diagram._pt_endpoint(constants, "triple", temperature_scale, temperature_offset, p_scale) if triple then table.insert(points, 1, triple) end end if append_critical then local critical = M.diagram._pt_endpoint(constants, "critical", temperature_scale, temperature_offset, p_scale) if critical then points[#points + 1] = critical end end return M.diagram.refresh_segments(points) end function M.diagram._pt_constant_property_curve(opts, input, value, family) opts = opts or {} M.diagram.prepare_fluid(opts) local fluid = opts.fluid or "R134a" local pmin = M._strict.require_decimal(opts.pressure_min, "pressure_min") local pmax = M._strict.require_decimal(opts.pressure_max, "pressure_max") value = number(value, family) if pmin <= 0 or pmax <= pmin or not finite_number(value) or (input == "Dmass" and value <= 0) then error("luacoolprop: invalid PT " .. family .. " arguments", 3) end local temperature_scale = to_number(opts.temperature_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local p_scale = to_number(opts.p_scale, 1e-5) return adaptive_logp_curve(log10(pmin), log10(pmax), function(t) return M.diagram._pt_point_from_logp_pair(fluid, input, value, t, temperature_scale, temperature_offset, p_scale) end, M._strict.domain_options(opts, { initial_intervals = to_number(opts[family .. "_initial_intervals"] or opts.initial_intervals, 18), max_depth = to_number(opts[family .. "_max_depth"] or opts.max_depth, 8), tolerance = to_number(opts[family .. "_tolerance"] or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 30), log_y = true, }, "PT " .. family)) end --- Generate one PT isentrope (`entropy` in J/(kg K)). function M.diagram.pt_isentrope_curve(opts, entropy) return M.diagram._pt_constant_property_curve(opts, "S", entropy, "isentrope") end --- Generate one PT isenthalp (`enthalpy` in J/kg). function M.diagram.pt_isenthalp_curve(opts, enthalpy) return M.diagram._pt_constant_property_curve(opts, "H", enthalpy, "isenthalp") end --- Generate one PT isochore (`volume` in m^3/kg). function M.diagram.pt_isochore_curve(opts, volume) volume = number(volume, "specific volume") if volume <= 0 then error("luacoolprop: PT isochore volume must be strictly positive", 2) end return M.diagram._pt_constant_property_curve(opts, "Dmass", 1 / volume, "isochore") end local function add_label(texout, id, style, point, text, coord_digits) if point then local name_style = id and ("name=" .. id .. "-label") or nil texout[#texout + 1] = "\\node[" .. join_style(name_style, style) .. "] at (axis cs:" .. fmt_num(point.x, coord_digits) .. "," .. fmt_num(point.y, coord_digits) .. ") {" .. text .. "}; " end end local function clamp01(x) if x < 0 then return 0 end if x > 1 then return 1 end return x end local function finite_plot_point(pt) return pt and finite_number(pt.x) and finite_number(pt.y) end --- Build an arc-length parametrisation of a thermodynamic curve. -- -- A scalar label position must not simply be a point index: adaptive sampling -- produces non-uniform point spacing. This helper builds a polyline measured in -- display-like coordinates, respecting the logarithmic modes of the diagram. -- It is used for manual labels, -- process labels, and default preferred positions forwarded to autonode. -- @tparam table points Curve points in PGFPlots coordinates. -- @tparam[opt] table bounds Axis-like bounds used for normalisation. -- @treturn table Arc-length path samples. local function make_label_path(points, bounds, log_x, log_y) local samples = {} local candidates, candidate = {}, {} for _, point in ipairs(points or {}) do if point._break_before and #candidate > 0 then candidates[#candidates + 1], candidate = candidate, {} end candidate[#candidate + 1] = point end if #candidate > 0 then candidates[#candidates + 1] = candidate end local source = {} for _, segment in ipairs(candidates) do if #segment > #source then source = segment end end local xspan, xmin_metric, yminlog, yspan log_x = to_bool(log_x or (bounds and bounds.log_x), false) log_y = to_bool(log_y == nil and (not bounds or bounds.log_y ~= false) or log_y or (bounds and bounds.log_y), true) if bounds then xmin_metric = log_x and log10(math.max(bounds.xmin, 1e-99)) or bounds.xmin local xmax_metric = log_x and log10(math.max(bounds.xmax, 1e-99)) or bounds.xmax xspan = math.max(1e-12, xmax_metric - xmin_metric) yminlog = log_y and log10(math.max(bounds.ymin, 1e-99)) or bounds.ymin local ymaxlog = log_y and log10(math.max(bounds.ymax, bounds.ymin * 1.0001)) or bounds.ymax yspan = math.max(1e-12, ymaxlog - yminlog) end for _, pt in ipairs(source) do if finite_plot_point(pt) then local lx = log_x and log10(pt.x) or pt.x local ly = log_y and log10(math.max(pt.y, 1e-99)) or pt.y local mx, my = lx, ly if bounds then mx = (lx - xmin_metric) / xspan my = (ly - yminlog) / yspan end samples[#samples + 1] = {x = pt.x, y = pt.y, lx = lx, ly = ly, mx = mx, my = my} end end if #samples == 0 then return samples end samples.log_x = log_x samples.log_y = log_y samples[1].s = 0 local total = 0 for i = 2, #samples do local dx = samples[i].mx - samples[i - 1].mx local dy = samples[i].my - samples[i - 1].my total = total + math.sqrt(dx * dx + dy * dy) samples[i].s = total end samples.total = total return samples end local function point_on_label_path(path, pos) if not path or #path == 0 then return nil end if #path == 1 or (path.total or 0) <= 1e-14 then local p = path[1] return {x = p.x, y = p.y, lx = p.lx, ly = p.ly, mx = p.mx, my = p.my, tx = 1, ty = 0} end pos = clamp01(pos) local target = pos * path.total local hi = 2 while hi <= #path and path[hi].s < target do hi = hi + 1 end if hi > #path then hi = #path end local lo = math.max(1, hi - 1) local a, b = path[lo], path[hi] local ds = math.max(1e-14, (b.s or 0) - (a.s or 0)) local r = (target - (a.s or 0)) / ds if r < 0 then r = 0 elseif r > 1 then r = 1 end local lx = a.lx + (b.lx - a.lx) * r local ly = a.ly + (b.ly - a.ly) * r local mx = a.mx + (b.mx - a.mx) * r local my = a.my + (b.my - a.my) * r local tx = b.mx - a.mx local ty = b.my - a.my local norm = math.sqrt(tx * tx + ty * ty) if norm <= 1e-14 then tx, ty = 1, 0 else tx, ty = tx / norm, ty / norm end return {x = path.log_x and 10 ^ lx or lx, y = path.log_y and 10 ^ ly or ly, lx = lx, ly = ly, mx = mx, my = my, tx = tx, ty = ty} end local function label_point(points, pos, log_x, log_y) local state = point_on_label_path( make_label_path(points, nil, log_x, log_y), pos) if not state then return nil end return {x = state.x, y = state.y} end local function label_segment(points, pos, log_x, log_y) if #points < 2 then return nil, nil end local delta = 0.004 local p1 = math.max(0, clamp01(pos) - delta) local p2 = math.min(1, clamp01(pos) + delta) if p1 == p2 then p1 = math.max(0, p1 - delta); p2 = math.min(1, p2 + delta) end local path = make_label_path(points, nil, log_x, log_y) local a = point_on_label_path(path, p1) local b = point_on_label_path(path, p2) if not a or not b then return nil, nil end return {x = a.x, y = a.y}, {x = b.x, y = b.y} end local function add_curve_label(texout, id, style, points, pos, text, coord_digits, sloped, allow_upside_down, log_x, log_y) if type(style) == "table" and type(points) ~= "table" then local old_style, old_points, old_pos, old_text = id, style, points, pos local old_coord_digits, old_sloped, old_allow = text, coord_digits, sloped id, style, points, pos, text, coord_digits, sloped, allow_upside_down = nil, old_style, old_points, old_pos, old_text, old_coord_digits, old_sloped, old_allow end if sloped and #points >= 2 then local a, b = label_segment(points, pos, log_x, log_y) if a and b then local name_style = id and ("name=" .. id .. "-label") or nil local slope_style = join_style(name_style, style, "sloped", allow_upside_down and "allow upside down" or nil) texout[#texout + 1] = "\\path (axis cs:" .. fmt_num(a.x, coord_digits) .. "," .. fmt_num(a.y, coord_digits) .. ") -- " .. "node[midway," .. slope_style .. "] {" .. text .. "} " .. "(axis cs:" .. fmt_num(b.x, coord_digits) .. "," .. fmt_num(b.y, coord_digits) .. "); " return end end add_label(texout, id, style, label_point(points, pos, log_x, log_y), text, coord_digits) end local function text_len_for_label(s) s = tostring(s or "") s = s:gsub("%b{}", function(x) return x:sub(2, -2) end) s = s:gsub("\\[a-zA-Z]+", "") :gsub("[%$%^_{}]", "") :gsub("\\,", "") return #s end local function normalized_label_position(points, pos, bounds) local state = point_on_label_path(make_label_path(points, bounds, bounds and bounds.log_x, bounds and bounds.log_y), pos) if not state then return nil end return {x = state.mx, y = state.my, point = {x = state.x, y = state.y}, tx = state.tx, ty = state.ty} end local function curve_bounds_from_records(records, opts) local xmin, xmax, ymin, ymax = math.huge, -math.huge, math.huge, -math.huge for _, rec in ipairs(records) do for _, pt in ipairs(rec.points or {}) do if finite_plot_point(pt) then if pt.x < xmin then xmin = pt.x end if pt.x > xmax then xmax = pt.x end if pt.y < ymin then ymin = pt.y end if pt.y > ymax then ymax = pt.y end end end end if xmin == math.huge or xmax <= xmin then xmin, xmax = 0, 1 end if ymin == math.huge or ymax <= ymin then ymin, ymax = 0.1, 10 end local xpad = (xmax - xmin) * to_number(opts.label_bounds_xpad, 0.03) xmin, xmax = xmin - xpad, xmax + xpad return {xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax} end local function make_candidate_positions(base, opts, ov) local manual = override_number(ov, "label_pos", nil) if manual then return {clamp01(manual)} end base = clamp01(base or 0.5) local n = math.max(1, math.floor(to_number(opts.label_candidate_count, 21))) local pmin = to_number(opts.label_min_pos, 0) local pmax = to_number(opts.label_max_pos, 1) local span = math.max(base - pmin, pmax - base) local out, seen = {}, {} local function addp(p) p = clamp01(math.max(pmin, math.min(pmax, p))) local k = string.format("%.4f", p) if not seen[k] then out[#out + 1] = p; seen[k] = true end end addp(base) for i = 1, n do local d = span * i / n addp(base + d) addp(base - d) end return out end local function label_box(rec, opts) local len = text_len_for_label(rec.text) local w = to_number(opts.label_box_width, nil) if not w then w = to_number(opts.label_box_base_width, 0.045) + len * to_number(opts.label_box_char_width, 0.0060) end local h = to_number(opts.label_box_height, 0.060) local mult = to_number(opts.label_box_sloped_multiplier, rec.sloped and 1.35 or 1.0) return w * mult, h * mult end local function overlap_penalty(candidate, placed, sep) local penalty = 0 for _, p in ipairs(placed) do local wx = 0.5 * (candidate.w + p.w) + sep local hy = 0.5 * (candidate.h + p.h) + sep local dx = math.abs(candidate.x - p.x) local dy = math.abs(candidate.y - p.y) if dx < wx and dy < hy then local ox = (wx - dx) / wx local oy = (hy - dy) / hy penalty = penalty + (ox * oy) ^ 2 end end return penalty end --- Emit label records without performing automatic placement. -- -- Render labels at their preferred positions without collision avoidance. -- -- Automatic label placement is delegated to pgfplots-autonode. This helper is -- used only when the user explicitly selects manual label placement and for -- simple process labels that are not routed through autonode. local function render_label_records(texout, records, opts) local list = {} for _, rec in ipairs(records or {}) do list[#list + 1] = rec end table.sort(list, function(a, b) return (a.serial or 0) < (b.serial or 0) end) for _, rec in ipairs(list) do local pos = override_number(rec.override or {}, "label_pos", rec.base_pos or 0.5) if not rec.skip_label then add_curve_label(texout, rec.id, rec.style, rec.points, pos, rec.text, rec.coord_digits or 6, rec.sloped, rec.allow_upside_down, opts and opts.log_x, opts and opts.log_y) end end end local function make_temperature_values(opts, pmin, _pmax, c) local unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) local offset = (unit == "celsius" or unit == "c" or unit == "degc") and 273.15 or 0 local mode = to_string(opts.temperature_mode or opts.t_mode, "linear"):lower() local explicit = opts.temperature_values or opts.t_values if explicit ~= nil and explicit ~= "" and explicit ~= "auto" then mode = "list" end local Tmin_auto = c.Ttriple and c.Ttriple + 2 or 250 local Tsat_min = nil if pmin and c.Tcrit and c.Ttriple then Tsat_min = safe_propsSI("T", "P", pmin, "Q", 0, opts.fluid or "R134a") end if Tsat_min then Tmin_auto = math.max(Tmin_auto, Tsat_min) end local Tmax_auto = c.Tcrit and math.min((c.Tmax or (c.Tcrit + 80)) - 1, c.Tcrit + 40) or 420 local Tmin_user = to_number(opts.temperature_min or opts.t_min, nil) local Tmax_user = to_number(opts.temperature_max or opts.t_max, nil) local Tstep_user = to_number(opts.temperature_step or opts.t_step, 10) local Tcount_user = to_number(opts.temperature_count or opts.t_count, nil) local preset_name = M._strict.choice_value(opts.temperature_preset or opts.t_preset, "default", "temperature preset", {"default", "refrigeration", "sparse", "dense"}) local Tmin_display = Tmin_user or (Tmin_auto - offset) local Tmax_display = Tmax_user or (Tmax_auto - offset) local preset_values if preset_name == "default" then preset_values = nil elseif preset_name == "refrigeration" then preset_values = {-40, -20, 0, 20, 40, 60, 80, 100, 120} elseif preset_name == "sparse" then preset_values = {-40, 0, 40, 80, 120} elseif preset_name == "dense" then preset_values = make_linear_values(-40, 120, 10) end local values_display = make_grid_values({ mode = mode, values = explicit, min = Tmin_display, max = Tmax_display, step = Tstep_user, count = Tcount_user, preset_values = preset_values, }) local out = {} for _, v in ipairs(values_display) do out[#out + 1] = v + offset end table.sort(out) return out end -- Keep formatting policy in the TeX frontend. Lua emits private dispatch -- hooks with both a siunitx unit expression and a generic-TeX fallback; -- luacoolprop.sty selects siunitx while plain TeX and ConTeXt ignore it. local function tex_format_number(value) return "\\csname LCP@format@number\\endcsname{" .. value .. "}" end local function tex_format_quantity(value, siunitx_unit, generic_unit) return "\\csname LCP@format@quantity\\endcsname{" .. value .. "}{" .. siunitx_unit .. "}{" .. generic_unit .. "}" end local function thermodynamic_symbol(opts, canonical_name, short_name, default) return to_string(opts[canonical_name] or opts[short_name], default) end local function format_temperature_label(T, opts) local unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) local symbol = thermodynamic_symbol(opts, "temperature_symbol", "t_symbol", "T") if unit == "kelvin" or unit == "k" then return "$" .. symbol .. "=" .. tex_format_quantity(fmt_t_value(T), "\\kelvin", "{\\rm K}") .. "$" end return "$" .. symbol .. "=" .. tex_format_quantity(fmt_t_value(T - 273.15), "\\degreeCelsius", "{}^{\\circ}{\\rm C}") .. "$" end local function make_entropy_values(opts, pmin, pmax, c) local fluid = opts.fluid or "R134a" local unit = M._strict.choice_value(opts.entropy_unit, "kjkgk", "entropy unit", {"kjkgk", "si", "jkgk", "j/kg/k"}) local factor = (unit == "si" or unit == "jkgk" or unit == "j/kg/k") and 1 or 1e-3 local mode = to_string(opts.entropy_mode or opts.s_mode, "linear"):lower() local explicit = opts.entropy_values or opts.s_values if explicit ~= nil and explicit ~= "" and explicit ~= "auto" then mode = "list" end local psat_ref = pmin if c.pcrit then psat_ref = math.min(psat_ref, c.pcrit * 0.999) end local sf = safe_propsSI("S", "P", psat_ref, "Q", 0, fluid) local sg = safe_propsSI("S", "P", psat_ref, "Q", 1, fluid) local Tref = c.Tcrit and math.min((c.Tmax or (c.Tcrit + 100)) - 1, c.Tcrit + 80) or nil local s_super = (Tref and pmin > 0) and safe_propsSI("S", "T", Tref, "P", pmin, fluid) or nil local smin_auto_si = sf or 800 local smax_auto_si = math.max(sg or smin_auto_si * 1.8, s_super or -1e99) if not finite_number(smin_auto_si) then smin_auto_si = 800 end if not finite_number(smax_auto_si) or smax_auto_si <= smin_auto_si then smax_auto_si = smin_auto_si + 800 end local smin_user = to_number(opts.entropy_min or opts.s_min, nil) local smax_user = to_number(opts.entropy_max or opts.s_max, nil) local sstep_user = to_number(opts.entropy_step or opts.s_step, factor == 1 and 100 or 0.1) local scount_user = to_number(opts.entropy_count or opts.s_count, nil) local preset_name = M._strict.choice_value(opts.entropy_preset or opts.s_preset, "default", "entropy preset", {"default", "refrigeration", "sparse", "dense"}) local smin_display = smin_user or (smin_auto_si * factor) local smax_display = smax_user or (smax_auto_si * factor) local preset_values if preset_name == "default" then preset_values = nil elseif preset_name == "sparse" then if factor == 1 then preset_values = make_linear_values(900, 1800, 200) else preset_values = make_linear_values(0.9, 1.8, 0.2) end elseif preset_name == "dense" then if factor == 1 then preset_values = make_linear_values(800, 1900, 100) else preset_values = make_linear_values(0.8, 1.9, 0.1) end elseif preset_name == "refrigeration" then if factor == 1 then preset_values = make_linear_values(900, 1800, 100) else preset_values = make_linear_values(0.9, 1.8, 0.1) end end local values_display = make_grid_values({ mode = mode, values = explicit, min = smin_display, max = smax_display, step = sstep_user, count = scount_user, preset_values = preset_values, }) local out = {} for _, v in ipairs(values_display) do out[#out + 1] = v / factor end table.sort(out) return out end local function format_entropy_label(s, opts) local unit = M._strict.choice_value(opts.entropy_unit, "kjkgk", "entropy unit", {"kjkgk", "si", "jkgk", "j/kg/k"}) local symbol = thermodynamic_symbol(opts, "entropy_symbol", "s_symbol", "s") if unit == "si" or unit == "jkgk" or unit == "j/kg/k" then return "$" .. symbol .. "=" .. tex_format_quantity(fmt_s_value(s), "\\joule\\per\\kilogram\\per\\kelvin", "{\\rm J\\,kg^{-1}\\,K^{-1}}") .. "$" end return "$" .. symbol .. "=" .. tex_format_quantity(fmt_s_value(s * 1e-3), "\\kilo\\joule\\per\\kilogram\\per\\kelvin", "{\\rm kJ\\,kg^{-1}\\,K^{-1}}") .. "$" end -- Select constant-enthalpy values for TS background curves. User-facing -- values default to kJ/kg while CoolProp always receives J/kg. Automatic -- limits are derived from actual fluid states, never from a refrigerant- -- specific table, so the same code also covers water, hydrocarbons, and -- future CoolProp fluids with a different enthalpy reference. local function make_enthalpy_values(opts, pmin, pmax, constants) local fluid = opts.fluid or "R134a" local unit = M._strict.choice_value(opts.enthalpy_unit or opts.h_unit, "kjkg", "enthalpy unit", {"kjkg", "si", "jkg", "j/kg"}) local factor = (unit == "si" or unit == "jkg" or unit == "j/kg") and 1 or 1e-3 local mode = to_string(opts.enthalpy_mode or opts.h_mode, "linear"):lower() local explicit = opts.enthalpy_values or opts.h_values if explicit ~= nil and explicit ~= "" and explicit ~= "auto" then mode = "list" end local saturation_pressure = pmin if constants.pcrit then saturation_pressure = math.min(saturation_pressure, constants.pcrit * 0.999) end local candidates = { safe_propsSI("H", "P", saturation_pressure, "Q", 0, fluid), safe_propsSI("H", "P", saturation_pressure, "Q", 1, fluid), } if constants.Tcrit then local upper_temperature = math.min( (constants.Tmax or constants.Tcrit + 100) - 1, constants.Tcrit + 60) candidates[#candidates + 1] = safe_propsSI("H", "P", pmin, "T", upper_temperature, fluid) candidates[#candidates + 1] = safe_propsSI("H", "P", pmax, "T", upper_temperature, fluid) end local hmin_si, hmax_si = math.huge, -math.huge for _, value in ipairs(candidates) do if finite_number(value) then hmin_si = math.min(hmin_si, value) hmax_si = math.max(hmax_si, value) end end if hmin_si == math.huge then hmin_si = 1e5 end if hmax_si == -math.huge or hmax_si <= hmin_si then hmax_si = hmin_si + 4e5 end local span = hmax_si - hmin_si hmin_si = hmin_si - 0.05 * span hmax_si = hmax_si + 0.05 * span local hmin = to_number(opts.enthalpy_min or opts.h_min, hmin_si * factor) local hmax = to_number(opts.enthalpy_max or opts.h_max, hmax_si * factor) local step = to_number(opts.enthalpy_step or opts.h_step, nil) local count = to_number(opts.enthalpy_count or opts.h_count, 8) local values = make_grid_values({ mode = mode, values = explicit, min = hmin, max = hmax, step = step, count = count, }) local out = {} for _, value in ipairs(values) do out[#out + 1] = value / factor end table.sort(out) return out end local function format_enthalpy_label(enthalpy, opts) local unit = M._strict.choice_value(opts.enthalpy_unit or opts.h_unit, "kjkg", "enthalpy unit", {"kjkg", "si", "jkg", "j/kg"}) local symbol = thermodynamic_symbol(opts, "enthalpy_symbol", "h_symbol", "h") if unit == "si" or unit == "jkg" or unit == "j/kg" then return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(enthalpy, 2), "\\joule\\per\\kilogram", "{\\rm J\\,kg^{-1}}") .. "$" end return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(enthalpy * 1e-3, 2), "\\kilo\\joule\\per\\kilogram", "{\\rm kJ\\,kg^{-1}}") .. "$" end local function make_isobar_values(opts, pmin, pmax) local unit = M._strict.choice_value(opts.isobar_unit, "bar", "isobar unit", {"bar", "si", "pa", "kpa", "mpa"}) local factor = 1e-5 if unit == "si" or unit == "pa" then factor = 1 elseif unit == "kpa" then factor = 1e-3 elseif unit == "mpa" then factor = 1e-6 end local explicit = opts.isobar_values local mode = to_string(opts.isobar_mode, "log"):lower() if explicit ~= nil and explicit ~= "" and explicit ~= "auto" then mode = "list" end local values = make_grid_values({ mode = mode, values = explicit, min = to_number(opts.isobar_min, pmin * factor), max = to_number(opts.isobar_max, pmax * factor), step = to_number(opts.isobar_step, nil), count = to_number(opts.isobar_count, 7), }) local out = {} for _, value in ipairs(values) do local pressure = value / factor if pressure <= 0 then error("luacoolprop: isobar values must be strictly positive", 3) end out[#out + 1] = pressure end table.sort(out) return out end local function format_pressure_label(pressure, opts) local unit = M._strict.choice_value(opts.isobar_unit, "bar", "isobar unit", {"bar", "si", "pa", "kpa", "mpa"}) local symbol = thermodynamic_symbol(opts, "pressure_symbol", "p_symbol", "p") if unit == "si" or unit == "pa" then return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(pressure, 3), "\\pascal", "{\\rm Pa}") .. "$" elseif unit == "kpa" then return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(pressure * 1e-3, 3), "\\kilo\\pascal", "{\\rm kPa}") .. "$" elseif unit == "mpa" then return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(pressure * 1e-6, 3), "\\mega\\pascal", "{\\rm MPa}") .. "$" end return "$" .. symbol .. "=" .. tex_format_quantity( fmt_num(pressure * 1e-5, 3), "\\bar", "{\\rm bar}") .. "$" end local function make_specific_volume_values(opts, pmin, pmax, c) local fluid = opts.fluid or "R134a" local unit = M._strict.choice_value( opts.specific_volume_unit or opts.volume_unit or opts.v_unit, "m3kg", "specific-volume unit", {"m3kg", "lkg", "l/kg", "dm3kg", "dm3/kg"}) local factor = 1 if unit == "lkg" or unit == "l/kg" or unit == "dm3kg" or unit == "dm3/kg" then factor = 1000 -- m^3/kg -> L/kg for display/input end local mode = to_string(opts.specific_volume_mode or opts.volume_mode or opts.v_mode, "log"):lower() local explicit = opts.specific_volume_values or opts.volume_values or opts.v_values if explicit ~= nil and explicit ~= "" and explicit ~= "auto" then mode = "list" end local psat_low = pmin if c.pcrit then psat_low = math.min(psat_low, c.pcrit * 0.999) end local rho_liq = safe_propsSI("Dmass", "P", psat_low, "Q", 0, fluid) local rho_vap = safe_propsSI("Dmass", "P", psat_low, "Q", 1, fluid) local vmin_auto_si = (rho_liq and rho_liq > 0) and 1 / rho_liq or 8e-4 local vmax_auto_si = (rho_vap and rho_vap > 0) and 1 / rho_vap or 2e-1 vmin_auto_si = vmin_auto_si * 0.95 vmax_auto_si = vmax_auto_si * 1.20 if not finite_number(vmin_auto_si) or vmin_auto_si <= 0 then vmin_auto_si = 8e-4 end if not finite_number(vmax_auto_si) or vmax_auto_si <= vmin_auto_si then vmax_auto_si = vmin_auto_si * 200 end local vmin_user = to_number(opts.specific_volume_min or opts.volume_min or opts.v_min, nil) local vmax_user = to_number(opts.specific_volume_max or opts.volume_max or opts.v_max, nil) local vstep_user = to_number(opts.specific_volume_step or opts.volume_step or opts.v_step, nil) local vcount_user = to_number(opts.specific_volume_count or opts.volume_count or opts.v_count, 8) local preset_name = M._strict.choice_value( opts.specific_volume_preset or opts.volume_preset or opts.v_preset, "default", "specific-volume preset", {"default", "refrigeration", "sparse", "dense"}) local vmin_display = vmin_user or (vmin_auto_si * factor) local vmax_display = vmax_user or (vmax_auto_si * factor) local preset_values if preset_name == "default" or preset_name == "refrigeration" then preset_values = (factor == 1) and {0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2} or {1,2,5,10,20,50,100,200} elseif preset_name == "sparse" then preset_values = (factor == 1) and {0.001,0.005,0.02,0.1} or {1,5,20,100} elseif preset_name == "dense" then preset_values = (factor == 1) and {0.0008,0.001,0.0015,0.002,0.003,0.005,0.007,0.01,0.015,0.02,0.03,0.05,0.07,0.1,0.15,0.2} or {0.8,1,1.5,2,3,5,7,10,15,20,30,50,70,100,150,200} end local values_display = make_grid_values({ mode = mode, values = explicit, min = vmin_display, max = vmax_display, step = vstep_user, count = vcount_user, preset_values = preset_values, }) local out = {} for _, v in ipairs(values_display) do local si = v / factor if si <= 0 then error("luacoolprop: specific-volume values must be strictly positive", 3) end out[#out + 1] = si end table.sort(out) return out end local function format_specific_volume_label(v, opts) local unit = M._strict.choice_value( opts.specific_volume_unit or opts.volume_unit or opts.v_unit, "m3kg", "specific-volume unit", {"m3kg", "lkg", "l/kg", "dm3kg", "dm3/kg"}) local symbol = thermodynamic_symbol(opts, "specific_volume_symbol", "v_symbol", "v") if unit == "lkg" or unit == "l/kg" or unit == "dm3kg" or unit == "dm3/kg" then return "$" .. symbol .. "=" .. tex_format_quantity(fmt_v_value(v * 1000), "\\litre\\per\\kilogram", "{\\rm L\\,kg^{-1}}") .. "$" end return "$" .. symbol .. "=" .. tex_format_quantity(fmt_v_value(v), "\\cubic\\metre\\per\\kilogram", "{\\rm m^3\\,kg^{-1}}") .. "$" end local function quality_curve_id(q, _opts) return curve_id("lcp-ph", "q", q, "", 4) end local function temperature_curve_id(T, opts) local unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) if unit == "kelvin" or unit == "k" then return curve_id("lcp-ph", "T", T, "K", 3) end return curve_id("lcp-ph", "T", T - 273.15, "C", 3) end local function entropy_curve_id(s, opts) local unit = M._strict.choice_value(opts.entropy_unit, "kjkgk", "entropy unit", {"kjkgk", "si", "jkgk", "j/kg/k"}) if unit == "si" or unit == "jkgk" or unit == "j/kg/k" then return curve_id("lcp-ph", "s", s, "JkgK", 3) end return curve_id("lcp-ph", "s", s * 1e-3, "kJkgK", 4) end local function specific_volume_curve_id(v, opts) local unit = M._strict.choice_value( opts.specific_volume_unit or opts.volume_unit or opts.v_unit, "m3kg", "specific-volume unit", {"m3kg", "lkg", "l/kg", "dm3kg", "dm3/kg"}) if unit == "lkg" or unit == "l/kg" or unit == "dm3kg" or unit == "dm3/kg" then return curve_id("lcp-ph", "v", v * 1000, "Lkg", 4) end return curve_id("lcp-ph", "v", v, "m3kg", 5) end local function pv_quality_curve_id(q, _opts) return curve_id("lcp-pv", "q", q, "", 4) end local function pv_temperature_curve_id(T, opts) local unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) if unit == "kelvin" or unit == "k" then return curve_id("lcp-pv", "T", T, "K", 3) end return curve_id("lcp-pv", "T", T - 273.15, "C", 3) end local function pv_entropy_curve_id(s, opts) local unit = M._strict.choice_value(opts.entropy_unit, "kjkgk", "entropy unit", {"kjkgk", "si", "jkgk", "j/kg/k"}) if unit == "si" or unit == "jkgk" or unit == "j/kg/k" then return curve_id("lcp-pv", "s", s, "JkgK", 3) end return curve_id("lcp-pv", "s", s * 1e-3, "kJkgK", 4) end local function ts_quality_curve_id(quality, _opts) return curve_id("lcp-ts", "q", quality, "", 4) end local function ts_enthalpy_curve_id(enthalpy, opts) local unit = M._strict.choice_value(opts.enthalpy_unit or opts.h_unit, "kjkg", "enthalpy unit", {"kjkg", "si", "jkg", "j/kg"}) if unit == "si" or unit == "jkg" or unit == "j/kg" then return curve_id("lcp-ts", "h", enthalpy, "Jkg", 2) end return curve_id("lcp-ts", "h", enthalpy * 1e-3, "kJkg", 3) end local function hs_quality_curve_id(quality, _opts) return curve_id("lcp-hs", "q", quality, "", 4) end local function hs_temperature_curve_id(temperature, opts) local unit = M._strict.choice_value(opts.temperature_unit, "celsius", "temperature unit", {"celsius", "c", "degc", "kelvin", "k"}) if unit == "kelvin" or unit == "k" then return curve_id("lcp-hs", "T", temperature, "K", 3) end return curve_id("lcp-hs", "T", temperature - 273.15, "C", 3) end local function hs_volume_curve_id(volume, opts) local unit = M._strict.choice_value(opts.specific_volume_unit or opts.v_unit, "m3kg", "specific-volume unit", {"m3kg", "lkg", "l/kg", "dm3kg", "dm3/kg"}) if unit == "lkg" or unit == "l/kg" then return curve_id("lcp-hs", "v", volume * 1000, "Lkg", 4) end return curve_id("lcp-hs", "v", volume, "m3kg", 5) end local function hs_pressure_curve_id(pressure, opts) local unit = M._strict.choice_value(opts.isobar_unit, "bar", "isobar unit", {"bar", "si", "pa", "kpa", "mpa"}) if unit == "si" or unit == "pa" then return curve_id("lcp-hs", "p", pressure, "Pa", 3) elseif unit == "kpa" then return curve_id("lcp-hs", "p", pressure * 1e-3, "kPa", 3) elseif unit == "mpa" then return curve_id("lcp-hs", "p", pressure * 1e-6, "MPa", 3) end return curve_id("lcp-hs", "p", pressure * 1e-5, "bar", 3) end function M.diagram._pt_entropy_curve_id(entropy, opts) local unit = M._strict.choice_value(opts.entropy_unit, "kjkgk", "entropy unit", {"kjkgk", "si", "jkgk", "j/kg/k"}) if unit == "si" or unit == "jkgk" or unit == "j/kg/k" then return curve_id("lcp-pt", "s", entropy, "JkgK", 3) end return curve_id("lcp-pt", "s", entropy * 1e-3, "kJkgK", 4) end function M.diagram._pt_enthalpy_curve_id(enthalpy, opts) local unit = M._strict.choice_value(opts.enthalpy_unit or opts.h_unit, "kjkg", "enthalpy unit", {"kjkg", "si", "jkg", "j/kg"}) if unit == "si" or unit == "jkg" or unit == "j/kg" then return curve_id("lcp-pt", "h", enthalpy, "Jkg", 2) end return curve_id("lcp-pt", "h", enthalpy * 1e-3, "kJkg", 3) end function M.diagram._pt_volume_curve_id(volume, opts) local unit = M._strict.choice_value(opts.specific_volume_unit or opts.v_unit, "m3kg", "specific-volume unit", {"m3kg", "lkg", "l/kg", "dm3kg", "dm3/kg"}) if unit == "lkg" or unit == "l/kg" or unit == "dm3kg" or unit == "dm3/kg" then return curve_id("lcp-pt", "v", volume * 1000, "Lkg", 4) end return curve_id("lcp-pt", "v", volume, "m3kg", 5) end --- Return true when label placement is delegated to pgfplots-autonode. -- In this mode the thermodynamic curve is still sampled by LuaCoolProp, but -- the label is emitted as a PGFPlots trailing path command so that -- pgfplots-autonode can sample the final PGFPlots path geometry. local function is_autonode_placement(opts) local placement = to_string(opts and opts.label_placement, "autonode"):lower() return not (placement == "none" or placement == "off" or placement == "manual") end --- Translate a LuaCoolProp label record into pgfplots-autonode keys. -- The function keeps LuaCoolProp's per-curve override semantics: manual -- positions become preferred positions, and fixed labels are represented by -- a degenerate admissible interval. local function autonode_bracket_options(rec, opts) opts = opts or {} local ov = rec.override or {} local preferred = override_number(ov, "label_pos", rec.base_pos or 0.5) local fixed = override_bool(ov, "label_fixed", false) local minpos = fixed and preferred or to_number(opts.label_min_pos, 0) local maxpos = fixed and preferred or to_number(opts.label_max_pos, 1) local node_options = rec.style or "luacoolprop label node" local samples = to_number(opts.autonode_candidates, to_number(opts.label_candidate_count, 51)) local strategy = to_string(opts.autonode_candidate_strategy, "around-preferred") local sloped = rec.sloped and "true" or "false" local allow_upside_down = rec.allow_upside_down and "true" or "false" local priority = rec.priority or 0 return join_style( "preferred pos=" .. fmt_num(clamp01(preferred), 5), "min pos=" .. fmt_num(clamp01(minpos), 5), "max pos=" .. fmt_num(clamp01(maxpos), 5), "samples=" .. tostring(math.max(1, math.floor(samples))), "candidate strategy=" .. strategy, "normal shift=" .. to_string(opts.autonode_normal_shift, "0pt"), "clearance=" .. to_string(opts.autonode_clearance, "1pt"), "inner sep=" .. to_string(opts.autonode_inner_sep, "1.5pt"), "node options={" .. node_options .. "}", "preferred weight=" .. to_string(opts.autonode_preferred_weight, "8"), "overlap weight=" .. to_string(opts.autonode_overlap_weight, "1000"), "priority=" .. tostring(priority), "sloped=" .. sloped, "allow upside down=" .. allow_upside_down ) end local function add_coords_plot(texout, id, style, points, coord_digits) texout[#texout + 1] = "\\addplot[" .. join_style(plot_name_style(id), style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(points) do append_coord(coords, pt, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " end --- Emit one curve plus an optional pgfplots-autonode trailing label command. -- This must be emitted immediately after the corresponding \addplot data so -- that PGFPlots can associate the autonode sampler with the current plot path. local function add_coords_plot_with_autonode(texout, rec, opts) texout[#texout + 1] = "\\addplot[" .. join_style(plot_name_style(rec.id), rec.style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(rec.points or {}) do append_coord(coords, pt, rec.coord_digits or 6) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}" if rec.label and rec.label_text then texout[#texout + 1] = "\\pgfplotsautonode[" .. autonode_bracket_options({ style = rec.label_style, base_pos = rec.base_pos, sloped = rec.sloped, allow_upside_down = rec.allow_upside_down, priority = rec.priority, override = rec.override, }, opts) .. "]{" .. rec.label_text .. "}" end texout[#texout + 1] = "; " end local function add_record(records, rec) rec.serial = #records + 1 records[#records + 1] = rec end local function copy_table(t) local out = {} for key, value in pairs(t or {}) do out[key] = value end return out end local function collect_ph_records(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin_axis, pmax_axis, c = resolve_pressure_limits(opts, fluid) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local curve_style = opts.curve_style or "" local records = {} local function final_curve_style(id, base) local ov = curve_override(opts, id) local style = join_style(base, override_string(ov, "curve_style", override_string(ov, "style", nil))) if forget_plot and not legend then style = join_style(style, "forget plot") end return style, ov end if to_bool(opts.isentropes, false) then local pmin, pmax = pmin_axis, pmax_axis if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isentrope_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isentrope_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isentrope_tolerance or opts.tolerance, 0.30), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local labels = to_bool(opts.isentrope_labels, to_bool(opts.labels, false)) local sloped = to_bool(opts.isentrope_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow = to_bool(opts.isentrope_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local every = math.max(1, math.floor(to_number(opts.isentrope_label_every or opts.label_every, 2))) local base_pos = to_number(opts.isentrope_label_pos or opts.label_pos, 0.45) local color = opts.isentrope_color or "green!50!black" local base_style = join_style(color, opts.isentrope_style or "line width=0.25pt,densely dotted") local label_style = opts.isentrope_label_style or opts.label_style or "luacoolprop isentrope label node" local values = make_entropy_values(opts, pmin, pmax, c) for i, s in ipairs(values) do local points = M.diagram.ph_isentrope_curve(gen, s) if #points >= 2 then local id = entropy_curve_id(s, opts) local style, ov = final_curve_style(id, base_style) local show = labels and (((i - 1) % every == 0) or i == #values) show = override_bool(ov, "label_show", override_bool(ov, "show_label", show)) add_record(records, { id = id, family = "isentrope", value = s, points = points, style = style, label = show, label_style = join_style(label_style, override_string(ov, "label_style", nil)), label_text = override_string(ov, "label_text", format_entropy_label(s, opts)), base_pos = override_number(ov, "label_pos", base_pos), sloped = override_bool(ov, "label_sloped", sloped), allow_upside_down = override_bool(ov, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isentrope_label_priority, 30), override = ov, }) end end end if to_bool(opts.isochores, false) then local pmin, pmax = pmin_axis, pmax_axis if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isochore_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isochore_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isochore_tolerance or opts.tolerance, 0.28), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local labels = to_bool(opts.isochore_labels, to_bool(opts.labels, false)) local sloped = to_bool(opts.isochore_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow = to_bool(opts.isochore_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local every = math.max(1, math.floor(to_number(opts.isochore_label_every or opts.label_every, 2))) local base_pos = to_number(opts.isochore_label_pos or opts.label_pos, 0.58) local color = opts.isochore_color or "orange!80!black" local base_style = join_style(color, opts.isochore_style or "line width=0.25pt,dash pattern=on 2pt off 1.4pt") local label_style = opts.isochore_label_style or opts.label_style or "luacoolprop isochore label node" local values = make_specific_volume_values(opts, pmin, pmax, c) for i, v in ipairs(values) do local points = M.diagram.ph_isochore_curve(gen, v) if #points >= 2 then local id = specific_volume_curve_id(v, opts) local style, ov = final_curve_style(id, base_style) local show = labels and (((i - 1) % every == 0) or i == #values) show = override_bool(ov, "label_show", override_bool(ov, "show_label", show)) add_record(records, { id = id, family = "isochore", value = v, points = points, style = style, label = show, label_style = join_style(label_style, override_string(ov, "label_style", nil)), label_text = override_string(ov, "label_text", format_specific_volume_label(v, opts)), base_pos = override_number(ov, "label_pos", base_pos), sloped = override_bool(ov, "label_sloped", sloped), allow_upside_down = override_bool(ov, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isochore_label_priority, 40), override = ov, }) end end end if to_bool(opts.isotherms, false) then local pmin, pmax = pmin_axis, pmax_axis if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.35), log_weight = to_number(opts.log_weight, 30), saturation_pressure_epsilon = to_number(opts.saturation_pressure_epsilon, 1e-5), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local labels = to_bool(opts.isotherm_labels, to_bool(opts.labels, false)) local sloped = to_bool(opts.isotherm_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow = to_bool(opts.isotherm_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local every = math.max(1, math.floor(to_number(opts.isotherm_label_every or opts.label_every, 2))) local base_pos = to_number(opts.isotherm_label_pos or opts.label_pos, 0.20) local color = opts.isotherm_color or "red!65!black" local base_style = join_style(color, opts.isotherm_style or "line width=0.25pt") local label_style = opts.isotherm_label_style or opts.label_style or "luacoolprop isotherm label node" local values = make_temperature_values(opts, pmin, pmax, c) for i, T in ipairs(values) do if (not c.Tmin or T >= c.Tmin - 1e-9) and (not c.Tmax or T <= c.Tmax + 1e-9) then local points = M.diagram.ph_isotherm_curve(gen, T) if #points >= 2 then local id = temperature_curve_id(T, opts) local style, ov = final_curve_style(id, base_style) local show = labels and (((i - 1) % every == 0) or i == #values) show = override_bool(ov, "label_show", override_bool(ov, "show_label", show)) add_record(records, { id = id, family = "isotherm", value = T, points = points, style = style, label = show, label_style = join_style(label_style, override_string(ov, "label_style", nil)), label_text = override_string(ov, "label_text", format_temperature_label(T, opts)), base_pos = override_number(ov, "label_pos", base_pos), sloped = override_bool(ov, "label_sloped", sloped), allow_upside_down = override_bool(ov, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isotherm_label_priority, 20), override = ov, }) end end end end if to_bool(opts.isoquality, true) then local pmax_quality = math.min(pmax_axis, c.pcrit or pmax_axis) local pmin = pmin_axis if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end if pmax_quality <= pmin then pmax_quality = pmin * 1.01 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax_quality, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.initial_intervals, 18), max_depth = to_number(opts.max_depth, 8), tolerance = to_number(opts.tolerance, 0.25), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local labels = to_bool(opts.quality_labels, to_bool(opts.labels, false)) local sloped = to_bool(opts.quality_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow = to_bool(opts.quality_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local every = math.max(1, math.floor(to_number(opts.quality_label_every or opts.label_every, 1))) local base_pos = to_number(opts.quality_label_pos or opts.label_pos, 0.35) local qc = opts.quality_color or "blue!60!black" local qbc = opts.quality_boundary_color or opts.boundary_color or "blue!80!black" local interior_style = join_style(qc, opts.interior_style or "line width=0.25pt") local boundary_style = join_style(qbc, opts.boundary_style or "line width=0.9pt") local label_style = opts.quality_label_style or opts.label_style or "luacoolprop quality label node" local quality_symbol = thermodynamic_symbol(opts, "quality_symbol", "q_symbol", "Q") local values = make_q_values(opts) for i, q in ipairs(values) do if q >= -1e-12 and q <= 1 + 1e-12 then local points = M.diagram.ph_isoquality_curve(gen, q) if #points >= 2 then local id = quality_curve_id(q, opts) local base = join_style(curve_style, (math.abs(q) < 1e-12 or math.abs(q - 1) < 1e-12) and boundary_style or interior_style) local style, ov = final_curve_style(id, base) local show = labels and (((i - 1) % every == 0) or i == #values) show = override_bool(ov, "label_show", override_bool(ov, "show_label", show)) add_record(records, { id = id, family = "quality", value = q, points = points, style = style, label = show, label_style = join_style(label_style, override_string(ov, "label_style", nil)), label_text = override_string(ov, "label_text", "$" .. quality_symbol .. "=" .. tex_format_number(fmt_q(q)) .. "$"), base_pos = override_number(ov, "label_pos", base_pos), sloped = override_bool(ov, "label_sloped", sloped), allow_upside_down = override_bool(ov, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.quality_label_priority, 10), override = ov, }) end end end end return records end -- Collect PV families into the format-independent curve-record model shared -- by every diagram renderer. local function collect_pv_records(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin_axis, pmax_axis, constants = resolve_pressure_limits(opts, fluid) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local curve_style = opts.curve_style or "" local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) local records = {} local function final_curve_style(id, base) local override = curve_override(opts, id) local style = join_style(base, override_string(override, "curve_style", override_string(override, "style", nil))) if forget_plot and not legend then style = join_style(style, "forget plot") end return style, override end local pmin_background = pmin_axis if constants.ptriple and pmin_background <= constants.ptriple then pmin_background = constants.ptriple * 1.05 end if to_bool(opts.isentropes, false) then local generator = { fluid = fluid, library = opts.library, pressure_min = pmin_background, pressure_max = pmax_axis, v_scale = v_scale, p_scale = p_scale, initial_intervals = to_number( opts.isentrope_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isentrope_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isentrope_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x_weight = to_number(opts.log_x_weight, 30), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local values = make_entropy_values(opts, pmin_background, pmax_axis, constants) local labels = to_bool(opts.isentrope_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isentrope_label_every or opts.label_every, 2))) local base_pos = to_number( opts.isentrope_label_pos or opts.label_pos, 0.45) local sloped = to_bool(opts.isentrope_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)) local allow = to_bool(opts.isentrope_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local base_style = join_style(curve_style, opts.isentrope_color or "green!50!black", opts.isentrope_style or "line width=0.25pt,densely dotted") local label_style = opts.isentrope_label_style or opts.label_style or "luacoolprop isentrope label node" for index, entropy in ipairs(values) do local points = M.diagram.pv_isentrope_curve(generator, entropy) if #points >= 2 then local id = pv_entropy_curve_id(entropy, opts) local style, override = final_curve_style(id, base_style) local show = labels and (((index - 1) % every == 0) or index == #values) show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = "isentrope", value = entropy, points = points, style = style, label = show, label_style = join_style(label_style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", format_entropy_label(entropy, opts)), base_pos = override_number(override, "label_pos", base_pos), sloped = override_bool(override, "label_sloped", sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isentrope_label_priority, 30), override = override, }) end end end if to_bool(opts.isotherms, false) then local generator = { fluid = fluid, library = opts.library, pressure_min = pmin_background, pressure_max = pmax_axis, v_scale = v_scale, p_scale = p_scale, initial_intervals = to_number( opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x_weight = to_number(opts.log_x_weight, 30), saturation_pressure_epsilon = to_number( opts.saturation_pressure_epsilon, 1e-5), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local values = make_temperature_values(opts, pmin_background, pmax_axis, constants) local labels = to_bool(opts.isotherm_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isotherm_label_every or opts.label_every, 2))) local base_pos = to_number( opts.isotherm_label_pos or opts.label_pos, 0.20) local sloped = to_bool(opts.isotherm_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)) local allow = to_bool(opts.isotherm_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local base_style = join_style(curve_style, opts.isotherm_color or "red!65!black", opts.isotherm_style or "line width=0.25pt") local label_style = opts.isotherm_label_style or opts.label_style or "luacoolprop isotherm label node" for index, temperature in ipairs(values) do if (not constants.Tmin or temperature >= constants.Tmin - 1e-9) and (not constants.Tmax or temperature <= constants.Tmax + 1e-9) then local points = M.diagram.pv_isotherm_curve(generator, temperature) if #points >= 2 then local id = pv_temperature_curve_id(temperature, opts) local style, override = final_curve_style(id, base_style) local show = labels and (((index - 1) % every == 0) or index == #values) show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = "isotherm", value = temperature, points = points, style = style, label = show, label_style = join_style(label_style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", format_temperature_label(temperature, opts)), base_pos = override_number(override, "label_pos", base_pos), sloped = override_bool(override, "label_sloped", sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isotherm_label_priority, 20), override = override, }) end end end end if to_bool(opts.isoquality, true) then local pmax_quality = math.min(pmax_axis, constants.pcrit or pmax_axis) local pmin_quality = pmin_background if pmax_quality <= pmin_quality then pmax_quality = pmin_quality * 1.01 end local generator = { fluid = fluid, library = opts.library, pressure_min = pmin_quality, pressure_max = pmax_quality, v_scale = v_scale, p_scale = p_scale, initial_intervals = to_number( opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x_weight = to_number(opts.log_x_weight, 30), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local values = make_q_values(opts) local labels = to_bool(opts.quality_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.quality_label_every or opts.label_every, 1))) local base_pos = to_number( opts.quality_label_pos or opts.label_pos, 0.35) local sloped = to_bool(opts.quality_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)) local allow = to_bool(opts.quality_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local interior_style = join_style(curve_style, opts.quality_color or "blue!60!black", opts.interior_style or "line width=0.25pt") local boundary_style = join_style(curve_style, opts.quality_boundary_color or opts.boundary_color or "blue!80!black", opts.boundary_style or "line width=0.9pt") local label_style = opts.quality_label_style or opts.label_style or "luacoolprop quality label node" local quality_symbol = thermodynamic_symbol(opts, "quality_symbol", "q_symbol", "Q") for index, quality in ipairs(values) do if quality >= -1e-12 and quality <= 1 + 1e-12 then local points = M.diagram.pv_quality_curve(generator, quality) if #points >= 2 then local id = pv_quality_curve_id(quality, opts) local base = (math.abs(quality) < 1e-12 or math.abs(quality - 1) < 1e-12) and boundary_style or interior_style local style, override = final_curve_style(id, base) local show = labels and (((index - 1) % every == 0) or index == #values) show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = "quality", value = quality, points = points, style = style, label = show, label_style = join_style(label_style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", "$" .. quality_symbol .. "=" .. tex_format_number(fmt_q(quality)) .. "$"), base_pos = override_number(override, "label_pos", base_pos), sloped = override_bool(override, "label_sloped", sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.quality_label_priority, 10), override = override, }) end end end end return records end -- Collect TS families into the format-independent curve-record model shared -- by every diagram renderer. Only sampling and coordinate conversion remain -- projection-specific. local function collect_ts_records(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin_axis, pmax_axis, constants = resolve_pressure_limits(opts, fluid) local pmin = pmin_axis if constants.ptriple and pmin <= constants.ptriple then pmin = constants.ptriple * 1.05 end local s_scale = to_number(opts.s_scale, 1e-3) local temperature_scale = to_number(opts.temperature_scale or opts.t_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local curve_style = opts.curve_style or "" local records = {} local function final_curve_style(id, base) local override = curve_override(opts, id) local style = join_style(base, override_string(override, "curve_style", override_string(override, "style", nil))) if forget_plot and not legend then style = join_style(style, "forget plot") end return style, override end if to_bool(opts.isenthalps, false) then local generator = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax_axis, s_scale = s_scale, temperature_scale = temperature_scale, temperature_offset = temperature_offset, initial_intervals = to_number( opts.isenthalp_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isenthalp_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isenthalp_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 8), temperature_weight = to_number(opts.temperature_weight, 0.01), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local values = make_enthalpy_values(opts, pmin, pmax_axis, constants) local labels = to_bool(opts.isenthalp_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isenthalp_label_every or opts.label_every, 2))) local base_pos = to_number( opts.isenthalp_label_pos or opts.label_pos, 0.65) local sloped = to_bool(opts.isenthalp_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)) local allow = to_bool(opts.isenthalp_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local base_style = join_style(curve_style, opts.isenthalp_color or "red!65!black", opts.isenthalp_style or "line width=0.25pt,densely dashed") local label_style = opts.isenthalp_label_style or opts.label_style or "luacoolprop isenthalp label node" for index, enthalpy in ipairs(values) do local points = M.diagram.ts_isenthalp_curve(generator, enthalpy) if #points >= 2 then local id = ts_enthalpy_curve_id(enthalpy, opts) local style, override = final_curve_style(id, base_style) local show = labels and (((index - 1) % every == 0) or index == #values) show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = "isenthalp", value = enthalpy, points = points, style = style, label = show, label_style = join_style(label_style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", format_enthalpy_label(enthalpy, opts)), base_pos = override_number(override, "label_pos", base_pos), sloped = override_bool(override, "label_sloped", sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.isenthalp_label_priority, 20), override = override, }) end end end if to_bool(opts.isoquality, true) then local pmax = math.min(pmax_axis, constants.pcrit or pmax_axis) if pmax <= pmin then pmax = pmin * 1.01 end local generator = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, s_scale = s_scale, temperature_scale = temperature_scale, temperature_offset = temperature_offset, initial_intervals = to_number( opts.quality_initial_intervals or opts.initial_intervals, 18), max_depth = to_number(opts.quality_max_depth or opts.max_depth, 8), tolerance = to_number(opts.quality_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 8), temperature_weight = to_number(opts.temperature_weight, 0.01), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local values = make_q_values(opts) local labels = to_bool(opts.quality_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.quality_label_every or opts.label_every, 1))) local base_pos = to_number( opts.quality_label_pos or opts.label_pos, 0.35) local sloped = to_bool(opts.quality_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)) local allow = to_bool(opts.quality_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local interior_style = join_style(curve_style, opts.quality_color or "blue!60!black", opts.interior_style or "line width=0.25pt") local boundary_style = join_style(curve_style, opts.quality_boundary_color or opts.boundary_color or "blue!80!black", opts.boundary_style or "line width=0.9pt") local label_style = opts.quality_label_style or opts.label_style or "luacoolprop quality label node" local quality_symbol = thermodynamic_symbol(opts, "quality_symbol", "q_symbol", "Q") for index, quality in ipairs(values) do if quality >= -1e-12 and quality <= 1 + 1e-12 then local points = M.diagram.ts_quality_curve(generator, quality) if #points >= 2 then local id = ts_quality_curve_id(quality, opts) local base = (math.abs(quality) < 1e-12 or math.abs(quality - 1) < 1e-12) and boundary_style or interior_style local style, override = final_curve_style(id, base) local show = labels and (((index - 1) % every == 0) or index == #values) show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = "quality", value = quality, points = points, style = style, label = show, label_style = join_style(label_style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", "$" .. quality_symbol .. "=" .. tex_format_number(fmt_q(quality)) .. "$"), base_pos = override_number(override, "label_pos", base_pos), sloped = override_bool(override, "label_sloped", sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", allow), coord_digits = coord_digits, priority = to_number(opts.quality_label_priority, 10), override = override, }) end end end end return records end local function collect_hs_records(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax, constants = resolve_pressure_limits(opts, fluid) if constants.ptriple and pmin <= constants.ptriple then pmin = constants.ptriple * 1.05 end local s_scale = to_number(opts.s_scale, 1e-3) local h_scale = to_number(opts.h_scale, 1e-3) local digits = math.floor(to_number(opts.coord_digits, 6)) local legend = to_bool(opts.legend, false) local forget = to_bool(opts.forget_plot, true) local curve_style = opts.curve_style or "" local records = {} local function styled(id, base) local override = curve_override(opts, id) local style = join_style(base, override_string(override, "curve_style", override_string(override, "style", nil))) if forget and not legend then style = join_style(style, "forget plot") end return style, override end local function record(family, value, points, id, base_style, label_spec) if #points < 2 then return end local style, override = styled(id, base_style) local show = label_spec.enabled and label_spec.selected show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = family, value = value, points = points, style = style, label = show, label_style = join_style(label_spec.style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", label_spec.text), base_pos = override_number(override, "label_pos", label_spec.pos), sloped = override_bool(override, "label_sloped", label_spec.sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", label_spec.allow), coord_digits = digits, priority = label_spec.priority, override = override, }) end local common = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, s_scale = s_scale, h_scale = h_scale, initial_intervals = to_number(opts.initial_intervals, 18), max_depth = to_number(opts.max_depth, 8), tolerance = to_number(opts.tolerance, 0.10), enthalpy_weight = to_number(opts.enthalpy_weight, 0.01), saturation_pressure_epsilon = to_number( opts.saturation_pressure_epsilon, 1e-5), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } if to_bool(opts.isobars, false) then local values = make_isobar_values(opts, pmin, pmax) local enabled = to_bool(opts.isobar_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isobar_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isobar_temperature_min = opts.isobar_temperature_min generator.isobar_temperature_max = opts.isobar_temperature_max generator.temperature_unit = opts.temperature_unit generator.isobar_initial_intervals = opts.isobar_initial_intervals generator.isobar_max_depth = opts.isobar_max_depth generator.isobar_tolerance = opts.isobar_tolerance for index, pressure in ipairs(values) do record("isobar", pressure, M.diagram.hs_isobar_curve(generator, pressure), hs_pressure_curve_id(pressure, opts), join_style(curve_style, opts.isobar_color or "purple!70!black", opts.isobar_style or "line width=0.25pt,dash pattern=on 3pt off 1.5pt"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isobar_label_style or opts.label_style or "luacoolprop isobar label node", text = format_pressure_label(pressure, opts), pos = to_number(opts.isobar_label_pos or opts.label_pos, 0.68), sloped = to_bool(opts.isobar_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isobar_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isobar_label_priority, 40), }) end end if to_bool(opts.isochores, false) then local values = make_specific_volume_values(opts, pmin, pmax, constants) local enabled = to_bool(opts.isochore_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isochore_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isochore_initial_intervals = opts.isochore_initial_intervals generator.isochore_max_depth = opts.isochore_max_depth generator.isochore_tolerance = opts.isochore_tolerance for index, volume in ipairs(values) do record("isochore", volume, M.diagram.hs_isochore_curve(generator, volume), hs_volume_curve_id(volume, opts), join_style(curve_style, opts.isochore_color or "orange!80!black", opts.isochore_style or "line width=0.25pt,densely dotted"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isochore_label_style or opts.label_style or "luacoolprop isochore label node", text = format_specific_volume_label(volume, opts), pos = to_number(opts.isochore_label_pos or opts.label_pos, 0.58), sloped = to_bool(opts.isochore_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isochore_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isochore_label_priority, 30), }) end end if to_bool(opts.isotherms, false) then local values = make_temperature_values(opts, pmin, pmax, constants) local enabled = to_bool(opts.isotherm_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isotherm_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isotherm_initial_intervals = opts.isotherm_initial_intervals generator.isotherm_max_depth = opts.isotherm_max_depth generator.isotherm_tolerance = opts.isotherm_tolerance for index, temperature in ipairs(values) do if (not constants.Tmin or temperature >= constants.Tmin) and (not constants.Tmax or temperature <= constants.Tmax) then record("isotherm", temperature, M.diagram.hs_isotherm_curve(generator, temperature), hs_temperature_curve_id(temperature, opts), join_style(curve_style, opts.isotherm_color or "red!65!black", opts.isotherm_style or "line width=0.25pt"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isotherm_label_style or opts.label_style or "luacoolprop isotherm label node", text = format_temperature_label(temperature, opts), pos = to_number(opts.isotherm_label_pos or opts.label_pos, 0.25), sloped = to_bool(opts.isotherm_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isotherm_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isotherm_label_priority, 20), }) end end end if to_bool(opts.isoquality, true) then local generator = copy_table(common) generator.pressure_max = math.min(pmax, constants.pcrit or pmax) local values = make_q_values(opts) local enabled = to_bool(opts.quality_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.quality_label_every or opts.label_every, 1))) local quality_symbol = thermodynamic_symbol(opts, "quality_symbol", "q_symbol", "Q") for index, quality in ipairs(values) do if quality >= 0 and quality <= 1 then local boundary = math.abs(quality) < 1e-12 or math.abs(quality - 1) < 1e-12 record("quality", quality, M.diagram.hs_quality_curve(generator, quality), hs_quality_curve_id(quality, opts), join_style(curve_style, boundary and (opts.quality_boundary_color or "blue!80!black") or (opts.quality_color or "blue!60!black"), boundary and (opts.boundary_style or "line width=0.9pt") or (opts.interior_style or "line width=0.25pt")), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.quality_label_style or opts.label_style or "luacoolprop quality label node", text = "$" .. quality_symbol .. "=" .. tex_format_number(fmt_q(quality)) .. "$", pos = to_number(opts.quality_label_pos or opts.label_pos, 0.35), sloped = to_bool(opts.quality_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.quality_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.quality_label_priority, 10), }) end end end return records end -- Collect PT families in the same renderer-neutral record format as every -- other built-in projection. The phase envelope is deliberately one curve: -- vapour quality is not an independent coordinate in the PT plane. function M.diagram._collect_pt_records(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax, constants = M.diagram._resolve_pt_pressure_limits(opts, fluid) local temperature_scale = to_number(opts.temperature_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local p_scale = to_number(opts.p_scale, 1e-5) local digits = math.floor(to_number(opts.coord_digits, 6)) local legend = to_bool(opts.legend, false) local forget = to_bool(opts.forget_plot, true) local curve_style = opts.curve_style or "" local records = {} local function styled(id, base) local override = curve_override(opts, id) local style = join_style(base, override_string(override, "curve_style", override_string(override, "style", nil))) if forget and not legend then style = join_style(style, "forget plot") end return style, override end local function record(family, value, points, id, base_style, spec) if #points < 2 then return end local style, override = styled(id, base_style) local show = spec.enabled and spec.selected show = override_bool(override, "label_show", override_bool(override, "show_label", show)) add_record(records, { id = id, family = family, value = value, points = points, style = style, label = show, label_style = join_style(spec.style, override_string(override, "label_style", nil)), label_text = override_string(override, "label_text", spec.text), base_pos = override_number(override, "label_pos", spec.pos), sloped = override_bool(override, "label_sloped", spec.sloped), allow_upside_down = override_bool(override, "label_allow_upside_down", spec.allow), coord_digits = digits, priority = spec.priority, override = override, }) end local common = { fluid = fluid, library = opts.library, reference_state = opts.reference_state, pressure_min = pmin, pressure_max = pmax, temperature_scale = temperature_scale, temperature_offset = temperature_offset, p_scale = p_scale, initial_intervals = to_number(opts.initial_intervals, 18), max_depth = to_number(opts.max_depth, 8), tolerance = to_number(opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 30), constants = constants, domain_policy = opts.domain_policy or opts.discontinuity_policy, } if to_bool(opts.isentropes, false) then local values = make_entropy_values(opts, pmin, pmax, constants) local enabled = to_bool(opts.isentrope_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isentrope_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isentrope_initial_intervals = opts.isentrope_initial_intervals generator.isentrope_max_depth = opts.isentrope_max_depth generator.isentrope_tolerance = opts.isentrope_tolerance for index, entropy in ipairs(values) do record("isentrope", entropy, M.diagram.pt_isentrope_curve(generator, entropy), M.diagram._pt_entropy_curve_id(entropy, opts), join_style(curve_style, opts.isentrope_color or "green!50!black", opts.isentrope_style or "line width=0.25pt,densely dotted"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isentrope_label_style or opts.label_style or "luacoolprop isentrope label node", text = format_entropy_label(entropy, opts), pos = to_number(opts.isentrope_label_pos or opts.label_pos, 0.45), sloped = to_bool(opts.isentrope_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isentrope_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isentrope_label_priority, 20), }) end end if to_bool(opts.isenthalps, false) then local values = make_enthalpy_values(opts, pmin, pmax, constants) local enabled = to_bool(opts.isenthalp_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isenthalp_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isenthalp_initial_intervals = opts.isenthalp_initial_intervals generator.isenthalp_max_depth = opts.isenthalp_max_depth generator.isenthalp_tolerance = opts.isenthalp_tolerance for index, enthalpy in ipairs(values) do record("isenthalp", enthalpy, M.diagram.pt_isenthalp_curve(generator, enthalpy), M.diagram._pt_enthalpy_curve_id(enthalpy, opts), join_style(curve_style, opts.isenthalp_color or "red!65!black", opts.isenthalp_style or "line width=0.25pt,densely dashed"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isenthalp_label_style or opts.label_style or "luacoolprop isenthalp label node", text = format_enthalpy_label(enthalpy, opts), pos = to_number(opts.isenthalp_label_pos or opts.label_pos, 0.62), sloped = to_bool(opts.isenthalp_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isenthalp_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isenthalp_label_priority, 30), }) end end if to_bool(opts.isochores, false) then local values = make_specific_volume_values(opts, pmin, pmax, constants) local enabled = to_bool(opts.isochore_labels, to_bool(opts.labels, false)) local every = math.max(1, math.floor(to_number( opts.isochore_label_every or opts.label_every, 2))) local generator = copy_table(common) generator.isochore_initial_intervals = opts.isochore_initial_intervals generator.isochore_max_depth = opts.isochore_max_depth generator.isochore_tolerance = opts.isochore_tolerance for index, volume in ipairs(values) do record("isochore", volume, M.diagram.pt_isochore_curve(generator, volume), M.diagram._pt_volume_curve_id(volume, opts), join_style(curve_style, opts.isochore_color or "orange!80!black", opts.isochore_style or "line width=0.25pt,dash pattern=on 2pt off 1.4pt"), { enabled = enabled, selected = ((index - 1) % every == 0) or index == #values, style = opts.isochore_label_style or opts.label_style or "luacoolprop isochore label node", text = format_specific_volume_label(volume, opts), pos = to_number(opts.isochore_label_pos or opts.label_pos, 0.58), sloped = to_bool(opts.isochore_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.isochore_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.isochore_label_priority, 40), }) end end if to_bool(opts.phase_envelope, true) then local generator = copy_table(common) generator.phase_envelope_initial_intervals = opts.phase_envelope_initial_intervals generator.phase_envelope_max_depth = opts.phase_envelope_max_depth generator.phase_envelope_tolerance = opts.phase_envelope_tolerance record("phase_envelope", nil, M.diagram.pt_phase_envelope_curve(generator), "lcp-pt-phase-envelope", join_style(curve_style, opts.phase_envelope_color or "blue!80!black", opts.phase_envelope_style or "line width=0.9pt"), { enabled = to_bool(opts.phase_envelope_labels, to_bool(opts.labels, false)), selected = true, style = opts.phase_envelope_label_style or opts.label_style or "luacoolprop phase envelope label node", text = opts.phase_envelope_label_text or "liquid--vapour equilibrium", pos = to_number(opts.phase_envelope_label_pos or opts.label_pos, 0.45), sloped = to_bool(opts.phase_envelope_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, true)), allow = to_bool(opts.phase_envelope_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)), priority = to_number(opts.phase_envelope_label_priority, 10), }) end return records end local function render_diagram_records(records, opts) local texout, labels = {}, {} if is_autonode_placement(opts) then for _, rec in ipairs(records) do add_coords_plot_with_autonode(texout, rec, opts) end return table.concat(texout) end -- Explicit manual path: no automatic collision avoidance is performed here. for _, rec in ipairs(records) do add_coords_plot(texout, rec.id, rec.style, rec.points, rec.coord_digits or 6) if rec.label then labels[#labels + 1] = { id = rec.id, points = rec.points, style = rec.label_style, text = rec.label_text, base_pos = rec.base_pos, sloped = rec.sloped, allow_upside_down = rec.allow_upside_down, coord_digits = rec.coord_digits, priority = rec.priority, serial = rec.serial, override = rec.override, } end end render_label_records(texout, labels, opts) return table.concat(texout) end local function label_record_from_curve_record(rec) return { id = rec.id, points = rec.points, style = rec.label_style, text = rec.label_text, base_pos = rec.base_pos, sloped = rec.sloped, allow_upside_down = rec.allow_upside_down, coord_digits = rec.coord_digits, priority = rec.priority, serial = rec.serial, override = rec.override, } end local function ph_family_options(opts, family) local out = copy_table(opts or {}) out.isoquality = false out.isotherms = false out.isentropes = false out.isochores = false if family == "quality" or family == "isoquality" or family == "q" then out.isoquality = true elseif family == "isotherm" or family == "isotherms" or family == "temperature" then out.isotherms = true elseif family == "isentrope" or family == "isentropes" or family == "entropy" then out.isentropes = true elseif family == "isochore" or family == "isochores" or family == "specific_volume" or family == "volume" then out.isochores = true else error("luacoolprop: unknown PH curve family '" .. tostring(family) .. "'", 3) end return out end local function render_diagram_curve_records_only(records, opts) local texout = {} if is_autonode_placement(opts) then for _, rec in ipairs(records) do add_coords_plot_with_autonode(texout, rec, opts) end else for _, rec in ipairs(records) do add_coords_plot(texout, rec.id, rec.style, rec.points, rec.coord_digits or 6) end end return table.concat(texout) end --- Serialize the selected PH quality curves as PGFPlots commands. -- This direct PH entry point accepts the compact option vocabulary. The -- registered `M.diagram.ph.quality_plots` entry point also accepts canonical -- option names and participates in generic diagram dispatch. -- @tparam[opt={}] table opts PH rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph_isoquality_plots(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin_axis, pmax_axis, c = resolve_pressure_limits(opts, fluid) local pmax_quality = math.min(pmax_axis, c.pcrit or pmax_axis) local pmin = pmin_axis if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end if pmax_quality <= pmin then pmax_quality = pmin * 1.01 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax_quality, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.initial_intervals, 18), max_depth = to_number(opts.max_depth, 8), tolerance = to_number(opts.tolerance, 0.25), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local labels = to_bool(opts.quality_labels, to_bool(opts.labels, false)) local sloped_labels = to_bool(opts.quality_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow_upside_down = to_bool(opts.quality_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local label_every = math.max(1, math.floor(to_number(opts.quality_label_every or opts.label_every, 1))) local label_pos = to_number(opts.quality_label_pos or opts.label_pos, 0.35) if label_pos < 0 then label_pos = 0 end if label_pos > 1 then label_pos = 1 end local curve_style = opts.curve_style or "" local quality_color = opts.quality_color or "blue!60!black" local quality_boundary_color = opts.quality_boundary_color or opts.boundary_color or "blue!80!black" local interior_style = join_style(quality_color, opts.interior_style or "line width=0.25pt") local boundary_style = join_style(quality_boundary_color, opts.boundary_style or "line width=0.9pt") local label_style = opts.quality_label_style or opts.label_style or "luacoolprop quality label node" local quality_symbol = thermodynamic_symbol(opts, "quality_symbol", "q_symbol", "Q") local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local texout = {} local q_values = make_q_values(opts) for q_index, q in ipairs(q_values) do if q >= -1e-12 and q <= 1 + 1e-12 then local points = M.diagram.ph_isoquality_curve(gen, q) if #points >= 2 then local is_boundary = math.abs(q) < 1e-12 or math.abs(q - 1) < 1e-12 local style = join_style(curve_style, is_boundary and boundary_style or interior_style) if forget_plot and not legend then style = join_style(style, "forget plot") end texout[#texout + 1] = "\\addplot[" .. join_style(style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(points) do append_coord(coords, pt, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " if labels and ((q_index - 1) % label_every == 0 or q_index == #q_values) then add_curve_label(texout, label_style, points, label_pos, "$" .. quality_symbol .. "=" .. tex_format_number(fmt_q(q)) .. "$", coord_digits, sloped_labels, allow_upside_down) end else texout[#texout + 1] = "\\relax " end end end return table.concat(texout) end --- Serialize the selected PH isotherms as PGFPlots commands. -- @tparam[opt={}] table opts PH rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph_isotherm_plots(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax, c = resolve_pressure_limits(opts, fluid) if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isotherm_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isotherm_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isotherm_tolerance or opts.tolerance, 0.35), log_weight = to_number(opts.log_weight, 30), saturation_pressure_epsilon = to_number(opts.saturation_pressure_epsilon, 1e-5), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local labels = to_bool(opts.isotherm_labels, to_bool(opts.labels, false)) local sloped_labels = to_bool(opts.isotherm_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow_upside_down = to_bool(opts.isotherm_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local label_every = math.max(1, math.floor(to_number(opts.isotherm_label_every or opts.label_every, 2))) local label_pos = to_number(opts.isotherm_label_pos or opts.label_pos, 0.20) local isotherm_color = opts.isotherm_color or "red!65!black" local style_base = join_style(isotherm_color, opts.isotherm_style or "line width=0.25pt") local label_style = opts.isotherm_label_style or opts.label_style or "luacoolprop isotherm label node" local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local t_values = make_temperature_values(opts, pmin, pmax, c) local texout = {} for i, T in ipairs(t_values) do if (not c.Tmin or T >= c.Tmin - 1e-9) and (not c.Tmax or T <= c.Tmax + 1e-9) then local points = M.diagram.ph_isotherm_curve(gen, T) if #points >= 2 then local style = style_base if forget_plot and not legend then style = join_style(style, "forget plot") end texout[#texout + 1] = "\\addplot[" .. join_style(style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(points) do append_coord(coords, pt, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " if labels and ((i - 1) % label_every == 0 or i == #t_values) then add_curve_label(texout, label_style, points, label_pos, format_temperature_label(T, opts), coord_digits, sloped_labels, allow_upside_down) end end end end return table.concat(texout) end --- Serialize the selected PH isentropes as PGFPlots commands. -- @tparam[opt={}] table opts PH rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph_isentrope_plots(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax, c = resolve_pressure_limits(opts, fluid) if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isentrope_initial_intervals or opts.initial_intervals, 14), max_depth = to_number(opts.isentrope_max_depth or opts.max_depth, 7), tolerance = to_number(opts.isentrope_tolerance or opts.tolerance, 0.30), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local labels = to_bool(opts.isentrope_labels, to_bool(opts.labels, false)) local sloped_labels = to_bool(opts.isentrope_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow_upside_down = to_bool(opts.isentrope_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local label_every = math.max(1, math.floor(to_number(opts.isentrope_label_every or opts.label_every, 2))) local label_pos = to_number(opts.isentrope_label_pos or opts.label_pos, 0.45) local isentrope_color = opts.isentrope_color or "green!50!black" local style_base = join_style(isentrope_color, opts.isentrope_style or "line width=0.25pt,densely dotted") local label_style = opts.isentrope_label_style or opts.label_style or "luacoolprop isentrope label node" local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local s_values = make_entropy_values(opts, pmin, pmax, c) local texout = {} for i, s in ipairs(s_values) do local points = M.diagram.ph_isentrope_curve(gen, s) if #points >= 2 then local style = style_base if forget_plot and not legend then style = join_style(style, "forget plot") end texout[#texout + 1] = "\\addplot[" .. join_style(style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(points) do append_coord(coords, pt, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " if labels and ((i - 1) % label_every == 0 or i == #s_values) then add_curve_label(texout, label_style, points, label_pos, format_entropy_label(s, opts), coord_digits, sloped_labels, allow_upside_down) end end end return table.concat(texout) end --- Serialize the selected PH isochores as PGFPlots commands. -- @tparam[opt={}] table opts PH rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph_isochore_plots(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax, c = resolve_pressure_limits(opts, fluid) if c.ptriple and pmin <= c.ptriple then pmin = c.ptriple * 1.05 end local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = to_number(opts.h_scale, 1e-3), p_scale = to_number(opts.p_scale, 1e-5), initial_intervals = to_number(opts.isochore_initial_intervals or opts.initial_intervals, 16), max_depth = to_number(opts.isochore_max_depth or opts.max_depth, 8), tolerance = to_number(opts.isochore_tolerance or opts.tolerance, 0.28), log_weight = to_number(opts.log_weight, 30), constants = c, domain_policy = opts.domain_policy or opts.discontinuity_policy, } local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local labels = to_bool(opts.isochore_labels, to_bool(opts.labels, false)) local sloped_labels = to_bool(opts.isochore_label_sloped, to_bool(opts.label_sloped or opts.sloped_labels, false)) local allow_upside_down = to_bool(opts.isochore_label_allow_upside_down, to_bool(opts.label_allow_upside_down, false)) local label_every = math.max(1, math.floor(to_number(opts.isochore_label_every or opts.label_every, 2))) local label_pos = to_number(opts.isochore_label_pos or opts.label_pos, 0.58) local isochore_color = opts.isochore_color or "orange!80!black" local style_base = join_style(isochore_color, opts.isochore_style or "line width=0.25pt,dash pattern=on 2pt off 1.4pt") local label_style = opts.isochore_label_style or opts.label_style or "luacoolprop isochore label node" local legend = to_bool(opts.legend, false) local forget_plot = to_bool(opts.forget_plot, true) local v_values = make_specific_volume_values(opts, pmin, pmax, c) local texout = {} for i, v in ipairs(v_values) do local points = M.diagram.ph_isochore_curve(gen, v) if #points >= 2 then local style = style_base if forget_plot and not legend then style = join_style(style, "forget plot") end texout[#texout + 1] = "\\addplot[" .. join_style(style, "unbounded coords=jump") .. "] coordinates {" local coords = {} for _, pt in ipairs(points) do append_coord(coords, pt, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " if labels and ((i - 1) % label_every == 0 or i == #v_values) then add_curve_label(texout, label_style, points, label_pos, format_specific_volume_label(v, opts), coord_digits, sloped_labels, allow_upside_down) end end end return table.concat(texout) end --- Serialize all enabled PH curve families as PGFPlots commands. -- Curves are collected in one pass so that names and autonode labels remain -- stable across family combinations. -- @tparam[opt={}] table opts PH rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph_plots(opts) opts = opts or {} -- The high-level plotter uses a unified pass. This gives all curves stable -- names and emits pgfplots-autonode labels immediately after their plots. return render_diagram_records(collect_ph_records(opts), opts) end --- Return data-dependent PGFPlots limits for a PH axis. -- @tparam[opt={}] table opts Options including `fluid`, `pressure_min`, -- `pressure_max`, `pressure_scale`, and `coord_digits`. -- @treturn string Comma-separated PGFPlots key list containing `ymin` and `ymax`. function M.diagram.ph_axis_style(opts) opts = opts or {} M.load_library(opts.library) local fluid = opts.fluid or "R134a" local pmin, pmax = resolve_pressure_limits(opts, fluid) local p_scale = to_number(opts.p_scale, 1e-5) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) return "ymin=" .. fmt_num(pmin * p_scale, coord_digits) .. ",ymax=" .. fmt_num(pmax * p_scale, coord_digits) end local function pv_family_options(opts, family) local out = copy_table(opts or {}) out.isoquality = false out.isotherms = false out.isentropes = false if family == "quality" or family == "isoquality" or family == "q" then out.isoquality = true elseif family == "isotherm" or family == "isotherms" or family == "temperature" then out.isotherms = true elseif family == "isentrope" or family == "isentropes" or family == "entropy" then out.isentropes = true else error("luacoolprop: unknown PV curve family '" .. tostring(family) .. "'", 3) end return out end --- Serialize all enabled PV curve families as PGFPlots commands. -- @tparam[opt={}] table opts PV rendering, sampling, style, and label options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.pv_plots(opts) opts = opts or {} local render_opts = copy_table(opts) render_opts.log_x = true return render_diagram_records(collect_pv_records(opts), render_opts) end --- Serialize one selected PV curve family as PGFPlots commands. -- @tparam[opt={}] table opts PV rendering, sampling, style, and label options. -- @tparam string family `quality`, `isotherm`, or `isentrope`. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.pv_family(opts, family) opts = opts or {} local family_opts = pv_family_options(opts, family) local render_opts = copy_table(opts) render_opts.log_x = true return render_diagram_curve_records_only(collect_pv_records(family_opts), render_opts) end --- Serialize selected PV quality curves. -- @tparam[opt={}] table opts PV options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.pv_quality_plots(opts) return M.diagram.pv_family(opts or {}, "quality") end --- Serialize selected PV isotherms. -- @tparam[opt={}] table opts PV options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.pv_isotherm_plots(opts) return M.diagram.pv_family(opts or {}, "isotherm") end --- Serialize selected PV isentropes. -- @tparam[opt={}] table opts PV options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.pv_isentrope_plots(opts) return M.diagram.pv_family(opts or {}, "isentrope") end --- Return data-dependent PGFPlots pressure limits for a PV axis. -- The generic TeX layer supplies `xmode=log,ymode=log`; this function only -- scales the physical pressure range. -- @tparam[opt={}] table opts Options including `fluid`, pressure limits, -- `p_scale`, and `coord_digits`. -- @treturn string Comma-separated PGFPlots axis keys. function M.diagram.pv_axis_style(opts) opts = opts or {} M.load_library(opts.library) local pmin, pmax = resolve_pressure_limits(opts, opts.fluid or "R134a") local p_scale = to_number(opts.p_scale, 1e-5) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) return "ymin=" .. fmt_num(pmin * p_scale, coord_digits) .. ",ymax=" .. fmt_num(pmax * p_scale, coord_digits) end local function ts_family_options(opts, family) local out = copy_table(opts or {}) out.isoquality = false out.isenthalps = false if family == "quality" or family == "isoquality" or family == "q" then out.isoquality = true elseif family == "isenthalp" or family == "isenthalps" or family == "enthalpy" then out.isenthalps = true else error("luacoolprop: unknown TS curve family '" .. tostring(family) .. "'", 3) end return out end --- Serialize all enabled TS curve families as PGFPlots commands. -- @tparam[opt={}] table opts TS rendering and sampling options. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ts_plots(opts, prepared_records) opts = opts or {} local render_opts = copy_table(opts) render_opts.log_y = false return render_diagram_records(prepared_records or collect_ts_records(opts), render_opts) end --- Serialize one selected TS family. -- @tparam[opt={}] table opts TS rendering and sampling options. -- @tparam string family `quality` or `isenthalp`. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ts_family(opts, family) opts = opts or {} local render_opts = copy_table(opts) render_opts.log_y = false return render_diagram_curve_records_only( collect_ts_records(ts_family_options(opts, family)), render_opts) end --- Serialize selected TS quality curves. -- @tparam[opt={}] table opts TS options. -- @treturn string PGFPlots commands. function M.diagram.ts_quality_plots(opts) return M.diagram.ts_family(opts or {}, "quality") end --- Serialize selected TS isenthalps. -- @tparam[opt={}] table opts TS options. -- @treturn string PGFPlots commands. function M.diagram.ts_isenthalp_plots(opts) return M.diagram.ts_family(opts or {}, "isenthalp") end --- Return data-dependent linear axis limits for a TS diagram. -- Explicit `entropy_axis_min`, `entropy_axis_max`, `temperature_axis_min`, and -- `temperature_axis_max` values override the automatically padded bounds. -- @tparam[opt={}] table opts TS options. -- @treturn string Comma-separated PGFPlots axis keys. function M.diagram.ts_axis_style(opts) opts = opts or {} local records = collect_ts_records(opts) local xmin, xmax, ymin, ymax = math.huge, -math.huge, math.huge, -math.huge for _, record in ipairs(records) do for _, point in ipairs(record.points or {}) do if finite_plot_point(point) then xmin = math.min(xmin, point.x) xmax = math.max(xmax, point.x) ymin = math.min(ymin, point.y) ymax = math.max(ymax, point.y) end end end if xmin == math.huge or xmax <= xmin then xmin, xmax = 0, 1 end if ymin == math.huge or ymax <= ymin then ymin, ymax = 200, 500 end local xpad = 0.04 * (xmax - xmin) local ypad = 0.04 * (ymax - ymin) xmin = to_number(opts.entropy_axis_min, xmin - xpad) xmax = to_number(opts.entropy_axis_max, xmax + xpad) ymin = to_number(opts.temperature_axis_min, ymin - ypad) ymax = to_number(opts.temperature_axis_max, ymax + ypad) local digits = math.floor(to_number(opts.coord_digits, 6)) return "xmin=" .. fmt_num(xmin, digits) .. ",xmax=" .. fmt_num(xmax, digits) .. ",ymin=" .. fmt_num(ymin, digits) .. ",ymax=" .. fmt_num(ymax, digits), records end local function hs_family_options(opts, family) local out = copy_table(opts or {}) out.isoquality, out.isochores, out.isotherms, out.isobars = false, false, false, false if family == "quality" or family == "isoquality" or family == "q" then out.isoquality = true elseif family == "isochore" or family == "isochores" or family == "volume" then out.isochores = true elseif family == "isotherm" or family == "isotherms" or family == "temperature" then out.isotherms = true elseif family == "isobar" or family == "isobars" or family == "pressure" then out.isobars = true else error("luacoolprop: unknown HS curve family '" .. tostring(family) .. "'", 3) end return out end --- Serialize all enabled HS curve families. -- @tparam[opt={}] table opts HS rendering and sampling options. -- @treturn string PGFPlots commands. function M.diagram.hs_plots(opts, prepared_records) opts = opts or {} local render_opts = copy_table(opts) render_opts.log_y = false return render_diagram_records(prepared_records or collect_hs_records(opts), render_opts) end --- Serialize one selected HS family. -- @tparam[opt={}] table opts HS options. -- @tparam string family `quality`, `isochore`, `isotherm`, or `isobar`. -- @treturn string PGFPlots commands. function M.diagram.hs_family(opts, family) opts = opts or {} local render_opts = copy_table(opts) render_opts.log_y = false return render_diagram_curve_records_only( collect_hs_records(hs_family_options(opts, family)), render_opts) end function M.diagram.hs_quality_plots(opts) return M.diagram.hs_family(opts or {}, "quality") end function M.diagram.hs_isochore_plots(opts) return M.diagram.hs_family(opts or {}, "isochore") end function M.diagram.hs_isotherm_plots(opts) return M.diagram.hs_family(opts or {}, "isotherm") end function M.diagram.hs_isobar_plots(opts) return M.diagram.hs_family(opts or {}, "isobar") end --- Return data-dependent linear axis limits for an HS diagram. -- @tparam[opt={}] table opts HS options and optional explicit axis bounds. -- @treturn string PGFPlots keys. function M.diagram.hs_axis_style(opts) opts = opts or {} local records = collect_hs_records(opts) local xmin, xmax, ymin, ymax = math.huge, -math.huge, math.huge, -math.huge for _, record in ipairs(records) do for _, point in ipairs(record.points or {}) do if finite_number(point.x) and finite_number(point.y) then xmin, xmax = math.min(xmin, point.x), math.max(xmax, point.x) ymin, ymax = math.min(ymin, point.y), math.max(ymax, point.y) end end end if xmin == math.huge or xmax <= xmin then xmin, xmax = 0, 1 end if ymin == math.huge or ymax <= ymin then ymin, ymax = 0, 1 end local xpad, ypad = 0.04 * (xmax - xmin), 0.04 * (ymax - ymin) xmin = to_number(opts.entropy_axis_min, xmin - xpad) xmax = to_number(opts.entropy_axis_max, xmax + xpad) ymin = to_number(opts.enthalpy_axis_min, ymin - ypad) ymax = to_number(opts.enthalpy_axis_max, ymax + ypad) local digits = math.floor(to_number(opts.coord_digits, 6)) return "xmin=" .. fmt_num(xmin, digits) .. ",xmax=" .. fmt_num(xmax, digits) .. ",ymin=" .. fmt_num(ymin, digits) .. ",ymax=" .. fmt_num(ymax, digits), records end function M.diagram._pt_family_options(opts, family) local out = copy_table(opts or {}) out.phase_envelope, out.isentropes, out.isenthalps, out.isochores = false, false, false, false local normalized = tostring(family):lower():gsub("[%s%-]+", "_") if normalized == "phase_envelope" or normalized == "saturation" or normalized == "equilibrium" then out.phase_envelope = true elseif normalized == "isentrope" or normalized == "isentropes" or normalized == "entropy" then out.isentropes = true elseif normalized == "isenthalp" or normalized == "isenthalps" or normalized == "enthalpy" then out.isenthalps = true elseif normalized == "isochore" or normalized == "isochores" or normalized == "volume" then out.isochores = true else error("luacoolprop: unknown PT curve family '" .. tostring(family) .. "'", 3) end return out end --- Serialize all enabled PT curve families. -- The pressure coordinate is logarithmic in the high-level TeX diagram. function M.diagram.pt_plots(opts, prepared_records) opts = opts or {} local render_options = copy_table(opts) render_options.log_y = true return render_diagram_records(prepared_records or M.diagram._collect_pt_records(opts), render_options) end --- Serialize one selected PT family. -- @tparam table opts PT options. -- @tparam string family `phase_envelope`, `isentrope`, `isenthalp`, or -- `isochore`. function M.diagram.pt_family(opts, family) opts = opts or {} local render_options = copy_table(opts) render_options.log_y = true return render_diagram_curve_records_only( M.diagram._collect_pt_records(M.diagram._pt_family_options(opts, family)), render_options) end function M.diagram.pt_phase_envelope_plots(opts) return M.diagram.pt_family(opts or {}, "phase_envelope") end function M.diagram.pt_isentrope_plots(opts) return M.diagram.pt_family(opts or {}, "isentrope") end function M.diagram.pt_isenthalp_plots(opts) return M.diagram.pt_family(opts or {}, "isenthalp") end function M.diagram.pt_isochore_plots(opts) return M.diagram.pt_family(opts or {}, "isochore") end --- Return data-dependent limits for a PT axis. -- Pressure limits are the requested physical interval. Temperature limits -- are derived from all enabled curves and may be overridden explicitly. function M.diagram.pt_axis_style(opts) opts = opts or {} local fluid = opts.fluid or "R134a" local pmin, pmax, constants = M.diagram._resolve_pt_pressure_limits(opts, fluid) local records = M.diagram._collect_pt_records(opts) local xmin, xmax = math.huge, -math.huge for _, record in ipairs(records) do for _, point in ipairs(record.points or {}) do if finite_plot_point(point) then xmin = math.min(xmin, point.x) xmax = math.max(xmax, point.x) end end end local temperature_scale = to_number(opts.temperature_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) if xmin == math.huge or xmax <= xmin then xmin = (constants.Ttriple - temperature_offset) * temperature_scale xmax = (constants.Tcrit - temperature_offset) * temperature_scale end local padding = 0.04 * math.max(1, xmax - xmin) xmin = to_number(opts.temperature_axis_min, xmin - padding) xmax = to_number(opts.temperature_axis_max, xmax + padding) if xmax <= xmin then error("luacoolprop: temperature_axis_max must be greater than temperature_axis_min", 2) end local p_scale = to_number(opts.p_scale, 1e-5) local digits = math.floor(to_number(opts.coord_digits, 6)) return "xmin=" .. fmt_num(xmin, digits) .. ",xmax=" .. fmt_num(xmax, digits) .. ",ymin=" .. fmt_num(pmin * p_scale, digits) .. ",ymax=" .. fmt_num(pmax * p_scale, digits), records end -- --------------------------------------------------------------------------- -- PH process segments for thermodynamic cycles -- --------------------------------------------------------------------------- local function trim(s) s = tostring(s or "") return (s:gsub("^%s+", ""):gsub("%s+$", "")) end local function clean_unit(u) u = tostring(u or ""):lower() u = u:gsub("%s+", "") :gsub("°", "deg") :gsub("\\,", "") :gsub("\\", "") :gsub("%{", "") :gsub("%}", "") :gsub("%/", "") :gsub("%*", "") :gsub("%-", "") :gsub("_", "") :gsub("%^", "") :gsub("[%(%)]", "") return u end local function parse_quantity_value(v, kind) if v == nil or v == "" or v == "auto" then return nil end if type(v) == "number" then if finite_number(v) then return v end error("luacoolprop: thermodynamic quantities must be finite", 3) end local s = trim(v) if s:find(",", 1, true) then error("luacoolprop: invalid quantity '" .. s .. "'; use '.' as the decimal separator", 3) end local x, unit = M._strict.decimal_number_prefix(s) if not x then error("luacoolprop: invalid quantity '" .. s .. "'; expected a finite decimal number, optionally followed by a unit " .. "(scientific notation such as 1.2E5Pa is accepted)", 3) end unit = clean_unit(unit) kind = tostring(kind or ""):lower() if kind == "p" or kind == "pressure" then if unit == "" or unit == "pa" or unit == "si" then return x end if unit == "kpa" then return x * 1e3 end if unit == "mpa" then return x * 1e6 end if unit == "bar" then return x * 1e5 end if unit == "mbar" then return x * 1e2 end error("luacoolprop: unknown pressure unit '" .. unit .. "'; expected Pa, kPa, MPa, bar, or mbar", 3) elseif kind == "t" or kind == "temperature" then if unit == "c" or unit == "degc" or unit == "celsius" then return x + 273.15 end if unit == "" or unit == "k" or unit == "kelvin" or unit == "si" then return x end error("luacoolprop: unknown temperature unit '" .. unit .. "'; expected K, kelvin, C, degC, or celsius", 3) elseif kind == "h" or kind == "enthalpy" then if unit == "kjkg" or unit == "kjkg1" then return x * 1e3 end if unit == "" or unit == "si" or unit == "jkg" or unit == "jkg1" then return x end error("luacoolprop: unknown specific-enthalpy unit '" .. unit .. "'; expected J/kg or kJ/kg", 3) elseif kind == "s" or kind == "entropy" then if unit == "kjkgk" or unit == "kjkg1k1" then return x * 1e3 end if unit == "" or unit == "si" or unit == "jkgk" or unit == "jkg1k1" then return x end error("luacoolprop: unknown specific-entropy unit '" .. unit .. "'; expected J/(kg K) or kJ/(kg K)", 3) elseif kind == "v" or kind == "specific_volume" then if unit == "" or unit == "si" or unit == "m3kg" or unit == "m3kg1" then return x end if unit == "lkg" or unit == "lkg1" or unit == "dm3kg" or unit == "dm3kg1" then return x * 1e-3 end error("luacoolprop: unknown specific-volume unit '" .. unit .. "'; expected m3/kg or L/kg", 3) elseif kind == "rho" or kind == "density" then if unit == "" or unit == "si" or unit == "kgm3" then return x end error("luacoolprop: unknown mass-density unit '" .. unit .. "'; expected kg/m3", 3) elseif kind == "q" or kind == "quality" then if unit == "" then return x end error("luacoolprop: vapour quality is dimensionless; unexpected unit '" .. unit .. "'", 3) end error("luacoolprop: internal error: unknown thermodynamic quantity kind '" .. kind .. "'", 3) end local function canonical_state_key(k) k = tostring(k or ""):lower():gsub("[%s%-]+", "_") if k == "pressure" then return "p" end if k == "temperature" or k == "t" then return "T" end if k == "quality" or k == "x" then return "q" end if k == "entropy" then return "s" end if k == "enthalpy" then return "h" end if k == "volume" or k == "specific_volume" or k == "specificvolume" then return "v" end if k == "density" or k == "rho" or k == "rhomass" then return "rho" end return k end local function parse_state_spec(spec) local st = {_raw = tostring(spec or ""), _props = {}} local body = tostring(spec or "") body = body:gsub("^%s*%{", ""):gsub("%}%s*$", "") for token in body:gmatch("[^,;]+") do local k, v = token:match("^%s*(.-)%s*[=:]%s*(.-)%s*$") if not k or k == "" or not v or v == "" then error("luacoolprop: malformed state field '" .. trim(token) .. "'; expected key=value. Decimal commas are not accepted; use '.'", 3) end local ck = canonical_state_key(k) local kinds = { T = "temperature", p = "pressure", h = "enthalpy", s = "entropy", q = "quality", v = "specific_volume", rho = "density", } local quantity_kind = kinds[ck] if not quantity_kind then error("luacoolprop: unknown state property '" .. trim(k) .. "'; expected p, T, h, s, q, v, or rho", 3) end if st._props[ck] then error("luacoolprop: duplicate state property '" .. trim(k) .. "'", 3) end st[ck] = parse_quantity_value(v, quantity_kind) st._props[ck] = true end if st.rho and not st.v and st.rho > 0 then st.v = 1 / st.rho end return st end local function state_pair(st) if st.p and st.q ~= nil then return "P", st.p, "Q", st.q end if st.T and st.q ~= nil then return "T", st.T, "Q", st.q end if st.p and st.T then return "P", st.p, "T", st.T end if st.p and st.h then return "P", st.p, "H", st.h end if st.p and st.s then return "P", st.p, "S", st.s end if st.p and st.v then return "P", st.p, "Dmass", 1 / st.v end if st.T and st.h then return "T", st.T, "H", st.h end if st.T and st.s then return "T", st.T, "S", st.s end if st.T and st.v then return "T", st.T, "Dmass", 1 / st.v end if st.h and st.s then return "H", st.h, "S", st.s end return nil end function M._strict.coolprop_output_from_pair(fluid, output, n1, v1, n2, v2) if not n1 then return nil end local key = output:lower() local cpout = ({ p = "P", t = "T", h = "H", s = "S", q = "Q", v = "Dmass", rho = "Dmass", })[key] if not cpout then return nil end local val = safe_propsSI(cpout, n1, v1, n2, v2, fluid) if not val then return nil end if key == "v" then if val == 0 then return nil end return 1 / val end return val end local function state_output_from_pair(fluid, st, output, n1, v1, n2, v2) local key = output:lower() if key == "p" and st.p then return st.p end if key == "t" and st.T then return st.T end if key == "h" and st.h then return st.h end if key == "s" and st.s then return st.s end if key == "q" and st.q ~= nil then return st.q end if key == "v" and st.v then return st.v end if key == "rho" and st.rho then return st.rho end return M._strict.coolprop_output_from_pair(fluid, output, n1, v1, n2, v2) end local function state_output(fluid, st, output) local n1, v1, n2, v2 = state_pair(st) return state_output_from_pair(fluid, st, output, n1, v1, n2, v2) end local function state_short(st) local parts = {} if st.p then parts[#parts + 1] = "p" end if st.T then parts[#parts + 1] = "T" end if st.h then parts[#parts + 1] = "h" end if st.s then parts[#parts + 1] = "s" end if st.q ~= nil then parts[#parts + 1] = "Q" end if st.v then parts[#parts + 1] = "v" end if st.rho then parts[#parts + 1] = "rho" end if #parts == 0 then return "no thermodynamic property" end return table.concat(parts, ", ") end local function state_pair_hint(st, role) role = role or "point" if st.q ~= nil and not st.p and not st.T then return "Point '" .. role .. "' contains Q but no saturation coordinate. Add p=... or T=..., for instance " .. role .. "={p=1bar,Q=" .. fmt_q(st.q) .. "}." end if st.p and not (st.T or st.h or st.s or st.v or st.q ~= nil) then return "Point '" .. role .. "' contains only p. Add one of T=..., h=..., s=..., v=... or Q=...." end if st.T and not (st.p or st.h or st.s or st.v or st.q ~= nil) then return "Point '" .. role .. "' contains only T. Add one of p=..., h=..., s=..., v=... or Q=...." end if st.h and not (st.p or st.T or st.s) then return "Point '" .. role .. "' contains only h. Add p=..., T=... or s=...." end if st.s and not (st.p or st.T or st.h) then return "Point '" .. role .. "' contains only s. Add p=..., T=... or h=...." end if st.v and not (st.p or st.T) then return "Point '" .. role .. "' contains only v. Add p=... or T=...." end return "Point '" .. role .. "' must contain two independent properties, for instance p=1bar,Q=1 or p=1bar,T=20C." end local function lcp_process_error(msg, suggestions) local out = {"luacoolprop: " .. msg} if suggestions and #suggestions > 0 then out[#out + 1] = "Suggestions:" for _, s in ipairs(suggestions) do out[#out + 1] = " - " .. s end end error(table.concat(out, "\n"), 3) end -- Tolerances are deliberately tighter than plotting precision but loose enough -- to absorb normal flash-solver roundoff when the same state is expressed by -- redundant input properties. The comparison uses the larger of the absolute -- floor and the relative term. M._strict.state_property_tolerances = { p = {absolute = 0.1, relative = 5e-7, name = "pressure p", unit = "Pa"}, T = {absolute = 1e-5, relative = 5e-7, name = "temperature T", unit = "K"}, h = {absolute = 1e-3, relative = 5e-7, name = "specific enthalpy h", unit = "J/kg"}, s = {absolute = 1e-5, relative = 5e-7, name = "specific entropy s", unit = "J/(kg K)"}, q = {absolute = 1e-9, relative = 5e-7, name = "vapour quality Q", unit = "1"}, v = {absolute = 1e-12, relative = 5e-7, name = "specific volume v", unit = "m3/kg"}, rho = {absolute = 1e-8, relative = 5e-7, name = "mass density rho", unit = "kg/m3"}, } function M._strict.property_difference(key, expected, actual) local specification = assert(M._strict.state_property_tolerances[key], "missing thermodynamic comparison tolerance") local difference = math.abs(expected - actual) local tolerance = math.max(specification.absolute, specification.relative * math.max(1, math.abs(expected), math.abs(actual))) return difference, tolerance, specification end function M._strict.validate_property_consistency(key, expected, actual, context) if not finite_number(actual) then lcp_process_error("CoolProp could not resolve " .. M._strict.state_property_tolerances[key].name .. " for " .. context, { "Check that the selected independent property pair is valid for this fluid.", }) end local difference, tolerance, specification = M._strict.property_difference(key, expected, actual) if difference > tolerance then lcp_process_error("inconsistent " .. specification.name .. " for " .. context, { string.format("Supplied value: %.12g %s.", expected, specification.unit), string.format("CoolProp value: %.12g %s.", actual, specification.unit), string.format("Absolute difference %.12g exceeds tolerance %.12g.", difference, tolerance), "Remove the redundant property or correct it so all properties describe one state.", }) end end local function validate_state_values(st, role) local suggestions = {} if st.p and st.p <= 0 then suggestions[#suggestions + 1] = "Pressure in point '" .. role .. "' must be strictly positive." end if st.T and st.T <= 0 then suggestions[#suggestions + 1] = "Temperature in point '" .. role .. "' must be in kelvin or given with C, for example T=20C." end if st.v and st.v <= 0 then suggestions[#suggestions + 1] = "Specific volume in point '" .. role .. "' must be strictly positive." end if st.rho and st.rho <= 0 then suggestions[#suggestions + 1] = "Density in point '" .. role .. "' must be strictly positive." end if st.q ~= nil and (st.q < 0 or st.q > 1) then suggestions[#suggestions + 1] = "Quality Q in point '" .. role .. "' must lie between 0 and 1." end if #suggestions > 0 then lcp_process_error("invalid state specification for point '" .. role .. "'", suggestions) end end local function property_name(kind) local names = {isobar = "pressure p", isotherm = "temperature T", isentropic = "entropy s", isenthalpic = "enthalpy h", isoquality = "quality Q", isochore = "specific volume v"} return names[kind] or "process constant" end local function process_constant_key(kind) return ({isobar = "p", isotherm = "T", isentropic = "s", isenthalpic = "h", isoquality = "q", isochore = "v"})[kind] end local function process_constant_kind(kind) return ({isobar = "pressure", isotherm = "temperature", isentropic = "entropy", isenthalpic = "enthalpy", isoquality = "quality", isochore = "specific_volume"})[kind] end local function explicit_constant_keys(kind) if kind == "isobar" then return {"pressure", "p", "value"}, "pressure=..." end if kind == "isotherm" then return {"temperature", "T", "t", "value"}, "temperature=..." end if kind == "isentropic" then return {"entropy", "s", "value"}, "entropy=..." end if kind == "isenthalpic" then return {"enthalpy", "h", "value"}, "enthalpy=..." end if kind == "isoquality" then return {"quality", "q", "Q", "value"}, "quality=..." end if kind == "isochore" then return {"specific_volume", "specific volume", "v", "value"}, "specific volume=..." end return {"value"}, "value=..." end local function complete_state(fluid, st, h_scale, p_scale, role, diagram_type, v_scale, s_scale, temperature_scale, temperature_offset) role = role or "point" if st.rho and not st.v and st.rho > 0 then st.v = 1 / st.rho end local n1, v1, n2, v2 = state_pair(st) if not n1 then lcp_process_error("cannot define point '" .. role .. "' from " .. state_short(st), {state_pair_hint(st, role)}) end -- Freeze the defining pair before adding derived fields. Re-selecting a -- pair after computing pressure can turn a valid two-phase T-rho input into -- the ambiguous saturation pair P-T. Resolve every property from this one -- pair, validate all explicitly supplied redundant properties, and finally -- replace the user spellings with CoolProp's canonical state values. local function original_pair_output(output) return M._strict.coolprop_output_from_pair(fluid, output, n1, v1, n2, v2) end local resolved = { p = original_pair_output("p"), T = original_pair_output("t"), h = original_pair_output("h"), s = original_pair_output("s"), v = original_pair_output("v"), rho = original_pair_output("rho"), } local quality = original_pair_output("q") if quality and quality >= 0 and quality <= 1 then resolved.q = quality end for key in pairs(st._props or {}) do if M._strict.state_property_tolerances[key] then M._strict.validate_property_consistency(key, st[key], resolved[key], "point '" .. role .. "'") end end st.p, st.T, st.h, st.s = resolved.p, resolved.T, resolved.h, resolved.s st.v, st.rho, st.q = resolved.v, resolved.rho, resolved.q if not st.p or not st.h or (diagram_type == "PV" and not st.v) or ((diagram_type == "TS" or diagram_type == "HS") and (not st.s or (diagram_type == "TS" and not st.T))) or (diagram_type == "PT" and not st.T) then lcp_process_error("CoolProp could not resolve point '" .. role .. "' from " .. state_short(st), { "Check that the chosen pair of properties is physically valid for the fluid.", "Check units: p defaults to Pa, T to K, h to J/kg, s to J/(kg K), v to m^3/kg." }) end if diagram_type == "PV" then st.x = st.v * (v_scale or 1) st.y = st.p * p_scale elseif diagram_type == "TS" then st.x = st.s * (s_scale or 1e-3) st.y = (st.T - (temperature_offset or 0)) * (temperature_scale or 1) elseif diagram_type == "HS" then st.x = st.s * (s_scale or 1e-3) st.y = st.h * h_scale elseif diagram_type == "PT" then st.x = (st.T - (temperature_offset or 0)) * (temperature_scale or 1) st.y = st.p * p_scale else st.x = st.h * h_scale st.y = st.p * p_scale end st.t = log10(st.p) return st end local function parse_process_constant(opts, names, kind) for _, name in ipairs(names) do local v = opts[name] if v ~= nil and v ~= "" and v ~= "auto" then return parse_quantity_value(v, kind) end end return nil end local function process_kind(opts) local raw = to_string(opts.type or opts.kind or opts.process, "isobar") local k = raw:lower():gsub("[%s%-_]+", "") local aliases = { isobar = "isobar", isobaric = "isobar", p = "isobar", pressure = "isobar", isotherm = "isotherm", isothermal = "isotherm", temperature = "isotherm", t = "isotherm", isentrope = "isentropic", isentropic = "isentropic", entropy = "isentropic", s = "isentropic", isenthalp = "isenthalpic", isenthalpic = "isenthalpic", enthalpy = "isenthalpic", h = "isenthalpic", throttle = "isenthalpic", isoquality = "isoquality", quality = "isoquality", q = "isoquality", isochore = "isochore", isochoric = "isochore", volume = "isochore", v = "isochore", } local kind = aliases[k] if not kind then lcp_process_error("unknown process type '" .. raw .. "'", { "Use one of: isobar, isotherm, isentropic, isenthalpic, quality, isochore.", "Aliases such as p, T, s, h, q and v are also accepted." }) end return kind end local function process_state_output(fluid, st, key) if key == "T" then return state_output(fluid, st, "t") end return state_output(fluid, st, key) end local function process_constant_from_states(kind, fluid, from, to, opts) local names, example = explicit_constant_keys(kind) local ckind = process_constant_kind(kind) local key = process_constant_key(kind) local explicit = parse_process_constant(opts, names, ckind) if explicit ~= nil then return explicit end if key and from[key] ~= nil then return from[key] end if key and to[key] ~= nil then return to[key] end local fv = key and process_state_output(fluid, from, key) or nil if fv ~= nil then return fv end local tv = key and process_state_output(fluid, to, key) or nil if tv ~= nil then return tv end local suggestions = { "Give the conserved property explicitly with " .. example .. ", or value=....", "Alternatively, make point 'from' or point 'to' sufficiently defined so that LuaCoolProp can compute " .. property_name(kind) .. "." } if not state_pair(from) then suggestions[#suggestions + 1] = state_pair_hint(from, "from") end if not state_pair(to) then suggestions[#suggestions + 1] = state_pair_hint(to, "to") end lcp_process_error("missing " .. property_name(kind) .. " for process type '" .. kind .. "'", suggestions) end local function validate_process_constant(kind, cval) local suggestions = {} if cval == nil or not finite_number(cval) then suggestions[#suggestions + 1] = "The conserved property must be a finite number." elseif kind == "isobar" and cval <= 0 then suggestions[#suggestions + 1] = "Pressure must be strictly positive. Example: p=10bar." elseif kind == "isotherm" and cval <= 0 then suggestions[#suggestions + 1] = "Temperature must be positive in kelvin, or use a Celsius unit. Example: T=20C." elseif kind == "isoquality" and (cval < 0 or cval > 1) then suggestions[#suggestions + 1] = "Quality Q must lie between 0 and 1." elseif kind == "isochore" and cval <= 0 then suggestions[#suggestions + 1] = "Specific volume must be strictly positive. Example: v=0.02m3kg." end if #suggestions > 0 then lcp_process_error("invalid conserved property for process type '" .. kind .. "'", suggestions) end end local function apply_process_constant(kind, st, cval) if kind == "isobar" and not st.p then st.p = cval end if kind == "isotherm" and not st.T then st.T = cval end if kind == "isentropic" and not st.s then st.s = cval end if kind == "isenthalpic" and not st.h then st.h = cval end if kind == "isoquality" and st.q == nil then st.q = cval end if kind == "isochore" and not st.v then st.v = cval end end function M._strict.validate_process_conservation(kind, cval, from, to) local key_by_kind = { isobar = "p", isotherm = "T", isentropic = "s", isenthalpic = "h", isoquality = "q", isochore = "v", } local key = key_by_kind[kind] if not key then return end M._strict.validate_property_consistency(key, cval, from[key], "process 'from' endpoint") M._strict.validate_property_consistency(key, cval, to[key], "process 'to' endpoint") end local function validate_process_endpoint(fluid, st, role) validate_state_values(st, role) if not state_pair(st) then lcp_process_error("cannot define point '" .. role .. "' from " .. state_short(st), {state_pair_hint(st, role)}) end end local function process_id(opts, kind, cval, from, to) if opts.name and opts.name ~= "" and opts.name ~= "auto" then return opts.name end local diagram_type = to_string(opts.diagram_type, "PH"):lower() local suffix = curve_id("lcp-" .. diagram_type .. "-process", kind, cval or 0, nil, 4) local fp = from and from.p and ("p" .. slug_number(from.p * 1e-5, 3) .. "bar") or "from" local tp = to and to.p and ("p" .. slug_number(to.p * 1e-5, 3) .. "bar") or "to" return suffix .. "-" .. fp .. "-to-" .. tp end local function reverse_points(points) -- Reconstruct the topology from the current flat array rather than trusting -- `points.segments`: a specialised sampler may have refined or extended the -- array after the adaptive pass. This keeps every late critical-point -- sample and reverses each connected component independently. local forward, current = {}, nil for _, point in ipairs(points) do if point._break_before or not current then current = {} forward[#forward + 1] = current end current[#current + 1] = point end local out = {segments = {}} for segment_index = #forward, 1, -1 do local segment = {} for point_index = #forward[segment_index], 1, -1 do local copy = {} for key, value in pairs(forward[segment_index][point_index]) do if key ~= "_break_before" then copy[key] = value end end segment[#segment + 1] = copy end if #out.segments > 0 then segment[1]._break_before = true end out.segments[#out.segments + 1] = segment for _, point in ipairs(segment) do out[#out + 1] = point end end out.domain_gap_count = math.max(0, #out.segments - 1) out.omitted_sample_count = points.omitted_sample_count out.partially_omitted = points.partially_omitted return out end local function two_point_curve(a, b) local function copy_point(state) return { x = state.x, y = state.y, p = state.p, h = state.h, T = state.T, s = state.s, q = state.q, v = state.v, rho = state.rho, t = state.t, } end return {copy_point(a), copy_point(b)} end --- Resolve and sample a thermodynamic process segment in a PH diagram. -- Supported process types are `isobar`, `isotherm`, `isentropic`, -- `isenthalpic`, `quality`, and `isochore`; `isoquality` is a compatibility -- alias for `quality`. -- @tparam table opts Process options, including `fluid`, `type`, `from`, and `to`. -- @treturn table Array of plot points. -- @treturn table Metadata (`id`, `kind`, conserved value, endpoints). function M.diagram.ph_process_points(opts) opts = opts or {} M.diagram.prepare_fluid(opts) M.load_library(opts.library) local fluid = opts.fluid or "R134a" local h_scale = to_number(opts.h_scale, 1e-3) local p_scale = to_number(opts.p_scale, 1e-5) local kind = process_kind(opts) local from = parse_state_spec(opts.from or opts.start or "") local to = parse_state_spec(opts.to or opts.stop or opts.finish or "") validate_state_values(from, "from") validate_state_values(to, "to") local cval = process_constant_from_states(kind, fluid, from, to, opts) validate_process_constant(kind, cval) apply_process_constant(kind, from, cval) apply_process_constant(kind, to, cval) validate_process_endpoint(fluid, from, "from") validate_process_endpoint(fluid, to, "to") from = complete_state(fluid, from, h_scale, p_scale, "from") to = complete_state(fluid, to, h_scale, p_scale, "to") M._strict.validate_process_conservation(kind, cval, from, to) local points if kind == "isobar" or kind == "isenthalpic" then points = two_point_curve(from, to) else local pmin, pmax = math.min(from.p, to.p), math.max(from.p, to.p) if math.abs(pmax - pmin) <= 1e-9 * math.max(1, pmax) then points = two_point_curve(from, to) else local gen = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, h_scale = h_scale, p_scale = p_scale, initial_intervals = to_number(opts.process_initial_intervals or opts.initial_intervals, 12), max_depth = to_number(opts.process_max_depth or opts.max_depth, 7), tolerance = to_number(opts.process_tolerance or opts.tolerance, 0.25), log_weight = to_number(opts.log_weight, 30), constants = M.fluid_constants(fluid, opts.library), domain_policy = opts.domain_policy or opts.discontinuity_policy, } if kind == "isotherm" then points = M.diagram.ph_isotherm_curve(gen, cval) elseif kind == "isentropic" then points = M.diagram.ph_isentrope_curve(gen, cval) elseif kind == "isoquality" then points = M.diagram.ph_isoquality_curve(gen, cval) elseif kind == "isochore" then points = M.diagram.ph_isochore_curve(gen, cval) else points = two_point_curve(from, to) end if from.p > to.p then points = reverse_points(points) end end end local id_opts = copy_table(opts) id_opts.diagram_type = "PH" return points, {kind = kind, constant = cval, from = from, to = to, id = process_id(id_opts, kind, cval, from, to), diagram_type = "PH", reference_state = M.reference_state(fluid)} end local function render_process_plot(points, meta, opts) local coord_digits = math.floor(to_number(opts.coord_digits, 6)) local color = opts.process_color or opts.color or "black" local style = join_style(plot_name_style(meta.id), color, opts.process_style or opts.style or "line width=1pt", "unbounded coords=jump") local texout = {} texout[#texout + 1] = "\\addplot[" .. style .. "] coordinates {" local coords = {} for _, point in ipairs(points) do append_coord(coords, point, coord_digits) end texout[#texout + 1] = table.concat(coords, " ") texout[#texout + 1] = "}; " if to_bool(opts.mark_endpoints or opts.markers or opts.marks, false) and #points >= 2 then local marker_style = join_style(opts.marker_style or "only marks,mark=*,mark size=1.2pt", opts.marker_color or color) texout[#texout + 1] = "\\addplot[" .. marker_style .. ",forget plot] coordinates {" .. "(" .. fmt_num(points[1].x, coord_digits) .. "," .. fmt_num(points[1].y, coord_digits) .. ") " .. "(" .. fmt_num(points[#points].x, coord_digits) .. "," .. fmt_num(points[#points].y, coord_digits) .. ")}; " end local label = opts.label or opts.label_text if label and label ~= "" and label ~= "none" then add_curve_label(texout, meta.id, opts.label_style or "luacoolprop process label node", points, to_number(opts.label_pos, 0.5), label, coord_digits, to_bool(opts.label_sloped or opts.sloped_label or opts.sloped, true), to_bool(opts.label_allow_upside_down, false), meta.diagram_type == "PV", meta.diagram_type == "PH" or meta.diagram_type == "PV" or meta.diagram_type == "PT") end return table.concat(texout), meta end --- Emit PGFPlots code for a thermodynamic process segment. -- @tparam table opts Process options. -- @treturn string PGFPlots/TikZ code. -- @treturn table Process metadata, including completed endpoint coordinates. function M.diagram.ph_process_plot(opts) opts = opts or {} M.load_library(opts.library) local points, meta = M.diagram.ph_process_points(opts) return render_process_plot(points, meta, opts) end --- Resolve and sample a thermodynamic process segment in a PV diagram. -- Isobars and isochores are exact straight segments. Isotherms, isentropes, -- qualities, and isenthalps use the same adaptive logarithmic sampling as the -- corresponding PV curves. -- @tparam table opts Process options, including `fluid`, `type`, `from`, `to`, -- `v_scale`, and `p_scale`. -- @treturn table Array of PV plot points. -- @treturn table Process metadata with completed endpoints. function M.diagram.pv_process_points(opts) opts = opts or {} M.diagram.prepare_fluid(opts) M.load_library(opts.library) local fluid = opts.fluid or "R134a" local h_scale = to_number(opts.h_scale, 1e-3) local v_scale = to_number(opts.v_scale, 1) local p_scale = to_number(opts.p_scale, 1e-5) local kind = process_kind(opts) local from = parse_state_spec(opts.from or opts.start or "") local to = parse_state_spec(opts.to or opts.stop or opts.finish or "") validate_state_values(from, "from") validate_state_values(to, "to") local constant = process_constant_from_states(kind, fluid, from, to, opts) validate_process_constant(kind, constant) apply_process_constant(kind, from, constant) apply_process_constant(kind, to, constant) validate_process_endpoint(fluid, from, "from") validate_process_endpoint(fluid, to, "to") from = complete_state(fluid, from, h_scale, p_scale, "from", "PV", v_scale) to = complete_state(fluid, to, h_scale, p_scale, "to", "PV", v_scale) M._strict.validate_process_conservation(kind, constant, from, to) local points if kind == "isobar" or kind == "isochore" then points = two_point_curve(from, to) else local pmin, pmax = math.min(from.p, to.p), math.max(from.p, to.p) if math.abs(pmax - pmin) <= 1e-9 * math.max(1, pmax) then points = two_point_curve(from, to) else local generator = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, v_scale = v_scale, p_scale = p_scale, initial_intervals = to_number( opts.process_initial_intervals or opts.initial_intervals, 12), max_depth = to_number(opts.process_max_depth or opts.max_depth, 7), tolerance = to_number(opts.process_tolerance or opts.tolerance, 0.18), log_weight = to_number(opts.log_weight, 30), log_x_weight = to_number(opts.log_x_weight, 30), constants = M.fluid_constants(fluid, opts.library), domain_policy = opts.domain_policy or opts.discontinuity_policy, } if kind == "isotherm" then points = M.diagram.pv_isotherm_curve(generator, constant) elseif kind == "isentropic" then points = M.diagram.pv_isentrope_curve(generator, constant) elseif kind == "isoquality" then points = M.diagram.pv_quality_curve(generator, constant) elseif kind == "isenthalpic" then points = pv_isenthalp_curve(generator, constant) else points = two_point_curve(from, to) end if from.p > to.p then points = reverse_points(points) end end end local id_opts = copy_table(opts) id_opts.diagram_type = "PV" return points, { kind = kind, constant = constant, from = from, to = to, id = process_id(id_opts, kind, constant, from, to), diagram_type = "PV", reference_state = M.reference_state(fluid), } end --- Emit PGFPlots code for a thermodynamic process in a PV diagram. -- @tparam table opts Process options. -- @treturn string PGFPlots/TikZ code. -- @treturn table Completed process metadata. function M.diagram.pv_process_plot(opts) opts = opts or {} M.load_library(opts.library) local points, metadata = M.diagram.pv_process_points(opts) return render_process_plot(points, metadata, opts) end --- Resolve and sample a thermodynamic process segment in a TS diagram. -- Isotherms and isentropes are exact straight lines in TS coordinates. -- Isobars are sampled with entropy as the independent variable so their -- two-phase plateau is preserved; isochores use temperature, while quality -- curves and isenthalps use adaptive logarithmic-pressure sampling. -- @tparam table opts Process options, including `fluid`, `type`, `from`, -- `to`, `s_scale`, and `temperature_scale`. -- @treturn table Array of TS plot points. -- @treturn table Process metadata with completed endpoints. function M.diagram.ts_process_points(opts) opts = opts or {} M.diagram.prepare_fluid(opts) M.load_library(opts.library) local fluid = opts.fluid or "R134a" local h_scale = to_number(opts.h_scale, 1e-3) local p_scale = to_number(opts.p_scale, 1e-5) local s_scale = to_number(opts.s_scale, 1e-3) local temperature_scale = to_number(opts.temperature_scale or opts.t_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local kind = process_kind(opts) local from = parse_state_spec(opts.from or opts.start or "") local to = parse_state_spec(opts.to or opts.stop or opts.finish or "") validate_state_values(from, "from") validate_state_values(to, "to") local constant = process_constant_from_states(kind, fluid, from, to, opts) validate_process_constant(kind, constant) apply_process_constant(kind, from, constant) apply_process_constant(kind, to, constant) validate_process_endpoint(fluid, from, "from") validate_process_endpoint(fluid, to, "to") from = complete_state(fluid, from, h_scale, p_scale, "from", "TS", nil, s_scale, temperature_scale, temperature_offset) to = complete_state(fluid, to, h_scale, p_scale, "to", "TS", nil, s_scale, temperature_scale, temperature_offset) M._strict.validate_process_conservation(kind, constant, from, to) local sampling = { initial_intervals = to_number( opts.process_initial_intervals or opts.initial_intervals, 12), max_depth = to_number(opts.process_max_depth or opts.max_depth, 7), tolerance = to_number(opts.process_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 8), log_y = false, y_weight = to_number(opts.temperature_weight, 0.01), } M._strict.domain_options(opts, sampling, "TS process") local points if kind == "isotherm" or kind == "isentropic" then points = two_point_curve(from, to) elseif kind == "isobar" then if math.abs(to.s - from.s) <= 1e-12 * math.max(1, math.abs(from.s)) then points = two_point_curve(from, to) else points = adaptive_logp_curve(from.s, to.s, function(entropy) return ts_point_from_pair(fluid, "P", constant, "S", entropy, entropy, s_scale, temperature_scale, temperature_offset) end, sampling) end elseif kind == "isochore" then local density = 1 / constant if math.abs(to.T - from.T) <= 1e-12 * math.max(1, from.T) then points = two_point_curve(from, to) else points = adaptive_logp_curve(from.T, to.T, function(temperature) return ts_point_from_pair(fluid, "T", temperature, "Dmass", density, temperature, s_scale, temperature_scale, temperature_offset) end, sampling) end else local pmin, pmax = math.min(from.p, to.p), math.max(from.p, to.p) if math.abs(pmax - pmin) <= 1e-9 * math.max(1, pmax) then points = two_point_curve(from, to) else local generator = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, s_scale = s_scale, temperature_scale = temperature_scale, temperature_offset = temperature_offset, initial_intervals = sampling.initial_intervals, max_depth = sampling.max_depth, tolerance = sampling.tolerance, log_weight = sampling.log_weight, temperature_weight = sampling.y_weight, constants = M.fluid_constants(fluid, opts.library), domain_policy = opts.domain_policy or opts.discontinuity_policy, } if kind == "isenthalpic" then points = M.diagram.ts_isenthalp_curve(generator, constant) elseif kind == "isoquality" then points = M.diagram.ts_quality_curve(generator, constant) else points = two_point_curve(from, to) end if from.p > to.p then points = reverse_points(points) end end end local id_opts = copy_table(opts) id_opts.diagram_type = "TS" return points, { kind = kind, constant = constant, from = from, to = to, id = process_id(id_opts, kind, constant, from, to), diagram_type = "TS", reference_state = M.reference_state(fluid), } end --- Emit PGFPlots code for a thermodynamic process in a TS diagram. -- @tparam table opts TS process options. -- @treturn string PGFPlots/TikZ code. -- @treturn table Completed process metadata. function M.diagram.ts_process_plot(opts) opts = opts or {} M.load_library(opts.library) local points, metadata = M.diagram.ts_process_points(opts) return render_process_plot(points, metadata, opts) end --- Resolve and sample a thermodynamic process in an HS diagram. -- Isentropes and isenthalps are exact vertical and horizontal segments. -- Other conserved-property paths use the same physical samplers as HS -- background families and retain complete endpoint metadata. -- @tparam table opts HS process options. -- @treturn table HS points. -- @treturn table Process metadata. function M.diagram.hs_process_points(opts) opts = opts or {} M.diagram.prepare_fluid(opts) M.load_library(opts.library) local fluid = opts.fluid or "R134a" local h_scale = to_number(opts.h_scale, 1e-3) local p_scale = to_number(opts.p_scale, 1e-5) local s_scale = to_number(opts.s_scale, 1e-3) local kind = process_kind(opts) local from = parse_state_spec(opts.from or opts.start or "") local to = parse_state_spec(opts.to or opts.stop or opts.finish or "") validate_state_values(from, "from") validate_state_values(to, "to") local constant = process_constant_from_states(kind, fluid, from, to, opts) validate_process_constant(kind, constant) apply_process_constant(kind, from, constant) apply_process_constant(kind, to, constant) validate_process_endpoint(fluid, from, "from") validate_process_endpoint(fluid, to, "to") from = complete_state(fluid, from, h_scale, p_scale, "from", "HS", nil, s_scale) to = complete_state(fluid, to, h_scale, p_scale, "to", "HS", nil, s_scale) M._strict.validate_process_conservation(kind, constant, from, to) local sampling = { initial_intervals = to_number( opts.process_initial_intervals or opts.initial_intervals, 12), max_depth = to_number(opts.process_max_depth or opts.max_depth, 7), tolerance = to_number(opts.process_tolerance or opts.tolerance, 0.10), log_weight = 1, log_y = false, y_weight = to_number(opts.enthalpy_weight, 0.01), } M._strict.domain_options(opts, sampling, "HS process") local points if kind == "isentropic" or kind == "isenthalpic" then points = two_point_curve(from, to) elseif kind == "isobar" then points = adaptive_logp_curve(from.s, to.s, function(entropy) return hs_point_from_pair(fluid, "P", constant, "S", entropy, entropy, s_scale, h_scale) end, sampling) else local pmin, pmax = math.min(from.p, to.p), math.max(from.p, to.p) if math.abs(pmax - pmin) <= 1e-9 * math.max(1, pmax) then points = two_point_curve(from, to) else local generator = { fluid = fluid, library = opts.library, pressure_min = pmin, pressure_max = pmax, s_scale = s_scale, h_scale = h_scale, initial_intervals = sampling.initial_intervals, max_depth = sampling.max_depth, tolerance = sampling.tolerance, enthalpy_weight = sampling.y_weight, constants = M.fluid_constants(fluid, opts.library), domain_policy = opts.domain_policy or opts.discontinuity_policy, } if kind == "isotherm" then points = M.diagram.hs_isotherm_curve(generator, constant) elseif kind == "isoquality" then points = M.diagram.hs_quality_curve(generator, constant) elseif kind == "isochore" then points = M.diagram.hs_isochore_curve(generator, constant) else points = two_point_curve(from, to) end if from.p > to.p then points = reverse_points(points) end end end local id_opts = copy_table(opts) id_opts.diagram_type = "HS" return points, { kind = kind, constant = constant, from = from, to = to, id = process_id(id_opts, kind, constant, from, to), diagram_type = "HS", reference_state = M.reference_state(fluid), } end --- Serialize one HS process path. -- @tparam table opts HS process options. -- @treturn string PGFPlots code. -- @treturn table Process metadata. function M.diagram.hs_process_plot(opts) opts = opts or {} M.load_library(opts.library) local points, metadata = M.diagram.hs_process_points(opts) return render_process_plot(points, metadata, opts) end --- Resolve and sample a thermodynamic process in a PT diagram. -- Isobars and isotherms are exact horizontal and vertical segments. The -- remaining conserved-property paths are sampled in logarithmic pressure; -- a constant-quality path lies on the liquid--vapour equilibrium curve -- because PT coordinates do not distinguish two-phase qualities. -- @tparam table opts PT process options. -- @treturn table PT points. -- @treturn table Process metadata with completed endpoints. function M.diagram.pt_process_points(opts) opts = opts or {} M.diagram.prepare_fluid(opts) M.load_library(opts.library) local fluid = opts.fluid or "R134a" local h_scale = to_number(opts.h_scale, 1e-3) local p_scale = to_number(opts.p_scale, 1e-5) local temperature_scale = to_number(opts.temperature_scale, 1) local temperature_offset = to_number(opts.temperature_offset, 0) local kind = process_kind(opts) local from = parse_state_spec(opts.from or opts.start or "") local to = parse_state_spec(opts.to or opts.stop or opts.finish or "") validate_state_values(from, "from") validate_state_values(to, "to") local constant = process_constant_from_states(kind, fluid, from, to, opts) validate_process_constant(kind, constant) apply_process_constant(kind, from, constant) apply_process_constant(kind, to, constant) validate_process_endpoint(fluid, from, "from") validate_process_endpoint(fluid, to, "to") from = complete_state(fluid, from, h_scale, p_scale, "from", "PT", nil, nil, temperature_scale, temperature_offset) to = complete_state(fluid, to, h_scale, p_scale, "to", "PT", nil, nil, temperature_scale, temperature_offset) M._strict.validate_process_conservation(kind, constant, from, to) local points if kind == "isobar" or kind == "isotherm" then points = two_point_curve(from, to) else local pmin, pmax = math.min(from.p, to.p), math.max(from.p, to.p) if math.abs(pmax - pmin) <= 1e-9 * math.max(1, pmax) then points = two_point_curve(from, to) else local generator = { fluid = fluid, library = opts.library, reference_state = opts.reference_state, pressure_min = pmin, pressure_max = pmax, temperature_scale = temperature_scale, temperature_offset = temperature_offset, p_scale = p_scale, initial_intervals = to_number( opts.process_initial_intervals or opts.initial_intervals, 12), max_depth = to_number(opts.process_max_depth or opts.max_depth, 7), tolerance = to_number(opts.process_tolerance or opts.tolerance, 0.10), log_weight = to_number(opts.log_weight, 30), constants = M.fluid_constants(fluid, opts.library), domain_policy = opts.domain_policy or opts.discontinuity_policy, } if kind == "isentropic" then points = M.diagram.pt_isentrope_curve(generator, constant) elseif kind == "isenthalpic" then points = M.diagram.pt_isenthalp_curve(generator, constant) elseif kind == "isochore" then points = M.diagram.pt_isochore_curve(generator, constant) elseif kind == "isoquality" then points = M.diagram._pt_constant_property_curve(generator, "Q", constant, "quality") else points = two_point_curve(from, to) end if from.p > to.p then points = reverse_points(points) end end end local id_options = copy_table(opts) id_options.diagram_type = "PT" return points, { kind = kind, constant = constant, from = from, to = to, id = process_id(id_options, kind, constant, from, to), diagram_type = "PT", reference_state = M.reference_state(fluid), } end --- Serialize one PT process path. function M.diagram.pt_process_plot(opts) opts = opts or {} M.load_library(opts.library) local points, metadata = M.diagram.pt_process_points(opts) return render_process_plot(points, metadata, opts) end -- Validate a user-supplied control-sequence stem before putting it in a -- \csname. Keeping this deliberately conservative makes the TeX bridge -- predictable and prevents process options from injecting arbitrary TeX. local function process_export_prefix(opts) local prefix = trim(opts.export_coordinates or "") if prefix == "" then return nil end if not prefix:match("^[A-Za-z@:_][A-Za-z0-9@:_]*$") then error("luacoolprop: 'export coordinates' must be a TeX control-sequence " .. "prefix containing only letters, digits, @, colon, or underscore", 3) end return prefix end -- Serialize process endpoints as global TeX definitions. The values use the -- same precision as the PGFPlots coordinates, so a published value always -- describes the path that TeX actually receives rather than a separately -- rounded thermodynamic result. local function process_coordinate_definitions(metadata, opts) local prefix = process_export_prefix(opts) if not prefix then return "" end local digits = math.floor(to_number(opts.coord_digits, 6)) local definitions = { {"FromX", metadata.from.x}, {"FromY", metadata.from.y}, {"ToX", metadata.to.x}, {"ToY", metadata.to.y}, } local state_properties = { {"PressureSI", "p"}, {"TemperatureK", "T"}, {"EnthalpySI", "h"}, {"EntropySI", "s"}, {"SpecificVolumeSI", "v"}, {"Quality", "q"}, } for _, endpoint in ipairs({{"From", metadata.from}, {"To", metadata.to}}) do for _, property in ipairs(state_properties) do local value = endpoint[2][property[2]] if finite_number(value) then definitions[#definitions + 1] = {endpoint[1] .. property[1], value} end end end local output = {} for _, definition in ipairs(definitions) do output[#output + 1] = "\\expandafter\\gdef\\csname " .. prefix .. definition[1] .. "\\endcsname{" .. fmt_num(definition[2], digits) .. "}" end return table.concat(output) end -- Logging is opt-in because large diagrams can contain many process segments. -- Use LuaTeX's terminal-and-log channel so the record is visible interactively -- and remains available in the reproducible build log. local function log_process_coordinates(metadata, opts) if not to_bool(opts.log_coordinates, false) then return end local digits = math.floor(to_number(opts.coord_digits, 6)) local message = "LuaCoolProp process coordinates [" .. metadata.id .. "]: " .. "from=(" .. fmt_num(metadata.from.x, digits) .. "," .. fmt_num(metadata.from.y, digits) .. "), to=(" .. fmt_num(metadata.to.x, digits) .. "," .. fmt_num(metadata.to.y, digits) .. ")" if texio and texio.write_nl then texio.write_nl("term and log", message) end end local function print_process(renderer, opts) local output, metadata = renderer(opts) tex.sprint(output .. process_coordinate_definitions(metadata, opts)) log_process_coordinates(metadata, opts) end --- Print PH quality plots into the active TeX input stream. -- @tparam[opt={}] table opts PH quality options. function M.tex.print_ph_isoquality(opts) tex.sprint(M.diagram.ph_isoquality_plots(opts or {})) end --- Print PH isotherm plots into the active TeX input stream. -- @tparam[opt={}] table opts PH isotherm options. function M.tex.print_ph_isotherms(opts) tex.sprint(M.diagram.ph_isotherm_plots(opts or {})) end --- Print PH isentrope plots into the active TeX input stream. -- @tparam[opt={}] table opts PH isentrope options. function M.tex.print_ph_isentropes(opts) tex.sprint(M.diagram.ph_isentrope_plots(opts or {})) end --- Print PH isochore plots into the active TeX input stream. -- @tparam[opt={}] table opts PH isochore options. function M.tex.print_ph_isochores(opts) tex.sprint(M.diagram.ph_isochore_plots(opts or {})) end --- Print one PH family with autonode-compatible serialization. -- This function exists for the generic TeX layer. Application Lua code should -- call the structured diagram renderer and handle the returned string instead. -- @tparam[opt={}] table opts PH family and label options. -- @tparam string family Canonical family name. function M.tex.print_ph_family_autonode(opts, family) opts = opts or {} local family_opts = ph_family_options(opts, family) local records = collect_ph_records(family_opts) -- Automatic labels are emitted directly through pgfplots-autonode trailing -- commands. LuaCoolProp does not run a separate collision-avoidance pass. tex.sprint(render_diagram_curve_records_only(records, opts)) end --- Print all enabled PH plots into the active TeX input stream. -- @tparam[opt={}] table opts PH diagram options. function M.tex.print_ph_plots(opts) tex.sprint(M.diagram.ph_plots(opts or {})) end --- Print a PGFPlots axis-style update for a PH request. -- @tparam[opt={}] table opts PH axis options. function M.tex.print_ph_axis_style(opts) local style = M.diagram.ph_axis_style(opts or {}) tex.sprint("\\pgfplotsset{every axis/.append style={" .. style .. "}}") end --- Define TeX macros containing commonly used fluid constants. -- This bridge writes `\LCPFluidName` and the available critical/triple-point -- macros into the active TeX input stream. -- @tparam[opt={}] table opts Options including `fluid`, `library`, `p_scale`, -- and `coord_digits`. function M.tex.print_fluid_constants(opts) opts = opts or {} local fluid = opts.fluid or "R134a" M.diagram.prepare_fluid(opts) local c = M.fluid_constants(fluid, opts.library) local pscale = to_number(opts.p_scale, 1e-5) local digits = math.floor(to_number(opts.coord_digits, 10)) local out = {} out[#out + 1] = "\\def\\LCPFluidName{" .. fluid .. "}" out[#out + 1] = "\\def\\LCPFluidReferenceState{" .. tostring(M.reference_state(fluid) or "DEF") .. "}" if c.pcrit then out[#out + 1] = "\\def\\LCPFluidPcritSI{" .. fmt_num(c.pcrit, digits) .. "}" end if c.pcrit then out[#out + 1] = "\\def\\LCPFluidPcritBar{" .. fmt_num(c.pcrit * 1e-5, digits) .. "}" end if c.pcrit then out[#out + 1] = "\\def\\LCPFluidPcritPlot{" .. fmt_num(c.pcrit * pscale, digits) .. "}" end if c.ptriple then out[#out + 1] = "\\def\\LCPFluidPtripleSI{" .. fmt_num(c.ptriple, digits) .. "}" end if c.ptriple then out[#out + 1] = "\\def\\LCPFluidPtripleBar{" .. fmt_num(c.ptriple * 1e-5, digits) .. "}" end if c.Tcrit then out[#out + 1] = "\\def\\LCPFluidTcritK{" .. fmt_num(c.Tcrit, digits) .. "}" end if c.Tcrit then out[#out + 1] = "\\def\\LCPFluidTcritC{" .. fmt_num(c.Tcrit - 273.15, digits) .. "}" end if c.Ttriple then out[#out + 1] = "\\def\\LCPFluidTtripleK{" .. fmt_num(c.Ttriple, digits) .. "}" end if c.Ttriple then out[#out + 1] = "\\def\\LCPFluidTtripleC{" .. fmt_num(c.Ttriple - 273.15, digits) .. "}" end if c.rhocrit then out[#out + 1] = "\\def\\LCPFluidRhocritSI{" .. fmt_num(c.rhocrit, digits) .. "}" end tex.sprint(table.concat(out)) end --- Print one PH process plot into the active TeX input stream. -- @tparam[opt={}] table opts PH process options. function M.tex.print_ph_process(opts) opts = opts or {} print_process(M.diagram.ph_process_plot, opts) end -- Structured PH API. Flat ph_* entry points are direct aliases, while every -- registered diagram implements the same table contract. local function canonical_ph_options(opts) return M.diagram.normalize_options("PH", opts, { quality = "isoquality", isotherm = "isotherms", isentrope = "isentropes", isochore = "isochores", enthalpy_scale = "h_scale", pressure_scale = "p_scale", }) end local function ph_api(renderer) return function(opts, ...) return renderer(canonical_ph_options(opts), ...) end end local ph_quality_curve = ph_api(M.diagram.ph_isoquality_curve) local ph_quality_plots = ph_api(M.diagram.ph_isoquality_plots) local ph_isotherm_curve = ph_api(M.diagram.ph_isotherm_curve) local ph_isotherm_plots = ph_api(M.diagram.ph_isotherm_plots) local ph_isentrope_curve = ph_api(M.diagram.ph_isentrope_curve) local ph_isentrope_plots = ph_api(M.diagram.ph_isentrope_plots) local ph_isochore_curve = ph_api(M.diagram.ph_isochore_curve) local ph_isochore_plots = ph_api(M.diagram.ph_isochore_plots) --- Canonical PH diagram implementation. -- Renderers return strings or Lua tables; they never print directly to TeX. -- The `families` field maps canonical family names to `{curve, plots}` tables. -- @table M.diagram.ph -- @field axis_style Return data-dependent PGFPlots axis keys. -- @field plots Serialize all enabled curve families. -- @field family Serialize one family selected by canonical name. -- @field process_points Resolve a process to point and metadata tables. -- @field process_plot Resolve and serialize a process. -- @field quality_curve Generate one constant-quality polyline. -- @field isotherm_curve Generate one constant-temperature polyline. -- @field isentrope_curve Generate one constant-entropy polyline. -- @field isochore_curve Generate one constant-specific-volume polyline. -- @field quality_plots Select and serialize quality curves. -- @field isotherm_plots Select and serialize isotherms. -- @field isentrope_plots Select and serialize isentropes. -- @field isochore_plots Select and serialize isochores. -- @field families Canonical family registry. M.diagram.ph = { axis_style = ph_api(M.diagram.ph_axis_style), plots = ph_api(M.diagram.ph_plots), process_points = ph_api(M.diagram.ph_process_points), process_plot = ph_api(M.diagram.ph_process_plot), quality_curve = ph_quality_curve, quality_plots = ph_quality_plots, isotherm_curve = ph_isotherm_curve, isotherm_plots = ph_isotherm_plots, isentrope_curve = ph_isentrope_curve, isentrope_plots = ph_isentrope_plots, isochore_curve = ph_isochore_curve, isochore_plots = ph_isochore_plots, families = { quality = { curve = ph_quality_curve, plots = ph_quality_plots, }, isotherm = { curve = ph_isotherm_curve, plots = ph_isotherm_plots, }, isentrope = { curve = ph_isentrope_curve, plots = ph_isentrope_plots, }, isochore = { curve = ph_isochore_curve, plots = ph_isochore_plots, }, }, } --- Serialize one canonical PH family. -- @tparam[opt={}] table opts PH diagram and family options. -- @tparam string family One of `quality`, `isotherm`, `isentrope`, or `isochore`. -- @treturn string PGFPlots `\addplot` commands. function M.diagram.ph.family(opts, family) opts = canonical_ph_options(opts) local family_opts = ph_family_options(opts, family) local records = collect_ph_records(family_opts) return render_diagram_curve_records_only(records, opts) end M.diagram.register_type("PH", M.diagram.ph) local function canonical_pv_options(opts) return M.diagram.normalize_options("PV", opts, { quality = "isoquality", isotherm = "isotherms", isentrope = "isentropes", specific_volume_scale = "v_scale", pressure_scale = "p_scale", }) end local function pv_api(renderer) return function(opts, ...) return renderer(canonical_pv_options(opts), ...) end end local pv_quality_curve = pv_api(M.diagram.pv_quality_curve) local pv_quality_plots = pv_api(M.diagram.pv_quality_plots) local pv_isotherm_curve = pv_api(M.diagram.pv_isotherm_curve) local pv_isotherm_plots = pv_api(M.diagram.pv_isotherm_plots) local pv_isentrope_curve = pv_api(M.diagram.pv_isentrope_curve) local pv_isentrope_plots = pv_api(M.diagram.pv_isentrope_plots) --- Canonical PV diagram implementation. -- The coordinate convention is `x=v*v_scale` and `y=p*p_scale`; raw SI -- values remain available as `point.v` and `point.p`. -- @table M.diagram.pv -- @field axis_style Return data-dependent PGFPlots pressure keys. -- @field plots Serialize all enabled curve families. -- @field family Serialize one canonical family. -- @field process_points Resolve a process to point and metadata tables. -- @field process_plot Resolve and serialize a process. -- @field quality_curve Generate one quality curve. -- @field isotherm_curve Generate one isotherm. -- @field isentrope_curve Generate one isentrope. -- @field quality_plots Select and serialize quality curves. -- @field isotherm_plots Select and serialize isotherms. -- @field isentrope_plots Select and serialize isentropes. -- @field families Canonical family registry. M.diagram.pv = { axis_style = pv_api(M.diagram.pv_axis_style), plots = pv_api(M.diagram.pv_plots), family = pv_api(M.diagram.pv_family), process_points = pv_api(M.diagram.pv_process_points), process_plot = pv_api(M.diagram.pv_process_plot), quality_curve = pv_quality_curve, quality_plots = pv_quality_plots, isotherm_curve = pv_isotherm_curve, isotherm_plots = pv_isotherm_plots, isentrope_curve = pv_isentrope_curve, isentrope_plots = pv_isentrope_plots, families = { quality = {curve = pv_quality_curve, plots = pv_quality_plots}, isotherm = {curve = pv_isotherm_curve, plots = pv_isotherm_plots}, isentrope = {curve = pv_isentrope_curve, plots = pv_isentrope_plots}, }, } M.diagram.register_type("PV", M.diagram.pv) local function canonical_ts_options(opts) return M.diagram.normalize_options("TS", opts, { quality = "isoquality", isenthalp = "isenthalps", entropy_scale = "s_scale", temperature_scale = "temperature_scale", }) end local function ts_api(renderer) return function(opts, ...) return renderer(canonical_ts_options(opts), ...) end end local ts_quality_curve = ts_api(M.diagram.ts_quality_curve) local ts_quality_plots = ts_api(M.diagram.ts_quality_plots) local ts_isenthalp_curve = ts_api(M.diagram.ts_isenthalp_curve) local ts_isenthalp_plots = ts_api(M.diagram.ts_isenthalp_plots) --- Canonical TS diagram implementation. -- The coordinate convention is `x=s*s_scale` and -- `y=T*temperature_scale`; the corresponding raw SI values remain available -- as `point.s` and `point.T`. -- @table M.diagram.ts -- @field axis_style Return data-dependent linear PGFPlots limits. -- @field plots Serialize all enabled TS families. -- @field family Serialize one selected TS family. -- @field process_points Resolve a process to point and metadata tables. -- @field process_plot Resolve and serialize a process. -- @field quality_curve Generate one constant-quality curve. -- @field isenthalp_curve Generate one constant-enthalpy curve. -- @field quality_plots Select and serialize quality curves. -- @field isenthalp_plots Select and serialize isenthalps. -- @field families Canonical family registry. M.diagram.ts = { axis_style = ts_api(M.diagram.ts_axis_style), plots = ts_api(M.diagram.ts_plots), family = ts_api(M.diagram.ts_family), process_points = ts_api(M.diagram.ts_process_points), process_plot = ts_api(M.diagram.ts_process_plot), quality_curve = ts_quality_curve, quality_plots = ts_quality_plots, isenthalp_curve = ts_isenthalp_curve, isenthalp_plots = ts_isenthalp_plots, families = { quality = {curve = ts_quality_curve, plots = ts_quality_plots}, isenthalp = {curve = ts_isenthalp_curve, plots = ts_isenthalp_plots}, }, } M.diagram.register_type("TS", M.diagram.ts) local function canonical_hs_options(opts) return M.diagram.normalize_options("HS", opts, { quality = "isoquality", isochore = "isochores", isotherm = "isotherms", isobar = "isobars", entropy_scale = "s_scale", enthalpy_scale = "h_scale", }) end local function hs_api(renderer) return function(opts, ...) return renderer(canonical_hs_options(opts), ...) end end local hs_quality_curve = hs_api(M.diagram.hs_quality_curve) local hs_isochore_curve = hs_api(M.diagram.hs_isochore_curve) local hs_isotherm_curve = hs_api(M.diagram.hs_isotherm_curve) local hs_isobar_curve = hs_api(M.diagram.hs_isobar_curve) --- Canonical HS diagram implementation. -- Coordinates are `x=s*entropy_scale` and `y=h*enthalpy_scale`. -- @table M.diagram.hs M.diagram.hs = { axis_style = hs_api(M.diagram.hs_axis_style), plots = hs_api(M.diagram.hs_plots), family = hs_api(M.diagram.hs_family), process_points = hs_api(M.diagram.hs_process_points), process_plot = hs_api(M.diagram.hs_process_plot), quality_curve = hs_quality_curve, quality_plots = hs_api(M.diagram.hs_quality_plots), isochore_curve = hs_isochore_curve, isochore_plots = hs_api(M.diagram.hs_isochore_plots), isotherm_curve = hs_isotherm_curve, isotherm_plots = hs_api(M.diagram.hs_isotherm_plots), isobar_curve = hs_isobar_curve, isobar_plots = hs_api(M.diagram.hs_isobar_plots), families = { quality = {curve = hs_quality_curve, plots = hs_api(M.diagram.hs_quality_plots)}, isochore = {curve = hs_isochore_curve, plots = hs_api(M.diagram.hs_isochore_plots)}, isotherm = {curve = hs_isotherm_curve, plots = hs_api(M.diagram.hs_isotherm_plots)}, isobar = {curve = hs_isobar_curve, plots = hs_api(M.diagram.hs_isobar_plots)}, }, } M.diagram.register_type("HS", M.diagram.hs) function M.diagram._canonical_pt_options(opts) return M.diagram.normalize_options("PT", opts, { phase_envelopes = "phase_envelope", isentrope = "isentropes", isenthalp = "isenthalps", isochore = "isochores", temperature_scale = "temperature_scale", pressure_scale = "p_scale", }) end function M.diagram._pt_api(renderer) return function(opts, ...) return renderer(M.diagram._canonical_pt_options(opts), ...) end end --- Canonical PT diagram implementation. -- Coordinates are `x=(T-temperature_offset)*temperature_scale` and -- `y=p*pressure_scale`. The high-level axis is linear in temperature and -- logarithmic in pressure. -- @table M.diagram.pt -- @field phase_envelope_curve Generate the liquid--vapour coexistence curve. -- @field isentrope_curve Generate one constant-entropy curve. -- @field isenthalp_curve Generate one constant-enthalpy curve. -- @field isochore_curve Generate one constant-specific-volume curve. M.diagram.pt = { axis_style = M.diagram._pt_api(M.diagram.pt_axis_style), plots = M.diagram._pt_api(M.diagram.pt_plots), family = M.diagram._pt_api(M.diagram.pt_family), process_points = M.diagram._pt_api(M.diagram.pt_process_points), process_plot = M.diagram._pt_api(M.diagram.pt_process_plot), phase_envelope_curve = M.diagram._pt_api(M.diagram.pt_phase_envelope_curve), phase_envelope_plots = M.diagram._pt_api(M.diagram.pt_phase_envelope_plots), isentrope_curve = M.diagram._pt_api(M.diagram.pt_isentrope_curve), isentrope_plots = M.diagram._pt_api(M.diagram.pt_isentrope_plots), isenthalp_curve = M.diagram._pt_api(M.diagram.pt_isenthalp_curve), isenthalp_plots = M.diagram._pt_api(M.diagram.pt_isenthalp_plots), isochore_curve = M.diagram._pt_api(M.diagram.pt_isochore_curve), isochore_plots = M.diagram._pt_api(M.diagram.pt_isochore_plots), families = { phase_envelope = { curve = M.diagram._pt_api(M.diagram.pt_phase_envelope_curve), plots = M.diagram._pt_api(M.diagram.pt_phase_envelope_plots)}, isentrope = {curve = M.diagram._pt_api(M.diagram.pt_isentrope_curve), plots = M.diagram._pt_api(M.diagram.pt_isentrope_plots)}, isenthalp = {curve = M.diagram._pt_api(M.diagram.pt_isenthalp_curve), plots = M.diagram._pt_api(M.diagram.pt_isenthalp_plots)}, isochore = {curve = M.diagram._pt_api(M.diagram.pt_isochore_curve), plots = M.diagram._pt_api(M.diagram.pt_isochore_plots)}, }, } M.diagram.register_type("PT", M.diagram.pt) -- The high-level TS, HS and PT TeX diagrams ask for data-dependent axis -- bounds immediately before asking for their plots. Their axis-style call -- necessarily samples every enabled curve. Reuse those *same records* for -- the following plot call, preserving every coordinate, disconnected-domain -- marker, named path and autonode path fraction. Never keep the cache beyond -- one matching call, and never reuse it after a reference-state change. local pending_axis_records local function snapshot_options(value, visiting) local kind = type(value) if kind ~= "table" then if kind == "string" or kind == "number" or kind == "boolean" or kind == "nil" then return value, true end return nil, false end if visiting[value] then return nil, false end visiting[value] = true local snapshot = {} for key, entry in pairs(value) do local key_kind = type(key) if key_kind ~= "string" and key_kind ~= "number" then visiting[value] = nil return nil, false end local copied, supported = snapshot_options(entry, visiting) if not supported then visiting[value] = nil return nil, false end snapshot[key] = copied end visiting[value] = nil return snapshot, true end local function same_options(first, second) if type(first) ~= type(second) then return false end if type(first) ~= "table" then return first == second end for key, value in pairs(first) do if not same_options(value, second[key]) then return false end end for key in pairs(second) do if first[key] == nil then return false end end return true end --- Print one operation from a registered diagram implementation for TeX. -- Supported operations are `plots`, `axis_style`, `family`, and `process`. -- This is the only structured renderer that is expected to write directly to -- TeX. It requires the LuaTeX global `tex` table and is not suitable for a -- standalone `texlua` program; call `M.diagram.get_type` there instead. -- @tparam string diagram_type Canonical uppercase diagram code. -- @tparam string operation Renderer operation. -- @tparam table opts Diagram or process options. For a process, -- `export_coordinates` publishes endpoint macros and `log_coordinates` -- records the endpoints through LuaTeX's terminal-and-log channel. -- @tparam[opt] string family Canonical family name for the `family` operation. function M.tex.print_diagram(diagram_type, operation, opts, family) local implementation = M.diagram.get_type(diagram_type) opts = opts or {} local renderer if operation == "process" then renderer = implementation.process_plot else renderer = implementation[operation] end if type(renderer) ~= "function" then error("luacoolprop: diagram type " .. diagram_type .. " does not implement operation '" .. tostring(operation) .. "'", 2) end if operation == "process" then pending_axis_records = nil print_process(renderer, opts) return end local output if operation == "axis_style" then local records pending_axis_records = nil output, records = renderer(opts) if records and (diagram_type == "TS" or diagram_type == "HS" or diagram_type == "PT") then local snapshot, supported = snapshot_options(opts, {}) if supported then pending_axis_records = { diagram_type = diagram_type, options = snapshot, reference_state = M.reference_state(opts.fluid or "R134a"), records = records, } end end output = "\\pgfplotsset{every axis/.append style={" .. output .. "}}" elseif operation == "plots" then local pending = pending_axis_records pending_axis_records = nil if pending and pending.diagram_type == diagram_type and pending.reference_state == M.reference_state(opts.fluid or "R134a") and same_options(pending.options, opts) then output = renderer(opts, pending.records) else output = renderer(opts) end elseif operation == "family" then pending_axis_records = nil output = renderer(opts, family) else pending_axis_records = nil output = renderer(opts) end tex.sprint(output) end return M