loading
%建立函数myfun求解参数的初始值 function F = myfun(x) clc
x1=[2883 3692 7 41.7] x2=[4742 7293 9 27.8] x3=[4948 7293 6 34.8] x4=[5006 8846 6 30.3] F =
[x(4)/(4*pi*((0.16*(x1(1)-x(1))*(1+0.0004*(x1(1)-x(1)))^(-0.5)))*(0.14*(x1(1)-x(1))*(1+0.003*(x1(1)-x(1)))^(-0.5)))...
*exp(-(x1(2)-x(2))^2/(2*((0.16*(x1(1)-x(1))*(1+0.0004*(x1(1)-x(1)))^(-0.5)))^2)-(x1(3)-x(3))^2/(2*(0.14*(x1(1)-x(1))*(1+0.003*(x1(1)-x(1)))^(-0.5))^2))-x1(4),
x(4)./(4*pi*((0.16*(x2(1)-x(1)).*(1+0.0004*(x2(1)-x(1))).^(-0.5))).*(0.14*(x2(1)-x(1)).*(1+0.003*(x2(1)-x(1))).^(-0.5)))...
.*exp(-(x2(2)-x(2)).^2./(2*((0.16*(x2(1)-x(1)).*(1+0.0004*(x2(1)-x(1))).^(-0.5)))^2)-(x2(3)-x(3)).^2./(2*(0.14*(x2(1)-x(1)).*(1+0.003*(x2(1)-x(1))).^(-0.5))^2))-x2(4),
x(4)./(4*pi*((0.16*(x3(1)-x(1)).*(1+0.0004*(x3(1)-x(1))).^(-0.5))).*(0.14*(x3(1)-x(1)).*(1+0.003*(x3(1)-x(1))).^(-0.5)))...
.*exp(-(0-x(2)).^2./(2*((0.16*(x3(1)-x(1)).*(1+0.0004*(x3(1)-x(1))).^(-0.5)))^2)-(16-x(3)).^2./(2*(0.14*(x3(1)-x(1)).*(1+0.003*(x3(1)-x(1))).^(-0.5))^2))-x3(4),
x(4)./(4*pi*((0.16*(x4(1)-x(1)).*(1+0.0004*(x4(1)-x(1))).^(-0.5))).*(0.14*(x4(1)-x(1)).*(1+0.003*(x4(1)-x(1))).^(-0.5)))...
.*exp(-(x4(2)-x(2)).^2./(2*((0.16*(x4(1)-x(1)).*(1+0.0004*(x4(1)-x(1))).^(-0.5)))^2)-(x4(3)-x(3)).^2./(2*(0.14*(x4(1)-x(1)).*(1+0.003*(x4(1)-x(1))).^(-0.5))^2))-x4(4)]
%x0=[-5 -5 -5 -5]
%options=optimset('Display','iter'); %[x,fval] = fsolve(@myfun,x0,options)
%对数据回归分析
function f=curvefun1(x,T) x1=T(:,1);y=T(:,2);z=T(:,3);
sigmay=0.16*(x1-x(1)).*(1+0.0004*(x1-x(1))).^(-0.5); sigmaz=0.14*(x1-x(1)).*(1+0.0003*(x1-x(1))).^(-0.5);
f=x(4)./(4*pi*sigmay.*sigmaz).*exp(-(y-x(2)).^2./(2*sigmay)-(z-x(3)).^2./(2*sigmaz))
%[beta，r，J]=nlinfit（T,y,’curvefun1’,x0）