fit_malus.nb

Least squares fit to the Malus' law.

Version of 12 Apr 2010 (Sasa)

1. Data

angles=((2π)/360) Range[0,180,10]//N

amplitudes={96.3,75.8,53.5,32.7,16.5,5.1,1.5,5.2,15.5,32.3,53.0,74.4,95.0,112.7,123.7,128.9,124.2,113.6,97.0};

data=Transpose@{angles,amplitudes};

n=Length[data]

2. Model

params={pperp,ppar,x0}

modelfn[x_,{pperp_,ppar_,x0_}]:=pperp+ppar ((Cos[x-x0])^2)

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

FindFit[data,{modelfn[x,params],pperp≥0,ppar≥0,0≤x0<2π},params,x]

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

{a,b,c}=PseudoInverse[{1,Cos[2#],Sin[2#]}&/@(First/@data)].Last/@data;
pfit={pperp->a-(Sqrt[(b^2)+(c^2)]),ppar->2(Sqrt[(b^2)+(c^2)]),x0->ArcTan[b,c]/2}

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}]];