311 lines
10 KiB
Rust
311 lines
10 KiB
Rust
use crate::domain::{
|
||
AnalysisResult, CalculationError, CalculationResult, Conclusion, DecorClass, IndexResult,
|
||
MaterialType, NuclideMeasurements, NuclideResult, SampleInput, Validity, Verdict,
|
||
};
|
||
use crate::mcm::run_monte_carlo;
|
||
|
||
/// 相对（扩展）不确定度可接受上限，用于 3.1 有效性判定与 legacy conclusion。
|
||
const ACCEPTANCE_LIMIT_PERCENT: f64 = 20.0;
|
||
/// 3.1 低活度豁免阈值：总比活度 ≤ 37 Bq/kg 时结果直接有效。
|
||
const TOTAL_ACTIVITY_EXEMPT: f64 = 37.0;
|
||
/// 指数的包含因子 k（GUM 2.2.4 / 2.2.6）。
|
||
const COVERAGE_FACTOR: f64 = 2.0;
|
||
|
||
pub fn calculate_sample(input: SampleInput) -> Result<CalculationResult, CalculationError> {
|
||
validate_input(&input)?;
|
||
|
||
let n = input.ra.measured_values.len();
|
||
let ra = calculate_nuclide(&input.ra)?;
|
||
let th = calculate_nuclide(&input.th)?;
|
||
let k = calculate_nuclide(&input.k)?;
|
||
|
||
let ira = calculate_ira(&ra);
|
||
let ir = calculate_ir(&ra, &th, &k);
|
||
let conclusion = if ira.relative_uncertainty_percent <= ACCEPTANCE_LIMIT_PERCENT
|
||
&& ir.relative_uncertainty_percent <= ACCEPTANCE_LIMIT_PERCENT
|
||
{
|
||
Conclusion::Ok
|
||
} else if n < 6 {
|
||
Conclusion::IncreaseMeasurementsToSix
|
||
} else {
|
||
Conclusion::RecalibrateInstrument
|
||
};
|
||
|
||
let analysis = analyze(input.material_type, &ra, &th, &k, &ira, &ir);
|
||
let mcm = run_monte_carlo(&ra, &th, &k, &input.material_type.primary_limits());
|
||
|
||
Ok(CalculationResult {
|
||
measurement_count: n,
|
||
ra,
|
||
th,
|
||
k,
|
||
ira,
|
||
ir,
|
||
conclusion,
|
||
analysis,
|
||
mcm,
|
||
})
|
||
}
|
||
|
||
fn validate_input(input: &SampleInput) -> Result<(), CalculationError> {
|
||
let counts = [
|
||
("Ra", input.ra.measured_values.len()),
|
||
("Th", input.th.measured_values.len()),
|
||
("K", input.k.measured_values.len()),
|
||
];
|
||
|
||
for (nuclide, count) in counts {
|
||
if count < 1 {
|
||
return Err(CalculationError::TooFewMeasurements { nuclide, count });
|
||
}
|
||
}
|
||
|
||
if input.ra.measured_values.len() != input.th.measured_values.len()
|
||
|| input.ra.measured_values.len() != input.k.measured_values.len()
|
||
{
|
||
return Err(CalculationError::MismatchedMeasurementCounts);
|
||
}
|
||
|
||
validate_nuclide("Ra", &input.ra)?;
|
||
validate_nuclide("Th", &input.th)?;
|
||
validate_nuclide("K", &input.k)?;
|
||
|
||
Ok(())
|
||
}
|
||
|
||
fn validate_nuclide(
|
||
nuclide: &'static str,
|
||
measurements: &NuclideMeasurements,
|
||
) -> Result<(), CalculationError> {
|
||
let calibration = measurements.calibration;
|
||
if !calibration.factor.is_finite()
|
||
|| !calibration.expanded_uncertainty_percent.is_finite()
|
||
|| !calibration.coverage_factor.is_finite()
|
||
|| calibration.coverage_factor == 0.0
|
||
{
|
||
return Err(CalculationError::InvalidCalibration { nuclide });
|
||
}
|
||
|
||
for value in &measurements.measured_values {
|
||
if !value.is_finite() {
|
||
return Err(CalculationError::NonFiniteValue {
|
||
field: "measured value",
|
||
});
|
||
}
|
||
}
|
||
|
||
Ok(())
|
||
}
|
||
|
||
fn calculate_nuclide(measurements: &NuclideMeasurements) -> Result<NuclideResult, CalculationError> {
|
||
let factor = measurements.calibration.factor;
|
||
let mean_measured = mean(&measurements.measured_values);
|
||
let mean_calibrated = mean_measured * factor;
|
||
let type_a_uncertainty = type_a_uncertainty(&measurements.measured_values)?;
|
||
let type_b_relative = measurements.calibration.expanded_uncertainty_percent
|
||
/ 100.0
|
||
/ measurements.calibration.coverage_factor;
|
||
let type_b_uncertainty = factor * type_b_relative;
|
||
let sensitivity_coefficient = mean_measured;
|
||
// 校准比活度 C = mean·a，对测量值 A 的灵敏系数为 a，故 A 类项为 a·uA。
|
||
let combined_uncertainty = ((factor * type_a_uncertainty).powi(2)
|
||
+ (sensitivity_coefficient * type_b_uncertainty).powi(2))
|
||
.sqrt();
|
||
|
||
Ok(NuclideResult {
|
||
mean_measured,
|
||
mean_calibrated,
|
||
type_a_uncertainty,
|
||
type_b_relative,
|
||
type_b_uncertainty,
|
||
sensitivity_coefficient,
|
||
combined_uncertainty,
|
||
})
|
||
}
|
||
|
||
fn calculate_ira(ra: &NuclideResult) -> IndexResult {
|
||
let value = ra.mean_calibrated / 200.0;
|
||
let standard_uncertainty = ra.combined_uncertainty / 200.0;
|
||
make_index(value, standard_uncertainty)
|
||
}
|
||
|
||
fn calculate_ir(ra: &NuclideResult, th: &NuclideResult, k: &NuclideResult) -> IndexResult {
|
||
let value = ra.mean_calibrated / 370.0 + th.mean_calibrated / 260.0 + k.mean_calibrated / 4200.0;
|
||
let standard_uncertainty = ((ra.combined_uncertainty / 370.0).powi(2)
|
||
+ (th.combined_uncertainty / 260.0).powi(2)
|
||
+ (k.combined_uncertainty / 4200.0).powi(2))
|
||
.sqrt();
|
||
make_index(value, standard_uncertainty)
|
||
}
|
||
|
||
/// 由指数值与标准不确定度构造完整的 `IndexResult`（含扩展不确定度与 GUM 真值区间）。
|
||
fn make_index(value: f64, standard_uncertainty: f64) -> IndexResult {
|
||
let expanded_uncertainty = standard_uncertainty * COVERAGE_FACTOR;
|
||
IndexResult {
|
||
value,
|
||
standard_uncertainty,
|
||
expanded_uncertainty,
|
||
relative_uncertainty_percent: relative_percent(standard_uncertainty, value),
|
||
relative_expanded_uncertainty_percent: relative_percent(expanded_uncertainty, value),
|
||
p2_5: value - expanded_uncertainty,
|
||
p97_5: value + expanded_uncertainty,
|
||
}
|
||
}
|
||
|
||
/// 3.1 有效性 + 3.2 临界值判定。
|
||
fn analyze(
|
||
material: MaterialType,
|
||
ra: &NuclideResult,
|
||
th: &NuclideResult,
|
||
k: &NuclideResult,
|
||
ira: &IndexResult,
|
||
ir: &IndexResult,
|
||
) -> AnalysisResult {
|
||
let total_calibrated_activity = ra.mean_calibrated + th.mean_calibrated + k.mean_calibrated;
|
||
|
||
// 3.1 有效性：低活度豁免，否则看 IRa 的相对扩展不确定度（k=2）。
|
||
let validity = if total_calibrated_activity <= TOTAL_ACTIVITY_EXEMPT {
|
||
Validity::LowActivityExempt
|
||
} else if ira.relative_expanded_uncertainty_percent <= ACCEPTANCE_LIMIT_PERCENT {
|
||
Validity::UncertaintyAcceptable
|
||
} else {
|
||
Validity::Invalid
|
||
};
|
||
|
||
// 3.2 临界值判定。
|
||
let verdict = if validity == Validity::Invalid {
|
||
Verdict::InvalidResult
|
||
} else {
|
||
match material {
|
||
MaterialType::DecorativeMaterial => judge_decorative(ira, ir),
|
||
single_tier => judge_single_tier(single_tier, ira, ir),
|
||
}
|
||
};
|
||
|
||
AnalysisResult {
|
||
total_calibrated_activity,
|
||
validity,
|
||
verdict,
|
||
}
|
||
}
|
||
|
||
/// 单个指数真值区间相对极限值的三态。
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||
enum TierCheck {
|
||
/// 区间整体在限值之下：合格。
|
||
Pass,
|
||
/// 区间整体在限值之上：超标。
|
||
Fail,
|
||
/// 区间跨越限值：需增加测量次数。
|
||
Straddle,
|
||
}
|
||
|
||
/// 判断指数真值区间相对极限值的位置；`None` 表示该指数无约束。
|
||
fn check_index(index: &IndexResult, limit: Option<f64>) -> TierCheck {
|
||
match limit {
|
||
None => TierCheck::Pass,
|
||
Some(limit) => {
|
||
if index.p97_5 < limit {
|
||
TierCheck::Pass
|
||
} else if index.p2_5 > limit {
|
||
TierCheck::Fail
|
||
} else {
|
||
TierCheck::Straddle
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/// 合并同一级别多个指数的判定：超标优先，其次跨越，全部合格才合格。
|
||
fn combine(checks: &[TierCheck]) -> TierCheck {
|
||
if checks.contains(&TierCheck::Fail) {
|
||
TierCheck::Fail
|
||
} else if checks.contains(&TierCheck::Straddle) {
|
||
TierCheck::Straddle
|
||
} else {
|
||
TierCheck::Pass
|
||
}
|
||
}
|
||
|
||
/// 主体/空心材料：单级判定。
|
||
fn judge_single_tier(material: MaterialType, ira: &IndexResult, ir: &IndexResult) -> Verdict {
|
||
let tier = &material.tiers()[0];
|
||
match combine(&[
|
||
check_index(ira, tier.ira_limit),
|
||
check_index(ir, tier.ir_limit),
|
||
]) {
|
||
TierCheck::Pass => Verdict::Qualified,
|
||
TierCheck::Straddle => Verdict::NeedMoreMeasurements,
|
||
TierCheck::Fail => Verdict::Unqualified,
|
||
}
|
||
}
|
||
|
||
/// 装饰装修材料：A→B→C 级联分类。
|
||
fn judge_decorative(ira: &IndexResult, ir: &IndexResult) -> Verdict {
|
||
let tiers = MaterialType::DecorativeMaterial.tiers();
|
||
let classes = [DecorClass::A, DecorClass::B, DecorClass::C];
|
||
for (tier, class) in tiers.iter().zip(classes) {
|
||
match combine(&[
|
||
check_index(ira, tier.ira_limit),
|
||
check_index(ir, tier.ir_limit),
|
||
]) {
|
||
TierCheck::Pass => return Verdict::DecorativeClass(class),
|
||
TierCheck::Straddle => return Verdict::NeedMoreMeasurements,
|
||
TierCheck::Fail => continue,
|
||
}
|
||
}
|
||
Verdict::DecorativeClass(DecorClass::Unqualified)
|
||
}
|
||
|
||
fn type_a_uncertainty(values: &[f64]) -> Result<f64, CalculationError> {
|
||
let n = values.len();
|
||
if n <= 1 {
|
||
// 单次测量：A 类不确定度为 0（PDF 2.2.1）。
|
||
Ok(0.0)
|
||
} else if n >= 6 {
|
||
Ok(sample_standard_deviation(values) / (n as f64).sqrt())
|
||
} else {
|
||
let coefficient =
|
||
range_coefficient(n).ok_or(CalculationError::UnsupportedRangeMethodCount { count: n })?;
|
||
let (min, max) = min_max(values);
|
||
Ok((max - min) / (coefficient * (n as f64).sqrt()))
|
||
}
|
||
}
|
||
|
||
fn sample_standard_deviation(values: &[f64]) -> f64 {
|
||
let avg = mean(values);
|
||
let sum_squared: f64 = values.iter().map(|value| (value - avg).powi(2)).sum();
|
||
(sum_squared / (values.len() as f64 - 1.0)).sqrt()
|
||
}
|
||
|
||
fn mean(values: &[f64]) -> f64 {
|
||
values.iter().sum::<f64>() / values.len() as f64
|
||
}
|
||
|
||
fn min_max(values: &[f64]) -> (f64, f64) {
|
||
values
|
||
.iter()
|
||
.fold((f64::INFINITY, f64::NEG_INFINITY), |(min, max), value| {
|
||
(min.min(*value), max.max(*value))
|
||
})
|
||
}
|
||
|
||
fn range_coefficient(n: usize) -> Option<f64> {
|
||
match n {
|
||
2 => Some(1.13),
|
||
3 => Some(1.69),
|
||
4 => Some(2.06),
|
||
5 => Some(2.33),
|
||
6 => Some(2.53),
|
||
_ => None,
|
||
}
|
||
}
|
||
|
||
fn relative_percent(uncertainty: f64, value: f64) -> f64 {
|
||
if value == 0.0 {
|
||
f64::INFINITY
|
||
} else {
|
||
uncertainty / value.abs() * 100.0
|
||
}
|
||
}
|