#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 31 19:43:42 2017

@author: Alain
"""

import matplotlib.pyplot as plt
import utilitaire as util

""" Titrage conductimétrique """


fichier = open("MgCl2HCl.txt")
lignes = fichier.readlines()
fichier.close()

V=[]
PH=[]
G=[]
first = True
for l in lignes :
    if first : 
        first = False
    else : 
        if len(l) >0 : 
            ls = l.split(";")
            V.append(float(ls[0]))
            PH.append(float(ls[1]))
            G.append(float(ls[2]))

fig, ax1 = plt.subplots()

ax1.plot(V, PH, "rx",markersize=3)
ax1.set_xlabel("V (mL)")
ax1.set_ylabel("pH", color="r")
ax1.tick_params("y", colors="r")
ax1.set_ylim(0,14)

ax2 = ax1.twinx()
ax2.plot(V, G, "b+")
ax2.set_ylabel("conductivité (mS)", color="b")
ax2.tick_params("y", colors="b")
ax2.set_ylim(0,10)
ax2.set_xlim(0,20)
ax2.set_xticks([i for i in range(0,22,2)])

plt.title("Titrage d'un mélange MgCl2 - HCl par NaOH")
plt.show()

V0 = 50
Gc=[]
for i in range(len(V)):
    Gc.append(G[i]*(V0+V[i])/V0)
    
plt.plot(V,G,"r+")
plt.plot(V,Gc,"bx")
plt.show()

a1,b1,r1 = util.regrelinEntre(V,Gc,0,8)
print(a1,b1)
a2,b2,r2 = util.regrelinEntre(V,Gc,10,12)
print(a2,b2)
a3,b3,r3 = util.regrelinEntre(V,Gc,15,20)
print(a3,b3)

Veq1 = (b1 - b2)/(a2 - a1)
Veq2 = (b2 - b3)/(a3 - a2)
print("Premier volume équivalent : ", round(Veq1,1), " mL")
print("Second volume équivalent : ", round(Veq2,1), " mL")


""" 

                    Couplage AMX 

"""

nuSpectro = 100 # fréquence du spectro en MHz
deltaA = 5 # déplacement chimique de A
nbM = 1 # nombre de protons de type M
JAM = 12 # constante de couplage AM
nbX = 2 # nombre de protons de type X
JAX = 8 # constante de couplage AX

def lorentz(delta, intensite, deltaRef, deltappm):
    return intensite / (1 + ((delta - deltaRef) / (5* deltappm) )**2)

def factorielle(n):
    if n == 0 : 
        return 1
    else : 
        return n*factorielle(n-1)

def combinaison(n,p):
    return factorielle(n)/(factorielle(p)*factorielle(n-p))

def genereFils(deltaRef, intensiteRef, nbFils, J):
    Delta = []
    I = []
    d = J/nuSpectro
    nI = 0
    for i in range(0,nbFils):
        nI = nI + combinaison(nbFils-1, i)
    if nbFils %2 == 0 : 
        d0 = (d/2) + ((nbFils/2)-1)*d
    else :
        d0 = ((nbFils-1)/2)*d
    for i in range(0, nbFils):
        Delta.append(deltaRef - d0 + i*d)
        I.append(intensiteRef*combinaison(nbFils-1, i))
    return Delta, I

""" programme principal """   
 
deltaPic = []
intensitePic = []
dp, ip = genereFils(deltaA, 1, nbM+1, JAM)
for i in range(len(dp)):
    dp2, ip2 = genereFils(dp[i], ip[i], nbX+1, JAX)
    for j in range(len(dp2)):
        deltaPic.append(dp2[j])
        intensitePic.append(ip2[j])

deltaMin = deltaA-0.5
deltaMax = deltaA+0.5
deltappm = 0.001
delta = deltaMin
D = []
I = []
while delta < deltaMax : 
    D.append(delta)
    intensite = 0
    for i in range(len(deltaPic)):
        intensite = intensite + lorentz(delta, intensitePic[i], deltaPic[i], deltappm)   
    I.append(intensite)
    delta = delta + deltappm

plt.plot(D, I)
plt.xlim(deltaMax,deltaMin)
plt.show()
