From b3265463a4cdf52c3993d909b8f87a47652afec9 Mon Sep 17 00:00:00 2001
From: Kirill Vasin <krikus.ms@gmail.com>
Date: Tue, 29 Oct 2024 19:48:06 +0100
Subject: [PATCH] manual phase unwrapping

---
 Kernel/Transmission.wl | 41 ++++++++++++++++++++++++++++++++++++++++-
 README.md              | 32 ++++++++++++++++++++++++++++++--
 2 files changed, 70 insertions(+), 3 deletions(-)

diff --git a/Kernel/Transmission.wl b/Kernel/Transmission.wl
index 04d20f3..c14cf2e 100644
--- a/Kernel/Transmission.wl
+++ b/Kernel/Transmission.wl
@@ -254,7 +254,7 @@ validateOptions[OptionsPattern[] ] := With[{},
   ]
 ]
 
-TransmissionUnwrap[t: TransmissionObject[a_], "Basic", OptionsPattern[]] := With[{
+TransmissionUnwrap[t: TransmissionObject[a_], "Basic" | Automatic, OptionsPattern[]] := With[{
   offset = offsetPhase[a],
   th = OptionValue["PhaseThreshold"]//N,
   phaseShift = OptionValue["PhaseShift"]
@@ -275,5 +275,44 @@ TransmissionUnwrap[t: TransmissionObject[a_], "Basic", OptionsPattern[]] := With
 
 Options[TransmissionUnwrap] = {"PhaseThreshold"->5.6, "PhaseShift"->0};
 
+
+TransmissionUnwrap[t: TransmissionObject[a_], "Held" | "Hold", OptionsPattern[]] := With[{
+  th = OptionValue["PhaseThreshold"]//N,
+  phaseShift = OptionValue["PhaseShift"]
+},
+  applyDefferedBranching[t, th, phaseShift]
+]
+
+
+splitPhase[tr_, th_:3.16] := With[{
+  phase = QuantityMagnitude[tr["Phase Features"], {1/"Centimeters", 1}]
+},
+  Split[phase, Not[Abs[(#2[[2]]-#1[[2]])] > th] &]
+]
+
+autoAdjust[joints_] := Map[Function[item, {item[[1]], item[[1]]}], joints];
+recombinePhase[sausages_, joints_] := With[{
+  accumulated = Accumulate[joints[[All,2]]] 
+}, 
+  Flatten[Join[{sausages[[1,All,2]]}, Table[
+    Map[Function[p, p[[2]] + 2Pi accumulated[[i-1]]], sausages[[i]]]
+  , {i, 2, Length[sausages]}]], 1]
+]
+
+applyDefferedBranching[tr_, th_:5.7, shift_:0] := Module[{sausages, joints, offset, freqs},
+  offset = 2 Pi (1/33.356) QuantityMagnitude[tr["\[Delta]t"], "Picoseconds"];
+  freqs = Normal[tr[[1]]["Frequencies"]];
+  sausages = splitPhase[tr, th];
+  joints = Table[{-Sign[sausages[[i+1, 1, 2]] - sausages[[i, -1, 2]]], 0}, {i, Length[sausages]-1}] // autoAdjust;
+
+
+  With[{sausages = sausages, joints = joints, off = offset freqs},
+    transmissionPhaseRecombine[tr, shift, off][{sausages, joints}] // Hold 
+  ]
+]
+
+transmissionPhaseRecombine[tr_, shift_, off_][{sausages_, joints_}] := (Append[tr, {"Phase"->NumericArray[recombinePhase[sausages, joints] + off], "PhaseShift"->shift}]);
+
+
 End[]
 EndPackage[]
\ No newline at end of file
diff --git a/README.md b/README.md
index 5655084..261c845 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@ A small library for high-precision material parameter extraction from time-domai
 
 ## Features
 - Fabry-Pérot deconvolution ⭐️
-- Informed phase unwrapping
+- Informed Automatic / semi-automatic phase unwrapping
 - High precision / various approximation methods for $n$, $\kappa$, and $\alpha$ solving
 - Works for both thin and thick samples
 - **GPU Acceleration** (OpenCL) ⭐️
@@ -197,7 +197,35 @@ TransmissionUnwrap[t_TransmissionObject, type_String, opts___] _TransmissionObje
 
 it performs phase-unwrapping procedure on `t` object and returns a new `TransmissionObject` with a modified phase. There is only one sheme specified by `type` for unwrapping is available for now 
 
-- `"Basic"` : uses informed phase unwrapping based on a time-delay between the sample and the reference signals
+##### type
+###### Automatic unwrapping
+- `"Basic"` or `Automatic` : uses informed phase unwrapping based on a time-delay between the sample and the reference signals
+
+###### Semi-automatic
+- `"Held"` : uses informed phase unwrapping based and returns a held expression with all phase shifts calculated for further modifications by a user.
+
+If `ReleaseHold` is applied the result will be the same as for `Automatic` type. To take advantage of held unwrapping, i.e. apply your adjustmenets to individual points you need to define a wrapper function. For example, this one will act like identity operation
+
+```mathematica
+myPhaseTransform[Hold[callback_[{parts_, joints_}]]] := Module[{
+  myJoints = joints
+},
+  (* modify points *)
+
+  callback[{parts, myJoints}]
+]
+
+newObject = TransmissionUnwrap[object, "Held"] // myPhaseTransform
+```
+
+the structure of `parts` and `joints` are
+```
+{part1  : {{wavenumber1, phase1}, ..}, part2 : {{wavenumberk, phasek}, ..}}
+{joint1 : {sign12, phaseJump12}, ..}
+```
+
+where `sign12` denotes the sign of a phase jump between parts, while `phaseJump12` is suggested by automatica unwrapper an integer to compenstate this discontinuity. 
+
 
 #### Options
 - `"PhaseThreshold"` : sets the treshold for a detector to remove `2Pi` jump. By the default is `5.6`