data133bis.py

# 20.03.2025

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.special import erf
from collections.abc import Iterable
import sys, os, time
from PIL import Image

save_dic = {}

plt.rcParams.update({
    "font.size": 20,              # Taille globale du texte
    "axes.labelsize": 25,         # Labels des axes
    "axes.titlesize": 25,         # Titre
    "xtick.labelsize": 20,        # Graduation axe x
    "ytick.labelsize": 20,        # Graduation axe y
    "legend.fontsize": 20,        # Légende
    "lines.linewidth": 4,       # Épaisseur des lignes
    "figure.figsize": (10, 7)     # Taille de la figure (en pouces)
})


### FUNCTIONS ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
### ___________________________________________________________________________________________________________________________________________________________


### tools

# error print on/off
def error_print_off():
    save_dic["original_stderr"] = sys.stderr
    class NullWriter:
        def write(self, *args, **kwargs): pass
    sys.stderr = NullWriter()
def error_print_on(): 
    try: sys.stderr = save_dic["original_stderr"]
    except: pass

# array check
def is_array(A):
    return isinstance(A, (list, tuple, np.ndarray)) or isinstance(A, Iterable)

# gif maker
def make_gif(output, img_del=True):
    image_filenames = []
    images = []
    for filename in sorted(os.listdir(out_folder)):
        if filename.endswith('.png'):
            print(filename)
            image_filenames.append(filename)
            images.append(Image.open(os.path.join(out_folder, filename)))
    print(" > bip-bouping in progress...")
    images[0].save(output, save_all=True, append_images=images[1:], duration=100, loop=0)
    print(f" > Done : {output}")
    if img_del:
        for filename in image_filenames: os.remove(os.path.join(out_folder, filename))

def pntply(para, var="c"):
    if len(para)==1: return f"{para[0]:.03}"
    if len(para)==2: return f"{para[0]:.03}c{para[1]:+.03}"
    c = ""
    for i,p in enumerate(para):
        if   i==0:           c += f"{p:.03}c{len(para)-1}"
        elif i==len(para)-1: c += f"{p:+.03}"
        elif i==len(para)-2: c += f"{p:+.03}c"
        else:                c += f"{p:+.03}c{len(para)-1}"
    return c

# read functions
def read(path):
    data = {"t":[],"z":[],"c":[],"m":[]}
    for line in list(open(path))[3:]:
        if line[:2] != "  ":
            data["t"].append(float(line.split(' ')[0]))
            for k in data: 
                if k!="t": data[k].append([])
        else:
            _,z,c,m = [float(x) for x in line.strip("\n").strip("  ").split(' ')]
            data["z"][-1].append(z); data["c"][-1].append(c); data["m"][-1].append(m)
    for k in data: data[k] = np.array(data[k]) # list -> nd.array
    return data
def read2(path1, path2):
    '''-> T, Z, Cs_solv, Cs_poly'''
    data = {"solvant":read(path1), "polymer":read(path2)}
    T = data["solvant"]["t"]    # same t for each types
    Z = data["solvant"]["z"][0] # same z for each t & types
    Cs_solv = data["solvant"]["c"]
    Cs_poly = data["polymer"]["c"]
    print(" > bip boup")
    return T, Z, Cs_solv, Cs_poly

# find zero
def zero(X,Y,i0=0):
    '''-> x0 / y(x0) = 0'''
    if Y[i0] == 0: return X[i0]
    sgn0 = Y[i0]/abs(Y[i0])
    for i,y in zip(range(i0,len(Y)),Y[i0:]):
        if Y[i] == 0:              return X[i]
        if Y[i]/abs(Y[i]) != sgn0: return X[i-1] + (X[i]-X[i-1])*Y[i-1]/(Y[i-1]-Y[i])

# find index
def get_index(X,x0s,i0=0):
    '''-> i / X[i] = x0'''
    tup = is_array(x0s)
    x0 = x0s[0] if tup else x0s
    I = None
    for i in range(i0,len(X)-1):
        if X[i]<=x0<X[i+1] or X[i]>=x0>X[i+1]: I = i; break
    if   tup and len(x0s) == 1: return [I]
    elif tup: return [I] + get_index(X, x0s[1:], I)
    else: return I

# curve fit models
def model_polynomial(x, *args): return sum([a*x**n for a,n in enumerate(args)])
def model_erf(xi, c0,D): return c0/2*(1+erf(xi/(2*D**0.5)))
def model_exp(x, A,k,x0): return A + np.exp(k*(x-x0))

def model_Bochard(c, D0, chi):
    return c*(1-c)*D0*(1/(1-c)-2*chi)

### Calculations

def get_c0(Z,Cs):
    return np.max(Cs)

def fix0_intersection(Z, Cs_solv, Cs_poly, i0=15):
    return np.array([Z-zero(Z,C,i0) for C in Cs_solv-Cs_poly]) # Zfixs

def fix0_halfc0(Z, Cs_solv, Cs_poly, c0_solv, c0_poly, i0=0):
    Zfixs_solv = np.array([Z-zero(Z,C,15) for C in Cs_solv-c0_solv/2])
    Zfixs_poly = np.array([Z-zero(Z,C,15) for C in Cs_poly-c0_poly/2])
    return Zfixs_solv, Zfixs_poly

def ZtoXi(Zs,T):
    if is_array(Zs[0]): return np.array([Z/t**.5 for Z,t in zip(Zs,T)])
    else: return np.array([Zs/t**.5 for t in T])

def get_mean(Xis_solv,Xis_poly,Cs_solv,Cs_poly, mti0,mti1,zi0,zi1, n=100):
    ximin = max(np.max(Xis_solv[mti0:mti1,zi0]), np.max(Xis_poly[mti0:mti1,zi0]))
    ximax = min(np.min(Xis_solv[mti0:mti1,zi1]), np.min(Xis_poly[mti0:mti1,zi1]))
    mXi = np.linspace(ximin,ximax,n)
    mC_solv = np.mean(np.array( [np.interp(mXi,Xi,C) for Xi,C in zip(Xis_solv[mti0:mti1],Cs_solv[mti0:mti1])] ),axis=0)
    mC_poly = np.mean(np.array( [np.interp(mXi,Xi,C) for Xi,C in zip(Xis_poly[mti0:mti1],Cs_poly[mti0:mti1])] ),axis=0)
    return mXi, mC_solv, mC_poly

def regress_Dcst(Xi,C,c0,sign,i0,i1):
    def model(xi,D): return model_erf(sign*xi,c0,D)
    (D,), _ = curve_fit(model,Xi[i0:i1+1],C[i0:i1+1])
    return D

def derive_D(Xi,C,polynomial_fix=0):
    '''integral_fix = "", "linear", "polynomial"'''
    Xip = np.gradient(Xi,C)
    Xii = np.zeros(len(Xi))
    for i in range(len(Xii)-1): Xii[i+1] = Xii[i] + (C[i+1]-C[i])*(Xi[i]+Xi[i+1])/2
    if   polynomial_fix == 1: Xii -= np.polyfit(C[:2],Xii[:2],1)[-1]
    elif polynomial_fix != 0: Xii -= np.polyfit(C[:int(len(C)*0.2)],Xii[:int(len(Xii)*0.2)],min(int(len(C)*0.2),polynomial_fix))[-1]
    D = -0.5*Xip*Xii
    return D, Xip, Xii


def plotWhite(fig,ax):
    fig.patch.set_facecolor("#355A6C")
    ax.set_facecolor("#628B9E")
    ax.grid(True, which='both', linestyle='--', linewidth=0.5, color='white', alpha=0.3)
    for spine in ax.spines.values():
        spine.set_color('white')
    ax.tick_params(colors='white', direction='in', length=6, width=1)
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')


### MAINS ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
### _______________________________________________________________________________________________________________________________________________________


path1 = "data133/dens_solvant_1_33__1_33.profile"
path2 = "data133/dens_polymer_1_33__1_33.profile"
paths = (path1,path2)
out_folder = "data_out"



def truc_0():
    T, Z, Cs_solv, Cs_poly = read2(*paths)

    ti0 = 0
    ti1 = len(T)-1
    zi0 = 0
    zi1 = 64

    c0_solv,  c0_poly  = get_c0(Z[ti0:ti1], Cs_solv[ti0:ti1]), get_c0(Z[ti0:ti1], Cs_poly[ti0:ti1])
    fig,ax = plt.subplots()
    for ti in range(ti0,ti1):
        c = (ti-ti0)/(ti1-ti0)
        ax.plot(Z[zi0:zi1], Cs_solv[ti,zi0:zi1], lw=3, color=(1-c,1,0,.5), zorder=10)
        ax.plot(Z[zi0:zi1], Cs_poly[ti,zi0:zi1], lw=3, color=(0,1-c,1,.5), zorder=10)
        ax.plot(Z[zi0:zi1], Cs_solv[ti,zi0:zi1]+Cs_poly[ti,zi0:zi1], lw=3, color=(.7,0,1-c,.5), zorder=10)
    ax.hlines((c0_solv,c0_poly), np.min(Z), np.max(Z), lw=1.5, ls=':', color="crimson", alpha=0.5, zorder=12)
    ax.grid()
    ax.set_xlabel("z")
    ax.set_ylabel("c")
    plt.show()
    plt.close() 


def truc_0bis():
    T, Z, Cs_solv, Cs_poly = read2(*paths)
    Cs_tot = Cs_solv + Cs_poly

    ti0 = 0
    ti1 = 199
    zi0 = 0
    zi1 = 64

    c0_solv,  c0_poly  = get_c0(Z[ti0:ti1], Cs_solv[ti0:ti1]), get_c0(Z[ti0:ti1], Cs_poly[ti0:ti1])
    fig,ax = plt.subplots()
    for ti in range(ti0,ti1):
        c = (ti-ti0)/(ti1-ti0)
        ax.plot(Z[zi0:zi1], Cs_solv[ti,zi0:zi1], lw=3, color=(1-c,1-0.5*c,0,.5), alpha=0.05, zorder=10)
        ax.plot(Z[zi0:zi1], Cs_poly[ti,zi0:zi1], lw=3, color=(0,1-c,1,.5),       alpha=0.05, zorder=10)
        if ti%1==0: 
            ax.plot(Z[zi0:zi1], Cs_tot[ti,zi0:zi1],  lw=3, color=(1-c,1-c,1-c,.5),   alpha=1, zorder=11)
        if ti%10==0:
            ax.plot(Z[zi0:zi1], Cs_tot[ti+1,zi0:zi1],  lw=4, color="black",   alpha=0.8, zorder=11)
    #ax.hlines((c0_solv,c0_poly), np.min(Z), np.max(Z), lw=1.5, ls=':', color="crimson", alpha=0.5, zorder=12)
    ax.set_xlabel("z")
    ax.set_ylabel("c")
    plotWhite(fig,ax)

    plt.show()
    plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# _____________________________________________________________________________________________________________________________


def truc_1():
    T, Z, Cs_solv, Cs_poly = read2(*paths)

    # time indices
    ti0 = 100
    ti1 = 200

    mti0 = 179
    mti1 = 199

    zi0 = 0
    zi1 = 64

    # regression xi intervals
    xi0_solv =  0.001
    xi1_solv =  0.006
    xi0_poly = -0.003
    xi1_poly =  0.001

    # c0 & ordinate fix
    c0_solv,  c0_poly  = get_c0(Z[mti0:mti1], Cs_solv[mti0:mti1]), get_c0(Z[mti0:mti1], Cs_poly[mti0:mti1])
    Zs_solv,  Zs_poly  = fix0_halfc0(Z, Cs_solv, Cs_poly, c0_solv, c0_poly, i0=15)
    Xis_solv, Xis_poly = ZtoXi(Zs_solv,T), ZtoXi(Zs_poly,T)

    # C(Xi) mean
    mXi, mC_solv, mC_poly = get_mean(Xis_solv,Xis_poly,Cs_solv,Cs_poly, mti0,mti1,zi0,zi1, zi1-zi0)

    do_regress = False

    # constant D regression
    if do_regress:
        ri0_solv, ri1_solv = get_index(mXi,(xi0_solv,xi1_solv))
        ri0_poly, ri1_poly = get_index(mXi,(xi0_poly,xi1_poly))
        D0_solv = regress_Dcst(mXi,mC_solv,c0_solv,  1, ri0_solv,ri1_solv)
        D0_poly = regress_Dcst(mXi,mC_poly,c0_poly, -1, ri0_poly,ri1_poly)
        rC_solv, rC_poly = model_erf(mXi,c0_solv,D0_solv), model_erf(-mXi,c0_poly,D0_poly)

    plot_C = False

    # plot C
    if plot_C:
        fig,ax = plt.subplots()
        for ti in range(ti0,ti1):
            c = (ti-ti0)/(ti1-ti0)
            ax.plot(Xis_solv[ti,zi0:zi1], Cs_solv[ti,zi0:zi1], lw=3, color=(1-c,1,0,.5), zorder=10)
            ax.plot(Xis_poly[ti,zi0:zi1], Cs_poly[ti,zi0:zi1], lw=3, color=(0,1-c,1,.5), zorder=10)
        ax.hlines((c0_solv,c0_poly), np.min(mXi), np.max(mXi), lw=1.5, ls=':', color="crimson", alpha=0.5, zorder=12)
        ax.plot(mXi, mC_solv, ls=':', lw=1.5, color="black", zorder=11)
        ax.plot(mXi, mC_poly, ls=':', lw=1.5, color="white", zorder=11)
        if do_regress: 
            ax.scatter((mXi[ri0_solv],mXi[ri1_solv],mXi[ri0_poly],mXi[ri1_poly]),
                       (mC_solv[ri0_solv],mC_solv[ri1_solv],mC_poly[ri0_poly],mC_poly[ri1_poly]), color="crimson", alpha=0.7, zorder=14)
            print(f"Regression results:\nSolvant: D ~ {D0_solv}\nPolymer: D ~ {D0_poly}")
            ax.plot(mXi, rC_solv, lw=2, color="crimson", alpha=0.7, zorder=13)
            ax.plot(mXi, rC_poly, lw=2, color="crimson", alpha=0.7, zorder=13)
        ax.grid()
        ax.set_xlabel("ξ = z/√t")
        ax.set_ylabel("c")
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~


    # new Xi lim
    xi0_solv = -0.008
    xi1_solv =  0.008
    xi0_poly = -0.007
    xi1_poly =  0.012

    Xi_solv = np.linspace(xi0_solv,xi1_solv,len(mXi))
    Xi_poly = np.linspace(xi1_poly,xi0_poly,len(mXi)) # !!
    mC_solv = np.interp(Xi_solv, mXi, mC_solv)
    mC_poly = np.interp(Xi_poly, mXi, mC_poly)
    if do_regress:
        rC_solv = np.interp(Xi_solv, mXi, rC_solv)
        rC_poly = np.interp(Xi_poly, mXi, rC_poly)

    poly_fix = 4

    # D derivation
    error_print_off()
    mD_solv, mXip_solv, mXii_solv = derive_D(Xi_solv, mC_solv, poly_fix)
    mD_poly, mXip_poly, mXii_poly = derive_D(Xi_poly, mC_poly, poly_fix)
    if do_regress:
        rD_solv, rXip_solv, rXii_solv = derive_D(Xi_solv, rC_solv, poly_fix)
        rD_poly, rXip_poly, rXii_poly = derive_D(Xi_poly, rC_poly, poly_fix)
    error_print_on()

    plot_D_exp = False
    plot_D_erf = False

    # plot D exp derivation
    if plot_D_exp:
        fig,((axs0,axs1,axs2,axs3),(axp0,axp1,axp2,axp3)) = plt.subplots(2,4)
        axs0.plot(mC_solv, Xi_solv,   lw=2, color="darkgreen",       label="ξ(c)")
        axs1.plot(mC_solv, mXip_solv, lw=2, color="yellowgreen",     label="dξ/dc")
        axs2.plot(mC_solv, mXii_solv, lw=2, color="orange",          label="∫dcξ")
        axs3.plot(mC_solv, mD_solv,   lw=2, color="crimson",         label="D(c)")
        axp0.plot(mC_poly, Xi_poly,   lw=2, color="darkblue",        label="ξ(c)")
        axp1.plot(mC_poly, mXip_poly, lw=2, color="rebeccapurple",   label="dξ/dc")
        axp2.plot(mC_poly, mXii_poly, lw=2, color="mediumvioletred", label="∫dcξ")
        axp3.plot(mC_poly, mD_poly,   lw=2, color="crimson",         label="D(c)")
        for ax in (axs0,axs1,axs2,axs3,axp0,axp1,axp2,axp3):
            ax.grid()
            ax.legend()
            ax.set_xlabel("c")
        axs3.set_ylim(0,None); axp3.set_ylim(0,None)
        fig.suptitle("Analysis from experimental data")

    # plot D erf derivation
    if plot_D_erf and do_regress:
        fig,((axs0,axs1,axs2,axs3),(axp0,axp1,axp2,axp3)) = plt.subplots(2,4)
        axs0.plot(rC_solv, Xi_solv,   lw=2, color="darkgreen",       label="ξ(c)")
        axs1.plot(rC_solv, rXip_solv, lw=2, color="yellowgreen",     label="dξ/dc")
        axs2.plot(rC_solv, rXii_solv, lw=2, color="orange",          label="∫dcξ")
        axs3.plot(rC_solv, rD_solv,   lw=2, color="crimson",         label="D(c)")
        axp0.plot(rC_poly, Xi_poly,   lw=2, color="darkblue",        label="ξ(c)")
        axp1.plot(rC_poly, rXip_poly, lw=2, color="rebeccapurple",   label="dξ/dc")
        axp2.plot(rC_poly, rXii_poly, lw=2, color="mediumvioletred", label="∫dcξ")
        axp3.plot(rC_poly, rD_poly,   lw=2, color="crimson",         label="D(c)")
        axs3.hlines((D0_solv),min(rC_solv),max(rC_solv), ls=':', color="darkred", alpha=0.6)
        axp3.hlines((D0_poly),min(rC_poly),max(rC_poly), ls=':', color="darkred", alpha=0.6)
        for ax in (axs0,axs1,axs2,axs3,axp0,axp1,axp2,axp3):
            ax.grid()
            ax.legend()
            ax.set_xlabel("c")
        axs3.set_ylim(0,None); axp3.set_ylim(0,None)
        axs2.set_title(f"polynomial fix at order {poly_fix}", fontsize=9)
        fig.suptitle("Analysis from the erf regression")

    if plot_D_erf or plot_D_exp:
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
    
    # fit D indices
    fit_i0_solv = 0
    fit_i1_solv = get_index(mC_solv, 4000)
    fit_i2_solv = -1
    fit_i0_poly = get_index(mC_poly,  200)
    fit_i1_poly = get_index(mC_poly, 5000)

    # fit X,Y
    fit_para1_solv = np.polyfit(mC_solv[fit_i0_solv:fit_i1_solv], mD_solv[fit_i0_solv:fit_i1_solv], 1)
    fit_para2_solv = np.polyfit(mC_solv[fit_i1_solv:fit_i2_solv], mD_solv[fit_i1_solv:fit_i2_solv], 1)
    fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly], np.log(mD_poly[fit_i0_poly:fit_i1_poly]), 2)
    fitX1_solv = np.linspace(0,                    mC_solv[fit_i1_solv], 100)
    fitX2_solv = np.linspace(mC_solv[fit_i1_solv], 5500,                 100)
    fitX1_poly = np.linspace(0,                    8000,                 100)
    fitY1_solv = np.polyval(fit_para1_solv, fitX1_solv)
    fitY2_solv = np.polyval(fit_para2_solv, fitX2_solv)
    fitY1_poly = np.exp(np.polyval(fit_para1_poly, fitX1_poly))

    plot_D = True

    # plot D
    if plot_D:
        fig,(ax0,ax1) = plt.subplots(1,2,figsize=(17,9))
        ax0.plot(mC_solv, mD_solv, lw=4, color=(0,0.5,0), zorder=8)
        ax0.plot(mC_solv, mD_solv, lw=8, color="white",   zorder=7)
        #ax0.scatter(*[[X[i] for i in (fit_i0_solv,fit_i1_solv,fit_i2_solv)] for X in (mC_solv,mD_solv)], s=100, color="darkgreen")
        ax1.plot(mC_poly, mD_poly, lw=4, color=(0,0,1), zorder=8)
        ax1.plot(mC_poly, mD_poly, lw=8, color="white", zorder=7)
        #ax1.scatter(*[[X[i] for i in (fit_i0_poly,fit_i1_poly)] for X in (mC_poly,mD_poly)], s=100, color="darkblue")
        for ax in (ax0,ax1):
            ax.grid()
            ax.set_xlabel("c")
            ax.set_ylabel("D")
            plotWhite(fig,ax)
        plt.subplots_adjust(left=0.05, right=0.95, wspace=0.2)
        #ax0.set_title("solvant")
        #ax1.set_title("polymer")
        ax0.set_xlim(0,5500)
        ax0.set_ylim(0,None)
        ax1.set_xlim(0,8000)
        ax1.set_ylim(0,None)
        #ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = {pntply(fit_para1_solv)}",      alpha=0.6)
        #ax0.plot(fitX2_solv, fitY2_solv, lw=3, color="orangered", label=f"D = {pntply(fit_para2_solv)}",      alpha=0.6)
        #ax1.plot(fitX1_poly, fitY1_poly, lw=3, color="red",       label=f"D = exp({pntply(fit_para1_poly)})", alpha=0.6)
        #ax0.legend(); ax1.legend()
        plt.show()
        plt.close()


def truc_1():
    T, Z, Cs_solv, Cs_poly = read2(*paths)

    # time indices
    ti0 = 100
    ti1 = 200

    mti0 = 179
    mti1 = 199

    zi0 = 0
    zi1 = 64

    # regression xi intervals
    xi0_solv = -0.008
    xi1_solv =  0.000
    xi0_poly = -0.003
    xi1_poly =  0.001

    # c0 & ordinate fix
    c0_solv,  c0_poly  = get_c0(Z[mti0:mti1], Cs_solv[mti0:mti1]), get_c0(Z[mti0:mti1], Cs_poly[mti0:mti1])
    Zs_solv,  Zs_poly  = fix0_halfc0(Z, Cs_solv, Cs_poly, c0_solv, c0_poly, i0=15)
    Xis_solv, Xis_poly = ZtoXi(Zs_solv,T), ZtoXi(Zs_poly,T)

    # C(Xi) mean
    mXi, mC_solv, mC_poly = get_mean(Xis_solv,Xis_poly,Cs_solv,Cs_poly, mti0,mti1,zi0,zi1, zi1-zi0)

    do_regress = False

    # constant D regression
    if do_regress:
        ri0_solv, ri1_solv = get_index(mXi,(xi0_solv,xi1_solv))
        ri0_poly, ri1_poly = get_index(mXi,(xi0_poly,xi1_poly))
        D0_solv = regress_Dcst(mXi,mC_solv,c0_solv,  1, ri0_solv,ri1_solv)
        D0_poly = regress_Dcst(mXi,mC_poly,c0_poly, -1, ri0_poly,ri1_poly)
        rC_solv, rC_poly = model_erf(mXi,c0_solv,D0_solv), model_erf(-mXi,c0_poly,D0_poly)

    plot_C = False

    # plot C
    if plot_C:
        fig,ax = plt.subplots()
        for ti in range(ti0,ti1):
            c = (ti-ti0)/(ti1-ti0)
            ax.plot(Xis_solv[ti,zi0:zi1], Cs_solv[ti,zi0:zi1], lw=3, color=(1-c,1-0.5*c,0,.5), zorder=10)
            ax.plot(Xis_poly[ti,zi0:zi1], Cs_poly[ti,zi0:zi1], lw=3, color=(0,1-c,1,.5), zorder=10)
        #ax.hlines((c0_solv,c0_poly), np.min(mXi), np.max(mXi), lw=1.5, ls=':', color="crimson", alpha=0.5, zorder=12)
        ax.plot(mXi, mC_solv, ls='--', lw=1.7, color="black", zorder=11)
        ax.plot(mXi, mC_poly, ls='--', lw=1.7, color="white", zorder=11)
        if do_regress: 
            #ax.scatter((mXi[ri0_solv],mXi[ri1_solv],mXi[ri0_poly],mXi[ri1_poly]),
            #           (mC_solv[ri0_solv],mC_solv[ri1_solv],mC_poly[ri0_poly],mC_poly[ri1_poly]), color="crimson", alpha=0.7, zorder=14)
            print(f"Regression results:\nSolvant: D ~ {D0_solv}\nPolymer: D ~ {D0_poly}")
            ax.plot(mXi, rC_solv, lw=7, color="crimson", alpha=0.7, zorder=13)
            ax.plot(mXi, rC_poly, lw=7, color="crimson", alpha=0.7, zorder=13)
        ax.grid()
        ax.set_xlabel("η = z/√t")
        ax.set_ylabel("c")
        plotWhite(fig,ax)
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~


    # new Xi lim
    xi0_solv = -0.008
    xi1_solv =  0.008
    xi0_poly = -0.007
    xi1_poly =  0.012

    Xi_solv = np.linspace(xi0_solv,xi1_solv,len(mXi))
    Xi_poly = np.linspace(xi1_poly,xi0_poly,len(mXi)) # !!
    mC_solv = np.interp(Xi_solv, mXi, mC_solv)
    mC_poly = np.interp(Xi_poly, mXi, mC_poly)
    if do_regress:
        rC_solv = np.interp(Xi_solv, mXi, rC_solv)
        rC_poly = np.interp(Xi_poly, mXi, rC_poly)

    poly_fix = 4

    # D derivation
    error_print_off()
    mD_solv, mXip_solv, mXii_solv = derive_D(Xi_solv, mC_solv, poly_fix)
    mD_poly, mXip_poly, mXii_poly = derive_D(Xi_poly, mC_poly, poly_fix)
    if do_regress:
        rD_solv, rXip_solv, rXii_solv = derive_D(Xi_solv, rC_solv, poly_fix)
        rD_poly, rXip_poly, rXii_poly = derive_D(Xi_poly, rC_poly, poly_fix)
    error_print_on()

    plot_D_exp = False
    plot_D_erf = False

    # plot D exp derivation
    if plot_D_exp:
        fig,((axs0,axs1,axs2,axs3),(axp0,axp1,axp2,axp3)) = plt.subplots(2,4)
        axs0.plot(mC_solv, Xi_solv,   lw=2, color="darkgreen",       label="ξ(c)")
        axs1.plot(mC_solv, mXip_solv, lw=2, color="yellowgreen",     label="dξ/dc")
        axs2.plot(mC_solv, mXii_solv, lw=2, color="orange",          label="∫dcξ")
        axs3.plot(mC_solv, mD_solv,   lw=2, color="crimson",         label="D(c)")
        axp0.plot(mC_poly, Xi_poly,   lw=2, color="darkblue",        label="ξ(c)")
        axp1.plot(mC_poly, mXip_poly, lw=2, color="rebeccapurple",   label="dξ/dc")
        axp2.plot(mC_poly, mXii_poly, lw=2, color="mediumvioletred", label="∫dcξ")
        axp3.plot(mC_poly, mD_poly,   lw=2, color="crimson",         label="D(c)")
        for ax in (axs0,axs1,axs2,axs3,axp0,axp1,axp2,axp3):
            ax.grid()
            ax.legend()
            ax.set_xlabel("c")
        axs3.set_ylim(0,None); axp3.set_ylim(0,None)
        fig.suptitle("Analysis from experimental data")

    # plot D erf derivation
    if plot_D_erf and do_regress:
        fig,((axs0,axs1,axs2,axs3),(axp0,axp1,axp2,axp3)) = plt.subplots(2,4)
        axs0.plot(rC_solv, Xi_solv,   lw=2, color="darkgreen",       label="ξ(c)")
        axs1.plot(rC_solv, rXip_solv, lw=2, color="yellowgreen",     label="dξ/dc")
        axs2.plot(rC_solv, rXii_solv, lw=2, color="orange",          label="∫dcξ")
        axs3.plot(rC_solv, rD_solv,   lw=2, color="crimson",         label="D(c)")
        axp0.plot(rC_poly, Xi_poly,   lw=2, color="darkblue",        label="ξ(c)")
        axp1.plot(rC_poly, rXip_poly, lw=2, color="rebeccapurple",   label="dξ/dc")
        axp2.plot(rC_poly, rXii_poly, lw=2, color="mediumvioletred", label="∫dcξ")
        axp3.plot(rC_poly, rD_poly,   lw=2, color="crimson",         label="D(c)")
        axs3.hlines((D0_solv),min(rC_solv),max(rC_solv), ls=':', color="darkred", alpha=0.6)
        axp3.hlines((D0_poly),min(rC_poly),max(rC_poly), ls=':', color="darkred", alpha=0.6)
        for ax in (axs0,axs1,axs2,axs3,axp0,axp1,axp2,axp3):
            ax.grid()
            ax.legend()
            ax.set_xlabel("c")
        axs3.set_ylim(0,None); axp3.set_ylim(0,None)
        axs2.set_title(f"polynomial fix at order {poly_fix}", fontsize=9)
        fig.suptitle("Analysis from the erf regression")

    if plot_D_erf or plot_D_exp:
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

    # fit D indices
    fit_i0_solv = 0
    fit_i1_solv = get_index(mC_solv, 4000)
    fit_i2_solv = -1
    fit_i0_poly = get_index(mC_poly,  200)
    fit_i1_poly = get_index(mC_poly, 5000)

    # fit X,Y
    fit_para1_solv = np.polyfit(mC_solv[fit_i0_solv:fit_i1_solv], mD_solv[fit_i0_solv:fit_i1_solv], 1)
    fit_para2_solv = np.polyfit(mC_solv[fit_i1_solv:fit_i2_solv], mD_solv[fit_i1_solv:fit_i2_solv], 1)
    fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly], np.log(mD_poly[fit_i0_poly:fit_i1_poly]), 2)
    fitX1_solv = np.linspace(0,                    mC_solv[fit_i1_solv], 100)
    fitX2_solv = np.linspace(mC_solv[fit_i1_solv], 5500,                 100)
    fitX1_poly = np.linspace(0,                    8000,                 100)
    fitY1_solv = np.polyval(fit_para1_solv, fitX1_solv)
    fitY2_solv = np.polyval(fit_para2_solv, fitX2_solv)
    fitY1_poly = np.exp(np.polyval(fit_para1_poly, fitX1_poly))
    
    fit_i0_solv, fit_i1_solv = 12, 28 # 12, 50
    fit_i0_poly, fit_i1_poly = get_index(mC_poly, (1000, 4000))

    nmC_solv, nmC_poly = mC_solv/c0_solv, mC_poly/c0_poly

    #(D0_solv, chi_solv), _ = curve_fit(model_Bochard, nmC_solv[fit_i0_solv:fit_i1_solv], mD_solv[fit_i0_solv:fit_i1_solv])
    #(D0_poly, chi_poly), _ = curve_fit(model_Bochard, nmC_poly[fit_i0_poly:fit_i1_poly], mD_poly[fit_i0_poly:fit_i1_poly])

    Dcc_solv = mD_solv/(nmC_solv*(1-nmC_solv))
    Dcc_poly = mD_solv/(nmC_solv*(1-nmC_solv))
    Ccc_solv = 1/(1-nmC_solv)
    Ccc_poly = 1/(1-nmC_poly)

    p_solv = np.polyfit(Ccc_solv[fit_i0_solv:fit_i1_solv], Dcc_solv[fit_i0_solv:fit_i1_solv], 1)
    p_poly = np.polyfit(Ccc_poly[fit_i0_poly:fit_i1_poly], Dcc_poly[fit_i0_poly:fit_i1_poly], 1)

    D0_solv, chi_solv = p_solv[0], p_solv[1]/(-2*p_solv[0])
    D0_poly, chi_poly = p_poly[0], p_poly[1]/(-2*p_poly[0])

    print(D0_solv, chi_solv)

    plot_D = True
    plot_D_type = 0

    print(fit_i0_solv,fit_i1_solv)

    # plot D
    if plot_D and plot_D_type == 0:
        fig,(ax0,ax1) = plt.subplots(1,2,figsize=(17,9))
        ax0.plot(mC_solv, mD_solv, ls="--", lw=4, color=(0,0.5,0), zorder=8)
        ax0.plot(mC_solv[fit_i0_solv:fit_i1_solv], mD_solv[fit_i0_solv:fit_i1_solv], lw=4, color=(0,0.5,0), zorder=9)
        ax0.plot(mC_solv, mD_solv, lw=8, color="white",   zorder=7)
        ax0.plot(mC_solv, model_Bochard(nmC_solv,D0_solv,chi_solv), color="crimson", lw=7, alpha=0.7, zorder=10)
        #ax0.scatter(*[[X[i] for i in (fit_i0_solv,fit_i1_solv,fit_i2_solv)] for X in (mC_solv,mD_solv)], s=100, color="darkgreen")
        ax1.plot(mC_poly, mD_poly, lw=4, color=(0,0,1), zorder=8)
        ax1.plot(mC_poly, mD_poly, lw=8, color="white", zorder=7)
        ax1.plot(mC_poly, model_Bochard(nmC_poly,D0_poly,chi_poly), color="crimson", lw=7, alpha=0.7, zorder=10)
        #ax1.scatter(*[[X[i] for i in (fit_i0_poly,fit_i1_poly)] for X in (mC_poly,mD_poly)], s=100, color="darkblue")
        for ax in (ax0,ax1):
            ax.grid()
            ax.set_xlabel("c")
            ax.set_ylabel("D")
            plotWhite(fig,ax)
        plt.subplots_adjust(left=0.05, right=0.95, wspace=0.2)
        #ax0.set_title("solvant")
        #ax1.set_title("polymer")
        ax0.set_xlim(0,5500)
        ax0.set_ylim(0,4.5e-5)
        ax1.set_xlim(0,8000)
        ax1.set_ylim(0,4e-5)
        #ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = {pntply(fit_para1_solv)}",      alpha=0.6)
        #ax0.plot(fitX2_solv, fitY2_solv, lw=3, color="orangered", label=f"D = {pntply(fit_para2_solv)}",      alpha=0.6)
        #ax1.plot(fitX1_poly, fitY1_poly, lw=3, color="red",       label=f"D = exp({pntply(fit_para1_poly)})", alpha=0.6)
        #ax0.legend(); ax1.legend()
        plt.show()
        plt.close() 


    if plot_D and plot_D_type == 1:
        nmC_solv, nmC_poly = mC_solv/c0_solv, mC_poly/c0_poly
        fig,(ax0,ax1) = plt.subplots(1,2,figsize=(20,9))
        ax0.plot(Ccc_solv, Dcc_solv, ls="--", lw=4, color=(0,0.5,0), zorder=8)
        ax0.plot(Ccc_solv[fit_i0_solv:fit_i1_solv], Dcc_solv[fit_i0_solv:fit_i1_solv], lw=4, color=(0,0.5,0), zorder=9)
        ax0.plot(Ccc_solv, Dcc_solv, lw=8, color="white",   zorder=7)
        ax0.plot(Ccc_solv, np.polyval(p_solv, Ccc_solv), color="crimson", lw=7, alpha=0.7, zorder=10)
        ax1.plot(Ccc_poly, Dcc_poly, lw=4, color=(0,0,1), zorder=8)
        ax1.plot(Ccc_poly, Dcc_poly, lw=8, color="white", zorder=7)
        ax1.plot(Ccc_poly, np.polyval(p_poly, Ccc_poly), color="crimson", lw=7, alpha=0.7, zorder=10)
        for ax in (ax0,ax1):
            ax.grid()
            ax.set_xlabel("1/(1-c/c0)")
            ax.set_ylabel("D/[(c/c0)(1-c/c0)]")
            plotWhite(fig,ax)
        plt.subplots_adjust(left=0.20, right=0.95, wspace=0.3)
        ax0.set_title("solvant")
        ax1.set_title("polymer")
        ax0.set_ylim(0,None)
        ax1.set_ylim(0,None)
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# _________________________________________________________________________________________________________________________________
               


def truc_2():
    T, Z, Cs_solv, Cs_poly = read2(*paths)

    # time indices
    mti0 = 5
    Dmti = 20
    interv = [(i,i+Dmti) for i in range(mti0,200-Dmti,1)]
    zi0 = 0
    zi1 = 64

    # new Xi lim
    xi0_solv = -0.008
    xi1_solv =  0.008
    xi0_poly = -0.007
    xi1_poly =  0.012

    poly_fix = 3

    for i, (mti0, mti1) in enumerate(interv):

        # c0 & ordinate fix
        c0_solv,  c0_poly  = get_c0(Z[mti0:mti1], Cs_solv[mti0:mti1]), get_c0(Z[mti0:mti1], Cs_poly[mti0:mti1])
        Zs_solv,  Zs_poly  = fix0_halfc0(Z, Cs_solv, Cs_poly, c0_solv, c0_poly, i0=15)
        Xis_solv, Xis_poly = ZtoXi(Zs_solv,T), ZtoXi(Zs_poly,T)

        # C(Xi) mean
        mXi, mC_solv, mC_poly = get_mean(Xis_solv,Xis_poly,Cs_solv,Cs_poly, mti0,mti1,zi0,zi1, zi1-zi0)

        Xi_solv = np.linspace(xi0_solv,xi1_solv,len(mXi))
        Xi_poly = np.linspace(xi1_poly,xi0_poly,len(mXi))
        mC_solv = np.interp(Xi_solv, mXi, mC_solv)
        mC_poly = np.interp(Xi_poly, mXi, mC_poly)

        # D derivation
        error_print_off()
        mD_solv = derive_D(Xi_solv, mC_solv, poly_fix)[0]
        mD_poly = derive_D(Xi_poly, mC_poly, poly_fix)[0]
        error_print_on()

        # fit D indices
        fit_i0_solv = 0
        fit_i1_solv = get_index(mC_solv, 4000)
        fit_i2_solv = -1
        fit_i0_poly = get_index(mC_poly,  200)
        fit_i1_poly = get_index(mC_poly, 5000)

        # fit X,Y
        try:
            fit_para1_solv = np.polyfit(mC_solv[fit_i0_solv:fit_i1_solv], mD_solv[fit_i0_solv:fit_i1_solv], 1)
            fit_para2_solv = np.polyfit(mC_solv[fit_i1_solv:fit_i2_solv], mD_solv[fit_i1_solv:fit_i2_solv], 1)
            fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly], np.log(mD_poly[fit_i0_poly:fit_i1_poly]), 2)
            fitX1_solv = np.linspace(0,                    mC_solv[fit_i1_solv], 100)
            fitX2_solv = np.linspace(mC_solv[fit_i1_solv], 5500,                 100)
            fitX1_poly = np.linspace(0,                    8000,                 100)
            fitY1_solv = np.polyval(fit_para1_solv, fitX1_solv)
            fitY2_solv = np.polyval(fit_para2_solv, fitX2_solv)
            fitY1_poly = np.exp(np.polyval(fit_para1_poly, fitX1_poly))
        except: continue

        # plot D
        fig,(ax0,ax1) = plt.subplots(1,2,figsize=(15,7))
        ax0.plot(mC_solv, mD_solv, lw=4, color="darkgreen", zorder=6)
        ax0.scatter(*[[X[i] for i in (fit_i0_solv,fit_i1_solv,fit_i2_solv)] for X in (mC_solv,mD_solv)], s=100, color="darkgreen", zorder=6)
        ax1.plot(mC_poly, mD_poly, lw=4, color="darkblue", zorder=6)
        ax1.scatter(*[[X[i] for i in (fit_i0_poly,fit_i1_poly)] for X in (mC_poly,mD_poly)], s=100, color="darkblue", zorder=6)
        for ax in (ax0,ax1):
            ax.grid()
            ax.set_xlabel("c")
            ax.set_ylabel("D")
        ax0.set_title("solvant")
        ax1.set_title("polymer")
        ax0.set_xlim(0,5500)
        ax0.set_ylim(0,4e-5)
        ax1.set_xlim(0,8000)
        ax1.set_ylim(0,5e-5)
        ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = {pntply(fit_para1_solv)}", alpha=0.6)
        ax0.plot(fitX2_solv, fitY2_solv, lw=3, color="orangered", label=f"D = {pntply(fit_para2_solv)}", alpha=0.6)
        ax1.plot(fitX1_poly, fitY1_poly, lw=3, color="red",       label=f"D = exp({pntply(fit_para1_poly)})", alpha=0.6)
        ax0.legend(); ax1.legend()
        plt.savefig(out_folder+f"/img_{i:03}.png")
        plt.close()
        print(f"{i+1:03}/{len(interv):03}:\t{mti0:03}-{mti1:03}")

    make_gif(out_folder+"/anim133_20.gif",True) # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# __________________________________________________________________________________________________________________________________



truc_1()