Harrison - 2 months ago 15

Python Question

I'm trying to create a line or scatter plot of this algorithm and it gives me the error

Traceback (most recent call last):

File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 73, in module

plt.plot(xr, P(xr))

File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 55, in P

if x > r:

ValueError: "The truth value of an array with more than one element is ambiguous. Useor`a.any()`

."`a.all()`

I have looked up possible solutions to this error but I don't think any apply to me.

`import numpy as np`

import scipy.integrate as integ

import matplotlib.pyplot as plt

rho = 1860

rhos = 250 #Assuming Nitrogen Ice

rhom = 1000 #Assuming Water

rhoc = 3500 #Assuming a mix of Olivine and Pyroxene

def rhof(x):

if x > r:

return "Point is outside of the planet"

elif x < r and x > rm:

return rhos

elif x < rm and x > rc:

return rhom

else:

return rhoc

r = 1.187e6

rc = 8.5e5 #Hypothesized

rm = 9.37e5 #Estimated based on crustal thickness of 250 km

Ts = 44

B = 0.15

G = 6.67e-11

m = 1.303e22

mc = (4*np.pi*rhoc*rc**3)/3

mm = (4*np.pi*rhom*((rm**3) - (rc**3)))/3

ms = (4*np.pi*rhos*((r**3) - (rm**3)))/3

Ic = 0.4*mc*rc**2

Im = 0.4*mm*rm**2

Is = 0.4*ms*r**2

Itot = Is + Im + Ic

def gi(x):

if x == r:

return G*m/r**2

elif x > r:

return "Point is outside of the planet"

elif x > rc and x < rm:

return (G*mc/rc**2) + (G*mm/((x**2) - (rc**2)))

elif x > rm and x < r:

return (G*mc/rc**2) + (G*mm/((rm**2) - (rc**2))) + (G*ms/((x**2) - (rm**2)))

else:

return G*((3*rhoc)/4*np.pi*x**3)/x**2

def Psmb(z):

return rhos*G*(4.0/3.0)*np.pi*(1/z**2)*(rhom*(rm**3) + rhos*(z - rm**3))

def Pmcb(z):

return rhom*G*(4.0/3.0)*np.pi*(1/z**2)*(rhoc*(rc**3) + rhom*(z - rc**3))

def P(x):

if x > r:

return "The point is outside of the planet"

elif x == r:

return 1

elif x > rm and x < r:

return (integ.quad(1000*gi(x), x, r))[0]

elif x == rm:

return (integ.quad(Psmb, x, r))[0]

elif x > rc and x < rm:

return (integ.quad(1000*gi(x), x, rm) + P(rm))[0]

elif x == rc:

return (integ.quad(Pmcb, x, rm) + P(rm))[0]

elif x < rc and x != 0:

return (integ.quad(1000*gi(x), x, rc) + P(rc))[0]

else:

return ((2.0/3.0)*np.pi*G*(rhoc**2)*r**2)

xr = np.linspace(0, 1187000, 1000)

plt.plot(xr, P(xr))

print("Mass = " + str(m) + " kg")

print("Radius = " + str(r) + " m")

print("Density = " + str(rho) + " kg/m^3")

print("Moment of Inertia = " + str(Itot) + " kgm^2")

print("Mean Surface Temperature = " + str(Ts) + " K")

print("Magnetic Field = " + str(B) + " nT")

print("Surface Gravity = " + str(gi(r)) + " m/sec^2")

print("Pressure at Surface = " + str(P(r)) + " Pa")

print("Pressure at Crust-Mantle Boundary = " + str(P(rm)) + " Pa")

print("Pressure at Mantle-Core Boundary = " + str(P(rc)) + " Pa")

print("Pressure at the Center = " + str(P(0)) + " Pa")

Is there a way to plot this function without separating each condition?

Answer

Numpy has a vectorize function that might be what you are looking for.

```
pFunc = np.vectorize(P)
plt.plot(xr, pFunc(xr))
```

Vectorize is essentially a *for* loop, so wont speed up your code.

If you are able to use Pandas then I reckon you could also try using the apply function:

```
import pandas as pd
xd = pd.Series(xr)
yr = xd.apply(lambda x: P(x))
plt.plot(xr, yr)
```

I believe the reason you were getting that specific error was because your condition was in effect asking if the array xr was greater than a specific number. Some items in the array were, some were not, so the result was ambiguous. The error message was asking if you wanted the question to be "is *anything* in xr greater than this number" OR " is *everything* in xr greater than this number".

Note: you should return np.nan rather than a string in your first condition. Also, the function integ.quad needs something that is callable as the first argument.