r_Gittergröße = 150
#r_Gittergröße = 90
r_Gittergröße = 100
phi_Gittergröße = 6*r_Gittergröße

assert phi_Gittergröße % 2 == 0

Geschwindigkeit = 4.5 # Verhältnis Grintsch : Weihnachtsmann

SchrittweiteW = 0.1
SchrittweiteG = SchrittweiteW * Geschwindigkeit


from math import sin,cos,pi

Masche_phi = 2*pi/phi_Gittergröße
Masche_r = 1/r_Gittergröße

SchrG = 10 # Grintsch bekommt genaue Gitterschrittweite
SchrittweiteG = SchrG * Masche_phi * 1.00000001
SchrittweiteW = SchrittweiteG / Geschwindigkeit

name = "Grintsch-%s-%sx%s-%s" % (Geschwindigkeit,r_Gittergröße,phi_Gittergröße,
                              SchrittweiteW)
#name = "".join("," if c=="." else c for c in name) # matplotlib mag keine Punkte
print("Name:",name)
print ("Grintschgeschwindigkeit =",Geschwindigkeit,
       "Gitter =",r_Gittergröße,"x",phi_Gittergröße)

Bereich_phi_G = int(SchrittweiteG//Masche_phi)
Bereich_r_W = int(SchrittweiteW//Masche_r)

"""
Punkt [i,j] entspricht (r*cos(phi), r*sin(phi))
mit r=(1-i/r_Gittergröße) und phi=(j*Masche_phi)
Grintsch ist auf (1,0).
"""

optG = {}
optW = {}

for j in range(phi_Gittergröße):
    for i in range(r_Gittergröße):
        optW[i,j]=0,0
        optG[i,j]=0,0
    phi = j*Masche_phi
    optW[0,j]=min(phi, 2*pi-phi),0

max_var_j = {}
for i in range(1,r_Gittergröße):
    Bereich_i = range(max(0,i-Bereich_r_W),min(r_Gittergröße,1+ i+Bereich_r_W))
    r1 = 1-i*Masche_r # p1 = (r1,0)
    "j2 kann von j um höchstens max_var_j[i,i2] in beiden Richtungen abweichen."
    for i2 in Bereich_i:
        r2 = 1-i2*Masche_r
        for j2 in range(phi_Gittergröße//2):
            phi = j2*Masche_phi
            p2x,p2y = r2*cos(phi),r2*sin(phi)
            Abst2 = (p2x-r1)**2 + p2y**2
            if Abst2 > SchrittweiteW**2:
                max_var_j[i,i2] = j2-1
                break
        else:
            max_var_j[i,i2] = phi_Gittergröße//2
#print("MV",max_var_j)

def backtrack_solution_W(i,j):
    Punkte_W = []
    Punkte_G = []
    j_G = 0
    _,num_Schritte = optW[i,j]
    while True:
        r = 1-i*Masche_r
        w,Schritte = optW[i,j]
        phi = (j+j_G)*Masche_phi
        Punkte_W.append((r*cos(phi),r*sin(phi)))
        phi = j_G*Masche_phi
        rr = 1+Schritte*0.1/num_Schritte
        Punkte_G.append((rr*cos(phi),rr*sin(phi)))
        if i==0: break
        i,j,_ = finde_max_W(i,j)
        if i==0: break
        j2,_ = finde_min_G(i,j)
        j_G -= (j2-j)
        j = j2
    return Punkte_W,Punkte_G
        
def finde_max_W(i,j):
    Wert,Schritte = -1,0
    for i2 in range(max(0,i-Bereich_r_W),min(r_Gittergröße,1+ i+Bereich_r_W)):
        j_Intervall = (j - max_var_j[i,i2], j + max_var_j[i,i2])
        if j - max_var_j[i,i2] <= phi_Gittergröße//2 <= j + max_var_j[i,i2]:
            j_Intervall = (phi_Gittergröße//2,)
#        for j2 in range(j - max_var_j[i,i2], j + max_var_j[i,i2]+1):
        for j2 in j_Intervall:
            jj2 = j2%phi_Gittergröße
            nWert,nSchritte = optG[i2,jj2]
            if (nWert,-nSchritte)>(Wert,-Schritte):
                Wert,Schritte = nWert,nSchritte
                i_opt,j_opt = i2,jj2
    return i_opt,j_opt, (Wert,Schritte+1)
        
def iterate_W():
    Änderung = False
    for i in range(1,r_Gittergröße):
        for j in range(phi_Gittergröße):
            i_opt,j_opt, neu = finde_max_W(i,j)
            if optW[i,j] != neu:
                optW[i,j] = neu
                Änderung = True
    return Änderung
            
def finde_min_G(i,j):
    if j<=Bereich_phi_G or j+Bereich_phi_G>=phi_Gittergröße:
        j_opt = 0
    elif 2*j>phi_Gittergröße:
        j_opt = j+Bereich_phi_G
    else:
        j_opt = j-Bereich_phi_G
    assert(0<=j_opt<phi_Gittergröße)
    Wert,Schritte = optW[i,j_opt]
    return j_opt,(Wert,Schritte+1)

    # altes langsameres Programm
    for j2 in range(j - Bereich_phi_G, j + Bereich_phi_G+1):
        jj2 = j2%phi_Gittergröße
        nWert,nSchritte = optW[i,jj2]
        if (nWert,-nSchritte)<(Wert,-Schritte):
            Wert,Schritte = nWert,nSchritte
            j_opt = jj2
    return j_opt,(Wert,Schritte+1)

def iterate_G():
    Änderung = False
    for i in range(r_Gittergröße):
        for j in range(phi_Gittergröße):
            _, neu = finde_min_G(i,j)
            if optG[i,j] != neu:
                optG[i,j] = neu
                Änderung = True
    return Änderung

eps = 1e-8  # sin und cos mathematisch gleich, aber nicht exakt gleich.
def prüfe_Monotonie(T):
    for j in range(phi_Gittergröße):
        for i in range(1,r_Gittergröße-1): # i=0 ist eine Ausnahme
            assert T[i,j][0]>=T[i+1,j][0],((i,j),(i+1,j),T[i,j],"C>",T[i+1,j])
    for i in range(r_Gittergröße):
        for j in range(phi_Gittergröße//2):
            assert T[i,j][0]<=T[i,j+1][0],(i,j,T[i,j],"A>",T[i,j+1])
            j2 = phi_Gittergröße-j-1
            j2p = (j2+1)%phi_Gittergröße
            assert T[i,j2p][0]<=T[i,j2][0],(i,j2,T[i,j2p],"B>",T[i,j2])

def printtable(s,T):
    return
    print(s)
    for i in range(r_Gittergröße):
        print(i,"".join(" --  " if T[i,j] is None else
                        "%5.2f/%d"% T[i,j] for j in range(phi_Gittergröße)))

def iterate():
    num_iter = 0
    while True:
        num_iter += 1
        print("It%d"%num_iter,"%.4f/%d"%optW[r_Gittergröße-1,0],end=" ",flush=True)
        ÄnderungW = iterate_W()
        #printtable("W",optW)
        ÄnderungG = iterate_G()
        #printtable("G",optG)
        if num_iter%10==0: print()
        if not (ÄnderungW or ÄnderungG):
            print()
            return num_iter

iterate()

Punkte_W,Punkte_G = backtrack_solution_W(r_Gittergröße-1,0)

print ("Grintschgeschwindigkeit =",Geschwindigkeit)
print ("#Schritte r =",r_Gittergröße, " #Schritte phi =",phi_Gittergröße)
print ("Masche_phi",Masche_phi, "   Masche_r",Masche_r)
print ("Bereich_phi_G",Bereich_phi_G, "  Bereich_r_W",Bereich_r_W)


printtable("W",optW)
prüfe_Monotonie(optW)
printtable("G",optG)
prüfe_Monotonie(optG)

print ("Lösung: Abstand %s, %d Schritte" % optW[r_Gittergröße-1,0])
print ("Fertig.")


def map_table(T,nam=""):
    pts = ((r*cos(phi),r*sin(phi))
           for i in range(r_Gittergröße)
           for j in range(phi_Gittergröße)
           for r,phi in [(1-i*Masche_r, j*Masche_phi)])
    x,y = zip(*pts)
    values = list(T[i,j][0]
                  for i in range(r_Gittergröße)
                  for j in range(phi_Gittergröße))
    colors = list(viridis (v/pi) for v in values)
#    pcm = ax.pcolor(x, y, values, cmap='viridis', shading='auto')
#    fig.colorbar(pcm, ax=ax[1], extend='max')

    fig, ax = plt.subplots(num=name+nam+".pdf")

    ax.axis('equal')
    ax.scatter(x,y, color=colors, # cmap="viridis") #,
               s=1)
    plt.show(block=False)

import matplotlib.pyplot as plt

from matplotlib import cm
viridis = cm.get_cmap('viridis')
#def color_map(v):
##    v /= pi
#    return viridis(v) #(1-v,0.1,v)

map_table(optW,"W")

fig, ax = plt.subplots(num=name+".pdf")

x,y=zip(*Punkte_W)
ax.plot(x, y, 'go-', markersize=2.3)
x,y=zip(*Punkte_G)
ax.plot(x, y, 'bo-', markersize=2.3,linewidth=0.5)
ax.grid()
ax.axis('equal')
plt.show()