Fill a region of an image based on inequality equations

Question

I want to plot a filled-in region on a 2D plot using matplotlib in python. The implicit region is defined by,

def f(x,y):
    return sqrt(2*x*y)-x-y

I want to add to an existing image a filled-in region in the XY plane defined by the two inequalities,

sqrt(x**2+y**2)-2<=f(x,y) && f(x,y)>=1

To be clear, I want to plot a region defined by all (x,y) values that satisfy the above inequalities. I have tried getting it to work with pyplot.contour but I can't get it to work. What is the correct approach?

EDIT: Below is an example image of my desired output:

This is an image of a planet casting a shadow on its rings based on its positioning with its star. I created that image in Mathematica, and now I need to create something equivalent in python. I have a plot of the planet and the rings already, and I am trying to add the shadow.

The shadow is defined by two inequalities and resides in the XY plane. The equations I gave above are highly simplified for convenience; I do not need help with the math, only implementing it in a plot.

In the python version, I also want the shadow to be a filled-in region that resides in the XY plane.

Answer 1

Directly using a contour plot will have problems because the function is only defined in two of the four quadrants. Therefore it is advisable to first transform the inequalities into a binary array over a grid.

To plot this you may use an imshow plot.

import numpy as np
import matplotlib.pyplot as plt

f = lambda x,y : np.sqrt(2*x*y)-x-y
g = lambda x,y : np.sqrt(x**2+y**2)-2

d = np.linspace(-2,2,2000)
x,y = np.meshgrid(d,d)

im = plt.imshow( ((g(x,y)<=f(x,y)) & (f(x,y)<=1)).astype(int) , 
                extent=(x.min(),x.max(),y.min(),y.max()),origin="lower", cmap="Greys")

plt.show()

You may equally use a contourf plot; here, in order to get the contour close to the binary values, it may make sense to use two levels close to 1. It'll also require to define an all-black colormap.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors

f = lambda x,y : np.sqrt(2*x*y)-x-y
g = lambda x,y : np.sqrt(x**2+y**2)-2

d = np.linspace(-2,2,2000)
x,y = np.meshgrid(d,d)
z = ((g(x,y)<=f(x,y)) & (f(x,y)<=1)).astype(int) 

cmap = matplotlib.colors.ListedColormap(["k","k"])
plt.contourf(x,y,z, [0.99999, 1.00001], cmap=cmap)

plt.show()