For a couple years now, I have been using OrdinaryDiffEq for developing a handful of models for describing the process of pharmaceutical lyophilization. The more time I spend, the more I learn about more aspects of Julia or packages that would have made my life easier if I had known about them sooner (e.g. Parameters.jl, ComponentArrays.jl). In that light, I am wondering if my current headache already has a good solution.

These models have anywhere between 10 and 30 physical properties that need to be defined for simulation. These range from universal physical properties (e.g. universal gas constant) to things that don’t often need tinkering (e.g. density of glass) to things that are different for every simulation (input temperature/pressure profiles). As I compare these models to experiment and refine them, I periodically find that I need to add temperature-dependence (or pressure-dependence, or frequency-dependence) to a parameter that I was previously treating as constant.

Up til now, my approach has been to rewrite the code in such a way that, instead of taking a constant value, the ODE function takes a callable (or anonymous function) which it evaluates everywhere it needs the property. This works, of course. But every time it happens I inevitably end up re-running the code several times finding all the places where I am now adding an anonymous function to a number, and I feel like it clutters the input script, e.g. turning

```
p = 10.0u"Pa"
```

into

```
p = t->10.0u"Pa"
```

which, in addition to being a little uglier, also makes error messages not as nice and confuses the non-Julia-users when I need to hand off my code to someone else.

I have also experimented with creating a custom type with pretty printing, with a single-argument constructor that essentially acts as a nicer-looking version of the anonymous function.

```
julia> p = RampedVariable(10.0u"Pa")
RampedVariable(10.0u"Pa")
julia> p(50u"s")
10.0u"Pa"
julia> p(Inf)
10.0u"Pa"
```

I implemented that type as below, for any who are curious. But I find that this method still ends up being more typing, and means I need to document more interfaces myself, etc.

My question is: is there a more elegant solution to this problem than, say, implementing every single parameter every time as a callable which may not actually depend on any of its inputs?

For reference, sometimes I am thinking of time-dependent external controls, weakly temperature-dependent parameters like thermal conductivity, or things like the dielectric loss coefficient which varies with temperature and electric field frequency.

```
struct RampedVariable
setpts::AbstractVector
ramprates::AbstractVector
holds::AbstractVector
timestops::AbstractVector
end
function RampedVariable(initial)
RampedVariable([initial], [], [], [0])
end
function RampedVariable(setpts, ramprate)
if length(ramprate) == 0 || length(setpts) == 1
@error "If no ramp necessary, construct RampedVariable with only one argument." ramprate
end
if length(ramprate) >= 2 || length(setpts) > 2
@error "For multiple ramps, need at least one hold time. Construct RampedVariable with three arguments." ramprate
end
if length(setpts) != 2
@error "Number of set points should be 1 more than ramps, since initial is included"
end
timestops = fill(0.0*setpts[1]/ramprate[1], 2)
timestops[2] = timestops[1] + (setpts[2]-setpts[1])/ramprate
RampedVariable(setpts, [ramprate], [], timestops)
end
function RampedVariable(setpts, ramprates, holds)
if (length(ramprates) != length(holds) + 1 ) || (length(ramprates)==0)
@error "Number of ramps should be zero or number of holds + 1"
end
if length(setpts) != length(ramprates) + 1
@error "Number of set points should be 1 more than ramps, since initial is included"
end
timestops = fill(0.0*setpts[1]/ramprates[1], length(ramprates) + length(holds) + 1)
(ramp, rest) = Iterators.peel(ramprates)
timestops[2] = timestops[1] + (setpts[2]-setpts[1])/ramp
for (i, ramp) in enumerate(rest)
timestops[2i+1] = timestops[2i] + holds[i]
timestops[2i+2] = timestops[2i+1] + (setpts[i+2]-setpts[i+1])/ramp
end
RampedVariable(setpts, ramprates, holds, timestops)
end
function (rv::RampedVariable)(t)
if length(rv.timestops) == 1
return rv.setpts[1]
end
im = findlast(rv.timestops .<= t)
if im == length(rv.timestops)
return rv.setpts[end]
elseif iseven(im)
return rv.setpts[im÷2+1]
else
ip = im+1
return (rv.setpts[ip÷2+1] - rv.setpts[ip÷2])/(rv.timestops[ip] - rv.timestops[im])*(t - rv.timestops[im]) + rv.setpts[ip÷2]
end
end
```