# -*- coding: utf-8 -*-
"""
Created on 04.11.2022

script to evaluate the viscosity of CA-acetone solutions

@author: csag4207
"""

import sys
sys.path.append('C:/Users/csag4207/Documents/forschung/python/scripts') # this has to be imported until the pythonpath is permanently changed to find the module
# from os import listdir
import matplotlib.pyplot as plt
import numpy as np
# from scipy import optimize
from datetime import datetime, timedelta

# from lmfit import Model, Parameters #, Parameter, report_fit
# # from PIL import Image
# import xlrd
import csv

#%%import the rheometer data

file_path1 = "C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/data/di_h2o.csv"

#convert file from unicode to utf-8 first via Editor - if time make suitable for Unicode
ws = []
buff = open(file_path1, 'r', encoding='utf-16')
#buff = unicodecsv.reader(file_path1, encoding = 'unicode')
csvfile = csv.reader(buff, delimiter='\t')
for row in csvfile:
    ws.append(row)
buff.close()
del buff, file_path1

# find different tests, number of measured intervals and header of each interval performed
testlines = [0]
inter = []
head = []
cuting = []
for i in np.arange(0, len(ws)):
    if 'Test:' in ws[i]:
        testlines.append(i)
    if 'Interval and data points:' in ws[i]: #anton paar changed the parameter names with the update from 01.07.19
        inter.append(i)
    if '[Pa]' in ws[i]:
        head.append(i)
    if ws[i] == []:
        cuting.append(i)

tests = np.unique([ws[i][1] for i in testlines])
interval = np.unique([ws[i][1] for i in inter])
header = np.squeeze(np.unique(np.array([ws[i] for i in head]), axis=0))
measurand = ws[7][1:-1]
# add a final line to the cuting array - update problem from 01.07.19
cuting.append(ws.__len__()) # cuting does not realise the last line

del testlines, inter

#make sure there has only been one type of measurement
if header.ndim > 1:
    print('Take care: more than one data matrix with different data is present')

tests = np.delete(tests, np.where(tests=='20221021_CA_viscosity_FW')[0])

print('possible data available:')
print(list(tests))
#%% delete faulty or unimportant data

use=[]

for i in tests:
    if 'h2o' in i:
        if '_new' in i:
            use.append(i)
    if 'h2o' not in i:
        if '_2' in i:
            use.append(i)

use.remove('8pc_CA_2 1')

#%%

#fill with information from the measurement script
indice = {}
for i in use:
    indice['_'.join(i.split('_')[:-1])] = [n for n, ws in enumerate(ws) if i in ws]

stop_indice = {}
for i in indice:
    stop_indice[i] = cuting[np.abs(np.array(cuting)+1-indice[i][0]).argmin()+5] # Anton paar has an empty line per interval after the update

# extract the data from the ws_matrix and throw all the lines that are of non-interest
data = {}
pop = {}
listing = {}
for i in indice.keys():
    data[i] = [ws[j] for j in np.arange(indice[i][0], stop_indice[i]) if ws[j] != []]
    listing[i]=[]
    for j in np.arange(0, data[i].__len__()):
        if 'Interval and data points:' in data[i][j] or 'Test:' in data[i][j]:
            listing[i].append(j)
        if 'Result:' in data[i][j] or 'Fundamental' in data[i][j]:
            listing[i].append(j)
        if '[1/s]' in data[i][j] or 'Interval data:' in data[i][j]:
            listing[i].append(j)
        if ''==data[i][j][1]:
            listing[i].append(j)
        

# delete the non-data from the dictionary

for i in data:
    data[i] = np.delete(np.array(data[i]), np.array(listing[i]))
    data[i] = [np.concatenate((np.array(j[1:3]).reshape((1, 2)).astype(float), np.array(j[4:8]).reshape((1, 4)).astype(float), np.array(j[10]).reshape((1, 1)).astype(float), np.array(j[12:14]).reshape((1, 2)).astype(float),  np.array(datetime.strptime(j[14], '%H:%M:%S')).reshape((1,1)), np.array(datetime.strptime(j[15], '%d/%m/%Y %X')).reshape((1,1))), axis=1) for j in data[i]]
    data[i] = np.concatenate(data[i], axis=0)

##count the ones in the first coloumn to find the different intervals
## the second interval is a frequency sweep, the fourth (linear from 1:20%) and fifth (linear from 20-500%) interval is a amplitude sweep, interval 1,3,6 are for 10Hz 0.3% gel equilibrium / background measurements (Chen2017)
ones={}
for i in data:
    ones[i] = []
    for j in np.arange(0, len(data[i])):
        if data[i][j, 0] == 1:
            ones[i].append(j)


#calculate a mean over interval 1 to have a strength approx. comparison
stats={}
ex = 20
for i in data:
    stats[i]={}
    stats[i]['mean1'] = np.mean(data[i][ones[i][0]+ex:ones[i][1]-1, :9], axis=0)
    stats[i]['std1'] = np.std(data[i][ones[i][0]+ex:ones[i][1]-1, :9].astype('float64'), axis=0)
    stats[i]['mean2'] = np.mean(data[i][ones[i][2]+ex:-1, :9], axis=0)
    stats[i]['std2'] = np.std(data[i][ones[i][2]+ex:-1, :9].astype('float64'), axis=0)


del ws, csvfile

#%%plot the viscosity values in comparison

selection = ['di_h2o', '2pc_CA', '4pc_CA', '6pc_CA', '8pc_CA', '10pc_CA']

plt.figure()
for i in selection:
    plt.bar(i, stats[i]['mean1'][6], yerr=stats[i]['std1'][6])
# plt.yscale('log')
plt.show()

#%% for the manuscript
for i in selection:
    print('visc /mPas ', stats[i]['mean1'][6], ' +- ', stats[i]['std1'][6])

#%% plot the important data

fig, axs = plt.subplots(2,1)

for i in data:
    if 'alginat' in i and '7.5' in i: 
        axs[0].plot(data[i][:,9], data[i][:,2], label = i)
axs[0].legend()
axs[0].set_yscale('log')

for i in data:
    if 'xanthan' in i and '7.5' in i:
        axs[1].plot(data[i][:,9], data[i][:,2], label = i)
axs[1].legend()
axs[1].set_yscale('log')

plt.xlabel('time /s')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.show()

#%% just plot the xanthan on an additional plot

plt.figure()
for i in data:
    if 'xanthan' in i and '7.5' in i:
        plt.plot(data[i][:,9], data[i][:,2], label = i)
plt.legend()
plt.yscale('log')

plt.xlabel('time /s')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.show()

#%% better use a boxplot to plot the stats

fig, axs = plt.subplots(1,2, figsize=(16, 10))
plt.rcParams.update({'font.size': 20})

cmap = plt.get_cmap('tab10')
colors = [cmap(i) for i in np.linspace(0.2, 0.8, 2)]

alg_list = []
xan_list = []

for i in data:
    if 'alginat' in i and '7.5' in i:
        alg_list.append([i.split('_')[0].replace('pc', ''), stats[i]['mean1'][2]])

alg_list = np.array(alg_list).astype('float')
alg_list = alg_list[alg_list[:,0].argsort()]


axs[0].bar(alg_list[:,0].astype('str'), alg_list[:,1], color=colors[0])
axs[0].tick_params(axis='x', labelrotation = 60)
axs[0].set_title('alginate')
axs[0].set_ylabel('storage modulus (10 Hz) /Pa')
axs[0].set_xlabel('glyoxal /%')
axs[0].set_ylim([0,16700])

for i in data:
    if 'xanthan' in i and '7.5' in i:
        xan_list.append([i.split('_')[0].replace('pc', ''), stats[i]['mean1'][2],  stats[i]['mean1'][3], stats[i]['mean1'][4], ])

xan_list = np.array(xan_list).astype('float')
xan_list = xan_list[xan_list[:,0].argsort()]

axs[1].bar(xan_list[:,0].astype('str'), xan_list[:,1], color=colors[1])
axs[1].tick_params(axis='x', labelrotation = 60)
axs[1].set_title('xanthan')
axs[1].set_xlabel('glyoxal /%')
axs[1].set_ylim([0,16700])
# plt.show()

throwing = np.array([-1, -3, 1])

xan_list[0] = np.mean([xan_list[1], xan_list[0]], axis=0)
xan_list[-4] = np.mean([xan_list[-4], xan_list[-3]], axis=0)
xan_list[-2] = np.mean([xan_list[-2], xan_list[-1]], axis=0)

xan_list = np.delete(xan_list, throwing, axis=0)

# plt.savefig('C:/Users/csag4207/Documents/forschung/konferenz/2021.10_nantes/pic/rheometry_both.png', dpi=500)


# %% change what is on the plot

fig, ax1 = plt.subplots()

ax1.bar(xan_list[:,0].astype('str'), xan_list[:,1], yerr=0.05*xan_list[:,1], alpha=0.5, color='red', capsize=5)
ax1.bar(xan_list[:,0].astype('str'), xan_list[:,2], yerr=0.05*xan_list[:,2], alpha=0.5, color='c', capsize=3)

ax1.tick_params(axis='x', labelrotation = 60)
ax1.set_title('xanthan')
ax1.set_xlabel('glyoxal /%')
ax1.set_ylim([0,18000])
ax1.set_ylabel('storage and loss modulus /Pa')


ax2 = ax1.twinx()

ax2.scatter(xan_list[:,0].astype('str'), xan_list[:,3], color='b')
ax2.set_ylabel('loss factor')


#%% check the frequency and amplitude response of the gels

fig, axs = plt.subplots(1, 2)

for i in data:
    if 'alginat' in i and '7.5' in i: 
        axs[0].plot(data[i][120:150,1], np.divide(data[i][120:150,2], data[i][118,2]), label = i)
        axs[0].plot(data[i][120:150,1], np.divide(data[i][120:150,3], data[i][118,3]), linestyle=':', label = i, color = 'b')
        axs[1].plot(data[i][120:150,1], data[i][120:150,4], label = i + '_loss_factor')
# axs[0].legend()
axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].set_ylim(0.01,7)

axs[1].legend()
axs[1].set_ylim(0.1,10)
axs[1].set_xscale('log')


plt.xlabel('time /s')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.show()

#%%


fig, axs = plt.subplots(1, 2)

for i in data:
    if 'xanthan' in i and '7.5' in i: 
        axs[0].plot(data[i][120:150,1], np.divide(data[i][120:150,2], data[i][118,2]), label = i)
        axs[0].plot(data[i][120:150,1], np.divide(data[i][120:150,3], data[i][118,3]), linestyle=':', label = i, color = 'b')
        axs[1].plot(data[i][120:150,1], data[i][120:150,4], label = i + '_loss_factor')
# axs[0].legend()
axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].set_ylim(0.1,3)

axs[1].legend()
axs[1].set_ylim(0,0.7)
axs[1].set_xscale('log')


plt.xlabel('time /s')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.show()

# no relative peak shift of the loss factor (thouch seems)
# no clear difference between high-low glyoxal samples

#%% amplitude test for alginate

fig, axs = plt.subplots(1, 2)

for i in data:
    if 'alginat' in i and '7.5' in i: 
        axs[0].plot(data[i][270:281,5], np.divide(data[i][270:281,2], data[i][118,2]), label = i)
        axs[0].plot(data[i][270:281,5], np.divide(data[i][270:281,3], data[i][118,3]), linestyle=':', label = i, color = 'b')
        axs[1].plot(data[i][270:281,5], data[i][270:281,4], label = i + '_loss_factor')
# axs[0].legend()
axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].set_ylim(0.01,7)

axs[1].legend()
axs[1].set_ylim(0.1,10)
axs[1].set_xscale('log')


plt.xlabel('amplitude /%')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.title('amplitude sweep')
plt.show()

#%% amplitude test for xanthan

fig, axs = plt.subplots(1, 2)

for i in data:
    if 'xanthan' in i and '7.5' in i: 
        axs[0].plot(data[i][270:281,5], np.divide(data[i][270:281,2], data[i][118,2]), label = i)
        axs[0].plot(data[i][270:281,5], np.divide(data[i][270:281,3], data[i][118,3]), linestyle=':', label = i, color = 'b')
        axs[1].plot(data[i][270:281,5], data[i][270:281,4], label = i + '_loss_factor')
# axs[0].legend()
axs[0].set_yscale('log')
axs[0].set_xscale('log')
axs[0].set_ylim(0.01,7)

axs[1].legend()
axs[1].set_ylim(0.1,10)
axs[1].set_xscale('log')


plt.xlabel('amplitude /%')
# plt.vlines(3331, 0, 0.35, color='r')
plt.ylabel('storage modulus /Pa')
plt.tight_layout()
plt.title('amplitude sweep')
plt.show()

# 0.75 and 1pc sample also showed low storage modulus/ 0.5 showed high (data artefact/sample preparation?)

#%% from xan-list[:,1] calculate the D_p between the cross-links
# from davidovich 2010 : G = 1/N * c/M_w * RT
# m_w = monomer molar mass
# R gas constant, T is temperature
# G shear/elastic modulus, which is similar to the storage modulus for low delta (loss factor) # https://www.researchgate.net/post/Is-there-a-relationship-between-Storage-modulus-and-elastic-modulus-for-a-solid
# otherwise calulate elastic_modulus = storage_modulus / cos(delta) 

c_xan = 0.03 # g/ml or kg/l
temp = 273 + 25 # deg C
m_w = 177.3 # g/mol from 2021.05/eval_hplc.py from the calcuation of the glyoxal to monomer ratio
r_gas = 8.31446261815324 # J⋅K−1⋅mol−1

def dp(x):
    return 1/x*c_xan/m_w*r_gas*temp

corrected_list = xan_list[:,1]/np.cos(xan_list[:,-1])

dp_glyo = np.fromiter(map(dp, corrected_list), dtype='float')*1000**2 # the 1000**2 comes from the calculation (check the units)

del_dp_glyo = np.sqrt(0.05**2 + (0.001/3)**2) # inclusion of the weighing error in case the paste has been produced incorrectly. However the 5% error remains unchained: assume a final error of 5%

#%%replot for the text

# selection = ['0pc_xanthan', '1pc_xanthan', '2pc_xanthan', '6pc_xanthan', '12pc_xanthan']

plt.rcParams.update({'font.size': 16})

cmap = plt.get_cmap('inferno') # prior it was vidris
colors = [cmap(i) for i in np.linspace(0.1, 0.8, 3)]

linest = [ (0, (3, 5, 1, 5)), ':', '-.', (0, (5, 5)), '--']

fig, ax1 = plt.subplots(figsize=(15, 10))

ax1.bar(xan_list[:,0].astype('str'), xan_list[:,1], yerr=0.05*xan_list[:,1], alpha=0.5, color=colors[0], capsize=5)
ax1.bar(xan_list[:,0].astype('str'), xan_list[:,2], yerr=0.05*xan_list[:,2], alpha=0.5, color=colors[1], capsize=3)

# ax1.plot(xan_list[:,0].astype('str'), xan_list[:,3], color='b') # somehow the axis does not fit anymore


ax1.tick_params(axis='x', labelrotation = 60)
ax1.set_title('xanthan')
ax1.set_xlabel('glyoxal /%')
ax1.set_ylim([0,18000])
ax1.set_ylabel('storage and loss modulus /Pa')


ax2 = ax1.twinx()
ax2.plot(xan_list[:,0].astype('str'), dp_glyo, color = colors[-1])
ax2.scatter(xan_list[:,0].astype('str'), np.multiply(xan_list[:,3], 1000), marker= 'o', s=80, color ='r')
ax2.set_ylabel('$d_p$ between cross-links \n as $10^3*$loss factor $\delta$')


plt.savefig('C:/Users/csag4207/Documents/forschung/arbeit/2022/Wundauflage/pic/rheometry.png', dpi=600)
# plt.savefig('C:/Users/csag4207/Documents/forschung/konferenz/2021.10_nantes/pic/rheometry_xanthan.png', dpi=500)


