Hi there, I'm trying to reproduce a model in the literature. Here is the code that causes the error (description+stack trace follows):
electron_chg = 1.6*10**(-19) # Coulombs
eps0 = 8.85418782e-12 # m^-3 kg^-1 s^4 A^2 = F/m
angs_per_meter = 1*10**10
T1 = 108
T11 = 120
T111 = 120
a10 = 3.34e5
a110 = 4.69e6
a1110 = -5.52e7
a12 = 3.23e8
a112 = 4.47e9
a123 = 4.91e9
m110 = -202e6
m1110 = 276e7
var('Px,Py,Pz, T')
a1(T) = a10*(T - T1)
a11(T) = a110*(T - T11)+m110
a111(T) = a1110*(T - T111)+m1110
Glan(Px,Py,Pz, T) = (a1(T))*(Px**2+Py**2+Pz**2) + (a11(T))*(Px**4+Py**4+Pz**4) + (a111(T))*(Px**6+Py**6+Pz**6) + a12*((Px**2)*(Py**2)+(Pz**2)*(Py**2)+(Px**2)*(Pz**2)) + a112*((Px**2)*(Py**4+Pz**4)+(Py**2)*(Px**4+Pz**4)+(Pz**2)*(Px**4+Py**4)) + a123*(Px**2)*(Py**2)*(Pz**2)
Ex(Px,Py,Pz, T) = derivative(Glan,Px)
Ey(Px,Py,Pz, T) = derivative(Glan,Py)
Ez(Px,Py,Pz, T) = derivative(Glan,Pz)
# this works
var('Ptet')
find_root(Ez(0,0,Ptet,32),0.0001,1.0)
# this doesn't, even for i = 0, which is the case that works above
var('Ptet')
Ptetx = zero_vector(101)
Ptety = zero_vector(101)
Ptetz = zero_vector(101)
for i in range(101):
print(i)
Ptetx[i] = find_root(Ex(0,0,Ptet,32+i),0.0001,1.0)
Ptety[i] = find_root(Ey(0,0,Ptet,32+i),0.0001,1.0)
Ptetz[i] = find_root(Ez(0,0,Ptet,32+i),0.0001,1.0)
Ptetx
Here is the full trace of the error:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-9-3d8c5561fbc3> in <module>
5 for i in range(Integer(101)):
6 print(i)
----> 7 Ptetx[i] = find_root(Ex(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
8 Ptety[i] = find_root(Ey(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
9 Ptetz[i] = find_root(Ez(Integer(0),Integer(0),Ptet,Integer(32)+i),RealNumber('0.0001'),RealNumber('1.0'))
/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/modules/free_module_element.pyx in sage.modules.free_module_element.FreeModuleElement.__setitem__ (build/cythonized/sage/modules/free_module_element.c:13856)()
1813 if n < 0 or n >= d:
1814 raise IndexError("vector index out of range")
-> 1815 self.set_unsafe(n, R(value))
1816
1817 cdef int set_unsafe(self, Py_ssize_t i, value) except -1:
/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9337)()
898 if mor is not None:
899 if no_extra_args:
--> 900 return mor._call_(x)
901 else:
902 return mor._call_with_args(x, args, kwds)
/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6043)()
285 raise TypeError("Cannot coerce {} to {}".format(x, C))
286 cdef Map m
--> 287 cdef Element e = method(C)
288 if e is None:
289 raise RuntimeError("BUG in coercion model: {} method of {} returned None".format(self.method_name, type(x)))
/opt/sagemath-9.2/local/lib/python3.7/site-packages/sage/rings/real_mpfr.pyx in sage.rings.real_mpfr.RealNumber._integer_ (build/cythonized/sage/rings/real_mpfr.c:17319)()
2227 return n
2228
-> 2229 raise TypeError("Attempt to coerce non-integral RealNumber to Integer")
2230
2231 def integer_part(self):
TypeError: Attempt to coerce non-integral RealNumber to Integer
The error describes itself as a type coercion, but I can't see what the origin is, given that the 0 case works outside the loop. Could it be related to how range() works, is it different than the python implementation?
