examp.m

这是在网上下的一个东东

%
%  This script runs the Backus-Gilbert computations for 
%
clear
%
% The data.
%
d=[5.973; 8.02];
sigma=[0.0005; 0.005];
Re=6370.8;
%
% Get the normalizing constraint.
%
q=[1.083221147; 1.757951474];
%
% First, try ri=0.15;
%
%
% Get the H matrix.  The following formulas for the integrals were done with
% Maple.
%
ri=1000/Re;
H=zeros(2,2);
H(1,1)=1.508616069 - 3.520104161*ri + 2.112062496*ri^2;
H(1,2)=3.173750352 - 7.140938293*ri + 4.080536168*ri^2;
H(2,1)=H(1,2);
H(2,2)=7.023621326 - 15.45196692*ri + 8.584426066*ri^2;
%
% Get the c coeffcients.
%
[c,lambda]=quadlin(H,q',[1.0]);
fprintf('Coefficients for ri=%f (%f km) are %f, %f\n',[ri,ri*Re,c(1),c(2)]);
fprintf('Estimate for ri=%f is %f \n',[ri,c'*d]);
fprintf('Standard deviation for ri=%f is %f\n',[ri,sqrt(c(1)^2*sigma(1)^2+c(2)^2*sigma(2)^2)]);
%
% Plot the averaging kernel.
%
figure(1);
bookfonts
r=(0.0:0.02:1.0)';
a=c(1)*g1(r)+c(2)*g2(r);
plot(r*Re,a,'k');
xlabel('radius (km)');
ylabel('A(r)');
H=text(100,1,'r=1000 km');
set(H,'FontSize',18);     
hold on
%
% Finally, try ri=0.8
%
%
% Get the H matrix.  The following formulas for the integrals were done with
% Maple.
%
ri=5000/Re;
H=zeros(2,2);
H(1,1)=1.508616069 - 3.520104161*ri + 2.112062496*ri^2;
H(1,2)=3.173750352 - 7.140938293*ri + 4.080536168*ri^2;
H(2,1)=H(1,2);
H(2,2)=7.023621326 - 15.45196692*ri + 8.584426066*ri^2;
%
% Get the c coeffcients.
%
[c,lambda]=quadlin(H,q',[1.0]);
fprintf('Coefficients for ri=%f (%f km) are %f, %f\n',[ri,ri*Re,c(1),c(2)]);
fprintf('Estimate for ri=%f is %f \n',[ri,c'*d]);
fprintf('Standard deviation for ri=%f is %f\n',[ri,sqrt(c(1)^2*sigma(1)^2+c(2)^2*sigma(2)^2)]);
%
% Plot the averaging kernel.
%
r=(0.0:0.02:1.0)';
a=c(1)*g1(r)+c(2)*g2(r);
plot(r*Re,a,'k');
H=text(3000,0.55,'r=5000 km');
set(H,'FontSize',18);     
%
%  Print out the results.
%
%print -deps akernels.eps