diff --git a/.agents/skills/ptf-extract/SKILL.md b/.agents/skills/ptf-extract/SKILL.md index f7d3cee..0499c79 100644 --- a/.agents/skills/ptf-extract/SKILL.md +++ b/.agents/skills/ptf-extract/SKILL.md @@ -26,8 +26,11 @@ input error. blockers and set affected functions to `blocked`; otherwise set reviewed, complete functions to `ready-for-implementation`. 4. Run `cargo run --manifest-path codegen/Cargo.toml -- validate` and fix - validation errors before finishing. Validation never justifies inferred - science. + validation errors before finishing. When a nontrivial formula expression is + repeated within one implementation, declare it once as an earlier local + implementation variable and reference that variable; retain the published + numeric lexemes and do not assign extra scientific semantics. Validation + never justifies inferred science. ## Output diff --git a/.agents/skills/ptf-generate/SKILL.md b/.agents/skills/ptf-generate/SKILL.md index 2aefca3..2423cc7 100644 --- a/.agents/skills/ptf-generate/SKILL.md +++ b/.agents/skills/ptf-generate/SKILL.md @@ -25,8 +25,12 @@ status when input validation fails; report the blocking input error. generated structures and classes; `$defs` keys only resolve local references. Do not infer missing science. 3. Validate, generate both retained targets, and run the required verification - gates. A generator capability gap is a blocker, never an invitation to - hand-write that computational target. + gates. Before validation, extract each repeated nontrivial formula + expression within a function into one earlier local implementation variable + and reference it thereafter; retain published numeric lexemes and do not + invent scientific semantics for the calculation intermediate. A generator + capability gap is a blocker, never an invitation to hand-write that + computational target. 4. After every required check passes, change the selected source functions from `ready-for-implementation` to `implemented`, then revalidate, regenerate, and prove the second generation pass is idempotent. diff --git a/codegen/src/formula.rs b/codegen/src/formula.rs index 8d9bbf6..9eb3cd1 100644 --- a/codegen/src/formula.rs +++ b/codegen/src/formula.rs @@ -45,13 +45,13 @@ pub(crate) enum ExprKind { Grouped(Box), } -#[derive(Clone, Copy, Debug, Eq, PartialEq)] +#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)] pub(crate) enum UnaryOp { Plus, Minus, } -#[derive(Clone, Copy, Debug, Eq, PartialEq)] +#[derive(Clone, Copy, Debug, Eq, Hash, PartialEq)] pub(crate) enum BinaryOp { Add, Subtract, diff --git a/codegen/src/model.rs b/codegen/src/model.rs index b3f7628..2939d8e 100644 --- a/codegen/src/model.rs +++ b/codegen/src/model.rs @@ -355,9 +355,16 @@ pub(crate) struct RawVariable { #[derive(Clone, Debug)] pub(crate) struct RawExpression { pub(crate) implementation_path: String, + pub(crate) source_location: SourceLocation, pub(crate) expression: crate::formula::Expr, } +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub(crate) struct SourceLocation { + pub(crate) line: usize, + pub(crate) column: usize, +} + #[cfg(test)] mod tests { use super::*; diff --git a/codegen/src/semantic.rs b/codegen/src/semantic.rs index 609ab07..b322ce7 100644 --- a/codegen/src/semantic.rs +++ b/codegen/src/semantic.rs @@ -1,8 +1,11 @@ -use std::{collections::BTreeMap, fmt}; +use std::{ + collections::{BTreeMap, HashMap}, + fmt, +}; use crate::{ formula::{self, Span}, - model::{RawExpression, RawFunction}, + model::{RawExpression, RawFunction, SourceLocation}, }; #[derive(Clone, Debug, PartialEq)] @@ -108,11 +111,23 @@ pub(crate) struct Error { pub(crate) function: String, pub(crate) implementation_path: String, pub(crate) span: Span, + source_location: Option>, message: String, } impl fmt::Display for Error { fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result { + if let Some(source_location) = self.source_location.as_deref() { + return write!( + formatter, + "{}:{}:{}..{}: {}", + self.specification_path, + source_location.line, + source_location.column + self.span.start, + source_location.column + self.span.end, + self.message, + ); + } write!( formatter, "{} -> function {} -> {}:{}..{}: {}", @@ -141,6 +156,8 @@ pub(crate) fn compile(raw: &RawFunction) -> Result { debug_assert_eq!(inputs[&input.name], index); } + validate_repeated_expressions(raw)?; + let variable_names: BTreeMap<_, _> = raw .variables .iter() @@ -179,6 +196,140 @@ pub(crate) fn compile(raw: &RawFunction) -> Result { }) } +#[derive(Clone, Debug, Eq, Hash, PartialEq)] +enum StructuralExpr { + Number(String), + Variable(String), + Unary { + op: formula::UnaryOp, + operand: Box, + }, + Binary { + op: formula::BinaryOp, + left: Box, + right: Box, + }, + Call { + name: String, + args: Vec, + }, + Grouped(Box), +} + +#[derive(Clone, Debug)] +struct Occurrence { + source_location: SourceLocation, + span: Span, +} + +fn validate_repeated_expressions(raw: &RawFunction) -> Result<(), Error> { + let mut occurrences = HashMap::new(); + for variable in &raw.variables { + let expression = &variable.expression; + if let Some((first, later)) = find_repeated_expression( + &expression.expression, + expression.source_location, + &mut occurrences, + ) { + return Err(source_error( + raw, + first.source_location, + first.span, + format!( + "expression is repeated at {}:{}:{}..{}; it must be extracted into an earlier implementation variable", + raw.specification_path.display(), + later.source_location.line, + later.source_location.column + later.span.start, + later.source_location.column + later.span.end, + ), + )); + } + } + Ok(()) +} + +fn find_repeated_expression( + expression: &formula::Expr, + source_location: SourceLocation, + occurrences: &mut HashMap, +) -> Option<(Occurrence, Occurrence)> { + if requires_extraction(expression) { + let structural = structural_expr(expression); + let occurrence = Occurrence { + source_location, + span: expression.span, + }; + if let Some(first) = occurrences.get(&structural) { + return Some((first.clone(), occurrence)); + } + occurrences.insert(structural, occurrence); + } + + match &expression.kind { + formula::ExprKind::Unary { operand, .. } | formula::ExprKind::Grouped(operand) => { + find_repeated_expression(operand, source_location, occurrences) + } + formula::ExprKind::Binary { left, right, .. } => { + find_repeated_expression(left, source_location, occurrences) + .or_else(|| find_repeated_expression(right, source_location, occurrences)) + } + formula::ExprKind::Call { args, .. } => args + .iter() + .find_map(|argument| find_repeated_expression(argument, source_location, occurrences)), + formula::ExprKind::Number(_) | formula::ExprKind::Variable(_) => None, + } +} + +fn requires_extraction(expression: &formula::Expr) -> bool { + match &expression.kind { + formula::ExprKind::Call { .. } => true, + formula::ExprKind::Binary { + op: formula::BinaryOp::Divide | formula::BinaryOp::Power, + .. + } => true, + formula::ExprKind::Binary { .. } => arithmetic_operation_count(expression) > 1, + formula::ExprKind::Unary { operand, .. } | formula::ExprKind::Grouped(operand) => { + arithmetic_operation_count(expression) > 1 || requires_extraction(operand) + } + formula::ExprKind::Number(_) | formula::ExprKind::Variable(_) => false, + } +} + +fn arithmetic_operation_count(expression: &formula::Expr) -> usize { + match &expression.kind { + formula::ExprKind::Number(_) | formula::ExprKind::Variable(_) => 0, + formula::ExprKind::Call { args, .. } => args.iter().map(arithmetic_operation_count).sum(), + formula::ExprKind::Grouped(operand) => arithmetic_operation_count(operand), + formula::ExprKind::Unary { operand, .. } => 1 + arithmetic_operation_count(operand), + formula::ExprKind::Binary { left, right, .. } => { + 1 + arithmetic_operation_count(left) + arithmetic_operation_count(right) + } + } +} + +fn structural_expr(expression: &formula::Expr) -> StructuralExpr { + match &expression.kind { + formula::ExprKind::Number(number) => StructuralExpr::Number(number.lexeme.clone()), + formula::ExprKind::Variable(name) => StructuralExpr::Variable(name.clone()), + formula::ExprKind::Unary { op, operand } => StructuralExpr::Unary { + op: *op, + operand: Box::new(structural_expr(operand)), + }, + formula::ExprKind::Binary { op, left, right } => StructuralExpr::Binary { + op: *op, + left: Box::new(structural_expr(left)), + right: Box::new(structural_expr(right)), + }, + formula::ExprKind::Call { name, args } => StructuralExpr::Call { + name: name.clone(), + args: args.iter().map(structural_expr).collect(), + }, + formula::ExprKind::Grouped(operand) => { + StructuralExpr::Grouped(Box::new(structural_expr(operand))) + } + } +} + fn insert_name( raw: &RawFunction, names: &mut BTreeMap, @@ -334,6 +485,23 @@ fn error( function: raw.name.clone(), implementation_path: implementation_path.to_owned(), span, + source_location: None, + message: message.into(), + } +} + +fn source_error( + raw: &RawFunction, + source_location: SourceLocation, + span: Span, + message: impl Into, +) -> Error { + Error { + specification_path: raw.specification_path.display().to_string(), + function: raw.name.clone(), + implementation_path: String::new(), + span, + source_location: Some(Box::new(source_location)), message: message.into(), } } @@ -353,7 +521,7 @@ mod tests { use crate::{ formula::parse, - model::{RawExpression, RawFunction, RawInput, RawVariable}, + model::{RawExpression, RawFunction, RawInput, RawVariable, SourceLocation}, }; use super::{BinaryOp, Expr, compile}; @@ -361,6 +529,7 @@ mod tests { fn expression(path: &str, source: &str) -> RawExpression { RawExpression { implementation_path: path.into(), + source_location: SourceLocation { line: 1, column: 1 }, expression: parse(path, source).unwrap(), } } @@ -420,6 +589,147 @@ mod tests { )); } + #[test] + fn permits_repeated_simple_multiplication() { + let raw = function( + &["x", "y"], + vec![ + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[0].expr", "x * y"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[1].expr", "x * y"), + }, + ], + ); + + assert!(compile(&raw).is_ok()); + } + + #[test] + fn rejects_repeated_power_with_both_expression_locations() { + let mut raw = function( + &["x"], + vec![ + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[0].expr", "x ^ 2"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[1].expr", "x ^ 2"), + }, + ], + ); + raw.variables[1].expression.source_location = SourceLocation { line: 2, column: 3 }; + + let error = compile(&raw).unwrap_err().to_string(); + assert!( + error.contains("specs/functions/example.md:1:1..6"), + "{error}" + ); + assert!( + error.contains("specs/functions/example.md:2:3..8"), + "{error}" + ); + assert!( + error.contains("must be extracted into an earlier implementation variable"), + "{error}" + ); + } + + #[test] + fn rejects_repeated_math_function_calls() { + let raw = function( + &["x"], + vec![ + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[0].expr", "sqrt(x)"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[1].expr", "sqrt(x)"), + }, + ], + ); + + assert!( + compile(&raw) + .unwrap_err() + .to_string() + .contains("must be extracted into an earlier implementation variable") + ); + } + + #[test] + fn rejects_repeated_compound_arithmetic_expressions() { + let raw = function( + &["x", "y", "z"], + vec![ + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[0].expr", "x + y * z"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[1].expr", "x + y * z"), + }, + ], + ); + + assert!( + compile(&raw) + .unwrap_err() + .to_string() + .contains("must be extracted into an earlier implementation variable") + ); + } + + #[test] + fn distinguishes_operand_order() { + let raw = function( + &["x", "y"], + vec![ + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[0].expr", "x / y"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[1].expr", "y / x"), + }, + ], + ); + + assert!(compile(&raw).is_ok()); + } + + #[test] + fn permits_references_to_an_earlier_extracted_expression() { + let raw = function( + &["x"], + vec![ + RawVariable { + name: "x_squared".into(), + expression: expression("implementation.variables[0].expr", "x ^ 2"), + }, + RawVariable { + name: "first".into(), + expression: expression("implementation.variables[1].expr", "x_squared + 1"), + }, + RawVariable { + name: "second".into(), + expression: expression("implementation.variables[2].expr", "x_squared + 2"), + }, + ], + ); + + assert!(compile(&raw).is_ok()); + } + #[test] fn rejects_unknown_forward_and_self_references_with_expression_paths() { for (source, expected) in [ diff --git a/codegen/src/specs.rs b/codegen/src/specs.rs index 56c629e..cddf040 100644 --- a/codegen/src/specs.rs +++ b/codegen/src/specs.rs @@ -86,11 +86,21 @@ pub(crate) fn load(root: &Path) -> Result> { )); } } + let expression_locations = match expression_locations(&text, &spec) { + Ok(locations) => locations, + Err(error) => { + errors.push(format!("{}:\n $:\n {error}", path.display())); + continue; + } + }; + let mut expression_locations = expression_locations.into_iter(); let implementations = spec .functions .iter() .map(|function| match &function.implementation { - Some(implementation) => compile(&path, function, implementation).map(Some), + Some(implementation) => { + compile(&path, function, implementation, &mut expression_locations).map(Some) + } None => Ok(None), }) .collect::, _>>(); @@ -138,6 +148,7 @@ fn compile( path: &Path, function: &crate::model::Function, implementation: &Implementation, + expression_locations: &mut impl Iterator, ) -> Result { let raw = RawFunction { specification_path: path.to_owned(), @@ -154,11 +165,19 @@ fn compile( .iter() .enumerate() .map(|(index, variable)| { + let source_location = expression_locations.next().ok_or_else(|| { + format!( + "{} -> function {} -> implementation.variables[{index}].expr: source location is unavailable", + path.display(), + function.name + ) + })?; expression( path, &function.name, format!("implementation.variables[{index}].expr"), &variable.expr, + source_location, ) .map(|expression| RawVariable { name: variable.name.clone(), @@ -176,6 +195,7 @@ fn expression( function: &str, implementation_path: String, source: &str, + source_location: crate::model::SourceLocation, ) -> Result { let location = format!( "{} -> function {function} -> {implementation_path}", @@ -184,11 +204,66 @@ fn expression( formula::parse(location, source) .map(|expression| RawExpression { implementation_path, + source_location, expression, }) .map_err(|error| error.to_string()) } +fn expression_locations( + text: &str, + spec: &Spec, +) -> Result, String> { + let expressions = spec + .functions + .iter() + .filter_map(|function| function.implementation.as_ref()) + .flat_map(|implementation| implementation.variables.iter()) + .map(|variable| variable.expr.as_str()); + let mut cursor = 0; + let mut locations = Vec::new(); + + for expression in expressions { + let Some((offset, next_cursor)) = find_expression(text, cursor, expression) else { + return Err(format!( + "could not locate formula expression `{expression}` in YAML source" + )); + }; + locations.push(location(text, offset)); + cursor = next_cursor; + } + Ok(locations) +} + +fn find_expression(text: &str, mut cursor: usize, expression: &str) -> Option<(usize, usize)> { + while let Some(relative) = text[cursor..].find("expr:") { + let field = cursor + relative; + let value_start = field + "expr:".len(); + let value_end = text[value_start..] + .find("expr:") + .map_or(text.len(), |next| value_start + next); + if let Some(relative) = text[value_start..value_end].find(expression) { + let value = value_start + relative; + return Some((value, value + expression.len())); + } + cursor = field + "expr:".len(); + } + None +} + +fn location(text: &str, offset: usize) -> crate::model::SourceLocation { + let prefix = &text[..offset]; + crate::model::SourceLocation { + line: prefix.bytes().filter(|byte| *byte == b'\n').count() + 1, + column: prefix + .rsplit_once('\n') + .map_or(prefix, |(_, line)| line) + .chars() + .count() + + 1, + } +} + fn validate_output( path: &Path, function: &crate::model::Function, @@ -238,7 +313,7 @@ mod tests { path::{Path, PathBuf}, }; - use super::load; + use super::{find_expression, load, location}; use crate::model::PythonGeneration; fn fixture_root(label: &str) -> PathBuf { @@ -455,4 +530,16 @@ mod tests { assert!(error.contains("filename stem must be an APA-style slug")); } + + #[test] + fn locates_block_and_quoted_formula_values_in_yaml_source() { + let source = "variables:\n - expr: >-\n x ^ 2\n - {expr: 'sqrt(x)'}\n"; + let (power, cursor) = find_expression(source, 0, "x ^ 2").unwrap(); + let (call, _) = find_expression(source, cursor, "sqrt(x)").unwrap(); + + assert_eq!(location(source, power).line, 3); + assert_eq!(location(source, power).column, 7); + assert_eq!(location(source, call).line, 4); + assert_eq!(location(source, call).column, 13); + } } diff --git a/specs/functions/hodnett2002.yaml b/specs/functions/hodnett2002.yaml index c91335b..5f2c177 100644 --- a/specs/functions/hodnett2002.yaml +++ b/specs/functions/hodnett2002.yaml @@ -144,18 +144,20 @@ functions: directly represented. implementation: variables: + - name: clay_squared + expr: clay ^ 2 - name: ln_alpha expr: (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity - 11.331 * ph + 0.019 * silt ^ 2) / 100 - name: alpha expr: exp(ln_alpha) - name: ln_n - expr: (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + 0.0070 * clay ^ 2 - 0.014 * sand * silt) / 100 + expr: (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + 0.0070 * clay_squared - 0.014 * sand * silt) / 100 - name: n expr: exp(ln_n) - name: theta_s expr: (81.799 + 0.099 * clay - 31.420 * bulk_density + 0.018 * cation_exchange_capacity + 0.451 * ph - 0.0005 * sand * clay) / 100 - name: theta_r - expr: (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - 0.831 * ph + 0.0018 * clay ^ 2 + 0.0026 * sand * clay) / 100 + expr: (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - 0.831 * ph + 0.0018 * clay_squared + 0.0026 * sand * clay) / 100 scientific_notes: | ## Supported models diff --git a/specs/functions/mayr1999.yaml b/specs/functions/mayr1999.yaml index 7096ca9..f607ae0 100644 --- a/specs/functions/mayr1999.yaml +++ b/specs/functions/mayr1999.yaml @@ -141,16 +141,22 @@ functions: great care; the paper provides that distribution only graphically. implementation: variables: + - name: silt_squared + expr: silt ^ 2 + - name: sand_squared + expr: sand ^ 2 + - name: organic_carbon_squared + expr: organic_carbon ^ 2 - name: log10_a_hc - expr: -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - 0.0021767278 * silt ^ 2 + 1.438224e-5 * silt ^ 3 + 8.040715e-4 * clay ^ 2 + 0.0044067117 * organic_carbon ^ 2 + expr: -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - 0.0021767278 * silt_squared + 1.438224e-5 * silt ^ 3 + 8.040715e-4 * clay ^ 2 + 0.0044067117 * organic_carbon_squared - name: log10_inv_b_hc - expr: -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density - 0.0497915563 * organic_carbon + 3.294687e-4 * sand ^ 2 - 1.689056e-6 * sand ^ 3 + 0.0011225373 * organic_carbon ^ 2 + expr: -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density - 0.0497915563 * organic_carbon + 3.294687e-4 * sand_squared - 1.689056e-6 * sand ^ 3 + 0.0011225373 * organic_carbon_squared - name: a_hc expr: 10 ^ log10_a_hc - name: b_hc expr: 10 ^ (-log10_inv_b_hc) - name: theta_s - expr: 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt + 0.0064338641 * clay - 0.3028160229 * bulk_density + 1.79762e-5 * sand ^ 2 - 3.134631e-5 * silt ^ 2 + expr: 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt + 0.0064338641 * clay - 0.3028160229 * bulk_density + 1.79762e-5 * sand_squared - 3.134631e-5 * silt_squared scientific_notes: | # Mayr and Jarvis (1999) diff --git a/specs/functions/saxton2006.yaml b/specs/functions/saxton2006.yaml index 1402473..6e10cfe 100644 --- a/specs/functions/saxton2006.yaml +++ b/specs/functions/saxton2006.yaml @@ -133,16 +133,18 @@ functions: - These statistical-average estimates should be calibrated to local measurements when available. implementation: variables: + - name: clay_organic_matter_term + expr: 0.027 * clay * organic_matter - name: theta_1500_preliminary expr: -0.024 * sand + 0.487 * clay + 0.006 * organic_matter + 0.005 * sand * organic_matter - 0.013 * clay * organic_matter + 0.068 * sand * clay + 0.031 - name: theta_1500 expr: theta_1500_preliminary + 0.14 * theta_1500_preliminary - 0.02 - name: theta_33_preliminary - expr: -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + 0.006 * sand * organic_matter - 0.027 * clay * organic_matter + 0.452 * sand * clay + 0.299 + expr: -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + 0.006 * sand * organic_matter - clay_organic_matter_term + 0.452 * sand * clay + 0.299 - name: theta_33 expr: theta_33_preliminary + 1.283 * theta_33_preliminary ^ 2 - 0.374 * theta_33_preliminary - 0.015 - name: theta_s_minus_33_preliminary - expr: 0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter - 0.027 * clay * organic_matter - 0.584 * sand * clay + 0.078 + expr: 0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter - clay_organic_matter_term - 0.584 * sand * clay + 0.078 - name: theta_s_minus_33 expr: theta_s_minus_33_preliminary + 0.636 * theta_s_minus_33_preliminary - 0.107 - name: theta_s @@ -153,10 +155,14 @@ functions: expr: -21.67 * sand - 27.93 * clay - 81.97 * theta_s_minus_33_preliminary + 71.12 * sand * theta_s_minus_33_preliminary + 8.29 * clay * theta_s_minus_33_preliminary + 14.05 * sand * clay + 27.16 - name: air_entry_tension expr: air_entry_preliminary + 0.02 * air_entry_preliminary ^ 2 - 0.113 * air_entry_preliminary - 0.70 + - name: ln_33 + expr: ln(33) + - name: ln_theta_33 + expr: ln(theta_33) - name: retention_b - expr: (ln(1500) - ln(33)) / (ln(theta_33) - ln(theta_1500)) + expr: (ln(1500) - ln_33) / (ln_theta_33 - ln(theta_1500)) - name: retention_a - expr: exp(ln(33) + retention_b * ln(theta_33)) + expr: exp(ln_33 + retention_b * ln_theta_33) - name: conductivity_lambda expr: 1 / retention_b - name: saturated_conductivity @@ -304,10 +310,14 @@ functions: - Use only for the 1500 to 33 kPa segment defined by the source. implementation: variables: + - name: ln_33 + expr: ln(33) + - name: ln_theta_33 + expr: ln(theta_33) - name: retention_b - expr: (ln(1500) - ln(33)) / (ln(theta_33) - ln(theta_1500)) + expr: (ln(1500) - ln_33) / (ln_theta_33 - ln(theta_1500)) - name: retention_a - expr: exp(ln(33) + retention_b * ln(theta_33)) + expr: exp(ln_33 + retention_b * ln_theta_33) - name: tension expr: retention_a * theta ^ (-retention_b) diff --git a/specs/functions/wang2012.yaml b/specs/functions/wang2012.yaml index a640bf4..2e6e2e6 100644 --- a/specs/functions/wang2012.yaml +++ b/specs/functions/wang2012.yaml @@ -115,14 +115,18 @@ functions: variables: - name: soil_organic_carbon_g_per_kg expr: 10.0 * soil_organic_carbon + - name: bulk_density_squared + expr: bulk_density ^ 2 + - name: log10_sand + expr: log10(sand) - name: log10_k_sat_cm_per_day - expr: 1.173 + 0.038 * silt + 0.690 * log10(sand) + 0.865 / sand - 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude + expr: 1.173 + 0.038 * silt + 0.690 * log10_sand + 0.865 / sand - 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude - name: k_sat_cm_per_day expr: 10.0 ^ log10_k_sat_cm_per_day - name: fc_percent - expr: 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + 19.198 / soil_organic_carbon_g_per_kg - 5.448 * bulk_density ^ 2 + 0.044 * soil_organic_carbon_g_per_kg ^ 2 + 1.975 * bulk_density * soil_organic_carbon_g_per_kg + expr: 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - 13.991 * log10_sand + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + 19.198 / soil_organic_carbon_g_per_kg - 5.448 * bulk_density_squared + 0.044 * soil_organic_carbon_g_per_kg ^ 2 + 1.975 * bulk_density * soil_organic_carbon_g_per_kg - name: sswc_percent - expr: 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand + 3.247 / soil_organic_carbon_g_per_kg - 17.096 * bulk_density ^ 2 + expr: 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand + 3.247 / soil_organic_carbon_g_per_kg - 17.096 * bulk_density_squared - name: theta_s expr: sswc_percent / 100.0 - name: theta_fc diff --git a/targets/ptfkit-native/cpp/hodnett2002.cppm b/targets/ptfkit-native/cpp/hodnett2002.cppm index b28060c..591c5c3 100644 --- a/targets/ptfkit-native/cpp/hodnett2002.cppm +++ b/targets/ptfkit-native/cpp/hodnett2002.cppm @@ -85,20 +85,21 @@ struct Hodnett2002PTFResult { inline Hodnett2002PTFResult calc_ptf_hodnett2002(double sand, double silt, double clay, double organic_carbon, double bulk_density, double cation_exchange_capacity, double ph) { + const double clay_squared = clay * clay; const double ln_alpha = (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity - 11.331 * ph + 0.019 * (silt * silt)) / 100.0; const double alpha = std::exp(ln_alpha); const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + - 0.0070 * (clay * clay) - 0.014 * sand * silt) / + 0.0070 * clay_squared - 0.014 * sand * silt) / 100.0; const double n = std::exp(ln_n); const double theta_s = (81.799 + 0.099 * clay - 31.420 * bulk_density + 0.018 * cation_exchange_capacity + 0.451 * ph - 0.0005 * sand * clay) / 100.0; const double theta_r = (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - 0.831 * ph + - 0.0018 * (clay * clay) + 0.0026 * sand * clay) / + 0.0018 * clay_squared + 0.0026 * sand * clay) / 100.0; return Hodnett2002PTFResult{alpha, n, theta_s, theta_r}; } diff --git a/targets/ptfkit-native/cpp/mayr1999.cppm b/targets/ptfkit-native/cpp/mayr1999.cppm index 21208bf..45793b2 100644 --- a/targets/ptfkit-native/cpp/mayr1999.cppm +++ b/targets/ptfkit-native/cpp/mayr1999.cppm @@ -71,20 +71,22 @@ struct Mayr1999PTFResult { [[nodiscard]] inline Mayr1999PTFResult calc_ptf_mayr1999(double sand, double silt, double clay, double bulk_density, double organic_carbon) { + const double silt_squared = silt * silt; + const double sand_squared = sand * sand; + const double organic_carbon_squared = organic_carbon * organic_carbon; const double log10_a_hc = -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - - 0.0021767278 * (silt * silt) + 1.438224e-5 * (silt * silt * silt) + - 8.040715e-4 * (clay * clay) + - 0.0044067117 * (organic_carbon * organic_carbon); + 0.0021767278 * silt_squared + 1.438224e-5 * (silt * silt * silt) + + 8.040715e-4 * (clay * clay) + 0.0044067117 * organic_carbon_squared; const double log10_inv_b_hc = -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density - 0.0497915563 * organic_carbon + - 3.294687e-4 * (sand * sand) - 1.689056e-6 * (sand * sand * sand) + - 0.0011225373 * (organic_carbon * organic_carbon); + 3.294687e-4 * sand_squared - 1.689056e-6 * (sand * sand * sand) + + 0.0011225373 * organic_carbon_squared; const double a_hc = std::pow(10.0, log10_a_hc); const double b_hc = std::pow(10.0, -log10_inv_b_hc); const double theta_s = 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt + 0.0064338641 * clay - 0.3028160229 * bulk_density + - 1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt); + 1.79762e-5 * sand_squared - 3.134631e-5 * silt_squared; return Mayr1999PTFResult{a_hc, b_hc, theta_s}; } diff --git a/targets/ptfkit-native/cpp/saxton2006.cppm b/targets/ptfkit-native/cpp/saxton2006.cppm index d9d87d3..47ff9ae 100644 --- a/targets/ptfkit-native/cpp/saxton2006.cppm +++ b/targets/ptfkit-native/cpp/saxton2006.cppm @@ -147,19 +147,20 @@ struct Saxton2006SalinityResult { */ [[nodiscard]] inline Saxton2006PTFResult calc_ptf_saxton2006(double sand, double clay, double organic_matter) { + const double clay_organic_matter_term = 0.027 * clay * organic_matter; const double theta_1500_preliminary = -0.024 * sand + 0.487 * clay + 0.006 * organic_matter + 0.005 * sand * organic_matter - 0.013 * clay * organic_matter + 0.068 * sand * clay + 0.031; const double theta_1500 = theta_1500_preliminary + 0.14 * theta_1500_preliminary - 0.02; const double theta_33_preliminary = -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + - 0.006 * sand * organic_matter - - 0.027 * clay * organic_matter + 0.452 * sand * clay + 0.299; + 0.006 * sand * organic_matter - clay_organic_matter_term + + 0.452 * sand * clay + 0.299; const double theta_33 = theta_33_preliminary + 1.283 * (theta_33_preliminary * theta_33_preliminary) - 0.374 * theta_33_preliminary - 0.015; const double theta_s_minus_33_preliminary = 0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter - - 0.027 * clay * organic_matter - 0.584 * sand * clay + 0.078; + clay_organic_matter_term - 0.584 * sand * clay + 0.078; const double theta_s_minus_33 = theta_s_minus_33_preliminary + 0.636 * theta_s_minus_33_preliminary - 0.107; const double theta_s = theta_33 + theta_s_minus_33 - 0.097 * sand + 0.043; @@ -171,9 +172,10 @@ inline Saxton2006PTFResult calc_ptf_saxton2006(double sand, double clay, double const double air_entry_tension = air_entry_preliminary + 0.02 * (air_entry_preliminary * air_entry_preliminary) - 0.113 * air_entry_preliminary - 0.70; - const double retention_b = - (std::log(1500.0) - std::log(33.0)) / (std::log(theta_33) - std::log(theta_1500)); - const double retention_a = std::exp(std::log(33.0) + retention_b * std::log(theta_33)); + const double ln_33 = std::log(33.0); + const double ln_theta_33 = std::log(theta_33); + const double retention_b = (std::log(1500.0) - ln_33) / (ln_theta_33 - std::log(theta_1500)); + const double retention_a = std::exp(ln_33 + retention_b * ln_theta_33); const double conductivity_lambda = 1.0 / retention_b; const double saturated_conductivity = 1930.0 * std::pow(theta_s - theta_33, 3.0 - conductivity_lambda); @@ -232,9 +234,10 @@ inline Saxton2006DensityResult calc_ptf_saxton2006_density(double normal_density */ [[nodiscard]] inline double calc_ptf_saxton2006_tension_dry(double theta, double theta_1500, double theta_33) { - const double retention_b = - (std::log(1500.0) - std::log(33.0)) / (std::log(theta_33) - std::log(theta_1500)); - const double retention_a = std::exp(std::log(33.0) + retention_b * std::log(theta_33)); + const double ln_33 = std::log(33.0); + const double ln_theta_33 = std::log(theta_33); + const double retention_b = (std::log(1500.0) - ln_33) / (ln_theta_33 - std::log(theta_1500)); + const double retention_a = std::exp(ln_33 + retention_b * ln_theta_33); return retention_a * std::pow(theta, -retention_b); } diff --git a/targets/ptfkit-native/cpp/wang2012.cppm b/targets/ptfkit-native/cpp/wang2012.cppm index 1198947..7aa7f91 100644 --- a/targets/ptfkit-native/cpp/wang2012.cppm +++ b/targets/ptfkit-native/cpp/wang2012.cppm @@ -70,20 +70,21 @@ inline Wang2012PTFResult calc_ptf_wang2012(double sand, double silt, double clay double bulk_density, double soil_organic_carbon, double altitude) { const double soil_organic_carbon_g_per_kg = 10.0 * soil_organic_carbon; - const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * std::log10(sand) + - 0.865 / sand - 0.030 * bulk_density * silt - + const double bulk_density_squared = bulk_density * bulk_density; + const double log10_sand = std::log10(sand); + const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * log10_sand + 0.865 / sand - + 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude; const double k_sat_cm_per_day = std::pow(10.0, log10_k_sat_cm_per_day); const double fc_percent = 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * std::log10(clay) - - 13.991 * std::log10(sand) + 42.261 * std::log10(soil_organic_carbon_g_per_kg) - - 11.763 / sand + 19.198 / soil_organic_carbon_g_per_kg - - 5.448 * (bulk_density * bulk_density) + + 13.991 * log10_sand + 42.261 * std::log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + + 19.198 / soil_organic_carbon_g_per_kg - 5.448 * bulk_density_squared + 0.044 * (soil_organic_carbon_g_per_kg * soil_organic_carbon_g_per_kg) + 1.975 * bulk_density * soil_organic_carbon_g_per_kg; const double sswc_percent = 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand + 3.247 / soil_organic_carbon_g_per_kg - - 17.096 * (bulk_density * bulk_density); + 17.096 * bulk_density_squared; const double theta_s = sswc_percent / 100.0; const double theta_fc = fc_percent / 100.0; const double k_sat = k_sat_cm_per_day / 8640000.0; diff --git a/targets/ptfkit-native/include/ptfkit/hodnett2002.h b/targets/ptfkit-native/include/ptfkit/hodnett2002.h index 54a3ceb..29c779d 100644 --- a/targets/ptfkit-native/include/ptfkit/hodnett2002.h +++ b/targets/ptfkit-native/include/ptfkit/hodnett2002.h @@ -82,20 +82,21 @@ typedef struct { static inline hodnett2002_ptf_result calc_ptf_hodnett2002(double sand, double silt, double clay, double organic_carbon, double bulk_density, double cation_exchange_capacity, double ph) { + const double clay_squared = clay * clay; const double ln_alpha = (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity - 11.331 * ph + 0.019 * (silt * silt)) / 100.0; const double alpha = exp(ln_alpha); const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + - 0.0070 * (clay * clay) - 0.014 * sand * silt) / + 0.0070 * clay_squared - 0.014 * sand * silt) / 100.0; const double n = exp(ln_n); const double theta_s = (81.799 + 0.099 * clay - 31.420 * bulk_density + 0.018 * cation_exchange_capacity + 0.451 * ph - 0.0005 * sand * clay) / 100.0; const double theta_r = (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - 0.831 * ph + - 0.0018 * (clay * clay) + 0.0026 * sand * clay) / + 0.0018 * clay_squared + 0.0026 * sand * clay) / 100.0; #ifdef __cplusplus return hodnett2002_ptf_result{alpha, n, theta_s, theta_r}; diff --git a/targets/ptfkit-native/include/ptfkit/mayr1999.h b/targets/ptfkit-native/include/ptfkit/mayr1999.h index 7077341..5407799 100644 --- a/targets/ptfkit-native/include/ptfkit/mayr1999.h +++ b/targets/ptfkit-native/include/ptfkit/mayr1999.h @@ -68,20 +68,22 @@ typedef struct { */ static inline mayr1999_ptf_result calc_ptf_mayr1999(double sand, double silt, double clay, double bulk_density, double organic_carbon) { + const double silt_squared = silt * silt; + const double sand_squared = sand * sand; + const double organic_carbon_squared = organic_carbon * organic_carbon; const double log10_a_hc = -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - - 0.0021767278 * (silt * silt) + 1.438224e-5 * (silt * silt * silt) + - 8.040715e-4 * (clay * clay) + - 0.0044067117 * (organic_carbon * organic_carbon); + 0.0021767278 * silt_squared + 1.438224e-5 * (silt * silt * silt) + + 8.040715e-4 * (clay * clay) + 0.0044067117 * organic_carbon_squared; const double log10_inv_b_hc = -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density - 0.0497915563 * organic_carbon + - 3.294687e-4 * (sand * sand) - 1.689056e-6 * (sand * sand * sand) + - 0.0011225373 * (organic_carbon * organic_carbon); + 3.294687e-4 * sand_squared - 1.689056e-6 * (sand * sand * sand) + + 0.0011225373 * organic_carbon_squared; const double a_hc = pow(10.0, log10_a_hc); const double b_hc = pow(10.0, -log10_inv_b_hc); const double theta_s = 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt + 0.0064338641 * clay - 0.3028160229 * bulk_density + - 1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt); + 1.79762e-5 * sand_squared - 3.134631e-5 * silt_squared; #ifdef __cplusplus return mayr1999_ptf_result{a_hc, b_hc, theta_s}; #else diff --git a/targets/ptfkit-native/include/ptfkit/saxton2006.h b/targets/ptfkit-native/include/ptfkit/saxton2006.h index 44b1d1e..daf43b0 100644 --- a/targets/ptfkit-native/include/ptfkit/saxton2006.h +++ b/targets/ptfkit-native/include/ptfkit/saxton2006.h @@ -145,19 +145,20 @@ typedef struct { */ static inline saxton2006_ptf_result calc_ptf_saxton2006(double sand, double clay, double organic_matter) { + const double clay_organic_matter_term = 0.027 * clay * organic_matter; const double theta_1500_preliminary = -0.024 * sand + 0.487 * clay + 0.006 * organic_matter + 0.005 * sand * organic_matter - 0.013 * clay * organic_matter + 0.068 * sand * clay + 0.031; const double theta_1500 = theta_1500_preliminary + 0.14 * theta_1500_preliminary - 0.02; const double theta_33_preliminary = -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + - 0.006 * sand * organic_matter - - 0.027 * clay * organic_matter + 0.452 * sand * clay + 0.299; + 0.006 * sand * organic_matter - clay_organic_matter_term + + 0.452 * sand * clay + 0.299; const double theta_33 = theta_33_preliminary + 1.283 * (theta_33_preliminary * theta_33_preliminary) - 0.374 * theta_33_preliminary - 0.015; const double theta_s_minus_33_preliminary = 0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter - - 0.027 * clay * organic_matter - 0.584 * sand * clay + 0.078; + clay_organic_matter_term - 0.584 * sand * clay + 0.078; const double theta_s_minus_33 = theta_s_minus_33_preliminary + 0.636 * theta_s_minus_33_preliminary - 0.107; const double theta_s = theta_33 + theta_s_minus_33 - 0.097 * sand + 0.043; @@ -169,8 +170,10 @@ static inline saxton2006_ptf_result calc_ptf_saxton2006(double sand, double clay const double air_entry_tension = air_entry_preliminary + 0.02 * (air_entry_preliminary * air_entry_preliminary) - 0.113 * air_entry_preliminary - 0.70; - const double retention_b = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); - const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); + const double ln_33 = log(33.0); + const double ln_theta_33 = log(theta_33); + const double retention_b = (log(1500.0) - ln_33) / (ln_theta_33 - log(theta_1500)); + const double retention_a = exp(ln_33 + retention_b * ln_theta_33); const double conductivity_lambda = 1.0 / retention_b; const double saturated_conductivity = 1930.0 * pow(theta_s - theta_33, 3.0 - conductivity_lambda); @@ -253,8 +256,10 @@ static inline saxton2006_density_result calc_ptf_saxton2006_density(double norma */ static inline double calc_ptf_saxton2006_tension_dry(double theta, double theta_1500, double theta_33) { - const double retention_b = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); - const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); + const double ln_33 = log(33.0); + const double ln_theta_33 = log(theta_33); + const double retention_b = (log(1500.0) - ln_33) / (ln_theta_33 - log(theta_1500)); + const double retention_a = exp(ln_33 + retention_b * ln_theta_33); return retention_a * pow(theta, -retention_b); } diff --git a/targets/ptfkit-native/include/ptfkit/wang2012.h b/targets/ptfkit-native/include/ptfkit/wang2012.h index cdabfe4..371d56d 100644 --- a/targets/ptfkit-native/include/ptfkit/wang2012.h +++ b/targets/ptfkit-native/include/ptfkit/wang2012.h @@ -67,19 +67,21 @@ static inline wang2012_ptf_result calc_ptf_wang2012(double sand, double silt, do double bulk_density, double soil_organic_carbon, double altitude) { const double soil_organic_carbon_g_per_kg = 10.0 * soil_organic_carbon; - const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * log10(sand) + - 0.865 / sand - 0.030 * bulk_density * silt - + const double bulk_density_squared = bulk_density * bulk_density; + const double log10_sand = log10(sand); + const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * log10_sand + 0.865 / sand - + 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude; const double k_sat_cm_per_day = pow(10.0, log10_k_sat_cm_per_day); const double fc_percent = - 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - - 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + - 19.198 / soil_organic_carbon_g_per_kg - 5.448 * (bulk_density * bulk_density) + + 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - 13.991 * log10_sand + + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + + 19.198 / soil_organic_carbon_g_per_kg - 5.448 * bulk_density_squared + 0.044 * (soil_organic_carbon_g_per_kg * soil_organic_carbon_g_per_kg) + 1.975 * bulk_density * soil_organic_carbon_g_per_kg; const double sswc_percent = 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand + 3.247 / soil_organic_carbon_g_per_kg - - 17.096 * (bulk_density * bulk_density); + 17.096 * bulk_density_squared; const double theta_s = sswc_percent / 100.0; const double theta_fc = fc_percent / 100.0; const double k_sat = k_sat_cm_per_day / 8640000.0; diff --git a/targets/ptfkit-py/src/ptfkit/hodnett2002.c b/targets/ptfkit-py/src/ptfkit/hodnett2002.c index be5f9e8..4ad1c9e 100644 --- a/targets/ptfkit-py/src/ptfkit/hodnett2002.c +++ b/targets/ptfkit-py/src/ptfkit/hodnett2002.c @@ -12,13 +12,14 @@ static void calc_ptf_hodnett2002_loop(char **args, const npy_intp *dimensions, const double bulk_density = *(const double *)args[4]; const double cation_exchange_capacity = *(const double *)args[5]; const double ph = *(const double *)args[6]; + const double clay_squared = clay * clay; const double ln_alpha = (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity - 11.331 * ph + 0.019 * (silt * silt)) / 100.0; const double alpha = exp(ln_alpha); const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + - 0.0070 * (clay * clay) - 0.014 * sand * silt) / + 0.0070 * clay_squared - 0.014 * sand * silt) / 100.0; const double n = exp(ln_n); const double theta_s = @@ -26,7 +27,7 @@ static void calc_ptf_hodnett2002_loop(char **args, const npy_intp *dimensions, 0.451 * ph - 0.0005 * sand * clay) / 100.0; const double theta_r = (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - - 0.831 * ph + 0.0018 * (clay * clay) + 0.0026 * sand * clay) / + 0.831 * ph + 0.0018 * clay_squared + 0.0026 * sand * clay) / 100.0; *(double *)args[7] = alpha; *(double *)args[8] = n; diff --git a/targets/ptfkit-py/src/ptfkit/mayr1999.c b/targets/ptfkit-py/src/ptfkit/mayr1999.c index 21ae5c3..b00530b 100644 --- a/targets/ptfkit-py/src/ptfkit/mayr1999.c +++ b/targets/ptfkit-py/src/ptfkit/mayr1999.c @@ -10,21 +10,24 @@ static void calc_ptf_mayr1999_loop(char **args, const npy_intp *dimensions, cons const double clay = *(const double *)args[2]; const double bulk_density = *(const double *)args[3]; const double organic_carbon = *(const double *)args[4]; + const double silt_squared = silt * silt; + const double sand_squared = sand * sand; + const double organic_carbon_squared = organic_carbon * organic_carbon; const double log10_a_hc = -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - - 0.0021767278 * (silt * silt) + - 1.438224e-5 * (silt * silt * silt) + 8.040715e-4 * (clay * clay) + - 0.0044067117 * (organic_carbon * organic_carbon); + 0.0021767278 * silt_squared + 1.438224e-5 * (silt * silt * silt) + + 8.040715e-4 * (clay * clay) + + 0.0044067117 * organic_carbon_squared; const double log10_inv_b_hc = -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density - 0.0497915563 * organic_carbon + - 3.294687e-4 * (sand * sand) - + 3.294687e-4 * sand_squared - 1.689056e-6 * (sand * sand * sand) + - 0.0011225373 * (organic_carbon * organic_carbon); + 0.0011225373 * organic_carbon_squared; const double a_hc = pow(10.0, log10_a_hc); const double b_hc = pow(10.0, -log10_inv_b_hc); const double theta_s = 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt + 0.0064338641 * clay - 0.3028160229 * bulk_density + - 1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt); + 1.79762e-5 * sand_squared - 3.134631e-5 * silt_squared; *(double *)args[5] = a_hc; *(double *)args[6] = b_hc; *(double *)args[7] = theta_s; diff --git a/targets/ptfkit-py/src/ptfkit/saxton2006.c b/targets/ptfkit-py/src/ptfkit/saxton2006.c index 83fcadc..3d10fa0 100644 --- a/targets/ptfkit-py/src/ptfkit/saxton2006.c +++ b/targets/ptfkit-py/src/ptfkit/saxton2006.c @@ -8,19 +8,20 @@ static void calc_ptf_saxton2006_loop(char **args, const npy_intp *dimensions, co const double sand = *(const double *)args[0]; const double clay = *(const double *)args[1]; const double organic_matter = *(const double *)args[2]; + const double clay_organic_matter_term = 0.027 * clay * organic_matter; const double theta_1500_preliminary = -0.024 * sand + 0.487 * clay + 0.006 * organic_matter + 0.005 * sand * organic_matter - 0.013 * clay * organic_matter + 0.068 * sand * clay + 0.031; const double theta_1500 = theta_1500_preliminary + 0.14 * theta_1500_preliminary - 0.02; - const double theta_33_preliminary = - -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + 0.006 * sand * organic_matter - - 0.027 * clay * organic_matter + 0.452 * sand * clay + 0.299; + const double theta_33_preliminary = -0.251 * sand + 0.195 * clay + 0.011 * organic_matter + + 0.006 * sand * organic_matter - + clay_organic_matter_term + 0.452 * sand * clay + 0.299; const double theta_33 = theta_33_preliminary + 1.283 * (theta_33_preliminary * theta_33_preliminary) - 0.374 * theta_33_preliminary - 0.015; const double theta_s_minus_33_preliminary = 0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter - - 0.027 * clay * organic_matter - 0.584 * sand * clay + 0.078; + clay_organic_matter_term - 0.584 * sand * clay + 0.078; const double theta_s_minus_33 = theta_s_minus_33_preliminary + 0.636 * theta_s_minus_33_preliminary - 0.107; const double theta_s = theta_33 + theta_s_minus_33 - 0.097 * sand + 0.043; @@ -32,8 +33,10 @@ static void calc_ptf_saxton2006_loop(char **args, const npy_intp *dimensions, co const double air_entry_tension = air_entry_preliminary + 0.02 * (air_entry_preliminary * air_entry_preliminary) - 0.113 * air_entry_preliminary - 0.70; - const double retention_b = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); - const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); + const double ln_33 = log(33.0); + const double ln_theta_33 = log(theta_33); + const double retention_b = (log(1500.0) - ln_33) / (ln_theta_33 - log(theta_1500)); + const double retention_a = exp(ln_33 + retention_b * ln_theta_33); const double conductivity_lambda = 1.0 / retention_b; const double saturated_conductivity = 1930.0 * pow(theta_s - theta_33, 3.0 - conductivity_lambda); @@ -89,8 +92,10 @@ static void calc_ptf_saxton2006_tension_dry_loop(char **args, const npy_intp *di const double theta = *(const double *)args[0]; const double theta_1500 = *(const double *)args[1]; const double theta_33 = *(const double *)args[2]; - const double retention_b = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); - const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); + const double ln_33 = log(33.0); + const double ln_theta_33 = log(theta_33); + const double retention_b = (log(1500.0) - ln_33) / (ln_theta_33 - log(theta_1500)); + const double retention_a = exp(ln_33 + retention_b * ln_theta_33); const double tension = retention_a * pow(theta, -retention_b); *(double *)args[3] = tension; for (int arg = 0; arg < 4; arg++) diff --git a/targets/ptfkit-py/src/ptfkit/wang2012.c b/targets/ptfkit-py/src/ptfkit/wang2012.c index 240164c..c49dc78 100644 --- a/targets/ptfkit-py/src/ptfkit/wang2012.c +++ b/targets/ptfkit-py/src/ptfkit/wang2012.c @@ -12,19 +12,21 @@ static void calc_ptf_wang2012_loop(char **args, const npy_intp *dimensions, cons const double soil_organic_carbon = *(const double *)args[4]; const double altitude = *(const double *)args[5]; const double soil_organic_carbon_g_per_kg = 10.0 * soil_organic_carbon; - const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * log10(sand) + + const double bulk_density_squared = bulk_density * bulk_density; + const double log10_sand = log10(sand); + const double log10_k_sat_cm_per_day = 1.173 + 0.038 * silt + 0.690 * log10_sand + 0.865 / sand - 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude; const double k_sat_cm_per_day = pow(10.0, log10_k_sat_cm_per_day); const double fc_percent = 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - - 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + - 19.198 / soil_organic_carbon_g_per_kg - 5.448 * (bulk_density * bulk_density) + + 13.991 * log10_sand + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + + 19.198 / soil_organic_carbon_g_per_kg - 5.448 * bulk_density_squared + 0.044 * (soil_organic_carbon_g_per_kg * soil_organic_carbon_g_per_kg) + 1.975 * bulk_density * soil_organic_carbon_g_per_kg; const double sswc_percent = 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand + 3.247 / soil_organic_carbon_g_per_kg - - 17.096 * (bulk_density * bulk_density); + 17.096 * bulk_density_squared; const double theta_s = sswc_percent / 100.0; const double theta_fc = fc_percent / 100.0; const double k_sat = k_sat_cm_per_day / 8640000.0; diff --git a/targets/ptfkit-rs/src/hodnett2002.rs b/targets/ptfkit-rs/src/hodnett2002.rs index 822d218..418f71f 100644 --- a/targets/ptfkit-rs/src/hodnett2002.rs +++ b/targets/ptfkit-rs/src/hodnett2002.rs @@ -85,6 +85,7 @@ pub fn calc_ptf_hodnett2002( cation_exchange_capacity: f64, ph: f64, ) -> Hodnett2002PTFResult { + let clay_squared = clay.powi(2); let ln_alpha = (-2.294f64 - 3.526f64 * silt + 2.440f64 * organic_carbon - 0.076f64 * cation_exchange_capacity - 11.331f64 * ph @@ -93,7 +94,7 @@ pub fn calc_ptf_hodnett2002( let alpha = ln_alpha.exp(); let ln_n = (62.986f64 - 0.833f64 * clay - 0.529f64 * organic_carbon + 0.593f64 * ph - + 0.0070f64 * clay.powi(2) + + 0.0070f64 * clay_squared - 0.014f64 * sand * silt) / 100.0f64; let n = ln_n.exp(); @@ -104,7 +105,7 @@ pub fn calc_ptf_hodnett2002( / 100.0f64; let theta_r = (22.733f64 - 0.164f64 * sand + 0.235f64 * cation_exchange_capacity - 0.831f64 * ph - + 0.0018f64 * clay.powi(2) + + 0.0018f64 * clay_squared + 0.0026f64 * sand * clay) / 100.0f64; Hodnett2002PTFResult { diff --git a/targets/ptfkit-rs/src/mayr1999.rs b/targets/ptfkit-rs/src/mayr1999.rs index 5053a76..5a324d2 100644 --- a/targets/ptfkit-rs/src/mayr1999.rs +++ b/targets/ptfkit-rs/src/mayr1999.rs @@ -72,28 +72,31 @@ pub fn calc_ptf_mayr1999( bulk_density: f64, organic_carbon: f64, ) -> Mayr1999PTFResult { + let silt_squared = silt.powi(2); + let sand_squared = sand.powi(2); + let organic_carbon_squared = organic_carbon.powi(2); let log10_a_hc = -4.9840297533f64 + 0.0509226283f64 * sand + 0.1575152771f64 * silt + 0.1240901644f64 * bulk_density - 0.1640033143f64 * organic_carbon - - 0.0021767278f64 * silt.powi(2) + - 0.0021767278f64 * silt_squared + 1.438224e-5f64 * silt.powi(3) + 8.040715e-4f64 * clay.powi(2) - + 0.0044067117f64 * organic_carbon.powi(2); + + 0.0044067117f64 * organic_carbon_squared; let log10_inv_b_hc = -0.8466880654f64 - 0.0046806123f64 * sand + 0.0092463819f64 * silt - 0.4542769707f64 * bulk_density - 0.0497915563f64 * organic_carbon - + 3.294687e-4f64 * sand.powi(2) + + 3.294687e-4f64 * sand_squared - 1.689056e-6f64 * sand.powi(3) - + 0.0011225373f64 * organic_carbon.powi(2); + + 0.0011225373f64 * organic_carbon_squared; let a_hc = 10.0f64.powf(log10_a_hc); let b_hc = 10.0f64.powf(-log10_inv_b_hc); let theta_s = 0.2345971971f64 + 0.0046614221f64 * sand + 0.0088163314f64 * silt + 0.0064338641f64 * clay - 0.3028160229f64 * bulk_density - + 1.79762e-5f64 * sand.powi(2) - - 3.134631e-5f64 * silt.powi(2); + + 1.79762e-5f64 * sand_squared + - 3.134631e-5f64 * silt_squared; Mayr1999PTFResult { a_hc, b_hc, diff --git a/targets/ptfkit-rs/src/saxton2006.rs b/targets/ptfkit-rs/src/saxton2006.rs index f1371ba..f464a3a 100644 --- a/targets/ptfkit-rs/src/saxton2006.rs +++ b/targets/ptfkit-rs/src/saxton2006.rs @@ -73,6 +73,7 @@ These statistical-average estimates should be calibrated to local measurements w #[cfg_attr(feature = "inline", inline)] #[must_use] pub fn calc_ptf_saxton2006(sand: f64, clay: f64, organic_matter: f64) -> Saxton2006PTFResult { + let clay_organic_matter_term = 0.027f64 * clay * organic_matter; let theta_1500_preliminary = -0.024f64 * sand + 0.487f64 * clay + 0.006f64 * organic_matter @@ -85,7 +86,7 @@ pub fn calc_ptf_saxton2006(sand: f64, clay: f64, organic_matter: f64) -> Saxton2 + 0.195f64 * clay + 0.011f64 * organic_matter + 0.006f64 * sand * organic_matter - - 0.027f64 * clay * organic_matter + - clay_organic_matter_term + 0.452f64 * sand * clay + 0.299f64; let theta_33 = theta_33_preliminary + 1.283f64 * theta_33_preliminary.powi(2) @@ -94,7 +95,7 @@ pub fn calc_ptf_saxton2006(sand: f64, clay: f64, organic_matter: f64) -> Saxton2 let theta_s_minus_33_preliminary = 0.278f64 * sand + 0.034f64 * clay + 0.022f64 * organic_matter - 0.018f64 * sand * organic_matter - - 0.027f64 * clay * organic_matter + - clay_organic_matter_term - 0.584f64 * sand * clay + 0.078f64; let theta_s_minus_33 = @@ -110,8 +111,10 @@ pub fn calc_ptf_saxton2006(sand: f64, clay: f64, organic_matter: f64) -> Saxton2 let air_entry_tension = air_entry_preliminary + 0.02f64 * air_entry_preliminary.powi(2) - 0.113f64 * air_entry_preliminary - 0.70f64; - let retention_b = (1500.0f64.ln() - 33.0f64.ln()) / (theta_33.ln() - theta_1500.ln()); - let retention_a = (33.0f64.ln() + retention_b * theta_33.ln()).exp(); + let ln_33 = 33.0f64.ln(); + let ln_theta_33 = theta_33.ln(); + let retention_b = (1500.0f64.ln() - ln_33) / (ln_theta_33 - theta_1500.ln()); + let retention_a = (ln_33 + retention_b * ln_theta_33).exp(); let conductivity_lambda = 1.0f64 / retention_b; let saturated_conductivity = 1930.0f64 * (theta_s - theta_33).powf(3.0f64 - conductivity_lambda); @@ -328,8 +331,10 @@ Use only for the 1500 to 33 kPa segment defined by the source."] #[cfg_attr(feature = "inline", inline)] #[must_use] pub fn calc_ptf_saxton2006_tension_dry(theta: f64, theta_1500: f64, theta_33: f64) -> f64 { - let retention_b = (1500.0f64.ln() - 33.0f64.ln()) / (theta_33.ln() - theta_1500.ln()); - let retention_a = (33.0f64.ln() + retention_b * theta_33.ln()).exp(); + let ln_33 = 33.0f64.ln(); + let ln_theta_33 = theta_33.ln(); + let retention_b = (1500.0f64.ln() - ln_33) / (ln_theta_33 - theta_1500.ln()); + let retention_a = (ln_33 + retention_b * ln_theta_33).exp(); retention_a * theta.powf(-retention_b) } #[cfg(test)] diff --git a/targets/ptfkit-rs/src/wang2012.rs b/targets/ptfkit-rs/src/wang2012.rs index b0c348c..cf6eb20 100644 --- a/targets/ptfkit-rs/src/wang2012.rs +++ b/targets/ptfkit-rs/src/wang2012.rs @@ -71,24 +71,26 @@ pub fn calc_ptf_wang2012( altitude: f64, ) -> Wang2012PTFResult { let soil_organic_carbon_g_per_kg = 10.0f64 * soil_organic_carbon; + let bulk_density_squared = bulk_density.powi(2); + let log10_sand = sand.log10(); let log10_k_sat_cm_per_day = - 1.173f64 + 0.038f64 * silt + 0.690f64 * sand.log10() + 0.865f64 / sand + 1.173f64 + 0.038f64 * silt + 0.690f64 * log10_sand + 0.865f64 / sand - 0.030f64 * bulk_density * silt - 0.00000995f64 * soil_organic_carbon_g_per_kg * altitude; let k_sat_cm_per_day = 10.0f64.powf(log10_k_sat_cm_per_day); let fc_percent = 46.481f64 - 4.757f64 * soil_organic_carbon_g_per_kg - 14.028f64 * clay.log10() - - 13.991f64 * sand.log10() + - 13.991f64 * log10_sand + 42.261f64 * soil_organic_carbon_g_per_kg.log10() - 11.763f64 / sand + 19.198f64 / soil_organic_carbon_g_per_kg - - 5.448f64 * bulk_density.powi(2) + - 5.448f64 * bulk_density_squared + 0.044f64 * soil_organic_carbon_g_per_kg.powi(2) + 1.975f64 * bulk_density * soil_organic_carbon_g_per_kg; let sswc_percent = 98.813f64 - 21.555f64 / bulk_density - 39.735f64 / silt - 2.091f64 / sand + 3.247f64 / soil_organic_carbon_g_per_kg - - 17.096f64 * bulk_density.powi(2); + - 17.096f64 * bulk_density_squared; let theta_s = sswc_percent / 100.0f64; let theta_fc = fc_percent / 100.0f64; let k_sat = k_sat_cm_per_day / 8640000.0f64;