Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
601 views
unlisted
ubuntu2404
Tweet
Share
Share
1
#!/usr/bin/env python3

2
# -*- coding: utf-8 -*-

3
"""

4
Created on Fri Nov  7 18:39:55 2025

5


6
@author: davide

7
Methods that calculates the orbit variation delta_C for an oblate particle given its aspect ratio

8
"""

9


10


11


12


13
import sympy as sym

14
import numpy as np

15


16


17
###kappa

18
k = sym.symbols('\kappa')

19
###\xi_0

20
xi = sym.symbols('xi', positive=True)

21
###\overline{\xi_0}

22
xib= sym.symbols("\overline{xi}", positive=True)

23
###orbit constant C

24
C = sym.symbols('C')

25


26
###inverse hyperbolic cotangent function

27
def invcoth(val):

28
    return 0.5 *  sym.log((val+1)/(val-1))

29
###aspect ratio depending functions for particle inertia

30
F1p_def = (4*xi**4-7*xi**2+3)*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2)

31
F2p_def = xib**4*(invcoth(xi)+xi*(xi*invcoth(xi)-1))/(40*xi*(1-2*xi**2)**2)

32
F3p_def = - xib**2*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2)

33
F4p_def = F3p_def

34
F5p_def = -F3p_def / 2

35
F6p_def = F5p_def

36
###

37
G1p_def = -(3*xi**4-5*xi**2+2)*((xi**2+1)*invcoth(xi)-xi) / (40*xi*(1-2*xi**2)**2)

38
G2p_def = (xi*(-xi**3+xi+(xi**4-1)*invcoth(xi))) / (40*(1-2*xi**2)**2)

39
G3p_def = G2p_def / xi**2

40
G4p_def = - G2p_def / xi**2

41
###aspect ratio depending functions for fluid inertia

42
F1f_def = (xi**(2)*(-648*xi**(12)+1350*xi**(10)-5571*xi**(8)+11841*xi**(6)-9269*xi**(4)+2263*xi**(2)+6) \

43
        -27*xi**(2)*(24*xi**(8)-14*xi**(6)-19*xi**(4)+16*xi**(2)-3)*xib**(8)*invcoth(xi)**(4)+9*xi \

44
         *(288*xi**(12)-564*xi**(10)-20*xi**(8)+799*xi**(6)-743*xi**(4)+261*xi**(2)-29)*xib**(4) \

45
         * invcoth(xi)**(3)+xi*(2592*xi**(14)-7020*xi**(12)+13932*xi**(10)-21123*xi**(8)+14255*xi**(6) \

46
        -577*xi**(4)-2711*xi**(2)+652)*invcoth(xi)-3*(1296*xi**(16)-4320*xi**(14)+5346*xi**(12) \

47
        -1477*xi**(10)-4260*xi**(8)+6116*xi**(6)-3492*xi**(4)+849*xi**(2)-58)*invcoth(xi)**(2)) \

48
         *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \

49
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)

50
F2f_def = -(xib**(2)*(-9*xi**(9)+30*xi**(7)-115*xi**(5)+90*xi**(3)-12*xi \

51
        +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) \

52
        +(27*xi**(10)-87*xi**(8)+133*xi**(6)-33*xi**(4)-52*xi**(2)+12)*invcoth(xi))) \

53
        *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \

54
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi)))**(-1)

55
F3f_def = -(xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \

56
        -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) \

57
        +273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3)+(-972*xi**(13)-324*xi**(11)+7365*xi**(9) \

58
        -10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi)*invcoth(xi)+3*(216*xi**(14)-378*xi**(12) \

59
        +109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4)+311*xi**(2)-22)*invcoth(xi)**(2)) \

60
        *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2)*(-3*xi**(3)+5*xi \

61
        +3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)

62
F4f_def = F3f_def

63
F5f_def = - F3f_def / 2.0

64
F6f_def = - F3f_def / 2.0

65
###

66
G1f_def = (xi**(2)*(81*xi**(10)-414*xi**(8)+1074*xi**(6)-1162*xi**(4)+479*xi**(2)-54) \

67
        +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) \

68
        +167*xi**(4)-129*xi**(2)+23)*xib**(4)*invcoth(xi)**(3)+(-324*xi**(13)+1566*xi**(11)-3309*xi**(9) \

69
        +3133*xi**(7)-1023*xi**(5)-79*xi**(3)+36*xi)*invcoth(xi) \

70
        +(486*xi**(14)-2214*xi**(12)+3819*xi**(10)-2568*xi**(8)-222*xi**(6)+1036*xi**(4)-355*xi**(2)+18) \

71
        *invcoth(xi)**(2))*(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \

72
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)

73
G2f_def = (-xi**(2)*(27*xi**(10)-180*xi**(8)+204*xi**(6)+68*xi**(4)-133*xi**(2)+18) \

74
        -9*xi**(4)*(3*xi**(4)+2*xi**(2)-1)*xib**(8)*invcoth(xi)**(4)+3*xi \

75
         *(36*xi**(10)-78*xi**(8)+73*xi**(6)-69*xi**(4)+35*xi**(2)-5)*xib**(4)*invcoth(xi)**(3) \

76
        +xi*(108*xi**(12)-630*xi**(10)+1041*xi**(8)-617*xi**(6)+115*xi**(4)-29*xi**(2)+12) \

77
         *invcoth(xi)+(-162*xi**(14)+810*xi**(12)-1551*xi**(10)+1600*xi**(8) \

78
        -1054*xi**(6)+448*xi**(4)-97*xi**(2)+6)*invcoth(xi)**(2)) \

79
         *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \

80
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)

81
G3f_def = (xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \

82
        -27*xi**(2)*(6*xi**(6)+xi**(4)-4*xi**(2)+1)*xib**(8)*invcoth(xi)**(4) \

83
        +9*xi*(12*xi**(10)+28*xi**(8)-201*xi**(6)+273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3) \

84
        +(-972*xi**(13)-324*xi**(11)+7365*xi**(9)-10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi) \

85
        *invcoth(xi)+3*(216*xi**(14)-378*xi**(12)+109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4) \

86
        +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) \

87
        *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1)

88
G4f_def = - G3f_def

89
###define the trigonometric integrals

90
###I1

91
I1=2*sym.pi

92
###I2

93
I2=2*sym.pi*(k-1)*(k+1)**(-1)

94
###I3

95
I3=2*sym.pi*(2*((C**(2)+1)*(C**(2)*k**(2)+1))**(-1/2)-1)

96
###I4

97
I4=2*sym.pi*(k-1)**(2)*(k+1)**(-2)

98
###I5+I6

99
I5I6 = -(4*sym.pi*(2*k**(2)*(3*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-8*C**(2)-6) \

100
    +4*k*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)) \

101
    +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)) \

102
    +k**(4)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-16*(C**(4)+C**(2))-2)-2)) \

103
    *((k**(2)-1)**(2)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)))**(-1)

104
###J1

105
J1=sym.pi*(k**(2)-1)*(k+1)**(-2)

106
###J2

107
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) 

108
###J3

109
J3=sym.pi*(k-1)**(2)*(k+1)**(-2)

110
###J4

111
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))) \

112
    -4*(k**(2)+1)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-1))*C**(-2)*(k**(2)-1)**(-3)

113
    

114
def shape_coeffs(r,e0):

115
    """returns the shape-dependent coefficients"""

116
    if(r>1.0):#prolates

117
        ###Particle inertia

118
        F1p_eff = F1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

119
        F2p_eff = F2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

120
        F3p_eff = F3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

121
        F4p_eff = F4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

122
        F5p_eff = F5p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

123
        F6p_eff = F6p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

124
        G1p_eff = G1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

125
        G2p_eff = G2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

126
        G3p_eff = G3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

127
        G4p_eff = G4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

128
        ###Fluid inertia

129
        F1f_eff = F1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

130
        F2f_eff = F2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

131
        F3f_eff = F3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

132
        F4f_eff = F4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

133
        F5f_eff = F5f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

134
        F6f_eff = F6f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

135
        G1f_eff = G1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

136
        G2f_eff = G2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

137
        G3f_eff = G3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

138
        G4f_eff = G4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf()

139
    else:#oblates

140
        ###Particle Inertia

141
        F1p_tmp = ((F1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

142
        F1p_eff = -sym.re(F1p_tmp.subs({xi:e0}).evalf())

143
        F2p_tmp = ((F2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

144
        F2p_eff = -sym.re(F2p_tmp.subs({xi:e0}).evalf())

145
        F3p_tmp = ((F3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

146
        F3p_eff = -sym.re(F3p_tmp.subs({xi:e0}).evalf())

147
        F4p_tmp = ((F4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

148
        F4p_eff = -sym.re(F4p_tmp.subs({xi:e0}).evalf())

149
        F5p_tmp = ((F5p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

150
        F5p_eff = -sym.re(F5p_tmp.subs({xi:e0}).evalf())

151
        F6p_tmp = ((F6p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

152
        F6p_eff = -sym.re(F6p_tmp.subs({xi:e0}).evalf())

153
        G1p_tmp = ((G1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

154
        G1p_eff = -sym.re(G1p_tmp.subs({xi:e0}).evalf())

155
        G2p_tmp = ((G2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

156
        G2p_eff = -sym.re(G2p_tmp.subs({xi:e0}).evalf())

157
        G3p_tmp = ((G3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

158
        G3p_eff = -sym.re(G3p_tmp.subs({xi:e0}).evalf())

159
        G4p_tmp = ((G4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

160
        G4p_eff = -sym.re(G4p_tmp.subs({xi:e0}).evalf())

161
        ###Fluid inertia

162
        F1f_tmp = ((F1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

163
        F1f_eff = -sym.re(F1f_tmp.subs({xi:e0}).evalf())

164
        F2f_tmp = ((F2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

165
        F2f_eff = -sym.re(F2f_tmp.subs({xi:e0}).evalf())

166
        F3f_tmp = ((F3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

167
        F3f_eff = -sym.re(F3f_tmp.subs({xi:e0}).evalf())

168
        F4f_tmp = ((F4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

169
        F4f_eff = -sym.re(F4f_tmp.subs({xi:e0}).evalf())

170
        F5f_tmp = ((F5f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

171
        F5f_eff = -sym.re(F5f_tmp.subs({xi:e0}).evalf())

172
        F6f_tmp = ((F6f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

173
        F6f_eff = -sym.re(F6f_tmp.subs({xi:e0}).evalf())

174
        G1f_tmp = ((G1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

175
        G1f_eff = -sym.re(G1f_tmp.subs({xi:e0}).evalf())

176
        G2f_tmp = ((G2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

177
        G2f_eff = -sym.re(G2f_tmp.subs({xi:e0}).evalf())

178
        G3f_tmp = ((G3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

179
        G3f_eff = -sym.re(G3f_tmp.subs({xi:e0}).evalf())

180
        G4f_tmp = ((G4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) 

181
        G4f_eff = -sym.re(G4f_tmp.subs({xi:e0}).evalf())

182
    return (F1p_eff,F2p_eff,F3p_eff,F4p_eff,F5p_eff,F6p_eff,G1p_eff,G2p_eff,G3p_eff,G4p_eff, 

183
           F1f_eff,F2f_eff,F3f_eff,F4f_eff,F5f_eff,F6f_eff,G1f_eff,G2f_eff,G3f_eff,G4f_eff)

184


185


186
def inertial_contributions(r,ee):

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

188
    particle_inertia_prolate =  (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \

189
                +(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee,k:r}).evalf()

190
    fluid_inertia_prolate = (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \

191
                +(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee,k:r}).evalf()

192
    particle_inertia_oblate = (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \

193
                +(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee,k:r}).evalf()

194
    fluid_inertia_oblate = (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \

195
                +(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee,k:r}).evalf()

196
    if r>1.0:

197
        return particle_inertia_prolate,fluid_inertia_prolate

198
    elif r>0.0 and r<1.0:

199
        return particle_inertia_oblate,fluid_inertia_oblate

200


201
def calculate_eccentricity(r):

202
    """calculate particle eccentricity from aspect ratio"""

203
    if r > 1.0:

204
        return r/(r**2-1)**0.5

205
    elif r < 1.0 and r > 0.0:

206
        return 1/(1-r**2)**0.5

207
    else:

208
        raise ValueError('Please pass a valid aspect ratio > 0.0 and != 1.0')

209
def deltaC_dabade16(r=0.56,Crange=None):

210
    """returns particle and fluid inertia orbit variations for the given aspect ratio and range of orbit costants"""

211
    #calculate the eccentricity

212
    ee = calculate_eccentricity(r)

213
    #get the inertial contributions: particle and fluid

214
    PINE,FINE = inertial_contributions(r,ee)

215
    if Crange is None:

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

217
        Crange = np.concatenate((Crange,np.linspace(4,100,100)))

218
    CDC_pine = []

219
    CDC_fine = []

220
    for j in range(len(Crange)): #loop on the C range

221
        CVAL = Crange[j] #set the orbit constant value

222
        ###particle inertia

223
        data = PINE.subs({C:CVAL}) #substitute the numerical values in the symbolic equation

224
        CDC_pine.append([CVAL,data]) #store the result in the list

225
        ###fluid inertia

226
        data = FINE.subs({C:CVAL}) #substitute the numerical values in the symbolic equation

227
        CDC_fine.append(data) #store the result in the list

228
    #convert to numpy arrays

229
    CDC_pine = np.array(CDC_pine)

230
    CDC_fine = np.array(CDC_fine)

231
    return np.column_stack((CDC_pine,CDC_fine))

232