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):

[schniefen]

This has been solved. The curves do match, it is just that, when the function is rotated, its central axes changes, and in the square, it is longer than that of the non-rotated one.