Dabade2016_drift.py

#!/usr/bin/env python
# coding: utf-8

# ![logo](JFM-notebooks-logo.jpg)

# ### $\mathit{Import \ useful \ modules}$

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
from pylab import *
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym


# ### $\mathit{Plot \ format}$

# In[2]:


SMALL_SIZE = 10
MEDIUM_SIZE = 20
BIGGER_SIZE = 24
plt.rc('font', size=BIGGER_SIZE)         # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('axes', labelpad=30)              # rotation
plt.rc('text', usetex=True) 
mpl.rcParams['figure.dpi'] = 300


# ### $\mathit{Load \ the \ reference \ for \ validation: \ Dabade \ et \ al., \ JFM, \ 2016}$

# In[3]:


###Particle inertia
#prolate: Figure 2.b, \xi_0 = 1.01
rod_pine = np.loadtxt('Rod_particle_ine_drift.txt')
#oblate: Figure 3, \xi_0 = 1.0001
disk_pine = np.loadtxt('Disk_particle_ine_drift.txt')
###Fluid inertia
#prolate: Figure 5.b, \xi_0 = 1.0001
rod_fine = np.loadtxt('Rod_fluid_ine_drift.txt')
#oblate: Figure 6.b, \xi_0 = 1.0001
disk_fine = np.loadtxt('Disk_fluid_ine_drift.txt')


# ### $\mathit{Define \ the \ symbols:}$

# In[4]:


###kappa
k = sym.symbols('\kappa')
###\xi_0
xi = sym.symbols('xi', positive=True)
###\overline{\xi_0}
xib= sym.symbols("\overline{xi}", positive=True)
###orbit constant C
C = sym.symbols('C')


# ### $\mathit{Define \ the \ inverse \ hyperbolic \ cotangent \ function}$

# In[5]:


###inverse hyperbolic cotangent function
def invcoth(val):
    return 0.5 *  sym.log((val+1)/(val-1))


# ### $\mathit{Define \ the \ aspect \ ratio \ depending \ functions \ F_i^p \ , G_i^p \ for \ particle \ inertia: \ equations \ 5.7 - 5.12}$

# In[6]:


F1p_def = (4*xi**4-7*xi**2+3)*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2)
F2p_def = xib**4*(invcoth(xi)+xi*(xi*invcoth(xi)-1))/(40*xi*(1-2*xi**2)**2)
F3p_def = - xib**2*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2)
F4p_def = F3p_def
F5p_def = -F3p_def / 2
F6p_def = F5p_def
###
G1p_def = -(3*xi**4-5*xi**2+2)*((xi**2+1)*invcoth(xi)-xi) / (40*xi*(1-2*xi**2)**2)
G2p_def = (xi*(-xi**3+xi+(xi**4-1)*invcoth(xi))) / (40*(1-2*xi**2)**2)
G3p_def = G2p_def / xi**2
G4p_def = - G2p_def / xi**2


# ### $\mathit{Define \ the \ aspect \ ratio \ depending \ functions \ F_i^f, \ G_i^f \ for \ fluid \ inertia: \ equations \ 6.1 - 6.8}$

# In[7]:


F1f_def = (xi**(2)*(-648*xi**(12)+1350*xi**(10)-5571*xi**(8)+11841*xi**(6)-9269*xi**(4)+2263*xi**(2)+6) \
        -27*xi**(2)*(24*xi**(8)-14*xi**(6)-19*xi**(4)+16*xi**(2)-3)*xib**(8)*invcoth(xi)**(4)+9*xi \
         *(288*xi**(12)-564*xi**(10)-20*xi**(8)+799*xi**(6)-743*xi**(4)+261*xi**(2)-29)*xib**(4) \
         * invcoth(xi)**(3)+xi*(2592*xi**(14)-7020*xi**(12)+13932*xi**(10)-21123*xi**(8)+14255*xi**(6) \
        -577*xi**(4)-2711*xi**(2)+652)*invcoth(xi)-3*(1296*xi**(16)-4320*xi**(14)+5346*xi**(12) \
        -1477*xi**(10)-4260*xi**(8)+6116*xi**(6)-3492*xi**(4)+849*xi**(2)-58)*invcoth(xi)**(2)) \
         *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)
F2f_def = -(xib**(2)*(-9*xi**(9)+30*xi**(7)-115*xi**(5)+90*xi**(3)-12*xi \
        +9*xib**(8)*(xi**(2)+1)*xi**(2)*invcoth(xi)**(3)-3*xib**(4)*(9*xi**(6)-10*xi**(4)-17*xi**(2)+14)*xi*invcoth(xi)**(2) \
        +(27*xi**(10)-87*xi**(8)+133*xi**(6)-33*xi**(4)-52*xi**(2)+12)*invcoth(xi))) \
        *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi)))**(-1)
F3f_def = -(xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \
        -27*xi**(2)*(6*xi**(6)+xi**(4)-4*xi**(2)+1)*xib**(8)*invcoth(xi)**(4)+9*xi*(12*xi**(10)+28*xi**(8)-201*xi**(6) \
        +273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3)+(-972*xi**(13)-324*xi**(11)+7365*xi**(9) \
        -10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi)*invcoth(xi)+3*(216*xi**(14)-378*xi**(12) \
        +109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4)+311*xi**(2)-22)*invcoth(xi)**(2)) \
        *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2)*(-3*xi**(3)+5*xi \
        +3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)
F4f_def = F3f_def
F5f_def = - F3f_def / 2.0
F6f_def = - F3f_def / 2.0
###
G1f_def = (xi**(2)*(81*xi**(10)-414*xi**(8)+1074*xi**(6)-1162*xi**(4)+479*xi**(2)-54) \
        +9*xi**(2)*(9*xi**(6)-7*xi**(2)+2)*xib**(8)*invcoth(xi)**(4)-3*xi*(108*xi**(10)-246*xi**(8)+69*xi**(6) \
        +167*xi**(4)-129*xi**(2)+23)*xib**(4)*invcoth(xi)**(3)+(-324*xi**(13)+1566*xi**(11)-3309*xi**(9) \
        +3133*xi**(7)-1023*xi**(5)-79*xi**(3)+36*xi)*invcoth(xi) \
        +(486*xi**(14)-2214*xi**(12)+3819*xi**(10)-2568*xi**(8)-222*xi**(6)+1036*xi**(4)-355*xi**(2)+18) \
        *invcoth(xi)**(2))*(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)
G2f_def = (-xi**(2)*(27*xi**(10)-180*xi**(8)+204*xi**(6)+68*xi**(4)-133*xi**(2)+18) \
        -9*xi**(4)*(3*xi**(4)+2*xi**(2)-1)*xib**(8)*invcoth(xi)**(4)+3*xi \
         *(36*xi**(10)-78*xi**(8)+73*xi**(6)-69*xi**(4)+35*xi**(2)-5)*xib**(4)*invcoth(xi)**(3) \
        +xi*(108*xi**(12)-630*xi**(10)+1041*xi**(8)-617*xi**(6)+115*xi**(4)-29*xi**(2)+12) \
         *invcoth(xi)+(-162*xi**(14)+810*xi**(12)-1551*xi**(10)+1600*xi**(8) \
        -1054*xi**(6)+448*xi**(4)-97*xi**(2)+6)*invcoth(xi)**(2)) \
         *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)
G3f_def = (xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \
        -27*xi**(2)*(6*xi**(6)+xi**(4)-4*xi**(2)+1)*xib**(8)*invcoth(xi)**(4) \
        +9*xi*(12*xi**(10)+28*xi**(8)-201*xi**(6)+273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3) \
        +(-972*xi**(13)-324*xi**(11)+7365*xi**(9)-10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi) \
        *invcoth(xi)+3*(216*xi**(14)-378*xi**(12)+109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4) \
        +311*xi**(2)-22)*invcoth(xi)**(2))*(120*xi**(2)*(2*xi**(2)-1)**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)
G4f_def = - G3f_def


# ### $\mathit{Define \ the \ particles}$

# In[8]:


###dictionary with the suggested values of \xi_0 to match the reference figures for validation
ee = {'pro_pine' : 1.01,
      'obl_pine' : 1.0001,
      'pro_fine' : 1.1,
      'obl_fine' : 1.0001}
###definition of particle aspect ratio according to prolate/oblate particle shape
rr = {'prolate' : lambda xi0 : xi0/(xi0**2-1)**0.5,
      'oblate'  : lambda xi0 : (xi0**2-1)**0.5/xi0}


# ### $\mathit{Function \ that \ evaluates \ the \ shape \ dependent \ coefficients \ for \ both \ particle \ and \ fluid \ inertia: the \ prolate - oblate \ transformation \ is \ automatically \ applied.}$

# In[9]:


def shape_coeffs(key,ee):
    if(key == 'prolate'):#prolates
        ###Particle inertia
        e0 = ee['pro_pine']
        F1p_eff = F1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F2p_eff = F2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F3p_eff = F3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F4p_eff = F4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F5p_eff = F5p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F6p_eff = F6p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G1p_eff = G1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G2p_eff = G2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G3p_eff = G3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G4p_eff = G4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        ###Fluid inertia
        e0 = ee['pro_fine']
        F1f_eff = F1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F2f_eff = F2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F3f_eff = F3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F4f_eff = F4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F5f_eff = F5f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        F6f_eff = F6f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G1f_eff = G1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G2f_eff = G2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G3f_eff = G3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
        G4f_eff = G4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()
    else:#oblates
        ###Particle Inertia
        e0 = ee['obl_pine']
        F1p_tmp = ((F1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F1p_eff = -sym.re(F1p_tmp.subs({xi:e0}).evalf())
        F2p_tmp = ((F2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F2p_eff = -sym.re(F2p_tmp.subs({xi:e0}).evalf())
        F3p_tmp = ((F3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F3p_eff = -sym.re(F3p_tmp.subs({xi:e0}).evalf())
        F4p_tmp = ((F4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F4p_eff = -sym.re(F4p_tmp.subs({xi:e0}).evalf())
        F5p_tmp = ((F5p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F5p_eff = -sym.re(F5p_tmp.subs({xi:e0}).evalf())
        F6p_tmp = ((F6p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F6p_eff = -sym.re(F6p_tmp.subs({xi:e0}).evalf())
        G1p_tmp = ((G1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G1p_eff = -sym.re(G1p_tmp.subs({xi:e0}).evalf())
        G2p_tmp = ((G2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G2p_eff = -sym.re(G2p_tmp.subs({xi:e0}).evalf())
        G3p_tmp = ((G3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G3p_eff = -sym.re(G3p_tmp.subs({xi:e0}).evalf())
        G4p_tmp = ((G4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G4p_eff = -sym.re(G4p_tmp.subs({xi:e0}).evalf())
        ###Fluid inertia
        e0 = ee['obl_fine']
        F1f_tmp = ((F1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F1f_eff = -sym.re(F1f_tmp.subs({xi:e0}).evalf())
        F2f_tmp = ((F2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F2f_eff = -sym.re(F2f_tmp.subs({xi:e0}).evalf())
        F3f_tmp = ((F3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F3f_eff = -sym.re(F3f_tmp.subs({xi:e0}).evalf())
        F4f_tmp = ((F4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F4f_eff = -sym.re(F4f_tmp.subs({xi:e0}).evalf())
        F5f_tmp = ((F5f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F5f_eff = -sym.re(F5f_tmp.subs({xi:e0}).evalf())
        F6f_tmp = ((F6f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        F6f_eff = -sym.re(F6f_tmp.subs({xi:e0}).evalf())
        G1f_tmp = ((G1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G1f_eff = -sym.re(G1f_tmp.subs({xi:e0}).evalf())
        G2f_tmp = ((G2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G2f_eff = -sym.re(G2f_tmp.subs({xi:e0}).evalf())
        G3f_tmp = ((G3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G3f_eff = -sym.re(G3f_tmp.subs({xi:e0}).evalf())
        G4f_tmp = ((G4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 
        G4f_eff = -sym.re(G4f_tmp.subs({xi:e0}).evalf())
    return (F1p_eff,F2p_eff,F3p_eff,F4p_eff,F5p_eff,F6p_eff,G1p_eff,G2p_eff,G3p_eff,G4p_eff, 
           F1f_eff,F2f_eff,F3f_eff,F4f_eff,F5f_eff,F6f_eff,G1f_eff,G2f_eff,G3f_eff,G4f_eff)


# ### $\mathit{Define \ the \ trigonometric \ integrals \ I_i \ and \ J_i \ introduced \ in \ Appendix C}$

# In[10]:


###I1
I1=2*sym.pi
###I2
I2=2*sym.pi*(k-1)*(k+1)**(-1)
###I3
I3=2*sym.pi*(2*((C**(2)+1)*(C**(2)*k**(2)+1))**(-1/2)-1)
###I4
I4=2*sym.pi*(k-1)**(2)*(k+1)**(-2)
###I5+I6
I5I6 = -(4*sym.pi*(2*k**(2)*(3*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-8*C**(2)-6) \
    +4*k*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)) \
    +sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))+4*(4*C**(2)+1)*k**(3)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)) \
    +k**(4)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-16*(C**(4)+C**(2))-2)-2)) \
    *((k**(2)-1)**(2)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)))**(-1)
###J1
J1=sym.pi*(k**(2)-1)*(k+1)**(-2)
###J2
J2 = -sym.pi*(-4*((C**2+1)*(C**2*k**2+1))**0.5+C**2*(k+1)**2+4)*C**(-2)*(k**2-1)**(-1) 
###J3
J3=sym.pi*(k-1)**(2)*(k+1)**(-2)
###J4
J4 = -sym.pi*(k**(2)-1)*(8*C**(4)*k**(3)+C**(2)*((k+1)**(4)-8*k**(2)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))) \
    -4*(k**(2)+1)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-1))*C**(-2)*(k**(2)-1)**(-3)


# ### $\mathit{Evaluate \ the \ shape \ dependent \ coefficients \ for \ a \ prolate \ particle}$

# In[11]:


F1p,F2p,F3p,F4p,F5p,F6p,G1p,G2p,G3p,G4p,F1f,F2f,F3f,F4f,F5f,F6f,G1f,G2f,G3f,G4f = shape_coeffs('prolate',ee)


# ### $\mathit{Evaluate \ the \ irreversible \ drift \ of \ a \ prolate \ particle \ due \ to \ particle \ inertia \ \Delta C_p ,\ as \ in \ eq. \ 5.19}$

# In[12]:


r = rr['prolate'](ee['pro_pine']) #select the correct aspect ratio to macth figure 2
PINE_pro = (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \
+(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee['pro_pine'],k:r}).evalf()


# ### $\mathit{Evaluate \ the \ irreversible \ drift \ of \ a \ prolate \ particle \ due \ to \ fluid \ inertia \ \Delta C_f, \ as \ in \ eq. 5.19}$

# In[13]:


r = rr['prolate'](ee['pro_fine']) #select the correct aspect ratio to match figure 5.b
FINE_pro = (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \
+(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee['pro_fine'],k:r}).evalf()


# ### $\mathit{Evaluate \ the \ shape \ dependent \ coefficients \ for \ an \ oblate \ particle}$

# In[14]:


F1p,F2p,F3p,F4p,F5p,F6p,G1p,G2p,G3p,G4p,F1f,F2f,F3f,F4f,F5f,F6f,G1f,G2f,G3f,G4f = shape_coeffs('oblate',ee)


# ### $\mathit{Evaluate \ the \ irreversible \ drift \ of \ an \ oblate \ particle \ due \ to \ particle \ inertia \ \Delta C_p ,\ as \ in \ eq. \ 5.19}$

# In[15]:


r = rr['oblate'](ee['obl_pine']) #select the correct aspect ratio to macth figure 3
PINE_obl = (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \
+(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee['obl_pine'],k:r}).evalf()


# ### $\mathit{Evaluate \ the \ irreversible \ drift \ of \ an \ oblate \ particle \ due \ to \ fluid \ inertia \ \Delta C_f, \ as \ in \ eq. 5.19}$

# In[16]:


r = rr['oblate'](ee['obl_fine']) #select the correct aspect ratio to match figure 6.b
FINE_obl= (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \
+(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee['obl_fine'],k:r}).evalf()


# ### $\mathit{Define \ a \ custom \ range \ of \ orbit \ constants \ C}$

# In[17]:


Crange = np.concatenate((np.linspace(0.00001,0.9,100),np.linspace(1,4,100)))
Crange = np.concatenate((Crange,np.linspace(4,100,100)))


# ### $\mathit{Calculate \ the \ drift \ over \ the \ C \ range}$

# In[18]:


CDC_pro_pine = [] #empty list to store the result for particle inertia on prolate
CDC_pro_fine = [] #empty list to store the result for fluid inertia on prolate
CDC_obl_pine = [] #empty list to store the result for particle inertia on oblate
CDC_obl_fine = [] #empty list to store the result for fluid inertia on oblate
for j in range(len(Crange)): #loop on the C range
    CVAL = Crange[j] #set the orbit constant value
    ###prolate, particle inertia
    data = PINE_pro.subs({C:CVAL}) #substitute the numerical values in the symbolic equation
    CDC_pro_pine.append([CVAL,data]) #store the result in the list
    ###prolate, fluid inertia
    data = FINE_pro.subs({C:CVAL}) #substitute the numerical values in the symbolic equation
    CDC_pro_fine.append([CVAL,data]) #store the result in the list
    ###oblate, particle inertia
    data = PINE_obl.subs({C:CVAL}) #substitute the numerical values in the symbolic equation
    CDC_obl_pine.append([CVAL,data]) #store the result in the list
    ###oblate, fluid inertia
    data = FINE_obl.subs({C:CVAL}) #substitute the numerical values in the symbolic equation
    CDC_obl_fine.append([CVAL,data]) #store the result in the list
#convert to numpy arrays
CDC_pro_pine = np.array(CDC_pro_pine)
CDC_pro_fine = np.array(CDC_pro_fine)
CDC_obl_pine = np.array(CDC_obl_pine)
CDC_obl_fine = np.array(CDC_obl_fine)


# ### $\mathit{Validate \ the \ method: \ fig. 2,b \ for \ the \ prolate \ particle \ inertia \ and \ fig. \ 3 \ for \ the \ oblate \ particle \ inertia}$

# In[19]:


###Create the figure with 4 sub-plots
fig, axs = plt.subplots(2,2,constrained_layout=True,sharex=True,figsize=(7,4))
###Figure 2.b of Dabade et al., 2016
ftr = -1/((ee['pro_pine']-1)**(3/2)*np.log(ee['pro_pine']-1))
axs[0,0].scatter(CDC_pro_pine[:,0]/(1+CDC_pro_pine[:,0]),CDC_pro_pine[:,1]/(1+CDC_pro_pine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.01}$')
axs[0,0].plot(rod_pine[:,0],rod_pine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--')
axs[0,0].set_ylim(-0.01,0.4)
axs[0,0].set_xlim(0.0,1.0)
axs[0,0].legend(fontsize=8,loc='lower center')
axs[0,0].title.set_text(r'$\mathit{Prolate, \ particle \ inertia}$')
axs[0,0].set_ylabel(r'$\mathit{ \frac{  {St}^{-1} \left( \Delta C_p / (C^2+1)  \right) }{ \left( (\xi_0^2-1)^{3/2} \log(\xi_0-1) \right)}  }$',fontsize=7,labelpad=5)
###Figure 5.b of Dabade et al., 2016
ftr = -((ee['pro_fine']-1)**(0.5)*np.log(ee['pro_fine']-1))
axs[0,1].scatter(CDC_pro_fine[:,0]/(1+CDC_pro_fine[:,0]),CDC_pro_fine[:,1]/(1+CDC_pro_fine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.1}$')
axs[0,1].plot(rod_fine[:,0],rod_fine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--')
axs[0,1].title.set_text(r'$\mathit{Prolate, \ fluid \ inertia}$')
axs[0,1].set_ylim(-0.01,0.2)
axs[0,1].set_ylabel(r'$\mathit{\frac{{Re_p}^{-1} \left( \Delta C_p / (C^2+1)  \right) }{ \left( (\xi_0-1)^{0.5} \log (\xi_0-1) \right)^{-1}} }$',fontsize=7,labelpad=5)
axs[0,1].legend(fontsize=8,loc='upper center')
###Figure 3 of Dabade et al., 2016
ftr = ee['obl_pine']**2
axs[1,0].scatter(CDC_obl_pine[:,0]/(1+CDC_obl_pine[:,0]),CDC_obl_pine[:,1]/(1+CDC_obl_pine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.0001}$')
axs[1,0].plot(disk_pine[:,0],disk_pine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--')
axs[1,0].title.set_text(r'$\mathit{Oblate, \ particle \ inertia}$')
axs[1,0].legend(fontsize=8,loc='upper center')
axs[1,0].set_ylim(-0.11,0.01)
axs[1,0].set_ylabel(r'$\mathit{{St}^{-1} \left( \Delta C_p / (C^2+1)  \right) \xi_0^2 }$',fontsize=7,labelpad=5)
###Figure 6.b of Dabade et al., 2016
ftr = (ee['obl_fine']-1)**0.5
axs[1,1].scatter(CDC_obl_fine[:,0]/(1+CDC_obl_fine[:,0]),CDC_obl_fine[:,1]/(1+CDC_obl_fine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.0001}$')
axs[1,1].plot(disk_fine[:,0],disk_fine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--')
axs[1,1].title.set_text(r'$\mathit{Oblate, \ fluid \ inertia}$')
axs[1,1].legend(fontsize=8,loc='upper center')
axs[1,1].set_ylim(-0.16,0.02)
axs[1,1].set_ylabel(r'$\mathit{{Re_p}^{-1} \left( \Delta C_f / (C^2+1)  \right) (\xi_0-1)^{0.5} }$',fontsize=6.5,labelpad=5)
###Format
for j in range(2):
    axs[1,j].set_xlabel(r'$\mathit{C/(C+1)}$',labelpad=5,fontsize=13)


# In[ ]: