ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Tue, 07 Mar 2023 11:49:36 +0100Switching coordinate charts mid-integration with solve_across_chartshttps://ask.sagemath.org/question/66826/switching-coordinate-charts-mid-integration-with-solve_across_charts/ I am trying to utilize `solve_across_charts` to avoid coordinate singularities by switching from spherical to cartesian coordinates when integrating geodesics in the Kerr spacetime. However the code just keeps infinitely running. Following the documentation, my code is (9.3 version),
M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
X.<t,r,th,ph>=M.chart(r"t r:(2.011,+infinity) th:(0,pi):\theta ph:(0,2*pi):periodic\phi")
C.<T,x,y,z>=M.chart()
C.add_restrictions(sqrt(x**2+y**2+z**2) > 2.011)
X_to_C=X.transition_map(C,[t,r*sin(th)*cos(ph),r*sin(th)*sin(ph), r*cos(th)])
C_to_X=C.transition_map(X,[T,sqrt(x**2+y**2+z**2),atan(sqrt(x**2+y**2)/z),atan(y/x)])
var('m,b,a', domain='real')
m=1
g = M.metric()
rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2
g[0,0] = -(1-2*m*r/rho2)
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2
def g00(m,r,th,a):return g[0,0](m,r,th,a)
def g03(m,r,th,a):return g[0,3](m,r,th,a)
def g11(m,r,th,a):return g[1,1](m,r,th,a)
def g22(m,r,th,a):return g[2,2](m,r,th,a)
def g33(m,r,th,a):return g[3,3](m,r,th,a)
def D(m,r,th,a): return (g03(m,r,th,a))^2-g00(m,r,th,a)*g33(m,r,th,a)
def initial_vector(r0, b,al,ph0=0, E=1, inward=False):
print(b)
t0,th0=0,pi/2
L = -b*E
vth0 = al/(r0*r0)
vt0=(E*g33(m,r0,th0,a)+L*g03(m,r0,th0,a))/(D(m,r0,th0,a))
vr0=sqrt((((1/D(m,r0,th0,a))*((L^2)*g00(m,r0,th0,a)+2*E*L*g03(m,r0,th0,a)+(E^2)*g33(m,r0,th0,a)))
-((vth0)^2)*g22(m,r0,th0,a))/g11(m,r0,th0,a))
if inward:
vr0 = - vr0
vph0 = -(1/D(m,r0,th0,a))*(E*g03(m,r0,th0,a)+L*(g00(m,r0,th0,a)))
p0 = M((t0, r0, th0, ph0), name='p_0')
return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')
r0 = 100
th0=pi/2
s = var('s')
ph0=0
v0 = initial_vector(r0, b=0,al=0, ph0=ph0, inward=True)
geod = M.integrated_geodesic(g, (s, 0, 150), v0, across_charts=True)
sol = geod.solve_across_charts(step=0.1, method='odeint', parameters_values={a:0.5},
verbose='False')
interp = geod.interpolate()
I think the issue is on the declaration of the transition maps, however I've tried everything I could think of with no luck. Any help would be much appreciated. BubbzTue, 07 Mar 2023 11:49:36 +0100https://ask.sagemath.org/question/66826/Evaluating the tangent vector of a previously saved and loaded integrated geodesichttps://ask.sagemath.org/question/62165/evaluating-the-tangent-vector-of-a-previously-saved-and-loaded-integrated-geodesic/I am trying to figure out how to load a previously computed and saved integrated geodesic and find the tangent vector using `geod.tangent_vector_eval_at` at a later time. The problem is that i can not seem to find a way to refer a `solution_key=` to the solution after i've loaded it in order to interpolate it. Is there a way to do that at a time different than the one that the solution is being computed? Thanks in advance.BubbzTue, 26 Apr 2022 11:33:08 +0200https://ask.sagemath.org/question/62165/