# -*- coding: utf-8 -*-
"""
Created on Mon May  3 13:48:06 2021

@author: Regina
"""


from matplotlib.legend_handler import HandlerTuple


import numpy as np
import matplotlib.pyplot as plt
from uncertainties import unumpy as unp
import statistics
import matplotlib.ticker as ticker
from scipy.stats import moment
plt.rcParams.update({
  "text.usetex": True,
})

size=20.0

plt.rcParams['figure.figsize'] = (6.4, 4.8) 
plt.rcParams['axes.labelsize'] =size
plt.rcParams['xtick.labelsize'] =size
plt.rcParams['ytick.labelsize'] =size
plt.rcParams['xtick.direction'] = "in"
plt.rcParams['ytick.direction'] = "in"
plt.rcParams['legend.fontsize'] =size
plt.rcParams['legend.title_fontsize'] =size
plt.rcParams['legend.frameon'] = False
plt.rcParams['legend.labelspacing'] = 0.05
plt.rcParams['lines.linewidth'] = 2.0
plt.rcParams['lines.markersize'] = 5.0
plt.rcParams['font.family'] = 'lmodern' 

def data_msd(data_path):

    f1 = open(data_path,'r')
    lines1 = f1.readlines()
    f1.close()

    i, j, count_block, time_here, msd_here, msd_here_nc, msd_here_nc_cc, msd_here_isf, msd_here_isf_nc, msd_here_isf_nc_erst, msd_here_isf_nc_kreuz = ([] for _ in range(11))
    
    for line in range(0,len(lines1[1:])):
            p=lines1[1:][line].split()
            if line%10!=0 or line==0:
                i.append(float(p[0]))
                j.append(float(p[1]))
                count_block.append(float(p[2]))
                time_here.append(float(p[3]))
                msd_here.append(float(p[4]))
                msd_here_nc.append(float(p[5]))
                msd_here_nc_cc.append(float(p[6]))
                msd_here_isf.append(float(p[7]))

    return i,j,count_block,time_here,msd_here,msd_here_nc,msd_here_nc_cc,msd_here_isf



def meanstd_msd(num,pfad,pfad2):
    #ddof=1 to get unbias std means with 1/n-1 *sum over sample (NOT the whole smple but just random subset)
    start=0

    y=np.array([data_msd(pfad+str(i)+pfad2) for i in range(start,num)])

    time_here_mean=np.mean( np.array([y[i][3] for i in range(0,num-start)]), axis=0 )

    arr=np.array([y[i][4] for i in range(0,num-start)])
    msd_here_mean=np.mean(arr, axis=0 )
    msd_here_std=np.std(arr, axis=0, ddof=1)/np.sqrt(num)


    arr_nc= np.array([y[i][5] for i in range(0,num-start)])
    msd_here_nc_mean=np.mean(arr_nc, axis=0 )
    msd_here_nc_std=np.std(arr_nc, axis=0, ddof=1)/np.sqrt(num)


    arr_cc=np.array([y[i][6] for i in range(0,num-start)])
    msd_here_nc_cc_mean=np.mean( arr_cc, axis=0 )
    msd_here_nc_cc_std=np.std( arr_cc, axis=0,ddof=1 )/np.sqrt(num)

    msd_here_isf_mean=np.mean( np.array([y[i][7] for i in range(0,num-start)]), axis=0)
    msd_here_isf_std=np.std( np.array([y[i][7] for i in range(0,num-start)]), axis=0, ddof=1 )/np.sqrt(num)


    return time_here_mean,msd_here_mean,msd_here_std,msd_here_nc_mean,msd_here_nc_std,msd_here_nc_cc_mean,msd_here_nc_cc_std,msd_here_isf_mean,msd_here_isf_std #specialfälle


samplenumber_5=100
time_01,msd_mean_01,msd_std_01,msd_nc_mean_01,msd_nc_std_01,msd_nc_cc_mean_01,msd_nc_cc_std_01,msd_isf_mean_01,msd_isf_std_01=meanstd_msd(samplenumber_5,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code/Data_observables_periodic_sinus_potential_V0_01/Data_','/data_observables_periodic_sinus_potential.dat')
time_1,msd_mean_1,msd_std_1,msd_nc_mean_1,msd_nc_std_1,msd_nc_cc_mean_1,msd_nc_cc_std_1,msd_isf_mean_1,msd_isf_std_1=meanstd_msd(samplenumber_5,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code/Data_observables_periodic_sinus_potential_V0_1/Data_','/data_observables_periodic_sinus_potential.dat')
time_3,msd_mean_3,msd_std_3,msd_nc_mean_3,msd_nc_std_3,msd_nc_cc_mean_3,msd_nc_cc_std_3,msd_isf_mean_3,msd_isf_std_3=meanstd_msd(samplenumber_5,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code/Data_observables_periodic_sinus_potential_V0_3/Data_','/data_observables_periodic_sinus_potential.dat')
time_5,msd_mean_5,msd_std_5,msd_nc_mean_5,msd_nc_std_5,msd_nc_cc_mean_5,msd_nc_cc_std_5,msd_isf_mean_5,msd_isf_std_5=meanstd_msd(samplenumber_5,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code/Data_observables_periodic_sinus_potential_V0_5/Data_','/data_observables_periodic_sinus_potential.dat')
time_10,msd_mean_10,msd_std_10,msd_nc_mean_10,msd_nc_std_10,msd_nc_cc_mean_10,msd_nc_cc_std_10,msd_isf_mean_10,msd_isf_std_10=meanstd_msd(samplenumber_5,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code/Data_observables_periodic_sinus_potential_V0_10/Data_','/data_observables_periodic_sinus_potential.dat')

D=1
#a=1/(2*np.pi)
a=1
tau=a**2/D
d=a


def Randomwalk1D(n): #n here is the no. of steps that we require
   x = 0.3
   y = 0.3
   xposition = [0.3] #starting from origin (0,0)
   yposition = [0.3] 
   dx=0.03
   for i in range (1,n+1):
       step = np.random.uniform(0,1)
       if step < 0.5: # if step is less than 0.5 we move up    
           x += dx
           y += dx  #moving up in u direction
       if step > 0.5: # if step is greater than 0.5 we move down  
           x += dx
           y += -dx #moving down in y direction
 
       xposition.append(x)
       yposition.append(y)
   return [xposition,yposition]
Randwalk = Randomwalk1D(15) #creating an object for the Randomwalk1D class and passing value of n as 100


def cosine_potential(x, U, a):
    return np.cos(x*2*np.pi)


#----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 10)) 

left, bottom, width, height = [0.65, 0.59, 0.28, 0.09]
ax3 = fig.add_axes([left, bottom, width, height])

x = np.linspace(-1, 1, 1000)
ax3.plot(x, cosine_potential(x, 1.0, 1.0),  color='black')
ax3.spines['bottom'].set_color('white')
ax3.spines['top'].set_color('white')
ax3.spines['left'].set_color('white')
ax3.spines['right'].set_color('white')
ax3.set_ylabel(r'$U(x)$')
ax3.set_xticks([], [])
ax3.set_yticks([], [])
dot = np.linspace(0.5, 0.5, 1000)
ax3.plot(0.5, 0.6,marker='o',markersize=12.0,color='b')
ax3.arrow(0.5, 0.6, 0.2, 0.0,head_width=0.2, head_length=0.1, width=0.01, fc='k', ec='k')
ax3.arrow(0.5, 0.6, -0.2, 0.0,head_width=0.2, head_length=0.1, width=0.01, fc='k', ec='k')

ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.set_xlabel(r'$t/\tau$')
ax1.set_ylabel(r'$\langle \Delta x(t)^2\rangle/ a^2$')

# Define custom ticks and tick labels for the x-axis
ax1.set_xticks([10**i for i in range(-5, 10)])
ax1.tick_params(axis='x', which='major', pad=5, length=4)
ax1.tick_params(axis='x', which='minor', length=2)
ax1.xaxis.set_minor_locator(ticker.LogLocator(subs=[2, 3, 4, 5, 6, 7, 8, 9],numticks=10))

# Define custom ticks and tick labels for the y-axis

ax1.set_yticks([10**i for i in range(-5, 5)])
ax1.tick_params(axis='y', which='minor', length=2)
ax1.yaxis.set_minor_locator(ticker.LogLocator(subs=[2, 3, 4, 5, 6, 7, 8, 9],numticks=10))
ax1.tick_params(axis='both', which='both', pad=6)


plot_options = [{"color": "#93A9EB"},
                {"color": "#6887E3"},
                {"color": "#3D65DB"},
                {"color": "blue"},
                {"color": "#193073"}]

for i, (time, msd_mean, msd_std) in enumerate(zip([time_01, time_1, time_3, time_5, time_10],
                                                  [msd_mean_01, msd_mean_1, msd_mean_3, msd_mean_5, msd_mean_10],
                                                  [msd_std_01, msd_std_1, msd_std_3, msd_std_5, msd_std_10])):
    plot_opts = plot_options[i]
    plot_opts.update({"label": f'${{{[0.1, 1, 3, 5, 10][i]}}}$'})
    ax1.plot(time / tau, msd_mean / d ** 2, **plot_opts)
    
# Add legend with custom font size and style
x_free = np.logspace(-3, -1)
x_eff = np.logspace(0, 1)

legend = ax1.legend(title=r"$\Delta U/k_\mathrm{B} T$",handlelength=1.0, handletextpad=0.5, framealpha=1)
plt.setp(legend.get_title(), fontsize=size)
plt.setp(legend.get_texts(), fontsize=size)
ax1.tick_params(axis='both', which='both', right='on',top='on',pad=6)
ax1.set_ylim(6*10**(-4), 5*10**(1))



ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_xlabel(r'$t/\tau$')

# Data to plot
msd_data = [
    (time_1, msd_nc_mean_1, msd_nc_std_1, r"$1.0 k_\mathrm{B} T$", 'solid', 'red', 0,20), #orange
    (time_01, msd_nc_mean_01, msd_nc_std_01, r"$0.1 k_\mathrm{B} T$", '-', '#ff7f7f', 13, 20), #darkbrown
    #(time_10, msd_nc_mean_10, msd_nc_std_10, "10.0", 'dotted', 'blue',20 , 20),
]

# Plot MSD data
for t, y, yerr, label, linestyle, color, position, rotation in msd_data:
    ax2.errorbar(t/ tau, y/ d**2, yerr=yerr/ d**2, linestyle=linestyle, color=color)
    ax2.text(t[position]/tau,2* y[position]/ d**2, r"$\Delta U$="+f"{label}", ha='left', va='bottom', rotation=rotation, color=color, fontsize=size)
 
# Plot complementary MSD data
complementary_data = [
    (time_1, np.abs(msd_nc_cc_mean_1), msd_nc_cc_std_1, r"$1.0k_\mathrm{B} T$",'solid', '#0f9200',0, 19,0.00000001), #pink
    (time_01, np.abs(msd_nc_cc_mean_01), msd_nc_cc_std_01,r"$0.1 k_\mathrm{B} T$", '-', '#30cb00',13, 20,0.0000001), #darkpurple
]

for t, y, yerr, label, linestyle, color, position, rotation, yval in complementary_data:
    ax2.errorbar(t/tau, np.abs(y)/ d**2, yerr=yerr/ d**2, linestyle=linestyle, color=color) #ACHTUNG zum Paper, hier kein faktor 2
    ax2.text(t[position]/tau, np.abs(y[position])/ d**2+yval, r"$\Delta U$="+f"{label}", ha='left', va='bottom', rotation=rotation, color=color, fontsize=size)

# Add custom legend
ax2.set_ylim(10**(-11),10**(1))

p1, = ax2.plot([0], [0], color="#ff7f7f", linestyle='-')
p2, = ax2.plot([0], [0], color='red',  linestyle='-')
p3, = ax2.plot([0], [0], linestyle='-', color='#30cb00')
p4, = ax2.plot([0], [0], linestyle='-', color='#0f9200')
legend2 = ax2.legend([(p3, p4),(p1, p2)], [r'$|\langle  \Delta x^{\mathrm{red}}(t) \Delta x(t)\rangle| / a ^2$  ',r'$\langle \Delta x^{\mathrm{red}}(t)^2\rangle / a^2$'], loc='upper left',handlelength=2.2,  handler_map={tuple: HandlerTuple(ndivide=None)})
ax2.tick_params(axis='both', which='both', right='on',top='on',pad=6)
ax2.add_artist(legend2)
ax1.text(-0.15, 0.95, '(a)', transform=ax1.transAxes, fontsize=20,
         verticalalignment='top')
ax2.text(-0.15, 0.95, '(b)', transform=ax2.transAxes, fontsize=20,
         verticalalignment='top')
plt.tight_layout()
#plt.savefig('../plots/sinus_potential/Fig_7.pdf')

plt.show()


#------------------
import matplotlib.ticker as ticker
from matplotlib.legend_handler import HandlerTuple


import matplotlib.colors as mcolors

fig, ax = plt.subplots()

# Set up plot
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'$t/\tau$')
ax.set_ylabel(r'$\Sigma (t)/a ^2$')



# Set the major and minor ticks on the y-axis
y_major_ticks = ticker.LogLocator(base=10, subs=[1.0], numticks=10)
y_minor_ticks = ticker.LogLocator(base=10, subs=np.arange(2, 10) * 0.1, numticks=10)
ax.yaxis.set_major_locator(y_major_ticks)
ax.yaxis.set_minor_locator(y_minor_ticks)

#make alpha values for error bars
errorbarcolor_darkorange = (*mcolors.hex2color('#fd8d3c') , 0.5) 
errorbarcolor_lightorange = (*mcolors.hex2color('#a63603') , 0.5) 
# Data to plot


msd_data = [
    
    (time_01, msd_std_01, np.zeros(len(time_10)), r"$1$ NS" , 'solid', "<", '#9ecae1','#6BAED6' ), #ultra lightblue
    (time_01, msd_nc_std_01, np.zeros(len(time_10)), r"$10$ NS", ':',">",'#ff7f7f','#ff7f7f' ), #ultra leightorange
    (time_1, msd_std_1, np.zeros(len(time_10)), r"$1$", 'solid',"o", 'blue','blue'), #leightblue
    (time_1, msd_nc_std_1, np.zeros(len(time_10)), r"$1$ NS" , ':', "s", 'red','red' ), #lightorange
    (time_10, msd_std_10, np.zeros(len(time_10)), r"$10$ ", '-', "^",  '#193073','#193073' ), #darkblue
    (time_10, msd_nc_std_10, np.zeros(len(time_10)), r"$10$ NS", ':',"D",'#6b0202','#6b0202' ), #darkorange

]

# Plot MSD data
for t, y, yerr, label, linestyle,marker, color,ecolor in msd_data:
    ax.errorbar(t/ tau, y/ d**2, yerr=(yerr/ d**2,yerr/ d**2) , linewidth=1, ecolor=ecolor,linestyle=linestyle, marker=marker, color=color, label=label)
 
p2, = ax.plot([0], [0], 'o', linestyle='-', color="blue") #leightblue
p1, = ax.plot([0], [0], 's', linestyle=':', color='red') #leightorange
p3, = ax.plot([0], [0], 'D', linestyle=':', color='#6b0202') #darkorange
p4, = ax.plot([0], [0], '^', linestyle='-', color='#193073') #darkblue
p5, = ax.plot([0], [0], '>', linestyle=':', color='#ff7f7f') #ultra leight orange
p6, = ax.plot([0], [0], '<', linestyle='-', color='#6BAED6') #ultra lightblue

legend3 =  ax.legend(handles=[(r'red \medspace std'), (p5,p6), ( p1, p2), ( p3, p4)], labels=[r"$\Delta U/k_\mathrm{B} T$" ,r'$0.1$',r'$1.0$', r'$10.0$'],handlelength=3.2,handler_map={tuple: HandlerTuple(ndivide=None)})
ax.add_artist(legend3)

ax.tick_params(axis='both', which='both', right='on',top='on',pad=6)
plt.tight_layout()
#plt.savefig('../plots/sinus_potential/Fig_8.pdf')
plt.show()


#----------------------------------------------------------------------------------------------------------------------------------------

fig, ax = plt.subplots()
# Set up plot
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'$t/\tau$')
ax.set_ylabel(r'$\Simga (t)/a ^2$')



# Set the major and minor ticks on the y-axis
y_major_ticks = ticker.LogLocator(base=10, subs=[1.0], numticks=10)
y_minor_ticks = ticker.LogLocator(base=10, subs=np.arange(2, 10) * 0.1, numticks=10)
ax.yaxis.set_major_locator(y_major_ticks)
ax.yaxis.set_minor_locator(y_minor_ticks)

#make alpha values for error bars
errorbarcolor_darkorange = (*mcolors.hex2color('#fd8d3c') , 0.5) 
errorbarcolor_lightorange = (*mcolors.hex2color('#a63603') , 0.5) 
# Data to plot


msd_data = [
    (time_01, msd_nc_std_01, np.zeros(len(time_10)), r"$10$ NS", ':',">",'#ff7f7f','#ff7f7f' ), #ultra leightorange
    (time_1, msd_nc_std_1, np.zeros(len(time_10)), r"$1$ NS" , ':', "s", 'red','red' ), #lightorange
    (time_10, msd_nc_std_10, np.zeros(len(time_10)), r"$10$ NS", ':',"D",'#6b0202','#6b0202' ), #darkorange
    (time_01, msd_nc_cc_std_01, np.zeros(len(time_10)), r"$10$ NS", ':',">",'black','black' ), #ultra leightorange
    (time_1, msd_nc_cc_std_1, np.zeros(len(time_10)), r"$1$ NS" , ':', "s", 'black','black' ), #lightorange
    (time_10, msd_nc_cc_std_10, np.zeros(len(time_10)), r"$10$ NS", ':',"D",'black','black' ), #darkorange

]

# Plot MSD data
for t, y, yerr, label, linestyle,marker, color,ecolor in msd_data:
    ax.errorbar(t/ tau, y/ d**2, yerr=(yerr/ d**2,yerr/ d**2) , linewidth=1, ecolor=ecolor,linestyle=linestyle, marker=marker, color=color, label=label)
 
p2, = ax.plot([0], [0], 'o', linestyle='-', color="blue") #leightblue
p1, = ax.plot([0], [0], 's', linestyle=':', color='red') #leightorange
p3, = ax.plot([0], [0], 'D', linestyle=':', color='#6b0202') #darkorange
p4, = ax.plot([0], [0], '^', linestyle='-', color='#193073') #darkblue
p5, = ax.plot([0], [0], '>', linestyle=':', color='#ff7f7f') #ultra leight orange
p6, = ax.plot([0], [0], '<', linestyle='-', color='#6BAED6') #ultra lightblue

legend3 =  ax.legend(handles=[(r'red \medspace std'), (p5,p6), ( p1, p2), ( p3, p4)], labels=[r"$\Delta U/k_\mathrm{B} T$" ,r'$0.1$',r'$1.0$', r'$10.0$'],handlelength=3.2,handler_map={tuple: HandlerTuple(ndivide=None)})
ax.add_artist(legend3)

ax.tick_params(axis='both', which='both', right='on',top='on',pad=6)
plt.tight_layout()
plt.savefig('../plots/sinus_potential/Fig_neu.pdf')
plt.show()


#2piperiodisch
samplenumber2=100
time_01,msd_mean_01,msd_std_01,msd_nc_mean_01,msd_nc_std_01,msd_nc_cc_mean_01,msd_nc_cc_std_01,msd_isf_mean_01,msd_isf_std_01=meanstd_msd(samplenumber2,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code_small_times/Data_observables_periodic_sinus_potential_V0_01_small/Data_','/data_observables_periodic_sinus_potential.dat')
time_1,msd_mean_1,msd_std_1,msd_nc_mean_1,msd_nc_std_1,msd_nc_cc_mean_1,msd_nc_cc_std_1,msd_isf_mean_1,msd_isf_std_1=meanstd_msd(samplenumber2,'../Data/Data_observables_periodic_sinus_potential_V0_modified_code_small_times/Data_observables_periodic_sinus_potential_V0_1_small/Data_','/data_observables_periodic_sinus_potential.dat')


D=1
a=1
#a=2*np.pi
tau=a**2/D
d=a


dimension=1

def numerisch(x, y, y_fehler):
    f4_der=[]
    x_f4=[]
    y_vacf = unp.uarray(y,y_fehler)
    for i in range(1,len(x)-1):
        if i%9!=0:
            h1 = x[i]-x[i-1]
            h2 = x[i+1]-x[i]
            f4_der.append((h1*y_vacf[i+1]-(h1+h2)*y_vacf[i]+h2*y_vacf[i-1])/(0.5*(h2*h1**2+h1*h2**2)))
            x_f4.append(x[i])

    x_f4=np.array(x_f4)
    f4_der=np.array(f4_der)

    return x_f4[0:27], unp.nominal_values(f4_der[0:27])/2./dimension, unp.std_devs(f4_der[0:27])/2./dimension




# Plot MSD data
fig, ax = plt.subplots()

# Data to plot
msd_data = [
    (time_1, -msd_mean_01,msd_std_01 , r'STD', 'solid', "o", 'blue', 'blue'), #blue
    (time_1, msd_nc_mean_01, msd_nc_std_01, r'NC', 'solid', "^", 'red','red' ), #orange
]


for t, y, yerr, label, linestyle, marker, color, ecolor in msd_data:
    x, y, y_err = numerisch(t/tau, y/a**2, yerr/a**2)
    ax.errorbar(x, y, yerr=y_err, linewidth=1, ecolor=ecolor, linestyle=linestyle, marker=marker, color=color, label=label)

# Set axes labels, limits, and scales
ax.set_xlabel(r'$t/\tau$')
ax.set_ylabel(r'$-Z(t)a^2/D^2$')
#ax.set_xlim([2*10**(0), 2*10**(2)])
ax.set_xscale('log')
ax.set_yscale('log')

# Add legend and grid
ax.legend(loc='best')
plt.tight_layout()
ax.tick_params(axis='both', which='both', right='on',top='on',pad=6)
#plt.savefig('../plots/sinus_potential/Fig_9.pdf')
# Show plot
plt.show()





# Plot MSD data
fig, ax = plt.subplots()

# Data to plot
msd_data = [
    (time_1, -msd_nc_cc_mean_1,msd_nc_cc_std_1 , r'Frenkel', 'solid', "o", 'blue', 'blue'), #blue
    (time_1, msd_nc_mean_1, msd_nc_std_1, r'NC', 'solid', "^", 'red','red' ), #orange
]


for t, y, yerr, label, linestyle, marker, color, ecolor in msd_data:
    x, y, y_err = numerisch(t/tau, y/a**2, yerr/a**2)
    ax.errorbar(x, y, yerr=y_err, linewidth=1, ecolor=ecolor, linestyle=linestyle, marker=marker, color=color, label=label)

# Set axes labels, limits, and scales
ax.set_xlabel(r'$t/\tau$')
ax.set_ylabel(r'$-Z(t)a^2/D^2$')
#ax.set_xlim([2*10**(0), 2*10**(2)])
ax.set_xscale('log')
ax.set_yscale('log')

# Add legend and grid
ax.legend(loc='best')
plt.tight_layout()
ax.tick_params(axis='both', which='both', right='on',top='on',pad=6)
#plt.savefig('../plots/sinus_potential/Fig_neu_vacf.pdf')
# Show plot
plt.show()









#-------------------------------------------------------------------------------------------------------------------------------------

#compute the sem of only r samples, because if i use semsem as my error of the sem, I need a new "mean"
def sem_subset(data,n):
    r=1
    l=int(n/r)
    # Divide data into r groups of l elements each
    data_split = np.split(data, r) #choose here for example 20=4*5 so 4 sets with 5 elements
    
    # Initialize an array to store the SEMs
    sems = np.empty([r,len(data[0])])
    
    # Loop over the groups and calculate the SEM for each
    for i in range(r):
        group = data_split[i]
        sems[i] = np.std(group, ddof=1,axis=0) / np.sqrt(l)  # we use 5 here because each group has 5 elements
    means=np.mean(sems,axis=0)
    semsem_subset = np.std(sems, ddof=1,axis=0) / np.sqrt(r) 
    return means, sems, semsem_subset
    

def meanstd_msd_subset(num,pfad,pfad2):
    #ddof=1 to get unbias std means with 1/n-1 *sum over sample (NOT the whole smple but just random subset)
    start=0

    y=np.array([data_msd(pfad+str(i)+pfad2) for i in range(start,num)])

    time_here_mean=np.mean( np.array([y[i][3] for i in range(0,num-start)]), axis=0 )

    arr=np.array([y[i][4] for i in range(0,num-start)])
    msd_here_std=sem_subset(arr,num)[0]
    msd_here_std_std=sem_subset(arr,num)[2]
            
    arr_nc= np.array([y[i][5] for i in range(0,num-start)])
    msd_here_nc_std=sem_subset(arr_nc,num)[0]
    msd_here_nc_std_std=sem_subset(arr_nc,num)[2]

    return time_here_mean,msd_here_std,sem_subset(arr,num)[1] ,msd_here_std_std,msd_here_nc_std,sem_subset(arr_nc,num)[1],msd_here_nc_std_std