Level curves don't match after rotation

[schniefen]

Consider the general 2D Gaussian function, centered at (0.5,0.5),

A*exp(-a*(-0.5 + x)**2-b*(-0.5 + x)*(-0.5 + y)-c*(-0.5 + y)**2)

where the covariance matrix can be written in terms of the coefficients a,b, and c as


2a & b
b & 2c


Rotating by 45 degrees counterclockwise gives

(a-b+c) & (a-c)
(a-c) & (a+b+c)


However, given a=1.25, b=0 and c=10000, and using Python to integrate over the unit square,

import numpy as np
import matplotlib.pyplot as plt
a=1.25
b=0   
c=10000
d=(a-b+c)/2
e=a-c
f=(a+b+c)/2
fig, ax = plt.subplots()
x,y=np.meshgrid(np.linspace(0,1,50),np.linspace(0,1,50))
z=3*np.exp(-a*(-0.5 + x)**2-b*(-0.5 + x)*(-0.5 + y)-c*(-0.5 + y)**2)
w=3*np.exp(-d*(-0.5 + x)**2-e*(-0.5 + x)*(-0.5 + y)-f*(-0.5 + y)**2) #rotated by 45 degrees counterclockwise
cs=ax.contour(x,y,z,levels=[0.8],colors='k',linestyles='dashed');
cs=ax.contour(x,y,w,levels=[0.8],colors='k',linestyles='dashed');

from scipy import integrate
h = lambda y, x: 3*np.exp(-a*(-0.5 + x)**2-b*(-0.5 + x)*(-0.5 + y)-c*(-0.5 + y)**2)
g = lambda y, x: 3*np.exp(-d*(-0.5 + x)**2-e*(-0.5 + x)*(-0.5 + y)-f*(-0.5 + y)**2)
print(integrate.dblquad(h, 0, 1, lambda x: 0, lambda x: 1))
print(integrate.dblquad(g, 0, 1, lambda x: 0, lambda x: 1))

And output:

(0.061757213121080706, 1.4742783672680448e-08)
(0.048117567144166894, 5.930455188853047e-12)

As well as the plot (where the one with coefficients a,b,c is the horizontal one, and the level curves are for C=z(x,y)=w(x,y)=0.8, both plotted over the unit square):