All Docstrings

A callback for use in simulating either the Pikal or RF model.

Terminates the time integration when the first component of the state vector, u[1], reaches 1e-10 (i.e. approaches zero).

Implemented very simply as endcond(u, t, integ) = u[1] - 1e-10 ContinuousCallback(endcond, terminate!, save_positions=(true, false))

ConstPhysProp(val)

Make a constant physical properties callable.

When called with any number of arguments, these objects simply return the constant value, e.g. cpp = ConstPhysProp(5u"W/m^2/K") can be called as cpp(x) or cpp() or cpp(x, y, z) to get 5u"W/m^2/K".

RampedVariable(constant_setpt)
RampedVariable(setpts, ramprate)
RampedVariable(setpts, ramprates, holds)

A convenience type for computing temperatures, pressures, etc. with multiple setpoints in sequence, and linear interpolation according to a fixed ramp rate between set points

Three main constructors are available: For a non-varying value, call with one argument:

RampedVariable(constant_setpt)

For one ramp from initial value to set point with indefinite hold, call with two arguments:

RampedVariable(setpts, ramprate)

And for multiple setpoints, call with three arguments:

RampedVariable(setpts, ramprates, holds)

With three arguments, setpts, ramprates, and holds should all be vectors, with lengths N+1, N, N-1 respectively.

The resulting RampedVariable rv = RampedVariable(...) can be called as rv(x) at any (dimensionally consistent) value of x, and will return the value at that time point along the ramp process.

A plot recipe is also provided for this type, e.g. plot(rv; tmax=10u"hr") where tmax indicates where to stop drawing the last setpoint hold.

A convenience type for dealing with the common functional form given to Rp and Kv.

An object Rp = RpFormFit(A, B, C) can be called as Rp(x), which simply computes A + B*x/(1 + C*x). Likewise, Kv = RpFormFit(Kc, Kp, Kd) can be called as Kv(p) to get Kc + Kp*p/(1 + Kd*p).

Be careful to pass dimensionally consistent values.

end_cond(u, t, integ)

Compute the end condition for primary drying (that mf or hf approaches zero).

extract_ts(rv; un)
extract_ts(a)

This extracts the time stops, which are typically not smooth, for controlled variables.

These time stops are used to ensure that the ODE solver treats these points as non-smooth, and does not attempt to interpolate across them.

This is defined for RampedVariable and DataInterpolations.AbstractInterpolation objects. For other objects, it falls back to returning [0.0].

get_tstops(controls)

Get the non-smooth time points for a set of time-varying controlled variables.

When called with a tuple, extract_ts is called on each element of the tuple, then those results are merged into a single vector of time points, dimensionless but given in hours.

Each element of the tuple should be a RampedVariable or a DataInterpolations.AbstractInterpolation object.

identify_pd_end(t, pch_pir, kind; window_width, tmin, tmax)

Identify the end of primary drying using second derivative of Pirani pressure.

Arguments

• t: time points
• pch_pir: Pirani pressure at given times
• kind: Val(:der2) for maximum of second derivative, Val(:onoff) for onset and offset.
• window_width: Width of the Savitzky-Golay filter window (default: 91)
• tmin: Minimum time for analysis (default: 0 hours, unitful)
• tmax: Maximum time for analysis (default: Inf hours, unitful)

If a single argument data is passed instead of t and pch_pir, the time and Pirani pressure are assumed to be accessible as either data.t and data.pch_pir or data[:t] and data[:pch_pir].

Approach

Use a Savitzky-Golay filter to smooth and take derivatives. For the onset-offset, find linear intersects of tangent at inflection point with maximum and minimum values of pressure. Those intersects are taken as onset and offset. For the second derivative, identify the maximum of the second derivative.

AbstractExpDatum

An abstract type for all types of experimental data used in fitting.

If you add a new type with this abstract type, you should define the following methods for your type:

• resid_name(::MyDatum)::Symbol: return a unique symbol identifying the type of residuals associated with your data type. This is used to map to a weight in the weights NamedTuple passed to obj_exp and err_exp.
• time_bound_data(::MyDatum): return true if your data type is associated with a time vector, and false otherwise. If true, you should also define:
  • has_timevec(md::MyDatum): return true if the instance md has its own time vector, and false if it uses the container's time vector.
  • nontrivial_t_range(md::MyDatum): return true if the instance md has a nontrivial set of time indices (i.e., not starting at 1 and incrementing by 1) into the container's or instance's time vector, and false otherwise.

For fitting with an optimization solver:

• obj_exp_datum(sol, st, dat::MyDatum; verbose=false) or obj_exp_datum(sol, dat::MyDatum; verbose=false), depending on if time_bound_data is true or false

For fitting with nonlinear least squares solver:

• err_exp_datum!(sol, st, dat::MyDatum, weight; verbose=false) or err_exp_datum!(sol, dat::MyDatum, weight; verbose=false), depending on if time_bound_data is true or false
• num_errs(md::MyDatum): return the number of residuals associated with your data type. This is used to determine the size of the residual vector for nonlinear least squares fitting.

EndTimeData(t_end)

An experimentally measured end of primary drying, for use in ExpFitData.

t_end is typically determined from non-temperature measurements (e.g. Pirani-CM convergence). It can be a single time, or a tuple of two times indicating a window of acceptable drying times: in the objective function, any model drying time inside that window is not penalized, while outside the window a squared error applies, as for the single-time case. If no EndTimeData object is present in the fit, the end of drying is ignored in the objective function.

ExpFitData(t, data)

A type for indicating how experimental data should be fit.

ExpFitData is a container holding a master time vector t and a tuple data of experimental data objects: any number of TfData, TvwSeriesData, TvwEndData, and EndTimeData objects. The order of objects in data determines the order of residuals in the fitting procedure.

Provided constructors:

ExpFitData(t, data)          # data is a tuple of pre-built data objectsExpFitData(t, obj1, obj2, ...)  # convenience: one or more pre-built data objects, in any order

The convenience constructors that accept raw temperature vectors with keyword arguments (Tvws=, t_end=) are only available through the deprecated PrimaryDryFit function.

Each data object is documented with its own fields and constructors. In particular, the temperature series (TfData, TvwSeriesData) carry a t_range selecting the time points they correspond to, and any data object may carry its own time vector t (see fit_t). An end of drying is supplied via an EndTimeData object and is ignored in the objective function if absent.

Currently-implemented data objects include:

• TfData: a single frozen-product temperature series
• TvwSeriesData: a single vial-wall temperature series
• TvwEndData: a single vial-wall temperature endpoint
• EndTimeData: a single end-of-drying time, or a tuple of two times indicating an acceptable window

Common Cases:

• Conventional, single thermocouple: ExpFitData(t, TfData(Tf))
• Conventional, multiple thermocouples: ExpFitData(t, (TfData(Tf1), TfData(Tf2), ...))
• Conventional with Pirani ending: ExpFitData(t, (TfData(Tf), EndTimeData(t_end)))
• RF with measured vial wall: ExpFitData(t, (TfData(Tf), TvwSeriesData(Tvws)))
• RF, matching model Tvw to experimental Tf[end] without measured vial wall: ExpFitData(t, (TfData(Tf), TvwEndData(Tvw_end)))

SolTrim

A struct to contain information about which model points overlap with experiment.

If preinterp is true, then the model solution is pre-interpolated to the experimental time points; the model's time indices 1:(length(sol.t)-1) correspond exactly to the experimental time indices ti_model_start:ti_model_end.

If preinterp is false, then the model solution is not pre-interpolated, and the model solution will need to be interpolated to the experimental time points. at time indices ti_m_start:ti_m_end.

TfData(Tf; t=missing)
TfData(Tf, t_range; t=missing)

A single experimental frozen-product (Tf) temperature series, for use in ExpFitData.

Fields:

• Tf: an AbstractVector of Temperature (one series).
• t_range: a range of integer indices into the time vector that Tf corresponds to; Tf[k] is the temperature measured at time t[t_range[k]].
• t: an optional AbstractVector of Time; when missing (the default), the container's time vector is used (see fit_t).

Constructors:

• TfData(Tf; t=missing): the time-index range defaults to 1:length(Tf), i.e. the series corresponds to the first length(Tf) points of the time vector.
• TfData(Tf, t_range; t=missing): use an explicit t_range (a range of the same length as Tf), e.g. to indicate the series was measured over a sub-window of the time vector.

TvwEndData(Tvw_end)

A single experimental vial-wall (Tvw) temperature endpoint, for use in ExpFitData.

TvwSeriesData(Tvws; t=missing)
TvwSeriesData(Tvws, t_range; t=missing)

A single experimental vial-wall (Tvw) temperature series, for use in ExpFitData.

Fields:

• Tvw: an AbstractVector of Temperature (one series).
• t_range: a range of integer indices into the time vector that Tvw corresponds to; Tvw[k] is the temperature measured at time t[t_range[k]].
• t: an optional AbstractVector of Time; when missing (the default), the container's time vector is used (see fit_t).

Constructors:

• TvwSeriesData(Tvws; t=missing): the time-index range defaults to 1:length(Tvws), i.e. the series corresponds to the first length(Tvws) points of the time vector.
• TvwSeriesData(Tvws, t_range; t=missing): use an explicit t_range (a range of the same length as Tvws), e.g. to indicate the series was measured over a sub-window of the time vector.

PrimaryDryFit(t, Tfs; Tvws, t_end)

Deprecated: this is now a constructor for ExpFitData.

err_exp!(errs, sol, efd; weights, verbose)

Evaluate the error between model solution sol and experimental data in efd.

The in-place version err_exp! fills the errs vector with the errors, while the non-in-place version err_exp returns a new vector of errors. In-place err_exp! thus doesn't allocate, but requires errs to have length num_errs(efd).

In contrast to obj_exp(), which sums all the squared residuals, this function fills the passed array errs with each separate residual, which is suited for least squares algorithms.

• errs is a vector of length num_errs(efd), which this function fills with the errors.

-sol is a solution to an appropriate model; see gen_sol_pd for a helper function.

• efd is an instance of ExpFitData, which contains some information about what to compare.
• weights is a NamedTuple mapping each experimental datum to an inverse-unit weight; the default is constructed with residual_weighting().

Each time series, plus the end time, is given equal weight by dividing by its length; error is given in K (but ustripped).

Note that if efd has vial wall temperatures (i.e. a TvwSeriesData or TvwEndData object in efd.data), the third-index variable in sol is assumed to be temperature, as is true for solutions with ParamObjRF.

If there are multiple series of Tf in efd, error is computed for each separately; likewise for Tvw.

err_exp(sol, efd; weights, verbose)

Evaluate the error between model solution sol and experimental data in efd.

The in-place version err_exp! fills the errs vector with the errors, while the non-in-place version err_exp returns a new vector of errors. In-place err_exp! thus doesn't allocate, but requires errs to have length num_errs(efd).

In contrast to obj_exp(), which sums all the squared residuals, this function fills the passed array errs with each separate residual, which is suited for least squares algorithms.

• errs is a vector of length num_errs(efd), which this function fills with the errors.

-sol is a solution to an appropriate model; see gen_sol_pd for a helper function.

• efd is an instance of ExpFitData, which contains some information about what to compare.
• weights is a NamedTuple mapping each experimental datum to an inverse-unit weight; the default is constructed with residual_weighting().

Each time series, plus the end time, is given equal weight by dividing by its length; error is given in K (but ustripped).

Note that if efd has vial wall temperatures (i.e. a TvwSeriesData or TvwEndData object in efd.data), the third-index variable in sol is assumed to be temperature, as is true for solutions with ParamObjRF.

If there are multiple series of Tf in efd, error is computed for each separately; likewise for Tvw.

exp_time_inds(st)

Return the range of time indices for the experimental time points that match up with the model solution.

fit_t(fitdat, obj)

Return the time vector to use for obj.

Three cases:

• obj.t[obj.t_range] if obj.t is set
• fitdat.t[obj.t_range] if t_range makes sense and obj.t is missing
• fitdat.t if obj is not tied to a time vector

Check if a given instance of an AbstractExpDatum subtype has its own self-contained time vector.

This is strictly false if time_bound_data is false

loss_weighting(; t, Tvw, Tf, kwargs...)

Construct a NamedTuple of weights per experiment type for use with obj_exp.

weights maps the names returned by resid_name to inverse-squared-unit weights, so even custom experimental data types can be weighted against each other without changing LyoPronto. t, Tf, and Tvw have built-in defaults.

The names of each kwarg must match the results of resid_name for each type of AbstractExpDatum, e.g. :Tf, :t, :Tvw for the builtin types.

The value of each kwarg should have inverse-square units, e.g. u"hr^-2" for :t.

model_result(sol, st, idx; unit, verbose)

Compute the model result for a given solution sol with time trimming st, variable idx.

If idx is specified, the function evaluates the solution at that index and gives it the units specified by the unit kwarg (defaulting to u"K").

If idx is not specified, the function evaluates calc_md_Q at all appropriate time points and returns a corresponding Table. If only a single column from that table is desired, pass the keyword argument var to select the column by name (e.g. var=:md).

If you have a choice between the two, it will be more efficient to use idx since it only has to interpolate the solution for that index, while the var option will evaluate the full model at all time points and then select the column.

model_time_inds(st)

Return the range of time indices for the model solution that match up with the experimental time points.

Check if a given instance of an AbstractExpDatum subtype has a nontrivial set of time indices.

"Nontrivial" in this case means that the indices do not start at 1 and/or do not increment by 1, i.e. the data series is not aligned with the start of the time vector or is not contiguous.

num_errs(efd)

Compute the number of data points available in efd for comparison to model solution.

This function is called on each object in efd.data and summed to give the total number of residuals. This is used for caching a residual vector for least-squares fitting, e.g. with err_exp! and nls_pd!.

obj_exp(sol, efd; weights, verbose)

Evaluate an objective function which compares model solution computed by sol to experimental data in efd.

• sol is a solution to an appropriate model; see gen_sol_pd for a helper.
• efd is an instance of ExpFitData, which contains some information about what to compare.
• weights is a NamedTuple maps each experimental datum to an inverse-squared-unit weight; the default is constructed with loss_weighting().

To add a custom type of experimental data, define resid_name(::MyDatum)::Symbol and map that symbol to a weight accessed in a NamedTuple, e.g. weights[:t] should be e.g. 1.0u"hr^-2".

Note that if efd has vial wall temperatures (i.e. a TvwSeriesData or TvwEndData object in efd.data), the third-index variable in sol is assumed to be temperature, as is true for the lumped capacitance model (see ParamObjRF.

If there are multiple series of Tf in efd, squared error is computed for each separately then summed; likewise for Tvw.

obj_expT(sol, efd; tweight, verbose, Tvw_weight)

A thin wrapper on obj_exp, for backwards compatibility.

resid_name(obj::AbstractExpDatum)::Symbol

Return a symbol identifying the type of residuals associated with obj. This is used to map to a weight in the weights NamedTuple passed to obj_exp and err_exp. The default names are:

• :Tf for TfData
• :Tvw for TvwSeriesData and TvwEndData
• :t for EndTimeData

If you define your own type of experimental data, you must define resid_name(::MyDatum)::Symbol to return a unique symbol for your data type, and then map that symbol to a weight, preferably with residual_weighting or loss_weighting, then pass the resulting NamedTuple to obj_exp and err_exp.

residual_weighting(; t, Tvw, Tf, kwargs...)

Construct a NamedTuple of weights per experiment type for use with err_exp.

weights maps the names returned by resid_name to inverse-unit weights, so even custom experimental data types can be weighted against each other without changing LyoPronto. t, Tf, and Tvw have built-in defaults.

The names of each kwarg must match the results of resid_name for each type of AbstractExpDatum, e.g. :Tf, :t, :Tvw for the builtin types.

The value of each kwarg should have inverse units, e.g. u"hr^-1" for :t.

time_bound_data(obj::AbstractExpDatum)::Bool

Trait for identifying whether a given type should be associated with a time vector

trim_sol(sol::ODESolution, t_exp)

Trim the solution sol to the experimental time grid t_exp, returning a SolTrim object.

As a first pass, this should be called on an ExpFitData. Thenshould be called separately for data objects which have their own time vector.

KBB_transform_basic(Kvwfg, Bfg, Bvwg)

Construct a typical transform for fitting Kvwf, Bf, and Bvw (as for a microwave cycle).

KBB_transform_bounded(
    Kvwfg,
    Bfg,
    Bvwg;
    Kvwf_scalefac,
    Bf_scalefac,
    Bvw_scalefac
)

Construct a bounded transform for fitting Kvwf, Bf, and Bvw (as for a microwave cycle).

Kvwf_scalefac, Bf_scalefac, and Bvw_scalefac are used to provide upper and lower bounds on the fitted parameter, as (Kvwfg/Kvwf_scalefac, Kvwfg*Kvwf_scalefac), etc. This is enforced with a logistic transform scaled and shifted appropriately.

KRp_transform_basic(Kshfg, R0g, A1g, A2g)

Construct a typical transform for fitting both Kshf and Rp.

K_transform_basic(Kshfg)

Construct a typical transform for fitting Kshf (a.k.a. Kv).

Rp_transform_basic(R0g, A1g, A2g)

Construct a typical transform for fitting Rp.

gen_nsol_pd(fitlog, tr, pos, fitdats; badprms, kwargs...)

Generate multiple solutions at once.

If the transformation tr makes something (e.g. NamedTuple) with properties separate and shared, then one each of separate is combined with shared, then they are matched up with each element of pos and fitdats. If pos is a single object, it is repeated. Further, if tr also has a field sep_inds, then those indices are used to map separate to the sets of pos and fitdats. This is useful for when e.g. 2 sets of separate parameters are to be applied across 5 different experiments.

gen_sol_pd(fitlog, tr, po; saveat, badprms, kwargs...)

Solve a model for the conditions po with tr(fitlog) mapped on top.

These models will typically be for primary drying, with a flat dimensionless vector fitlog for fitted parameter values, a mapping tr from fitlog to named coefficients, and other parameters in po.

If given, fitdat is used to set saveat for the ODE solution.

The equations used are determined by the type of po, which (with the magic of dispatch) is used to set up an ODE system.

tr should be a TransformTuple object, from TransformVariables, which maps a plain Vector{Float64} to a NamedTuple with fields matching the properties and units of the po object. To construct callable functions, have the transform map to callable structs.

This small function runs

fitprm = transform(tr, fitlog)new_params = setproperties(po, fitprm)!isnothing(badprms) && badprms(new_params) && return Val(NaN)prob = ODEProblem(new_params; tspan=(0.0, 1000.0))sol = solve(prob, Rodas4(autodiff=AutoForwardDiff(chunksize=2)); saveat, kwargs...)

which is wrapped to avoid code duplication.

So, to choose which parameters to fitting, all that is necessary is to provide an appropriate transform tr. Therefore this function can be used for both K-Rp fitting, or just Rp, or just a subset of the 3 Rp coefficients.

Other kwargs are passed directly (as is) to the ODE solve call.

nls_pd!(errs, fitlog, tpf; badprms, weights, verbose)

Calculate the errors for fitting parameters to primary drying data. This directly calls gen_sol_pd, then err_exp!, so see those docstrings.

nls_pd(fitlog, tpf; badprms, weights, verbose)

Calculate the errors for fitting parameters to primary drying data. This directly calls gen_sol_pd, then err_exp, so see those docstrings.

obj_pd(fitlog, tpf; weights, badprms, verbose)

Calculate the sum of squared error (objective function) for fitting parameters to primary drying data. This directly calls gen_sol_pd, then obj_exp, so see those docstrings.

objn_pd(fitlog, tpf; weights, badprms, verbose)

Calculate the sum of squared error (objective function) for fitting parameters to primary drying data. This directly calls gen_nsol_pd, then obj_exp, so see those docstrings.

cycledataplot(tablelike, (Tname1, ...), Tsh_name, (p_name1, ...))
cycledataplot!(tablelike, (Tname1, ...), Tsh_name, (p_name1, ...))

Plot recipe for plotting both temperatures and pressures of a lyophilization cycle.

This requires at least two subplots, one for temperature and one for pressure. Either call twinx(plot()) then cycledataplot!(...), or give a layout like cycledataplot(..., layout=(2,1))"

cycledataplot(tablelike, (Tname1, ...), Tsh_name, (p_name1, ...))
cycledataplot!(tablelike, (Tname1, ...), Tsh_name, (p_name1, ...))

Plot recipe for plotting both temperatures and pressures of a lyophilization cycle.

This requires at least two subplots, one for temperature and one for pressure. Either call twinx(plot()) then cycledataplot!(...), or give a layout like cycledataplot(..., layout=(2,1))"

exppplot(t, p1, [p2, ...,], (name1, [name2, ...]))
exppplot!(t, p1, [p2, ...,], (name1, [name2, ...]))

Plot recipe for plotting pressures of a lyophilization cycle.

exppplot(t, p1, [p2, ...,], (name1, [name2, ...]))
exppplot!(t, p1, [p2, ...,], (name1, [name2, ...]))

Plot recipe for plotting pressures of a lyophilization cycle.

exptfplot(time, T1, [T2, ...]; nmarks=Inf, showline=false, labsuffix=", exp.")
exptfplot!(time, T1, [T2, ...]; nmarks=Inf, showline=false, labsuffix=", exp.")

Plot recipe for one or more experimentally measured product temperatures, all at same times. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{fi}"*labsuffix.

nmarks is an integer indicating how many points to plot with markers, so that the number of markers is reasonable even with many data points. To plot all points (default), pass nmarks=Inf; to plot no points, pass nmarks=0.

This recipe defaults to showing data points without a line, but a line can be added by passing showline=true or specifying a linewidth (or lw).

exptfplot(time, T1, [T2, ...]; nmarks=Inf, showline=false, labsuffix=", exp.")
exptfplot!(time, T1, [T2, ...]; nmarks=Inf, showline=false, labsuffix=", exp.")

Plot recipe for one or more experimentally measured product temperatures, all at same times. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{fi}"*labsuffix.

nmarks is an integer indicating how many points to plot with markers, so that the number of markers is reasonable even with many data points. To plot all points (default), pass nmarks=Inf; to plot no points, pass nmarks=0.

This recipe defaults to showing data points without a line, but a line can be added by passing showline=true or specifying a linewidth (or lw).

exptvwplot(time, T1, [T2, ...]; skip=1, nmarks=Inf, showline=false, labsuffix=", exp.")
exptvwplot!(time, T1, [T2, ...]; skip=1, nmarks=Inf, showline=false, labsuffix=", exp.")

Plot recipe for a set of experimentally measured vial wall temperatures. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{vwi}"*labsuffix. skip is an integer, indicating how many points to skip at a time, so that the dotted line looks dotted even with noisy data.

nmarks is an integer indicating how many points to plot with markers, so that the number of markers is reasonable even with many data points. To plot all points (default), pass nmarks=Inf; to plot no points, pass nmarks=0.

This recipe defaults to showing data points without a line, but a line can be added by passing showline=true or specifying a linewidth (or lw).

exptvwplot(time, T1, [T2, ...]; skip=1, nmarks=Inf, showline=false, labsuffix=", exp.")
exptvwplot!(time, T1, [T2, ...]; skip=1, nmarks=Inf, showline=false, labsuffix=", exp.")

Plot recipe for a set of experimentally measured vial wall temperatures. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{vwi}"*labsuffix. skip is an integer, indicating how many points to skip at a time, so that the dotted line looks dotted even with noisy data.

nmarks is an integer indicating how many points to plot with markers, so that the number of markers is reasonable even with many data points. To plot all points (default), pass nmarks=Inf; to plot no points, pass nmarks=0.

This recipe defaults to showing data points without a line, but a line can be added by passing showline=true or specifying a linewidth (or lw).

modconvtplot(sols; showmarks=false, labsuffix = ", model")
modconvtplot!(sols; showmarks=false, labsuffix = ", model")

Plot recipe for one or multiple solutions to the Pikal model, e.g. the output of gen_sol_pd. This adds a series to the plot for each passed solution, with labels defaulting to "T_{fi}"*labsuffix.

This recipe defaults to showing a line with no markers, but markers can be added by passing showmarks=true or specifying a markershape (or ms).

modconvtplot(sols; showmarks=false, labsuffix = ", model")
modconvtplot!(sols; showmarks=false, labsuffix = ", model")

Plot recipe for one or multiple solutions to the Pikal model, e.g. the output of gen_sol_pd. This adds a series to the plot for each passed solution, with labels defaulting to "T_{fi}"*labsuffix.

This recipe defaults to showing a line with no markers, but markers can be added by passing showmarks=true or specifying a markershape (or ms).

modrftplot(sol, labsuffix=", model", showmarks=false, trimend=0)
modrftplot!(sol, labsuffix=", model", showmarks=false, trimend=0)

Plot recipe for one solution to the lumped capacitance model.

This adds two series to the plot, with labels defaulting to ["T_f" "T_{vw}"] .* labsuffix. The optional argument trimend controls how many time points to trim from the end (which is helpful if temperature shoots up as mf -> 0).

Since this is a recipe, any Plots.jl keyword arguments can be passed to modify the plot.

This recipe defaults to showing a line with no markers, but markers can be added by passing showmarks=true or specifying a markershape (or ms).

modrftplot(sol, labsuffix=", model", showmarks=false, trimend=0)
modrftplot!(sol, labsuffix=", model", showmarks=false, trimend=0)

Plot recipe for one solution to the lumped capacitance model.

This adds two series to the plot, with labels defaulting to ["T_f" "T_{vw}"] .* labsuffix. The optional argument trimend controls how many time points to trim from the end (which is helpful if temperature shoots up as mf -> 0).

Since this is a recipe, any Plots.jl keyword arguments can be passed to modify the plot.

This recipe defaults to showing a line with no markers, but markers can be added by passing showmarks=true or specifying a markershape (or ms).

pressurenames(name)

Attempt to catch common names or shorthands for Pirani and capacitance manometer, then give back a nice label.

summary_md_Q(sol)

This function makes a Table of the results of calc_md_Q at all of the given ODESolution's time points.

For more performance-sensitive uses (e.g. in parameter fitting), see model_result.

tendplot(t_end)
tendplot!(t_end)
tendplot(t_end1, t_end2)
tendplot!(t_end1, t_end2)

Plot recipe that adds a labeled vertical line to the plot at time t_end. A default label and styling are applied, but these can be modified by keyword arguments as usual If two time points are passed, a light shading is applied between instead of a vertical line.

tendplot(t_end)
tendplot!(t_end)
tendplot(t_end1, t_end2)
tendplot!(t_end1, t_end2)

Plot recipe that adds a labeled vertical line to the plot at time t_end. A default label and styling are applied, but these can be modified by keyword arguments as usual If two time points are passed, a light shading is applied between instead of a vertical line.

get_vial_mass(vialsize)

Return vial mass for given ISO vial size.

Uses a table from a SCHOTT manual, stored internally in a CSV.

get_vial_radii(vialsize)

Return inner and outer radius for passed ISO vial size.

Uses a table from a SCHOTT manual, stored internally in a CSV.

get_vial_shape(vialsize::String)

Return a NamedTuple with a slew of vial dimensions, useful for drawing the shape of the vial with make_outlines.

Uses a table from a SCHOTT manual, stored internally in a CSV.

get_vial_thickness(vialsize)

Return vial wall thickness for given ISO vial size.

Uses a table from a SCHOTT manual, stored internally in a CSV.

make_outlines(dims, Vfill)

Return a sequence of points (ready to be made into Plots.Shapes for the vial and fill volume, with Unitful dimensions, for given vial dimensions.

This is a convenience function for making figures illustrating fill depth.

lyo_1d_dae_f = ODEFunction(lyo_1d_dae!, mass_matrix=Diagonal([1.0, 0.0]))

Compute the right hand side function for the Pikal model.

The DAE system which is the Pikal model (1 ODE, one nonlinear algebraic equation for pseudosteady conditions) is here treated as a constant-mass-matrix implicit ODE system. The implementation is in lyo_1d_dae! and calc_md_Q.

The initial conditions u0 = [h_f, Tf] should be unitless, but are internally assigned to be in [cm, K]. The unitless time is taken to be in hours, so derivatives are given in unitless [cm/hr, K/hr].

params is a ParamObjPikal, which can be constructed with the following form (helping with readability):

params = ParamObjPikal((    (Rp, hf0, csolid, ρsolution),    (Kshf, Av, Ap),    (pch, Tsh) ,))

where those listed following are callables returning Quantitys, and the rest are Quantitys. See RpFormFit and RampedVariable for convenience types that can help with the callables.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)
• Kshf(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).
• Tsh(t), pch(t) return shelf temperature and chamber pressure respectively at time t.

struct ParamObjPikal{__T_Rp, __T_hf0, __T_csolid, __T_ρsolution, __T_Kshf, __T_Av, __T_Ap, __T_pch, __T_Tsh} <: LyoPronto.ParamObj

The ParamObjPikal type is a container for the parameters used in the Pikal model.

params is a ParamObjPikal, which can be constructed with the following form (helping with readability):

params = ParamObjPikal((    (Rp, hf0, csolid, ρsolution),    (Kshf, Av, Ap),    (pch, Tsh) ,))

where those listed following are callables returning Quantitys, and the rest are Quantitys. See RpFormFit and RampedVariable for convenience types that can help with the callables.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)
• Kshf(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).
• Tsh(t), pch(t) return shelf temperature and chamber pressure respectively at time t.

struct ParamObjRF{__T_Rp, __T_hf0, __T_csolid, __T_ρsolution, __T_Kshf, __T_Av, __T_Ap, __T_pch, __T_Tsh, __T_P_per_vial, __T_mf0, __T_cpf, __T_mv, __T_cpv, __T_f_RF, __T_eppf, __T_eppvw, __T_Kvwf, __T_Bf, __T_Bvw} <: LyoPronto.ParamObj

The ParamObjRF type is a container for the parameters used in the RF model.

Since it has many fields, the recommended constructor accepts a tuple of tuples, as follows, to help avoid ordering mistakes:

params = ParamObjRF((       (Rp, hf0, cSolid, ρsolution),    (Kshf_f, Av, Ap),    (pch, Tsh, P_per_vial),    (mf0, cpf, mv, cpv),    (f_RF, eppf, eppvw),    (Kvwf, Bf, Bvw),))

The parameters should all be Unitful quantities with appropriate dimensions, with some exceptions which are callables returning quantities. See RpFormFit and RampedVariable for convenience types that can help with these cases.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)

• Kshf_f(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).

• Tsh(t), pch(t), P_per_vial(t) return shelf temperature, chamber pressure, and microwave power respectively at time t.

• Arad and alpha were used in a prior version of the model, and are not used in the current version; they will be removed in a future version.

calc_hRp_T(po, efd; i)

For experimental conditions given by a po and experimental data given by efd, compute the effective $R_p$ and dry layer height $h_d$ over time.

Since po is a a ParamObjPikal, you will need to construct that object–the value of $R_p$ will not be used here, so set it to any dummy value.

If efd has multiple temperature series, pass an index i to select which series to use. Otherwise, the first series will be used.

calc_md_Q(u, po, t)

With the Pikal model, compute model quantities at u=[hf, Tf], time t, and conditions po.

po should be a ParamObjPikal; the return will be a named tuple with Unitful quantities in the following fields:

• md: mass flow rate (g/hr)
• Q_shf: heat transfer from shelf to product (W)

This allows assessment of the model's outputs without needing to rewrite the model equations.

calc_md_Q(u, po, tn)

Compute the mass flow and heat transfer terms for the lumped-capacitance microwave-assisted model.

Returns a named tuple with the following fields, all as Unitful quantities:

• md: mass flow rate (g/hr)
• Q_shf: heat transfer from shelf to product (W)
• Q_vwf: heat transfer from vial wall to product (W)
• Q_RF_f: volumetric heating of product (W)
• Q_RF_vw: volumetric heating of vial wall (W)
• Q_shw: heat transfer from shelf to vial wall (W)

lumped_cap_rf!(du, u, params, tn)

Compute the right-hand-side function for the ODEs making up the lumped-capacitance microwave-assisted model.

To access the values of the various heat transfer terms, use [calc_md_Q](@ref) to compute them; that function is used internally by this function.

du refers to [dmf/dt, dTf/dt, dTvw/dt], with u = [mf, Tf, Tvw]. u is taken without units but assumed to have the units of [g, K, K] (which is internally added). tn is assumed to be in hours (internally added), so dudt is returned with assumed units [g/hr, K/hr, K/hr] to be consistent.

Use the ParamObjRF type to hold the parameters.

params = ParamObjRF((       (Rp, hf0, cSolid, ρsolution),    (Kshf_f, Av, Ap),    (pch, Tsh, P_per_vial),    (mf0, cpf, mv, cpv),    (f_RF, eppf, eppvw),    (Kvwf, Bf, Bvw),))

The parameters should all be Unitful quantities with appropriate dimensions, with some exceptions which are callables returning quantities. See RpFormFit and RampedVariable for convenience types that can help with these cases.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)

• Kshf_f(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).

• Tsh(t), pch(t), P_per_vial(t) return shelf temperature, chamber pressure, and microwave power respectively at time t.

• Arad and alpha were used in a prior version of the model, and are not used in the current version; they will be removed in a future version.

lyo_1d_dae!(du, u, params, t)

Internal implementation of the Pikal model. See lyo_1d_dae_f for the wrapped version, which is more fully documented.

qrf_integrate(sol, RF_params)

Compute the integral over time of each heat transfer mode in the lumped capacitance model.

RF_params should represent the same parameters used to generate the solution sol, which (if OrdinaryDiffEq doesn't change) can likely be accessed as sol.prob.p.

Returns a NamedTuple with keys which match the result of calc_md_Q (Q_sub, Q_shf, Q_vwf, Q_RF_f, Q_RF_vw, Q_shw).

rf_lumcap_EM_violate(po)

Compute whether the given RF parameters will violate conservation of energy.

Specifically: calls calc_md_Q(calc_u0(po), po, 0.0), and returns true if the Q_RF_f + Q_RF_vw > po.P_per_vial(0.0). To be used with the badprms keyword argument to obj_pd and similar.

Molecular weight of water

Glass heat capacity Rough estimate from Bansal and Doremus 1985, "Handbook of Glass Properties"

Ice heat capacity IAPWS for ice: use triple point value for simplicity

Permittivity of free space

Glass thermal conductivity From Bansal and Doremus 1985, rough estimate based an a sodium borosilicate glass

Thermal conductivity of ice From Slack, 1980 200K : .032 W/cmK 250K : .024 W/cmK 273K : .0214 W/cmK

Thermal conductivity of sucrose, Caster grade powder From McCarthy and Fabre, 1989 book chapter

Density of ice at triple point.

Heat of sublimation of ice Approximate: coresponds to 254K or 225K From Feistel and Wagner, 2006

Dielectric loss coefficient of borosilicate glass, per Schott's testing

Water vapor viscosity in dilute limit

From Hellmann and Vogel, 2015 250K: 8.054 μPas 260K: 8.383 μPas 270K: 8.714 μPa*s

Glass density SCHOTT measurement, posted online

Density of sucrose, Caster grade powder From McCarthy and Fabre, 1989 book chapter

Stefan-Boltzmann Constant

Using the same Arrhenius fit as calc_psub, compute the sublimation temperature at a given pressure.

calc_psub(T) 
calc_psub(T::Q) where Q<:Quantity

Compute pressure (in Pascals) of sublimation at temperature T in Kelvin.

This is essentially an Arrhenius fit, where we compute: psub = pref * exp(-ΔHsub*Mw/R/T)

Single-purpose module for computing the dielectric loss coefficient of ice, as a function of temperature and frequency. This module exports a function ϵppf(T, f) for doing so. (Set as a separate module just to keep namespaces clean.) This code is a nearly-direct implementation of the correlation, Eqs. 3-6 presented in: Takeshi Matsuoka, Shuji Fujita, Shinji Mae; Effect of temperature on dielectric properties of ice in the range 5–39 GHz. J. Appl. Phys. 15 November 1996; 80 (10): 5884–5890. https://doi.org/10.1063/1.363582

ϵppf(T, f)

Compute the dielectric loss of ice as a function of temperature and frequency. Expects temperature in Kelvin and frequency in Hz, with Unitful units.

eq_cap_line(d, vt, l, volume)

Compute the equipment capability line for given geometry parameters. The resulting object can be called with a pressure in mTorr to get the corresponding mass flow rate in kg/hr. The line's slope and intercept can be accessed with line.k and line.b, respectively.

If provided with Unitful.Quantity inputs, the geometry parameters will be converted to the expected units (mm for lengths, m^3 for volume) before calculation, and the result will have Unitful units as well. If provided with plain floats, plain floats will be returned.

The input values are:

• d: Diameter of the duct [$mm$].
• vt: Thickness of the butterfly valve inside the duct [$mm$]
• l: Length of the duct [$mm$].
• v: Volume of the product chamber [$m^3$].

eq_cap_line_new(d, vt, l, volume)

Compute the equipment capability line for given geometry parameters using interpolation on the pressures at the four mass flow rate sample points.

The resulting object can be evaluated at pressure p, as a float in mTorr (if dimensions were given as plain floats) or with Unitful pressure (if dimensions were given as quantities) to get the corresponding mass flow rate in kg/hr. To invert a resulting line, access its slope with line.k (in kg/hr/mTorr) and its intercept with line.b (in kg/hr).

If provided with Unitful.Quantity inputs, the geometry parameters will be converted to the expected units (mm for lengths, m^3 for volume) before calculation, and the result will have Unitful units as well. If provided with plain floats, plain floats will be returned.

The input values are:

• d: Diameter of the duct [$mm$].
• vt: Thickness of the butterfly valve inside the duct [$mm$]
• l: Length of the duct [$mm$].
• v: Volume of the product chamber [$m^3$].

eq_cap_pressure(m, D, valve_thickness, L, volume)

Call eq_cap_line to get the equipment capability line for the given geometry parameters, then evaluate that line with mass flow rate m in kg/hr to get minimum controllable pressure in mTorr.

If provided with Unitful.Quantity inputs, the geometry parameters will be converted to the expected units (mm for lengths, m^3 for volume) before calculation, and the result will have Unitful units as well. If provided with plain floats, plain floats will be returned.

The input values are:

• d: Diameter of the duct [$mm$].
• vt: Thickness of the butterfly valve inside the duct [$mm$]
• l: Length of the duct [$mm$].
• v: Volume of the product chamber [$m^3$].

eq_cap_pressures_new(d, vt, l, volume)

Compute the pressures at the four mass flow rate sample points for given geometry parameters. The resulting vector is in the same order as M_DOT, and has units of mTorr

If provided with Unitful.Quantity inputs, the geometry parameters will be converted to the expected units (mm for lengths, m^3 for volume) before calculation, and the result will have Unitful units as well. If provided with plain floats, plain floats will be returned.

The input values are:

• d: Diameter of the duct [$mm$].
• vt: Thickness of the butterfly valve inside the duct [$mm$]
• l: Length of the duct [$mm$].
• v: Volume of the product chamber [$m^3$].

is_in_model_range(d, vt, l, volume)

Check whether the given geometry parameters are within the range of the original data used to generate the equipment capability curves.