# -*- coding: utf-8 -*-
"""
Created on Fri Jul  8 08:58:35 2016

@author: softmatter
"""

import numpy as np
import matplotlib.pyplot as plt

"""
Declare paths to txt documents to find values of:
- beta_des for monolayer (beta_d_m) ;
- beta _des for bilayer (beta_d_b) ;
- regimes encountered (regime).
"""
path_beta_d_b=
path_beta_d_m=
path_regime=
beta_d_b=np.genfromtxt(path_beta_d_b) # Array of values of beta_des for the bilayer
beta_d_m=np.genfromtxt(path_beta_d_m) # Array of values of beta_des for the monolayers
regime=np.genfromtxt(path_regime) # Array of indexes indicating the regime of evolution
n_d_b=np.shape(beta_d_b)[0] # Number of simulated beta_des for the bilayer
n_d_m=np.shape(beta_d_b)[1] # Number of simulated beta_des for the monolayer

""" 
We only keep values of beta_des which are sufficiently far not to make the plot too crowdy.
The value "step" can be modified.
"""
step=0.1*np.log(10)
I=[0]
for i in range(n_d_b):
    if np.absolute(np.log(beta_d_b[i,0])-np.log(beta_d_b[I[-1],0]))>step:
        I+=[i]
    elif regime[i,-1]!=regime[I[-1],-1]:
        I+=[i]
J=[0]
for i in range(n_d_m):
    if np.absolute(np.log(beta_d_m[0,i])-np.log(beta_d_m[0,J[-1]]))>step:
        J+=[i]
    elif regime[0,i]!=regime[0,I[-1]]:
        J+=[i]
    
regime=regime[I,:][:,J]
beta_d_m=beta_d_m[I,:][:,J]
beta_d_b=beta_d_b[I,:][:,J]

"""
We gather indexes of arrays in other arrays according to the regime 
encountered in the simulation.
"""
# We gather crossover behaviours
regime[regime==6]=4
regime[regime==1]=2
regime[regime==7]=8
I_eq, J_eq=np.where(regime==-2) # Droplets are stable but the DIB is unstable
I_eq_DIB,J_eq_DIB=np.where(regime==-1) # The DIB is stable.
I_0, J_0=np.where(regime==0) # Bilayer expansion
I_2, J_2=np.where(regime==2) # Up-Down crossover
I_3, J_3=np.where(regime==3) # Crossover behaviour
I_4, J_4=np.where(regime==4) # Regime in which the bilayer radius decreases and polar angles do not have monotonic behaviour
I_5, J_5=np.where(regime==5) # Unzipping
I_8, J_8=np.where(regime==8) # Down-Up crossover

"""
Plot.
"""
fig,ax=plt.subplots() # The markersize (ms), meaning the size of the points representing data in the plot can be modified, so as their shape and their color.
ax.loglog(beta_d_m[I_0,J_0],beta_d_b[I_0,J_0],marker='o',ms=10,linestyle='',color='Crimson',label='Bilayer expansion')
ax.loglog(beta_d_m[I_2,J_2],beta_d_b[I_2,J_2],marker='s',ms=10,linestyle='',color='LightGreen')
ax.loglog(beta_d_m[I_3,J_3],beta_d_b[I_3,J_3],marker='s',ms=10,linestyle='',color='LightGreen',label='Crossover')
ax.loglog(beta_d_m[I_4,J_4],beta_d_b[I_4,J_4],marker='s',ms=10,linestyle='',color='LightGreen')
ax.loglog(beta_d_m[I_5,J_5],beta_d_b[I_5,J_5],marker='^',ms=10,linestyle='',color='DarkSlateBlue',label='Unzipping')
ax.loglog(beta_d_m[I_8,J_8],beta_d_b[I_8,J_8],marker='s',ms=10,ls='',color='LightGreen')
ax.set_xlabel(r'$\beta_{des}^{m}$',fontsize=40) # fontsize can be modified.
ax.set_ylabel(r'$\beta_{des}^{b}$',fontsize=40) # fontsize can be modified.
ax.tick_params(labelsize=35) # The labelsize of the numbers on axes can be modified.
ax.set_xlim(7e-4,1000) # The x-limits of the plot can be modified.
ax.set_ylim(7e-4,1000) # The y limits of the plot can be modified.
# Plot of the legend in ncol=3 columns, the position is set by the loc argument and the bbox_to_anchor argument, the size of characters can be modified too.
ax.legend(loc='upper center',bbox_to_anchor=(0.5,1.5),prop={'size':24},fancybox=True,ncol=3)
fig.show()