# -*- coding: utf-8 -*- """ Created on Mon Mar 1 11:17:37 2021 @author: Ehresmann """ #import sys import numpy as np #import pandas as pd import matplotlib.pyplot as plt #import matplotlib.animation as ani #import glob # Videobibliothek importieren import cv2 import os import h5py # Berechnet die effektive Feldlandschaft für sphärische Partikel # in einem definierten Abstand zum Substrat über einem parallele- # Streifen Domänenmuster in EB-Schichtsystemen und daraus die # dynamisch variierende Landschaft der potenziellen Energie für # superparamagnetische Partikel. Es wird die trigonometrische # Näherung genutzt, in der die Streufeldlandschaft des Substrates # und die Form der von außen angelegten Magnetfeldpulse als # trigonometrische Funktionen (sin, cos) angenähert werden. # ---------------------------------------------------- # -----------------Dateneingabe----------------------- # ---------------------------------------------------- # ------------Substrat und Domänenmuster-------------- # Magnetische Feldkonstante in Vs/(Am) mu_0 = 1.257e-6 # Dicke des Ferromagneten in Mikrometern t_FM = 10.0e-3 # Sättigungsmagnetisierung des Ferromagneten in kA/m hier Co70Fe30 # (nicht Co30Fe70) nach Zingsem et al. JPD 50, 495006 (2017) und # nach Bardos et al. JAP 40, 1371 (1969) Ms_FM = 1800 # Periode des Musters aus parallelen Streifen in Mikrometern # Für unsere parallelen Streifen mit Breite 5 Mikrometern wären das # 10 Mikrometer d_dom = 10.0 k = 2*np.pi/d_dom # Bei ungleichmäßiger Breite der Parallelstreifen-Domänen durch den # lithografischen Prozess: Aufteilung der Domänenbreiten bei gegebener # Periodenlänge d_dom der in -x und in +x magnetisierten Domänen # Eingabe der Domänenbreite der in -x magnetisierten Domänen in Mikrometern b_links = 5 if b_links < d_dom: b_rechts = d_dom - b_links else: print('WARNING: domain width too broad for given periodicity') # Durchnummerierung und Anzahl der Perioden des head-to-head/ tail-to-tail # Domänenmusters für das gerechnet werden soll Periode_d_dom = np.linspace(-3, 5, 9, endpoint = True) #------------------Partikeleigenschaften-------------- # Partikelradius in Mikrometer. R_SPP = 1.0 # Höhe der Partikelmitte über dem FM für die gerechnet werden soll # (größer als der Radius) in Mikrometern z0 = 1.6 # Gebe Warnung aus, wenn der Partikelradius größer als der Abstand des # Partikelzentrums zur Oberfläche ist if z0 <= R_SPP: print('WARNING: particle radius larger than particle-center distance to surface') # Effektive Magnetische Suszeptibiität der Partikel bei kleinen Feldern # nach der Erörterung von Markus Gusenbauer et al...Effektiv bedeutet hier, # dass die Differenz zwischen Partikelsuszeptibilität und der des umgebenden # Mediums (hier: Wasser) benutzt wird. Da dieser Wert mit einem # Experiment (Partikel in Wasser) verglichen wurde, beinhaltet er # automatisch auch den Effekt des Entmagnetisierungsfeldes für kugelförmige # Objekte (Entmagnetisierungsfaktor dort 1/3). Sollte nur die Suszeptibililität # eines unendlich ausgedehnten Materials gegeben sein, dann ist der Zusammenhang # zwischen dieser Suszeptibilität und der Suszeptibilität einer Kugel aus dem # gleichen Material gegeben durch: chi_Kugel = 3*chi_Mat/(3 + chi_Mat) chi_Kugel = 0.383 #--------------Eigenschaften des externen Feldes----------------------------- # Nach der Erörterung im eigenen Manuskript ist die Amplitude der z-Komponente # des Streufeldes # für ein SPP mit Radius 1 Mikrometer in einer Höhe von 1,6 Mikrometern # oberhalb der ferromagnetischen Schicht etwa 0.55 kA/m. Wir wählen daher # die Amplitude der x-Komponente des externen Feldes in kA/m in etwa in der # gleichen Größenordnung Hextmax_x = 1.5 # Amplitude der z-Komponente des externen Feldes in kA/m Hextmax_z = 0.5 # Frequenz der x-Komponente des externen Feldes in Hz nu_x = 12 # Frequenz der z-Komponente des externen Feldes in Hz nu_z = 12 # Phasendifferenz der Pulse zum Zeitpunkt t=0 dphi = np.pi/2 #-------Optimierung der Framerate in der Simulation--------- # Für wieviele Perioden des äußeren Feldes soll gerechnet werden? # Es soll hier die Anzahl der Perioden des langsamer gepulsten Feldes # angegeben werden, für die gerechnet werden soll N_Per_langsam = 6 # Ermittlung der Anzahl N_Per der Perioden für die gerechnet werden soll # und Ermittlung der Periodendauern T_langsam und T_schnell der langsamen und # schnellen Pulsung if nu_x >= nu_z: nu = nu_x N_Per_schnell = N_Per_langsam*nu_x/nu_z T_langsam = 1/nu_z T_schnell = 1/nu_x else: N_Per_schnell = N_Per_langsam*nu_z/nu_x nu = nu_z T_langsam = 1/nu_x T_schnell = 1/nu_z # Zahl der Stroboskopbilder je Sekunde (Framerate): # Hier wird eine dynamische Frame_Rate benutzt: Es sollen je Periode der # schnelleren Pulsung 72 Frames (ohne letzten Punkt) generiert werden # (alle 5° ein Bild, sinnvoll für Frequenzen im niedrigen zweistelligen # Bereich). Damit ergibt sich die Frame_Rate für die Rechnung zu # 72/T_schnell oder 72*nu Frame_Rate = 72*nu # Gesamtzahl der Frames für die Rechnung: # Die Gesamtzeit für die zu berechnenden Bilder ergibt sich zu # N_Per_langsam*T_langsam. Die Gesamtzahl der Bilder ergibt sich zu # N_Per_langsam*T_langsam*Frame_Rate Frames_Zahl = int(N_Per_langsam*T_langsam*Frame_Rate) print(T_langsam) print(T_schnell) print(N_Per_langsam) print(N_Per_schnell) print(Frame_Rate) print('Echtzeit für Rechnung', N_Per_langsam*T_langsam, 's') print('Framesanzahl', Frames_Zahl) # Liste für die Durchnummerierung der Stroboskopbilder für das gesamte Video: # Die Liste hat Frames_Zahl Elemente, die Indizierung startet bei 0 und hört # bei (Frames_Zahl-1) auf framenoArray = np.linspace(0, Frames_Zahl, endpoint = False, num=int(Frames_Zahl), dtype=int) print(framenoArray[0],framenoArray[1],framenoArray[2],framenoArray[-1]) # Initialisierung des Bildzählers, zeigt auf erstes Bild mit Index 0 framno_ctr = 0 # Array der echten Zeitpunkte, für die Frames erstellt werden sollen, # also für die die x- und z-Komponenten der MFL und die # dynamische Potenziallandschaft berechnet werden sollen tArray = 1/Frame_Rate*framenoArray print(tArray[0],tArray[1],tArray[-1]) # ----Generierung der Berechnungsdaten und Abspeichern im h5-file-------- # Wichtungsfaktor für die laterale Streufeldmittelung über das Partikelvolumen # wird für diese Rechnung numerisch aus dem Ergebnis für die gesonderte Rechnung # genommen. Für den hier gerechneten Fall mit R=1mum und d=10mum ist der # als konstant angenommene Reduktionsfaktor 15% (also 0.85). # Nun wird der Parameter A berechnet, der (vermutlich: siehe eigene # Ausarbeitung) die Mittelung der Streufelder über das Partikelvolumen # grob nähert. Für diese Berechnung wird der A-Faktor nach Yellen et al # Lab Chip 7, 1681 (2007) berechnet, modifiziert mit der in Zingsem # et al. JPD 50, 495006 (2017) angegebenen magnetischen Poldichte, sowie der # dort ebenfalls angegebenen Näherung für die Überlagerung der # Streufelder aus benachbarten Domänenwänden. Es wird # dazu das Magnetfeld als Funktion von z über einer head-to-head Domänenwand # mit gleichen Sättigungsmagnetisierungen rechts und links der Wand # angenommen. Dies kann später noch geändert werden, da die # Sättigungsmagnetisierungen links und rechts der Wand wegen des # Ionenbeschusses wahrscheinlich unterschiedlich sind. # Die Domänenwand wird als infinitesimal schmal angenommen. A_part = 0.85*(Ms_FM*t_FM)/(2*np.pi)*(1/z0-2/(np.sqrt(z0**2+(d_dom/2)**2))+2/np.sqrt(z0**2+4*(d_dom/2)**2)-4/d_dom*(np.log(2)-0.5)) # Definition der Funktionen für die beiden Feldpulskomponenten zum Zeitpunkt # t und als Funktion des Ortes x, hier im trigonometrischen Modell def Hzeff(x,t): return A_part*np.cos(2*np.pi/d_dom*x)+Hextmax_z*np.sin(2*np.pi*nu_z*t+dphi) def Hxeff(x,t): return A_part*np.sin(2*np.pi/d_dom*x)+Hextmax_x*np.sin(2*np.pi*nu_x*t) def Ueff(x,t): return -chi_Kugel*4/3*np.pi*R_SPP**3*mu_0*(Hzeff(x,t)**2+Hxeff(x,t)**2)/1e-6 # Wertearray für die Stützstellen des relevanten x-Bereichs für die # Berechnung der Felder. Stützstellen ergeben sich als 100/ Mikrometer # für den lateralen Bereich der durch Per_d*d_dom beschrieben wird. # Zum späteren Durchzählen wird das x_Bereichsarray als Integer-Zahlen # definiert x_Bereich = np.linspace(Periode_d_dom[0]*d_dom*100, Periode_d_dom[-1]*d_dom*100, int((Periode_d_dom[-1]*d_dom-Periode_d_dom[0]*d_dom))*100+1, dtype=int, endpoint = True) print(x_Bereich[0:6],x_Bereich[2998:3001],x_Bereich[7998:8001]) # Initialisierung des Resultatearrays: Es besteht aus der Anzahl an # Stützstellen an Zeilen und den Spalten x-Position, Hxeff, # Hzeff, Ueff und dUeff/dx. Das ResultateArray # wird für jedes Frame überschrieben und zuvor für das jeweilige Frame # im h5-file abgespeichert. ResultateArray = np.zeros((x_Bereich[-1]-x_Bereich[0]+1, 5)) print(np.shape(ResultateArray)) # Generiere ein h5 File für den vorliegenden Parametersatz. # Für jeden Frame wird ein Datensatz abgespeichert (also die Zahl der # Datensätze entspricht der Zahl der Frames und diese werden nach der # Nummer des jeweiligen Stroboskopbilds (Frames) benannt. Jeder # Datensatz besteht aus einem 2d-array mit einer Anzahl von Zeilen, # die der Anzahl von Stützstellen entspricht und 4 Spalten, welche # die Position der Stützstelle, die Hx-Komponente, die Hz-Komponente # und das Potential U enthalten. # Öffne zunächst die Datei mit dem Attribut append # und lege sie in der Progammvariable h5f ab. with h5py.File('./TWMP_PD_Hxmax='+str(round(Hextmax_x,1))+'_Hzmax='+str(round(Hextmax_z,1))+'_fx='+str(round(nu_x,2))+'_fz='+str(round(nu_z,2))+'_sinsin.h5', 'a') as h5f: for frame in framenoArray: t0 = 1/Frame_Rate*frame #Berechne den echten Zeitpunkt des jeweiligen Frames print('FrameNo:', frame) print('Echtzeit:', t0) # Für jedes Frame berechne die Werte für alle Stützstellen. x0 läuft von -3000 # bis +5000 (zB), aber die Indizierung in den Arrays startet immer bei 0, # daher muss der Zähler xctr (läuft zB. von 0 bis 8000) vom # x-Wert x0 unterschieden werden. Die Funktionen werden an den Orten mit # der Einheit Mikrometern berechnet, sodass im Argument der # Funktionen der Stützstellenwert x0 immer durch # 100 dividiert werden muss. for x0 in x_Bereich: xctr = int(x0)-int(x_Bereich[0]) ResultateArray[xctr,0] = x0/100 ResultateArray[xctr,1] = Hxeff(x0/100,t0) ResultateArray[xctr,2] = Hzeff(x0/100,t0) ResultateArray[xctr,3] = Ueff(x0/100,t0) ResultateArray[:,4] = -np.gradient(Ueff(x_Bereich/100,t0),ResultateArray[:,0]) print('x0 = ', ResultateArray[0:6,0],ResultateArray[2998:3001,0],ResultateArray[7998:8001,0]) print('t0 = ', t0) print('Hx = ', ResultateArray[2998,1],ResultateArray[2999,1],ResultateArray[3000,1],ResultateArray[3001,1]) print('Hz = ', ResultateArray[100,2]) print('Ueff = ', ResultateArray[100,3]) print('dUeff/dx = ', ResultateArray[100,4]) print(ResultateArray.shape) h5f.create_dataset('FrameNo'+str(frame), data=ResultateArray) h5f.close() # ----------------------------------------------------- # Ende Dateneingabe und Datengenerierung # ----------------------------------------------------- # ----------------------------------------------------- # Auswertung der Simulationen und Darstellung der Ergebnisse # (1): Laterale Geschwindigkeit eines Minimums der dPEL # ----------------------------------------------------- # Lade die Resultate Arrays für jedes frame with h5py.File('./TWMP_PD_Hxmax='+str(round(Hextmax_x,1))+'_Hzmax='+str(round(Hextmax_z,1))+'_fx='+str(round(nu_x,2))+'_fz='+str(round(nu_z,2))+'_sinsin.h5', 'r') as h5f: # for frame in framenoArray: # loaded_ResultateArray = h5f['FrameNo'+str(frame)][:] loaded_ResultateArray = h5f['FrameNo'+'0'][:] # Bestimme für das jeweilige frame den Wert des tiefsten Minimums für # das Potential und den höchsten Wert für den Gradienten des # Potentials MinimumPotValue = np.min(loaded_ResultateArray[:,3]) MaximumPotGradValue = np.max(loaded_ResultateArray[:,4]) # Bestimme für das jeweilige frame den Index für das Minimum des Potentials # mit dem kleinsten x-Wert und den Index für das Maximum des Gradienten; # ist notwendig, da mehrere Minima/ Maxima mit gleichem Wert zu erwarten sind MinimumPotValueIndex = np.min(np.where(loaded_ResultateArray[:,3] == MinimumPotValue)) MaximumPotGradValueIndex = np.min(np.where(loaded_ResultateArray[:,4] == MaximumPotGradValue)) # Bestimme die x-Position dieses Potentialminimums in Mikrometern und des # Gradientenmaximums ebenfalls in Mikrometern MinimumPotValuePosition = (-3000+MinimumPotValueIndex)/100 MaximumPotGradValuePosition = (-3000+MaximumPotGradValueIndex)/100 # Erstelle ein Array mit den Positionen des print('Potentialminimum = ', MinimumPotValue) print('Potentialgradientenmaximum = ', MaximumPotGradValue) print('Minimum Potential Value Index = ', MinimumPotValueIndex) print('Maximum Potential Gradient Value Index = ', MaximumPotGradValueIndex) print('Minimum Potential Value Position = ', MinimumPotValuePosition) print('Maximum Potential Gradient Value Position = ', MaximumPotGradValuePosition) plt.plot(loaded_ResultateArray[:,0],loaded_ResultateArray[:,3], label='Ueff') plt.plot(loaded_ResultateArray[:,0],loaded_ResultateArray[:,4], label='dUeff/dx') plt.legend(loc='upper left') #plt.ylim(0, 4.5) plt.show() h5f.close() # ---------------------------------------------------- # Diagrammeigenschaften # ---------------------------------------------------- with h5py.File('./TWMP_PD_Hxmax='+str(round(Hextmax_x,1))+'_Hzmax='+str(round(Hextmax_z,1))+'_fx='+str(round(nu_x,2))+'_fz='+str(round(nu_z,2))+'_sinsin.h5', 'r') as h5f: for frame in framenoArray: t0 = 1/Frame_Rate*frame #Berechne die echte Zeit des jeweiligen Frames print('FrameNo:', frame) print('Echtzeit:', t0) # lade für jedes Frame das ResultateArray in loaded_ResultateArray loaded_ResultateArray = h5f['FrameNo'+str(frame)][:] # check ob richtig geladen print('x0 = ', loaded_ResultateArray[0:6,0],loaded_ResultateArray[2998:3001,0],loaded_ResultateArray[7998:8001,0]) print('t0 = ', t0) print('Hx = ', loaded_ResultateArray[100,1]) print('Hz = ', loaded_ResultateArray[100,2]) print('Ueff = ', loaded_ResultateArray[100,3]) print(loaded_ResultateArray.shape) # Definition der Diagrammgröße in Zoll und Auflösung in dpi # Breite w = 8 # Höhe h = 6 # Auflösung aufl = 200 # Erstelle ein neues Matplot-Figure-Objekt mit w x h Zoll und # einer Auflösung von aufl dpi: fig = plt.figure(figsize=(w, h), dpi=aufl) # Titel der Abbildung: # ax.set_title('Sinus und Cosinus', fontsize=20, color='black') # Festlegung der Anzahl der Panels in der Abbildung # In dieser Abbildung ein 1x1 großes Panel-Gitter erstellen; # Als aktuelles Panel wird das erste dieses Gitters ausgewählt: sub1 = fig.add_subplot(211) sub2 = fig.add_subplot(212) fig.subplots_adjust(hspace=0) # ---------------------------------------------------- # Ende Diagrammeigenschaften # ---------------------------------------------------- # ----------------------------------------------------- # Achsen, Gitter und Tickparameter # ----------------------------------------------------- # y-Achsenbeschriftung: sub1.set_ylabel(r'H$_{eff,x,z}$ [kA/m]', fontsize=18) sub2.set_ylabel(r'U$_{eff,x,z}$ [$10^{-17}$ J]', fontsize=18) # secax = ax.secondary_yaxis('right') # secax.set_ylabel('angular rotation frequency [1/s]', fontsize=16) # x-Achsenbeschriftung: sub2.set_xlabel('x [µm]', fontsize=18) # Zeige Minor Ticks? # ax.minorticks_on() # Richtung und Größe der Ticks # ax.tick_params(which='minor', length=5, width=2, direction='in') sub1.tick_params(which='major', length=10, width=2, direction='in', bottom=True, top=True, left=True, right=True, color='black') sub2.tick_params(which='major', length=10, width=2, direction='in', bottom=True, top=True, left=True, right=True, color='black') # Schalte x-Achsenbeschriftung für das obere Panel aus sub1.set_xticklabels([]) # Setze die y-Achsenbeschriftungsgröße auf 12 pt sub1.tick_params(axis='both', which='major', labelsize=16) sub2.tick_params(axis='both', which='major', labelsize=16) # Soll ein Gitternetz mit eingezeichnet werden? sub1.grid(which='both', ls='--') sub2.grid(which='both', ls='--') # Achsen-Wertebereiche festlegen; wenn nicht festgelegt wird der Bereich # automatisch generiert # Syntax: plt.axis([xmin, xmax, ymin, ymax]) sub1.set_xlim(-30,50) sub2.set_xlim(-30,50) sub1.set_ylim(-4.99,4.99) sub2.set_ylim(-3.29,-0.0001) # Drucke die Achsen in dickerer Strichstärke mit # mpl.spines.Spine(ax, 'top', linewidth=2) sub1.spines['top'].set_linewidth(2) sub2.spines['top'].set_linewidth(2) sub1.spines['bottom'].set_linewidth(2) sub2.spines['bottom'].set_linewidth(2) sub1.spines['left'].set_linewidth(2) sub2.spines['left'].set_linewidth(2) sub1.spines['right'].set_linewidth(2) sub2.spines['right'].set_linewidth(2) # ----------------------------------------------------- # Ende Achsen, Gitter und Tickparameter # ----------------------------------------------------- # ------------------------------------------------------ # Darstellung der Funktion und speichern # ------------------------------------------------------ # Zeige die Echtzeit des Frames im Bild an sub1.text(-28, 4.99, 't ='+'{0:.3f}'.format(t0)+' s', backgroundcolor='w', fontsize=14) # arrowprops=dict(arrowstyle="->", relpos=(0., 0.)), # horizontalalignment='right', verticalalignment='top', # ) # Markiere mit einem Pfeil eine Stelle eines Anfangsmaximums des # Potentials um die Dynamik in Bildern zu veranschaulichen sub2.annotate('', xytext=(-20, 0), xy=(-20, -0.5), textcoords='data', arrowprops=dict(arrowstyle="simple", relpos=(0., 0.), color='b'), horizontalalignment='right', verticalalignment='top', ) #ax.annotate('r = 2.0 $\\mu$m', xy=(1.7, 1.5), xycoords='data', # xytext=(6, 2.5), textcoords='data', fontsize=16, # arrowprops=dict(arrowstyle="->", relpos=(0., 0.)), # horizontalalignment='right', verticalalignment='top', # ) #ax.annotate('r = 4.0 $\\mu$m', xy=(2.4, 1.65), xycoords='data', # xytext=(6, 2.25), textcoords='data', fontsize=16, # arrowprops=dict(arrowstyle="->", relpos=(0., 0.)), # horizontalalignment='right', verticalalignment='top', # ) #ax.annotate('r = 8.0 $\\mu$m', xy=(3.3, 1.8), xycoords='data', # xytext=(6, 2.0), textcoords='data', fontsize=16, # arrowprops=dict(arrowstyle="->", relpos=(0., 0.)), # horizontalalignment='right', verticalalignment='top', # ) # Funktion in blauer Farbe mit durchgezogener Linie und 2 Pixel Breite # plotten: plt.plot(x, cos_x, color="blue", linewidth=1.0, linestyle="-") # Plotte die x-Feldkomponente in schwarz gestrichelt, die z-Feldkomponente # in blau durchgezogen und das effektive POtential in blau durchgezogen # in separatem Panel unter den Feldkomponenten mit Strichstärke 2: sub1.plot(loaded_ResultateArray[:,0], loaded_ResultateArray[:,2], color="blue", linewidth=2.0, linestyle="-") sub1.plot(loaded_ResultateArray[:,0], loaded_ResultateArray[:,1], color="black", linewidth=2.0, linestyle="--") sub2.plot(loaded_ResultateArray[:,0], loaded_ResultateArray[:,3]/10, color="red", linewidth=2.0, linestyle="-.") plt.yticks(fontsize=16) plt.xticks(fontsize=16) # Speichern als PNG-Datei: fig.savefig('TWMP_PD_Hxmax='+str(round(Hextmax_x,1))+'_Hzmax='+str(round(Hextmax_z,1))+'_fx='+str(round(nu_x,2))+'_fz='+str(round(nu_z,2))+'_sinsin_t='+'{0:.3f}'.format(t0)+'.png', dpi=200) plt.show() h5f.close() # ------------------------------------------------------------ # Erstelle das Video # ------------------------------------------------------------ # Eingabe- und Ausgabepfade für Bilder und Video definieren #image_folder = './images/' # Ersetzen Sie durch den Pfad zu Ihren Bildern image_folder = './' # Name der Ausgabevideodatei video_name = 'TMWP_PD_Hxmax='+str(round(Hextmax_x,1))+'_Hzmax='+str(round(Hextmax_z,1))+'_fx='+str(round(nu_x,2))+'_fz='+str(round(nu_z,2))+'_sinsin.mp4' # Liste der Bilder abrufen images = [img for img in os.listdir(image_folder) if img.endswith(".png")] # Lesen Sie das erste Bild, um die Größe zu erhalten first_image = cv2.imread(os.path.join(image_folder, images[0])) height, width, layers = first_image.shape # Video-Writer initialisieren # 'mp4v' für das MP4-Format verwenden fourcc = cv2.VideoWriter_fourcc(*'mp4v') # Video erstellen mit Framerate 10 und Bilder bestimmter Breite und Höhe Video = cv2.VideoWriter(video_name, fourcc, 10, (width, height)) # Bild für Bild zum Video hinzufügen for image in images: img_path = os.path.join(image_folder, image) frame = cv2.imread(img_path) Video.write(frame) # Video-Writer freigeben Video.release() cv2.destroyAllWindows()