data11bis.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

# decorator for Tuck models
def deco_model_Xi(model_Xi):
    @np.vectorize
    def new_model(c,n):
        if 0 <= c <= 0.5: return  model_Xi(2*c,n)
        else:             return -model_Xi(2*(1-c),n)
    return new_model


### Models

# curve fit models
def model_erf(xi, c0,D): return c0/2*(1+erf(xi/(2*D**0.5)))
def model_exp(x, k,x0): return np.exp(k*(x-x0))

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

# Tuck models
@deco_model_Xi
def model_Xi1(c,n): return np.exp(n)-np.exp(n*c)
@deco_model_Xi
def model_Xi2(c,n): return 1-c**n
@deco_model_Xi
def model_Xi3(c,n): return (1-c)**n
@deco_model_Xi
def model_Xi4(c,n): return np.cos(1/2*np.pi*c)
@deco_model_Xi
def model_Xi5(c,n): return 1-np.sin(1/2*np.pi*c)

def model_D1(c,n): return 1/2*np.exp(n*c)*(1-np.exp(n*c)+n*c*np.exp(n))
def model_D2(c,n): return 1/2*n*c**2*(1-c**n/(1+n))
def model_D3(c,n): return n/2*(n+1)*((1-c)**(n-1)-(1-c)**(2*n))
def model_D4(c,n): return 1/2*np.sin(1/2*np.pi*c)**2
def model_D5(c,n): return 1/2*np.cos(1/2*np.pi*c)*(np.cos(1/2*np.pi*c)+1/2*np.pi*c-1)


### 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:
        fix = np.polyfit(C[:len(C)//5],Xii[:len(Xi)//5],min(len(C)//5,polynomial_fix))
        Xii -= np.polyfit(C[:len(C)//5],Xii[:len(Xi)//5],min(len(C)//5,polynomial_fix))[-1]
    D = -0.5*Xip*Xii
    return D, Xip, Xii, fix


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 = "data11/dens_solvant_1_1.profile"
path2 = "data11/dens_polymer_1_1.profile"
paths = (path1,path2)
out_folder = "data_out"


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

    ti0 = 50
    ti1 = 149
    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), zorder=10)
        ax.plot(Z[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(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_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 = 50
    ti1 = 149

    mti0 = 50 #50
    mti1 = 70 #70

    zi0 = 0
    zi1 = 64

    # regression xi intervals
    xi0_solv =  -0.01
    xi1_solv =   0.003
    xi0_poly =  -0.005
    xi1_poly =   0.005

    # 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 = True

    # 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_Xi = 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=1, zorder=13)
            #ax.plot(mXi, rC_poly, lw=7, color="crimson", alpha=0.7, zorder=13)
        ax.set_xlabel("η = z/√t")
        ax.set_ylabel("c")
        plotWhite(fig,ax)
        plt.show()
        plt.close()
        
    # plot Xi
    if plot_Xi:
        fig,ax = plt.subplots()
        for ti in range(ti0,ti1):
            c = (ti-ti0)/(ti1-ti0)
            ax.plot(Cs_solv[ti,zi0:zi1], Xis_solv[ti,zi0:zi1], lw=3, color=(1-c,1-0.5*c,0,.5), zorder=10)
            ax.plot(Cs_poly[ti,zi0:zi1], Xis_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(mC_solv, mXi, ls='--', lw=1.7, color="black", zorder=11)
        ax.plot(mC_poly, mXi, ls='--', lw=1.7, color="white", zorder=11)
        ax.set_ylabel("η")
        ax.set_xlabel("c")
        ax.set_ylim(-0.02,0.02)
        plotWhite(fig,ax)
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~


    # new Xi lim
    xi0_solv = -0.0123
    xi1_solv =  0.0060
    xi0_poly = -0.0050
    xi1_poly =  0.0100

    old_mC_solv = mC_solv

    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)

    poly_fix = 4

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

    X = np.linspace(0,6000,200)
    mY = np.polyval(parafixms,X) - parafixms[-1]
    if do_regress: rY = np.polyval(parafixrs,X) - parafixrs[-1] #                           !!!!!!!!!!!!!!!!!!!!!!!!

    plot_D_exp = False

    # plot D exp derivation
    if plot_D_exp:
        fig,((axs0,axs1,axs2,axs3)) = plt.subplots(1,4,figsize=(25,6))
        axs0.plot(old_mC_solv, mXi, lw=3.5, ls="--", color="white")
        axs0.plot(mC_solv, Xi_solv,   lw=4, color="white", label="η(c)" )
        axs1.plot(mC_solv, mXip_solv, lw=4, color="white", label="dη/dc")
        axs2.plot(mC_solv, mXii_solv, lw=4, color="white", label="∫dcη" )
        axs3.plot(mC_solv, mD_solv,   lw=4, color="white", label="D(c)" )
        plt.subplots_adjust(left=0.05, right=0.95, wspace=0.5)
        for ax in (axs0,axs1,axs2,axs3):
            ax.grid()
            ax.set_xlabel("c")
            plotWhite(fig,ax)
        axs0.set_ylabel("η(c)"); axs1.set_ylabel("dη/dc"); axs2.set_ylabel("∫dcη"); axs3.set_ylabel("D(c)"); 
        axs3.set_ylim(0,None)
        axs2.plot(X,mY, color="crimson", ls=':', lw=3, alpha=0.7)
        
        plt.show()
        plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
    
    # fit D indices
    fit_i0_solv = 0
    fit_i1_solv = get_index(mC_solv, 3500)
    fit_i2_solv = get_index(mC_solv, 5500)
    fit_i0_poly = get_index(mC_poly,  500)
    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],   0)
    #fit_para2_solv = np.polyfit(mC_solv[fit_i1_solv:fit_i2_solv+1], np.log(mD_solv[fit_i1_solv:fit_i2_solv+1]), 2)
    #fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly+1], np.log(mD_poly[fit_i0_poly:fit_i1_poly+1]), 2)
    #print(fit_para1_poly)
    #fitX1_solv = np.linspace(0, mC_solv[fit_i1_solv],100)
    #fitX2_solv = np.linspace(mC_solv[fit_i1_solv],6000,100)
    #fitX1_poly = np.linspace(0,7000,100)
    #fitY1_solv = np.polyval(fit_para1_solv, fitX1_solv)
    #fitY2_solv = np.exp(np.polyval(fit_para2_solv, fitX2_solv))
    #fitY1_poly = np.exp(np.polyval(fit_para1_poly, fitX1_poly)) # np.polyval(fit_para1_poly, fitX1_poly)

    fit_i0_solv, fit_i1_solv = 25, 53
    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,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,199-Dmti,1)]
    zi0 = 0
    zi1 = 64

    # new Xi lim
    xi0_solv = -0.0123
    xi1_solv =  0.0060
    xi0_poly = -0.0050
    xi1_poly =  0.0100

    poly_fix = 3

    Dflat1_solv = []
    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, 3500)
        fit_i2_solv = get_index(mC_solv, 5000)
        fit_i0_poly = get_index(mC_poly,  500)
        fit_i1_poly = get_index(mC_poly, 5000)
        #print(fit_i0_solv, fit_i1_solv, fit_i2_solv, fit_i0_poly, fit_i1_poly)

        # 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],   0); Dflat1_solv.append(fit_para1_solv[0])
            fit_para2_solv = np.polyfit(mC_solv[fit_i1_solv:fit_i2_solv+1], np.log(mD_solv[fit_i1_solv:fit_i2_solv+1]), 2)
            fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly+1], np.log(mD_poly[fit_i0_poly:fit_i1_poly+1]), 2)
            fitX1_solv = np.linspace(0,                    mC_solv[fit_i1_solv], 100)
            fitX2_solv = np.linspace(mC_solv[fit_i1_solv], 6000,                 100)
            fitX1_poly = np.linspace(0,                    7000,                 100)
            fitY1_solv = np.polyval(fit_para1_solv, fitX1_solv); flatY1_solv = np.polyval([np.mean(Dflat1_solv)], fitX1_solv)
            fitY2_solv = np.exp(np.polyval(fit_para2_solv, fitX2_solv))
            fitY1_poly = np.exp(np.polyval(fit_para1_poly, fitX1_poly))
        except: continue

        # plot
        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,6000)
        ax0.set_ylim(0,7e-5)
        ax1.set_xlim(0,7000)
        ax1.set_ylim(0,1.75e-5)
        ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = {pntply(fit_para1_solv)}",            alpha=0.6)
        ax0.plot(fitX1_solv, flatY1_solv,lw=3, color="salmon",    label=f"Dmoy = {pntply([np.mean(Dflat1_solv)])}", alpha=0.5, zorder=5)
        ax0.plot(fitX2_solv, fitY2_solv, lw=3, color="orangered", label=f"D = exp({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+"/anim11_20.gif",True) # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# __________________________________________________________________________________________________________________________________



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

    # time indices
    ti0 = 99
    ti1 = 199

    mti0 = 179
    mti1 = 199

    zi0 = 0
    zi1 = 64

    # regression xi intervals
    xi0_solv = -0.012
    xi1_solv =  0.000
    xi0_poly =  0.005
    xi1_poly =  0.010

    # 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)
    Cs_solv  /= c0_solv; Cs_poly  /= c0_poly # normalization
    Xis_solv /= c0_solv; Xis_poly /= c0_poly # normalization

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

    # Tuck models regress


    # plot C
    fig,ax = plt.subplots(figsize=(15,7))
    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.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)
    ax.grid()
    ax.set_xlabel("η /c0")
    ax.set_ylabel("c/c0")
    plt.show()
    plt.close() # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# ___________________________________________________________________________________________________________________________________


truc_1()