clc
clear all
close all
z=0;
for xz=0:1:3
for yz=0:1:3
z=z+1;
sigmax=1;
sigmay=1;
orderx=xz;
ordery=yz;
sigma2szeRatiox = 3 + 0.25 * orderx - 2.5/((orderx-6)^2+(orderx-9)^2);
sigma2szeRatioy = 3 + 0.25 * ordery - 2.5/((ordery-6)^2+(ordery-9)^2);
sxe = ceil(sigma2szeRatiox*sigmax);
sye = ceil(sigma2szeRatioy*sigmay);
t1 = -sxe:0.1:sxe; % row vector
t1 = t1(:); % column vector
t2 = -sye:0.1:sye;
t2=t2(:);
% precalculate some orders of sigma
sigmax2 = sigmax^2;
sigmay2 = sigmay^2;
sqrt2  = sqrt(2);

% Calculate the gaussian

basegaussx =(1/sqrt(2*pi*sigmax2))* exp(- t1.^2 / (2*sigmax2) );
basegaussy =(1/sqrt(2*pi*sigmay2))* exp(- t2.^2 / (2*sigmay2) );

%  H( t/(sigma*sqrt2) ), where H() is the
% Hermite polynomial. 
x = t1 / (sigmax*sqrt2);
y = t2 / (sigmay*sqrt2);

% Depending on the order, calculate the Hermite polynomial .

switch orderx
    case 0,  part1 = 1;
    case 1,  part1 = 2*x;
    case 2,  part1 = -2 + 4*x.^2;
    case 3,  part1 = -12*x + 8*x.^3;
    case 4,  part1 = 12 - 48*x.^2 + 16*x.^4;
    case 5,  part1 = 120*x - 160*x.^3 + 32*x.^5;
end

switch ordery
    case 0,  part2 = 1;
    case 1,  part2 = 2*y;
    case 2,  part2 = -2 + 4*y;
    case 3,  part2 = -12*y + 8*y.^3;
    case 4,  part2 = 12 - 48*y.^2 + 16*y.^4;
    case 5,  part2 = 120*y - 160*y.^3 + 32*y.^5;
end

% apply Hermite polynomial to gauss
k1 = (-1).^orderx .* part1 .* basegaussx;
k2 = (-1).^ordery .* part2 .* basegaussy;


% NORMALIZE


norm_defaultx = 1 / sum(sum( basegaussx ));
norm_defaulty = 1 / sum(sum( basegaussy ));

norm_hermitex = 1/(sigmax*sqrt2)^orderx;
norm_hermitey = 1/(sigmay*sqrt2)^orderx;


k1  = k1 .* ( norm_defaultx * norm_hermitex );
k2  = k2 .* ( norm_defaulty * norm_hermitey );
k=k2*k1';

[X1,X2]=meshgrid(t1',t2');
X=[X1(:) X2(:)];
subplot(4,4,z);
contour(X1,X2,k,'LineWidth',1.5);
title(['Order=',int2str(xz+yz),'  with m=',int2str(xz),', n=',int2str(yz)]);
end
end