Add files via upload

LB.py

@@ -0,0 +1,165 @@
+"""
+@author: EAK
+"""
+
+from matplotlib.ticker import MaxNLocator
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes, mark_inset
+from matplotlib import cm
+import matplotlib.pyplot as plt
+import numpy as np
+import matplotlib.ticker as ticker
+import os
+import mainantonpaar
+
+#### data read-in and preprocessing ##########################################
+
+# change to directory where data are stored
+dir='/Users/Emile/Documents/Graduate/AntonPaarData/Porifera/summary/shear data'
+os.chdir(dir)
+
+# first round LAOS
+fn1    = '02082020_sponge12_samp3_strainsweep_8mm_3mm 0'
+run1   = mainantonpaar.LAOS(fn1, r = 0.006)
+t1     = run1.Timep
+gam0_1 = run1.gam0
+sig0_1 = run1.sig0
+gamma1 = run1.Strain
+sigma1 = run1.Stress
+GpM1   = run1.Gpm
+GpL1   = run1.Gpl
+Gstar1 = run1.Gstarnorm
+n1     = run1.n
+# second round LAOS
+fn2    = '02082020_sponge12_samp3_strainsweep_8mm_3mm 1'
+run2   = mainantonpaar.LAOS(fn2, r = 0.006)
+t2     = run2.Timep
+sig0_2 = run2.sig0
+gamma2 = run2.Strain
+sigma2 = run2.Stress
+GpM2   = run2.Gpm
+GpL2   = run2.Gpl
+Gstar2 = run2.Gstarnorm
+n2     = run2.n
+# number of amplitudes
+Nf  = np.size(gam0_1)
+# amplitude array for looping
+Nf_ = np.arange(0, Nf, 1)
+# strain amplitudes, change to percent
+gam0 = [100*gam0_1[i] for i in Nf_]
+
+#### plotting ################################################################
+
+# color scheme
+cwsubsec = np.linspace(0.7, 0.1, Nf)   # first round, blues and reds
+clrscw   = [cm.pink(x) for x in cwsubsec]
+grsubsec = np.linspace(0.2, 0.7, Nf) # second round, transparent grays
+clrsgr   = [cm.gray(x) for x in grsubsec]
+# strain amplitude labels
+lbls = ['$\gamma_{0}=$ %.2f' % gam0[i] +'\%' for i in Nf_]
+
+# time series
+fig0, axs0 = plt.subplots(nrows = 1, ncols = 2, figsize = (22, 11))
+[axs0[0].plot(t1[i], gamma1[i], linewidth = 3, color = clrscw[i]) 
+ for i in Nf_[1::2]] # strain vs time, 1st run
+[axs0[0].plot(t2[i], gamma2[i], linewidth = 3, color = clrsgr[i]) 
+ for i in Nf_[1::2]] # strain vs time, 2nd run
+[axs0[1].plot(t1[i], sigma1[i]/1000, linewidth = 3, color = clrscw[i]) # kPa
+ for i in Nf_[1::2]] # stress vs time, 1st run
+[axs0[1].plot(t2[i], sigma2[i]/1000, linewidth = 3, color = clrsgr[i]) # kPa
+ for i in Nf_[1::2]] # stress vs time, 2nd run
+axs0[0].set_ylabel(r'$\gamma$ (\%)', fontsize = 48, labelpad = -30)
+axs0[1].set_ylabel(r'$\tau$ (kPa)', fontsize = 48, labelpad = -30)
+maxx =  np.max(np.abs(axs0[0].get_xlim()))
+maxy = [np.max(np.abs(axs0[i].get_ylim())) for i in [0,1]]
+[axs0[i].set_ylim([-maxy[i]+0.06*maxy[i],maxy[i]-0.06*maxy[i]]) for i in [0,1]]
+[axs0[j].set_xlabel(r'time (s)', fontsize = 48) for j in [0,1]]
+[axs0[j].tick_params(axis = 'both', which = 'major', labelsize = 38) 
+ for j in [0,1]]
+# insets
+axins = [inset_axes(axs0[i], width = '40%', height = '40%', loc = 1) 
+         for i in [0,1]]
+[axins[0].plot(t1[i], gamma1[i], lw = 3, c = clrscw[i]) for i in Nf_[1::2]]
+[axins[0].plot(t2[i], gamma2[i], lw = 3, c = clrsgr[i]) for i in Nf_[1::2]]
+[axins[i].tick_params(axis = 'both', labelsize = 28) for i in [0,1]]
+[axins[i].set_xlim(0.5*maxx, maxx-0.1*maxx) for i in [0,1]]
+[axins[i].set_ylim(-0.2*maxy[i], -0.001*maxy[i]) for i in [0,1]]
+[axins[1].plot(t1[i], sigma1[i]/1000, lw = 3, c = clrscw[i]) 
+ for i in Nf_[1::2]]
+[axins[1].plot(t2[i], sigma2[i]/1000, lw = 3, c = clrsgr[i]) 
+ for i in Nf_[1::2]]
+[mark_inset(axs0[i], axins[i], loc1 = 4, loc2 = 3, fc = 'none', ec = '0')
+ for i in [0,1]]
+fig0.subplots_adjust(wspace = 0.2)
+fig0.savefig('plots/'+fn1[:-2]+'_timeseries.svg', bbox_inches = 'tight', 
+            dpi = 1200)
+
+# Lissajous-Bowditch curves
+fig1, ax1 = plt.subplots(figsize = (12,12))
+x1 = np.asarray(gamma1, dtype = object)
+y1 = np.asarray(sigma1, dtype = object)
+x2 = np.asarray(gamma2, dtype = object)
+y2 = np.asarray(sigma2, dtype = object)
+# wrap-around so curves are totally closed
+x1 = [np.append(x1[i][:], x1[i][0]) for i in Nf_]
+y1 = [np.append(y1[i][:], y1[i][0]) for i in Nf_]
+x2 = [np.append(x2[i][:], x2[i][0]) for i in Nf_]
+y2 = [np.append(y2[i][:], y2[i][0]) for i in Nf_]
+# plot
+[ax1.plot(x2[i], y2[i]/1000, lw = 3, c = clrsgr[i]) for i in Nf_[1::2]]
+[ax1.plot(x1[i], y1[i]/1000, lw = 3, c = clrscw[i], label = lbls[i])
+ for i in Nf_[1::2]]
+# lines illustrating minimum and maximum moduli
+x  = np.array([-maxy[0]+0.06*maxy[0], 1.04*maxy[0]])
+ym = GpM1[-1]*x/100
+yl = GpL1[-1]*x/100
+ax1.plot(x, ym/1000, linestyle = '--', lw = 3, c = 'black', marker = None)
+ax1.plot(x, yl/1000, linestyle = '-.', lw = 3, c = 'black', marker = None)
+ax1.tick_params(axis = 'both', which = 'major', labelsize = 48)
+ax1.set_xlabel(r'$\gamma$ (\%)',  fontsize = 58)
+ax1.set_ylabel(r'$\tau$ (kPa)', fontsize = 58, labelpad = -30)
+ax1.set_xlim([-1.008*maxy[0], 1.008*maxy[0]])
+ax1.set_ylim([-1.008*maxy[1], 1.008*maxy[1]])
+ax1.xaxis.set_major_locator(ticker.MaxNLocator(5))
+ax1.yaxis.set_major_locator(ticker.MaxNLocator(6))
+hndls, lbls = ax1.get_legend_handles_labels()
+fig1.legend(hndls, lbls, loc = 'center right', bbox_to_anchor = (1.32,0.32), 
+            ncol = 1, fancybox = True, shadow = True, fontsize = 38)
+# small strain inset
+axins1 = inset_axes(ax1, width = '38%', height = '38%', loc = 4)
+[axins1.plot(x2[i], y2[i]/1000, lw = 3, c = clrsgr[i]) for i in Nf_[1::2]]
+[axins1.plot(x1[i], y1[i]/1000, lw = 3, c = clrscw[i]) for i in Nf_[1::2]]
+axins1.set_xlim(-0.077*maxy[0], 0.077*maxy[0])
+axins1.set_ylim(-0.077*maxy[1], 0.077*maxy[1])
+axins1.xaxis.set_major_locator(ticker.MaxNLocator(3))
+axins1.yaxis.set_major_locator(ticker.MaxNLocator(3))
+axins1.xaxis.set_ticks_position('top') 
+axins1.tick_params(axis = 'both', labelsize = 28)
+mark_inset(ax1, axins1, loc1 = 2, loc2 = 1, fc = 'none', ec = '0')
+# large strain
+axins2 = inset_axes(ax1, width = '38%', height = '38%', loc = 2)
+[axins2.plot(x2[i], y2[i]/1000, lw = 3, c = clrsgr[i]) for i in Nf_[1::2]]
+[axins2.plot(x1[i], y1[i]/1000, lw = 3, c = clrscw[i]) for i in Nf_[1::2]]
+axins2.set_xlim(0.21*maxy[0], maxy[0]-0.077*maxy[0])
+axins2.set_ylim(0.21*maxy[1], maxy[1]-0.077*maxy[1])
+axins2.xaxis.set_major_locator(ticker.MaxNLocator(2))
+axins2.yaxis.set_major_locator(ticker.MaxNLocator(2))
+axins2.yaxis.set_ticks_position('right') 
+axins2.tick_params(axis = 'both', labelsize = 28)
+mark_inset(ax1, axins2, loc1 = 4, loc2 = 1, fc = 'none', ec = '0')
+# save
+fig1.savefig('plots/'+fn1[:-2]+'_LB.svg', bbox_inches = 'tight', dpi = 1200)
+
+# power spectra
+fig2, ax2 = plt.subplots(figsize = (12,12))
+[ax2.plot(n1[i], Gstar1[i], linewidth = 3, color = clrscw[i]) 
+ for i in Nf_[1::2]]
+[ax2.plot(n2[i], Gstar2[i], linewidth = 3, color = clrsgr[i]) 
+ for i in Nf_[1::2]]
+ax2.set_xlabel('$n$', fontsize = 42)
+ax2.set_ylabel('$|G^{*}_{n}|$/$|G^{*}_{1}|$', fontsize = 42)
+ax2.xaxis.set_major_locator(MaxNLocator(integer = True))     
+ax2.tick_params(axis = 'both', which = 'major', labelsize = 48)
+ax2.set_xlim([0, 12])
+ax2.set_ylim([0.0002, 1])
+ax2.set_yscale('log')
+fig2.savefig('plots/'+fn1[:-2]+'_power.svg', bbox_inches = 'tight', dpi = 1200)

frequencysweepsummary.py

@@ -0,0 +1,184 @@
+"""
+@author: EAK
+"""
+
+from matplotlib import cm
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import pandas as pd
+import numpy as np
+import os
+import rheofuncs
+
+#### Data Read-In ############################################################
+
+# directory
+dir = '/Users/Emile/Documents/Graduate/AntonPaarData/Porifera/summary/'+\
+      'frequency data'
+os.chdir(dir)
+# file
+datstr = 'freqsweepfullsummary'
+dat    = pd.read_csv(datstr + '.csv').values 
+n      = np.arange(0, np.int64(np.size(dat, 1)/8), 1) # number of species array
+# columns
+phi   = [dat[0, 8*i+1] for i in n] # mass porosity
+phi_e = [dat[0, 8*i+2] for i in n] # mass porosity uncertainty
+w     = [dat[:, 8*i+3].astype('float64') for i in n] # frequency, rad/s
+Gp    = [dat[:, 8*i+4].astype('float64') for i in n] # storage modulus, Pa
+Gpp   = [dat[:, 8*i+5].astype('float64') for i in n] # loss modulus, Pa
+Gp_e  = [dat[:, 8*i+6].astype('float64') for i in n] # stored uncertainty, Pa
+Gpp_e = [dat[:, 8*i+7].astype('float64') for i in n] # loss uncertainty, Pa
+
+#### Calculations ############################################################
+
+# loss tangent
+tand = [Gpp[i]/Gp[i] for i in n]
+# propagate its uncertainty from G' and G" uncertainties
+tand_e = [np.sqrt((Gp_e[i]*Gpp[i]/Gp[i]**2)**2+(Gpp_e[i]/Gp[i])**2) for i in n]
+# linear viscoelastic fit parameters file
+fitstr, mdlstr = 'fitparams'+datstr[9:13], '_Zener' # choose which one
+fit = pd.read_csv(fitstr + mdlstr + '.csv').values
+params = [fit[:, i + 1].astype('float64') for i in n]
+# get best fit values by inputting fit parameters into the model
+Gpmodel  = getattr(rheofuncs, mdlstr[1:] + 'Gp')
+Gppmodel = getattr(rheofuncs, mdlstr[1:] + 'Gpp')
+Gpbest   = [ Gpmodel(w[i], *params[i][:-1]) for i in n]
+Gppbest  = [Gppmodel(w[i], *params[i][:-1]) for i in n]
+tandbest = [Gppbest[i]/Gpbest[i] for i in n]
+# power law exponent of G'
+m     = [np.polyfit(np.log10(w[i]), np.log10(Gp[i]), 1, cov = True) for i in n]
+alpha = [m[i][0][0] for i in n]
+alpha_e = [np.sqrt(np.diag(m[i][1])[0]) for i in n]
+
+#### Plotting ################################################################
+# font
+mpl.rc('font', family = 'Courier New')
+mpl.rc('text', usetex = True)
+# full species list
+spp   = ['Axinella polycapella','Callyspongia sp.','Cinachyrella apion',
+         'Cliona celata','Tectitethya keyensis','Tethya aurantium', 
+         'Aplysina fulva','Igernella notabilis']
+spp   = [r'\textit{%s}' % spp[i] for i in n[:-1]] # italicize
+spp.insert(len(spp), 'synthetic sponge')  # add synthetic sponge label
+# colors
+spixclrs   = [cm.cool(x) for x in np.linspace(0, 1, 6)]
+nospixclrs = [cm.spring(x) for x in np.linspace(0.6, 0.89, 2)]
+clrs       = [spixclrs[5], spixclrs[4], spixclrs[3], spixclrs[2], spixclrs[1],
+              spixclrs[0], nospixclrs[0], nospixclrs[1]]
+clrs.insert(len(spp), 'xkcd:dried blood')
+# markers
+mrkrs = ['s','s', 's','s','s','s','o','o','*']
+
+# Gp vs omega
+fig0, axs0 = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 10))
+[axs0.errorbar(w[i], Gp[i], yerr = Gp_e[i], ms = 20, fmt = mrkrs[i], mew = 3,
+               mfc = clrs[i], mec = 'black', ecolor = 'dimgray', 
+               label = spp[i], capsize = 6) for i in n]
+[axs0.plot(w[i], Gpbest[i], ls = '--', lw = 3, c = clrs[i]) for i in n]
+axs0.set_xscale('log')
+axs0.set_yscale('log')
+axs0.set_ylim([35000, 2500000])
+axs0.tick_params(axis = 'both', which = 'major', labelsize = 44)
+axs0.set_xlabel('$\omega$ (rad/s)', fontsize = 44, labelpad = -6)
+axs0.set_ylabel('$G^{\prime}$ (Pa)', fontsize = 44, labelpad = -10)
+hndls, lbls = axs0.get_legend_handles_labels() # create legend
+fig0.legend(hndls, lbls, loc = 'center right', bbox_to_anchor = (-0.075,0.5), 
+            ncol = 1, fancybox = True, shadow = True, fontsize = 40)
+fig0.savefig('plots/Gpvsomega'+datstr[9:13]+mdlstr+'.svg', dpi = 1200, 
+             bbox_inches = 'tight')
+
+# Gpp vs omega
+fig1, axs1 = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 10))
+[axs1.errorbar(w[i], Gpp[i], yerr = Gpp_e[i], ms = 18, fmt = mrkrs[i], 
+               mew = 3, mfc = clrs[i], mec = 'black', ecolor = 'dimgray', 
+               label = spp[i], capsize = 6) for i in n]
+[axs1.plot(w[i], Gppbest[i], ls = '--', lw = 3, c = clrs[i]) for i in n]
+axs1.plot(w[8], Gppbest[8], ls = '--', lw = 3, c = clrs[8])
+axs1.set_xscale('log')
+axs1.set_yscale('log')
+axs1.set_ylim([2000, 400000])
+axs1.tick_params(axis = 'both', which = 'major', labelsize = 44)
+axs1.set_xlabel('$\omega$ (rad/s)', fontsize = 44, labelpad = -6)
+axs1.set_ylabel('$G^{\prime \prime}$ (Pa)', fontsize = 44, labelpad = -10)
+fig1.savefig('plots/Gppvsomega'+datstr[9:13]+mdlstr+'.svg', dpi = 1200, 
+             bbox_inches = 'tight')
+
+# loss tangent vs omega
+fig2, axs2 = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 10))
+[axs2.errorbar(w[i], tand[i], yerr = tand_e[i], ms = 18, fmt = mrkrs[i], 
+               mew = 3, mfc = clrs[i], mec = 'black', ecolor = 'dimgray', 
+               label = spp[i], capsize = 6) for i in n]
+[axs2.plot(w[i], tandbest[i], ls = '--', lw = 3, c = clrs[i]) for i in n]
+max2, min2 = np.max(axs2.get_ylim()), np.min(axs2.get_ylim())
+axs2.set_xscale('log')
+axs2.set_yscale('log')
+axs2.set_ylim([0.03, 0.3])
+axs2.tick_params(axis = 'both', which = 'major', labelsize = 44)
+axs2.tick_params(axis = 'y', which = 'minor', labelsize = 28)
+axs2.set_xlabel('$\omega$ (rad/s)', fontsize = 44, labelpad = -6)
+axs2.set_ylabel(r'$\tan\delta$', fontsize = 44, labelpad = -10)
+fig2.savefig('plots/tandelvsomega'+datstr[9:13]+mdlstr+'.svg', dpi = 1200, 
+             bbox_inches = 'tight')
+
+# Gpp vs Gp
+fig3, axs3 = plt.subplots(nrows = 1, ncols = 1, figsize = (10, 10))
+[axs3.errorbar(Gp[i][6], Gpp[i][6], xerr = Gp_e[i][6], yerr = Gpp_e[i][6], 
+               ms = 18, fmt = mrkrs[i], mew = 3, mfc = clrs[i], 
+               mec = 'black', ecolor = 'dimgray', capsize = 6) for i in n]
+x = [10**4, 10**5, 10**6, 10**7]
+axs3.plot(x,[0.06*x[i] for i in [0,1,2,3]], lw = 3, ls = '-.', c = 'silver',
+          label = '$G^{\prime \prime}=0.06G^{\prime}$')
+axs3.plot(x,[0.18*x[i] for i in [0,1,2,3]], lw = 3, ls = ':', c = 'silver',
+          label = '$G^{\prime \prime}=0.18G^{\prime}$')
+axs3.set_xscale('log')
+axs3.set_yscale('log')
+axs3.set_xlim([4*10**4, 2*10**6])
+axs3.set_ylim([2*10**3, 3*10**5])
+axs3.tick_params(axis = 'both', which = 'major', labelsize = 44)
+axs3.set_xlabel('$G^{\prime}(\omega=\pi)$ (Pa)', fontsize = 44, labelpad = -6)
+axs3.set_ylabel('$G^{\prime \prime}(\omega=\pi)$ (Pa)', fontsize = 44)
+hndls, lbls = axs3.get_legend_handles_labels() # create legend
+axs3.legend(hndls, lbls, loc = 'upper left', fancybox = True, shadow = True, 
+            fontsize = 40)
+fig3.savefig('plots/GppvsGp.svg', dpi = 1200, bbox_inches = 'tight')
+
+# alpha vs. specimen
+fig4, axs4 = plt.subplots(nrows = 1, ncols = 1, figsize = (10, 10))
+X = np.argsort(-np.asarray(alpha))
+[axs4.bar(n[i], alpha[X[i]], 0.65, yerr = alpha_e[X[i]], color = clrs[X[i]], 
+          capsize = 6, ecolor = 'dimgray') for i in n]
+axs4.set_xticks(n)
+axs4.set_xticklabels([])
+axs4.set_ylim(0, 0.076)
+axs4.set_yticks([0.01, 0.03, 0.05, 0.07])
+axs4.tick_params(axis = 'y', which = 'major', labelsize = 44)
+axs4.set_ylabel(r'$\alpha$', fontsize = 44, labelpad = 20, rotation = 0)
+fig4.savefig('plots/alphabarplot.svg', dpi = 300, bbox_inches = 'tight')
+
+# porosity vs. specimen
+fig5, axs5 = plt.subplots(nrows = 1, ncols = 1, figsize = (10, 10))
+X = np.argsort(-np.asarray(phi))
+[axs5.bar(n[i], phi[X[i]], 0.65, yerr = phi_e[X[i]], color = clrs[X[i]], 
+          capsize = 6, ecolor = 'dimgray') for i in n]
+spp_sort = [spp[X[i]] for i in n]
+axs5.set_xticks(n)
+axs5.set_xticklabels(spp_sort, fontdict = {'fontsize': 33}, rotation = 90)
+axs5.set_ylim(0, 0.96)
+axs5.set_yticks([0.2, 0.4, 0.6, 0.8])
+axs5.tick_params(axis = 'y', which = 'major', labelsize = 44)
+axs5.set_ylabel(r'$\phi_{w}$', fontsize = 44, labelpad = 20, rotation = 0)
+fig5.savefig('plots/phibarplot.svg', dpi = 300, bbox_inches = 'tight')
+
+# spicules bar plot
+fig6, axs6 = plt.subplots(nrows = 1, ncols = 1, figsize = (10, 10))
+spix   = [43900, 1540, 4020, 2000000, 2780, 5960, 0, 0, 0]  # number per cm^3
+spix_e = [20000,  400, 2000, 1750000, 2000, 2000, 0, 0, 0]  # uncertainty
+[axs6.bar(n[i], spix[X[i]], 0.65, yerr = spix_e[X[i]], color = clrs[X[i]], 
+          capsize = 6, ecolor = 'dimgray', rasterized = True) for i in n]
+axs6.set_xticks(n)
+axs6.set_xticklabels([])
+axs6.set_yticks([1000, 10000, 100000, 1000000])
+axs6.set_yscale('log')
+axs6.tick_params(axis = 'y', which = 'major', labelsize = 44)
+axs6.set_ylabel('$N_{spicules}$/cm$^{3}$', fontsize = 44, labelpad = 10)
+fig6.savefig('plots/spiculesbarplot.svg', dpi = 300)