0

I'm working on this code that requires the use of elliptic integrals to solve an equation as a function of time. Most of the code is pretty straight forward but it runs into an error on the final equation containing the elliptic integrals, reading "cannot create mpf from array". Not quite sure if there is an easy fix but any insight would be greatly appreciated. The code can be found below:

import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
from mpmath import *

G = 6.6743*10**(-11) #Gravitational Constant
M = 10*(1.988*10**30) #Mass of the Black Hole (kg)
m = 1*(1.988*10**30) #Mass of the Companion Star (kg)
Mt = M + m #Total Mass(kg)
q = m/M #Mass Ratio
c = 2.99792458*10**8 #Speed of Light (m/s)
Period = 10 #Orbital Period (Days)
P = Period*86400 #Orbital Period (Seconds)
phi = 0.01*(np.pi/(180)) #Inclination from Line-of-Sight
t0 = Period/2 #Pulse Mid-Point
a = ((P**2*G*(Mt))/(4*np.pi**2))**(1/3) #Semi-Major Axis
omega = ((G*(m + M))/a**3)**0.5 #Angular Velocity
vtran = a*omega #Transverse Velocity
Rs = (2*G*M)/c**2 #Schwarzschild Radius
Rein = (2*Rs*a)**0.5 #Einstein Radius
te = (Rein/vtran)/86400 #Einstein Time
u0 = (a/Rein)*phi #Closest Angular Approach
pstar = ((1*(6.957*10**8))/Rein)
t = np.linspace(0,P,2500000) #Time Vector
u = ((u0**2)+((((t/86400)-t0)/te))**2)**0.5 #Angular Separation

n = ((4*u*pstar)/(u + pstar)**2)
k = ((4*n)/(4 + (u - pstar)**2))**0.5

Amp = (1/(2*(np.pi)))*(((((u+pstar)/(pstar**2))*np.sqrt(4+(u-pstar)**2)*ellipe(k**2))-(((u-pstar)/(pstar**2))*(((8+u**2-pstar**2)/(np.sqrt(4+(u-pstar)**2)))*ellipk(k**2)))+(((4*(u-pstar)**2)/(pstar**2*(u+pstar)))*(((1+pstar**2)/(np.sqrt(4+(u-pstar)**2)))*ellippi(n,k**2)))))
D.L
  • 4,339
  • 5
  • 22
  • 45
sorabella91
  • 1
  • 1
  • `mpmath` does not work with `numpy`. A `mpf` is a scalar, a single number. as a coding style isue, that `amp` ,ine that requires scroling is ugle and hard to read. I can't quickly tell where you are using `mpmath` – hpaulj Sep 13 '22 at 05:56

1 Answers1

0

The problem is specifically with mpmath.ellipe() and mpmath.ellipk() and mpath.ellippi().

The docs are here: https://mpmath.org/doc/current/functions/elliptic.html#ellipe (and similarly for ellipk and ellippi).

But in summary, as the comment has pointed out, mpmath.ellipe(*args) Called with a single argument m, evaluates the Legendre complete elliptic integral of the second kind, E(m)

You have passed in a list.

D.L
  • 4,339
  • 5
  • 22
  • 45
  • Thank you both for your answers! I have now circumvented the issue by doing a For loop for t instead of using linspace. This seems to have fixed everything. – sorabella91 Sep 13 '22 at 17:26