% PROGRAM runBG 
%-----------------------------------------------------------------------
% The program demonstrates the use of the Backus-Gilbert method for
% inverting the Laplace transform from noisy data.
%-----------------------------------------------------------------------
% CALLS TO : BG.m
% 25/1/94 Erkki Somersalo 
%-----------------------------------------------------------------------

% Data points s:

N = 20;
alpha = -2; beta = 3;
s = exp(alpha + (beta - alpha)/(N-1)*(0:N-1));

% Reconstruction points t:

M = 10;
a = 0.01; b = 18;
t = a + (b - a)/M*(0:M-1);

% Plot points tplot

Mplot = 100;
tplot = a + (b - a)/(Mplot-1)*(0:Mplot-1);
PP = exp(-tplot'*s);

% Matrices S, C and g

nn = 2;  %order of the weight function
aux = 1 ./(s'*ones(1,N) + ones(N,1)*s);
g = 1 ./s';
sigma2 = input('Anna kohinan varianssi: ');
sigma = sqrt(sigma2);
C = sigma2*eye(N);
gamma2 = input('Anna rekonstruktion varianssi: ');

Phi = [];
for m = 1:M
  S = aux;
  for ell = 1:nn
   S = (-t(m))^ell*aux + ell*aux.*S;
  end
  phi = BG(S,C,g,gamma2);
  Phi = [Phi,phi];
end

% Plotting the approximating kernels

A = PP*Phi; A = A(:,2:length(A(1,:)));
figure(1)
plot(tplot,A,'r'), grid, title('Approksimoivat ytimet')

% Generating the data

choice = input('Linear (1) or hat (2) function?  ');

if choice == 2
 a = 4; bb = 8;
 g1 = 1/a; g2 = -1/(bb-a);
 c = bb/(bb-a);
 g = g1*(1 - exp(-a*s))./s./s - g1*a*exp(-a*s)./s;
 g = g + g2*(a*exp(-a*s) - bb*exp(-bb*s))./s;
 g = g + (g2+c*s).*(exp(-a*s) - exp(-bb*s))./s./s;
 ftrue = g1*tplot.*(tplot<=a) + (g2*tplot +c).*(tplot>a & tplot<=bb);
 fftrue = g1*t.*(t<=a) + (g2*t +c).*(t>a & t<=bb);
else
 g = 1 ./(s.*s);
 ftrue = tplot;
 fftrue = t;
end

figure(2)
plot(tplot,ftrue,'r'),hold on
plot(t,fftrue,'ro')
S1 = sprintf('Kohinan varianssi = %5.3f, ',sigma^2);
S2 = sprintf('rekonstruktion varianssi = %5.3f',gamma2);
for nn = 1:10
  gg = g + sigma*randn(size(g));
  f = gg*Phi;
  plot(t,f,'r'), hold on
end
grid, title([S1,S2]);
hold off




