Add files via upload

Pot_MultiLayers.m

@@ -0,0 +1,53 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vb=0.3;                     % potential barrier [eV]
+
+%zt = [  10 10 10 9 10 8 10 7 10 6 10 ]*1E-9;  % thickness of each layer [nm]
+%Vt = [  1 0 1 0 1 0 1 0 1 0 1   ]*Vb;        % barrier height of each layer [eV]
+
+zt = [ 10 10 5 3 10 ]*1E-9;  % thickness of each layer [nm]
+Vt = [  1  0 1 0  1  ]*Vb;        % barrier height of each layer [eV]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Discretisation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% here, I descretize the grid z and the potential V0
+
+for i=1:length(zt)
+  if i==1
+    zv{1} = 0:dz:zt(1);
+    z=zv{1};
+    V0=zv{1}*0+Vt(1);
+  else
+    zv{i} = (z(end)+dz):dz:(z(end) + zt(i));
+    z  = [ z  zv{i} ];
+    V0 = [ V0   (zv{i}*0+1) * Vt(i)  ];
+  end
+end
+
+% here, I kind of patch the grid for the TMM
+
+if TM_Method==1
+  clear zv
+  for i=1:length(zt)
+    zzt((i-1)*Ns+1:i*Ns)=zt(i)/Ns;
+    VVt((i-1)*Ns+1:i*Ns)=Vt(i);  
+  end
+  zzt=[1e-12  zzt 1e-12];
+  VVt=[10 VVt 10];
+
+  for i=1:length(zzt)
+    if i==1
+      ZZ(1) = zzt(1);
+      zv{1} = 0:dz:zzt(1);
+      zz    = zv{1};
+    else
+      ZZ(i) = ZZ(end)+zzt(i);
+      zv{i} = (zz(end)+dz):dz:ZZ(end);
+      zz    = [zz zv{i}];
+    end
+  end
+
+end

Pot_Parabolic.m

@@ -0,0 +1,20 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vb=0.3;                 % potential barrier [eV]
+a=1e15;
+Ltot=50e-9;
+
+z=0:dz:Ltot;
+V0=a*(z-Ltot/2).^2;
+
+V0(V0>Vb)=Vb;
+
+if TM_Method==1
+  display('ERROR: TM_Method not valid for that potential')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pot_Sinus.m

@@ -0,0 +1,21 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vb=2;                 % potential barrier [eV]
+a=1e15;
+Ltot=50e-9;
+
+z=0:dz:Ltot;
+V0=Vb/2*(sin(7*pi/Ltot*z) +1);
+V0(1)=Vb;
+V0(end)=Vb;
+%V0(V0>Vb)=Vb;
+
+if TM_Method==1
+  display('ERROR: TM_Method not valid for that potential')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Schroed1D_Euler_f.m

@@ -0,0 +1,87 @@
+function[E,psi]=Schroed1D_Euler_f(z,V0,Mass,n,dE,precision)
+
+method=2;         % method "2" is much more acurate than "1" but take also more time...
+e=min(V0);
+Emax=max(V0)+0.1;
+
+C=0; N=0; E=[];
+psi=z*0+1; psi_old=psi;
+
+
+while e<Emax && length(E)<n
+
+    C = C+1;
+
+    psi_old = psi;
+    psi=Schroed1D_Euler_Eval(z,V0,Mass,e,method);
+
+    if (sign( psi(end) ) ~= sign( psi_old(end) ) ) && C>1 % here, I catch a quantum state because the last point of psi change sign
+
+        N=N+1; de=dE;
+
+        while abs(de)>precision
+            if sign( psi(end) ) ~= sign( psi_old(end) )
+                 de = -de/2;
+            end
+            e=e+de;
+            psi_old=psi;
+            psi=Schroed1D_Euler_Eval(z,V0,Mass,e,method);
+        end
+
+        E(N,:)=e; PSI(:,N)= psi;
+        C=0;
+    end	
+
+    e=e+dE;
+end
+
+psi=PSI;
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function[psi]=Schroed1D_Euler_Eval(z,V0,Mass,ee,method)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h=6.62606896E-34;               %% Planck constant [J.s]
+hbar=h/(2*pi);
+e=1.602176487E-19;              %% electron charge [C]
+me=9.10938188E-31;              %% electron mass [kg]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+dz=z(2)-z(1);
+y = [0 1];
+
+if method==1
+
+    for i = 1:length(z)                                  %% number of z steps to take
+        dy(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * y(1) ; %% Equation for dv/dz
+        dy(1) =  Mass*me*y(2);                           %% Equation for dx/dz
+        y = y + dz*dy ;                                  %% integrate both equations with Euler
+        psi(i)= y(1);
+    end
+
+elseif method==2
+
+    for i = 1:length(z)                                  %% number of z steps to take
+
+        dy(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * y(1) ; %% Equation for dv/dz
+        dy(1) =  Mass*me*y(2);                           %% Equation for dx/dz
+
+        K = y + 0.5*dz*dy;
+        dK(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * K(1) ; %% Equation for dv/dz
+        dK(1) =  Mass*me*K(2);                           %% Equation for dx/dz
+
+        y =  y + dz*dK;
+        psi(i)= y(1);
+    end
+end
+
+end

Schroed1D_FEM_f.m

@@ -0,0 +1,50 @@
+function[E,psi]=Schroed1D_FEM_f(z,V0,Mass,n)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h=6.62606896E-34;               %% Planck constant [J.s]
+hbar=h/(2*pi);
+e=1.602176487E-19;              %% electron charge [C]
+me=9.10938188E-31;              %% electron mass [kg]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nz=length(z);
+dz=z(2)-z(1);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%% Building of the operators %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DZ2 =(-2)*diag(ones(1,Nz)) + (1)*diag(ones(1,Nz-1),-1) + (1)*diag(ones(1,Nz-1),1);
+DZ2=DZ2/dz^2;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%% Building of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+H = -hbar^2/(2*me*Mass) * DZ2  +  diag(V0*e)  ;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% Diagonalisation of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+H = sparse(H);
+[psi,Energy] = eigs(H,n,'SM');
+E = diag(Energy)/e ;
+E=real(E);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Normalization of the Wavefunction %%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for i=1:n
+    psi(:,i)=psi(:,i)/sqrt(trapz(z',abs(psi(:,i)).^2));  % normalisation at 1
+end
+
+psi=psi(:,end:-1:1);
+E=E(end:-1:1);
+
+end

Schroed1D_PWE_f.m

@@ -0,0 +1,91 @@
+function[E,psi]=Schroed1D_PWE_f(z,V0,Mass,n,Nz,NG)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h=6.62606896E-34;               %% Planck constant [J.s]
+hbar=h/(2*pi);
+e=1.602176487E-19;              %% electron charge [C]
+me=9.10938188E-31;              %% electron mass [kg]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%% Interpolation on a grid that have 2^N points %%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+NG = 2*floor(NG/2);             % round to lower even number
+
+zz=linspace(z(1),z(end),Nz);
+V=interp1(z,V0,zz);
+
+dz=z(2)-z(1);
+dzz=zz(2)-zz(1);
+Ltot=zz(end)-zz(1);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Building Potential in Fourier space %%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vk = fftshift(fft(V))*dzz/Ltot;
+Vk =Vk(Nz/2-NG+1:Nz/2+NG+1);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% Reciprocal lattice vectors %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+G = (-NG/2:NG/2)*2*pi/Ltot;
+NG=length(G);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%% Building Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+idx = ones(NG,1)*(1:NG);        % matrix with elements equal to the column index idx(i,j) = j
+idx = idx' - idx;               % matrix with elements equal to row - column idx(i,j) = i - j;
+idx = idx + NG ;                % idx(i,j) = i - j + N
+
+D2=diag(G(1:NG)).^2;            % Operator second derivative
+
+H = hbar^2/(2*me*Mass)*D2 + Vk(idx)*e;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% Solving Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+H = sparse(H);
+[psik, Ek] = eigs(H,n,'SM');
+E = diag(Ek)  / e;
+E=abs(E);
+%E=real(E);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% Transforming & Scaling the waves functions %%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for i=1:n
+    psizz    = invFFT1D(psik(:,i).',Nz)/dzz ; % take care that the transpose does not take the complex conjugate
+    psi_temp = interp1(zz,psizz,z);
+    psi(:,i) = psi_temp/sqrt(trapz(z,abs(psi_temp).^2));  % normalisation at 1
+end
+
+% in Octave, the order of the eigen values are reversed...
+%psi=psi(:,end:-1:1);
+%E=E(end:-1:1);
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [Vz] = invFFT1D(Vk,Nz)
+
+Nk=length(Vk);
+N00=Nz/2-ceil(Nk/2);
+
+Vk00= [zeros(1,N00+1)  Vk  zeros(1,N00)] ; % adding zeros on the missing Fourier components
+Vz=ifft(ifftshift(Vk00));
+
+end