function qout = qfun(parameters,t)
%Default water supply rate at nodes
%Input arguemtns:
%   parameters: this paramteers structure passed to network_percolation.
%   qfun uses fields
%       qscaled:    flag that determines which of two functional forms are
%                   used. If qscaled is false, qout = q_in_0. If qscaled is
%                   true,
%                   q_in = q_in_0*(q_0_old+(q_0-q_0_old*min((t-t_in)/t_q,1)))
%                   where t_start is the start time of the simulation
%       q_0, q_0_old, t_in, t_q: parameters activated by qscaled as above
%       oscillate:  if qscaled is set to 'true', a second flag 'oscillate'
%                   may be set whe set, an additional oscillatory term is
%                   added
%                   q_in = q_in_0*(q_0_old+(q_0-q_0_old*min((t-t_in)/t_q,1))) + 
%                           q_in_0*delta_q*cos(omega*t)*min((t-t_in)/t_q,1);
%       delta_q, omega, t_q: additional parameters for oscillatory water
%                   supply
%       moulins:    optional substructure, used only if qscaled is set to
%                   true. Allows localized water supply to be defined and
%                   phased in from a spatially smooth background defined
%                   through q_in_0 (qscaled could be repurposed for this by
%                   using q_0 and q_0_old as vectors rather than q_in_0
%                   'moulins' must have the following fields
%           active: flag that determines whether moulin input is set. If
%                   false, no other fields are needed. If true, need the
%                   following
%           base:   n_nodes-by-1 list of mean water supply levels through
%                   moulins
%           base_old: n_nodes-by-1 list of previous mean water supply levels
%                   through moulins
%           t_moulin: time scale for transition
%           oscillate: flag to add oscillatory term to moulin water supply
%           amplitude: amplitude of oscillatory supply
%           phase_moulin:      phase shift
%           envelope:   amplitude of envelope oscillation
%           phase_envelope: phase of envelope oscillation
%           omega_moulin:   angular frequency of oscillation
%           omega_envelope: angular frequency of envelope oscillation
%                   Water supply then works as
%                       q_in= q_in_background + base_old+(base-base_old)*min((t-t_in)/t_moulin,1)
%                       + sum(amplitude*(1+envelope*cos(omega_envelope*t+phase_envelope))*cos(omega_moulin*t+phase_moulin)*min((t-t_in)/t_moulin,1));
%                   where q_in_background is computed as described under
%                   qscaled/qoscillate above. All parameters from oscillate
%                   to phase_envelope above are n_nodes-by-m arrays (to
%                   account for multiple frequencies) and the sum is taken
%                   over the second index
%Output
%   qout:   n_nodes-by-1 vector giving instantaneous water input to each
%           node at time t


if parameters.qscaled
    qout = parameters.q_in_0*(parameters.q_0_old+(parameters.q_0-parameters.q_0_old)*min((t-parameters.t_in)/parameters.t_q,1));
    %water supply ramp-up
    if isfield(parameters,'oscillate') && parameters.oscillate
        qout = qout + parameters.q_in_0*parameters.delta_q*cos(parameters.omega*t)*min((t-parameters.t_in)/parameters.t_q,1);
    end
    if isfield(parameters,'moulins') && parameters.moulins.active
        qout = qout + parameters.moulins.base_old + (parameters.moulins.base-parameters.moulins.base_old)*min((t-parameters.t_in)/parameters.moulins.t_moulin,1);
        %moulin supply ramp-up
        if parameters.moulins.oscillate
            qout = qout + sum(parameters.moulins.amplitude.*(1+parameters.moulins.envelope.*cos(parameters.moulins.omega_envelope*t+parameters.moulins.phase_envelope))...
                .*cos(parameters.moulins.omega_moulin*t+parameters.moulins.phase_moulin),2)*min((t-parameters.t_in)/parameters.moulins.t_moulin,1);
        end
        %oscillatory supply to moulins
    end
else
    qout = parameters.q_in_0;
end

end

