data22.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 = {}


### 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, k,x0): return np.exp(k*(x-x0))

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


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


path1 = "data22/dens_solvant_2_2.profile"
path2 = "data22/dens_polymer_2_2.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.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_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.002
    xi1_solv = 0.005
    xi0_poly = 0.000
    xi1_poly = 0.006

    # 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.007
    xi1_solv =  0.005
    xi0_poly = -0.005
    xi1_poly =  0.013

    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 = get_index(mC_solv,  500)
    fit_i1_solv = get_index(mC_solv, 3700)
    fit_i2_solv = -1
    fit_i0_poly = get_index(mC_poly,  200)
    fit_i1_poly = get_index(mC_poly, 3700)
    fit_i2_poly = get_index(mC_poly, 5500)
    fit_i3_poly = -1

    # fit X,Y
    fit_para1_solv = np.polyfit(mC_solv[fit_i0_solv:fit_i1_solv], np.log(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], 0)
    fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly], np.log(mD_poly[fit_i0_poly:fit_i1_poly]), 2)
    fit_para2_poly = np.polyfit(mC_poly[fit_i1_poly:fit_i2_poly], mD_poly[fit_i1_poly:fit_i2_poly], 0)
    fit_para3_poly = np.polyfit(mC_poly[fit_i2_poly:fit_i3_poly], np.log(mD_poly[fit_i2_poly:fit_i3_poly]), 1)
    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,                    mC_poly[fit_i1_poly], 100)
    fitX2_poly = np.linspace(mC_poly[fit_i1_poly], mC_poly[fit_i2_poly], 100)
    fitX3_poly = np.linspace(mC_poly[fit_i2_poly], 9500,                 100)
    fitY1_solv = np.exp(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))
    fitY2_poly = np.polyval(fit_para2_poly, fitX2_poly)
    fitY3_poly = np.exp(np.polyval(fit_para3_poly, fitX3_poly))

    plot_D = True

    # plot D
    if plot_D:
        fig,(ax0,ax1) = plt.subplots(1,2)
        ax0.plot(mC_solv, mD_solv, lw=4, color="darkgreen")
        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="darkblue")
        ax1.scatter(*[[X[i] for i in (fit_i0_poly,fit_i1_poly,fit_i2_poly,fit_i3_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")
        ax0.set_title("solvant")
        ax1.set_title("polymer")
        ax0.set_xlim(0,5500)
        ax0.set_ylim(0,None)
        ax1.set_xlim(0,9500)
        ax1.set_ylim(0,None)
        ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = exp({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)
        ax1.plot(fitX2_poly, fitY2_poly, lw=3, color="orangered", label=f"D = {pntply(fit_para2_poly)}",      alpha=0.6)
        ax1.plot(fitX3_poly, fitY3_poly, lw=3, color="goldenrod", label=f"D = exp({pntply(fit_para3_poly)})", alpha=0.6)
        ax0.legend(); ax1.legend()
        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.007
    xi1_solv =  0.005
    xi0_poly = -0.005
    xi1_poly =  0.013

    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 = get_index(mC_solv,  500)
        fit_i1_solv = get_index(mC_solv, 3900)
        fit_i2_solv = -1
        fit_i0_poly = get_index(mC_poly,  500)
        fit_i1_poly = get_index(mC_poly, 4000)
        fit_i2_poly = get_index(mC_poly, 5500)
        fit_i3_poly = -1

        # fit X,Y
        try:
            fit_para1_solv = np.polyfit(mC_solv[fit_i0_solv:fit_i1_solv], np.log(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], 0)
            fit_para1_poly = np.polyfit(mC_poly[fit_i0_poly:fit_i1_poly], np.log(mD_poly[fit_i0_poly:fit_i1_poly]), 2)
            fit_para2_poly = np.polyfit(mC_poly[fit_i1_poly:fit_i2_poly], mD_poly[fit_i1_poly:fit_i2_poly], 0)
            fit_para3_poly = np.polyfit(mC_poly[fit_i2_poly:fit_i3_poly], np.log(mD_poly[fit_i2_poly:fit_i3_poly]), 1)
            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,                    mC_poly[fit_i1_poly], 100)
            fitX2_poly = np.linspace(mC_poly[fit_i1_poly], mC_poly[fit_i2_poly], 100)
            fitX3_poly = np.linspace(mC_poly[fit_i2_poly], 9500,                 100)
            fitY1_solv = np.exp(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))
            fitY2_poly = np.polyval(fit_para2_poly, fitX2_poly)
            fitY3_poly = np.exp(np.polyval(fit_para3_poly, fitX3_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,fit_i2_poly,fit_i3_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,1.2e-5)
        ax1.set_xlim(0,9500)
        ax1.set_ylim(0,8e-5)
        ax0.plot(fitX1_solv, fitY1_solv, lw=3, color="red",       label=f"D = exp({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)
        ax1.plot(fitX2_poly, fitY2_poly, lw=3, color="orangered", label=f"D = {pntply(fit_para2_poly)}",      alpha=0.6)
        ax1.plot(fitX3_poly, fitY3_poly, lw=3, color="goldenrod", label=f"D = exp({pntply(fit_para3_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+"/anim22_20.gif",True) # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
# __________________________________________________________________________________________________________________________________



truc_1()