Auswertung Versuch O8 mit GnuPlot

Mit Hilfe der Fraunhoferschen Näherung kann das Kirchhoffsche Beugungsintegral auf einfache Weise gelöst werden. Dazu werden zwei Annahmen gemacht:

1. Das auf das beugende Objekt fallende Licht ist parallelen und kann daher mit ebenen Wellenfronten beschrieben werden.
2. Die Beugungserscheinungen an diesem Objekt, dessen beugende Strukturen sich innerhalb eines Kreises mit dem Radius $\rho$ befinden, werden im Fernfeld betrachtet. Dies bedeutet, dass der Beobachtungsabstand $L$ zum beugenden Objekt groß ist, dass die Fresnelzahl $\frac{\rho^2}{\lambda\,L}$ sehr viel kleiner als 1 ist.

Im weiteren wird auch davon ausgegangen, dass die beugenden Blenden und der Schirm zur Beobachtung der Beugungsbilder senkrecht zur Ausbreitungsrichtung der einfallenden Lichtwellen angeordnet sind.

Unter diesen Voraussetzungen können die für die Beugung relevanten Abmessungen, wie die Breite eines Spalts oder der Durchmesser einer Lochblende, auf einfache Weise aus den beobachteten Beugungsbilder bestimmt werden. Die Messunsicherheit, die dabei erreicht werden kann, ist deutlich kleiner als die, die mit anderen direkten optischen Messmethoden, zum Beispiel der Bestimmung mit Hilfe eines Lichtmikroskops möglich ist.

Einzelspalt

Die Intensität des an einem Einzelspalt gebeugten Lichts wird unter den oben genannten Annahmen durch die Proportionalität

$$I \propto \left(\frac{\sin(\Theta)}{\Theta}\right)^2 \Theta = \frac{ \pi B \sin \alpha}{\lambda}$$ (1)

bestimmt (Lipson, Lipson, Tannhauser, Optik (1997), Seite 175). Hierbei bezeichnet $B$ die Breite des Spalts, $\lambda$ die Wellenlänge des Lichts und $\alpha$ den Winkel zur Ausbreitungsrichtung der einfallenden ebenen Lichtwellen. Die Minima dieser Funktion liegen bei

$$\Theta = k\,\pi k = 1,2,3,\ldots$$

$$\sin \alpha = k \frac{\lambda}{B}$$

Mit Hilfe der doppelten Kleinwinkelnäherung,

$$\sin \alpha \approx \alpha \approx\tan \alpha = \frac{a}{L}$$ (2)

ergibt zwischen der Ordnung $k$ und dem Abstand $a_k$ der Minima vom zentralen Hauptmaximum ein linearer Zusammenhang

$$a_k = \frac{L\,\lambda}{B}\,k$$ (3)

Aus dem Anstieg dieser Geraden und den bekannten Größen $\lambda$ und $L$ kann die Spaltbreite berechnet werden. Diese Vorgehensweise ist immer dann zulässig, wenn der relative Fehler der Kleinwinkelnäherung deutlich kleiner ist als die relative Messunsicherheit des Abstandes

$$\frac{ \tan \alpha - \sin \alpha}{\sin \alpha} =\sqrt{1\,+\,\left... ... \approx \frac{1}{2}\left(\frac{a_{\mathrm{max}}}{L}\right)^2 < \frac{u_L}{L}$$

Diese Bedingung ist erfüllt, wenn für die Messunsicherheit von $L$ gilt:

$$u_L > \frac{a_{\mathrm{max}}^2}{2\,L}$$ (4)

wobei $a_{\mathrm{max}}$ der größte beobachtete Minimaabstand ist.

Auswertung der Minimaabstände

Auf Grund der hohen Intensität und der Breite des Hauptmaximums ist dessen genaue Lage auf dem Beobachtungsschirm nur sehr schwer zu bestimmen. Unter Ausnutzung der Symmetrie des Beugungsbilds kann dieses Problem umgangen werden. Werden die Positionen der Minima $x_{-k}$ und $x_{+k}$ auf beiden Seiten des Hauptmaximums auf einer willkürlich positionierten Längenskala bestimmt, so gilt:

$$a_k = \frac{\vert x_{+k} - x_{-k}\vert}{2}$$

Die Lage des Hauptmaximums ist für alle Werte von $k$ durch

$$x_{0_k} = \frac{x_{+k} + x_{-k}}{2}$$

bestimmt. Eine Betrachtung mit dem Fehlerfortpflanzungsgesetz ergibt, dass die Unsicherheiten der so bestimmten Werte von $a_k$ und $x_{0_k}$ identisch sind. Das bedeutet, das die Standardabweichung des Mittelwertes $\overline{x_0}$ mit der Standardabweichung der einzelnen $a_k$ Werte übereinstimmt und daher als Maß für deren Messunsicherheit benutzt werden kann. Bei dieser Art der Auswertung muss die Gerade (Gleichung 3) zwingend durch den Nullpunkt gehen. Es kann keinen physikalisch oder messtechnisch begründeten Offset geben.

Da bei dieser Messreihe nur sehr wenige Daten mit GnuPlot ausgewertet werden, kann die Tabelle mit den Werten für $k$, $x_{-k}$ und $x_{+k}$ als inline Daten direkt in die Datei mit den GnuPlot-Befehlen aufgenommen werden.

#k xlinks xrechts
$data « EOD
10 47 143
9 53 138
8 57 134
7 63 129
6 67 125
5 72 120
4 77 116
3 82 111
2 87 106
1 92 101
EOD

Die Standardabweichung der Mittelpunktslage kann neben vielen weiteren statistischen Parametern aus diesen Werten direkt berechnet werden.

stats $data using (($2+$3)/2)
n = STATS_records
mw = STATS_mean
std = STATS_ssd

Die im weiteren noch benötigten Werte, die Anzahl der Datenpunkte, der Mittelwert der Lage des Hauptmaximums $x_0$ sowie dessen Standardabweichung, werden dabei gleich in entsprechend benannten Variablen gespeichert. Nach der Definition der sehr einfachen Fitfunktion und dem sinnvollen Setzen des Anfangswertes für den einzigen Parameter $\theta_1$ kann die Fitroutine aufgerufen werden.

a(x) = theta1*x
theta1=5
fit a(x) $data using 1:(($3-$2)/2) via theta1

Dieser Fit erfolgte ohne Gewichtung, da bei dieser Auswertemethode die Messunsicherheiten aller Abstände $a_k$ gleich der Standardabweichung von $x_0$ sind. Die Berücksichtigung des systematischen Restfehlers der zur Messung der Minimalagen verwendeten Längenskala ergibt nur sehr geringe Unterschiede zwischen den einzelnen Messunsicherheiten $u_{a_k}$, so dass auch in diesem Fall auf eine Gewichtung verzichtet werden kann.

Nach dem Setzen sinnvoller Werte für die Beschriftung und den Plotbereich kann die Grafik mit den Messdaten und der berechneten Ausgleichsgeraden gezeichnet werden. Der Wert der Standardabweichung von $x_0$ wird zur Darstellung der Fehlerintervalle genutzt.

set xlabel "Ordnung"
set xrange [0:n+1]
set ylabel "a [mm]"
plot $data using 1:(($3-$2)/2):(std) with errorbars pt 7 ps 0.75\
title "measured data", a(x) title "linear fit"

Abbildung 1: Beugung am Spalt: gemessene Abstände der Minima vom Hauptmaximum und daraus berechnete Ausgleichsgerade

Aus dem Wert des Fitparameters $\theta_1$ und dessen Fehler, beide Werte werden im Terminalfenster von GnuPlot angezeigt und können der Datei fit.log entnommen werden, kann bei bekanntem Abstand $L$ und bekannter Wellenlänge $\lambda$ die Breite des beugenden Spalts

$$B = \frac{L\,\lambda}{\theta_1}$$

und durch Anwendung des Fehlerfortpflanzungsgesetzes die zugehörige Unsicherheit berechnet werden.

Auswertung der Minimalpositionen

Für die Lagen der Minima des Beugungsbilds des Spalts auf einer willkürlich positionierten Skala ergibt sich unter Anwendung der Kleinwinkelnäherung (Gleichung 2) der Zusammenhang:

$$x_k = x_0 + \frac{L\,\lambda}{B}\,k$$ (5)

Hierbei bedeutet $x_0$ die Lage des Hauptmaximums. Bei dieser Art der Auswertung gibt es keine statistisch begründete Möglichkeit, die Messunsicherheit der einzelnen $x_k$ Werte zu bestimmen, so dass diese nur unter Berücksichtigung des systematische Restfehlers der verwendeten Skala und der konkreten Messbedingungen abgeschätzt werden kann. Auch in diesem Fall kann auf eine Gewichtung verzichtet werden, da sich die ergebenden Messunsicherheiten $u_{x_k}$nur wenig unterscheiden werden.

Der für der Auswertung notwendige Datensatz enthält nur die zwei Spalten für k und $x_k$

$data « EOD
-10 47
-9 53
-8 57

$\vdots$

8 134
9 138
10 143
EOD

Wenn eine Gewichtung nicht notwendig ist, können für den Fall eines einfachen Geradenausgleiches

$$y(x) = \theta_1 + \theta_2 x$$

die Parameter $\theta_1$ und $\theta_2$ und deren Fehler mit dem Befehls

stats $data

berechnet und im Terminalfenster ausgegeben werden. Um diese Werte auch in weiteren Rechnungen zu benutzen sollten sie in geeignet benannte Variablen übernommen werden.

theta1=STATS_intercept
utheta1=STATS_intercept_err
theta2=STATS_slope
utheta2=STATS_slope_err

Für die Darstellung der erhaltene Ausgleichsgerade, der Messdaten und deren abgeschätzten Messunsicherheiten, sind noch einige weitere Größen zu definieren.

uy=0.5
y(x) = theta1+theta2*x
set xlabel "Ordnung"
set ylabel "Minima Position [mm]"

Danach kann die Graphik erstellt werden.

plot $data using 1:2:(uy) with errorbars pt 7 ps 0.75 \
title "measured data", y(x) title "linear fit"

Abbildung 2: Beugung am Spalt: gemessene Lagen der Minima und daraus berechnete Ausgleichsgerade

Wie im vorigen Abschnitt kann bei bekanntem Abstand $L$ und bekannter Wellenlänge $\lambda$ aus dem Wert des Anstiegs $\theta_2$ die Spaltbreite

$$B = \frac{L\,\lambda}{\theta_2}$$

berechnet werden. Mit Hilfe des Fehlerfortpflanzungsgesetzes ergibt sich aus den Unsicherheiten $u_L$, $u_\lambda$ und $u_{\theta_2}$ die Unsicherheit der Spaltbreite.

Lochblende

Die Intensitätsverteilung des an einer Lochblende gebeugten Lichts lässt sich in der Fraunhoferschen Näherung mit der Gleichung

$$I = I_0 \left(\frac{2\, J_1(\Theta)}{\Theta}\right)^2 \Theta = \frac{\pi\,D\,\sin \alpha}{\lambda}$$ (6)

beschreiben (Lipson, Lipson, Tannhauser, Optik (1997), Seite 178). $J_1$ bezeichnet hierbei die Besselfunktion 1.Ordnung. $D$ ist der Durchmesser der Lochblende, $\lambda$ die Wellenlänge des verwendeten Lichts und $\alpha$ den Beobachtungswinkel bezogen auf die optischen Achse.

Eine erste Abschätzung des Durchmessers der beugenden Lochblende kann mit Hilfe des Winkels des ersten Minimas $\alpha_1$ vorgenommen werden, für das gilt:

$$\sin \alpha_1 = 1.22\;\frac{\lambda}{D}$$

Wenn die Bedingung der Gleichung 4 erfüllt ist, kann auch hier wieder die Kleinwinkelnäherung

$$\sin \alpha \approx \alpha \approx \tan \alpha = \frac{x - x_0}{L}$$

angewendet werden. Der Wert $x_0$ bezeichnet die Lage des 0.Maximums auf der Messskala $x$. $L$ ist der Abstand zwischen der Lochblende und dem Eintrittsspalt der verwendeten Photodiode. Damit ergibt sich für $\Theta$ der vereinfachte Ausdruck:

$$\Theta = \frac{\pi\;D\;(x-x_0)}{\lambda \; L}$$

Die Größen $D$, $\lambda$ und $L$ können zu einem Geometrieparameter zusammengefasst werden, in den zur Vereinfachung der nächsten Schritte auch noch die Konstante $\pi$ aufgenommen wird.

$$G = \pi \frac{D}{\lambda \; L}$$

Mit dem sich ergebenden Ausdruck für $\Theta$

$$\Theta = G\;(x-x_0)$$

kann die Fitfunktion für die gemessenen Intensitätsverteilung $I(x)$ aufgestellt werden,

$$I(x) = 4 \,I_0 \left(\frac{J_1\left(G\;(x-x_0)\right)}{G\;(x-x_0)}\right)^2 + I_B$$ (7)

die als anzupassende Parameter die Größen $I_0$, $x_0$, und $G$ enthält. $I_B$ wurde als zusätzliche, von $x$ unabhängige Größe eingeführt, um das nicht zu vermeidende Restlicht zu berücksichtigen. Dieser ,,Dunkelstrom`` kann einfach durch Unterbrechung des Strahlenganges vor der Lochblende gemessen werden und liegt meistens mehr als eine Größenordnung unter dem kleinsten in der Intensitätsverteilung gemessenen Strom.

Um die Daten mit GnuPlot auszuwerten, ist als erstes die Erstellung einer Datei mit den Messdaten in tabellarischer Form erforderlich. Die erste Spalte enthält die gemessene x-Position (z.B. in mm), die zweite den gemessenen Strom (z.B. in nA). Als Beispiel soll im weiteren die Datei Messdaten-20170317-Platz4.txt verwendet werden. Das Einlesen und Darstellen der Daten aus der vorbereiteten Datei erfolgt durch Eingabe des Befehls:

set logscale y
plot "Messdaten-20170317-Platz4.txt" using 1:2

Dabei wird eine logarithmische y-Skala verwendet, damit auch die Nebenmaxima dargestellt werden. Die Definition der Intensitätsverteilung hinter der Lochblende erfolgt durch:

I(x) = (x==x0) ? I0+Ib : \
I0*(2*besj1(G*(x-x0))/(G*(x-x0)))**2+Ib

Die Stelle $x = x0$ erfordert eine besondere Behandlung, da an dieser Stelle eine Division durch 0 auftritt, die einen Abbruch der Fitroutine zur Folge hätte. Der Grenzwert

$$\lim_{x\to 0} \, \frac{J_1(x)}{x} = \frac{1}{2}$$

ist eindeutig definiert. Als nächstes werden die notwendigen numerischen Größen, die aus dem ersten Plot abgeschätzt werden können, festgelegt.

D=0.3
lambda=632.8*10**-6
l=1591
G=pi*D/(lambda*l)
x0=6.5
I0=7300
Ib=0.007

Alle Längen sind hier in $mm$ und alle Ströme in $nA$ angegeben. Dies kann gleich in die Beschriftung der Achsen aufgenommen werden.

set xlabel "x [mm]"
set ylabel "I [nA]"

Damit kann I(x) zusammen mit den Messdaten angezeigt werden.

set logscale
plot "Messdaten-20170317-Platz4.txt" using 1:2, I(x)

Durch manuelle Veränderung des unbekannten Wertes für den Blendendurchmesser kann man versuchen, die Anpassung an die Messdaten manuell weiter zu verbessern und einen guten Startwerte für den nichtlinearen Fitalgorithmus zu finden. Danach kann der erste ungewichtete Fit erfolgen.

fit I(x) "Messdaten-20170317-Platz4.txt" using 1:2 via x0,I0,G

Da es sich bei $I_B$, im Unterschied zu $x_o$ und $I_0$, um eine direkt gemessene, von $x$ unabhängige Größe handelt, sollte dieser Wert nicht als Fitparameter verwendet werden. Das Ergebnis kann mit

set logscale y
set samples 2500
plot "Messdaten-20170317-Platz4.txt" using 1:2 \
title "measured data", I(x) title "fitted curve I(x)"

dargestellt werden. Zur richtigen Darstellung der Minima wurde die Zahl der Stützstellen für die Berechnung der Fitfunktion gegenüber dem Standardwert 100 deutlich erhöht. Dafür sind in diesem Fall 2500 Stützstellen notwendig.

Abbildung 3: Beugung an einer Lochblende: gemessene Intensitätsverteilung

In dieser Darstellung fällt auf, dass die Abweichungen zwischen den Messdaten und der Fitkurve im Bereich der Minima besonders ausgeprägt sind. Ursache dafür ist die Breite des Eintrittsspaltes vor der Photodiode, die größer als Null sein muss. Dies führt dazu, dass der gemessene Strom dem Integral über einen durch die Breite diese Eintrittsspaltes bestimmten Bereich der theoretischen Intensitätsverteilung entspricht.

Berücksichtigung der Auflösungsfunktion

Um die gemessenen Intensitätsverteilung richtig zu beschreiben, ist die Auswirkung der Auflösungsfunktion $r(x)$ des Eintrittsspaltes zu berücksichtigen. Die Auflösungsfunktion $r(x)$ kann durch eine Rechteckfunktion $\Pi(\tau)$ beschrieben werden.

$$r(x) = \frac{\Pi \left(\frac{x}{2 b}\right)}{2 b} = \begin{cases... ...ght \vert \leq 1 \\ 0 & \left \vert \frac{x}{b} \right \vert > 1 \end{cases}$$

$b$ ist dabei die halbe Breite des Eintrittsspaltes. Für diese Auflösungsfunktion gilt mit $b > 0$

$$\int_{-\infty }^{\infty } r(x) \, dx = 1$$

Die gemessene Intensität $I_M$ ergibt sich aus der Faltung der durch Gleichung 7 gegebenen Intensitätsverteilung $I(x)$ mit der Auflösungsfunktion $r(x)$ des Eintrittsspaltes.

$$I_M(x) = ( I \star r)(x) = \int I(\tau) r(x - \tau) d\tau$$

Dieses Faltungsintegral kann analytisch gelöst werden. Unter Ausnutzung der Distributivität der Faltung $r \star ( g + h ) = r \star g + r \star h$ und der Assoziativität mit der skalaren Multiplikation $r \star ( a g ) = a ( r \star g )$ sowie durch Ersetzen des Terms $(x - x0)$ durch $\tilde{x}$ lässt sich obige Gleichung einfach umformen (siehe hierzu : Mathematica Notebook O8-LoesungFaltungsintegral.nb und O8-ConvolutionSpaltPhotodiode.nb). Man erhält für den gemessenen Strom :

$$I_M(x) = I_0 f(x-x_0) + I_B$$ (8)

mit

$$f(\tilde{x}) = \frac{1}{2 b}\int_{-b}^{b}\left(\frac{2 J_1(G (\tilde{x}-\tau))}{G (\tilde{x}-\tau)}\right)^2 d \tau$$

$$= \frac{4}{3 b} \left( \frac{\left( G^2 (b-\tilde{x})^2-\frac{1}{2}... ...de{x})^2-\frac{1}{2}\right) J_1(G (b+\tilde{x})){}^2}{G^2(b+\tilde{x})} \right.$$

$$+ (b-\tilde{x}) J_0(G (b-\tilde{x})){}^2 + (b+\tilde{x}) J_0(G (b+\tilde{x})){}^2$$

$$- \left.\frac{J_0(G (b-\tilde{x})) J_1(G (b-\tilde{x}))}{G} -\frac{J_0(G (b+\tilde{x})) J_1(G (b+\tilde{x}))}{G}\right)$$

Für $\vert\tilde{x}\vert = b$ kann der Wert von $f(\tilde{x})$ nicht numerisch berechnet werden. Der Grenzwert

$$\lim_{\vert\tilde{x}\vert\to b}\,f(\tilde{x}) = \frac{4}{3 b} \le... ...){}^2}{b G^2} + 2 b J_0(2 b G){}^2 - \frac{ J_1(2 b G) J_0(2 b G)}{G}\right)$$

ist definiert und kann anstelle von $f(\tilde{x})$ genutzt werden. Der Wert von $f(\tilde{x})$ an der Stelle $\tilde{x}=0$

$$f_0 = f(0) = \frac{4}{3 b} \left(\frac{ \left(2 b^2 G^2-1\right) J_1(b G){}^2}{b G^2}+2 b J_0(b G){}^2-\frac{2 J_0(b G) J_1(b G)}{G}\right)$$

ist sowohl von $G$ als auch von $b$ abhängig. Dies führt zu einer unnötigen, beim Fitten störenden Korrelation zwischen $I_0$ und $G$ sowie zwischen $I_0$ und $b$. Um dies zu vermeiden, sollte anstelle von Gleichung 8 der Ausdruck

$$I_M(x) = I_0\; \frac{f(x-x_0)}{f_0} + I_B$$ (9)

verwendet werden.

Diese Funktion wird in GnuPlot am besten aus mehrere Teilfunktionen zusammengesetzt. Zuerst werden die drei Hilfsfunktionen f0(b,G), fb(b,G) und f(x)

f0(b,G) = (2*b**2*G**2-1)*(besj1(b*G))**2/(b*G**2) \
+ 2*b*(besj0(b*G))**2-2*besj0(b*G)*besj1(b*G)/G

fb(b,G) = (2*b**2*G**2-0.25)*(besj1(2*b*G))**2/(b*G**2) \
+ 2*b*(besj0(2*b*G))**2-besj0(2*b*G)*besj1(2*b*G)/G

f(x) = (abs(x)==b)? fb(b,G) : \
(G**2*(b-x)**2-0.5)*(besj1(G*(b-x)))**2/(G**2*(b-x)) \
+ (G**2*(b+x)**2-0.5)*(besj1(G*(b+x)))**2/(G**2*(b+x)) \
+ (b-x)*(besj0(G*(b-x)))**2 \
+ (b+x)*(besj0(G*(b+x)))**2 \
- (besj0(G*(b-x)))*(besj1(G*(b-x)))/G \
- (besj0(G*(b+x)))*(besj1(G*(b+x)))/G

definiert, mit deren Hilfe dann die eigentliche Fitfunktion festgelegt wird.

I(x) = (x==x0)? I0+Ib : I0*f(x-x0)/f0(b,G)+Ib

Die weiteren Schritte sind analog der oben beschriebenen Vorgehensweise. Die dort erhaltenen Ergebnisse können als Startwerte verwendet werden. Für die Breite des Eintrittsspaltes wird ein Wert von etwa einem Millimeter angenommen.

b = 0.5
Ib=0.007
fit I(x) "Messdaten-20170317-Platz4.txt" using 1:2 via x0,I0,G,b

Zur richtigen Darstellung der Minima beim Zeichnen ist die Zahl der Stützstellen wieder deutlich zu erhöhen.

set samples 2500
set logscale y
plot "Messdaten-20170317-Platz4.txt" using 1:2 \
title "measured data", I(x) title "fitted curve I(x)"

Die graphische Darstellung zeigt, dass die Anpassung an die Daten deutlich besser ist.

Abbildung 4: Messdaten und berechnete Fitfunktion unter Berücksichtigung der Breite des Eintrittsspaltes

Dennoch weicht die berechnete Fitfunktion im Bereich der Minima immer noch deutlich von den Messwerten ab. Diese Differenzen werden bei höheren Ordnungen noch auffälliger und betreffen fast alle Punkte mit kleinen Werten des gemessenen Stromes.

Der durchgeführten Ausgleichsrechnung lag die Annahme zu Grunde, dass alle Stromwerte mit der gleichen Messunsicherheit bestimmt wurden. Daher wurde auf eine Gewichtung verzichtet. Im Anbetracht der Tatsache, dass der Wertebereich der Strommessung fünf Größenordnungen umfasst, ist diese Annahme nicht gerechtfertigt und es muss eine Gewichtung durchgeführt werden.

Gewichtung mittels Poissonverteilung

Die Intensität des an der Lochblende gebeugten Lichtes wird mit einer als Stromquelle betriebenen Silizium-Photodiode gemessen, die im Quasi-Kurzschluss betrieben wird. In dieser Betriebsart ist der Strom in Sperrrichtung über viele Größenordnungen linear von der Bestrahlungsstärke abhängig. Es gilt:

$$I = \frac{Q}{t}$$

$$= \frac{\eta\;q_e\;N_{\mathrm{Photon}}}{t}$$

$$= \frac{\eta\;q_e}{t} N_{\mathrm{Photon}}$$

$Q$ ist die Ladungsmenge, die innerhalb der Zeitspanne $t$ durch die auftreffenden Photonen in der Photodiode generiert wird. Der Wert des Quantenwirkungsgrades der Photodiode $\eta$, die Elementarladung $q_e$ und die interne Zeitkonstante $t$ des verwendeten Picoamperemeters werden im weiteren zu der Konstanten $C_I$ zusammengefasst. Damit kann man für I schreiben:

$$I = C_I\; N_{\mathrm{Photon}}$$

Die Anzahl der Photonen ist poissonverteilt. Damit gilt für die Standardabweichung des gemessenen Stromes

$$u_I = C_I\;\sqrt{N_{\mathrm{Photon}}} = \sqrt{C_I} \sqrt{I} \propto \sqrt{I}$$

wobei der Wert der Größe $\sqrt{C_I}$ nicht bekannt ist. Für die Festlegung von Gewichtsfaktoren, ist es nicht notwendig den genauen Wert der einzelnen Messunsicherheiten zu kennen, die gegebene Proportionalität ist ausreichend. Da es sich bei der hier diskutierten Gewichtung zwar um eine physikalisch begründete Abschätzung der Messunsicherheit handelt, die darin enthaltene Proportionalitätskonstante $\sqrt{C_I}$ jedoch unbekannt ist, muss die Option errorscaling gesetzt sein. Dies entspricht der Grundeinstellung von GnuPlot.

set fit errorscaling

Nun ist noch die Funktion zur Berechnung von $u_I$ zu definieren. Dabei wird für die Proportionalitätskonstante der Einfachheit halber der Wert 1 festgelegt. Damit ergibt sich:

u(I)=sqrt(I)

Die Ergebnisse des ungewichteten Fits können als sinnvolle Startwerte benutzt werden.

x0=6.45069
I0=7249.12
G=0.845147
b=0.628305
Ib=0.007

Danach kann der Fit mit den vier Parametern x0, I0, G und b gestartet und das Ergebnis zusammen mit den Messdaten dargestellt werden. Die zur Gewichtung benötigten Messunsicherheiten werden dabei aus den Messwerten berechnet.

fit I(x) "Messdaten-20170317-Platz4.txt" \
using 1:2:(u($2)) yerr via x0,I0,G,b
set logscale y
set samples 2500
plot "Messdaten-20170317-Platz4.txt" using 1:2 \
title "measured data", I(x) title "fitted curve I(x)"

Ein Plotten der Fehlerintervalle der Messdaten ist auf Grund der Tatsache, dass nur die Proportionalität zu den aus der Poisson-Verteilung abgeschätzten Unsicherheiten genutzt wurde, nicht sinnvoll.

Abbildung 5: Messdaten und berechnete Fitfunktion unter Berücksichtigung der Breite des Eintrittsspaltes mit Gewichtung nach Poisson-Verteilung

Allein schon der Vergleich mit dem ungewichteten Fit (Abbildung 4) zeigt, dass sich die Anpassung der Fitkurve an die Messdaten noch einmal deutlich verbessert hat. Ein Vergleich der Werte $\chi^2/\mathrm{DoF}$ (DoF steht für die Anzahl der Freiheitsgrade) ist hier nicht möglich, da einmal ohne und einmal mit Gewichtung gefitted wurde. Im ersten Fall entspricht der Wert $\chi^2/\mathrm{DoF}$ der Restvarianz, dem Teil der Gesamtvarianz, der nicht durch den Fit erklärt wird. Im zweiten Fall gibt der Wert $\chi^2/\mathrm{DoF}$ die Varianz der Gewichtseinheit an, die von dem bei der Gewichtung frei gewählten Proportionalitätsfaktor abhängt.

Instrumentelle Gewichtung

Wenn Angaben zum systematischen Restfehler des verwendeten Strommessgerätes bekannt sind, können diese zur Gewichtung benutzt werden. Das hier verwendete Picoamperemeter MV40 verfügt über insgesamt 16 Messbereiche (300$\mu$A, 100$\mu$A, 30$\mu$A, 10$\mu$A, 3$\mu$A, 1$\mu$A, 300nA, 100nA, 30nA, 10nA, 3nA, 1nA, 300pA, 100pA, 30pA, 10pA), die sich über mehr als 7 Größenordnungen erstrecken. Der Hersteller gibt den Gerätefehler mit 1,5% des Messbereichsendwertes an. In den Picoampere-Bereichen kommt noch eine Unsicherheit von 3pA hinzu. Geht man davon aus, das immer im jeweils kleinst möglichen Messbereich gemessen wurde, lässt sich die Unsicherheit $u_I$ der Strommessung aus dem jeweiligen Messwert ableiten.

GK = 0.015
ZF = 0.003
u(I) = (I>100000)? GK*300000 : (I>30000)? GK*100000 :\
(I>10000)? GK*30000 : (I>3000)? GK*10000 : (I>1000)? GK*3000 :\
(I>300)? GK*1000 : (I>100)? GK*300 : (I>30)? GK*100 :\
(I>10)? GK*30 : (I>3)? GK*10 : (I>1)? GK*3 : (I>0.3)? GK*1.0:\
(I>0.1)? GK*0.3+ZF : (I>0.03)? GK*0.1+ZF :\
(I>0.01)? GK*0.03+ZF : GK*0.01+ZF

Die Ergebnisse des ungewichteten Fits können auch in diesem Fall als sinnvolle Startwerte genutzt werden. Zum Starten der Fitroutine wird der gleiche GnuPlot Befehl genutzt wie im Fall der Gewichtung mittels Poisson-Verteilung.

x0=6.45069
I0=7249.12
G=0.845147
b=0.628305
Ib=0.007
fit I(x) "Messdaten-20170317-Platz4.txt" \
using 1:2:(u($2)) yerr via x0,I0,G,b

Das Ergebnis dieses Fits hängt sehr empfindlich von der Wahl der Startwerte für die Fitparameter ab. Schon kleine Änderungen können dazu führen, dass das angewendete Marquardt-Levenberg-Verfahren in einem lokalen und nicht im globalen Minimum konvergiert. In der graphischen Darstellung wird dies meist sofort deutlich. Hier helfen oftmals nur das Wissen über die physikalischen Zusammenhänge und über die genutzten statistischen Verfahren sowie Erfahrung und Inspiration weiter. Nicht um sonst heißt es im GnuPlot-Manual unter Verweis auf das fudgit Projekt: ``Nonlinear fitting is an art!''.

Da in diesem Fall abgeschätzte Messunsicherheiten genutzt wurden, können in der graphischen Darstellung diese als Fehlerbalken mit eingezeichnet werden.

set logscale y
set samples 2500
plot "Messdaten-20170317-Platz4.txt" using 1:2:(u($2)) \
with errorbars pt 7 ps 0.5 title "measured data", \
I(x) title "fitted curve I(x)"

Abbildung 6: Messdaten und berechnete Fitfunktion unter Berücksichtigung der Breite des Eintrittsspaltes mit instrumenteller Gewichtung

Auf Grund der verwendeten unterschiedlichen Methoden zur Berechnung der Gewichte, einmal auf Grundlage der erwarteten statistischen Verteilung, einmal auf Grundlage der Herstellerangaben zur Messunsicherheit des verwendeten Messgeräts, ist ein Vergleich der beiden $\chi^2/\mathrm{DoF}$ Werte nicht zulässig.

Ähnlich wie bei der Auswertung der Beugung am Spalt kann aus dem letztendlich erhaltenen Wert des Parameters $G$ bei bekanntem Abstand $L$ und bekannter Wellenlänge $\lambda$ der gesuchte Blendendurchmesser $D$ und dessen Unsicherheit berechnet werden.

Als letztes ist noch GnuPlot zu beenden:

quit

Literatur

Lipson, Stephen G., Lipson, Henry S., and Tannhauser, David S.
Optik, Mit 125 Aufgaben und vollständigen Lösungen. Optical Physics <dt.>
Springer-Lehrbuch. Berlin [u.a.]: Springer, 1997

http://www.springer.com/de/book/9783540619123 abgerufen am 11.4.2017 12:06 Uhr