Matrizen-Funktion Fitten (Thema Ellipsometrie)

HecMac · Mittwoch 27. März 2024, 17:33

Hallo,
für eine Experiment muss ich eine Funktion fitten, die sich aus der Multiplikation mehrere Matrizen ergibt. (Hintergrund: Ich berechne den Brechungsindex einer Probe). Um das fitten einer solchen Funktion besser zu verstehen habe ich erstmal eine einfacheres Code-Beispiel geschrieben, dass Messdaten aus einen ähnlichen Zusammenhang generiert, wo die zu fittenden Parameter bekannt sind (Matrize A). Ziel ist es jetzt die Fit Funktion so anzuwenden, das ich auf die Parameter in etwa komme die ich als Beispiel selbst gewählt habe. Dann kann ich das Prinzip auf mein eigentliches Problem übertragen.

Hier der Code:

Code: Alles auswählen

import numpy as np
import math
import matplotlib.pyplot as plt
import cmath
from scipy.optimize import curve_fit

#Eingangsintensität in Watt
I0=1

#Erstelle Jonesvektor, linear 45°: ξ = −45°, χ = 180° Gemäß Seite 3
Jin = np.array([[1],
                [1]])*(1/np.sqrt(2))

#Winkel des rotierenden Bauteils Alpha
Alpha_deg = 10*np.arange((360/10)+1) #[0,45,90] #np.arange(3)

#Analysatormatrize A
A = np.array([[1.8, 2.2],
              [3.7, 4.5]])

#Intensitätsverlauf Array erstellt
Iv = np.array([])

for Alpha in Alpha_deg:
    #Drehmatrize erstellen für Transformation in den Analysator hinein
    D = np.array([[ np.cos(np.radians(Alpha)), np.sin(np.radians(Alpha))],
                  [-np.sin(np.radians(Alpha)), np.cos(np.radians(Alpha))]])
    
    #Drehmatrize erstellen für Transformation aus den Analysator heraus
    D2 = np.array([[ np.cos(np.radians(-Alpha)), np.sin(np.radians(-Alpha))],
                   [-np.sin(np.radians(-Alpha)), np.cos(np.radians(-Alpha))]])
    
    #Reslutierender Jonesvektor
    Jout=D2@A@D@Jin
    
    #Resultierende Intensität = Betragsquadrat von Jout * Eingangsintensität I0
    I = (np.sum(np.abs(Jout)**2))*I0 
    
    #Intensitätsverlaufs Array ein Element hinzufügen
    Iv = np.append(Iv, I)  
    

# Plot des Intensitätsverlaufs
#plt.plot(Alpha_deg, Iv)
#plt.title('Intensitätsverlauf in Abhängigkeit von Alpha')
#plt.xlabel('Alpha (Grad)')
#plt.ylabel('Intensität')
#plt.grid(True)
#plt.show()

    
Messdaten=np.array([[0,0], 
             [0,0], 
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0],
             [0,0]])
     
for i in range(0,37,1):
     Messdaten[i][0]=Alpha_deg[i]*1000
     Messdaten[i][1]=Iv[i]*1000
Messdaten=Messdaten*0.001     
    
#Ab hierbeginnt mein Versuch dieParameter a, b, c, d "her zu fitten" die in der oberen Matrize A enthalten sind    
def Funktion(Alpha, a, b, c, d):
    #Eingangsintensität in Watt
    I0=1

    #Erstelle Jonesvektor, linear 45°: ξ = −45°, χ = 180° Gemäß Seite 3
    Jin = np.array([[1],
                    [1]])*(1/np.sqrt(2))

    #Analysatormatrize A
    A = np.array([[a, b],
                  [c, d]])
    
    #Drehmatrize erstellen für Transformation in den Analysator hinein
    D = np.array([[ np.cos(np.radians(Alpha)), np.sin(np.radians(Alpha))],
                  [-np.sin(np.radians(Alpha)), np.cos(np.radians(Alpha))]])
    
    #Drehmatrize erstellen für Transformation aus den Analysator heraus
    D2 = np.array([[ np.cos(np.radians(-Alpha)), np.sin(np.radians(-Alpha))],
                   [-np.sin(np.radians(-Alpha)), np.cos(np.radians(-Alpha))]])
    
    #Reslutierender Jonesvektor
    Jout=D2@A@D@Jin
    
    #Resultierende Intensität = Betragsquadrat von Jout * Eingangsintensität I0
    I = (np.sum(np.abs(Jout)**2))*I0 
    
    return I


#Dieser Fit funktioniert noch nicht...
Fit=curve_fit(Funktion, Messdaten[:,0], Messdaten[:,1], [1.8,2.2,3.7,4.5])

Die Parameter die ich finden möchte sind in der oben in der Matrize A enthalten. Wenn ich es schaffe mit der Fit Funktion diese wieder zu erhalten ist mein Ziel erreicht. Leider erhalte ich immer folgenden Fehler:

...
File c:\...\untitled10.py:117 in Funktion
Jout=D2@A@D@Jin

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 37)

Was muss ich am Code anpassen?

Danke !

Sirius3 · Donnerstag 28. März 2024, 08:50

`math`, `cmath` und `matplotlib` werden importiert, aber nicht benutzt. Das kann also weg.
Zuerst einmal solltest Du sprechendere Variablennamen benutzen, Variabelnamen werden Python generell komplett klein geschrieben. Jin wäre besser jones_vector und Alpha_deg besser angels_in_degrees, wobei hier das `s` wichtig ist, denn es sind ja mehrere!
Wenn man einen Zahlenbereich abdecken will, dann benutzt man die passende Step-Size: `np.arange(0, 360.1, 10.0)` oder gleich `np.linspace(0, 360, 37, endpoint=True)`
Rotationsmatrizen sind Unitäre Matrizen, die Inverse ist einfach die transponierte, so dass nur eine `rotation_matrix` nötig ist und nicht `D` und `D2`.
np.append darf man nicht verwenden, weil das viel Speicher hin und herkopiert. Statt dessen benutzt man eine einfache Liste und erzeugt nach der Schleife daraus ein numpy-Array.

Code: Alles auswählen

intensities = []
for angle in angles:
    #Drehmatrize erstellen für Transformation in den Analysator hinein
    cos_alpha = np.cos(np.radians(angle))
    sin_alpha = np.sin(np.radians(angle))
    rotation_matrix = np.array([[cos_alpha, sin_alpha], [-sin_alpha, cos_alpha]])
    
    #Reslutierender Jonesvektor
    jones_out = rotation_matrix.T @ analysator @ rotation_matrix @ jones_in

    #Resultierende Intensität = Betragsquadrat von Jout * Eingangsintensität I0
    intensity = (abs(jones_out) ** 2).sum() * initial_intensity
    intensities.append(intensity)

intensities = np.array(intensities)

Wenn man bei numpy mit Schleifen über Array-Einträge arbeitet, dann macht man etwas falsch:

Code: Alles auswählen

measurement_data = np.array([intensities, angles]).T

Was ist der Sinn, erst etwas mit 1000 zu multiplizieren um es gleich danach wieder durch 1000 zu teilen???

Eine Funktion wird wie Variablen komplett klein geschrieben, und hat einen sprechenden Namen, `Funktion` erfüllt keine der Anforderungen.

Zum eigentlichen Problem: hast Du die Dimensionen der Drehmatrix angeschaut? Hast Du verstanden, warum die Dimensionen so sind? Hast Du Dir auf dem Papier die Operation mal aufgeschrieben, die Durchgeführt werden soll?