Equation of a flat top gaussian I'm trying to map a set of data points that seems to form a kind of a flat top gaussian in one dimension. I'm used to fitting data to a gaussian using non-linear least squares. I'm trying to understand how to do the same, for a flat top gaussian?
The approach is the same.  A Gaussian distribution has two parameters, the mean and standard deviation.  You add up the square of the error of each point to get an error function, which depends on those two parameters, then minimize the error function to find the best fit set of parameters.  For your case, you need to define the "flat topped Gaussian" as a function of some number of parameters.  You might just chop off the head of the Gaussian at $\pm n \sigma$, leading to something like $$M(\mu , \sigma , n,x)=\begin {cases} \frac K{\sqrt {2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma ^2}}&amp;|x| \gt n \sigma \\ \frac K{\sqrt {2 \pi \sigma^2}} e^{-\frac{(n \sigma-\mu)^2}{2\sigma ^2}}&amp;|x| \le n \sigma\end {cases}$$ where you choose $K$ as a function of $n$ so the integral over the distribution remains $1$  Now compute the error of each point from $M$, square it, add them all up, and that is your error function to minimize over the parameters.
This is my implementation in Matlab. 
I use FWHM (full width half maximum) which is the location where the Gaussian is equal to 0.5.
FunctionResult

function[y] = flatTopGaussian(sizex,FWHM,plot_result)
% Produces a flat top Gaussian.
%   Equal to 1 between -FWHM/2,FWHM/2
%   Gaussian function for the other values
% Inputs
%   sizex:  length of y
%   FWHM:   Full width half maximum
%   plot_result: should we plot the result

%%
% sizex=200;
% FWHM=2*20;
x_vec=-sizex/2+1:sizex/2;
sigma=FWHM/2.355;
[y] = 2*createGaussianPeak_1D(sizex, sigma);
y(abs(x_vec)&lt;FWHM/2)=1;

%% Plotting
if plot_result
    flat_win=zeros(size(x_vec)); flat_win(abs(x_vec)&lt;FWHM/2)=1;
    figure(1); plot(x_vec,flat_win,x_vec,y);
end
end

function [y] = createGaussianPeak_1D(sizex,output_sigma)
sxh = ceil((sizex) / 2);
mult = -0.5 / (output_sigma * output_sigma);
y=zeros(1,sizex);
for j=1:sizex
    jh = j - sxh;
    y(j) = exp(mult * (jh * jh));
end
end