Python - Returning number of elements in array above a curved function

Question

I'm trying to extract the number of values above the dashed, dotted-dashed, and dotted curves in my plot.

I can always eyeball the number of points that exceed the extent of these curves, but I do not want to risk any errors in this method.

First off, these curves are defined by a mass variable and a distance variable I call in my code Escape_Velocity. It is part of a class profile I have made.

import numpy as np
import NFW_Profile as NFW

h0 = 0.704
NFWc12 = NFW.Profile(c=12.5, z=0.0, h=h0)
nfwEsc = NFWc12.Escape_Velocity(Mass, Radius)

On my plot, the curves are based on this function, where I pick an arbitrary Mass value and Radius is fed a np.linspace value.

The values on my plot are from my analysis from data sets, where we can call RadVel to be the y-axis values "Radial Velocity" and Dist to be the x-axis values "Galactocentric Distance".

The key part is that both of these arrays of elements are of the same dimension and that they can be indexed with one another. For example the 20th indexed value RadVel corresponds to the 20th indexed value of Dist.

---

So my plan from is to

1. Pick an arbitrary value of Mass
2. Take an element from Dist
3. Feed that value to Escape_Velocity
4. If the RadVel value, corresponding to the element fed to Escape_Velocity, is greater than the Escape_Velocity value, it counts.

5. If the Escape_Velocity is greater than the corresponding RadVel value, it does not count.

   import numpy as np
   import NFW_Profile as NFW
   
   h0 = 0.704
   NFWc12 = NFW.Profile(c=12.5, z=0.0, h=h0)
   
   
   def UnboundSubhalos(mass, Dist):
       Vesc = NFWc12.Escape_Velocity(mass, Dist)
   
       PosValues = [i for i in RadVel if i => Vesc.all()]
       NegValues = [i for i in RadVel if i <= -Vesc.all()]
   
       return  len(PosValues) + len(NegValues)
   

Answer 1

As far as I got this, you want to find the number of red dots that are outside the region defined by a non-linear function you call Escape_Velocity and the red dots have x-positions and y-positions that are stored in the arrays Dist and RedVal, respectively, ordered so that a red dot has a position (x,y)=(Dist[n], RedVel[n]). Also, it looks like (and your definition of UnboundSubhalos suggests it) the dotted lines defined by Escape_Velocity are symmetric around RedVel = 0, i.e. you actually want to find the number of red dots with an absolute value of RedVel larger than Escape_Velocity.

In this case you can simply do the following: Assuming Escape_Velocity is a function you can pass an array to, you find the Escape_Velocity for every element of your array Dist and store it in the array esc_vel = Escape_Velocity(Mass, Dist). Else you will have to calculate these values with a loop over the elements of Dist. From your code I can assume that Escape_Velocity gives a positive value, i.e. the dotted line at positive RedVal values. Then the array outside_dots = np.abs(RedVel) >= esc_vel contains True for every dot outside the region enclosed by the dotted curves and False for every dot inside. Then np.sum(outside_dots) gives you the desired number of dots outside.

The following code works like described:

import numpy as np
import matplotlib.pyplot as plt

number_dots = 100
x = np.random.random(number_dots)  # random numbers in the range [0, 1]
y = (np.random.random(number_dots) - 0.5) * 2  # random numbers in the range [-1, 1]

curved_function = lambda z: 1 - np.sqrt(z)/3

plt.plot(x, y, 'ro')
plt.xlabel('x')
plt.ylabel('y')

def count_dots_outside(x, y):
        outside_dots = np.abs(y) > curved_function(x)
        return np.sum(outside_dots)

plt.title('Number of dots outside the curved function: '+str(count_dots_outside(x,y)))
xx = np.linspace(0,1,100)
plt.plot(xx, curved_function(xx), 'black', linestyle='dashed')
plt.plot(xx, -curved_function(xx), 'black', linestyle='dashed')
plt.show()