function yy = gc_deconv(data_y)
% deconvolute glowcurve

threshold = 0.1;    % stop deconvoluting when the peak intensity is this
                    % fraction of the maximum peak
E = 1.671;          % trap activation energy in eV

% assume data_y is a wellformatted source of data with length 701

yy = zeros(1,701); % this will hold the sum of the deconvoluted functions

axis([350 700]);
hold off;
plot(data_y);
hold on;

% assume peak 5 is the largest, unannealed peaks will not be subtracted
% first, therefore the low temperature part is deconvolutet frst.

%data_lt = data_y(1:533)
a = zeros(1, 168);
data_lt = [data_y(1:533)', a]; % expand data_lt to full size
data_lt = data_lt';

a = max(data_lt); % this tells how much signal there is in the LT spectrum
threshold = a * threshold; % the threshold level is based on the LT
                           % region alone.


while (a > threshold) && (a > 1)
  %  disp("-------");
  [a,b] = max(data_lt);        % a holds intensity, b temperature
                                 %  disp ("Peak (T,I):"), disp (b), disp(a);
                                 % but we do not want to include the data
                                 % from the end point at 533
  if b < 533
    [x,y] = gc_kernel(a,b,E);   % guess the peak
    data_lt = data_lt - y';       % subtract the peak
    plot(y,'b--');
    printf(" %.2f %.2f %.2f \n", b, a, sum(y));
    yy = yy + y;
    a = 0; % escape this loop
  end
end
  
% now deconvolve the HT part.

data_y = data_y - yy'; % subtract allread found peaks from data set
a = max(data_y);

while (a > threshold) && (a > 1)
%  disp("-------");
  [a,b] = max(data_y);        % a holds intensity, b temperature
%  disp ("Peak (T,I):"), disp (b), disp(a);
  [x,y] = gc_kernel(a,b,E);   % guess the peak
  data_y = data_y - y';       % subtract the peak
  plot(y,'b--');
  printf(" %.2f %.2f %.2f \n", b, a, sum(y));
  yy = yy + y;
end
hold off;


%plot(yy);
%hold on;
%data_y = data_y +yy';
%plot(data_y);


chi = norm(yy'-data_y);

