propNavMsg should handle {multiple PRN/consts} x {multiple epochs}

[ammtc]

Right now, propNavMsg requires the user to manually replicate a single timepoint of interest, epoch0, to equal the number of {PRN / constellation index} pairs:

prns = [1 2 10 19 1 3]'; consts = [1 1 2 3 2 1]';
epoch0 = navsu.time.cal2epochs(2020,2,15,13,30,0);
epochs = epoch0*ones(size(prns));

datai = navsu.geo.propNavMsg(beph,prns,consts,epochs);

At a minimum, we should be able to handle two obvious degenerate cases of the above:

1. Call propNavMsg with multiple PRN / constellation pairs, as above, but with a single epoch argument. In other words, the next-to-last line in the code block above should be automatically handled if size(epochs) == size(epoch0) == [1 1].

2. Call propNavMsg with a single PRN / constellation pair, but with a vector of epochs at which to calculate. In this case, we should calculate the results for the PRN in question at each requested timestep.

The third, most general case, would also be very useful — calculate the results for each PRN / constellation pair at each of the times in the epochs vector. However, this will require an extra loop in the propagation code, and a bit of careful thinking about to how to organize the outputs so as to maintain backwards compatibility, if at all possible, for those who have been using earlier versions of propNavMsg.