# -*- coding: utf-8 -*-
"""
Created on Tue Mar 11 13:56:52 2025

@author: regin
"""



import numpy as np
from scipy.linalg import eig
from scipy.special import iv,i0  # Modified Bessel function of the first kind
import matplotlib.pyplot as plt
from nice_figures import *
from uncertainties import ufloat, unumpy as unp


import pandas as pd

# Define constants (set to 1 for now)
#D_x = 1
beta=1
#a=2*np.pi
period=2*np.pi
Q_1 =1#2*np.pi/period
epsilon=10**(-5)
# Matrix size (truncation) 10 is quite low, 250 is quite high
num = 40
figsize = load_style()
cols = load_cols()
colors=[cols['blue'],cols['orange'],cols['green'],cols['red'],cols['purple'],cols['brown'],cols['pink'],cols['olive']]
w, h = figsize[('APS', '1-column')]


def L_matrix_Fokker_planck(num,q,u, f,D_x):
    
    L_q = np.zeros((2*num+1, 2*num+1), dtype=complex)
    for mu in range(-int(num), int(num)+1):
         for nu in range(-int(num), int(num)+1):
             delta = 1 if mu == nu else 0
             delta_p1 = 1 if mu == nu +  1 else 0
             delta_m1 = 1 if mu == nu - 1 else 0
             #L_q[mu+num,nu+num] =  -( D_x * mu**2 +  D_x *1j*f*mu )*delta-u/2*mu*(delta_p1-delta_m1)+ delta*(-1j*q*f-2*q*D_x*mu-D_x*q**2)-q*u/2*(delta_p1-delta_m1)
             #L_q[mu+num,nu+num] =  -( D_x * mu**2 +  D_x *1j*f*mu )*delta-D_x*u/2*mu*(delta_p1-delta_m1)+ D_x*delta*(-1j*q*f-2*q*Q_1*mu-q**2)-D_x*q*u/2*Q_1*(delta_p1-delta_m1)
             L_q[mu+num,nu+num] =  -Q_1**2*D_x*(u*mu/2*(delta_p1-delta_m1) + delta*(1j*mu*f+mu**2)) - D_x*Q_1*q*(u/2*(delta_p1-delta_m1)+ delta*(1j*f+2*mu))-D_x*q**2*delta

    D, W0,V0  = eig(L_q, left=True) #eigenvalue, left eigenvector, right eigenvector
    
    W = np.zeros((2*num+1, 2*num+1), dtype=complex)
    V = np.zeros((2*num+1, 2*num+1), dtype=complex)

    for i in range(-num,num+1):
        norm=np.sum(np.conjugate(W0[:,i])*V0[:,i])
        W[:,i]=W0[:,i]/np.conjugate(np.sqrt(norm))
        V[:,i]=V0[:,i]/np.sqrt(norm)
       # print( np.dot(np.conjugate(W[:,i]),  V[:,i]))
    
    D=-D    #(-1.) because of our definition: check np.dot(L,V0[:,50])+D[50]*V0[:,50] =0 (+ and not - works!, as in our definition of the eigenvalue equation)

    return D, W,V



 
def L_matrix_FP_kk(num,q,u, f,D_x):
    L_q = np.zeros((2*num+1, 2*num+1), dtype=complex)
    for mu in range(-int(num), int(num)+1):
         for nu in range(-int(num), int(num)+1):
             delta = 1 if mu == nu else 0
             L_q[mu+num,nu+num] = delta*(-D_x*q**2)

    return L_q


def L_matrix_FP_k(num,q,u, f,D_x):
    L_q = np.zeros((2*num+1, 2*num+1), dtype=complex)
    for mu in range(-int(num), int(num)+1):
         for nu in range(-int(num), int(num)+1):
             delta = 1 if mu == nu else 0
             delta_p1 = 1 if mu == nu +  1 else 0

             delta_m1 = 1 if mu == nu - 1 else 0
             L_q[mu+num,nu+num] =  - D_x*Q_1*q*(u/2*(delta_p1-delta_m1)+ delta*(1j*f+2*mu))  #+ delta*(1j*f-2*mu))-D_x*q**2*delta

    return L_q


def L_matrix_FP(num,q,u, f,D_x):
    L_q = np.zeros((2*num+1, 2*num+1), dtype=complex)

    for mu in range(-int(num), int(num)+1):
         for nu in range(-int(num), int(num)+1):

             delta = 1 if mu == nu else 0
             delta_p1 = 1 if mu == nu +  1 else 0
             delta_m1 = 1 if mu == nu - 1 else 0
             L_q[mu+num,nu+num] =  -Q_1**2*D_x*(u*mu/2*(delta_p1-delta_m1) + delta*(1j*mu*f+mu**2))
             #L_q[mu+num,nu+num] =  (- D_x * mu**2 - 1j*f*mu )*delta+u/2*mu*(delta_p1-delta_m1)

    return L_q


def bralLketr(l, n, M, r, m):
    return np.dot(np.conj(l[:,n]), M @ r[:,m])
def braLket_fill(matrix,u,f,D_x):
    braLket=np.zeros((2*num+1,2*num+1),dtype=np.complex_) 
    D, W, V=L_matrix_Fokker_planck(num,0,u, f,D_x)
    for index_l in range(-num,num+1):
        for index_r in range(-num,num+1):
            braLket[index_l][index_r] = bralLketr(W, index_l, matrix, V, index_r)
    return braLket

def kappa_1(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    n0=np.where(np.abs(D)<10**(-5))[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)
    #print(n0[0])
    braLket_x=braLket_fill(tot_q,u,f,D_x)
    return  np.real(1j*t*braLket_x[n0[0],n0[0]]) # the 2 comes from the fact 

def mu_1(t,u,f,D_x):
    return kappa_1(t,u,f,D_x)
def mu_tilde_1(t,u,f,D_x):
    return kappa_1(t,u,f,D_x)
def kappa_tilde_1(t,u,f,D_x):
    return kappa_1(t,u,f,D_x)

def kappa_2(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    n0=np.where(np.abs(D)<10**(-5))[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)
    braLket_x=braLket_fill(tot_q,u,f,D_x)
    varianz_s=0.0+0.0j

    for index,lamb in enumerate(D):
        if abs(lamb)!=0.0 :  
            varianz_s+=braLket_x[n0[0],index]*braLket_x[index,n0[0]]*(-1+np.exp((-lamb)*t)+(lamb)*t)/(lamb)**2 #only the terms linear in k contruibute, the others become k^3 or higher
    varianz_s-=t*D_x #?+t**2/2* bralLketr(W,nullentry,tot_k,V,nullentry)**2  #the quadratic k term 
    return  np.real(-2.*varianz_s) # the 2 comes from the fact 

def eff(u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    n0=np.where(np.abs(D)<10**(-5))[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)
    braLket_x=braLket_fill(tot_q,u,f,D_x)
    varianz_s=0.0+0.0j

    for index,lamb in enumerate(D):
        if abs(lamb)!=0.0 :  
            varianz_s+=braLket_x[n0[0],index]*braLket_x[index,n0[0]]*(1)/(lamb)**2 #only the terms linear in k contruibute, the others become k^3 or higher
    return  np.real(2.*varianz_s) # the 2 comes from the fact 


def D_infinity(u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    n0=np.where(np.abs(D)<10**(-5))[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)
    braLket_x=braLket_fill(tot_q,u,f,D_x)
    D_value=0.0+0.0j

    for index,lamb in enumerate(D):
        if abs(lamb)!=0.0 :  
            D_value-=braLket_x[n0[0],index]*braLket_x[index,n0[0]]/(lamb)#only the terms linear in k contruibute, the others become k^3 or higher
    D_value+=D_x #?+t**2/2* bralLketr(W,nullentry,tot_k,V,nullentry)**2  #the quadratic k term 
    return  np.real(D_value) # the 2 comes from the fact 


def kappa_tilde_2(t,u,f,D_x):
    return kappa_2(t,u,f,D_x)-2*D_x*t

def mu_2(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    tot_k=L_matrix_FP_k(num,1,u, f,D_x)
    #tot_kk=L_matrix_FP_kk(num,1,u, f,D_x)
    varianz_s=0.0+0.0j
    nullentry=np.where(abs(D)<epsilon)[0][0]
    for index,lamb in enumerate(D):
        if abs(lamb)!=0.0 :
            varianz_s+=bralLketr(W,nullentry,tot_k,V,index)*bralLketr(W,index,tot_k,V,nullentry)*(-1+np.exp((-lamb)*t)+(lamb)*t)/(lamb)**2 #only the terms linear in k contruibute, the others become k^3 or higher
    varianz_s+=-t*D_x  +t**2*bralLketr(W,nullentry,tot_k,V,nullentry)**2 /2 #the quadratic k term 
    return  -2.*varianz_s # the 2 comes from the fact 

def mu_tilde_2(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    tot_k=L_matrix_FP_k(num,1,u, f,D_x)
    tot_kk=L_matrix_FP_kk(num,1,u, f,D_x)
    
    varianz_s=0.0+0.0j
    nullentry=np.where(abs(D)<epsilon)[0][0]
    for index,lamb in enumerate(D):
        if abs(lamb)!=0.0 :
            varianz_s+=bralLketr(W,nullentry,tot_k,V,index)*bralLketr(W,index,tot_k,V,nullentry)*(-1+np.exp((-lamb)*t)+(lamb)*t)/(lamb)**2 #only the terms linear in k contruibute, the others become k^3 or higher
    varianz_s+=t**2*bralLketr(W,nullentry,tot_k,V,nullentry)**2 /2 #the quadratic k term 
    return  -2.*varianz_s # the 2 comes from the fact 

def compute_HA(q, Qmu, Qnu, t, u,f):
    D_x=1
    tau=1/(D_x*beta*np.sqrt(u**2-f**2)*Q_1**2)
    nu=(Qnu+q)
    mu=(Qmu+q)
    exponent =  -0.5 * D_x * (mu**2 - 2 * np.exp(-t / tau) * mu *nu + nu**2) * tau- 1j*period/2*(mu-nu)
    return np.exp(exponent)
def harm_variance(t,u,f,D_x):
    D_x=1
    tau=1/(D_x*beta*np.sqrt(u**2-f**2)*Q_1**2)
    return 2*tau*D_x*(1-np.exp(-t/tau))


def harm_D(t,u,f,D_x):
    D_x=1
    tau=1/(D_x*beta*np.sqrt(u**2-f**2)*Q_1**2)
    return D_x*(np.exp(-t/tau))
def deterministic_v(u,f,D_x):
    D_x=1
    return D_x*np.sqrt(f**2-u**2)*Q_1




def bands_real():
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$5$',r'$7$',r'$12$',r'$15$',r'$20$']
    #us=[0.0,2,4,6,8,10]
    mark=['s','o','^','v','d']
    lines=['-','--',':']
    fs=[0.1,0.6]    
    
    fig, ax = plt.subplots(figsize=(1.2* w, h))
    ax.text( 0.07, 0.95, "(a)", transform=ax.transAxes, ha="right", va="top", )
    for i, (mark_i, color,f,line) in enumerate(
        zip(mark, colors,fs,lines)):
        
        bandsre=[]

        sp=0.005
        for q in np.arange(-0.5,0.5,sp):
            D,W,V=L_matrix_Fokker_planck(num,q,1, f,1)
            D_sorted = D[np.argsort(D.imag)]
            bandsre.append(np.real((D_sorted)))

        if i==0:
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+2] for a in bandsre]),linestyle=line,label="2",color=colors[0])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+1] for a in bandsre]),linestyle=line,label="1",color=colors[1])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num] for a in bandsre]),linestyle=line,label="0",color=colors[2])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-1] for a in bandsre]),linestyle=line,label="-1",color=colors[3])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-2] for a in bandsre]),linestyle=line,label="-2",color=colors[4])
        else :
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+2] for a in bandsre]),linestyle=line,color=colors[0])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+1] for a in bandsre]),linestyle=line,color=colors[1])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num] for a in bandsre]),linestyle=line,color=colors[2])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-1] for a in bandsre]),linestyle=line,color=colors[3])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-2] for a in bandsre]),linestyle=line,color=colors[4])
        ticks = np.linspace(-np.pi, np.pi, 5)
        ax.set_xticks(ticks)
        ax.set_xticklabels([r"$-\pi$", r"$-\pi/2$", r"$0$", r"$\pi/2$", r"$\pi$"])
        ax.legend(title='n',frameon=True)
        ax.set_ylabel(r'Re$\lambda_{nq}L^2/D$')
        ax.set_xlabel(r'$qL$')
bands_real()
plt.tight_layout()
plt.savefig("figures/bands_re.pdf")
plt.show()  

def bands_imag():
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$5$',r'$7$',r'$12$',r'$15$',r'$20$']
    #us=[0.0,2,4,6,8,10]
    mark=['s','o','^','v','d']
    lines=['-','--',':']
    fs=[0.1,0.6]    
    
    fig, ax = plt.subplots(figsize=(1.2* w, h))
    ax.text( 0.07, 0.95, "(b)", transform=ax.transAxes, ha="right", va="top", )
    for i, (mark_i, color,f,line) in enumerate(
        zip(mark, colors,fs,lines)):
        
        bandsre=[]

        sp=0.005
        for q in np.arange(-0.5,0.5,sp):
            D,W,V=L_matrix_Fokker_planck(num,q,1, f,1)
            D_sorted = D[np.argsort(D.imag)]
            bandsre.append(np.imag((D_sorted)))

        if i==0:
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+2] for a in bandsre]),linestyle=line,label="2",color=colors[0])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+1] for a in bandsre]),linestyle=line,label="1",color=colors[1])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num] for a in bandsre]),linestyle=line,label="0",color=colors[2])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-1] for a in bandsre]),linestyle=line,label="-1",color=colors[3])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-2] for a in bandsre]),linestyle=line,label="-2",color=colors[4])
        else :
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+2] for a in bandsre]),linestyle=line,color=colors[0])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num+1] for a in bandsre]),linestyle=line,color=colors[1])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num] for a in bandsre]),linestyle=line,color=colors[2])    
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-1] for a in bandsre]),linestyle=line,color=colors[3])
            ax.plot(np.arange(-0.5,0.5,sp)*2*np.pi,4*np.pi**2*np.array([a[num-2] for a in bandsre]),linestyle=line,color=colors[4])
        ticks = np.linspace(-np.pi, np.pi, 5)
        ax.set_xticks(ticks)
        ax.set_xticklabels([r"$-\pi$", r"$-\pi/2$", r"$0$", r"$\pi/2$", r"$\pi$"])
        ax.legend(title='n')
        ax.set_ylabel(r'Im$\lambda_{nq}L^2/D$')
        ax.set_xlabel(r'$qL$')
bands_imag()
plt.tight_layout()
plt.savefig("figures/bands_im.pdf")
plt.show()  






import glob
from pathlib import Path
files = glob.glob("data_observables_variable_name/msd_u*_f*.txt")
import re

def parse_params(fname):
    m = re.search(r"u([0-9.]+)_f([0-9.]+)", fname)
    return float(m.group(1)), float(m.group(2))

parse_params('data_observables_variable_name/msd_u100_f201.txt')

BASE = Path("data_observables_variable_name")
def filename(u, f,path):
    return BASE/Path(path) / f"msd_u{u:.2f}_f{f:.2f}.txt"
#obs  = np.loadtxt(filename(100, 100,'plot_variance_force')) 


def plot_velocity_force():

    D_x=1
    kbTs=[0.01,0.1,1.0]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$100$',r'$10$',r'$1$',r'$12$',r'$15$',r'$20$']
   
    #mark=['s','o','^','v','d','p']
    mark=['s','o','^']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    
    
    us=1
    f= np.logspace(0, 2, 100)
    plotbox2=[]
    for j in range(len(f)):
        det_v_here=deterministic_v(us,f[j],D_x)
        plotbox2.append(det_v_here/(D_x/period))
    ax.plot(f/us,np.array(plotbox2)/us, color='black',linestyle=':', label=r'$v_{\mathrm{det}}$')   

    for i, (mark_i, color) in enumerate(zip(mark, colors)):
        kbT=kbTs[i]
        #f=np.linspace(0.5, 1.5,60)/kbT
        us=1/kbT
        f= np.logspace(-1, 2, 100)*us
        
        plotbox=[]
        for j in range(len(f)):
            #t=np.array([1000,1001])
            #Dhere=np.gradient(kappa_2(t,us,f[j],D_x) ,t)[0] / 2
            Kappa_here=kappa_1(1,us,f[j],D_x)
            plotbox.append(Kappa_here/(D_x/period))
        #print(plotbox)
        ax.plot(f/us,np.array(plotbox)/us, 
        color=color,
        #marker=mark_i,
        markerfacecolor='None',
        markeredgecolor=color,
        label=labels[i]
        )

        f_sim=np.array([0.2,0.56,0.8,1.0,3.16,17.78,80.])*us
        plotbox3=[]
        for j in range(len(f_sim)):
            obs  = np.loadtxt(filename(us, f_sim[j],'plot_velocity_force')) 
            det_v_sim=np.mean((np.array(obs[:,1])/obs[:,0])[-13:-4])
            plotbox3.append(det_v_sim)
        ax.plot(f_sim/us, np.array(plotbox3)/us,marker='.',linestyle='None',color='black',markersize=4)


                            
    f_help= np.logspace(0.5, 1.5, 50)
    ax.plot(f_help, 4*f_help,color='black' )
    ax.text(0.65, 0.6, r"$\propto f$", transform=ax.transAxes, fontsize=14)
        
    ax.set_xlabel(r"$f/u$")
    ax.set_ylabel(r"$(v/u)L/D $")
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_ylim(0.05, 1000)
    ax.legend(title=r"$u$",frameon=False)
'''
plot_velocity_force()
plt.tight_layout()
plt.savefig("figures/plot_velocity.pdf")

plt.show()
'''
def plot_variance_potential():
    #t=np.logspace(-4, 3, 400) #define the times for the Variance
    #D_x=1
    #f=100
    #us=[0.0,40,100,113,117,130]
    

           

    t=np.logspace(-3, 4, 400) #define the times for the Variance 
    D_x=1
    f=0.001
    us=[5]
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$5$',r'$7$',r'$12$',r'$15$',r'$20$']
    #us=[0.0,2,4,6,8,10]
    mark=['s','o','^','v','d']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
        #ax.axhline(eff(us[i],f,D_x)/period**2,color=color,markeredgecolor=color)
        ax.plot(t*D_x/period**2,harm_variance(t,us[i],f,D_x)/period**2)
        ax.plot(t*D_x/period**2,
                kappa_2(t,us[i],f,D_x)/period**2,
                #mu_2(t,us[i],f,D_x), 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i],
                linestyle='-'
                )

        ax.plot(t[:200]*D_x/period**2,2*(D_x)*t[:200]/period**2,color='black',linestyle='-')
        if us[i]<7:
            ax.plot(t[200:]*D_x/period**2,2*(D_infinity(us[i],f,D_x))*t[200:]/period**2,color='black',linestyle=':')
        ax.plot(t*D_x/period**2,2*(D_infinity(us[i],f,D_x))*t/period**2+eff(us[i],f,D_x)/period**2,color='black',linestyle='--')
        ax.axhline(eff(us[i],f,D_x)/period**2,color='black',linestyle=':')      
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$\mathrm{Var}(t)/L^2$")
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$U_1/k_bT$",frameon=False)
'''
plot_variance_potential()
plt.tight_layout()
#plt.savefig("figures/plot_variance_potential_eff_f=0.pdf")
plt.show()
'''
def plot_variance_potential_eff():
    #t=np.logspace(-4, 3, 400) #define the times for the Variance
    #D_x=1
    #f=100
    #us=[0.0,40,100,113,117,130]
    
    t=np.logspace(-3, 3, 400) #define the times for the Variance 
    D_x=1
    f=3
    us=[5,8,10]
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$5$',r'$7$',r'$12$',r'$15$',r'$20$']
    #us=[0.0,2,4,6,8,10]
    mark=['s','o','^']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
        #ax.axhline(eff(us[i],f,D_x)/period**2,color=color,markeredgecolor=color)

        ax.plot(t*D_x/period**2,
                kappa_2(t,us[i],f,D_x)/period**2,
                #mu_2(t,us[i],f,D_x), 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i],
                linestyle='-'
                )
        ax.plot(t*D_x/period**2,harm_variance(t,us[i],f,D_x)/period**2, color=color, linestyle=':')
        ax.plot(t*D_x/period**2,(eff(us[i],f,D_x)*t/t)/period**2, color=color, linestyle='-.')        
        #ax.plot(t[:200]*D_x/period**2,2*(D_x)*t[:200]/period**2,color='black',linestyle='-')
        #if us[i]<7:
        #    ax.plot(t[200:]*D_x/period**2,2*(D_infinity(us[i],f,D_x))*t[200:]/period**2,color='black',linestyle=':')
        ##ax.plot(t*D_x/period**2,2*(harm_variance(t,us[i],f,D_x)+D_infinity(us[i],f,D_x))*t/period**2,color=color,linestyle='--')
        #ax.axhline(eff(us[i],f,D_x)/period**2,color='black',linestyle=':')      
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$\mathrm{Var}(t)/L^2$")
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$U_1/k_bT$",frameon=False)

'''
plot_variance_potential_eff()
plt.tight_layout()
#plt.savefig("figures/plot_variance_potential_eff_f_neq0.pdf")
plt.show()
'''
def kramer(f,u,Dx):
    return np.exp(2*u*np.sqrt(1 - (f/u)**2) + 2*f*np.arcsin(f/u) - f*np.pi) \
            / ((2 * np.pi) * u * np.sqrt(1 - (f/u)**2)) \
            * (1 / (1 + np.exp(-2 * np.pi * f)))
def kramer2a(f,u,Dx):        
    return (1 / (2 * np.pi * u * np.sqrt(1 - (f/u)**2)) \
            * np.exp((np.sqrt(2)*4/3) * u * (1 - (f/u))**(3/2)))
def kramer2b(f,u,Dx):        
    return (1 / (2 * np.pi * u * np.sqrt(1 - (f/u)**2)) \
            * np.exp((2/3) * u * (1 - (f/u)**2)**(3/2)))
                
        
        
def periodtime(f,u,Dx):        
    return period/Dx/Q_1/np.sqrt(f**2-u**2)/period**2*Dx    
def tautime(f,u,Dx):
    return 1 /(Dx*np.sqrt(u**2-f**2)* Q_1**2)/period**2*Dx

def plot_variance_force():
    t=np.logspace(-5, 3, 400) #define the times for the Variance
    Dx=1

    #u=1
    #fs=[0,0.1,0.5,1,1.5,4]
    
    u=100 # U/kbT=1/0.01
    
    #fs=np.array([0,8,9.6,10,15,20])
    fs=np.array([0,85,90,100,110,200])
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$0.0$',r'$0.85$',r'$0.9$',r'$1.0$',r'$1.1$',r'$2.0$']
    #labels=[r'$0$',r'$0.8$',r'$0.96$',r'$1.0$',r'$1.5$',r'$2.0$']
    mark=['s','o','^','v','d','p']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    ax.text( 0.07, 0.97, "(a)", transform=ax.transAxes, ha="right", va="top", )
   
    for i, (mark_i, color) in enumerate(zip(mark, colors)):

        ax.loglog(t/period**2*Dx,
                kappa_2(t,u,fs[i],Dx)/period**2,
                #scnd_mom(t,u,fs[i],D_x), 
                color=color,linestyle='-',
                #marker=mark_i,
                markerfacecolor='None',markeredgecolor=color,label=labels[i]
                )

        obs  = np.loadtxt(filename(u, fs[i],'plot_variance_force')) 
        if fs[i]==0:
            ax.plot(obs[:,0][::2][:-9]  , (obs[:,2]-obs[:,1]**2)[::2][:-9],marker='.',linestyle='None',color='black',markersize=4)
        else:
            ax.plot(obs[:,0][::2][:-1]  , (obs[:,2]-obs[:,1]**2)[::2][:-1],marker='.',linestyle='None',color='black',markersize=4)
        #ax.axvline(kramer(fs[i],u,Dx),color=color,linestyle='-.')
        #ax.axvline(tautime(fs[i],u,Dx),color=color, linestyle=':')
    for i, (mark_i, color) in enumerate(zip(mark, colors)):
        if fs[i]<u:
            ax.plot(t*Dx/period**2,harm_variance(t,u,fs[i],Dx)/period**2, color='black', linestyle=':')
            
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$\mathrm{Var}(t)/L^2$")
    #ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$f/u$",frameon=False)
'''
plot_variance_force()
plt.tight_layout()
#plt.savefig("figures/plot_variance_force.pdf")
plt.show()
'''
#1 / (2* np.pi * u * np.sqrt(1 - (f/u)**2)) \        * np.exp((1/3) * u * (1 - (f/u)**2)**(3/2))

#np.exp(2*u*np.sqrt(1 - (f/u)**2) + 2*f*np.arcsin(f/u) - f*np.pi) \  / ((2 * np.pi) * u * np.sqrt(1 - (f/u)**2)) \ * (1 / (1 + np.exp(-2 * np.pi * f)))
        
def plot_diffusion_potential():
    t=np.logspace(-3, 2, 400) #define the times for the Variance
    Dx=1
    #f=1
    #us=[0.1,0.5,1,2,2.5,3]
    f=10
    us=[0,6.0,10,12,15,20.0]
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$0$',r'$6$',r'$10$',r'$12$',r'$15$',r'$20$']
   
    mark=['s','o','^','v','d','p']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
                Dhere=np.gradient(kappa_2(t,us[i],f,Dx) /(period)**2,t*D_x/(period)**2) / 2
                ax.plot(t*D_x/(period)**2,Dhere, 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i]
                )
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$D(t)/D$")
    #ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$u$",frameon=False)
''''
plot_diffusion_potential()
plt.tight_layout()
plt.savefig("figures/diffusion_potentials.pdf")
plt.show()

'''


def plot_diffusion():
    t=np.logspace(-3, 2, 400) #define the times for the Variance
    D_x=1.0
    u=10
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$0$',r'$5$',r'$7$',r'$10$',r'$15$',r'$20$']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    #colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    fs=[0,5,7,10,15,20]
    mark=['s','o','^','v','d','p']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
        Dhere=np.gradient(kappa_2(t,u,fs[i],D_x) /(period)**2,t*D_x/(period)**2) / 2
        ax.plot(t*D_x/(period)**2,
                Dhere/D_x, 
                #scnd_mom(t,u,fs[i],D_x), 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i]
                )
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$D(t)/D$")
    #ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$f$",frameon=False)
'''
plot_diffusion()
plt.tight_layout()
plt.savefig("figures/diffusion.pdf")
plt.show()

'''

def mean_derivative(t, var, discard_last=5, use_before=10):

    t = np.asarray(t)
    var = np.asarray(var)

    N = len(t)
    # Fenster definieren
    start = N - discard_last - use_before - 1
    end = N - discard_last
    t_win = t[start:end]
    var_win = var[start:end]

    # diskrete Ableitung
    dvar_dt = np.diff(var_win) / np.diff(t_win)

    # Mittelwert der Ableitung
    return np.mean(dvar_dt)


def plot_giant_diffusion2():
    
    #f=1
    #us=[0.1,0.5,,2,2.5,3]
    D_x=1

    #us=1
    #us=[0.001,6.0,10,12,15,20.0]

    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    kbT=1
    us=[100,10,1.0]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$100$',r'$10$',r'$1$',r'$12$',r'$15$',r'$20$']
   
    #mark=['s','o','^','v','d','p']
    mark=['s','o','^']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    ax.text( 0.07, 0.97, "(b)", transform=ax.transAxes, ha="right", va="top", )
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
        
                f=np.linspace(0.5, 1.5,200)*us[i]
                #f= np.logspace(-1.5, 0.5, 50)*us[i]

                plotbox=[]
                for j in range(len(f)):
                    #t=np.array([1000,1001])
                    #Dhere=np.gradient(kappa_2(t,us,f[j],D_x) ,t)[0] / 2
                    Dhere=D_infinity(us[i],f[j],D_x)
                    #Dhere=1
                    plotbox.append(Dhere/D_x)
                #print(plotbox)
                ax.plot(f/us[i],np.array(plotbox), 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                linestyle='-',
                label=labels[i]
                )
                f=np.linspace(0.5, 1.5,5)*us[i]
                plotbox2=[]
                for j in range(len(f)):

                    obs  = np.loadtxt(filename(us[i], f[j],'giant_diffusion')) 
                    giantdiff=mean_derivative((obs[:,0]), (obs[:,2]-obs[:,1]**2), 5, 20)/2
                    plotbox2.append(giantdiff/D_x)
                ax.plot(f/us[i] , np.array(plotbox2),marker='.',linestyle='None',color='black',markersize=4)

                #print(i)
    ax.set_xlabel(r"$f/u$")
    ax.set_ylabel(r"$D_\infty/D$")
    ax.set_ylim(-0.1,19)
    #ax.set_yscale('log')
    #ax.set_xscale('log')
    ax.legend(title=r"$u$",frameon=False)

plot_giant_diffusion2()
plt.tight_layout()
plt.savefig("figures/giant_diffusion.pdf")
plt.show()

def kappa_3(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0,u, f,D_x)

    n0=np.where(np.abs(D)<10**(-5))[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)

    braLket_x_q=braLket_fill(tot_q,u,f,D_x)
  
    moment3=0+0j
    for index_l,lamb in enumerate(D):
        if abs(lamb)>epsilon:
            function3_0=((np.exp(-t* lamb)*2.0+np.exp(-t* lamb)*t * lamb+ (-2.0+t *lamb)))/lamb**3
            #function3_1=(-2.0-2.0*np.exp(-t*lamb)-2.0* t*lamb+t**2 *lamb**2)/(2.0*lamb**3)
            moment3+=function3_0*braLket_x_q[n0,index_l]*braLket_x_q[index_l,n0]*(braLket_x_q[index_l,index_l]-braLket_x_q[n0,n0])            
            #moment3+=2*function3_1*braLket_x_q[n0,index_l]*braLket_x_q[index_l,n0]*braLket_x_q[n0,n0]
            #print(braLket_x_q[index_l,index_l])
            for index_e,eta in enumerate(D):
                if abs(eta)>epsilon and abs(eta-lamb)>epsilon:
                    function3_3=(np.exp(-t*lamb)+t* lamb-1) /((lamb**2)* (eta-lamb) )+(np.exp(-t*eta)+t *eta-1)/((eta**2)* (lamb-eta))
                    moment3+=function3_3*braLket_x_q[n0,index_l]*braLket_x_q[index_l,index_e]*braLket_x_q[index_e,n0]
    return (-6j)*moment3


def mu_3(t,u,f,D_x):
    k1=kappa_1(t,u,f,D_x)
    moment3=kappa_3(t,u,f,D_x)+3*kappa_2(t,u,f,D_x)*k1+k1**3         
    return moment3

 
def plot_skewness():
    t=np.logspace(-3, 2.5, 400) #define the times for the Variance
    D_x=1.0
    u=10
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$0$',r'$0.8$',r'$0.95$',r'$1.0$',r'$1.5$',r'$2.0$']
    fs=[0.0,8,9.5,10,15,20]
    #u=1
    #fs=[0.0001,0.5,1,2,3,10]
    mark=['s','o','^','v','d','p']
    fig, ax = plt.subplots(figsize=(1.2 * w, h))

    
    ax.text( 0.07, 0.97, "(a)", transform=ax.transAxes, ha="right", va="top", )
    for i, (mark_i, color) in enumerate(zip(mark, colors)):
        
        ax.semilogx(t/period**2*D_x,                
                kappa_3(t,u,fs[i],D_x)/kappa_2(t,u,fs[i],D_x)**(3/2), 
                #mu_2(t,u,fs[i],D_x), 
                color=color,
                #marker=mark_i,
                linestyle='-',
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i]
                )
    for i, (mark_i, color) in enumerate(zip(mark, colors)):
        
        obs  = np.loadtxt(filename(u, fs[i],'skewness')) 
        obs_k3=obs[:,3]-3*obs[:,2]*obs[:,1]+2*obs[:,1]**3
        obs_k2=obs[:,2]-obs[:,1]**2
        if fs[i]<8:
            ax.plot(obs[:,0][::][1:-20]  , (obs_k3/obs_k2**(3/2))[::][1:-20],marker='.',linestyle='None',color='black',markersize=4)
        else:
            ax.plot(obs[:,0][::][:-10]  , (obs_k3/obs_k2**(3/2))[::][:-10],marker='.',linestyle='None',color='black',markersize=4)
            

    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$\mathrm{Skew}(t)$")
    #ax.set_yscale('log')
    ax.set_xscale('log')
    ax.legend(title=r"$f/u$",frameon=False)
    ax.set_ylim(-0.9,2.1)

plot_skewness()
plt.tight_layout()
#plt.savefig("figures/skewness.pdf")
plt.show()


def function4_(t,lamb,eta,xi):
    return (-1+np.exp(-lamb*t)+(lamb)*t)/(lamb)**2/(lamb-eta)/(lamb-xi)+(-1+np.exp((-eta)*t)+(eta)*t)/(eta)**2/(eta-lamb)/(eta-xi)+(-1+np.exp((-xi)*t)+(xi)*t)/(xi)**2/(xi-eta)/(xi-lamb)


def mu_tilde_4(t,u,f,D_x):
    D, W, V=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    n0=np.where(D==0)[0]
    tot_q=L_matrix_FP_k(num,1,u, f,D_x)
    
    braLket_x=braLket_fill(tot_q,u,f,D_x)
    
    moment4=t**4/24*braLket_x[n0[0],n0[0]]**4
    print('here', moment4)
    #moment4=0+0j
    for index_l,lamb in enumerate(D):
        if index_l!=n0[0]:
            function4_0=(-6+6 *np.exp(-t* lamb)+6* t* lamb-3 *t**2 *lamb**2+t**3* lamb**3)/(6 *lamb**4)
            function4_1=(6.0- 4 *t *lamb + t**2* lamb**2 - 2* np.exp(-t *lamb)* (3+t* lamb))/(2 *lamb**4)
            function4_2=((np.exp(-t* lamb)*6.0+np.exp(-t* lamb)*4 *t* lamb+np.exp(-t* lamb)*t**2* lamb**2+2 *(-3+t* lamb)))/(2* lamb**4)
           
            moment4+=3*function4_0*braLket_x[n0[0],index_l]*braLket_x[index_l,n0[0]]*braLket_x[n0[0],n0[0]]*braLket_x[n0[0],n0[0]]
            moment4+=2*function4_1*braLket_x[n0[0],index_l]*braLket_x[index_l,index_l]*braLket_x[index_l,n0[0]]*braLket_x[n0[0],n0[0]]
            moment4+=function4_1*braLket_x[n0[0],index_l]*braLket_x[index_l,n0[0]]*braLket_x[n0[0],index_l]*braLket_x[index_l,n0[0]]
            #print(braLket_x[index_l,n0[0]])
            moment4+=function4_2*braLket_x[n0[0],index_l]*braLket_x[index_l,index_l]*braLket_x[index_l,index_l]*braLket_x[index_l,n0[0]]
            for index_e,eta in enumerate(D):
                if abs(eta)>epsilon and abs(eta-lamb)>epsilon:
                    function4_3=(-1 + np.exp(-t*eta) +   t*eta)/(eta**3*(eta - lamb)) + t**2/( 2*eta*lamb) + (-1 + np.exp(-t* lamb) +   t* lamb)/(lamb**3* (-eta + lamb))
                    function4_6a=((np.exp(-t*(eta))*lamb**3+(eta-lamb)**2*(-lamb+eta*(-2+t*lamb))+np.exp(-t*(lamb))*eta**2*(eta*(2+t*lamb)-lamb*(3+t*lamb))))/(eta**2*(eta-lamb)**2*lamb**3)
                    function4_6b=((np.exp(-t*(lamb))*eta**3+(lamb-eta)**2*(-eta+lamb*(-2+t*eta))+np.exp(-t*(eta))*lamb**2*(lamb*(2+t*eta)-eta*(3+t*eta))))/(lamb**2*(lamb-eta)**2*eta**3)
             
                    moment4+=2*function4_3*braLket_x[n0[0],index_l]*braLket_x[index_l,index_e]*braLket_x[index_e,n0[0]]*braLket_x[n0[0],n0[0]]
                    moment4+=function4_3*braLket_x[n0[0],index_l]*braLket_x[index_l,n0[0]]*braLket_x[n0[0],index_e]*braLket_x[index_e,n0[0]]
                    moment4+=function4_6b*braLket_x[n0[0],index_l]*braLket_x[index_l,index_e]*braLket_x[index_e,index_e]*braLket_x[index_e,n0[0]]
                    moment4+=function4_6a*braLket_x[n0[0],index_l]*braLket_x[index_l,index_e]*braLket_x[index_e,index_l]*braLket_x[index_l,n0[0]]
                    moment4+=function4_6a*braLket_x[n0[0],index_l]*braLket_x[index_l,index_l]*braLket_x[index_l,index_e]*braLket_x[index_e,n0[0]]
                    
                    for index_x,xi in enumerate(D):
                        if abs(xi)>epsilon and abs(xi-lamb)>epsilon and abs(eta-xi)>epsilon:
                            moment4+=function4_(t,lamb,eta,xi)*braLket_x[n0[0],index_l]*braLket_x[index_l,index_e]*braLket_x[index_e,index_x]*braLket_x[index_x,n0[0]]

    return 24*moment4

def mu_4(t,u,f,D_x):
    return mu_tilde_4(t,u,f,D_x)+12*D_x*t*mu_tilde_2(t,u,f,D_x)+12*D_x**2*t**2

def kappa_4(t,u,f,D_x):
    mu1=kappa_1(t,u,f,D_x)
    mu2=np.real(mu_2(t,u,f,D_x))
    mu3=np.real(mu_3(t,u,f,D_x))
    mu4=np.real(mu_4(t,u,f,D_x))
    return mu4-4*mu3*mu1-3*mu2**2+12*mu2*mu1**2-6*mu1**4

def plot_kumulant_3():
    t=np.logspace(-3, 3, 400) #define the times for the Variance
    D_x=1.0
    u=1
    #labels=[r'$1$',r'$2$',r'$3$']
    #mark=['s','o','^']
    #colors=['r', '#009900', 'b']
    #fs=[1,10,15,20.0,25.0,45.0,35.0,50]
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    labels=[r'$0$',r'$10$',r'$15$',r'$20$',r'$30$',r'$50$']
    fs=[0.001,1,2,3,4,10]
    mark=['s','o','^','v','d','p']
    Fig, ax = plt.subplots()
    for i, (mark_i, color) in enumerate(
        zip(mark, colors)):
        ax.semilogx(t,
                kappa_3(t,u,fs[i],D_x)/kappa_2(t,u,fs[i],D_x)**(3/2), 
                #mu_2(t,u,fs[i],D_x), 
                color=color,
                #marker=mark_i,
                markerfacecolor='None',
                markeredgecolor=color,
                label=labels[i]
                )
    
'''
plot_kumulant_3()
plt.legend()    
plt.show()
'''


'''

f_values=[0.1,5,8,10,15,20]
data = []
header = "t"  # Header for the file
#f_values = [0.001,0.1,0.5,1,2]
f_values=[0.1,5,8,10,15,20]
t=np.logspace(-3, 2, 300)
for f in f_values:
    u=10
    D_x=1
    NGP=kappa_4(t,u,f,D_x)/3/kappa_2(t,u,f,D_x)**2
    data.append(np.real(NGP))
    header += f"\tf={f}"  # Add column header for each V_0
    plt.semilogx(t,NGP)
# Combine time and ISF data
plt.show()

data = np.column_stack([t] + data)
# Save to .dat file
np.savetxt("NG_data.dat", data, header=header, fmt="%.6e", delimiter="\t")
# Save to .txt file (same content)
np.savetxt("NG_data.txt", data, header=header, fmt="%.6e", delimiter="\t")
'''
#analytic data
data_analytic = np.loadtxt("NG_data.txt", skiprows=1)  # Load data, skipping header
time_analytic, *NG_analytic = data_analytic.T 


def plot_kurtosis():
    u=10
    f_values=[0.1,5,8,10,15,20]
    fig, ax = plt.subplots(figsize=(1.2 * w, h))
    ax.text( 0.07, 0.97, "(b)", transform=ax.transAxes, ha="right", va="top", )
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    for i, (f, NG,color) in enumerate(zip(f_values, NG_analytic,colors)):    
        if i>1:
            ax.plot(time_analytic/period**2,
                    NG, 
                    #mu_2(t,u,fs[i],D_x), 
                    color=color,
                    #marker=mark_i,
                    linestyle='-',
                    markerfacecolor='None',
                    markeredgecolor=color,
                    label=f'{f*0.1:.1f}'
                    )
    for i, (f, NG,color) in enumerate(zip(f_values, NG_analytic,colors)):    
        if i>1:       
            obs  = np.loadtxt(filename(u, f,'skewness')) 
            
            obs_k4=obs[:,4]-4*obs[:,3]*obs[:,1]-3*obs[:,2]**2+12*obs[:,2]*obs[:,1]**2-6*obs[:,1]**4
            obs_k2=obs[:,2]-obs[:,1]**2
            
            if f==15:
                ax.plot(obs[:,0][::][0:-22]  ,(obs_k4/3/obs_k2**(2))[::][0:-22],marker='.',linestyle='None',color='black',markersize=4)
            else:
                ax.plot(obs[:,0][::][1:-12]  ,(obs_k4/3/obs_k2**(2))[::][1:-12],marker='.',linestyle='None',color='black',markersize=4)
                
    ax.set_xlabel(r"$tD/L^2$")
    ax.set_ylabel(r"$\alpha_2(t)$")
    #ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_ylim(-0.5,1.7)
    ax.legend(title=r"$f/u$",frameon=False)

plot_kurtosis()
plt.tight_layout()
#plt.savefig("figures/kurtosis.pdf")
plt.show()

def compute_ISF(mu, nu, t,q,u,f):

    D_x=1
    q=q*2*np.pi/period
    
    eigenvalues, eigenvectorsl, eigenvectorsr=L_matrix_Fokker_planck(num,q,u, f,D_x)
    #eigenvalues0, eigenvectorsl0, eigenvectorsr0=L_matrix_Fokker_planck(num, 0.00001,u, f,D_x)
    eigenvalues0, eigenvectorsl0, eigenvectorsr0=L_matrix_Fokker_planck(num,0.0,u, f,D_x)
    null=np.where(np.abs(eigenvalues0)<0.00001)[0]
    result = 0
    for n, lambda_nq in enumerate(eigenvalues): #!! inserted minus sign
        exp_term = np.exp(-lambda_nq * t)
        sum_chi = 0
        sum_xi = 0
        #print(n)
        #np.dot(eigenvectors[:, n],eigenvectors0[:,null[0]])

        for chi in range(-num,num+1):
            if nu>0:
                sum_chi += (eigenvectorsr[chi-nu+num, n] *np.conjugate(eigenvectorsl0[chi+num,null[0]]))
            else:
                sum_chi += (eigenvectorsr[chi+num, n] *np.conjugate(eigenvectorsl0[chi+nu+num,null[0]]))

        for xi in range(-num,num+1):
            if mu>0:
                sum_xi += (np.conjugate(eigenvectorsl[xi-mu+num, n]) *eigenvectorsr0[xi+num,null[0]])
            else:
                sum_xi += (np.conjugate(eigenvectorsl[xi+num, n]) *eigenvectorsr0[xi+mu+num,null[0]])
        result += exp_term * sum_chi * sum_xi
    #vt=kappa_1(t,u,f,D_x)
    vt=0
    #np.exp(1j*q*vt)
    #np.exp(1j*period/2*(mu-nu))
    #np.sqrt((mu+q)**2+(nu+q)**2)
    return (result)





def plot_ISF_real_imag_grid(
    cases=((1, 1, 0.0), (1, 1, 0.5), (3, 3, 0.5)),
    u=10,
    f_ratios=(0.0, 0.8, 1.0, 1.5,2.0),
):
    t = np.logspace(-3, 3, 400)
    t_scaled = t / period**2
    Dx=1
    fig, axes = plt.subplots(nrows=2,ncols=3,figsize=(3 * 1.2 * w, 2 * h),sharex=True,)
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    panel_labels = ["(a)", "(b)", "(c)", "(d)", "(e)", "(f)"]
    for col, (mu, nu, q) in enumerate(cases):
        ax_re = axes[0, col]
        ax_im = axes[1, col]

        for f, c in zip(f_ratios, colors):
            ISF = compute_ISF(-mu, -nu, t, q, u, f * u)
            ax_re.plot(t_scaled, np.real(ISF), color=c, label=rf"${f}$",linestyle='-')
            ax_im.plot(t_scaled, np.imag(ISF), color=c, label=rf"${f}$",linestyle='-')
            

        ax_re.axvline(kramer( 0.8 * u,u,Dx),color=colors[1],linestyle='--')
        ax_re.axvline(tautime(0.8 * u,u,Dx),color='darkgrey', linestyle='--')
        ax_im.axvline(kramer( 0.8 * u,u,Dx),color=colors[1],linestyle='--')
        ax_im.axvline(tautime(0.8 * u,u,Dx),color='darkgrey', linestyle='--')
        # Harmonic approximation (real part only)
        for f in (0.0, 0.8):
            ax_re.plot( t_scaled,np.real(compute_HA(q, mu, nu, t, u, f * u)),color="black",linestyle=":",)

        ax_re.set_xscale("log")
        ax_im.set_xscale("log")
        ax_re.set_ylim(-0.7,1.1)
        ax_im.set_ylim(-0.85,0.6)
        ax_re.set_xlim(0.00002,1)
        ax_im.set_xlim(0.00002,1)
        if q==0.0:
            ax_re.set_title(rf"$\mu=\nu={mu},\,qL=0$")
        if q==0.5:
            ax_re.set_title(rf"$\mu=\nu={mu},\,qL=\pi$")


    label_fs = 13
    title_fs = 13
    tick_fs = 12
    legend_fs = 11
    
    # axis labels
    axes[0, 0].set_ylabel(r"$\mathrm{Re}\,F_{\mu,\mu}(q,t)$", fontsize=label_fs)
    axes[1, 0].set_ylabel(r"$\mathrm{Im}\,F_{\mu,\mu}(q,t)$", fontsize=label_fs)
    
    for ax in axes[1, :]:
        ax.set_xlabel(r"$tD/L^2$", fontsize=label_fs)
    
    # titles, ticks, legends
    for ax in axes.flatten():
        ax.title.set_fontsize(title_fs)
        ax.tick_params(axis="both", labelsize=tick_fs)
        ax.legend(title=r"$f/u$", frameon=False, fontsize=legend_fs, title_fontsize=legend_fs)
    


    for ax, label in zip(axes.flatten(), panel_labels):
        ax.text( 0.1, 0.90, label, transform=ax.transAxes, ha="right", va="top", fontsize=label_fs )

        
    for o in range(3):
        for f, c in zip(f_ratios, colors):      
            obs  = np.loadtxt(filename(u, f * u,'ISF')) 
            axes[0, o].plot(obs[:,0][::][0:-6]  , obs[:,5+o][::][0:-6]  ,marker='.',linestyle='None',color='black',markersize=4)
        for f, c in zip(f_ratios, colors): 
            obs  = np.loadtxt(filename(u, f * u,'ISF')) 
            axes[1, o].plot(obs[:,0][::][0:-6]   , obs[:,11+o][::][0:-6]  ,marker='.',linestyle='None',color='black',markersize=4)
                   
    
    fig.tight_layout()
    return fig, axes
'''
fig, axes = plot_ISF_real_imag_grid()
#plt.savefig("figures/ISF_diagonal.pdf")
plt.show()

'''

def plot_ISF_real_imag_grid_offdiag(
    #cases=((2, 1, 0.0), (0, 1, 0.5), (0, -1, 0.5)),
    cases=((2, 1, 0.0), (0, 1, 0.5), (0, -1, 0.5)),
    u=10,
    f_ratios=(0.0, 0.8, 1.0, 1.5, 2.0),
):
    t = np.logspace(-3, 3, 400)
    t_scaled = t / period**2
    Dx=1
    fig, axes = plt.subplots(nrows=2,ncols=3,figsize=(3 * 1.2 * w, 2 * h),sharex=True,)
    colors=['r', '#009900', 'b', '#bf00ff','#ffc40c','#e6550d']
    panel_labels = ["(a)", "(b)", "(c)", "(d)", "(e)", "(f)"]
    for col, (mu, nu, q) in enumerate(cases):
        ax_re = axes[0, col]
        ax_im = axes[1, col]

        for f, c in zip(f_ratios, colors):
            ISF = compute_ISF(-nu, -mu, t, q, u, f * u)
            ax_re.plot(t_scaled, np.real(ISF), color=c, label=rf"${f}$",linestyle='-')
            ax_im.plot(t_scaled, np.imag(ISF), color=c, label=rf"${f}$",linestyle='-')
        ax_re.axvline(kramer( 0.8 * u,u,Dx),color=colors[1],linestyle='--')
        ax_re.axvline(tautime(0.8 * u,u,Dx),color='darkgrey', linestyle='--')
        ax_im.axvline(kramer( 0.8 * u,u,Dx),color=colors[1],linestyle='--')
        ax_im.axvline(tautime(0.8 * u,u,Dx),color='darkgrey', linestyle='--')
        # Harmonic approximation (real part only)
        for f in [0.0, 0.8]:
            ax_re.plot( t_scaled,np.real(compute_HA(q, mu, nu, t, u, f * u)),color="black",linestyle=":",)

        ax_re.set_xscale("log")
        ax_im.set_xscale("log")
        ax_re.set_ylim(-1.1,0.5)
        ax_im.set_ylim(-1,0.8)
        ax_re.set_xlim(0.0001,1)
        ax_im.set_xlim(0.0001,1)
        if q==0.0:
            ax_re.set_title(rf"$\mu=\nu={mu},\,qL=0$")
        if q==0.5:
            ax_re.set_title(rf"$\mu=\nu={mu},\,qL=\pi$")







    label_fs = 13
    title_fs = 13
    tick_fs = 12
    legend_fs = 11
    
  
    # titles, ticks, legends
    for ax in axes.flatten():
        ax.title.set_fontsize(title_fs)
        ax.tick_params(axis="both", labelsize=tick_fs)
    

    # Axis labels
    axes[0, 0].set_ylabel(r"$\mathrm{Re}\,F_{\mu,\nu}(q,t)$", fontsize=label_fs)
    axes[1, 0].set_ylabel(r"$\mathrm{Im}\,F_{\mu,\nu}(q,t)$", fontsize=label_fs)

    for ax in axes[1, :]:
        ax.set_xlabel(r"$tD/L^2$", fontsize=label_fs)

    for ax, label in zip(axes.flatten(), panel_labels):
        ax.text( 0.95, 0.95, label, transform=ax.transAxes, ha="right", va="top", fontsize=label_fs )
        
    axes[1, 0].legend(title=r"$f/u$", frameon=False,loc=3, fontsize=legend_fs, title_fontsize=legend_fs)
    axes[1, 1].legend(title=r"$f/u$", frameon=False,loc=4, fontsize=legend_fs, title_fontsize=legend_fs)
    axes[1, 2].legend(title=r"$f/u$", frameon=False,loc=3, fontsize=legend_fs, title_fontsize=legend_fs)

    for o in range(3):
        for f, c in zip(f_ratios, colors):      
            obs  = np.loadtxt(filename(u, f * u,'ISF')) 
            axes[0, o].plot(obs[:,0][::][0:-6]  , obs[:,8+o][::][0:-6]  ,marker='.',linestyle='None',color='black',markersize=4)
        for f, c in zip(f_ratios, colors): 
            obs  = np.loadtxt(filename(u, f * u,'ISF')) 
            axes[1, o].plot(obs[:,0][::][0:-6]   , obs[:,14+o][::][0:-6]  ,marker='.',linestyle='None',color='black',markersize=4)
              
    fig.tight_layout()
    return fig, axes
'''
fig, axes = plot_ISF_real_imag_grid_offdiag()
plt.savefig("figures/ISF_offdiagonal.pdf")
plt.show()
'''