# -*- coding: utf-8 -*-
"""
Created on Tue Jun 21 13:03:22 2016

@author: Benjamin
"""
from __future__ import division
import numpy as np
import scipy.integrate as integrate
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import math

"""
The first part of the code is the same as the one to compute Fig2.
More information can be found in the dedicated file.
"""

# Geometrical functions

def f(theta):
    return (2*np.pi/3.)*(1-np.cos(theta))**2*(2+np.cos(theta))/((np.sin(theta)
    )**3)

def g(theta):
    return 4*np.pi*(1.-np.cos(theta))/((np.sin(theta))**2)

def f_prime(theta):
    return 2*np.pi/(1.+np.cos(theta))**2
    
def g_prime(theta):
    return 4*np.pi*np.sin(theta)/(1.+np.cos(theta))**2

    
# Experimental data

time=np.genfromtxt('../Data_DIBs/bilayer_area_data.txt')[:,0]
time_star=np.genfromtxt('../Data_DIBs/a_star.txt')[:,0]
area_star=np.genfromtxt('../Data_DIBs/a_star.txt')[:,1]
bilayer=np.genfromtxt('../Data_DIBs/bilayer_area_data.txt')[:,1]
monolayer_left=np.genfromtxt('../Data_DIBs/left_monolayer_area_data.txt')[:,1]
monolayer_right=np.genfromtxt('../Data_DIBs/right_monolayer_area_data.txt'
)[:,1]

a_exp=np.sqrt(bilayer/np.pi)
theta_left_exp=np.arccos(4*np.pi/(2*monolayer_left/a_exp**2)-1)
theta_right_exp=np.arccos(4*np.pi/(2*monolayer_right/a_exp**2)-1)

# Fixed parameters

V = 0.18e-4*(0.1e-3*np.pi*0.1e-3) # Total volume of the setup [m^3]
N=6.8e11 # Total number of lipids in the setup
N_a=6.02e23 # Avogadro number [mol^(-1)]
T=298 # Temperature [K]
V_w=18e-6 # Molar volume of water [m^3.mol^(-1)]
k_B=1.38e-23 # Boltzmann constant [J.K^(-1)]

class Simulation(object):
    
    def __init__(self, exc_s, gam_0, P_f, xi, k_on, k_off_m, k_off_b,
                 C_out,C_in,mu_b):
        
        """
        Create a framework for the simulation with free parameters.
        
        Arguments:
        - exc_s: total number of available site on each monolayer per unit 
        surface [/m^2]
        - gam_0: surface tension of o/w interface for a clean interface [N/m]
        - P_f: permeation coefficient of water in DOPC monolayer [m/s]
        - xi: diffusion constant of lipids from the monolayer to the bilayer
        [/m/s]
        - k_on: adsorption constant of lipids [/s]
        - k_off_m: monolayer desorption constant of lipids [/s]
        - k_off_b: bilayer desorption constant of lipids [/s]
        - C_out: outside osmolarity [mol/m³]
        - C_in: inside osmolarity [mol/m³]
        
        Raises:
        - TypeError if parameters are not floats 
        """
        
        try:
            self.k_off_m=float(k_off_m)
            self.k_off_b=float(k_off_b)
            self.k_on=float(k_on)
            self.exc_s=float(exc_s)
            self.gam_0=float(gam_0)
            self.P_f=float(P_f)
            self.xi=float(xi)
            self.C_out=float(C_out)
            self.C_in=float(C_in)
            self.mu_b=float(mu_b)
        except (TypeError, ValueError,AssertionError):
            raise TypeError ('''Simulation parameters must be floats, except m 
            which needs to be an integer.''')
        
    def values_arrays_init(self):
        
        """
        Create arrays for each variable:
        - a: radius of the bilayer [m]
        - theta_l: left droplet co-latitude [rad]
        - theta_r: right droplet co-latitude [rad]
        - R_l: left droplet radius [m]
        - R_r: right droplet radius [m]
        - S_l: left droplet area [m^2]
        - S_r: right droplet area [m^2]
        - A: bilayer area [m^2]
        - V_l: left droplet volume [m^3]
        - V_r: right droplet volume [m^3]
        - exc_l: left monolayer surface excess [/m^2]
        - exc_r: right monolayer surface excess [/m^2]
        - exc_b: bilayer surface excess [/m^2]
        """
        
        self.a=np.array([])
        self.theta_l=np.array([])
        self.theta_r=np.array([])
        self.R_l=np.array([])
        self.R_r=np.array([])
        self.S_l=np.array([])
        self.S_r=np.array([])
        self.A=np.array([])
        self.exc_l=np.array([])
        self.exc_r=np.array([])
        self.exc_b=np.array([])
        
    def initialization(self,exc_l0,exc_r0,exc_b0):
        
        """
        Create initial conditions for the simulation.
        
        Arguments:
        - exc_l0: initial surface excess of left monolayer [/m^2]
        - exc_r0: initial surface excess of right monolayer [/m^2]
        - exc_b0: initial surface excess of the bilayer [/m^2]
        
        Return: 
        - Array containing initial conditions:
            - a0: initial bilayer radius [m]
            - theta_l0: initial left droplet polar angle [rad]
            - theta_r0: initial right droplet polar angle [rad]
            - exc_l0: initial left monolayer surface excess [/m^2]
            - exc_r0: initial right monolayer surface excess [/m^2]
            - exc_b0: initial bilayer surface excess [/m^2]
        
        Raise:
        - TypeError if initial surface excesses are not floats
        - ValueError if the programm cannot compute the initial value of 
        lipids density because simulation parameters are not sensible or
        because surface excesses are negative
        """
        
        try:
            exc_l0=float(exc_l0)
            exc_r0=float(exc_r0)
            exc_b0=float(exc_b0)
        except (TypeError,ValueError,AssertionError):
            raise TypeError ('Initial surface excesses must be floats.')
        if (exc_l0<0 or exc_r0<0 or exc_b0<0):
            raise ValueError ('Surface excesses cannot be negative.')
        a0= a_exp[0]
        theta_l0=theta_left_exp[0]
        theta_r0=theta_right_exp[0]
        self.V_l0=0.5*a0**3*f(theta_l0)
        self.V_r0=0.5*a0**3*f(theta_r0)
        return np.array([a0,theta_l0,theta_r0,exc_l0,exc_r0,exc_b0])
        
    def evolution_volume(self,a,theta,exc,side):
        
        """
        Compute the time derivative of droplet volume.
        
        Arguments:
        - a: bilayer radius [m]
        - theta: droplet polar angle [rad]
        - exc: monolayer surface excess [/m^2]
        - side: 'l' for the left droplet and 'r' for the right droplet
    
        Return:
        - Time derivative of droplet volume
    
        Raises:
        - TypeError if side is not a string
        - ValueError if side is not 'l' or 'r' or if at least one surface
        excess is larger than maximal surface excess
        """
        
        try:
            side=str(side)
        except (TypeError,ValueError,AssertionError):
            raise TypeError('side must be a string.')
        if side=='l':
            return -self.P_f*0.5*a**2*g(theta)*(
            self.C_out*V_w-self.C_in*V_w*self.V_l0/(0.5*a**3*f(theta)
            )-2*V_w*np.sin(theta)*(self.gam_0+k_B*T*self.exc_s*np.log(
            1-exc/self.exc_s))/(a*N_a*k_B*T))
        if side=='r':
            return -self.P_f*0.5*a**2*g(theta)*(
            self.C_out*V_w-self.C_in*V_w*self.V_r0/(0.5*a**3*f(theta)
            )-2*V_w*np.sin(theta)*(self.gam_0+k_B*T*self.exc_s*np.log(
            1-exc/self.exc_s))/(a*N_a*k_B*T))
        else:
            raise ValueError('''side must be 'l' or 'r'.''')
        
    def evolution_a(self,a,exc_l,exc_r,exc_b):
        
        """
        Compute the time derivative of bilayer radius.
        
        Arguments:
        - a: bilayer radius [m]
        - exc_l: left monolayer surface excess [/m^2]
        - exc_r: right monolayer surface excess [/m^2]
        - exc_b: bilayer surface excess [/m^2]
    
        Return:
        - Time derivative of bilayer radius
        
        Raise:
        - ValueError if at least one surface excess is negative or if one 
        surface excess is larger than maximum surface excess
        """
        if (math.isnan(np.log(exc_b)) or math.isnan(np.log(exc_b))):
            raise ValueError('Bilayer surface excess must be positive.')
        if (math.isinf(np.log(exc_l)) or math.isinf(np.log(exc_l))):
            raise ValueError('Left monolayer surface excess must be positive.')
        if (math.isinf(np.log(exc_r)) or math.isinf(np.log(exc_r))):
            raise ValueError('Right monolayer surface excess must be positive.'
            )
        if (math.isnan(np.log(self.exc_s-exc_b)) or math.isnan(np.log(
        self.exc_s-exc_b))):
            raise ValueError('''Bilayer surface excess must be lower than 
            maximum surface excess.''')
        if (math.isinf(np.log(self.exc_s-exc_l)) or math.isinf(np.log(
        self.exc_s-exc_l))):
            raise ValueError('''Left monolayer surface excess must be must be 
            lower than maximum surface excess.''')
        if (math.isinf(np.log(self.exc_s-exc_r)) or math.isinf(np.log(
        self.exc_s-exc_r))):
            raise ValueError('''Left monolayer surface excess must be must be 
            lower than maximum surface excess.''')
        else:
            return self.xi*(np.log(exc_l*exc_r/((self.exc_s-exc_l)*(
            self.exc_s-exc_r)))-2*self.mu_b)/(2*exc_b
            )-self.k_off_b*a/2+self.k_on*(self.exc_s-exc_b)*a/(2*exc_b)
        
    def evolution_theta(self,a,theta_l,theta_r,exc_l,exc_r,exc_b,side):
        
        """
        Compute the time derivative of droplet polar angle.
        
        Arguments:
        - a: bilayer radius [m]
        - theta_l: left droplet polar angle [rad]
        - theta_r: right droplet polar angle [rad]
        - exc_l: left monolayer surface excess [/m^2]
        - exc_r: right monolayer surface excess [/m^2]
        - side: 'l' for the left droplet and 'r' for the right droplet
    
        Return:
        - Time derivative of droplet polar angle
    
        Raises:
        - TypeError if side is not a string
        - ValueError if side is not 'l' or 'r', if at least one surface 
        excess is negative or if at least one surface excess is larger than
        maximum surface excess
        """ 
        try:
            side=str(side)
        except (TypeError,ValueError,AssertionError):
            raise TypeError('side must be a string.')
        if side=='l':
            return (1/(a**3*f_prime(theta_l)))*(2*self.evolution_volume(
            a,theta_l,exc_l,side)-3*a**2*self.evolution_a(a,exc_l,exc_r,exc_b)*f(theta_l))
        if side=='r':
            return (1/(a**3*f_prime(theta_r)))*(2*self.evolution_volume(
            a,theta_r,exc_r,side)-3*a**2*self.evolution_a(a,exc_l,exc_r,exc_b)*f(theta_r))
        else:
            raise ValueError('''side must be 'l' or 'r'.''')
                
    def evolution_exc_m_diffusion(self,a,theta_l,theta_r,exc_l,exc_r,exc_b,side
    ):
         
        """
        Compute the time derivative of monolayer surface excess due to 
        adsorption and desorption of lipids from oil (mainly reverse micelles).
        
        Arguments:
        - a: bilayer radius [m]
        - theta_l: left droplet polar angle [rad]
        - theta_r: right droplet polar angle [rad]
        - exc_l: left monolayer surface excess [/m^2]
        - exc_r: right monolayer surface excess [/m^2]
        - side: 'l' for the left droplet and 'r' for the right droplet
    
        Return:
        - Time derivative of monolayer surface excess due to adsorption and 
        desorption of free lipids in oil
    
        Raises:
        - TypeError if side is not a string
        - ValueError if side is not 'l' or 'r', if at least one surface 
        excess is negative or if at least one surface excess is larger than
        maximum surface excess.
        """
        
        try:
            side=str(side)
        except (TypeError,ValueError,AssertionError):
            raise TypeError('side must be a string.')
        if side=='l':
            return self.k_on*(self.exc_s-exc_l)-self.k_off_m*exc_l-exc_l*(
            2*self.evolution_a(a,exc_l,exc_r,exc_b)/a+(g_prime(theta_l)/g(
            theta_l))*self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,exc_b,side))
        if side=='r':
            return self.k_on*(self.exc_s-exc_r)-self.k_off_m*exc_r-exc_r*(
            2*self.evolution_a(a,exc_l,exc_r,exc_b)/a+(g_prime(theta_r)/g(
            theta_r))*self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,exc_b,side))
        else:
            raise ValueError('''side must be 'l' or 'r'.''')
        
    def evolution_exc_m_flow_bilayer(self,a,theta,exc_m,exc_b):
        
        """
        Compute the time derivative of monolayer surface excess due to lipid 
        flux from the monolayer to the bilayer.
        
        Arguments:
        - a: bilayer radius [m]
        - theta: droplet polar angle [rad]
        - exc_m: monolayer surface excess [/m^2]
        - exc_b: bilayer surface excess [/m^2]
    
        Return:
        - Time derivative of monolayer surface excess due to lipid flux from 
        the monolayer to the bilayer
        
        Raise:
        - ValueError if at least one surface excess is negative or if at least
        one surface excess is larger than maximum surface excess
        """
        
        if (math.isnan(np.log(exc_b)) or math.isnan(np.log(exc_b))):
            raise ValueError('Bilayer surface excess must be positive.')
        if (math.isinf(np.log(exc_m)) or math.isinf(np.log(exc_m))):
            raise ValueError('Monolayer surface excess must be positive.')
        if (math.isnan(np.log(self.exc_s-exc_b)) or math.isnan(np.log(
        self.exc_s-exc_b))):
            raise ValueError('''Bilayer surface excess must be lower than 
            maximum surface excess.''')
        if (math.isnan(np.log(self.exc_s-exc_m)) or math.isnan(np.log(
        self.exc_s-exc_m))):
            raise ValueError('''Bilayer surface excess must be lower than 
            maximum surface excess.''')
        return 4*np.pi*self.xi*(self.mu_b-np.log(exc_m/(self.exc_s-exc_m)))/(a*g(theta))
           
    def evolution(self,t_f,dt,y0):
            
        """
        Compute the dynamics of the system from t=0 to t_f using odeint 
        function of scipy.integrate (LSODA Fortran algorithm) with time step 
        dt.
            
        Arguments:
        - t_f: computation duration [s]
        - dt: simulation time step [s]
        - y0: array of initial conditions 
        
        Raises:
        - TypeError if parameters are not float
        - ValueError if parameters are not positive or if dt is larger than t_f
        - ValueError if at least one surface excess becomes negative
        """
            
        try:
            t_f=float(t_f)
            self.dt=float(dt)
        except (ValueError,TypeError,AssertionError):
            raise TypeError('Simulation time parameters must be floats.')
        if t_f<0 or dt<0 or dt==0.:
            raise ValueError('Simulation time parameters must be positive.')
        if t_f<dt:
            raise ValueError('''Simulation duration must be larger than 
            simulation time step.''')
        def lindiff(y,t):
            a, theta_l, theta_r,exc_l, exc_r, exc_b = y # unpack current values of y
            derivs = [self.evolution_a(a,exc_l,exc_r,exc_b
            ),self.evolution_theta(a,theta_l,theta_r,exc_l,exc_r,exc_b,'l'
            ),self.evolution_theta(a,theta_l,theta_r,exc_l,exc_r,exc_b,'r'
            ),self.evolution_exc_m_diffusion(
            a,theta_l,theta_r,exc_l,exc_r,exc_b,'l'
            )+self.evolution_exc_m_flow_bilayer(a,theta_l,exc_l,exc_b
            ),self.evolution_exc_m_diffusion(
            a,theta_l,theta_r,exc_l,exc_r,exc_b,'r'
            )+self.evolution_exc_m_flow_bilayer(a,theta_r,exc_r,exc_b),0]
            return derivs  
            
        # Make time array for solution
        self.t = np.arange(0., t_f, dt)
        
        # Resolution and storage
        psoln=integrate.odeint(lindiff, y0,self.t)
        self.a=psoln[:,0]
        self.theta_l=psoln[:,1]
        self.theta_r=psoln[:,2]
        self.exc_l=psoln[:,3]
        self.exc_r=psoln[:,4]
        self.exc_b=psoln[:,5]
        self.R_l=self.a/np.sin(self.theta_l)
        self.R_r=self.a/np.sin(self.theta_r)
        self.S_l=self.a**2*0.5*g(self.theta_l)
        self.S_r=0.5*self.a**2*g(self.theta_r)
        self.A=np.pi*self.a**2

    def plot_solutions(self):
            
        """
        Plot simulation results.
        
        Returns:
        - Plot of bilayer radius a [micrometers] as a function of time [s]
        - Plot of droplets polar angles [deg] as a function of time [s]
        - Plot of areas [square micrometers] and area ratio as a function of 
        time [s]
        - Plot surface excesses [in units of exc_s] as a function of time [s]
        - Plot lipid numbers as a function of time [s]
        - Plot alpha coefficients [m/s] as a function of time [s]
        - Plot surface tensions [mN/m] as a function of time [s]
        - Plot volumes [fL] as a function of time [s]
        - Plot surface tensions [mN/m] as a function of surface excesses [in 
        units of exc_s]
        """
            
        fig1,ax1=plt.subplots()
        ax1.plot(self.t,self.a*1e6,linestyle='--',linewidth=5,color='darkred',
                 label=r'$a_{th}$')
        ax1.set_xlabel('Evaporation time (s)',fontsize=40)
        ax1.set_ylabel(r'Bilayer radius ($\mu m$)',fontsize=40)
        ax1.legend(loc=0,prop={'size':50})
        ax1.tick_params(labelsize=30)
        fig1.show()
            
        fig2,ax2=plt.subplots()
        ax2.plot(self.t,self.theta_l*180/np.pi,color='forestgreen',
                 linestyle='-.',linewidth=5,label=r'$\theta_{l}^{th}$')
        ax2.plot(self.t,self.theta_r*180/np.pi,linestyle=':',
                 color='darkblue',linewidth=5,label=r'$\theta_{r}^{th}$')
        ax2.set_xlabel('Evaporation time (s)',fontsize=40)
        ax2.set_ylabel(r'Angle ($^{\circ}$)',fontsize=40)
        ax2.tick_params(labelsize=30)     
        ax2.legend(loc=0,prop={'size':35})
        fig2.show()
            
        fig3,ax3=plt.subplots()
        ax4=ax3.twinx()
        fig3.set_size_inches(15,22,forward=True)
        ax3.plot(self.t,self.S_l*1e12,color='forestgreen',linestyle='-.',
                 linewidth=3,label='Left monolayer')
        ax3.plot(self.t,self.S_r*1e12,linestyle=':',color='darkblue',
                 linewidth=3,label='Right monolayer')
        ax3.plot(self.t,self.A*1e12,linestyle='dashed',color='darkred',
                 linewidth=3,label='Bilayer')
        ax4.plot(self.t,2*np.pi/(0.5*g(self.theta_l)+0.5*g(self.theta_r)),
                 color='slategray',linestyle='-',linewidth=2,
                 label=r'$A_{th}^{*}$')
        ax3.set_xlabel('Evaporation time (s)',fontsize=23)
        ax3.set_ylabel(r'Surface area ($\mu m^{2}$)',fontsize=25)
        ax4.set_ylabel(r'$A^{*}=\frac{2A}{S_{l}+S_{r}}$',fontsize=35)
        ax4.tick_params(labelsize=30)
        ax3.tick_params(labelsize=30)
        ax4.legend(loc=1,prop={'size':27})
        ax3.legend(loc=6,prop={'size':18})
        fig3.show()
                    
        fig5,ax5=plt.subplots()
        fig5.set_size_inches(15,22,forward=True)
        ax5.plot(self.t,self.exc_l/self.exc_s,color='forestgreen',
                 linestyle='-.',linewidth=5,label=r'$\Gamma_{m}^{l}$')
        ax5.plot(self.t,self.exc_r/self.exc_s,linestyle=':',
                 color='darkblue',linewidth=5,label=r'$\Gamma_{m}^{r}$')
        ax5.plot(self.t,self.exc_b/self.exc_s,linewidth=5,linestyle='dashed',
                 color='darkred',label=r'$\Gamma_{b}$')
        ax5.set_xlabel('Evaporation time (s)',fontsize=40)
        ax5.set_ylabel(r'Surface excess ($\Gamma_{\infty}$)',fontsize=40)
        ax5.legend(loc=0,prop={'size':43})
        ax5.tick_params(labelsize=30)
        fig5.show()
        
        fig6,ax6=plt.subplots()
        fig6.set_size_inches(15,22,forward=True)
        ax6.plot(self.t,self.exc_l*self.S_l*1e-8,color='forestgreen',
                 linestyle='-.',linewidth=5,label=r'$N_{m}^{l}$')
        ax6.plot(self.t,self.exc_r*self.S_r*1e-8,linestyle=':',
                 color='darkblue',linewidth=5,label=r'$N_{m}^{r}$')
        ax6.plot(self.t,self.exc_b*2*self.A*1e-8,linestyle='dashed',
                 color='darkred',linewidth=5,label=r'$N_{b}$')
        ax6.set_xlabel('Evaporation time (s)',fontsize=40)
        ax6.set_ylabel(r'Lipid number $\left(\times\ 10^{8}\right)$',
                       fontsize=40)
        ax6.tick_params(labelsize=30)
        ax6.legend(loc=0,prop={'size':35})
        fig6.show()
        
        fig8,ax8=plt.subplots()
        ax8.plot(self.t,1e3*(self.gam_0+k_B*T*self.exc_s*np.log(
        1-self.exc_l/self.exc_s)),color='forestgreen',linestyle='-.',
        linewidth=5,label=r'$\gamma_{m}^{l}$')
        ax8.plot(self.t,1e3*(self.gam_0+k_B*T*self.exc_s*np.log(
        1-self.exc_r/self.exc_s)),linestyle=':',color='darkblue',linewidth=5,
        label=r'$\gamma_{m}^{r}$')
        gam_b0=-np.cos(theta_left_exp[0])*(self.gam_0+k_B*T*self.exc_s*np.log(
        1-self.exc_l[0]/self.exc_s))-np.cos(theta_right_exp[0])*(
        self.gam_0+k_B*T*self.exc_s*np.log(1-self.exc_r[0]/self.exc_s))
        ax8.plot(self.t,1e3*(gam_b0+k_B*T*self.exc_s*np.log((
        1-self.exc_b/self.exc_s)/(1-self.exc_b[0]/self.exc_s))),
        linestyle='dashed',color='darkred',linewidth=5,label=r'$\gamma_{b}$')
        ax8.plot(self.t,((self.gam_0+k_B*T*self.exc_s*np.log(
        1-self.exc_l/self.exc_s))*np.cos(self.theta_l)+(
        self.gam_0+k_B*T*self.exc_s*np.log(1-self.exc_r/self.exc_s))*np.cos(
        self.theta_r)+gam_b0+k_B*T*self.exc_s*np.log((1-self.exc_b/self.exc_s
        )/(1-self.exc_b[0]/self.exc_s)))*1e3,linestyle='-',linewidth=2,
        color='PaleVioletRed',
        label=r'$\gamma_{b}+\gamma_{m}^{l}\cos\theta_{l}+\gamma_{m}^{r}\cos\theta_{r}$')
        ax8.set_xlabel('Evaporation time (s)',fontsize=40)
        ax8.set_ylabel('Surface tension (mN/m)',fontsize=40)
        ax8.tick_params(labelsize=30)
        ax8.legend(loc=0,prop={'size':35})
        fig8.show()
        
        fig9,ax9=plt.subplots()
        ax9.plot(self.t,1e18*0.5*self.a**3*f(self.theta_l),
                 color='forestgreen',linestyle='-.',linewidth=5,
                 label=r'$V_{l}$')
        ax9.plot(self.t,1e18*0.5*self.a**3*f(self.theta_r),linestyle=':',
                 color='darkblue',linewidth=5,label=r'$V_{r}$')
        ax9.set_xlabel('Evaporation time (s)',fontsize=40)
        ax9.set_ylabel('Droplet volume (fL)',fontsize=40)
        ax9.tick_params(labelsize=30)
        ax9.legend(loc=0,prop={'size':43})
        fig9.show()
        
    def fixed_points(self,exc_b0,guess):
        
        """
        Determine the fixed points of the differential system.
        
        Argument:
        - exc_b0: initial bilayer surface excess [/m^2]
        - guess: first guess for fixed point
        
        Returns:
        - Array containing fixed point values for:
            - a: bilayer radius [m]
            - theta_l: left droplet polar angle [rad]
            - theta_r: right droplet polar angle [rad]
            - exc_l: left monolayer surface excess [/m^2]
            - exc_r: right monolayer surface excess [/m^2]
            - exc_b: bilayer surface excess [/m^2]
        - 'DIB is unstable.' if there is not any fixed point
        """
        
        def lindiff(y,t):
            a, theta_l, theta_r,exc_l, exc_r, = y # unpack current values of y
            derivs = [self.evolution_a(a,exc_l,exc_r,exc_b0
            ),self.evolution_theta(a,theta_l,theta_r,exc_l,exc_r,exc_b0,'l'
            ),self.evolution_theta(a,theta_l,theta_r,exc_l,exc_r,exc_b0,'r'
            ),self.evolution_exc_m_diffusion(
            a,theta_l,theta_r,exc_l,exc_r,exc_b0,'l'
            )+self.evolution_exc_m_flow_bilayer(a,theta_l,exc_l,exc_b0
            ),self.evolution_exc_m_diffusion(
            a,theta_l,theta_r,exc_l,exc_r,exc_b0,'r'
            )+self.evolution_exc_m_flow_bilayer(a,theta_r,exc_r,exc_b0)]
            return derivs
        try:
            return optimize.anderson(lindiff,guess)
        except Exception:
            return 'DIB is unstable.'

"""
Values of parameters that can be modified like in the code for Fig2.
"""
# Fixed parameters
exc_s=
gam_0=
k_on=
k_off_b=
exc_l0=
exc_r0=
exc_b0=
C_out=
C_in=
P_f=
k_off_m=
xi=
mu_b=
# Radius of the droplet
R=2*(gam_0+k_B*T*exc_s*np.log(1-exc_b0/exc_s))/(N_a*k_B*T*(C_out-C_in))
# K is the array containing the values of desorption constants (can be modified)
K=np.concatenate((np.mgrid[1e-5:1e-4:1e-5],np.mgrid[1e-4:1e-3:1e-4],np.mgrid[1e-3:1e-2:1e-3],np.mgrid[1e-2:1e-1:1e-2],np.mgrid[1e-1:1:1e-1],np.mgrid[1:10:1])) 
tau_d=1/K # Desorption timescales arrays
tau_e=R/(P_f*C_out*V_w) # Evaporation timescale
tau_f=exc_s**0.5/xi # Flow timescale
n_d=len(K) # Number of desorption timescales
# Array that will contain an integer, to characterize the different ways the DIB becomes unstable
regime=np.ones((n_d,n_d))*(-3)
for i in range(n_d):
    for j in range(n_d):
        eq_d=0
        eq_DIB=0
        non_phys=0
        stop=np.zeros(9) # stop is an array full of zero except at one position where it is equal to 1, this position is a reference to the regime the DIB has experienced.
        k_off_b=K[i]
        k_off_m=K[j]
        sim=Simulation(exc_s,gam_0,P_f,xi,k_on,k_off_m,k_off_b,C_out,C_in,mu_b) 
        sim.values_arrays_init()
        ICs=sim.initialization(exc_l0,exc_r0,exc_b0)
        """
        We look if the system has fixed points, if not the result is a string that says that the system is unstable.
        """
        x=sim.fixed_points(exc_b0,ICs) 
        if isinstance(x,str):
            t_f=1 # Initial duration of the simulation
            v_l=0.5*a_exp[0]**3*f(theta_left_exp[0]) # Initial volume of left droplet
            v_r=0.5*a_exp[0]**3*f(theta_right_exp[0]) # Initial volume of right droplet
            a=a_exp[0] # Initial value of bilayer radius
            theta_l=theta_left_exp[0] # Initial value of left polar angle 
            theta_r=theta_right_exp[0] # Initial value of right polar angle
            # All the folllowing values can be modified.
            theta_min=95*np.pi/180 # Minimum value of angle, if achieved, stop the simulation
            theta_max=150*np.pi/180 # Maximum value of angle, if achieved, stop the simulation
            v_l_min=50e-18 # Minimum volume of left droplet, if achieved stop the simulation
            v_r_min=50e-18 # Minimum volume of right droplet, if achieved stop the simulation
            a_lim=0.2e-6 # Minimum volume of bilayer radius, if achieved stop the simulation
            while (v_l>v_l_min and v_r>v_r_min and a>a_lim and theta_l<theta_max and theta_r<theta_max and theta_l> theta_min and theta_r>theta_min and eq_d<1 and non_phys<1):
                # While some observabes keep values in a certain range, we simulate the system.                
                sim.values_arrays_init()
                try:
                    sim.evolution(t_f,1e-2,ICs)
                    if (np.absolute(0.5*sim.a[-1]**3*f(sim.theta_l[-1])/v_l-1
                    )<1e-6 and np.absolute(0.5*sim.a[-1]**3*f(sim.theta_r[-1])/v_r-1
                    )<1e-6 and 0.5*sim.a[-1]**3*f(sim.theta_l[-1]
                    )>v_l_min and 0.5*sim.a[-1]**3*f(sim.theta_r[-1]
                    )>v_r_min and np.absolute(sim.a[-1]/a-1)<1e-6):
                        # Eventhough we said the system has no fixed point, if the system is really slightly perturbed, we can consider it is at equilibrium
                        eq_d=1 # If the volume of the droplets and the radius of the bilaer are not modified by 1e-6 in relative value (can be changed) the droplets are said at equilibrium
                        if sim.a[-1]>a_lim:
                            eq_DIB=1 # BUT the DIB is only at equilibrium if the bilayer radius is large enought
                    v_l=0.5*sim.a[-1]**3*f(sim.theta_l[-1])
                    v_r=0.5*sim.a[-1]**3*f(sim.theta_r[-1])
                    a=sim.a[-1]
                    theta_l=sim.theta_l[-1]
                    theta_r=sim.theta_r[-1]
                except ValueError:
                    non_phys=1
                t_f+=1 # We increase the total duration of the simulation and start again.
            if eq_d>0:
                if eq_DIB>0:
                    regime[i,j]=-1 # -1: integer for stable DIBs
                else:
                    regime[i,j]=-2 # -2: integer for stable droplets but unstable DIB
            else:
                if non_phys>0:
                    # If a Value Error was raised during the simulation, it means we simulate the system too long.
                    sim.values_arrays_init()
                    sim.evolution(t_f-2,1e-2,ICs)
                while regime[i,j]==-3:
                    # We compute first derivatives of some arrays.
                    var_a=np.gradient(sim.a)
                    var_theta_l=np.gradient(sim.theta_l)
                    var_theta_r=np.gradient(sim.theta_r)
                    # If the minimum of a is close to the initial value (in a sense that can be modified), we consider it has only increased.
                    if np.min(sim.a)>3.56e-6:
                        stop[0]=1 # 1: regime of bilayer expansion.
                    # If the maximum of a is close to its initial value (in a sense that can be modified), we consider it has only decreased.
                    elif np.max(sim.a)<3.8e-6:
                        # If the maximum of angles is not too far from the initial value, we consider angles have only decreased.
                        if max(np.max(sim.theta_l),np.max(sim.theta_r))<2.10:
                            stop[3]=1 # regime 3: both angles and bialyer radius decrease, this is the crossover behaviour.
                        # If the minimum of angles is not to far from the initial value; we consider angles have only increased.
                        elif min(np.min(sim.theta_l),np.min(sim.theta_r))>1.85:
                            stop[5]=1 # 5: unzipping
                        else:
                            if (var_theta_l[0]>0 and var_theta_r[0]>0):
                                stop[6]=1 #6: bilayer radius decreases and angles first increase and then decrease.
                            else:
                                stop[4]=1 #4: bilayer radius decreases and angles first decrease and finally increase.
                    else:
                        if var_a[0]>0:
                            if max(np.max(sim.theta_l),np.max(sim.theta_r))<2.10:
                                stop[1]=1 #1: Up/Down crossover with decreasing angles
                            else:
                                stop[2]=1 #2: Up/Down crossover with increasing angles
                        else:
                            if max(np.max(sim.theta_l),np.max(sim.theta_r))<2.10:
                                stop[8]=1 #8: Down/Up crossover with decreasing angles
                            else:
                                stop[7]=1 #7: Down/Up crossover with increasing angles
                    regime[i,j]=np.where(stop>0)[0] # We save the integer in the regime array
        else:
            if x[0]>a_lim:
                regime[i,j]=-1 # If we had a fixed point, it means that we have stable DIB or droplets, depending of the size of the bilayer.
            else:
                regime[i,j]=-2
            
"""
We create 2D arrays for the plot after: in (i,j), beta_d_b is the array of beta_des for the bilayer
and beta_d_m is the array of beta_des for the monolayers.
"""
beta_d_m=np.zeros((n_d,n_d))
beta_d_b=np.zeros((n_d,n_d))
for i in range(n_d):
    for j in range(n_d):
        beta_d_b[i,j]=tau_d[i]/tau_e
        beta_d_m[i,j]=tau_d[j]/tau_e

"""
We save data. The filenames can be modified in variables path...
"""
beta_d_b_path=
beta_d_m_path=
regime_path=
np.savetxt(beta_d_b_path,beta_d_b)
np.savetxt(beta_d_m_path,beta_d_m)
np.savetxt(regime_path,regime)