# -*- coding: utf-8 -*-
"""
Created on Tue Aug 09 11:12:47 2016

@author: Benjamin GUISELIN
"""

from __future__ import division
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import math

# Geometrical functions:
"""
f is the function so that the volume of the two droplets with bilayer radius a and common polar angle theta is V=a**3*f(theta).
f_prime is the derivative of f. 
g is the function so that the area of the two droplets with bilayer radius a and common polar angle theta is S=a**2*g(theta).
g_prime is the derivative of g.
"""

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
    
# Fixed parameters:
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_m, k_on_b, k_off_m, k_off_b,
                 C_out_l, C_out_r, C_in_l,C_in_r, mu_b, exc_b):
        
        """
        Create a framework for the simulation with free parameters.
        
        Arguments:
        - k_off_m: monolayer desorption constant of lipids [/s].
        - k_off_b: bilayer desorption constant of lipids [/s].
        - k_on_m: monolayer adsorption constant of lipids [/s].
        - k_on_b: blayer adsorption constant of lipids [/s].
        - exc_s: total number of available sites 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]
        - C_out_l: left outside osmolarity [mol/m^3], related to the osmotic pressure
          by C_out=Pi_out/(N_a*k_B*T)
        - C_out_r: right outside osmolarity [mol/m^3]
        - C_in_l: left inside osmolarity [mol/m^3]
        - C_in_r: right inside osmolarity [mol/m^3]
        - mu_b: chemical potential of lipids on the bilayer [k_{B}T]
        - exc_b: bilayer lipid density [/m^2]
        
        Raises:
        - TypeError if parameters are not floats 
        - ValueError if parameters are negative (except C_out or mu_b) or if 
        bilayer surface excess is above the maximum packing density
        """
        
        try:
            self.k_off_m=float(k_off_m)
            self.k_off_b=float(k_off_b)
            self.k_on_m=float(k_on_m)
            self.k_on_b=float(k_on_b)
            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_l=float(C_out_l)
            self.C_out_r=float(C_out_r)
            self.C_in_l=float(C_in_l)
            self.C_in_r=float(C_in_r)
            self.mu_b=float(mu_b)
            self.exc_b=float(exc_b)
        except (TypeError, ValueError,AssertionError):
            raise TypeError ('Simulation parameters must be floats.')
        if (self.k_off_m<0 or self.k_off_b<0):
            raise ValueError ('Desorption constants must be positive.')
        if (self.k_on_m <0 or self.k_on_b <0):
            raise ValueError ('Adsorption constants must be positive.')
        if (self.C_in_l<0 or self.C_in_r <0):
            raise ValueError ('Inside osmolarity must be positive.')
        if self.exc_s <0:
            raise ValueError ('Langmuir parameter must be positive.')
        if self.gam_0<0:
            raise ValueError ('Surface tension must be positive.')
        if self.P_f <0:
            raise ValueError ('Water permeation coefficient must be positive.')
        if self.xi <0:
            raise ValueError ('''Diffusion constant of lipids on layers must
            be positive.''')
        if (self.exc_b<0 or self.exc_b>self.exc_s):
            raise ValueError ('''Bilayer lipid density must be positive and 
            lower than the maximum packing density.''')

    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].
        - exc_l: left monolayer lipid density [/m^2].
        - exc_r: right monolayer lipid density [/m^2].
        - V_l: left droplet volume [m^3].
        - V_r: right droplet volume [m^3].
        - N_l: left monolayer lipid number.
        - N_r: right monolayer lipid number.
        - N_b: bilayer lipid number.
        - N_des_l: number of lipids desorbing from the left monolayer.
        - N_des_r: number of lipids desorbing from the right monolayer.
        - N_des_b: number of lipids desorbing from the bilayer.
        - N_flow_l_b: number of flowing lipids from the left monolayer to the
        bilayer.
        - N_flow_r_b: number of flowing lipids from the right monolayer to the
        bilayer.
        - data: experimental data.
        """
        
        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.V_l=np.array([])
        self.V_r=np.array([])
        self.N_l=np.array([])
        self.N_r=np.array([])
        self.N_b=np.array([])
        self.N_des_l=np.array([])
        self.N_des_r=np.array([])
        self.N_des_b=np.array([])
        self.N_flow_l_b=np.array([])
        self.N_flow_r_b=np.array([])
        self.data=np.array([])
        
    def initialization(self,a0,theta_l0,theta_r0,exc_l0,exc_r0):
        
        """
        Create initial conditions for the simulation.
        
        Arguments:
        - a0: initial bilayer radius [m].
        - theta_l0: initial left droplet polar angle [rad].
        - theta_r0: initial right droplet polar angle [rad].
        - exc_l0: initial lipid density on the left monolayer [/m^2].
        - exc_r0: initial lipid density on the right monolayer [/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 lipid density [/m^2].
            - exc_r0: initial right monolayer lipid density [/m^2].
        
        Raise:
        - TypeError if initial conditions are not floats.
        - ValueError if initial conditions are negative, or if initial lipid 
        densities are above the maximum packing density or if polar angles are
        not in the range [90°,180°].
        """
        
        try:
            a0=float(a0)
            theta_l0=float(theta_l0)
            theta_r0=float(theta_r0)
            exc_l0=float(exc_l0)
            exc_r0=float(exc_r0)
        except (TypeError,ValueError,AssertionError):
            raise TypeError ('Initial conditions must be floats.')
        if (exc_l0<0 or exc_r0<0 or a0<0 or theta_l0<0 or theta_r0<0):
            raise ValueError ('Initial conditions cannot be negative.')
        if max(exc_l0,exc_r0)>self.exc_s:
            raise ValueError ('''Initial lipid densities cannot be above the 
            maximum packing density.''')
        if (min(theta_l0,theta_r0)<np.pi/2. or max(theta_l0,theta_r0)>np.pi):
            raise ValueError ('''Droplet polar angles must be in the range 
            [90°,180°].''')
        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])
        
    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 lipid density [/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'.
        """
        
        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_l*V_w-self.C_in_l*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_r*V_w-self.C_in_r*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):
        
        """
        Compute the time derivative of bilayer radius.
        
        Arguments:
        - a: bilayer radius [m].
        - exc_l: left monolayer lipid density [/m^2].
        - exc_r: right monolayer lipid density [/m^2].
    
        Return:
        - Time derivative of bilayer radius.
        
        Raise:
        - ValueError if at least one lipid density is negative or above the 
        maximum packing density.
        """
        if (math.isinf(np.log(exc_l/(self.exc_s-exc_l))) or math.isnan(np.log(
        exc_l/(self.exc_s-exc_l)))):
            raise ValueError('''Left monolayer lipid density must be positive 
            and lower than the maximum packing density.''')
        if (math.isinf(np.log(exc_r/(self.exc_s-exc_r))) or math.isnan(np.log(
        exc_r/(self.exc_s-exc_r)))):
            raise ValueError('''Right monolayer lipid density must be positive 
            and lower than the maximum packing density.''')
        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*self.exc_b
            )-self.k_off_b*a/2+self.k_on_b*(a/2.)*(self.exc_s-self.exc_b
            )/self.exc_b
        
    def evolution_theta(self,a,theta_l,theta_r,exc_l,exc_r,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 lipid density [/m^2].
        - exc_r: right monolayer lipid density [/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' or if at least one surface 
        lipid density is negative or above the maxmimum packing density.
        """ 
        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*f(theta_l)*self.evolution_a(
            a,exc_l,exc_r))
        if side=='r':
            return (1/(a**3*f_prime(theta_r)))*(2*self.evolution_volume(
            a,theta_r,exc_r,side)-3*a**2*f(theta_r)*self.evolution_a(
            a,exc_l,exc_r))
        else:
            raise ValueError('''side must be 'l' or 'r'.''')
                
    def evolution_exc_m(self,a,theta_l,theta_r,exc_l,exc_r,side):
         
        """
        Compute the time derivative of monolayer surface excess.
        
        Arguments:
        - a: bilayer radius [m].
        - theta_l: left droplet polar angle [rad].
        - theta_r: right droplet polar angle [rad].
        - exc_l: left monolayer lipid density [/m^2].
        - exc_r: right monolayer lipid density [/m^2].
        - side: 'l' for the left droplet and 'r' for the right droplet.
    
        Return:
        - Time derivative of monolayer lipid density.
    
        Raises:
        - TypeError if side is not a string.
        - ValueError if side is not 'l' or 'r' or if at least one surface 
        lipid density is negative.
        """
        
        try:
            side=str(side)
        except (TypeError,ValueError,AssertionError):
            raise TypeError('side must be a string.')
        if side=='l':
            return self.k_on_m*(self.exc_s-exc_l
            )-self.k_off_m*exc_l-exc_l*(2*self.evolution_a(a,exc_l,exc_r)/a+(
            g_prime(theta_l)/g(theta_l))*self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,side))+4*np.pi*self.xi*(
            self.mu_b-np.log(exc_l/(self.exc_s-exc_l)))/(a*g(theta_l))
        if side=='r':
            return self.k_on_m*(self.exc_s-exc_r
            )-self.k_off_m*exc_r-exc_r*(2*self.evolution_a(a,exc_l,exc_r)/a+(
            g_prime(theta_r)/g(theta_r))*self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,side))+4*np.pi*self.xi*(
            self.mu_b-np.log(exc_r/(self.exc_s-exc_r)))/(a*g(theta_r))
        else:
            raise ValueError('''side must be 'l' or 'r'.''')
        
    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 for the bilayer radius, droplet polar
        angles (left/right) and initial monolayer lipid densities (left/right).
        
        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 lipid density 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.''')
        
        # Evolution operator of the system
        def lindiff(y,t):
            a, theta_l, theta_r,exc_l, exc_r = y
            derivs = [self.evolution_a(a,exc_l,exc_r),self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,'l'),self.evolution_theta(
            a,theta_l,theta_r,exc_l,exc_r,'r'),self.evolution_exc_m(
            a,theta_l,theta_r,exc_l,exc_r,'l'),self.evolution_exc_m(
            a,theta_l,theta_r,exc_l,exc_r,'r')]
            return derivs  
            
        # 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.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
        self.V_l=0.5*self.a**3*f(self.theta_l)
        self.V_r=0.5*self.a**3*f(self.theta_r)
        self.N_l=self.exc_l*self.S_l
        self.N_r=self.exc_r*self.S_r
        self.N_b=self.exc_b*self.A*2
        
        """
        Number of desorbing lipids from the layers computed by a numerical 
        integration (trapezoidal rule) of the desorption term. At t=0, the 
        number of desorbing lipids is equal to 0.
        """
        
        self.N_des_l=integrate.cumtrapz(self.k_off_m*self.exc_l*self.S_l,
                                        self.t,initial=0) 
        self.N_des_r=integrate.cumtrapz(self.k_off_m*self.exc_r*self.S_r,
                                        self.t,initial=0)
        self.N_des_b=integrate.cumtrapz(self.k_off_b*self.exc_b*self.A*2,
                                        self.t,initial=0)
        
        """
        Number of flowing lipids from the monolayers to the bilayer computed by
        a numerical integration (trapezoidal rule) of the flow term. At t=0, 
        the number of flowing lipids is equal to 0.
        """
        
        self.N_flow_l_b=integrate.cumtrapz(-2*np.pi*self.a*self.xi*(
        self.mu_b-np.log(self.exc_l/(self.exc_s-self.exc_l))),self.t,initial=0)
        self.N_flow_r_b=integrate.cumtrapz(-2*np.pi*self.a*self.xi*(
        self.mu_b-np.log(self.exc_r/(self.exc_s-self.exc_r))),self.t,initial=0)
                            
    def data_importation(self,data_name):
        
        """
        Import data simulations from document called data_name (extension 
        required). The document must be in the same folder as the code (otherwise please
        include the path to the file) and must be an array of 4 columns: 
        (time, left monolayer, right monolayer, bilayer). 
        Areas must be given in m^2 and time in seconds.
        
        Argument:
        - data_name: name of the document containing experimental data.
        
        Returns:
        - TypeError if data_name is not a string.
        - ValueError if the document is not in the folder or if the data array
        does not have the right dimensions.
        """
        
        try:
            data_name=str(data_name)
        except (TypeError,ValueError,AssertionError):
            raise TypeError ('''The name of the document containing data must 
            be a string.''')
        try:
            data=np.genfromtxt(str(data_name))
            self.data=np.zeros((data.shape[0],))
        except IOError:
            raise ValueError ('Data cannot be found. Check the name and the path again.')
        if np.absolute(data.shape[1]-4)>0:
            raise ValueError ('Data array does not have the right dimensions.')
        self.data[:,-1]=self.data[:,-2]*2/(self.data[:,1]+self.data[:,2])
    
    def error(self,sigma=5e-12):
            
        """
        Compute simulation error with respect to experimental data.
        
        Argument:
        - sigma: standard error of experimental areas (5e-12 by default).
        
        Return:
        - Simulation error.
        
        Raise:
        - ValueError if simulation duration is less than experiment duration 
        t_exp [s].
        """
        
        try:
            time=self.data[:,0]
            monolayer_left=self.data[:,1]
            monolayer_right=self.data[:,2]
            bilayer=self.data[:,3]
            area_star=self.data[:,4]
        except IndexError:
            raise ValueError ('Experimental data have not been imported yet.')
            
        # Gathering of simulation points corresponding to experimental points        
        J=[]
        for i in range(len(time)):
            j=np.argmin(np.absolute(time[i]-self.t))
            if np.absolute(time[i]-self.t[j])<=self.dt/2:
                J+=[j]   
        
        # Error computation 
        if len(J)<len(time):
            raise ValueError('''Simulation duration is less than experiment
            duration.''')
        else:
            E=(np.pi*self.a[J]**2-bilayer)**2/(2*sigma**2) # error on bilayer area
            F=(self.a[J]**2*0.5*g(self.theta_l[J])-monolayer_left)**2/(
            2*sigma**2) # error on left monolayer area
            G=(self.a[J]**2*0.5*g(self.theta_r[J])-monolayer_right)**2/(
            2*sigma**2) # error on right monolayer area
            sigma_star=area_star*np.sqrt((sigma/bilayer)**2+2*(sigma/(
            monolayer_left+monolayer_right))**2) # uncertainty on aspect ratio
            E_star=(2*np.pi/(0.5*g(self.theta_l[J])+0.5*g(self.theta_r[J])
            )-area_star)**2/(2*sigma_star**2) # error on aspect ratio
            return np.sum(E)+np.sum(F)+np.sum(G)+np.sum(E_star)

    def plot_solutions(self,every=100,lw=4,ms=20,ls=35,fs=40):
            
        """
        Plot simulation results.
        
        Arguments:
        - every: spacing between markers on line plots (100 by default).
        - lw: linewidth for plots (2 by default).
        - ms: markersize for plots (20 by default).
        - ls: size for the ticks of the plots (35 by default).
        - fs: fontsiwe for the titles of the axes (40 by default).
               
        Returns:
        - Plot of bilayer radius a [microns] 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 lipid densities [in units of exc_s] as a function of time [s].
        - Plot lipid numbers on layers as a function of time [s].
        - Plot droplets' radii [microns] as a function of time [s].
        - Plot volumes [fL] as a function of time [s].
        - Plot exchanged lipid numbers as a function of time [s].
        
        Raises:
        - TypeError if parameters are not integers.
        - ValueError if parameters are not positive.
        """
        
        try:
            every=int(every)
            lw=int(l)
            ms=int(ms)
            ls=int(ls)
            fs=int(fs)
        except (TypeError,ValueError,AssertionError):
            raise TypeError ('Parameters must be integers.')
        if (every<0 or lw<0 or ms<0 or ls<0 or fs<0):
            raise ValueError ('Parameters must be positive.')
        
        # Bilayer radius:
        fig1,ax1=plt.subplots()
        ax1.plot(self.t,self.a*1e6,marker='d',linestyle='-',markersize=ms,
                 linewidth=lw,markevery=every,color='lightcoral') # 'd' = diamond
        ax1.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax1.set_ylabel(r'Bilayer radius ($\mu m$)',fontsize=fs)
        ax1.tick_params(labelsize=ls)
        fig1.show()
        
        # Polar angles: 
        fig2,ax2=plt.subplots()
        ax2.plot(self.t,self.theta_l*180/np.pi,markevery=every,marker='v',
                 markersize=ms,color='mediumseagreen',linewidth=lw,
                 linestyle='-',label=r'$\theta_{l}^{th}$') # left polar angle, 'v' = down triangle
        ax2.plot(self.t,self.theta_r*180/np.pi,linestyle='-',marker='p',
                 markersize=ms,color='steelblue',linewidth=lw,markevery=every,
                 label=r'$\theta_{r}^{th}$') # right polar angle, 'p' = pentagon
        ax2.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax2.set_ylabel(r'Angle ($^{\circ}$)',fontsize=fs)
        ax2.tick_params(labelsize=ls)
        ax2.legend(loc='upper center',bbox_to_anchor=(0.5,1.3),ncol=2,
                   prop={'size':45},fancybox=True)
        fig2.show()
        
        # Surface areas and surface aspect ratio:
        fig3,ax3=plt.subplots()
        ax4=ax3.twinx() # Double y-axis to plot the aspect ratio on the same plot
        ax3.plot(self.t,self.S_l*1e12,marker='v',markevery=every,markersize=ms,
                 color='mediumseagreen',linestyle='-',linewidth=lw,
                 label='Left monolayer')
        ax3.plot(self.t,self.S_r*1e12,marker='p',linestyle='-',markevery=every,
                 markersize=ms,color='steelblue',linewidth=lw,
                 label='Right monolayer')
        ax3.plot(self.t,self.A*1e12,marker='d',linestyle='-',markevery=every,
                 markersize=ms,color='lightcoral',linewidth=lw,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=lw,
                 label=r'$A_{th}^{*}$')
        ax3.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax3.set_ylabel(r'Surface area ($\mu m^{2}$)',fontsize=fs)
        ax4.set_ylabel(r'$A^{*}=\frac{2A}{S_{l}+S_{r}}$',fontsize=fs)
        ax4.tick_params(labelsize=ls)
        ax3.tick_params(labelsize=ls)
        ax4.set_yticks([0.0,0.2,0.4,0.6,0.8,1.0])
        try:
            time=self.data[:,0]
            monolayer_left=self.data[:,1]
            monolayer_right=self.data[:,2]
            bilayer=self.data[:,3]
            area_star=self.data[:,4]
            ax3.plot(time,monolayer_left*1e12,linestyle='-',color='olive',
                     marker='^',linewidth=lw,markersize=ms,fillstyle='none',
                     markeredgewidth=ms/3,label='Left monolayer (exp)')
            ax3.plot(time,monolayer_right*1e12,linestyle='-',fillstyle='none',
                     color='mediumturquoise',linewidth=lw,markersize=ms,
                     marker='o',markeredgewidth=4,
                     label='Right monolayer (exp)')
            ax3.plot(time,bilayer*1e12,linestyle='-',fillstyle='none',
                     marker='s',color='tomato',linewidth=lw,markeredgewidth=ms/3,
                     markersize=ms,label='Bilayer (exp)')
            ax4.plot(time,area_star,color='black',linestyle='dashed',
                     linewidth=lw,label=r'$A_{exp}^{*}$')
            ax4.legend(loc='upper left',prop={'size':32},fancybox=True,ncol=2,
                       bbox_to_anchor=(0.1,1.2))
            ax3.legend(loc='upper right',prop={'size':18},ncol=3,fancybox=True,
                       bbox_to_anchor=(1,1.2))
            fig3.show()
        except IndexError:
            ax4.legend(loc='upper left',prop={'size':32},fancybox=True,
                       bbox_to_anchor=(0.1,1.2))
            ax3.legend(loc='upper right',prop={'size':18},ncol=3,fancybox=True,
                       bbox_to_anchor=(1,1.15))
            fig3.show()
              
        # Lipid densities:
        fig5,ax5=plt.subplots()
        ax5.plot(self.t,self.exc_l/self.exc_s,color='mediumseagreen',
                 marker='v',markersize=ms,markevery=every,linestyle='-',
                 linewidth=lw,label=r'$\Gamma_{m}^{l}$') # left monolayer lipid density
        ax5.plot(self.t,self.exc_r/self.exc_s,marker='p',linestyle='-',
                 markersize=ms,color='steelblue',linewidth=lw,
                 markevery=(int(every/2),every),label=r'$\Gamma_{m}^{r}$') # right monolayer lipid density
        ax5.plot(self.t,np.ones(len(self.t))*self.exc_b/self.exc_s,marker='d',
                 linestyle='-',markersize=ms,markevery=every,linewidth=lw,
                 color='lightcoral',label=r'$\Gamma_{b}$') # bilayer lipid density
        ax5.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax5.set_ylabel(r'Surface density ($\Gamma_{\infty}$)',fontsize=fs)
        ax5.legend(loc='upper center',bbox_to_anchor=(0.5,1.35),ncol=3,
                   prop={'size':48},fancybox=True)
        ax5.tick_params(labelsize=ls)
        fig5.show()
        
        # Lipid numbers:
        k=int(np.log(self.N_l[0])/np.log(10)) # Power of lipid numbers in scientific notation
        fig6,ax6=plt.subplots()
        ax6.plot(self.t,self.N_l/10**k,color='mediumseagreen',marker='v',
                 markersize=ms,markevery=every,linestyle='-',linewidth=lw,
                 label=r'$N_{m}^{l}$') # left monolayer lipid number
        ax6.plot(self.t,self.N_r/10**k,marker='p',linestyle='-',linewidth=lw,
                 markersize=ms,markevery=every,color='steelblue',
                 label=r'$N_{m}^{r}$') # right monolayer lipid number
        ax6.plot(self.t,self.N_b/10**k,marker='d',markevery=every,linewidth=lw,
                 linestyle='-',markersize=ms,color='lightcoral',
                 label=r'$N_{b}$') # bilayer lipid number
        ax6.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax6.set_ylabel(r'Lipid number ($\times \ 10^{}$)'.format(k),
                       fontsize=fs)
        ax6.tick_params(labelsize=ls)
        ax6.legend(loc='upper center',prop={'size':40},fancybox=True,ncol=3,
                   bbox_to_anchor=(0.5,1.3))
        fig6.show()
        
        # Radii:
        fig7,ax7=plt.subplots()
        ax7.plot(self.t,self.R_l*1e6,linestyle='-',marker='^',markevery=every,
                 color='mediumseagreen',markersize=ms,linewidth=lw,
                 label=r'$R_{l}$') # left droplet radius
        ax7.plot(self.t,self.R_r*1e6,linestyle='-',marker='p',color='steelblue',
                markersize=ms,linewidth=lw,markevery=every,label=r'$R_{r}$') # right droplet radius
        ax7.tick_params(labelsize=ls) 
        ax7.legend(loc='upper center',prop={'size':40},fancybox=True,ncol=2,
                   bbox_to_anchor=(0.5,1.25))
        ax7.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax7.set_ylabel(r'Droplet radius ($\mu$m)',fontsize=fs)
        fig7.show()
        
        # Volumes:
        fig8,ax8=plt.subplots()
        ax8.plot(self.t,1e18*self.V_l,color='mediumseagreen',marker='v',
                 markersize=ms,markevery=every,linestyle='-',linewidth=lw,
                 label=r'$V_{l}$') # left droplet volume
        ax8.plot(self.t,1e18*self.V_r,marker='p',linestyle='-',markersize=ms,
                 markevery=every,color='steelblue',linewidth=lw,
                 label=r'$V_{r}$') # right droplet volume
        ax8.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax8.set_ylabel('Droplet volume (fL)',fontsize=fs)
        ax8.tick_params(labelsize=ls)
        ax8.legend(loc='upper center',prop={'size':48},fancybox=True,ncol=2,
                   bbox_to_anchor=(0.5,1.3))
        fig8.show()
        
        # Exchanged lipid numbers:
        fig9, ax9=plt.subplots()
        ax9.plot(self.t,self.N_des_l/10**k,marker='v',color='mediumseagreen',
                 linestyle='-',linewidth=lw,markersize=ms,markevery=every,
                 label=r'$N_{des}^{l}$') # lipids desorbing from the left monolayer
        ax9.plot(self.t,self.N_des_r/10**k,marker='p',color='steelblue',
                 linewidth=lw,linestyle='-',markersize=ms,
                 markevery=(int(every/2),every),label=r'$N_{des}^{r}$') # lipids desorbing from the right monolayer
        ax9.plot(self.t,self.N_des_b/10**k,marker='d',color='lightcoral',
                 linewidth=lw,linestyle='-',markersize=ms,markevery=every,
                 label=r'$N_{des}^{b}$') # lipids desorbing from the bilayer monolayer
        ax9.plot(self.t,self.N_flow_l_b/10**k,marker='<',color='darkorange',
                 linewidth=lw,linestyle='-',markersize=ms,markevery=every,
                 label=r'$N_{flow}^{l\rightarrow b}$') # '<' = left triangle, lipids flowing from the left monolayer to the bilayer
        ax9.plot(self.t,self.N_flow_r_b/10**k,marker='>',color='mediumorchid',
                 linewidth=lw,linestyle='-',markersize=ms,
                 markevery=(int(every/2),every),
                 label=r'$N_{flow}^{r\rightarrow b}$') # '>' = right triangle, lipids flowing from the right monolayer to the bilayer
        ax9.set_xlabel('Evaporation time (s)',fontsize=fs)
        ax9.set_ylabel(r'Lipid number ($\times \ 10^{}$)'.format(k),
                       fontsize=fs)
        ax9.tick_params(labelsize=ls)
        ax9.legend(loc='upper left',prop={'size':36},fancybox=True,ncol=5,
                   bbox_to_anchor=(-0.15,1.35))
        fig9.show()

# Enter the numerical values of the parameters:         
exc_s=
gam_0=
k_off_m=
k_off_b=
k_on_m=
k_on_b=
C_in_l=
C_in_r=
xi=
P_f=
C_out_l=
C_out_r=
exc_b=
mu_b=
sim=Simulation(exc_s, gam_0, P_f, xi, k_on_m, k_on_b, k_off_m, k_off_b, 
               C_out_l, C_out_r, C_in_l,C_in_r, mu_b, exc_b)
sim.values_arrays_init()
# Enter the initial conditions:
a0=
theta_l0=
theta_r0=
exc_l0=
exc_r0=
ICs=sim.initialization(a0,theta_l0,theta_r0,exc_l0,exc_r0)
# Enter the simulation duration and the timestep:
t_f=
dt=
sim.evolution(t_f,dt,ICs)
sim.plot_solutions()
