Revisions

doubleegrip.py

@@ -8,6 +8,7 @@
 import agedepth
 from matplotlib import pyplot as plt
 import mcfab as mc
+import specfabfuns as sff
 with open('path2dSGdt10.pkl', 'rb') as f:
     path2d = pickle.load(f)
 
@@ -59,8 +60,8 @@
 for i in tqdm(range(len(paths))):
     paths[i].Temperature()
     paths2[i].Temperature()
-    paths[i].fabric_mc(npoints)
-    paths2[i].fabric_mc(npoints,x='Reduced')
+    paths[i].fabric_sf(L)
+    paths2[i].fabric_sf(L,x='Reduced')
     depths[i] = paths[i].d[-1]
 
 
@@ -96,7 +97,7 @@
 import cartopy.crs as ccrs
 import cartopy.mpl.geoaxes as geoaxes
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
-from meso_fab_mc import BuildHarmonics
+from mcfab import BuildHarmonics
 import matplotlib.patheffects as path_effects
 
 plt.rcParams.update({
@@ -207,7 +208,7 @@ def loadEGRIP(loc):
 #Custom legend
 from matplotlib.lines import Line2D
 legend_elements = [Line2D([0], [0], color='k', lw=2, label='SpecCAF'),
-                     Line2D([0], [0], color='k', lw=2, linestyle='--', label=r"$\tilde{\lambda}' = \tilde{\lambda}/4$"),
+                     Line2D([0], [0], color='k', lw=2, linestyle='--', label=r"$\tilde{\lambda}' = 0.25\tilde{\lambda}$"),
                    Line2D([0], [0], marker='o', color='w', label='EGRIP',
                           markerfacecolor='k', markersize=4),
                           Line2D([0], [0], marker='x', color='w', label='Pole figures',
@@ -244,7 +245,7 @@ def loadEGRIP(loc):
 
 fig,axs = plt.subplots(1, 3, sharey=True,figsize=(6,4))
 
-colors = sns.color_palette("cmo.ice", 3)
+colors = sns.color_palette("deep", 3)
 
 
 
@@ -291,7 +292,6 @@ def loadEGRIP(loc):
 
 
 #%%
-L=8
 mmax = 8
 vmax=2
 
@@ -301,7 +301,7 @@ def J(odf):
     for l in range(0,odf.L+1,2):
         Sff = 0*Sff
         for m in range(0,l+1,1):
-            Sff=Sff+np.abs(odf.f[odf.sh.idx(l,abs(m))])**2
+            Sff=Sff+np.abs(odf.f[odf.idx(l,abs(m))])**2
         J=J+Sff
     return J
 
@@ -348,25 +348,24 @@ def angle_correction(xyz):
 espec = np.sort(np.array([ev_n[j],ev_z[j],ev_s[j]]))
 espec2 = np.sort(np.array([ev_n2[j],ev_z2[j],ev_s2[j]]))
 
-n1 = paths[j].n[-1,...]
-m1 = paths[j].m[-1,...]
+f1 = paths[j].f[-1,...]
+
+f2 = paths2[j].f[-1,...]
 
-n2 = paths2[j].n[-1,...]
-m2 = paths2[j].m[-1,...]
 
 fig,ax = plt.subplots(1,3,figsize=(6,3.3),subplot_kw=\
                       {'projection':ccrs.AzimuthalEquidistant(90,90)})
-odf1 = BuildHarmonics(n1,m1,L,mmax)
-odf2 = BuildHarmonics(n2,m2,L,mmax)
+odf1 = sff.Plotting(L,f1)
+odf2 = sff.Plotting(L,f2)
 odf_exp = BuildHarmonics(xyz,m,L,mmax)
 
 odf1.plot(fig,ax[0],hemisphere=True)
 odf2.plot(fig,ax[1],hemisphere=True)
 odf_exp.plot(fig,ax[2],hemisphere=True)
 
 
-J1 = J(odf1)
-J2 = J(odf2)
+J1 = odf1.J()
+J2 = odf2.J()
 J_exp = J(odf_exp)
 
 # ax[0].set_title('(a) SpecCAF'\
@@ -413,17 +412,17 @@ def angle_correction(xyz):
     espec = np.sort(np.array([ev_n[j],ev_z[j],ev_s[j]]))
     espec2 = np.sort(np.array([ev_n2[j],ev_z2[j],ev_s2[j]]))
 
-    n1 = paths[j].n[-1,...]
-    m1 = paths[j].m[-1,...]
 
-    n2 = paths2[j].n[-1,...]
-    m2 = paths2[j].m[-1,...]
+    f1 = paths[j].f[-1,...]
+
+    f2 = paths2[j].f[-1,...]
+
 
     # ax = subfig.subplots(nrows=1, ncols=3,subplot_kw=\
     #                         {'projection':ccrs.AzimuthalEquidistant(90,90)})    
     ax = axs[row,:]
-    odf1 = BuildHarmonics(n1,m1,L,mmax)
-    odf2 = BuildHarmonics(n2,m2,L,mmax)
+    odf1 = sff.Plotting(L,f1)
+    odf2 = sff.Plotting(L,f2)
     odf_exp = BuildHarmonics(xyz,m,L,mmax)
 
     pcol1 = odf1.plot(fig,ax[0],hemisphere=True,colorbar=True,pad=0.05)
@@ -432,8 +431,8 @@ def angle_correction(xyz):
 
 
 
-    J1 = J(odf1)
-    J2 = J(odf2)
+    J1 = odf1.J()
+    J2 = odf2.J()
     J_exp = J(odf_exp)
 
     fmax1 = pcol1.get_clim()[1]
@@ -447,7 +446,7 @@ def titlestr(J,fmax):
     colnumeral  = ['i','ii','iii']
 
     ax[0].set_title('(' +colnumeral[0] + ') SpecCAF' + titlestr(J1,fmax1))
-    ax[1].set_title('(' +colnumeral[1] + r") $\tilde{\lambda}' = \tilde{\lambda}/4$" + titlestr(J2,fmax2))
+    ax[1].set_title('(' +colnumeral[1] + r") $\tilde{\lambda}' = 0.25\tilde{\lambda}$" + titlestr(J2,fmax2))
     ax[2].set_title('(' +colnumeral[2] + ') EGRIP ice core data' + titlestr(J_exp,fmax_exp))
 
     #subfig.suptitle(rowletter[row] + f' Pole figures at {depth:.0f} m',y=0.96)

parameters.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+
+
+
+
+def Richards2021(T):
+    T[T<-30] = -30
+    iotaD = 0.0259733*T + 1.95268104
+    iotaS = np.zeros_like(T)
+    lambtilde = (0.00251776*T + 0.41244777)/np.sqrt(2) #correction to agree with SpecCAF
+    betatilde = 5*(0.35182521*T + 12.17066493)/np.sqrt(2)
+
+    Ecc = np.ones_like(T)
+    Eca = np.ones_like(T)
+    power = np.ones_like(T)
+
+
+    x = np.array([iotaD,iotaS,lambtilde,betatilde,Ecc,Eca,power])
+
+    return x.T
+
+
+def Richards2021Reduced(T,reduce=0.25):
+    T[T<-30] = -30
+    iotaD = 0.0259733*T + 1.95268104
+    iotaS = np.zeros_like(T)
+    lambtilde = reduce*(0.00251776*T + 0.41244777)/np.sqrt(2) #correction to agree with SpecCAF
+    betatilde = 5*(0.35182521*T + 12.17066493)/np.sqrt(2)
+
+    Ecc = np.ones_like(T)
+    Eca = np.ones_like(T)
+    power = np.ones_like(T)
+
+
+    x = np.array([iotaD,iotaS,lambtilde,betatilde,Ecc,Eca,power])
+
+    return x.T
+
+def Elmer(T):
+    iotaD = 0.94*np.ones_like(T)
+    iotaS = 0.6*np.ones_like(T)
+    lambtilde = np.sqrt(2)*2e-3 * np.exp(np.log(10)*T/10)
+    betatilde = np.zeros_like(T)
+    Ecc = np.ones_like(T)
+    Eca = 25*np.ones_like(T)
+    power = np.ones_like(T)
+
+    x = np.array([iotaD,iotaS,lambtilde,betatilde,Ecc,Eca,power])
+
+    return x.T

plotdivide.py

@@ -22,6 +22,7 @@
 colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
 
 npoints=10000
+L=20
 location = 'GRIP'#,'DomeF','DomeC','Talos']
 #locations = ['DomeF']
 
@@ -45,9 +46,9 @@
 
 
 
-a2,a4,n,m = track.fabric(p,npoints)
+a2,a4,f = track.fabric_sf(p,L)
 
-a2r,a4r,nr,mr = track.fabric(p,npoints,x='Reduced')
+a2r,a4r,fr = track.fabric_sf(p,L,x='Reduced')
 
 eigvals = np.linalg.eigvals(a2)
 
@@ -205,7 +206,8 @@ def inset(fig,ax,insx,insy,j,r=0.1,hemisphere=True,vmax=0.7):
 fig,ax = plt.subplots(1,3,figsize=(6,4),sharey=True)
 
 import cmocean
-colors = sns.color_palette("cmo.ice", 3)
+colors = sns.color_palette("deep", 3)
+
 
 #plot smallest eigenvalue in left plot, etc.
 ax[0].plot(eigvals[:,0],p.d,color=colors[0],linewidth=2,label='SpecCAF')

quickplot.py

@@ -52,7 +52,7 @@
 ev_z = np.zeros(len(paths))
 for i in tqdm(range(len(paths))):
     paths[i].Temperature()
-    paths[i].fabric_mc(npoints,)
+    paths[i].fabric_sf(L)
     depths[i] = paths[i].d[-1]
 
     ## Due to lack of vertical shear we know one eigenvector is (0,0,1),
@@ -240,6 +240,8 @@ def inset(fig,ax,insx,insy,j,tind=-1,r=0.1,lon=90,lat=90,hemisphere=True,vmax=0.
 
 
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+import specfabfuns as sff
+vmax = 0.7
 def inset(fig,ax,insx,insy,j,tind=-1,r=0.1,lon=90,lat=90,hemisphere=False,vmax=0.7):
     centre=ax.transAxes.inverted().transform(ax.transData.transform((insx, insy)))
     long_s = np.arctan2(paths[j].vy_s[-1],paths[j].vx_s[-1])*180/np.pi
@@ -252,7 +254,8 @@ def inset(fig,ax,insx,insy,j,tind=-1,r=0.1,lon=90,lat=90,hemisphere=False,vmax=0
                                  ))
 
     #odf = Reconstruct(paths[j].mmc.n,paths[j].mmc.m)
-    odf = BuildHarmonics(paths[j].n[tind,...],paths[j].m[tind,...],L,mmax)
+    # odf = BuildHarmonics(paths[j].n[tind,...],paths[j].m[tind,...],L,mmax)
+    odf = sff.Plotting(L,paths[j].f[tind,...])
     pcol = odf.plot(fig,ins,hemisphere,vmax=vmax,colorbar=False)
 
     geo = ccrs.RotatedPole()