# -*- coding: utf-8 -*-
"""
Created on Tue Oct 24 08:42:55 2017

@author: Benjamin
"""

from __future__ import division
from collections import Counter
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as col
import scipy.interpolate
import numpy as np
import numpy.linalg as linalg
from matplotlib.ticker import FormatStrFormatter


def create_proxy(label):
    line = matplotlib.lines.Line2D([0], [0], linestyle='',
                mec='none', marker=label[0], color=col.hsv_to_rgb((label[1],label[2],0.7)))
    return line

L=np.mgrid[0:1:80j]
L=np.delete(L,0)
L=np.delete(L,1) # and get rid 1 in the array
exc_s=2.3e18
gam_0=25e-3
k_on=0
k_off=5e-3
k_off_m=k_off
k_off_b=k_off_m
exc_l0=1.42e18
exc_r0=1.42e18
exc_b0=1.42e18
G=exc_b0
k_on=k_off_m*exc_b0/(exc_s-exc_b0) # Equilibrium value of k_on_m
k_on_m=k_on
k_on_b=k_on
C_in=0
xi=9.85e10
k_B=1.38e-23
T=298
N_a=6.02e23
P_f=80e-6
mu_b=np.log(k_on/k_off_m) # Equilibrium value of mu_b
C=np.array([1e-4,1.5e-4,2e-4,3e-4,5e-4,8e-4,1e-3,1.5e-3,2e-3,3e-3,5e-3,8e-3,1e-2,1.5e-2,
            2e-2,3e-2,5e-2,8e-4,1e-1,2e-1,3e-1,5e-1,8e-1,1,2,3,5,8,10,13,15,20]) # Outside osmolarities that are checked
C=np.flipud(C)
C_grid=C

fp = open("EU_tuples.txt","r")
adsorption = fp.read()
fp.close()
absorption = [map(lambda x: [int(num) for num in x.split(" ")[0:-1]],line) 
               for line in map(lambda ln : ln.split(",")[0:-1],adsorption.splitlines() )]
absorption_strength = [map(lambda x: max(x)/sum(x) ,line) for line in absorption]
absorption_value = [map(lambda x: x.index(max(x)), line) for line in absorption]
a=np.array(absorption_value)
print(len(absorption_value))
"""
print((absorption_value[:][10]))
fromTop = map(lambda x: Counter((absorption_value[:][10])[:x])[0],range(1,len(L)))
fromBase = map(lambda x: Counter((absorption_value[:][10])[x+1:])[2],range(len(L-1)))
print(fromTop)
print(fromBase)
print([x == y for (x,y) in zip(fromTop,fromBase)].index(True))

points = []

for i in range(len(C)):
	
	fromTop = map (lambda x: Counter( absorption_value[i][:x] )[0],range(1,len(L)))
	fromBase = map(lambda x: Counter( absorption_value[i][x+1:])[2],range(len(L)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorption_value[i][crossover.index(True)] == 0 or absorption_value[i][crossover.index(True)] == 2 :
			points.append( (crossover.index(True),i) )
	

	fromTop = map (lambda x: Counter( absorption_value[i][:x] )[1],range(1,len(L)))
	fromBase = map(lambda x: Counter( absorption_value[i][x+1:])[2],range(len(L)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorption_value[i][crossover.index(True)] == 1 or absorption_value[i][crossover.index(True)] == 2 :
			points.append( (crossover.index(True),i) )

	fromTop = map (lambda x: Counter( absorption_value[i][:x] )[0],range(1,len(L)))
	fromBase = map(lambda x: Counter( absorption_value[i][x+1:])[1],range(len(L)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorption_value[i][crossover.index(True)] == 0 or absorption_value[i][crossover.index(True)] == 1 :
			points.append( (crossover.index(True),i) )

absorptionT = zip(*absorption_value)

for i in range(len(L)):
	
	fromTop = map (lambda x: Counter( absorptionT[i][:x] )[0],range(1,len(C)))
	fromBase = map(lambda x: Counter( absorptionT[i][x+1:])[2],range(len(C)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorptionT[i][crossover.index(True)] == 0 or absorptionT[i][crossover.index(True)] == 2 :
			points.append( (i,crossover.index(True)) )
	

	fromTop = map (lambda x: Counter( absorptionT[i][:x] )[1],range(1,len(C)))
	fromBase = map(lambda x: Counter( absorptionT[i][x+1:])[2],range(len(C)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorptionT[i][crossover.index(True)] == 1 or absorptionT[i][crossover.index(True)] == 2 :
			points.append( (i,crossover.index(True)) )

	fromTop = map (lambda x: Counter( absorptionT[i][:x] )[0],range(1,len(C)))
	fromBase = map(lambda x: Counter( absorptionT[i][x+1:])[1],range(len(C)-1))
	crossover = [x == y for (x,y) in zip(fromTop,fromBase)]
	if True in crossover:
		if absorptionT[i][crossover.index(True)] == 0 or absorptionT[i][crossover.index(True)] == 1 :
			points.append( (i,crossover.index(True)) )
		
"""
R=2*(gam_0+k_B*T*exc_s*np.log(1-G/exc_s))/(C_grid*N_a*k_B*T)
beta=exc_s*R*k_off/xi
R*=1e6
TH=np.pi - np.arcsin(L)
fig,ax=plt.subplots()
fig.canvas.draw()
ax.set_ylim(1,2e5)
ax.set_xlim(87,183)
ax.set_yticks([1,1e1,1e2,1e3,1e4,1e5])
ax.set_xticks([90,120,150,180])
m=ax.get_yticks()
ax.set_yticks(m)

ax2=ax.twinx()

ax2.set_ylim(1,2e5)
ax2.set_xlim(87,183)


for i in range(len(C)):
    for j in range(len(L)):
        ste = 2.0*(absorption_strength[i][j]-0.5)
        val = absorption_value[i][j]
        if val == 0:
           val = 1.0/3.0 #green stable
	   if absorption_strength[i][j] > 0.55:
	 	shape = '^'
	   else:
		shape = 'o'
        elif val == 1:
           val = 0.0 # red expand
	   if absorption_strength[i][j] > 0.55:
		shape = 'D'
	   else:
		shape = 'o'
        elif val == 2:
           val = 2.0/3.0 #blue unzip
	   if absorption_strength[i][j] > 0.55:
		shape = 'p'
	   else:
		shape = 'o'
	elif val == 3:
           val = 1.0/6.0
	   if absorption_strength[i][j] > 0.55:
		shape = '*'
	   else:
		shape = 'o'           
	if absorption_strength[i][j] < 0.551:
		ste = 0.0
	if j%2==0:
      		ax.semilogy([TH[j]*180/np.pi],[R[i]],marker=shape,ls='',color=col.hsv_to_rgb((val,ste,0.7)),ms=12,label='unstable')
		ax2.semilogy([TH[j]*180/np.pi],[R[i]],marker='',ls='')

"""
for point in points:
	ax.semilogy([TH[point[0]]*180/np.pi],[R[point[1]]],marker='o',color="black",ms=9,label="guide")

fp = open("guideEU/EU_guide_out.csv","r")
guide = fp.read()
fp.close()
guide = [float(num) for num in guide.splitlines()]

xes = [90 + x/10.0 for x in range(0,901)]
yes = [pow(10,x) for x in guide]
ax.semilogy(xes,yes,'--',color="Black",lw='5')


guide = map( lambda (x,y): (TH[x]*180/np.pi,R[y]),list(set(points)))
print (guide)
fp = open('EU_guide.txt','w')
for item in guide:
	fp.write("%f,%f\n" % (item[0], item[1]))
fp.close()
"""

ax.tick_params(labelsize=40)
ax.set_xlabel(r'Polar angle $\theta$ ($^{\circ}$)',fontsize=40)
ax.set_ylabel(r'$R^{*}$ ($\mu m$)',fontsize=40)  

ax2_t=m*exc_s*1e-6*k_off/xi
ax2.set_yticks(m)
labels = [item.get_text() for item in ax2.get_yticklabels()]
labels =  [r'1.2x10$^{-1}$']+[r'{:.1f}x10$^{}$'.format(ax2_t[i]/10**(i-1),i-1) for i in range(1,len(ax2_t))]
ax2.set_yticklabels(labels,rotation=0)
ax2.tick_params(labelsize=40)
ax2.set_ylabel(r'$\beta_{flow}$',fontsize=40)
print(ax2_t)

'''
ax.legend(loc='upper center',prop={'size':30},fancybox=True,ncol=2,bbox_to_anchor=(0.5,1.6))
'''
labels = [['^',1.0/3.0,1.0],['D',0.0,1.0],['p',2.0/3.0,1.0],['o',0.0,0.0]]
descriptions = ['stable', 'bilayer expansion', 'unzipping','no dominant behaviour']
proxies = [create_proxy(item) for item in labels]
ax.legend(proxies,descriptions,numpoints=1,markerscale=2.5,fontsize=35,loc='upper left',fancybox=True)


ax2.set_ylim(1,2e5)
fig.show()
raw_input("");


      
