Find points that lie on a specific Gaussian contour

Question

I have a 2-D distribution of points (dark green points in the image) that I fit with a Bivariate Gaussian distribution and I have taken only the points with $f(x,y) > 0.1$ (yellow points in the image).

In the image are also present: the contour line (blue line) generated with the Python built-in function; the contour points (red points) calculated with my function. I did something wrong during the calculation because the points should line over the displayed Gaussian contour.

link to image since I can't post

My reasoning so far has been the following: I started from the equation of the Bivariate Gaussian Function with $\mu_x=0,\mu_y=0,\rho=0$: $$ f(x,y)=\frac{1}{2\pi\sigma_x\sigma_y} \exp \Bigg( -\frac{1}{2}(\frac{x^2}{\sigma^2_x}+\frac{y^2}{\sigma^2_y}) \Bigg) $$ $$ f(x,y)=c $$ In my specific case $c=0.1$. $$ \frac{1}{2\pi\sigma_x\sigma_y} \exp \Bigg( -\frac{1}{2}(\frac{x^2}{\sigma^2_x}+\frac{y^2}{\sigma^2_y}) \Bigg)=c $$ $$ \frac{x^2}{\sigma^2_x}+\frac{y^2}{\sigma^2_y} =-2\ln(2\pi\sigma_x\sigma_yc) $$ $$ \frac{x^2}{-2\ln(2\pi\sigma_x\sigma_yc)\sigma^2_x}+\frac{y^2}{-2\ln(2\pi\sigma_x\sigma_yc)\sigma^2_y} =1 $$ Thus the major and minor semi-axes are: $$ a=\sigma_x\sqrt{-2\ln(2\pi\sigma_x\sigma_yc)} $$ $$ b=\sigma_y\sqrt{-2\ln(2\pi\sigma_x\sigma_yc)} $$ And finally the $x$ and $y$ coordinates of the Ellipse: $$ 0\leq\theta<2\pi $$ $$ x=a\cos(\theta) $$ $$ y=b\sin(\theta) $$ Then I shifted and rotated the points using the $\mu$ vector and the rotation matrix $R$.

I used the following Python code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

def eigsorted(cov):
    vals, vecs = np.linalg.eigh(cov)
    order = vals.argsort()[::-1]
    return vals[order], vecs[:,order]

def createConfidenceEllipsePoints(M, p):
    T = M[:, [1,2]]
    mu = np.mean(T, axis=0)
    covariance = np.cov(T,rowvar=False)
    vals,vecs = eigsorted(covariance)
    theta = np.arctan2(*vecs[:,0][::-1])
    c = -2*np.log(2*np.pi*np.sqrt(covariance[0,0])*np.sqrt(covariance[1,1])*p)
    width, height = np.sqrt(c) *np.sqrt(np.array([covariance[0,0],covariance[1,1]]))
    rotation = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]])
    thetaGrid = np.linspace(0, 2*np.pi,100)
    P = np.column_stack((width*np.cos(thetaGrid),height*np.sin(thetaGrid)))
    P = np.dot(P,rotation)
    P = mu + P
    return P 

def gaussianContour(mvn, maxVal):
    x, y = np.mgrid[0:maxVal:.01, 0:maxVal:.01]
    pos = np.empty(x.shape + (2,))
    pos[:, :, 0] = x; pos[:, :, 1] = y
    return x, y, mvn.pdf(pos)

def computeGaussinaPDF(M):
    T = M[:, [1,2]]
    mu = np.mean(T, axis=0)
    covariance = np.cov(T,rowvar=False)
    mvn = multivariate_normal(mu,covariance) #create a multivariate Gaussian object with specified mean and covariance matrix
    return  np.column_stack((M, mvn.pdf(T))), mvn


# G1 is a matrix 1700x3, first colum id, second column x coordinates, third column y coordinates

plt.figure()
plt.hold(True)
plt.scatter(G1[:,1],G1[:,2],facecolors='none', edgecolors='g')
plt.xlabel("mean_intensity_fitc")
plt.ylabel("mean_intensity_apc")
plt.grid();

pdfG1, pfG1 = computeGaussinaPDF(G1)
G1Pdf = pdfG1[pdfG1[:,3] > 0.1,:]
x1,y1,contour1 = gaussianContour(pfG1, np.amax(G1[:,1]))
plt.contour(x1, y1, contour1, [0.1])
plt.scatter(G1Pdf[:,1],G1Pdf[:,2],facecolors='none', edgecolors='y')
points = createConfidenceEllipsePoints(G1, 0.1)
plt.scatter(points[:,0],points[:,1],color='r')
plt.show(block=True)

What am I doing wrong? Any help is greatly appreciated!