from math import* from pylab import* import numpy as np import datetime from scipy.stats import norm import sys from scipy import linalg import time import timeit T0=time.time() ########################################################################## '''Dziala robienie widma dla dowolnego hamiltonianu jesli sie poda macierze przejscia Dodano mozliwosc roznienia mapy w funkcji kata w plaszyczyznie studni oraz obracania kierunkiem E''' ########################################################################## ########################### STALE LICZBOWE ############################## hbar=1. # [eVs] = 1 bo juz w magnetonie bohra siedzi hbar ha=4.135667516e-10 # ueV*s - Planck constant miu = 5.7883818066e-5 # [eV/T] magneton Bohra AA=(680E-9) # A*IS [eV] i tu tez jest hbar gJon=2.0 # bezwymiarowe aa=(320E-9)*1 # [eV] no i tu tez jest hbar kb = 8.617333262E-5 # Boltzman constant eV/K #mub*g*dB=dE=0.5788 meV/T*2.002*dB=1.1587meV/T*dB def Energia_na_GHz(Energia): '''Energia w ueV na GHz E/h=f''' return Energia/ha/10e9 def GHz_na_energie(GHz): '''Energia w Gz na ueV E=f*h''' return GHz*10e9*ha #ha jest w ueV*s #ustawienie szerokosc wyswietlania np.core.arrayprint._line_width=300 ########################################################################## ########################################################################## SpinJadrowyJonu=5./2. # trzeba wpisac z kropka aby dzielil float SpinElektronowyJonu=5./2. # trzeba wpisac z kropka aby dzielil float JonI_wymiar=int(2*SpinJadrowyJonu+1) JonS_wymiar=int(2*SpinElektronowyJonu+1) #kron ########## WARTOSCI LICZBOWE PIERWIASTKOW PRZY OPERATORACH ############### def Pierwiastek_Operator_pl( Spin , Spin_z): return hbar*np.sqrt( Spin*(Spin + 1 ) - Spin_z*(Spin_z + 1) ) ##################################### def Pierwiastek_Operator_mi( Spin , Spin_z): return hbar*np.sqrt( Spin*(Spin + 1 ) - Spin_z*(Spin_z - 1) ) ############## SPIN PLUS & MINUS & Z ##################################### ########################################################################## def OperatorPlus(Spin): '''S+|Spin,Spin_z>=sqrt(S(S+1)-Sz(Sz+1))|Spin,Spin_z+1> w osi pionowej a, poziomo b''' Operator = np.zeros( (int(Spin*2+1),int(Spin*2+1)), complex ) for a in np.arange(-Spin,Spin+1,1): for b in np.arange(-Spin,Spin+1,1): if (a-b)==1 and b<=(Spin-1): index_b=int(b+Spin) index_a=int(a+Spin) Operator[index_b,index_a]+=Pierwiastek_Operator_pl( Spin , b) return Operator ########################################################################## def OperatorMinus(Spin): '''S+|Spin,Spin_z>=sqrt(S(S+1)-Sz(Sz-1))|Spin,Spin_z+1> w osi pionowej a, poziomo b''' Operator = np.zeros( (int(Spin*2+1),int(Spin*2+1)), complex ) for a in np.arange(-Spin,Spin+1,1): for b in np.arange(-Spin,Spin+1,1): if (a-b)==-1 and b>=-(Spin-1): index_b=b+Spin index_a=a+Spin Operator[int(index_b),int(index_a)]+=Pierwiastek_Operator_mi( Spin , b) return Operator ########################################################################## def OperatorZ(Spin): '''SZ|Spin,Spin_z>=Sz|Spin,Spin_z> w osi pionowej a, poziomo b''' Operator = np.zeros( (int(Spin*2+1),int(Spin*2+1)), complex ) for a in np.arange(-Spin,Spin+1,1): for b in np.arange(-Spin,Spin+1,1): if (a==b): index_b=b+Spin index_a=a+Spin Operator[int(index_b),int(index_a)]+=b return Operator ########################################################################## def Rotation(Matrix,alphaDeg=0,ProperJz=''): '''Obrot wokol osi Z : A'=R*A*R^-1=exp(-i*alpha*Jz)*A*exp(i*alpha*Jz) ProperJz to odpowiadajaca wymiarem macierzy Matrix macierz zetowa alphaDeg to kat o ktory obracamy w stopniach Sx ->Sy dla alpha=-90deg''' alpha=alphaDeg*np.pi/180. NewMatrix=np.dot(linalg.expm(-1.j*alpha*ProperJz),np.dot(Matrix,linalg.expm(1.j*alpha*ProperJz))) return NewMatrix ########################################################################## ########################################################################## IJonPl=OperatorPlus(SpinJadrowyJonu) IJonMi=OperatorMinus(SpinJadrowyJonu) SJonPl=OperatorPlus(SpinElektronowyJonu) SJonMi=OperatorMinus(SpinElektronowyJonu) SJonX = (SJonPl + SJonMi) / 2.0 SJonY = (SJonPl - SJonMi) / 2.0j SJonZ = OperatorZ(SpinElektronowyJonu) IJonX = (IJonPl + IJonMi) / 2.0 IJonY = (IJonPl - IJonMi) / 2.0j IJonZ = OperatorZ(SpinJadrowyJonu) OneSJon=np.eye(int(2*SpinElektronowyJonu+1)) OneIJon=np.eye(int(2*SpinJadrowyJonu+1)) OneIS =np.kron(OneIJon,OneSJon) #KOLEJNOSC W ILOCZYNACH TENSOROWYCH: # kron(elekron,kron(dziura,kron(SpinJadroJonu,spinElektrJonu))) # wtedy mamy # e0,dz0,I0,S0 # e0,dz0,I0,S1 # ...# e0,dz0,I1,S0 # e0,dz0,I1,S1# itd #zgodnie z zasadami mnozenia tensorowego (kronecker product): #dla A=[[0,1][-1,0]] #b = [[b00,b01][b10,b11]) #C = [[C00,C01][C10,C11]) # kron (A,b,C) # # | 0 0 | 1*b00*[C] 1*b01*[C] | # | 0 0 | 1*b10*[C] 1*b11*[C] | # | -1*b00*[C] -1*b01*[C] | 0 0 | # | -1*b00*[C] -1*b01*[C] | 0 0 | # #gdzie [C] to cala macierz , wiec ostatecznie wymiar to 2*2*2=8 # # wymiar hamiltonianu to ElS_wymiar*DzS_wymiar*JonI_wymiar*JonS_wymiar ########################################################################## ########################################################################## def HamiltonianJonuZeeman(B=0.6,ThetaDeg=0,PhiDeg=0): '''Zwraca hamiltoniann dla zadanego kata pola podanego w stopniach dla wartosci Pola B w eV Kat Theta liczony od kierunku Z do plaszczyzny Kat Phi liczony od osi X do przeciwnie do ruchu wskazowek zegara parzac od gory na probke''' ThetaRad=np.pi*ThetaDeg/180. PhiRad=np.pi*PhiDeg/180. Bx= B*np.sin(ThetaRad)*np.cos(PhiRad) #*cos(PhixRad) By= B*np.sin(ThetaRad)*np.sin(PhiRad) #*sin(PhixRad) Bz= B*np.cos(ThetaRad) return gJon*miu*(Bx*np.kron(OneIJon,SJonX)+By*np.kron(OneIJon,SJonY)+Bz*np.kron(OneIJon,SJonZ)) ########################################################################## ########################################################################## def HamiltonianJonuIS(DD=8.75e-7,EE=0,alphaDeg=0): ''' Hamiltonian nie zalezny od czasu jonu bez pola magnetycznego w eV alpha to kat obrocenia HEstrain''' HHyperfine = AA *( np.kron(IJonX,SJonX)+np.kron(IJonY,SJonY)+np.kron(IJonZ,SJonZ) ) HDstrain = DD *( np.kron(OneIJon,np.linalg.matrix_power(SJonZ,2))-OneIS*SpinElektronowyJonu*(SpinElektronowyJonu+1)/3. ) #HEstrain = EE *(np.kron(OneIJon,np.linalg.matrix_power(SJonX,2))-np.kron(OneIJon,np.linalg.matrix_power(SJonY,2))) HEstrainRot = EE *(np.kron(OneIJon,Rotation(np.linalg.matrix_power(SJonX,2),alphaDeg,SJonZ))-np.kron(OneIJon,Rotation(np.linalg.matrix_power(SJonY,2),alphaDeg,SJonZ))) Hcrystal =(aa/6.)*(np.kron(OneIJon,np.linalg.matrix_power(SJonX,4))+np.kron(OneIJon,np.linalg.matrix_power(SJonY,4))+np.kron(OneIJon,np.linalg.matrix_power(SJonZ,4))-OneIS*SpinElektronowyJonu*(SpinElektronowyJonu+1.)*(3.*SpinElektronowyJonu**2+3.*SpinElektronowyJonu-1)*(1./5.) ) H=HHyperfine*1.+HDstrain*1.+HEstrainRot*1.+Hcrystal*1. return H ########################################################################## ########################################################################## def PrzejsciaRzadkie(P,tol=1e-8): '''Bierze rzadka macierz przejscia i zwraca indeksy elementow niezerowych ; P=[Pplus,Pminus]''' P0pocz=[] P0kon=[] P0wart=[] P0=P[0] P1pocz=[] P1kon=[] P1wart=[] P1=P[1] print( "Oblicznie niezerowych przejsc, p<",tol," sa pomijane!") for NoPoczatkowyWekt in range(np.shape(P[0])[0]): for NoKoncowyWekt in range(np.shape(P[0])[1]): if abs(P0[NoPoczatkowyWekt,NoKoncowyWekt])>tol: P0pocz.append(NoPoczatkowyWekt) P0kon.append(NoKoncowyWekt) P0wart.append(P0[NoPoczatkowyWekt,NoKoncowyWekt]) if abs(P1[NoPoczatkowyWekt,NoKoncowyWekt])>tol: P1pocz.append(NoPoczatkowyWekt) P1kon.append(NoKoncowyWekt) P1wart.append(P1[NoPoczatkowyWekt,NoKoncowyWekt]) print( " P0 ma ", len(P0pocz)," niezerowych przejsc") print( " P1 ma ", len(P1pocz)," niezerowych przejsc") return P0pocz,P0kon,P0wart,P1pocz,P1kon,P1wart ########################################################################## def plot_energy_vs_B(B_min,B_max,B_step,ifSave=False): fig = plt.figure(figsize=(10,10)) #a = fig.subplots(2,2) #rows, cols grid = plt.GridSpec(2, 3, wspace=0.4, hspace=0.3) leftup = fig.add_subplot(grid[:,:2]) #leftdown = fig.add_subplot(grid[1,:2]) textplot = fig.add_subplot(grid[:,2],frameon = True) textplot.set_xticks([]) textplot.set_yticks([]) textplot.set_ylim(0, 100) textplot.set_xlim(0, 100) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") t = textplot.text(50, 107, "Info", ha="left", va="bottom", color="k", family="sans-serif", fontweight="bold", fontsize=16) t = textplot.text(60, 104, timestamp, ha="center", va="bottom", color="k", family="sans-serif", fontweight="normal", fontsize=12) leftup.set_title('Energy spectrum') leftup.set_ylabel('Energy (meV)') leftup.set_xlabel('Magnetic field (B)') Es = [] for b in np.arange(B_min, B_max, B_step): Hmag = HamiltonianJonuZeeman(b, ThetaDeg=0, PhiDeg=45) E,v = np.linalg.eigh(H1 + Hmag) Es.append(E) Es = np.array(Es) for i in range(len(E)): leftup.plot(np.arange(B_min, B_max, B_step),Es[:,]*1000,'-')#,color='black' #leftup.set_xlim(min(ys1D),max(ys1D)) plt.show() ### Uwzglednienie temperatury, czynnik Boltzmanna ####################### def Boltzman_factor_calc(EeV,T_K): p = np.exp(-EeV/(kb*T_K)) return p def Boltzman_factor_find(EeV,T_K,dic): if (EeV,T_K) in dic: return dic[(EeV,T_K)] else: f = Boltzman_factor_calc(EeV,T_K) dic[(EeV, T_K)] = f dic['calc'] += 1 return f ####################################################################### def MakeSpektrum(H,P=0,dEbin=0,ifPlotSpectrum=False,ifPlotTransitions=False,ifPlotSave=False,MinTransitionStrength=1e-4,temperature_K=1000): '''Zwraca spektrum w meV zbinowane o szerokosci binu dEbin w meV, P to macierz przejscia podanie dEbin=0 nie binuje tylko kazde przejscie zostaje zapisane linie NIE sa rozmazane gaussem!!! ifPlotSpectrum wyswietla pojedyncze widmo Prawd(E) ifPlotTransitions wyswietla przejscia i poziomy energetyczne ifPlotSave wlacza automatyczne zapisywanie do pdf obu plikow jako Transitions_%Y%m%d-%H%M%S.pdf Energy_spectrum_%Y%m%d-%H%M%S.pdf MinTransitionStrength minimalna sila przejcia dla ktorej jest dodawane przejcie do listy UWAGA!:aby utworzyc widmo ODMR dEbin musi byc rowny zero!!! ''' Energie=[] Prawdopodobienstwa=[] WartosciW,WektoryW= np.linalg.eigh(H) #WartosciW,WektoryW = zip(*sorted(zip(WartosciW,WektoryW))) #WartosciW,WektoryW = np.array(WartosciW),np.array(WektoryW) #print WartosciW,WektoryW if np.all(np.shape(np.array(P))[0]==0): raise ValueError('P not detected!!!') # if ifPlotTransitions==True: figure() for NoPoczatkowyWekt in range(np.shape(H)[0]): for NoKoncowyWekt in range(NoPoczatkowyWekt*1,np.shape(H)[0]): dEmeV=(WartosciW[NoKoncowyWekt]-WartosciW[NoPoczatkowyWekt])/1e-3 if dEmeV>0 and NoPoczatkowyWekt!=NoKoncowyWekt: #dEmeV=(WartosciW[NoKoncowyWekt]-WartosciW[NoPoczatkowyWekt])/1e-3 #7.36s exp #5.81s dic #5.27s fixed PwekPocz=np.dot(P[0],WektoryW[:,NoPoczatkowyWekt]) SilaOsc=np.vdot(WektoryW[:,NoKoncowyWekt],PwekPocz) #Boltz_fac = 1 Boltz_fac_1 = Boltzman_factor_find(WartosciW[NoPoczatkowyWekt],temperature_K,Boltzmann_dic) Boltz_fac_2 = Boltzman_factor_find(WartosciW[NoKoncowyWekt],temperature_K,Boltzmann_dic) Boltz_fac = Boltz_fac_1 - Boltz_fac_2 #Boltz_fac = Boltzman_factor_find(WartosciW[NoPoczatkowyWekt],temperature_K,Boltzmann_dic) #Boltz_fac_sum += Boltz_fac #Boltz_fac = np.exp(-WartosciW[NoPoczatkowyWekt]/(kb*temperature_K)) PrawdJedn=1*abs(SilaOsc)**2*Boltz_fac #tu jest wszytko #PwekPocz=np.dot(P[1],WektoryW[:,NoPoczatkowyWekt]) #SilaOsc=np.vdot(WektoryW[:,NoKoncowyWekt],PwekPocz) #PrawdJedn+=1*abs(SilaOsc)**2 *np.exp(-WartosciW[NoPoczatkowyWekt]/(kb*temperature_K))#tu nie dziala #DEBUGOWANIE #oklist=range(1,36)+range(115,118) #####################33 #if NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt in oklist: if 0: SzIn,IzIn,SzKo,IzKo=0,0,0,0 for npo in range(len(WektoryW[NoPoczatkowyWekt])): if abs(WektoryW[NoPoczatkowyWekt][npo])>=0.1: SzIn=np.real(np.vdot(WektoryW[:,npo],np.dot(np.kron(OneIJon,SJonZ),WektoryW[:,npo])))#float(imag(T[npo][npo])) #S IzIn=np.real(np.vdot(WektoryW[:,npo],np.dot(np.kron(IJonZ,OneSJon),WektoryW[:,npo])))#float(real(T[npo][npo])) #I print( WektoryW[NoPoczatkowyWekt][npo],'IzIn:',IzIn,", SzIn: ",SzIn) for nko in range(len(WektoryW[NoKoncowyWekt])): if abs(WektoryW[NoKoncowyWekt][nko])>=0.1: SzKo=np.real(np.vdot(WektoryW[:,nko],np.dot(np.kron(OneIJon,SJonZ),WektoryW[:,nko])))#float(imag(T[nko][nko])) #S IzKo=np.real(np.vdot(WektoryW[:,nko],np.dot(np.kron(IJonZ,OneSJon),WektoryW[:,nko])))#float(real(T[nko][nko])) #I print( WektoryW[NoKoncowyWekt][nko],', Izko:',IzKo,", SzKo: ",SzKo) if True or (IzKo-IzIn==0 and abs(SzKo-SzIn)==1): print( "Para NO", NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt," $$$$$$$$$$$$ f=", PrawdJedn, " dE=",dEmeV) print( "dI=",IzKo-IzIn," , dS=",SzKo-SzIn, (" ........................................ PRZEJSCIE ..." if (IzKo-IzIn<=0.1 and abs(SzKo-SzIn)>=0.9) else " .................")) print( '_______________________________________________') #KONIEC DEBUGOWANIA if PrawdJedn>=MinTransitionStrength: Energie.append(dEmeV) Prawdopodobienstwa.append(PrawdJedn) #print(PrawdJedn,dEmeV,WartosciW[NoPoczatkowyWekt]) if ifPlotTransitions==True: plot(NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt,WartosciW[NoPoczatkowyWekt]/1e-3,lw=1,c='blue',marker='o') plot(NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt,WartosciW[NoKoncowyWekt]/1e-3,lw=1,c='black',marker='.') #if PrawdJedn>=MinTransitionStrength: plot([NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt,NoPoczatkowyWekt*np.shape(H)[0]+NoKoncowyWekt], [WartosciW[NoPoczatkowyWekt]/1e-3,WartosciW[NoKoncowyWekt]/1e-3], 'r-', alpha=erf(PrawdJedn),linewidth=2) #annotate('Text', xy=(... strzalki... if 1: Boltz_fac_sum = np.sum([abs(Boltzman_factor_find(Ei,temperature_K,Boltzmann_dic)) for Ei in WartosciW]) xxx = np.array(PrawdJedn) PrawdJedn = xxx/Boltz_fac_sum PrawdJedn = PrawdJedn.tolist() #print('Boltzmann!!!!!!!!',Boltz_fac_sum) if ifPlotTransitions==True: title('Energy transitions (blue = initial state)') xlabel('Transition number (initNo*total+finalNo)') ylabel('Energy [meV]') timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") if ifPlotSave==True: savefig('Transitions_'+timestamp+'.pdf') show() close() show() if dEbin==0: EnergieListaKoncowa=Energie PrawdopodobienstwaListaKoncowa=Prawdopodobienstwa PlotType='.' else: Emin=min(Energie)-5*dEbin Emax=max(Energie)+5*dEbin ileE=int((Emax-Emin)/dEbin) if ileE<1.1e5: print( "Spektrum sklada sie z ", ileE, " pikseli") else: raise ValueError("Za duzo pikseli!!!! zwieksz wartosc dEbin znacznie powyzej ", (Emax-Emin)/1e5, " meV") dEbin=(Emax-Emin)/1e6 Emin=min(Energie)-5*dEbin Emax=max(Energie)+5*dEbin NewProbs=np.zeros(int((Emax-Emin)/dEbin)) #to bedzie widmo NewEnergies=np.arange(Emin,Emax-dEbin,dEbin) wsp=ileE/(Emax-Emin) for i in range(len(Energie)): E=Energie[i] kanalE=int(wsp*E-wsp*Emin+0.5) NewProbs[kanalE]+=Prawdopodobienstwa[i] EnergieListaKoncowa=NewEnergies PrawdopodobienstwaListaKoncowa=NewProbs PlotType='-' if ifPlotSpectrum==True: plot(EnergieListaKoncowa,PrawdopodobienstwaListaKoncowa,PlotType) xlabel('Energy [meV]') ylabel('Probability [arb. units]') title('Energy spectrum') if ifPlotSave==True: timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") savefig('Energy_spectrum_'+timestamp+'.pdf') show() close() show() # #linalg.eig(a) Compute the eigenvalues and right eigenvectors of a square array. #linalg.eigh(a[, UPLO]) Return the eigenvalues and eigenvectors of a Hermitian or symmetric matrix. #linalg.eigvals(a) Compute the eigenvalues of a general matrix. #linalg.eigvalsh(a[, UPLO]) Compute the eigenvalues of a Hermitian or real symmetric matrix. # A.getH() #Wart, Wekt = np.linalg.eigh(H) return EnergieListaKoncowa,PrawdopodobienstwaListaKoncowa ########################################################################## def Gauss(ilePixGauss=111): gauss=[] #(1/(sqrt(2*pi)*sigma)*exp(-0.5*((x-x0)/sigma)**2) x0=(ilePixGauss-1)/2. sigma=ilePixGauss/550. for i in range(0,ilePixGauss,1): gauss.append(np.exp(-0.5*((i-x0)/sigma)**2)) return gauss ilePixGauss=50001 gauss=Gauss(ilePixGauss) def WartoscPrzejscDlaMikrofali(EnergieListaKoncowa,PrawdopodobienstwaListaKoncowa,MWfGHz=13.1,dfHz=0.25): # Emin...dErozmaz....E_MWfGHz....dErozmaz...Emax # dErozmaz=szerokosc gaussa tzn dfGHz*N #szerokosc realnie powinna miec okolo 0.01meV tzn 2GHz # poniewaz szerokosc gaussa jest bardzo mala #Tstart=time.time() N=1e9 dfHz=dfHz*1 dfGHz=1e-9*dfHz dErozmaz=GHz_na_energie(dfGHz*N)/1000. # meV=ueV*1000 E_MWfGHz=GHz_na_energie(MWfGHz)/1000. # meV=ueV*1000 #print " Ecentr: ",E_MWfGHz," dErozmaz: ",dErozmaz WartoscPrawdDlaMWfGHz=0 OffsetPix=int((ilePixGauss-1)/2.) dEbinGaussa=3.*dErozmaz/OffsetPix #n=0 for i in range(len(EnergieListaKoncowa)): Ei=EnergieListaKoncowa[i] #print "dErozmaz: ",dErozmaz,"E_MWfGHz: ",E_MWfGHz," P: ",WartoscPrawdDlaMWfGHz, ",Ei: ",Ei,abs(Ei-E_MWfGHz)<=dErozmaz dEi=abs(Ei-E_MWfGHz) if (dEi<=dErozmaz): #s=sqrt(GHz_na_energie(dfGHz)) #WartoscPrawdDlaMWfGHz+=PrawdopodobienstwaListaKoncowa[i]*norm.pdf(E_MWfGHz,scale=s,loc=Ei)/norm.pdf(Ei,scale=s,loc=Ei) WartoscPrawdDlaMWfGHz+=PrawdopodobienstwaListaKoncowa[i]*gauss[OffsetPix+int(dEi/dEbinGaussa)] #n+=1 #Tend=time.time() #print "WartoscPrzejscDlaMikrofali: dfHz=" , str(dfHz)," #E=",len(EnergieListaKoncowa)," ; time=", (Tend-Tstart)/10 #print EnergieListaKoncowa #print PrawdopodobienstwaListaKoncowa #print WartoscPrawdDlaMWfGHz return WartoscPrawdDlaMWfGHz ########################################################################## def nic(): pass def MakeSingleODMRspectrum(Bmin=0.3,Bmax=0.85,Bstep=0.0002,BthetaDeg=0,BphiDeg=0,MWfGHz=15.6,dfHz=0.7,ifPlot=False,ifSave=True,fname_additon='',temperature=1000): '''Zwraca liste pol magnetycznych i prawdopodobienstw dla widma ODMR dla zadanej czestosci mikrofal MWfGHz=15.6 w GHz, szerokosci linii gauss dfHz=0.7, kat pola wzgledem osi zet podany w stopniach BthetaDeg=0''' #T0=time.time() bx=[] px=[] Tham=0 Tspec=0 Tgauss=0 n=0 Bproc=(Bmax-Bmin)/Bstep/100. print( "ODMR spectr.for theta="+str(BthetaDeg),'deg',' phi =',str(BphiDeg),'deg') for b in np.arange(Bmin,Bmax,Bstep): T1=time.time() H=HamiltonianJonuZeeman(b,BthetaDeg,BphiDeg) T2=time.time() EL,PL=MakeSpektrum(H1+H,P=(P1,P2),ifPlotSpectrum=False,ifPlotTransitions=False,ifPlotSave=False,dEbin=0,temperature_K=temperature) T3=time.time() p=WartoscPrzejscDlaMikrofali(EL,PL,MWfGHz,dfHz) T4=time.time() Tham+=T2-T1 Tspec+=T3-T2 Tgauss+=T4-T3 n+=1. if False: print( "\rODMR spectr.for theta="+str(BthetaDeg)+" deg in progress: ",int(round((b-Bmin)/(Bmax-Bmin-Bstep)*100)),'%') bx.append(b) px.append(p) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") if ifSave else nic() if ifPlot: plot(bx,px,'.-',lw=0.5,ms=2) xlabel('Magnetic field B [T]') ylabel('Probability [arbitrary units]') title('ODMR spectrum for f='+str(MWfGHz)+" GHz") savefig('ODMR_Spectrum_'+str(MWfGHz)+"_GHz_"+timestamp+'.pdf') if ifSave else nic() show() if ifSave: SaveTwoColumn(bx,px,'ODMR_Spectrum_'+str(MWfGHz)+"_GHz_"+timestamp+fname_additon+'.txt') print( "\n-------------------------") print( "In MakeSingleODMRspectrum: HamiltonianJonuZeeman",Tham/n,"s; MakeSpektrum:",Tspec/n,"s ; WartoscPrzejscDlaMikrofali:",Tgauss/n,"s") return bx,px #----------------------------------------- def SaveTwoColumn(X,Y,fname): fout=open(fname,'w') for i in range(len(X)): fout.write(str(X[i])+'\t'+str(Y[i])+'\n') fout.close() ########################################################################## def PlotEnergyLevelsVsB(Bmin=0,Bmax=0.1,Bstep=0.00001,BthetaDeg=0,BphiDeg=0,ifPlot=True,ifSave=True): wartWiele=[] bxWiele=[] for b in np.arange(Bmin,Bmax,Bstep): H=HamiltonianJonuZeeman(b,BthetaDeg,BphiDeg) Wart, Wekt = np.linalg.eigh(H+1*H1) #np.vdot(Wekt[4],np.dot(np.kron(IJonZ,OneSJon),Wekt[4])) #np.vdot(Wekt[4],np.dot(np.kron(OneIJon,SJonZ),Wekt[4])) bx=b*np.ones(len(Wart)) wartWiele.append(np.sort(Wart*1e6)) bxWiele.append(bx) wartWiele=np.array(wartWiele) bxWiele=np.array(bxWiele) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") if ifSave else nic() def nic(): pass if ifPlot: for i in range(len(bxWiele[1,:])): xlabel('Magnetic field B [T]') ylabel('Energy [ueV]') title('Energy levels vs B') plot(bxWiele[:,i],wartWiele[:,i],'k-',lw=1,alpha=1) savefig('Energy_levels_ueV_vs_B_'+timestamp+'.pdf') if ifSave else nic() show() if ifSave: SaveMultiColumn(bxWiele,wartWiele,'Energy_levels_ueV_vs_B_'+timestamp+'.txt') return bxWiele,wartWiele #----------------------------------------- def SaveMultiColumn(Xs,Ys,fname): ''' niech x=[1,2,3,4], a y sa dwie kolumny y1 y2 format dotej funkcji musi wynoscic: xs=[[1,1],[2,2],[3,3],[4,4]] ys=[[y11,y21],[y12,y22],[y13,y23],[y14,y24]] dedykowane funkcji PlotEnergyLevelsVsB''' fout=open(fname,'w') for wiersz in range(len(Xs[:,0])): fout.write(str(Xs[wiersz,0])+'\t') #we wszytkich jest ten sam x dla tego mozna braz z [0] for kolumna in range(len(Ys[0])): fout.write(str(Ys[wiersz,kolumna])+'\t') fout.write('\n') fout.close() ########################################################################## def MakeODMRangleMap(Bmin=0.3,Bmax=0.85,Bstep=0.0002,BthetaDegMin=0,BthetaDegMax=0,BthetaDegStep=0,BphiDeg=0,MWfGHz=15.6,dfHz=0.7,ifPlot=False,ifSave=True,fname_additon=''): BXs=[] PZs=[] TYs=[] for th in np.arange(BthetaDegMin,BthetaDegMax+BthetaDegStep,BthetaDegStep): bx,pz=MakeSingleODMRspectrum(Bmin,Bmax,Bstep,th,BphiDeg,MWfGHz,dfHz,False,False) BXs.append(bx) PZs.append(pz) TYs.append(th) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") if ifSave else nic() if ifPlot: xlabel('Magnetic field B [T]') ylabel('Magnetic filed angle [deg]') title('ODMR angular map') imshow(PZs,aspect='auto',extent=[Bmin,Bmax,BthetaDegMin,BthetaDegMax],origin='lower', interpolation='none')#interpolation='spline36' savefig('ODMR_map_'+timestamp+'.pdf') if ifSave else nic() show() if ifSave: SaveThreeColumn(BXs,TYs,PZs,'Do_publikacji_ODMR_map_112300B_'+timestamp+'_'+fname_additon+'.txt') return BXs,TYs,PZs #----------------------------------------- def SaveThreeColumn(BXs,TYs,PZs,fname): osX=BXs[0] fout=open(fname,'w') fout.write('-\t') for x in osX: fout.write(str(x)+'\t') fout.write('\n') for t in range(len(TYs)): fout.write(str(TYs[t])) for i in range(len(osX)): fout.write('\t'+str(PZs[t][i])) fout.write('\n') fout.close() def GHz_na_T(GHz,g=2.0): return 0.071448*GHz/g def T_na_GHz(T,g=2.0): return T*g/(0.071448) ########################################################################## H1=HamiltonianJonuIS(DD=600.0e-9,EE=-0e-6) P1=np.kron(OneIJon,1*np.sign(SJonMi)) #dziala P1=np.kron(OneIJon,1*SJonMi) #dziala P2=np.kron(OneIJon,1*np.sign(SJonPl)) #nie dziala P2=np.kron(OneIJon,1*SJonPl) #nie dziala Boltzmann_dic = {} Boltzmann_dic['calc']=0 #H=HamiltonianJonuZeeman(0.5445,0,0) #MakeSpektrum(H+H1,P=(P1,P2),ifPlotSpectrum=True,ifPlotTransitions=True,ifPlotSave=True,dEbin=0.0001) #PlotEnergyLevelsVsB(Bmin=0.5,Bmax=0.6,Bstep=0.00001,BthetaDeg=0,BphiDeg=45,ifPlot=True,ifSave=True) #P1pocz,P1kon,P1wart,P2pocz,P2kon,P2wart=PrzejsciaRzadkie((P1,P2),tol=1e-8) #T=np.kron(IJonZ,OneSJon)+np.kron(OneIJon,SJonZ*1j)#####################3 #H=HamiltonianJonuZeeman(0.5445,0,0) #MakeSpektrum(H+H1,P=(P1,P2),ifPlotSpectrum=True,ifPlotTransitions=True,ifPlotSave=True,dEbin=0.0001) #BX,PX=MakeSingleODMRspectrum(Bmin=0.41,Bmax=0.51,Bstep=0.0001,BthetaDeg=0,BphiDeg=45,MWfGHz=12.9,dfHz=2,ifPlot=1) #bx,ey=PlotEnergyLevelsVsB(Bmin=0,Bmax=0.9,Bstep=0.0003,BthetaDeg=0,BphiDeg=0,ifPlot=True,ifSave=True) #BXs,TYs,PZs=MakeODMRangleMap(Bmin=0.0,Bmax=0.9,Bstep=0.0003,BthetaDegMin=0,BthetaDegMax=90,BthetaDegStep=2,BphiDeg=0,MWfGHz=15.6,dfHz=0.7,ifPlot=True,ifSave=True) #Box,Poy=MakeSingleODMRspectrum(Bmin=0.54,Bmax=0.58,Bstep=0.0002,BthetaDeg=0,MWfGHz=15.6,dfHz=0.7,ifPlot=True,ifSave=False) #H=HamiltonianJonuZeeman(0.08,0) #MakeSpektrum(H+H1,P=(P1,P2),ifPlotSpectrum=True,ifPlotTransitions=True,ifPlotSave=False,dEbin=0.0) #P=SJonPl #numbers=arange(1,9.99,0.05) #numbers=array(numbers) #Dec0=list(numbers*1e-9) #Dec1=list(numbers*1e-8) #Dec2=list(numbers*1e-7) #Decs=Dec0+Dec1+Dec2 #Elist=[150e-9,175e-9,200e-9,300e-9,100e-9,1e-9,50e-9,400e-9,600e-9,1000e-9] #Dlist=[1e-9,50e-9,100e-9,150e-9,200e-9,300e-9,400e-9,500e-9,800e-9,900e-9,1000e-9,2000e-9] #PhiList=[0,5,10,15,20,25,30,35,40,45] PhiList=[45]#[0,15,30,45,60,75,90,105,120,135,150,165,180]#arange(45,137.5,0.5) Dlist=list(array(arange(0,2752.5,2.5))*1e-9)#[50e-9,480e-9,1000e-9] Elist=[0]#[50e-9,150e-9,480e-9] AlphaList=[45]#[0,15,30,45] Tlist = list(array(arange(1.5,10.5,0.5)))+list(array(arange(10,110,10)))+list(array(arange(100,1100,100))) #========================= Bmin=0.2#0.28 Bmax=0.55#0.4748 Bstep=0.0002 BthetaDeg=0 MWfGHz=10.564 dfHz=3.4#1.7 T=2.75# #========================= def temp_dep_D550(T_K): '''skrajna prawa do skrajnej lewej UW1028 1.51/0.85=1.776 0.937/0.52=1.81 2.88 do 1.82 stosunek wniosek: T=1.95K do 3.55K''' r =(-2-1.3698*T_K)/(1-1.35423*T_K) return r if 1: Bstep=0.0002 ifPlot=False ifSave=False DD = 2800e-9#550e-9 EE = 0 alph = 45 BphiDeg = 45 T = 2.75#2.13*1000 fname_additon = '' H1=HamiltonianJonuIS(DD,EE,alph) T1=time.time() n=10 for i in range(n): BX,PX=MakeSingleODMRspectrum(Bmin,Bmax,Bstep,BthetaDeg,BphiDeg,MWfGHz,dfHz,ifPlot,ifSave,fname_additon,temperature =T) T2=time.time() print('-------------------\n\nMakeSingleODMRspectrum time secs',(T2-T1)/n) elif 0: for alph in AlphaList: for PH in PhiList: for EE in Elist: DDs=[] BXs=[] PZs=[] for DD in Dlist: print( "EE=",EE, " eV for D=", DD, "eV , Phi=",PH,"deg ; Theta=",str(BthetaDeg),"deg ; alphaE=",alph," deg") H1=HamiltonianJonuIS(DD,EE,alph) fnameadd='f'+str(MWfGHz)+'GHz'+'_dfHz'+str(dfHz)+'_DD'+str(DD)+'_EE'+str(EE)+'_eV_Phi'+str(PH)+'_deg_Theta'+str(BthetaDeg)+'_deg'+'_alph'+str(alph)+'_deg' BphiDeg=PH ifPlot=False ifSave=False fname_additon=fnameadd BX,PX=MakeSingleODMRspectrum(Bmin,Bmax,Bstep,BthetaDeg,BphiDeg,MWfGHz,dfHz,ifPlot,ifSave,fname_additon,T) BXs.append(BX) PZs.append(PX) DDs.append(DD) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") SaveThreeColumn(BXs,DDs,PZs,"W_funkcji_D_sigma_T_2.75_K_"+timestamp+'_'+fnameadd+'.txt') else: for alph in AlphaList: for PH in PhiList: for EE in Elist: EE = Elist[0] Ts=[] BXs=[] PZs=[] DD = 520e-9 for T in Tlist: print( "EE=",EE, " eV for D=", DD, "eV , Phi=",PH,"deg ; Theta=",str(BthetaDeg),"deg ; alphaE=",alph," deg") H1=HamiltonianJonuIS(DD,EE,alph) fnameadd='f'+str(MWfGHz)+'GHz'+'_dfHz'+str(dfHz)+'_DD'+str(DD)+'_EE'+str(EE)+'_eV_Phi'+str(PH)+'_deg_Theta'+str(BthetaDeg)+'_deg'+'_alph'+str(alph)+'_deg' BphiDeg=PH ifPlot=False ifSave=False fname_additon=fnameadd BX,PX=MakeSingleODMRspectrum(Bmin,Bmax,Bstep,BthetaDeg,BphiDeg,MWfGHz,dfHz,ifPlot,ifSave,fname_additon,T) BXs.append(BX) PZs.append(PX) Ts.append(T) timestamp=datetime.datetime.now().strftime("%Y%m%d-%H%M%S") SaveThreeColumn(BXs,Ts,PZs,"W_funkcji_T_Sigma_"+timestamp+'_'+fnameadd+'.txt')