|
| 1 | +/// Permanent tide systems |
| 2 | +use crate::authoring::*; |
| 3 | + |
| 4 | +// ----- F O R W A R D ----------------------------------------------------------------- |
| 5 | + |
| 6 | +fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { |
| 7 | + let mut successes = 0_usize; |
| 8 | + let n = operands.len(); |
| 9 | + let ellps = op.params.ellps(0); |
| 10 | + let Ok(coefficient) = op.params.real("coefficient") else { |
| 11 | + return successes; |
| 12 | + }; |
| 13 | + |
| 14 | + for i in 0..n { |
| 15 | + let mut coord = operands.get_coord(i); |
| 16 | + let phibar = ellps.latitude_geographic_to_geocentric(coord[1]); |
| 17 | + let s = phibar.sin(); |
| 18 | + coord[2] += coefficient * (-0.198) * (1.5 * s * s - 0.5); |
| 19 | + operands.set_coord(i, &coord); |
| 20 | + successes += 1; |
| 21 | + } |
| 22 | + |
| 23 | + successes |
| 24 | +} |
| 25 | + |
| 26 | +// ----- I N V E R S E ----------------------------------------------------------------- |
| 27 | + |
| 28 | +fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { |
| 29 | + let mut successes = 0_usize; |
| 30 | + let n = operands.len(); |
| 31 | + let ellps = op.params.ellps(0); |
| 32 | + let Ok(coefficient) = op.params.real("coefficient") else { |
| 33 | + return successes; |
| 34 | + }; |
| 35 | + |
| 36 | + for i in 0..n { |
| 37 | + let mut coord = operands.get_coord(i); |
| 38 | + let phibar = ellps.latitude_geographic_to_geocentric(coord[1]); |
| 39 | + let s = phibar.sin(); |
| 40 | + coord[2] -= coefficient * (-0.198) * (1.5 * s * s - 0.5); |
| 41 | + operands.set_coord(i, &coord); |
| 42 | + successes += 1; |
| 43 | + } |
| 44 | + |
| 45 | + successes |
| 46 | +} |
| 47 | + |
| 48 | +// ----- C O N S T R U C T O R --------------------------------------------------------- |
| 49 | + |
| 50 | +// Example... |
| 51 | +#[rustfmt::skip] |
| 52 | +pub const GAMUT: [OpParameter; 5] = [ |
| 53 | + OpParameter::Flag { key: "inv" }, |
| 54 | + OpParameter::Real { key: "k", default: Some(0.3) }, |
| 55 | + OpParameter::Text { key: "ellps", default: Some("GRS80") }, |
| 56 | + OpParameter::Text { key: "from", default: None }, |
| 57 | + OpParameter::Text { key: "to", default: None } |
| 58 | +]; |
| 59 | + |
| 60 | +pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result<Op, Error> { |
| 61 | + let mut op = Op::plain(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT, ctx)?; |
| 62 | + let k = op.params.real("k")?; |
| 63 | + |
| 64 | + let Ok(to) = op.params.text("to") else { |
| 65 | + return Err(Error::MissingParam( |
| 66 | + "permtide: must specify 'to=' as exactly one of {'mean', 'zero', 'free'}".to_string(), |
| 67 | + )); |
| 68 | + }; |
| 69 | + let Ok(from) = op.params.text("from") else { |
| 70 | + return Err(Error::MissingParam( |
| 71 | + "permtide: must specify 'from=' as exactly one of {'mean', 'zero', 'free'}".to_string(), |
| 72 | + )); |
| 73 | + }; |
| 74 | + |
| 75 | + let coefficient = match (to.as_str(), from.as_str()) { |
| 76 | + ("mean", "mean") => 0.0, |
| 77 | + ("mean", "zero") => 1.0, |
| 78 | + ("mean", "free") => 1.0 + k, |
| 79 | + ("zero", "zero") => 0.0, |
| 80 | + ("zero", "mean") => -1.0, |
| 81 | + ("zero", "free") => k, |
| 82 | + ("free", "free") => 0.0, |
| 83 | + ("free", "mean") => -(1.0 + k), |
| 84 | + ("free", "zero") => -k, |
| 85 | + _ => f64::NAN, |
| 86 | + }; |
| 87 | + |
| 88 | + if coefficient.is_nan() { |
| 89 | + return Err(Error::BadParam( |
| 90 | + "'to=' or 'from='".to_string(), |
| 91 | + "must be one of {'mean', 'zero', 'free'}".to_string(), |
| 92 | + )); |
| 93 | + } |
| 94 | + |
| 95 | + op.params.real.insert("coefficient", coefficient); |
| 96 | + Ok(op) |
| 97 | +} |
| 98 | + |
| 99 | +// ----- T E S T S --------------------------------------------------------------------- |
| 100 | + |
| 101 | +#[cfg(test)] |
| 102 | +mod tests { |
| 103 | + use super::*; |
| 104 | + use float_eq::assert_float_eq; |
| 105 | + |
| 106 | + #[test] |
| 107 | + fn permtide() -> Result<(), Error> { |
| 108 | + let mut ctx = Minimal::default(); |
| 109 | + |
| 110 | + // A test point near Copenhagen |
| 111 | + let pnt = Coor4D::geo(55., 12., 0., 0.); |
| 112 | + |
| 113 | + // Mean -> zero |
| 114 | + let op = ctx.op("permtide from=mean to=zero ellps=GRS80")?; |
| 115 | + let mut operands = [pnt]; |
| 116 | + ctx.apply(op, Fwd, &mut operands)?; |
| 117 | + assert_float_eq!(operands[0][2], 0.099407199, abs_all <= 1e-8); |
| 118 | + ctx.apply(op, Inv, &mut operands)?; |
| 119 | + assert_float_eq!(operands[0][2], pnt[2], abs_all <= 1e-12); |
| 120 | + |
| 121 | + // Mean -> free |
| 122 | + let op = ctx.op("permtide from=mean to=free ellps=GRS80")?; |
| 123 | + let mut operands = [pnt]; |
| 124 | + ctx.apply(op, Fwd, &mut operands)?; |
| 125 | + assert_float_eq!(operands[0][2], 0.1292293587824579, abs_all <= 1e-8); |
| 126 | + ctx.apply(op, Inv, &mut operands)?; |
| 127 | + assert_float_eq!(operands[0][2], pnt[2], abs_all <= 1e-12); |
| 128 | + |
| 129 | + // Inversion |
| 130 | + let fwd_op = ctx.op("permtide from=mean to=zero ellps=GRS80")?; |
| 131 | + let inv_op = ctx.op("permtide from=zero to=mean ellps=GRS80 inv")?; |
| 132 | + let mut operands = [pnt]; |
| 133 | + ctx.apply(fwd_op, Fwd, &mut operands)?; |
| 134 | + let fwd_h = operands[0][2]; |
| 135 | + |
| 136 | + let mut operands = [pnt]; |
| 137 | + ctx.apply(inv_op, Fwd, &mut operands)?; |
| 138 | + let inv_h = operands[0][2]; |
| 139 | + assert_float_eq!(fwd_h, inv_h, abs_all <= 1e-20); |
| 140 | + |
| 141 | + // Bad tide system names |
| 142 | + assert!(matches!( |
| 143 | + ctx.op("permtide from=cheese to=zero ellps=GRS80"), |
| 144 | + Err(Error::BadParam(_, _)) |
| 145 | + )); |
| 146 | + |
| 147 | + // Missing tide system names |
| 148 | + assert!(matches!(ctx.op("permtide"), Err(Error::MissingParam(_)))); |
| 149 | + assert!(matches!( |
| 150 | + ctx.op("permtide to=zero ellps=GRS80"), |
| 151 | + Err(Error::MissingParam(_)) |
| 152 | + )); |
| 153 | + assert!(matches!( |
| 154 | + ctx.op("permtide from=mean ellps=GRS80"), |
| 155 | + Err(Error::MissingParam(_)) |
| 156 | + )); |
| 157 | + Ok(()) |
| 158 | + } |
| 159 | +} |
0 commit comments