How to identify a bimodal distribution? I understand that once we plot the values as a chart, we can identify a bimodal distribution by observing the twin-peaks, but how does one find it programmatically? (I am looking for an algorithm.)
 A: Identifying a mode for a continuous distribution requires smoothing or binning the data.
Binning is typically too procrustean: the results often depend on where you place the bin cutpoints.
Kernel smoothing (specifically, in the form of kernel density estimation) is a good choice.  Although many kernel shapes are possible, typically the result does not depend much on the shape.  It depends on the kernel bandwidth.  Thus, people either use an adaptive kernel smooth or conduct a sequence of kernel smooths for varying fixed bandwidths in order to check the stability of the modes that are identified.  Although using an adaptive or "optimal" smoother is attractive, be aware that most (all?) of these are designed to achieve a balance between precision and average accuracy: they are not designed to optimize estimation of the location of modes.
As far as implementation goes, kernel smoothers locally shift and scale a predetermined function to fit the data.  Provided that this basic function is differentiable--Gaussians are a good choice because you can differentiate them as many times as you like--then all you have to do is replace it by its derivative to obtain the derivative of the smooth.  Then it's simply a matter of applying a standard zero-finding procedure to detect and test the critical points.  (Brent's method works well.)  Of course you can do the same trick with the second derivative to get a quick test of whether any critical point is a local maximum--that is, a mode.
For more details and working code (in R) please see https://stats.stackexchange.com/a/428083/919.
A: I came late to the party, but if you are just interested in whether it is multimodal or not, meaning you are not interested in the number of modes, you should look at diptest.
In R the package is called diptest.
A: There is a well-known paper by Silverman that deals with this issue. It employs kernel-density estimation. See

B. W. Silverman, Using kernel density estimates to investigate multimodality, J. Royal
  Stat. Soc. B, vol. 43, no. 1, 1981, pp. 97-99.

Note that there are some errors in the tables of the paper. This is just a starting point, but a pretty good one. It provides a well-defined algorithm to use, in the event that's what you're most looking for. You might look on Google Scholar at papers that cite it for more "modern" approaches.
A: The definition in wiki is slightly confusing to me.  The probability of a continous data set having just one mode is zero.  A simple way to program a bimodal distrubiton is with two seperate normal distributions centered differently.  This creates two peaks or what wiki calls modes.  You can actually use almost any two distributions, but one of the harder statistical opportunities is to find how the data set was formed after combining the two random data distributions.
