def plot_tr(vol_tr,st_point_dep1,st_point_ada1,dep_vec1_2,ada_vec1_2,bet,gam,alpha3,n_sp,r,c):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load

    import pyabf
    cell_num = '95810005'
    abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
    abf.setSweep(15)
    vol_base = -62.86698717948718
    abf.sweepY = abf.sweepY + vol_base

    t = sym.Symbol('t')
    alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load()
    #V=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*((-2)*sym.exp(1)**(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta*((-1)+gamma)+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(1/2)*((-1)+alpha+beta)*((-1)+beta+gamma**2)+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)*V0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+alpha*((-1)+(1+(-1)*beta)**(1/2)+beta)*((-1)+beta+gamma**2)+(-1)*((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+((-1)+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(-1)*((-1)+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2)))+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(alpha*(1+(1+(-1)*beta)**(1/2)+(-1)*beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+(-1)*(1+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(1+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2))))
    #Iadap=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*(2*sym.exp(1) **(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*alpha*(1+(-1)*beta)**(1/2)*((-1)+beta+gamma**2)+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*alpha*(1+(-1)*beta)**(1/2)+(-2)*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+IdA0*beta*gamma+(-1)*IdA0*beta**2*gamma+(-1)*gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+beta*gamma**2+(-1)*IaA0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2)))+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)+(-1)*alpha*(1+(-1)*beta)**(1/2)+2*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+(-1)*beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+(-1)*IdA0*beta*gamma+IdA0*beta**2*gamma+gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+(-1)*beta*gamma**2+(-1)*IaA0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2))));
    #Idep=sym.exp(1)**(((-1)*t+t0)*gamma)*IdA0;





    t0_val=0
    vtm=np.mean(vol_tr)/66.35
    vrm=-1.0
    d_in=st_point_dep1[r]
    a_in=st_point_ada1[r]
    aux = V.subs(alpha, alpha3[1]).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0,d_in)
    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = abf.sweepX * 1000 / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]

    ada = Iadap.subs(alpha, alpha3[1]).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)
    plt.figure()
    plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


    vrm=-0.7349960338647149
    d_inc=dep_vec1_2[r,c]
    a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep ****'+str(i)+' *****spike****')
        print(dep + d_inc)
        print('***** ada ****'+str(i)+' *****spike****')
        print(ada + a_inc)


        aux = V.subs(alpha, alpha3[1]).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = abf.sweepX[(abf.sweepX * 1000 / 15.58)>t_init]* 1000 / 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:
            t_next = x_vals[y_vals > vtm][0]


            t_next=x_vals[y_vals>vtm][0]
            ada = Iadap.subs(alpha, alpha3[1]).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        else:
            i=n_sp+1

def plot_tr2(vtm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/0.252587)+10
    #cell_num = '95810005'
    abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
    abf.setSweep(tr)
    vol_base = -62.86698717948718
    abf.sweepY = abf.sweepY + vol_base

    t = sym.Symbol('t')
    alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load()
    #V=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*((-2)*sym.exp(1)**(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta*((-1)+gamma)+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(1/2)*((-1)+alpha+beta)*((-1)+beta+gamma**2)+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)*V0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+alpha*((-1)+(1+(-1)*beta)**(1/2)+beta)*((-1)+beta+gamma**2)+(-1)*((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+((-1)+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(-1)*((-1)+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2)))+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(alpha*(1+(1+(-1)*beta)**(1/2)+(-1)*beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+(-1)*(1+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(1+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2))))
    #Iadap=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*(2*sym.exp(1) **(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*alpha*(1+(-1)*beta)**(1/2)*((-1)+beta+gamma**2)+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*alpha*(1+(-1)*beta)**(1/2)+(-2)*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+IdA0*beta*gamma+(-1)*IdA0*beta**2*gamma+(-1)*gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+beta*gamma**2+(-1)*IaA0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2)))+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)+(-1)*alpha*(1+(-1)*beta)**(1/2)+2*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+(-1)*beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+(-1)*IdA0*beta*gamma+IdA0*beta**2*gamma+gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+(-1)*beta*gamma**2+(-1)*IaA0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2))));
    #Idep=sym.exp(1)**(((-1)*t+t0)*gamma)*IdA0;





    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0,d_in)
    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = abf.sweepX * 1000 / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]
    t_out.append(t_next)

    ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)
    plt.figure()
    plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


    vrm=-0.7349960338647149
    #d_inc=dep_vec1_2[r,c]
    #a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep ****'+str(i)+' *****spike****')
        print(dep + d_inc)
        print('***** ada ****'+str(i)+' *****spike****')
        print(ada + a_inc)


        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = abf.sweepX[(abf.sweepX * 1000 / 15.58)>t_init]* 1000 / 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)
            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        else:
            i=n_sp+1

    plt.plot(abf.sweepX * 1000, abf.sweepY)
    return t_out


def plot_tr3(vtm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    tr=min(int(alpha3/0.252587)+10,15)
    #cell_num = '95810005'
    abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
    abf.setSweep(tr)
    vol_base = -62.86698717948718
    abf.sweepY = abf.sweepY + vol_base

    t = sym.Symbol('t')
    alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load()
    #V=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*((-2)*sym.exp(1)**(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta*((-1)+gamma)+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(1/2)*((-1)+alpha+beta)*((-1)+beta+gamma**2)+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)*V0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+alpha*((-1)+(1+(-1)*beta)**(1/2)+beta)*((-1)+beta+gamma**2)+(-1)*((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+((-1)+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(-1)*((-1)+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2)))+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(alpha*(1+(1+(-1)*beta)**(1/2)+(-1)*beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+(-1)*(1+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(1+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2))))
    #Iadap=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*(2*sym.exp(1) **(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*alpha*(1+(-1)*beta)**(1/2)*((-1)+beta+gamma**2)+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*alpha*(1+(-1)*beta)**(1/2)+(-2)*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+IdA0*beta*gamma+(-1)*IdA0*beta**2*gamma+(-1)*gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+beta*gamma**2+(-1)*IaA0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2)))+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)+(-1)*alpha*(1+(-1)*beta)**(1/2)+2*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+(-1)*beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+(-1)*IdA0*beta*gamma+IdA0*beta**2*gamma+gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+(-1)*beta*gamma**2+(-1)*IaA0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2))));
    #Idep=sym.exp(1)**(((-1)*t+t0)*gamma)*IdA0;





    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0,d_in)
    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = abf.sweepX * 1000 / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_out=[]
    if len(np.nonzero(y_vals > vtm)[0]) > 0:
        t_next=x_vals[y_vals>vtm][0]
        t_out.append(t_next)
        ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, vrm).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
        dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


        #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
        #    spt = pickle.load(f)

        #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
        #    Iadap_mat = pickle.load(f)

        #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
        #    Idep_mat = pickle.load(f)

        #auxa = (IaA - a_in) ** 2
        #auxd = (IdA - d_in) ** 2
        #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
        #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
        x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


        x_vals = x_sector[np.nonzero(x_sector)]
        y_vals = lam_x(x_vals)
        plt.figure()
        plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        vrm=-0.7349960338647149
        #d_inc=dep_vec1_2[r,c]
        #a_inc=ada_vec1_2[r,c]
        for i in range(n_sp):
            t_init=t_next+2/15.58
            print('******dep ****'+str(i)+' *****spike****')
            print(dep + d_inc)
            print('***** ada ****'+str(i)+' *****spike****')
            print(ada + a_inc)


            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
            #print('ada comp')
            #print( ada+a_inc)
            #print('dep comp')
            #print(dep + d_inc)

    #        print('ada')
    #        print(ada_nc)
    #        print('dep')
    #        print(dep_nc)
            lam_x = sym.lambdify(t, aux, modules=['numpy'])
            x_vals = abf.sweepX[(abf.sweepX * 1000 / 15.58)>t_init]* 1000 / 15.58
            y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
            t_init = t_next + 2 / 15.58

            if len(np.nonzero(y_vals > vtm)[0]) > 0:

                t_next=x_vals[y_vals>vtm][0]
                t_out.append(t_next)
                ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
                dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
    #        auxa = (IaA - a_in) ** 2
    #        auxd = (IdA - d_in) ** 2
    #        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
    #        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
                x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


                x_vals = x_sector[np.nonzero(x_sector)]
                y_vals = lam_x(x_vals)
                #plt.figure()
                plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
            else:
                i=n_sp+1

        plt.plot(abf.sweepX * 1000, abf.sweepY)
    else:
        plt.figure()
        plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        plt.plot(abf.sweepX * 1000, abf.sweepY)
    return t_out




def plot_trnew(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/0.252587)+10
    #cell_num = '95810005'
    abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
    abf.setSweep(tr)
    vol_base = -62.86698717948718
    abf.sweepY = abf.sweepY + vol_base

    t = sym.Symbol('t')
    alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load()
    #V=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*((-2)*sym.exp(1)**(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta*((-1)+gamma)+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(1/2)*((-1)+alpha+beta)*((-1)+beta+gamma**2)+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)*V0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+alpha*((-1)+(1+(-1)*beta)**(1/2)+beta)*((-1)+beta+gamma**2)+(-1)*((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+((-1)+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(-1)*((-1)+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2)))+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(alpha*(1+(1+(-1)*beta)**(1/2)+(-1)*beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+(-1)*(1+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(1+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2))))
    #Iadap=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*(2*sym.exp(1) **(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*alpha*(1+(-1)*beta)**(1/2)*((-1)+beta+gamma**2)+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*alpha*(1+(-1)*beta)**(1/2)+(-2)*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+IdA0*beta*gamma+(-1)*IdA0*beta**2*gamma+(-1)*gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+beta*gamma**2+(-1)*IaA0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2)))+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)+(-1)*alpha*(1+(-1)*beta)**(1/2)+2*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+(-1)*beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+(-1)*IdA0*beta*gamma+IdA0*beta**2*gamma+gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+(-1)*beta*gamma**2+(-1)*IaA0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2))));
    #Idep=sym.exp(1)**(((-1)*t+t0)*gamma)*IdA0;





    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = abf.sweepX * 1000 / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]
    t_out.append(t_next)

    ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)
    plt.figure()
    plt.plot(31.1 + x_vals * 15.58, y_vals * Vconvfact)


    #vrm=-0.7349960338647149
    #d_inc=dep_vec1_2[r,c]
    #a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep ****'+str(i)+' *****spike****')
        print(dep + d_inc)
        print('***** ada ****'+str(i)+' *****spike****')
        print(ada + a_inc)


        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = abf.sweepX[(abf.sweepX * 1000 / 15.58)>t_init]* 1000 / 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)
            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(31.1 + x_vals * 15.58, y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)

            i=n_sp+1

    #plt.plot(abf.sweepX * 1000, abf.sweepY)
    return t_out



def plot_trnew2(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load

    st_sign = 531.1
    corr=int(alpha3/0.252587)*0.2
    f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

    aux=f.readline()
    tim=[]
    vol=[]
    trovato = False
    for aux in f:
        tim.append(np.float64(aux.split('\t')[0]))
        vol.append(np.float64(aux.split('\t')[1]))
    tim = np.array(tim)
    vol = np.array(vol)

    tim_aux=np.linspace(tim.min(), tim.max(), 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/0.252587)+10
    #cell_num = '95810005'


    t = sym.Symbol('t')
    alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load()
    #V=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*((-2)*sym.exp(1)**(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta*((-1)+gamma)+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(1/2)*((-1)+alpha+beta)*((-1)+beta+gamma**2)+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)*V0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+alpha*((-1)+(1+(-1)*beta)**(1/2)+beta)*((-1)+beta+gamma**2)+(-1)*((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+((-1)+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(-1)*((-1)+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2)))+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(alpha*(1+(1+(-1)*beta)**(1/2)+(-1)*beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+((-1)+beta)*((IaA0+(-1)*IdA0)*beta**2+(-1)*(1+(1+(-1)*beta)**(1/2))*beta*((-1)+IdA0*((-1)+gamma))+(1+(1+(-1)*beta)**(1/2))*((-1)+gamma**2)+IaA0*beta*((-1)+gamma**2))))
    #Iadap=(1/2)*sym.exp(1)**((-1)*t0*(1+(-1)*beta)**(1/2)+(-1)*t*((1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*beta)**(-3/2)*((-1)+beta+gamma**2)**(-1)*(2*sym.exp(1) **(t*(1+(-1)*beta)**(1/2)+t0*((1+(-1)*beta)**(1/2)+gamma))*IdA0*(1+(-1)*beta)**(3/2)*beta+(-2)*sym.exp(1)**(t0*(1+(-1)*beta)**(1/2)+t*((1+(-1)*beta)**(1/2)+gamma))*alpha*(1+(-1)*beta)**(1/2)*((-1)+beta+gamma**2)+sym.exp(1)**(t*(2*(1+(-1)*beta)**(1/2)+gamma))*(1+(-1)*alpha*(1+(-1)*beta)**(1/2)+(-2)*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+IdA0*beta*gamma+(-1)*IdA0*beta**2*gamma+(-1)*gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+beta*gamma**2+(-1)*IaA0*((-1)+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+(-1)*V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2)))+sym.exp(1)**(2*t0*(1+(-1)*beta)**(1/2)+t*gamma)*((-1)+(-1)*alpha*(1+(-1)*beta)**(1/2)+2*beta+(-1)*IdA0*(1+(-1)*beta)**(1/2)*beta+alpha*(1+(-1)*beta)**(1/2)*beta+(-1)*beta**2+IdA0*(1+(-1)*beta)**(1/2)*beta**2+(-1)*IdA0*beta*gamma+IdA0*beta**2*gamma+gamma**2+alpha*(1+(-1)*beta)**(1/2)*gamma**2+(-1)*beta*gamma**2+(-1)*IaA0*(1+(1+(-1)*beta)**(1/2))*((-1)+beta)*((-1)+beta+gamma**2)+V0*(1+beta**2+(-1)*gamma**2+beta*((-2)+gamma**2))));
    #Idep=sym.exp(1)**(((-1)*t+t0)*gamma)*IdA0;





    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = tim_aux / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]
    t_out.append(t_next)

    ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)

    plt.figure()
    plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)


    #vrm=-0.7349960338647149
    #d_inc=dep_vec1_2[r,c]
    #a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep ****'+str(i)+' *****spike****')
        print(dep + d_inc)
        print('***** ada ****'+str(i)+' *****spike****')
        print(ada + a_inc)


        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)
            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)

            i=n_sp+1

    plt.plot(tim, vol)
    return t_out


def plot_tr_inter(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3,intervalli):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load

    tim_aux = np.linspace(0.0, intervalli[len(intervalli)-1][1], 10000)
    plt.figure()

    for ind_int in range(len(intervalli)):
        st_sign = intervalli[ind_int][0]
        len_sign = (intervalli[ind_int][1]-intervalli[ind_int][0])/15.58
        corr=int(alpha3/0.252587)*0.2


        #alpha3 = [0.252587 * 4, 0.252587 * 5]
        t_out = []
        tr=int(alpha3/0.252587)+10
        #cell_num = '95810005'


        t = sym.Symbol('t')
        alpha,beta,gamma,IaA0,IdA0,t0,V0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0')

        [V, Iadap, Idep] = load()


        if ind_int==0:
            t0_val = 0
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
        else:
            t0_val = t_next#(intervalli[ind_int][1]-intervalli[ind_int-1][1])/15.58
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, y_vals[len(y_vals) - 1]).subs(IaA0,a_in).subs(IdA0, d_in)

        print('***** dep init ****')
        print(d_in)
        print('***** ada init ****')
        print(a_in)

        lam_x = sym.lambdify(t, aux, modules=['numpy'])

        x_vals = tim_aux[(tim_aux / 15.58) > t0_val] / 15.58

        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        if len(x_vals[y_vals>vtm])>0:
            t_next=t0_val+min(x_vals[y_vals>vtm][0]-t0_val,len_sign)
        else:
            t_next=t0_val+len_sign
        t_out.append(t_next)

        ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
        dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


        x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


        x_vals = x_sector[np.nonzero(x_sector)]
        y_vals = lam_x(x_vals)

        #plt.figure()
        if ind_int==0:
            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            plt.plot( x_vals * 15.58, y_vals * Vconvfact)

        stop=False
        i=0
        while(not(stop)):
            t_init = t_next + 2 / 15.58

            print('******dep ****'+str(i)+' *****spike****')
            print(dep + d_inc)
            print('***** ada ****'+str(i)+' *****spike****')
            print(ada + a_inc)
            if ind_int==0:
                aux1=0
            else:
                aux1=(st_sign/15.58)

            if t_init < aux1+len_sign:


                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
                lam_x = sym.lambdify(t, aux, modules=['numpy'])
                x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
                y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
                t_init = t_next + 2 / 15.58

                if len(x_vals[y_vals > vtm])>0:
                    t_next=t0_val+min(x_vals[y_vals>vtm][0],len_sign)
                else:
                    t_next =t0_val+ len_sign
                t_out.append(t_next)
                ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
                dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
                x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


                x_vals = x_sector[np.nonzero(x_sector)]
                y_vals = lam_x(x_vals)
            #plt.figure()

                if ind_int == 0:
                    plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
                else:
                    plt.plot(x_vals * 15.58, y_vals * Vconvfact)

                i=i+1
            else:
                #tim_aux = np.linspace(tim.min(), tim.max(), 10000)
                if (ind_int<len(intervalli)-1):
                    t_init=intervalli[ind_int][1]/15.58
                    t_next=intervalli[ind_int+1][0]/15.58
                    print('si riparte da................')
                    print(y_vals[len(y_vals)-1])
                    print('....................................')
                    aux = V.subs(alpha, 0.0).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, y_vals[len(y_vals)-1]).subs(IaA0,ada).subs(IdA0, dep)
                    ada = Iadap.subs(alpha, 0).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, y_vals[len(y_vals) - 1]).subs(IaA0, ada).subs(IdA0, dep).subs(t, t_next)
                    dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep).subs(t, t_next)

                    lam_x = sym.lambdify(t, aux, modules=['numpy'])
                    x_vals = tim_aux / 15.58
                    x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

                    x_vals = x_sector[np.nonzero(x_sector)]

                    y_vals = lam_x(x_vals)
                # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
                    #t_next = (intervalli[ind_int+1][0]-st_sign)/15.58
                    t_out.append(t_next)


                    #plt.figure()
                    plt.plot( x_vals * 15.58, y_vals * Vconvfact)

                    print('aoa')
                    print(n_sp)
                    print(len(np.nonzero(y_vals > vtm)[0]))
                    print(vtm)
                i = i + 1
                stop = True

#        plt.plot(tim, vol)
    return t_out


def plot_tr_inter2(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,gam,n_sp,cell_num,alpha3,intervalli,tc):
    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load

    tim_aux = np.linspace(0.0, intervalli[len(intervalli)-1][1], 10000)
    plt.figure()
    t_out = []

    for ind_int in range(len(intervalli)):
        st_sign = intervalli[ind_int][0]
        len_sign = (intervalli[ind_int][1]-intervalli[ind_int][0])/15.58
        corr=int(alpha3/0.252587)*0.2


        #alpha3 = [0.252587 * 4, 0.252587 * 5]

        tr=int(alpha3/0.252587)+10
        #cell_num = '95810005'


        t = sym.Symbol('t')
        alpha,beta,gamma,IaA0,IdA0,t0,V0,ada0,dep0= sym.symbols('alpha,beta,gamma,IaA0,IdA0,t0,V0,ada0,dep0')

        [V, Iadap, Idep] = load()

        V2 = -1 + (-(V0) - 1) * sym.exp(1) ** (-300 * (t-t0))

        if ind_int==0:
            t0_val = 0
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
            #V2.subs(t,t0).subs(V0,-1) -Vconvfact + (-(V0 * Vconvfact) - Vconvfact) * sym.exp(1) ** (-300 * t)
            #aux=V2.subs(t0, t0_val).subs(V0, -1)
        else:
            t0_val = t_next#(intervalli[ind_int][1]-intervalli[ind_int-1][1])/15.58
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, y_vals[len(y_vals) - 1]).subs(IaA0,a_in).subs(IdA0, d_in)

        print('***** dep init ****')
        print(d_in)
        print('***** ada init ****')
        print(a_in)

        lam_x = sym.lambdify(t, aux, modules=['numpy'])

        x_vals = tim_aux[(tim_aux / 15.58) > t0_val] / 15.58

        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        if len(x_vals[y_vals>vtm])>0:
            t_next=t0_val+min(x_vals[y_vals>vtm][0]-t0_val,len_sign)
        else:
            t_next=t0_val+len_sign
        t_out.append(t_next)

        ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
        dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)


        x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


        x_vals = x_sector[np.nonzero(x_sector)]
        y_vals = lam_x(x_vals)

        #plt.figure()
        if ind_int==0:
            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            plt.plot( x_vals * 15.58, y_vals * Vconvfact)

        stop=False
        i=0
        while(not(stop)):
            t_init = t_next + 2 / 15.58

            print('******dep ****'+str(i)+' *****spike****')
            print(dep + d_inc)
            print('***** ada ****'+str(i)+' *****spike****')
            print(ada + a_inc)
            if ind_int==0:
                aux1=0
            else:
                aux1=(st_sign/15.58)

            if t_init < aux1+len_sign:


                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
                lam_x = sym.lambdify(t, aux, modules=['numpy'])
                x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
                y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
                t_init = t_next + 2 / 15.58

                if len(x_vals[y_vals > vtm])>0:
                    t_next=t0_val+min(x_vals[y_vals>vtm][0]-t0_val,len_sign)
                else:
                    t_next =t0_val+ len_sign
                t_out.append(t_next)
                ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
                dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
                x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


                x_vals = x_sector[np.nonzero(x_sector)]
                y_vals = lam_x(x_vals)
            #plt.figure()

                if ind_int == 0:
                    plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
                else:
                    plt.plot(x_vals * 15.58, y_vals * Vconvfact)

                i=i+1
            else:
                #tim_aux = np.linspace(tim.min(), tim.max(), 10000)
                if (ind_int<len(intervalli)-1):
                    t_init=intervalli[ind_int][1]/15.58
                    t_next=intervalli[ind_int+1][0]/15.58
                    print('si riparte da................')
                    print(y_vals[len(y_vals)-1])
                    print('....................................')
                    aux = V.subs(alpha, 0.0).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, y_vals[len(y_vals)-1]).subs(IaA0,ada).subs(IdA0, dep)
                    V2 = -1 + (V0 + 1) * sym.exp(1) ** (- (t - t0)/tc)
                    Iadap2=a_in+(ada0-a_in) * sym.exp(1) ** (- (t - t0)/tc)
                    Idep2=d_in+(dep0-d_in) * sym.exp(1) ** (- (t - t0)/tc)

                    aux = V2.subs(t0, t_init).subs(V0, y_vals[len(y_vals) - 1])
                    ada = Iadap2.subs(t0, t_init).subs(ada0,ada).subs(t, t_next)
                    dep = Idep2.subs(t0, t_init).subs(dep0, dep).subs(t, t_next)

                    #ada = Iadap.subs(alpha, 0).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, y_vals[len(y_vals) - 1]).subs(IaA0, ada).subs(IdA0, dep).subs(t, t_next)
                    #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep).subs(t, t_next)

                    lam_x = sym.lambdify(t, aux, modules=['numpy'])
                    x_vals = tim_aux / 15.58
                    x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

                    x_vals = x_sector[np.nonzero(x_sector)]
                    y_vals = np.zeros(len(x_vals))

                    for indic in range(len(x_vals)):
                        y_vals[indic]=aux.subs(t,x_vals[indic])


                # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
                    #t_next = (intervalli[ind_int+1][0]-st_sign)/15.58
                    t_out.append(t_next)


                    #plt.figure()
                    if (np.size(y_vals) == 1):
                        y_vals = np.ones(np.size(x_vals)) * y_vals

                    plt.plot( x_vals * 15.58, y_vals * Vconvfact)

                    print('aoa')
                    print(n_sp)
                    print(len(np.nonzero(y_vals > vtm)[0]))
                    print(vtm)
                i = i + 1
                stop = True

#        plt.plot(tim, vol)
    return t_out



def plot_tr_v3(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,delta1,n_sp,cell_num,alpha3,sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    st_sign = 531.1
    corr=alpha3/sc
    f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

    aux=f.readline()
    tim=[]
    vol=[]
    trovato = False
    for aux in f:
        tim.append(np.float64(aux.split('\t')[0]))
        vol.append(np.float64(aux.split('\t')[1]))
    tim = np.array(tim)
    vol = np.array(vol)

    tim_aux=np.linspace(tim.min(), tim.max(), 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0,d_in).subs(Psi, psi1)

    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = tim_aux / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]
    t_out.append(t_next)
################################################################################
    #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    #dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)
#################################################################################

    ada=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0, d_in).subs(Psi, psi1).subs(t, t_next)
    dep= Idep.subs(beta, bet).subs(IdA0, d_in).subs(t, t_next).subs(t0, t0_val)

    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)

    plt.figure()
    plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)


    #vrm=-0.7349960338647149
    #d_inc=dep_vec1_2[r,c]
    #a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep ****'+str(i)+' *****spike****')
        print(dep + d_inc)
        print('***** ada ****'+str(i)+' *****spike****')
        print(ada + a_inc)


        #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0, dep+d_inc).subs(Psi, psi1)

        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)
            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(Psi, psi1).subs(t, t_next)
            dep = Idep.subs(beta, bet).subs(IdA0, dep+d_inc).subs(t, t_next).subs(t0, t_init)
            #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)

            i=n_sp+1

    plt.plot(tim, vol)
    return t_out


def plot_tr_v3_vect(Vconvfact,vtm,vrm,a_in,d_in,a_inc,d_inc,bet,delta1,alpha3,sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    n_sp=np.size(a_inc)
    st_sign = 531.1
    corr=alpha3/sc
    f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

    aux=f.readline()
    tim=[]
    vol=[]
    trovato = False
    for aux in f:
        tim.append(np.float64(aux.split('\t')[0]))
        vol.append(np.float64(aux.split('\t')[1]))
    tim = np.array(tim)
    vol = np.array(vol)

    tim_aux=np.linspace(tim.min(), tim.max(), 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0,d_in).subs(Psi, psi1)

    print('***** dep init ****')
    print(d_in)
    print('***** ada init ****')
    print(a_in)

    lam_x = sym.lambdify(t, aux, modules=['numpy'])
    x_vals = tim_aux / 15.58
    y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
    t_next=x_vals[y_vals>vtm][0]
    t_out.append(t_next)
################################################################################
    #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0, d_in).subs(t,t_next)
    #dep = Idep.subs(gamma, gam).subs(t0, t0_val).subs(IdA0, d_in).subs(t,t_next)
#################################################################################

    ada=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0, d_in).subs(Psi, psi1).subs(t, t_next)
    dep= Idep.subs(beta, bet).subs(IdA0, d_in).subs(t, t_next).subs(t0, t0_val)


    #with open('spt_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    spt = pickle.load(f)

    #with open('Iadap_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + '_v0_' + str(vrm) + 'bet_' + str(bet) + 'gam_' + str(gam) + '_.pkl',"rb") as f:
    #    Iadap_mat = pickle.load(f)

    #with open('Idep_surf_IaA_' + str(IaA_max) + '_' + str(IaA_min) + '_IdA_' + str(IdA_max) + '_' + str(IdA_min) + '_bin_' + str(bin_x_unt) + 'gam_' + str(gam) + '_.pkl', "rb") as f:
    #    Idep_mat = pickle.load(f)

    #auxa = (IaA - a_in) ** 2
    #auxd = (IdA - d_in) ** 2
    #ada_nc=Iadap_mat[4,auxa.argmin(),auxd.argmin()]
    #dep_nc=Idep_mat[4, auxa.argmin(),auxd.argmin()]
    x_sector = np.logical_and(t0_val < x_vals, x_vals < t_next) * x_vals


    x_vals = x_sector[np.nonzero(x_sector)]
    y_vals = lam_x(x_vals)

    plt.figure()
    plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)


    #vrm=-0.7349960338647149
    #d_inc=dep_vec1_2[r,c]
    #a_inc=ada_vec1_2[r,c]
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('******dep ****'+str(i)+' *****spike****')
        print(dep )
        print('******dep inc ****' + str(i) + ' *****spike****')
        print(d_inc[i])
        print('******dep tot ****' + str(i) + ' *****spike****')
        print(dep + d_inc[i])

        #print('***** ada ****'+str(i)+' *****spike****')
        #print(ada)
        #print('***** ada inc ****' + str(i) + ' *****spike****')
        #print( a_inc[i])
        print('******ada tot ****' + str(i) + ' *****spike****')
        print(ada + a_inc[i])

        #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, ada+a_inc[i]).subs(IdA0, dep+d_inc[i]).subs(Psi, psi1)

        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)
            print('***** t_next****')
            print(t_next)

            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc[i]).subs(IdA0, dep+d_inc[i]).subs(Psi, psi1).subs(t, t_next)
            dep = Idep.subs(beta, bet).subs(IdA0, dep+d_inc[i]).subs(t, t_next).subs(t0, t_init)
            #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            #plt.figure()
            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)

            i=n_sp+1

    plt.plot(tim, vol)
    plt.title('trace (' + str(sc * alpha3) + 'nA)')
    return t_out


def plot_tr_v3_vect2(Vconvfact,vtm,vrm,a_inc,d_inc,bet,delta1,alpha3,sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    n_sp=np.size(a_inc)
    st_sign = 531.1
    corr=alpha3/sc
    f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

    aux=f.readline()
    tim=[]
    vol=[]
    trovato = False
    for aux in f:
        tim.append(np.float64(aux.split('\t')[0]))
        vol.append(np.float64(aux.split('\t')[1]))
    tim = np.array(tim)
    vol = np.array(vol)

    tim_aux=np.linspace(tim.min(), tim.max(), 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0,d_in).subs(Psi, psi1)
    plt.figure()
    t_next=0
    for i in range(n_sp):
        t_init=t_next+2/15.58
        print('******dep inc ****' + str(i) + ' *****spike****')
        print(d_inc[i])

        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada inc ****' + str(i) + ' *****spike****')
        print( a_inc[i])

        #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, a_inc[i]).subs(IdA0, d_inc[i]).subs(Psi, psi1)

        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58)>t_init]/ 15.58
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)
        t_init = t_next + 2 / 15.58

        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next)

#            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc[i]).subs(IdA0, dep+d_inc[i]).subs(Psi, psi1).subs(t, t_next)
#            dep = Idep.subs(beta, bet).subs(IdA0, dep+d_inc[i]).subs(t, t_next).subs(t0, t_init)
            #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/15.58
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)

            i=n_sp+1

    plt.plot(tim, vol)
    plt.title('trace ('+ str(sc*alpha3)+'nA)')
    return t_out



def plot_tr_v3_vect3_dep_ini(Vconvfact,vtm,vrm,a_inc,bet,delta1,alpha3,sc, tao_m, Idep_ini_vr, st_sign, dur_sign,Idep_ini):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    plotta=False
    n_sp=np.size(a_inc)
    corr=alpha3

    tim_aux=np.linspace(0, 1000, 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0,d_in).subs(Psi, psi1)
    if plotta:
        plt.figure()
    time=[]
    voltage=[]
    t_next=0
    for i in range(n_sp+1):
        if i>0:
            t_init=t_next+2/tao_m
        else:
            t_init=t_next


        #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, Idep_ini).subs(Psi, psi1)
        else:
            print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
            print('***** ada inc ****' + str(i) + ' *****spike****')
            print(a_inc[i-1])
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, a_inc[i-1]).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)

        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao_m)>t_init]/ tao_m
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next*tao_m)

#            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc[i]).subs(IdA0, dep+d_inc[i]).subs(Psi, psi1).subs(t, t_next)
#            dep = Idep.subs(beta, bet).subs(IdA0, dep+d_inc[i]).subs(t, t_next).subs(t0, t_init)
            #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            if plotta:
                plt.plot(st_sign + x_vals * tao_m, y_vals * Vconvfact)
            time.append(st_sign + x_vals * tao_m)
            voltage.append(y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao_m
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            if plotta:
                plt.plot(st_sign + x_vals * tao_m, y_vals * Vconvfact)
            time.append(st_sign + x_vals * tao_m)
            voltage.append(y_vals * Vconvfact)

            i=n_sp+1
        t_init = t_next + 2 / tao_m
    #plt.plot(tim, vol)
    #plt.title('trace ('+ str(sc*alpha3)+'nA)')
    return t_out,time,voltage

def plot_tr_v3_vect3_dep_ini2(Vconvfact,vtm,vrm,a_inc,bet,delta1,alpha3,sc, tao_m, Idep_ini_vr, st_sign, dur_sign,Idep_ini,Ith):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    plotta=False
    n_sp=np.size(a_inc)
    corr = round(alpha3 * sc) / 1000

    tim_aux=np.linspace(0, 1000, 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t0_val=0
    #vtm=np.mean(vol_tr)/66.35
    #vrm=-1.0
    #d_in=st_point_dep1[r]
    #a_in=st_point_ada1[r]
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t0_val).subs(V0, -1).subs(IaA0, a_in).subs(IdA0,d_in)
    #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t0_val).subs(V0, -1).subs(IaA0,a_in).subs(IdA0,d_in).subs(Psi, psi1)
    if plotta:
        plt.figure()
    time=[]
    voltage=[]
    t_next=0
    for i in range(n_sp+1):
        if i>0:
            t_init=t_next+2/tao_m
        else:
            t_init=t_next


        #aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, Idep_ini*(corr-Ith)*(corr>Ith)).subs(Psi, psi1)
        else:
            print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
            print('***** ada inc ****' + str(i) + ' *****spike****')
            print(a_inc[i-1])
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, a_inc[i-1]).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)

        #print('ada comp')
        #print( ada+a_inc)
        #print('dep comp')
        #print(dep + d_inc)

#        print('ada')
#        print(ada_nc)
#        print('dep')
#        print(dep_nc)
        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao_m)>t_init]/ tao_m
        y_vals = lam_x(x_vals)
    # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:



            t_next=x_vals[y_vals>vtm][0]
            t_out.append(t_next*tao_m)

#            ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc[i]).subs(IdA0, dep+d_inc[i]).subs(Psi, psi1).subs(t, t_next)
#            dep = Idep.subs(beta, bet).subs(IdA0, dep+d_inc[i]).subs(t, t_next).subs(t0, t_init)
            #ada = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0, t_init).subs(V0, vrm).subs(IaA0,ada+a_inc).subs(IdA0, dep+d_inc).subs(t, t_next)
            #dep = Idep.subs(gamma, gam).subs(t0, t_init).subs(IdA0, dep+d_inc).subs(t, t_next)
#        auxa = (IaA - a_in) ** 2
#        auxd = (IdA - d_in) ** 2
#        ada_nc = Iadap_mat[4, auxa.argmin(), auxd.argmin()]
#        dep_nc = Idep_mat[4, auxa.argmin(), auxd.argmin()]
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals


            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            if plotta:
                plt.plot(st_sign + x_vals * tao_m, y_vals * Vconvfact)
            time.append(st_sign + x_vals * tao_m)
            voltage.append(y_vals * Vconvfact)
        else:
            print('aoa')
            print(n_sp)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao_m
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            if plotta:
                plt.plot(st_sign + x_vals * tao_m, y_vals * Vconvfact)
            time.append(st_sign + x_vals * tao_m)
            voltage.append(y_vals * Vconvfact)

            i=n_sp+1
        t_init = t_next + 2 / tao_m
    #plt.plot(tim, vol)
    #plt.title('trace ('+ str(sc*alpha3)+'nA)')
    return t_out,time,voltage


def plot_tr_from_fit(Vconvfact, vtm, vrm, funzione,popt, bet, delta1,alpha3, sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    #n_sp=np.size(a_inc)
    st_sign = 531.1
    corr=alpha3/sc
    f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

    aux=f.readline()
    tim=[]
    vol=[]
    trovato = False
    for aux in f:
        tim.append(np.float64(aux.split('\t')[0]))
        vol.append(np.float64(aux.split('\t')[1]))
    tim = np.array(tim)
    vol = np.array(vol)

    tim_aux=np.linspace(tim.min(), tim.max(), 10000)

    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    plt.figure()
    i=0
    while t_init*15.58<500:

    #while i < 20:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada inc ****' + str(i) + ' *****spike****')
        print(funzione(t_init*15.58,*popt))
        print('t_init')
        print(t_init)
        print(t_init * 15.58)
        print('bet')
        print(bet)
        print('delta1')
        print(delta1)
        print('vrm')
        print(vrm)
        print('psi1')
        print(psi1)
        print(i)

        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,funzione(t_init * 15.58,*popt)).subs(IdA0, 0).subs(Psi, psi1)

        else:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, funzione(t_init*15.58,*popt[0])).subs(IdA0, 0).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58) > t_init] / 15.58
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            t_out.append(t_next)
            print("t_next .....")
            print(t_next)
            print(t_next*15.58)



            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next=500
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)

        i = i + 1
        print('t_next')
        print(t_next)
        t_init = t_next + 2 / 15.58

    plt.plot(tim, vol)
    plt.title('trace (' + str(alpha3/sc) + 'nA)')
    return t_out


def plot_tr_from_fit2(Vconvfact, vtm, vrm, exp_cum,popt,popt_dep, bet, delta1,alpha3, sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    #n_sp=np.size(a_inc)
    st_sign = 531.1
    corr=alpha3/sc
    try:
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

        aux=f.readline()
        tim=[]
        vol=[]
        trovato = False
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
            vol.append(np.float64(aux.split('\t')[1]))
        tim = np.array(tim)
        vol = np.array(vol)

        tim_aux=np.linspace(tim.min(), tim.max(), 10000)
    except:
        print('no trace')
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        aux = f.readline()
        tim = []
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
        tim = np.array(tim)

        tim_aux = np.linspace(tim.min(), tim.max(), 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    plt.figure()
    i=0
    #while i<20:
    while t_init * 15.58 < 500:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        print(exp_cum(t_init*15.58,popt[0],popt[1],popt[2]))
        print('***** dep  ****' + str(i) + ' *****spike****')
        print(exp_cum(t_init * 15.58, popt_dep[0], popt_dep[1],popt_dep[2]))

        popt_dep
        print('t_init')
        print(t_init)
        print(t_init * 15.58)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,exp_cum(t_init * 15.58,popt[0],popt[1],popt[2])).subs(IdA0, exp_cum(t_init * 15.58, popt_dep[0], popt_dep[1], popt_dep[2])).subs(Psi, psi1)

        else:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, exp_cum(t_init*15.58,popt[0],popt[1],popt[2])).subs(IdA0, exp_cum(t_init * 15.58, popt_dep[0], popt_dep[1], popt_dep[2])).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58) > t_init] / 15.58
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            t_out.append(t_next)
            print("t_next .....")
            print(t_next)
            print(t_next*15.58)


            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/15.58
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
            except:
                plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        t_init = t_next + 2 / 15.58
    try:
        plt.plot(tim, vol)
    except:
        print('no trace')
    plt.title('trace (' + str(alpha3/sc) + 'nA)')
    return t_out


def plot_tr_from_fit3(Vconvfact, vtm, vrm, exp_cum,popt, bet, delta1,alpha3, sc):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    #n_sp=np.size(a_inc)
    st_sign = 531.1
    corr=alpha3/sc
    try:
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

        aux=f.readline()
        tim=[]
        vol=[]
        trovato = False
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
            vol.append(np.float64(aux.split('\t')[1]))
        tim = np.array(tim)
        vol = np.array(vol)

        tim_aux=np.linspace(tim.min(), tim.max(), 10000)
    except:
        print('no trace')
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        aux = f.readline()
        tim = []
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
        tim = np.array(tim)

        tim_aux = np.linspace(tim.min(), tim.max(), 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    plt.figure()
    i=0
    #while i<20:
    while t_init * 15.58 < 500:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        print(exp_cum(t_init*15.58,popt[0],popt[1],popt[2]))

        print('t_init')
        print(t_init)
        print(t_init * 15.58)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,exp_cum(t_init * 15.58,popt[0],popt[1],popt[2])).subs(IdA0, 0).subs(Psi, psi1)

        else:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, exp_cum(t_init*15.58,popt[0],popt[1],popt[2])).subs(IdA0, 0).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / 15.58) > t_init] / 15.58
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            t_out.append(t_next)
            print("t_next .....")
            print(t_next)
            print(t_next*15.58)


            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/15.58
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact)
            except:
                plt.plot(st_sign + x_vals * 15.58, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        t_init = t_next + 2 / 15.58
    try:
        plt.plot(tim, vol)
    except:
        print('no trace')
    plt.title('trace (' + str(alpha3/sc) + 'nA)')
    return t_out


def plot_tr_from_fit4(Vconvfact, vtm, vrm, exp_cum,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    try:
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

        aux=f.readline()
        tim=[]
        vol=[]
        trovato = False
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
            vol.append(np.float64(aux.split('\t')[1]))
        tim = np.array(tim)
        vol = np.array(vol)

        tim_aux=np.linspace(tim.min(), tim.max(), 100000)
    except:
        print('no trace')
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        aux = f.readline()
        tim = []
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
        tim = np.array(tim)

        tim_aux = np.linspace(tim.min(), tim.max(), 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    plt.figure()
    i=0
    #while i<20:
    while t_init * tao < 400:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(exp_cum(t_init * tao, popt[0], popt[1], popt[2]))
        else:
            print(exp_cum(t_next*tao,popt[0],popt[1],popt[2]))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, 0).subs(Psi, psi1)

        else:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, exp_cum(t_next*tao,popt[0],popt[1],popt[2])).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if (t_next*tao<400):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)


            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        t_init = t_next + 2 / tao
    try:
        plt.plot(tim, vol)
    except:
        print('no trace')
    plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out

def plot_tr_from_fit5(Vconvfact, vtm, vrm, exp_cum,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,tempo_soglia):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3
    plotta=False
    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            cell_num = '95810010'
            abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            abf.setSweep(tr)
            print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            print(tr)
            vol_base = 0#-62.86698717948718
            vol = abf.sweepY + vol_base
            tim=abf.sweepX * 1000
            tim_aux = np.linspace(tim.min(), tim.max(), 100000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        aux = f.readline()
        tim = []
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
        tim = np.array(tim)

        tim_aux = np.linspace(tim.min(), tim.max(), 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    if plotta:
        plt.figure()
    i=0
    #while i<20:
    while t_init * tao < dur_sign:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(exp_cum(t_init * tao, popt[0], popt[1]))
        else:
            new_ada=exp_cum(t_next*tao,popt[0],popt[1])*(((t_next*tao)>tempo_soglia)*15+1)
            print(new_ada)
            print(exp_cum(t_next*tao,popt[0],popt[1]))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if t_init * tao== dur_sign:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,adaf).subs(IdA0, depf).subs(Psi, psi1)

        else:
            if i==0:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, 0).subs(Psi, psi1)

            else:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if i == 0:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, 0).subs(IdA0, 0).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, 0).subs(t0, t_init).subs(t, t_next)
            else:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0,vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t, t_next)
            print('.................adap....')
            print(adaf)
            print('..................dep....')
            print(depf)
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            if plotta:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            adaf=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,new_ada ).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t,dur_sign)
            depf=Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t,dur_sign)
            t_init=t_next
            vrm=y_vals[len(y_vals)-1]
            alpha3=0

        else:
            t_init = t_next + 2 / tao
    try:
        if plotta:
            plt.plot(tim, vol)
    except:
        print('no trace')
    if plotta:
        plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out

def plot_tr_from_fit6(Vconvfact, vtm, vrm, exp_cum,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,tempo_soglia):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3

    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            #cell_num = '95810010'
            #abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            #abf.setSweep(tr)
            #print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            #print(tr)
            #vol_base = 0#-62.86698717948718
            #vol = abf.sweepY + vol_base
            #tim=abf.sweepX * 1000
            tim_aux = np.linspace(0, 1000, 10000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        #f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        aux = f.readline()
        tim = []
        for aux in f:
            tim.append(np.float64(aux.split('\t')[0]))
        tim = np.array(tim)

        tim_aux = np.linspace(0, 1000, 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    plt.figure()
    i=0
    #while i<20:
    while t_init * tao < dur_sign:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(exp_cum(t_init * tao, popt[0], popt[1]))
        else:
            print(exp_cum(t_next*tao,popt[0],popt[1])*(((t_next*tao)>tempo_soglia)*15+1))
            print(exp_cum(t_next*tao,popt[0],popt[1]))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if i==0:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, 0).subs(Psi, psi1)

        else:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, exp_cum(t_next*tao,popt[0],popt[1])*(((t_next*tao)>tempo_soglia)*10+1)).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)

            plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = 500/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            t_init=t_next
        else:
            t_init = t_next + 2 / tao
    try:
        plt.plot(tim, vol)
    except:
        print('no trace')
    plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out

def plot_tr_from_fit_neuron_dep_ini(Vconvfact, vtm, vrm, funzione,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,tempo_soglia,dep_ini):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3
    plotta=False
    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            cell_num = '95810010'
            abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            abf.setSweep(tr)
            print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            print(tr)
            vol_base = 0#-62.86698717948718
            vol = abf.sweepY + vol_base
            tim=abf.sweepX * 1000
            tim_aux = np.linspace(tim.min(), tim.max(), 100000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        #f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        #aux = f.readline()
        tim = []
        #for aux in f:
         #   tim.append(np.float64(aux.split('\t')[0]))
        #tim = np.array(tim)

        tim_aux = np.linspace(0, 1000, 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    if plotta:
        plt.figure()
    i=0
    ada_vec=[]
    time=[]
    voltage=[]
    #while i<20:
    while t_init * tao < dur_sign:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(funzione(t_init * tao, *popt))
        else:
            new_ada=funzione(t_next*tao,*popt)*(((t_next*tao)>tempo_soglia)*15+1)
            ada_vec.append(new_ada)
            print(new_ada)
            print(funzione(t_next*tao,*popt))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if t_init * tao== dur_sign:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,adaf).subs(IdA0, depf).subs(Psi, psi1)

        else:
            if i==0:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, dep_ini).subs(Psi, psi1)

            else:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if i == 0:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, 0).subs(IdA0, dep_ini).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, dep_ini).subs(t0, t_init).subs(t, t_next)
            else:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0,vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t, t_next)
            print('.................adap....')
            print(adaf)
            print('..................dep....')
            print(depf)
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            time.append(st_sign + x_vals * tao)
            voltage.append(y_vals * Vconvfact)
            if plotta:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact)
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact*np.ones(x_vals.size))
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            adaf=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,new_ada ).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t,dur_sign)
            depf=Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t,dur_sign)
            t_init=t_next
            vrm=y_vals[len(y_vals)-1]
            alpha3=0

        else:
            t_init = t_next + 2 / tao
    try:
        if plotta:
            plt.plot(tim, vol)
    except:
        print('no trace')
    if plotta:
        plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out,ada_vec,time,voltage

def plot_tr_from_fit_neuron_dep_ini2(Vconvfact, vtm, vrm, funzione,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,tempo_soglia,dep_ini,Ith):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3
    plotta=False
    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            cell_num = '95810010'
            abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            abf.setSweep(tr)
            print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            print(tr)
            vol_base = 0#-62.86698717948718
            vol = abf.sweepY + vol_base
            tim=abf.sweepX * 1000
            tim_aux = np.linspace(tim.min(), tim.max(), 100000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        #f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        #aux = f.readline()
        tim = []
        #for aux in f:
         #   tim.append(np.float64(aux.split('\t')[0]))
        #tim = np.array(tim)

        tim_aux = np.linspace(0, 1000, 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    if plotta:
        plt.figure()
    i=0
    ada_vec=[]
    time=[]
    voltage=[]
    #while i<20:
    while t_init * tao < dur_sign:
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(funzione(t_init * tao, *popt))
        else:
            new_ada=funzione(t_next*tao,*popt)*(((t_next*tao)>tempo_soglia)*15+1)
            ada_vec.append(new_ada)
            print(new_ada)
            print(funzione(t_next*tao,*popt))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if t_init * tao== dur_sign:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,adaf).subs(IdA0, depf).subs(Psi, psi1)

        else:
            if i==0:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0, dep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(Psi, psi1)

            else:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if i == 0:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, 0).subs(IdA0, dep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, dep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(t0, t_init).subs(t, t_next)
            else:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0,vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, Idep_ini_vr).subs(t0, t_init).subs(t, t_next)
            print('.................adap....')
            print(adaf)
            print('..................dep....')
            print(depf)
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            time.append(st_sign + x_vals * tao)
            voltage.append(y_vals * Vconvfact)
            if plotta:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact)
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact*np.ones(x_vals.size))
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            adaf=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,new_ada ).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t,dur_sign)
            depf=Idep.subs(beta, bet).subs(IdA0, Idep_ini_vr).subs(t0, t_init).subs(t,dur_sign)
            t_init=t_next
            vrm=y_vals[len(y_vals)-1]
            alpha3=0

        else:
            t_init = t_next + 2 / tao
    try:
        if plotta:
            plt.plot(tim, vol)
    except:
        print('no trace')
    if plotta:
        plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out,ada_vec,time,voltage





def plot_tr_from_fit_neuron_rette_dep_ini(Vconvfact, vtm, vrm, funzione,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,lin_func_inf,vinc_inf,lin_func_sup,vinc_sup,Idep_ini):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3
    import sympy as sym
    I = sym.Symbol('I')
    plotta=False
    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    if corr*1000<vinc_inf:
        dur_sign=min(dur_sign,lin_func_inf.subs(I,corr*1000))

    if corr*1000>vinc_sup:
        dur_sign=min(dur_sign,lin_func_sup.subs(I,corr*1000))
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            cell_num = '95810010'
            abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            abf.setSweep(tr)
            print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            print(tr)
            vol_base = 0#-62.86698717948718
            vol = abf.sweepY + vol_base
            tim=abf.sweepX * 1000
            tim_aux = np.linspace(tim.min(), tim.max(), 100000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        #f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        #aux = f.readline()
        tim = []
        #for aux in f:
         #   tim.append(np.float64(aux.split('\t')[0]))
        #tim = np.array(tim)

        tim_aux = np.linspace(0, 1000, 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    if plotta:
        plt.figure()
    i=0
    ada_vec=[]
    time=[]
    voltage=[]
    #while i<20:
    while t_init * tao < dur_sign :
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(funzione(t_init * tao, *popt))
        else:
            new_ada=funzione(t_next*tao,*popt)
            ada_vec.append(new_ada)
            print(new_ada)
            print(funzione(t_next*tao,*popt))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if t_init * tao== dur_sign:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,adaf).subs(IdA0, depf).subs(Psi, psi1)

        else:
            if i==0:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0,Idep_ini).subs(Psi, psi1)

            else:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if i == 0:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, 0).subs(IdA0, Idep_ini).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, Idep_ini).subs(t0, t_init).subs(t, t_next)
            else:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0,vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t, t_next)
            print('.................adap....')
            print(adaf)
            print('..................dep....')
            print(depf)
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            time.append(st_sign + x_vals * tao)
            voltage.append(y_vals * Vconvfact)
            if plotta:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact)
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact*np.ones(x_vals.size))
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            adaf=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,new_ada ).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t,dur_sign)
            depf=Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t,dur_sign)
            t_init=t_next
            vrm=y_vals[len(y_vals)-1]
            alpha3=0

        else:
            t_init = t_next + 2 / tao
    try:
        if plotta:
            plt.plot(tim, vol)
    except:
        print('no trace')
    if plotta:
        plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out,ada_vec,time,voltage

def plot_tr_from_fit_neuron_rette_dep_ini2(Vconvfact, vtm, vrm, funzione,popt, bet, delta1,alpha3, sc,tao,Idep_ini_vr,st_sign,dur_sign,Neuron_tr,lin_func_inf,vinc_inf,lin_func_sup,vinc_sup,Idep_ini,Ith):

    import pyabf
    import matplotlib.pyplot as plt
    import sympy as sym
    import numpy as np
    from load_eq import load_v3
    import sympy as sym
    I = sym.Symbol('I')
    plotta=False
    #n_sp=np.size(a_inc)

    corr=round(alpha3*sc)/1000
    if corr*1000<vinc_inf:
        dur_sign=min(dur_sign,lin_func_inf.subs(I,corr*1000))

    if corr*1000>vinc_sup:
        dur_sign=min(dur_sign,lin_func_sup.subs(I,corr*1000))
    try:
        if Neuron_tr:
            f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_"+str(corr)[0:3]+".soma.v.txt", "r")
    # f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\CA1_pyr\\r_seed2_0\\step_" + str(corr)[0:3] + ".soma.v.txt","r")

            aux=f.readline()
            tim=[]
            vol=[]
            trovato = False
            for aux in f:
                tim.append(np.float64(aux.split('\t')[0]))
                vol.append(np.float64(aux.split('\t')[1]))
            tim = np.array(tim)
            vol = np.array(vol)

            tim_aux=np.linspace(tim.min(), tim.max(), 100000)
        else:
            tr = int(round(alpha3*sc)/200) + 10
            cell_num = '95810010'
            abf = pyabf.ABF('C:\\Users\\INA\\Downloads\\PC-cAC\\PC-cAC\\tutti\\' + cell_num + '.abf')
            abf.setSweep(tr)
            print('trrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr')
            print(tr)
            vol_base = 0#-62.86698717948718
            vol = abf.sweepY + vol_base
            tim=abf.sweepX * 1000
            tim_aux = np.linspace(tim.min(), tim.max(), 100000)
        #plt.plot(, abf.sweepY)

    except:
        print('no trace')
        #f = open("C:\\Users\\INA\\PycharmProjects\\neuron_model\\CA1_pyr2\\r_seed2_0\\step_1.0.soma.v.txt","r")
        #aux = f.readline()
        tim = []
        #for aux in f:
         #   tim.append(np.float64(aux.split('\t')[0]))
        #tim = np.array(tim)

        tim_aux = np.linspace(0, 1000, 10000)
    #alpha3 = [0.252587 * 4, 0.252587 * 5]
    t_out = []
    tr=int(alpha3/sc)+10
    #cell_num = '95810005'

    t = sym.Symbol('t')
    delta, Psi, alpha, beta, gamma, IaA0, IdA0, t0, V0 = sym.symbols('delta,Psi,alpha,beta,gamma,IaA0,IdA0,t0,V0')

    [V, Iadap, Idep] = load_v3()
    psi1 = ((-4) * bet + ((1 + delta1) ** 2)) ** (0.5)

    t_init=0
    if plotta:
        plt.figure()
    i=0
    ada_vec=[]
    time=[]
    voltage=[]
    #while i<20:
    while t_init * tao < dur_sign :
        print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
        print('***** ada  ****' + str(i) + ' *****spike****')
        if i==0:
            print(funzione(t_init * tao, *popt))
        else:
            new_ada=funzione(t_next*tao,*popt)
            ada_vec.append(new_ada)
            print(new_ada)
            print(funzione(t_next*tao,*popt))

        print('t_init')
        print(t_init)
        print(t_init * tao)


        # aux = V.subs(alpha, alpha3).subs(beta, bet).subs(gamma, gam).subs(t0,t_init ).subs(V0, vrm).subs(IaA0, ada+a_inc).subs(IdA0,dep+d_inc)
        if t_init * tao== dur_sign:
            aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,adaf).subs(IdA0, depf).subs(Psi, psi1)

        else:
            if i==0:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, -1).subs(IaA0,0).subs(IdA0,Idep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(Psi, psi1)

            else:
                aux = V.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1)


        lam_x = sym.lambdify(t, aux, modules=['numpy'])
        x_vals = tim_aux[(tim_aux / tao) > t_init] / tao
        y_vals = lam_x(x_vals)
        # plt.plot(31.1 + x_vals * 15.58, y_vals * 66.35)


        if len(np.nonzero(y_vals > vtm)[0]) > 0:

            t_next = x_vals[y_vals > vtm][0]#+t_init
            if i == 0:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0, 0).subs(IdA0, Idep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, Idep_ini*(corr*1000-Ith)*(corr*1000>Ith)).subs(t0, t_init).subs(t, t_next)
            else:
                adaf = Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0,vrm).subs(IaA0, new_ada).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t, t_next)
                depf = Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t, t_next)
            print('.................adap....')
            print(adaf)
            print('..................dep....')
            print(depf)
            if (t_next*tao<dur_sign):
                t_out.append(t_next*tao)
                print("t_next .....")
                print(t_next)
                print(t_next*tao)
            else:
                t_next=dur_sign/tao

            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals

            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            time.append(st_sign + x_vals * tao)
            voltage.append(y_vals * Vconvfact)
            if plotta:
                plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            #st_sign=st_sign#+t_next*15.58
            print("st_sign .....")
            print(st_sign)
        else:
            print('aoa')
            print(i)
            print(len(np.nonzero(y_vals > vtm)[0]))
            print(vtm)
            t_next = dur_sign/tao
            x_sector = np.logical_and(t_init < x_vals, x_vals < t_next) * x_vals
            x_vals = x_sector[np.nonzero(x_sector)]
            y_vals = lam_x(x_vals)
            try:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact)
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact)
            except:
                time.append(st_sign + x_vals * tao)
                voltage.append(y_vals * Vconvfact*np.ones(x_vals.size))
                if plotta:
                    plt.plot(st_sign + x_vals * tao, y_vals * Vconvfact*np.ones(x_vals.size))

        i = i + 1
        print('t_next')
        print(t_next)
        if t_next==dur_sign:
            adaf=Iadap.subs(alpha, alpha3).subs(beta, bet).subs(delta, delta1).subs(t0, t_init).subs(V0, vrm).subs(IaA0,new_ada ).subs(IdA0, Idep_ini_vr).subs(Psi, psi1).subs(t,dur_sign)
            depf=Idep.subs(beta, bet).subs(IdA0, new_ada).subs(t0, t_init).subs(t,dur_sign)
            t_init=t_next
            vrm=y_vals[len(y_vals)-1]
            alpha3=0

        else:
            t_init = t_next + 2 / tao
    try:
        if plotta:
            plt.plot(tim, vol)
    except:
        print('no trace')
    if plotta:
        plt.title('trace (' + str(alpha3*sc/1000) + 'nA)')
    return t_out,ada_vec,time,voltage