Self-Organized Map Gaussian Neighborhood Function

I typically use linkdist when running a SOM in MATLAB, however, the research shows that a Gaussian neighborhood function may often be the better choice. Does anyone know how to implement a Gaussian neighborhood function rather than using a built-in function? I've attached the function I use to run a SOM.

function [classescol,means,SOMx,PCAx] = bsompca(data,Hdim,Vdim,O_Lrate,T_Lrate,topology,distfunc,Zoption,retainpc,varargin)

clear net
olrate=O_Lrate;
tlrate=T_Lrate;
somrows=Hdim;
somcols=Vdim;
alldata=data;
topo=topology;
distfunction=distfunc;

retain=retainpc;
opt=Zoption;
[ostepsx,datacols]=size(alldata);

if olrate<=tlrate
    warning('O_Lrate should be greater than T_Lrate. Proceeding nonetheless...');
end

if nargin==9
    iterations=min(ostepsx*2,10000);
      osteps=round(iterations/2);
else
    iterations=varargin{1};
      osteps=round(iterations/2);
end

if strcmp(retain,'none')==1;
    PCAx{1,1}='No PCA completed';
    x = double(transpose(alldata));
    warning('to make sense of Sammon map, do subsequent PCA and use plotgrid3d function to map mean PC score by classescol');
else
    
    if opt==1
        [loadings,score,latent,~,expl] = pca(zscore(alldata));
    else
        [loadings,score,latent,~,expl] = pca(alldata);
    end

    if retain==1
        eigen=latent>1;
        retainPCnum=sum(eigen);
        retainpcs=score(:,1:retainPCnum);
        explvar=sum(expl(1:retainPCnum));
        x = double(transpose(retainpcs));
    elseif retain==0
        retainPCnum=datacols;
        retainpcs=score;
        explvar=sum(expl(1:retainPCnum));
        x = double(transpose(score));
    elseif retain>0 && retain<1
        retainPCnum=sum(cumsum(expl)<retain*100);
        retainpcs=score(:,1:retainPCnum);
        explvar=sum(expl(1:retainPCnum));
        x = double(transpose(retainpcs));
    else
        error('retainpc was entered incorrectly')
    end

    PCAx{1,1}='explvar Retained';
    PCAx{2,1}='retainPCnum';
    PCAx{3,1}='retainpcs';
    PCAx{4,1}='eigenvalues';
    PCAx{5,1}='Explaind Var all';
    PCAx{6,1}='Loadings';

    PCAx{1,2}=explvar;
    PCAx{2,2}=retainPCnum;
    PCAx{3,2}=retainpcs;
    PCAx{4,2}=latent;
    PCAx{5,2}=expl;
    PCAx{6,2}=loadings;
end

%% BEGIN SOMMing
net=selforgmap([somrows somcols],osteps,1,topo,distfunction);
net.layerWeights{1,1}.learnParam.order_lr = olrate;
net.layerWeights{1,1}.learnParam.tune_lr = tlrate;
net.trainParam.epochs = iterations;


[net2,tr] = train(net,x);
y = net2(x);
classes = vec2ind(y);

classescol=transpose(classes);
means=grpstats(alldata,classescol);
plotsompos(net2);
plotsomhits(net2,x);

somweights=net2.IW{1,1};
CPdist=dist(transpose(somweights));
DailyDist=transpose(dist(somweights,x));

posx = eval(['@' topo]);
pos=posx(somrows,somcols);

ITz=eval(['@' distfunction]);
distall=ITz(pos);
initdist=max(distall(:));

SOMx{1,1}='SOM neural network';
SOMx{2,1}='tr';
SOMx{3,1}='somweights';
SOMx{4,1}='CPdist';
SOMx{5,1}='DailyDist';
SOMx{6,1}='Settings';
SOMx{7,1}='initdist';

Settings{1,1}='Dimensions';
Settings{2,1}='O_Lrate';
Settings{3,1}='T_Lrate';
Settings{4,1}='topology';
Settings{5,1}='distfunc';
Settings{6,1}='iterations';
Settings{7,1}='PCA Settings for retain';

Settings{1,2}=[Hdim,Vdim];
Settings{2,2}=O_Lrate;
Settings{3,2}=T_Lrate;
Settings{4,2}=topology;
Settings{5,2}=distfunction;
Settings{6,2}=iterations;
Settings{7,2}=retain;

SOMx{1,2}=net2;
SOMx{2,2}=tr;
SOMx{3,2}=somweights;
SOMx{4,2}=CPdist;
SOMx{5,2}=DailyDist;
SOMx{6,2}=Settings;
SOMx{7,2}=initdist;





end