public variable Norm, Center, Sigma;
Norm = 100.0;
Center = 12.0;  % Angstrom
Sigma = 0.025;  % Angstrom

% Generate fake data using a Gaussian line profile

public define init_data (lo, hi)
{
   variable x = struct 
     {
	bin_lo, bin_hi, value, err
     };
   
   x.bin_lo = @lo;
   x.bin_hi = @hi;
   
   fit_fun ("gauss(1)");
   set_par (1, Norm);
   set_par (2, Center, 0, 6, 24);
   set_par (3, Sigma, 0, 0.0025, 0.25);
   
   x.value = eval_fun (lo, hi);

   variable i;
   
   foreach ([0:length(lo)-1])
     {
	i = ();
	x.value[i] = prand(x.value[i]);
     }

   x.err = sqrt(x.value);
   
   return x;
}

    % Generate random datasets until we find one
    % whose best-fit center is extremely close to the
    % center of the underlying Gaussian.
    
public define generate_ideal_random_dataset (lo, hi)
{
   variable d, k = 0;
   
   delete_data (all_data);
   
   forever
     {
	d = init_data (lo, hi);
	
	if (k == 0)
	  () = define_counts (d);
	else
	  put_data_counts (1, d);
	
	() = fit_counts();
	
	k++;
	
	if (abs(get_par(2) - Center) < 1.e-5)
	  break;
     }
}

    % Generate N random datasets, fitting a Gaussian to each
    % then plot the distribution of best-fit center values.

public define fit_random_datasets (lo, hi, N)
{
   variable k, d, center= Double_Type[N];

   delete_data (all_data);
   
   _for (0, N-1, 1)
     {
	k = ();

	d = init_data (lo, hi);
	
	if (k == 0)
	  () = define_counts (d);
	else
	  put_data_counts (1, d);
	
	() = fit_counts();
	center[k] = get_par(2);
     }

   variable id = plot_open ("ex_distribution.ps/cps");
   label (latex2pg ("\\lambda_{center} [\\A]"), "N", "Distribution of Best-Fit Values");

   variable xlo, xhi;
   (xlo, xhi) = linear_grid (11.98, 12.02, 128);
   limits;
   hplot (xlo, xhi, histogram (center,xlo,xhi));
   
   plot_close (id);
}

% Generate a wavelength grid and a fake dataset:

public variable Lo, Hi;
(Lo, Hi) = linear_grid (10,14,256);
() = define_counts (init_data(Lo,Hi));