fit_spotsize.nb

Least squares fit to determine the gaussian beam width (spotsize)

Version of 12 Apr 2010 (Sasa)

1. Data

powers={56.7,56.6,56.6,56.5,56.4,56.3,56.1,55.7,55.0,54.1,52.7,51.2,49.1,46.3,42.9,39.0,34.8,30.7,26.2,21.8,17.8,14.8,11.5,8.2,6.2,4.4,3.1,2.2,1.6,1.1,0.8,0.6,0.5,0.4,0.3,0.3, 0.3,0.3};

knifehights=50Range@Length@powers//N;

data=Transpose@{knifehights,powers};

n=Length[data]

2. Model

params={p0,p1,h0,w}

modelfn[x_,{p0_,p1_,h0_,w_}]:=p0+(p1/2)(1-Erf[(Sqrt[2])(x-h0)/w])

s2fn[data_,params_]:=Module[{m,d,x=First/@data,y=Last/@data},
m=modelfn[#,params]&/@x;
d=y-m;
d.d
]

m=Length[params]

3. Fit

pfit=FindFit[data,{modelfn[x,params],
0<p0<10,0<p1<100,500<h0<1500,100<w<1000
},params,x]

D[s2fn[data,params],#]&/@params/.pfit//N

Plot[modelfn[x,params/.pfit],{x,Min[First/@data],Max[First/@data]},PlotStyle->Thick,Frame->True,Epilog->{Red,AbsolutePointSize[5],Point/@Transpose@{First/@data,Last/@data}}]

4. Parameter errors

s2min=s2fn[data,params/.pfit]

rms=(Sqrt[(s2min/(n-m))])

atamat=(1/2)Outer[D[((n-m)/s2min) s2fn[data,params],#1,#2]&,params,params]/.pfit//N

covmat=Inverse[atamat]

pstderr=Table[covmat[[i,i]],{i,m}]//Sqrt

5. Summary

MapThread[Print[StringJoin[#1," = ",#2," +/- ",#3]]&,Map[ToString,{params,params/.pfit,pstderr},{2}]];