Could it not be improved by, say, something like this:
newphase = phase + timestep * frequency
Where newphase is a local variable, probably the same type as phase, a member variable of the particle emitter, where timestep is the amount of time that has passed since the last physics update step, and where frequency is the number of particles emitted per second.
steps = floor(newphase) - floor(phase)
Where steps is a local integer variable.
oldpos = pos
pos = somefunc(oldpos)...
That is, compute the new position (pos, but don't throw away the old position (oldpos) just yet.
for each integer i on the range (phase, newphase]
{
frac = (i - phase) / (newphase - phase)
ppos = oldpos * (1 - frac) + frac * pos
emit_particle_at_pos(ppos)
}
Where emit_particle_at_pos is a stand-in for the particle emission call, using the argument ppos as the position of the emitted aprticle. Note that i may be equal to newphase, but cannot be equal to phase, thus ensuring that particles don't get emitted twice.
phase = newphase
Last but not least, update the phase.
If I've thought this through right, it would appear to evenly time the emission of the particles. However, using this technique might cause rapidly increasing numbers of particles, which could be bad for framerates.
The technique could also be fairly easily modified to create evenly spaced particles,
id est the distance between each particle, rather than the time, would be equal.