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
\begin{pmatrix} 2a & b\\ b & 2c \end{pmatrix}
Rotating by 45 degrees counterclockwise gives
$$\frac{2}{4}\left( \begin{array}{rr} 1 & -1\\ 1 & 1\end{array}\right)\left(\begin{array}{cc} 2a & b\\ b & 2c\end{array}\right)\left( \begin{array}{rr} 1 & 1\\ -1 & 1\end{array}\right)= \cdots =\left(\begin{array} ((a-b+c) & (a-c)\\(a-c) & (a+b+c)\end{array}\right)$$
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 (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):
