Plot the Gaussian Distribution in 3D In probability theory, the normal (or Gaussian) distribution is a very common continuous probability distribution. Normal distributions are important in statistics and are often used in the natural and social sciences to represent real-valued random variables whose distributions are not known.


The challenge

Your challenge is to plot the probability density of the Gaussian Distribution on a 3-dimensional plane. This function is defined as:



Where:




A = 1, σx = σy = σ

Rules


Your program must take one input σ, the standard deviation.
Your program must print a 3D plot of the Gaussian Distribution in the highest quality as your language/system allows.
Your program may not use a direct Gaussian Distribution or probability density builtin.
Your program does not have to terminate.
Your plot may be in black and white or color.
Your plot must have grid lines on the bottom. Grid lines on the sides (as shown in the examples) are unnecessary.
Your plot does not need to have line numbers next to the grid lines.


Scoring

As usual in code-golf, the submission with the least bytes wins! I may never "accept" an answer using the button, unless one is incredibly small and intuitive.

Example output

Your output could look something like this:



Or it could look like this:



More valid outputs. Invalid outputs.
C++, 3477 3344 bytes

Byte count does not include the unnecessary newlines.
MD XF golfed off 133 bytes.

There's no way C++ can compete for this, but I thought it would be fun to write a software renderer for the challenge.  I tore out and golfed some chunks of GLM for the 3D math and used Xiaolin Wu's line algorithm for rasterization.  The program outputs the result to a PGM file named g.



#include&lt;array&gt;
#include&lt;cmath&gt;
#include&lt;vector&gt;
#include&lt;string&gt;
#include&lt;fstream&gt;
#include&lt;algorithm&gt;
#include&lt;functional&gt;
#define L for
#define A auto
#define E swap
#define F float
#define U using
U namespace std;
#define K vector
#define N &lt;&lt;"\n"
#define Z size_t
#define R return
#define B uint8_t
#define I uint32_t
#define P operator
#define W(V)&lt;&lt;V&lt;&lt;' '
#define Y template&lt;Z C&gt;
#define G(O)Y vc&lt;C&gt;P O(vc&lt;C&gt;v,F s){vc&lt;C&gt;o;L(Z i=0;i&lt;C;++i){o\
[i]=v[i]O s;}R o;}Y vc&lt;C&gt;P O(vc&lt;C&gt;l, vc&lt;C&gt;r){vc&lt;C&gt;o;L(Z i=0;i&lt;C;++i){o[i]=l[i]O r[i];}R o;}
Y U vc=array&lt;F,C&gt;;U v2=vc&lt;2&gt;;U v3=vc&lt;3&gt;;U v4=vc&lt;4&gt;;U m4=array&lt;v4,4&gt;;G(+)G(-)G(*)G(/)Y F d(
vc&lt;C&gt;a,vc&lt;C&gt;b){F o=0;L(Z i=0;i&lt;C;++i){o+=a[i]*b[i];}R o;}Y vc&lt;C&gt;n(vc&lt;C&gt;v){R v/sqrt(d(v,v));
}v3 cr(v3 a,v3 b){R v3{a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]};}m4 P*(
m4 l,m4 r){R{l[0]*r[0][0]+l[1]*r[0][1]+l[2]*r[0][2]+l[3]*r[0][3],l[0]*r[1][0]+l[1]*r[1][1]+
l[2]*r[1][2]+l[3]*r[1][3],l[0]*r[2][0]+l[1]*r[2][1]+l[2]*r[2][2]+l[3]*r[2][3],l[0]*r[3][0]+
l[1]*r[3][1]+l[2]*r[3][2]+l[3]*r[3][3]};}v4 P*(m4 m,v4 v){R v4{m[0][0]*v[0]+m[1][0]*v[1]+m[
2][0]*v[2]+m[3][0]*v[3],m[0][1]*v[0]+m[1][1]*v[1]+m[2][1]*v[2]+m[3][1]*v[3],m[0][2]*v[0]+m[
1][2]*v[1]+m[2][2]*v[2]+m[3][2]*v[3],m[0][3]*v[0]+m[1][3]*v[1]+m[2][3]*v[2]+m[3][3]*v[3]};}
m4 at(v3 a,v3 b,v3 c){A f=n(b-a);A s=n(cr(f,c));A u=cr(s,f);A o=m4{1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1};o[0][0]=s[0];o[1][0]=s[1];o[2][0]=s[2];o[0][1]=u[0];o[1][1]=u[1];o[2][1]=u[2];o[0]
[2]=-f[0];o[1][2]=-f[1];o[2][2]=-f[2];o[3][0]=-d(s,a);o[3][1]=-d(u,a);o[3][2]=d(f,a);R o;}
m4 pr(F f,F a,F b,F c){F t=tan(f*.5f);m4 o{};o[0][0]=1.f/(t*a);o[1][1]=1.f/t;o[2][3]=-1;o[2
][2]=c/(b-c);o[3][2]=-(c*b)/(c-b);R o;}F lr(F a,F b,F t){R fma(t,b,fma(-t,a,a));}F fp(F f){
R f&lt;0?1-(f-floor(f)):f-floor(f);}F rf(F f){R 1-fp(f);}struct S{I w,h; K&lt;F&gt; f;S(I w,I h):w{w
},h{h},f(w*h){}F&amp;P[](pair&lt;I,I&gt;c){static F z;z=0;Z i=c.first*w+c.second;R i&lt;f.size()?f[i]:z;
}F*b(){R f.data();}Y vc&lt;C&gt;n(vc&lt;C&gt;v){v[0]=lr((F)w*.5f,(F)w,v[0]);v[1]=lr((F)h*.5f,(F)h,-v[1]
);R v;}};I xe(S&amp;f,v2 v,bool s,F g,F c,F*q=0){I p=(I)round(v[0]);A ye=v[1]+g*(p-v[0]);A xd=
rf(v[0]+.5f);A x=p;A y=(I)ye;(s?f[{y,x}]:f[{x,y}])+=(rf(ye)*xd)*c;(s?f[{y+1,x}]:f[{x,y+1}])
+=(fp(ye)*xd)*c;if(q){*q=ye+g;}R x;}K&lt;v4&gt; g(F i,I r,function&lt;v4(F,F)&gt;f){K&lt;v4&gt;g;F p=i*.5f;F
q=1.f/r;L(Z zi=0;zi&lt;r;++zi){F z=lr(-p,p,zi*q);L(Z h=0;h&lt;r;++h){F x=lr(-p,p,h*q);g.push_back
(f(x,z));}}R g;}B xw(S&amp;f,v2 b,v2 e,F c){E(b[0],b[1]);E(e[0],e[1]);A s=abs(e[1]-b[1])&gt;abs
(e[0]-b[0]);if(s){E(b[0],b[1]);E(e[0],e[1]);}if(b[0]&gt;e[0]){E(b[0],e[0]);E(b[1],e[1]);}F yi=
0;A d=e-b;A g=d[0]?d[1]/d[0]:1;A xB=xe(f,b,s,g,c,&amp;yi);A xE=xe(f,e,s,g,c);L(I x=xB+1;x&lt;xE;++
x){(s?f[{(I)yi,x}]:f[{x,(I)yi}])+=rf(yi)*c;(s?f[{(I)yi+1,x}]:f[{x,(I)yi+1}])+=fp(yi)*c;yi+=
g;}}v4 tp(S&amp;s,m4 m,v4 v){v=m*v;R s.n(v/v[3]);}main(){F l=6;Z c=64;A J=g(l,c,[](F x,F z){R
v4{x,exp(-(pow(x,2)+pow(z,2))/(2*pow(0.75f,2))),z,1};});I w=1024;I h=w;S s(w,h);m4 m=pr(
1.0472f,(F)w/(F)h,3.5f,11.4f)*at({4.8f,3,4.8f},{0,0,0},{0,1,0});L(Z j=0;j&lt;c;++j){L(Z i=0;i&lt;
c;++i){Z id=j*c+i;A p=tp(s,m,J[id]);A dp=[&amp;](Z o){A e=tp(s,m,J[id+o]);F v=(p[2]+e[2])*0.5f;
xw(s,{p[0],p[1]},{e[0],e[1]},1.f-v);};if(i&lt;c-1){dp(1);}if(j&lt;c-1){dp(c);}}}K&lt;B&gt; b(w*h);L(Z i
=0;i&lt;b.size();++i){b[i]=(B)round((1-min(max(s.b()[i],0.f),1.f))*255);}ofstream f("g");f 
W("P2")N;f W(w)W(h)N;f W(255)N;L(I y=0;y&lt;h;++y){L(I x=0;x&lt;w;++x)f W((I)b[y*w+x]);f N;}R 0;}



l is the length of one side of the grid in world space.
c is the number of vertices along each edge of the grid.
The function that creates the grid is called with a function that takes two inputs, the x and z (+y goes up) world space coordinates of the vertex, and returns the world space position of the vertex.
w is the width of the pgm
h is the height of the pgm
m is the view/projection matrix. The arguments used to create m are...


field of view in radians
aspect ratio of the pgm
near clip plane
far clip plane
camera position
camera target
up vector



The renderer could easily have more features, better performance, and be better golfed, but I've had my fun!