import numpy as np
import math
import matplotlib.pyplot as plt

# - Described at: https://mathworld.wolfram.com/CircleMap.html
#   - Both the formula and modulo 1 trace to the above.
# - Described https://en.wikipedia.org/wiki/Arnold_tongue
#   - The above uses "+ K/2*pi" while MathWorld uses "- K/2*pi"
#   - The choice of omega = 1/3 traces to the above.
#   - The choice of kInterval of (0, 4 * math.pi) traces to the above, except that we start at -2
#     to get a more complete picture.
#   - Issue: the choice of omega and kInterval trace only to Wikipedia, not to other source.
# - The choice of starting at 0.5 is not based on anything; it is a guess taken over
#   from the logistic map.

preset = 1

saveFigure = True
markerSize = 0.001
markerStyle = "," # . - point; , - pixel
if preset == 1:
  kInterval = (-2, 4 * math.pi)
  kResolution = 0.002
  mapIterationCount = 1000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0.5
  yLim = (-0.1, 1.1)
  omega = 1.0/3.0
  fileName = "CircleMapOrbitDiag.png"

lims = np.zeros(mapIterationCount)

fig, biax = plt.subplots()
fig.set_size_inches(16, 9)

itemCount = 0
for K in np.arange(kInterval[0], kInterval[1], kResolution):
  itemCount += 1
  lims[0] = startingIterationValue
  for i in xrange(mapIterationCount - 1):
    lims[i + 1] = (lims[i] + omega + K / (2 * math.pi) * math.sin(2 * math.pi * lims[i])) % 1
  if itemCount % 1000 == 0:
    print("Parameter values processed:", itemCount)

  biax.plot([K] * lastMapIterationCountToPlot,
            lims[(mapIterationCount - lastMapIterationCountToPlot):], "b.",
            markersize=markerSize,
            marker=markerStyle)

biax.set(xlabel="K", ylabel="x",
         title="Circle map - (x + omega - K / (2 * pi) * sin(2 * pi * x)) mod 1 - bifurcation diagram/orbit diagram" + "\n" +
               "K resolution: " + str(kResolution) +
               ", iteration count: " + str(mapIterationCount) +
               ", last values to plot count: " + str(lastMapIterationCountToPlot) +
               (", omega: %.3f" % omega))
               
biax.set_ylim(yLim[0], yLim[1])
if saveFigure:
  plt.savefig(fileName, dpi=600)
else:
  plt.show()