Plotting a Graph from Data

I'm new to MatLab and I'm modelling a semiconductor laser with the aim of plotting 'L' (quantum well width) against 'Band_Gap'. I've previously tried 'plot(L,Band_Gap);' in the command window but it just gives me a blank graph. How would i go about plotting this?

My current code

Answer 1

[Star Strider]

Try:

figure
plot(L,Band_Gap)
grid

in your script that calculates those values instead.

[SKC1997]

Hi,

I've just tried this and it is still returning a blank graph as shown. What am I doing wrong?

Cheers

[Star Strider]

I have no idea.

I cannot run your code, because for some reason, MATLAB will not run images of code. If you copy your code from your Editor and paste it to a Comment here, I will probably be able to run it, and I might be able to help.

[SKC1997]

This my code in it's current from. I basically want to be able to plot Band Gap against L for values of L between 0 and 11*10^-9.

%% Constants with all value containing Sb and all AlAs values removed
clear all; close all
%effective masses 
m_GaAs_SO = 0.17;
m_GaAs_lh = 0.073;
m_GaAs_hh001 = 0.35;
m_GaAs_hh111 = 0.87;
m_GaAs_c = 0.07;

m_InAs_SO = 0.15;
m_InAs_lh = 0.055;
m_InAs_hh001 = 0.39;
m_InAs_hh111 = 0.98;
m_InAs_c = 0.031;

%Added by Simon
m_InP_SO = 0.19;
m_InP_lh = 0.14;
m_InP_hh001 = 0.48;
m_InP_hh111 = 1.13;
m_InP_c = 0.064;

m_GaP_SO = 0.23;
m_GaP_lh = 0.16;
m_GaP_hh001 = 0.41;
m_GaP_hh111 = 1.01;
m_GaP_c = 0.15;

%Lattice Constants with InP and GaP
a_GaAs = 5.653;
a_InAs = 6.058;
a_InP = 5.896 ;
a_GaP = 5.451 ;

%Hydrostatic deformation potentials with InP and GaP
av_GaAs = 1.16 ;
av_InAs = 1.00 ;
av_InP = 1.27 ;
av_GaP = 1.70;


ac_GaAs = -7.17 ;
ac_InAs = -5.08 ;
ac_InP = -5.04 ;
ac_GaP = -7.14;

%Elastic Constants with InP and GaP
C11_GaAs = 1.18;
C12_GaAs = 0.54;
C44_GaAs = 0.59;

C11_InAs = 0.83;
C12_InAs = 0.45;
C44_InAs = 0.40;

C11_InP = 1.02;
C12_InP = 0.58;
C44_InP = 0.46;

C11_GaP = 1.41;
C12_GaP = 0.62;
C44_GaP = 0.70;

%Spin-orbit Splitting with InP and GaP
Delta_GaAs = 0.341 ;
Delta_InAs = 0.39 ;
Delta_InP = 0.108 ;
Delta_GaP = 0.08 ;

%Valence Band Energy (Valence Band Offset) with InP and GaP
EVBGaAs = -0.8 ;
EVBInAs = -0.59 ;
EVBInP = -0.94 ;
EVBGaP = -1.27 ;

%Gamma Band Gaps with InP and GaP
bGaAs = 1.519 ;
bInAs = 0.417 ;
bInP = 1.4236 ;
bGaP = 2.886 ;

%Conduction Band Energy (Sum of Valence band offset & Bandgap)
ECBGaAs = bGaAs+EVBGaAs ;
ECBInAs = bInAs+EVBInAs ;
ECBInP = bInP+EVBInP ;
ECBGaP = bGaP+EVBGaP ;

%Split-off Band Energy (Valence band energy minus Delta SO)
ESOGaAs = EVBGaAs - Delta_GaAs;
ESOInAs = EVBInAs - Delta_InAs;
ESOInP = EVBInP - Delta_InP;
ESOGaP = EVBGaP - Delta_GaP;
%% Parameters UNSURE
y = 1 ; % from page 5852 '... or Ga(0.47)In(0.53)As'
t = 1 ; % 'GaInP Crossover composition' pg 5840 - 0.7
k = 0 ; % 'alloy becomes indirect above 0.45' pg 5844
p = 1 ; % Arbitrary
%% Ternaries - All found on 1991 Krijn paper

ECBGaInAs = (1-y)*ECBGaAs + y*ECBInAs ;
ECBGaInP = (1-t)*ECBGaP + t*ECBInP ; %added by SC
ECBGaAsP = (1-k)*ECBGaAs + k*ECBGaP ; %SC
ECBInAsP = p*ECBInAs + (1-p)*ECBInP; %SC

EVBGaInAs = (1-y)*EVBGaAs + y*EVBInAs ;
EVBGaInP = (1-t)*EVBGaP + t*EVBInP ;%SC
EVBGaAsP = (1-k)*EVBGaAs + k*EVBGaP ; %SC
EVBInAsP = p*EVBInAs + (1-p)*EVBInP; %SC

ESOGaInAs = (1-y)*ESOGaAs + y*ESOInAs ;
ESOGaInP = (1-t)*ESOGaP + t*ESOInP ;%SC
ESOGaAsP = (1-k)*ESOGaAs + k*ESOGaP ; %SC
ESOInAsP = p*ESOInAs + (1-p)*ESOInP; %SC

a_GaInAs = (1-y)*a_GaAs + y*a_InAs ;
a_GaInP = (1-t)*a_GaP + t*a_InP ;%SC
a_GaAsP = (1-k)*a_GaAs + k*a_GaP ; %SC
a_InAsP = p*a_InAs + (1-p)*a_InP; %SC

ac_GaInAs = (1-y)*ac_GaAs + y*ac_InAs ;
ac_GaInP = (1-t)*ac_GaP + t*ac_InP ;%SC
ac_GaAsP = (1-k)*ac_GaAs + k*ac_GaP ; %SC
ac_InAsP = p*ac_InAs + (1-p)*ac_InP; %SC

av_GaInAs = (1-y)*av_GaAs + y*av_InAs ;
av_GaInP = (1-t)*av_GaP + t*av_InP ;%SC
av_GaAsP = (1-k)*av_GaAs + k*av_GaP ; %SC
av_InAsP = p*av_InAs + (1-p)*av_InP; %SC

C11_GaInAs = (1-y)*C11_GaAs + y*C11_InAs ;
C11_GaInP = (1-t)*C11_GaP + t*C11_InP ;%SC
C11_GaAsP = (1-k)*C11_GaAs + k*C11_GaP ; %SC
C11_InAsP = p*C11_InAs + (1-p)*C11_InP; %SC

C12_GaInAs = (1-y)*C12_GaAs + y*C12_InAs ;
C12_GaInP = (1-t)*C12_GaP + t*C12_InP ;%SC
C12_GaAsP = (1-k)*C12_GaAs + k*C12_GaP ; %SC
C12_InAsP = p*C12_InAs + (1-p)*C12_InP; %SC

C44_GaInAs = (1-y)*C44_GaAs + y*C44_InAs ;
C44_GaInP = (1-t)*C44_GaP + t*C44_InP ;%SC
C44_GaAsP = (1-k)*C44_GaAs + k*C44_GaP ; %SC
C44_InAsP = p*C44_InAs + (1-p)*C44_InP; %SC

%% GaInPAs Barrier Parameters 
w = 0.47 ; % Taken from 'laser structure PDF'
v = 0.24 ; % Taken from 'laser structure PDF'


%% Ga_v_In_1-v_P_w_As_1-w_ Barrier
%Rough attempt using GaInAsSb as a template for GaInPAs - Interpolation equation from Band Parameters paper

ECBGaInPAs = (v*(1-v)*(w*ECBGaInP + (1-w)*ECBGaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*ECBGaAsP + (1-v)*ECBInAsP))/(v*(1-v) + w*(1-w)) ;
% ^ Conduction band edge
EVBGaInPAs = (v*(1-v)*(w*EVBGaInP + (1-w)*EVBGaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*EVBGaAsP + (1-v)*EVBInAsP))/(v*(1-v) + w*(1-w)) ; 
% ^ Valence band edge
ESOGaInPAs = (v*(1-v)*(w*ESOGaInP + (1-w)*ESOGaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*ESOGaAsP + (1-v)*ESOInAsP))/(v*(1-v) + w*(1-w)) ;
% ^ Spin-orbit split-off band edge

ECBGaInAsVEC = ECBGaInAs ; %repmat(ECBAlGaAsSb,1,501) ;%going to use GaInAs in for calculation of band gaps NOT SURE
EVBGaInAsVEC = EVBGaInAs; %repmat(EVBAlGaAsSb,1,501) ;
ESOGaInAsVEC = ESOGaInAs ; %repmat(ESOAlGaAsSb,1,501) ;

a_GaInPAs = (v*(1-v)*(w*a_GaInP + (1-w)*a_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*a_GaAsP + (1-v)*a_InAsP))/(v*(1-v) + w*(1-w)) ;

ac_GaInPAs = (v*(1-v)*(w*ac_GaInP + (1-w)*ac_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*ac_GaAsP + (1-v)*ac_InAsP))/(v*(1-v) + w*(1-w)) ;
av_GaInPAs = (v*(1-v)*(w*av_GaInP + (1-w)*av_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*av_GaAsP + (1-v)*av_InAsP))/(v*(1-v) + w*(1-w)) ;

C11_GaInPAs = (v*(1-v)*(w*C11_GaInP + (1-w)*C11_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*C11_GaAsP + (1-v)*C11_InAsP))/(v*(1-v) + w*(1-w)) ;
C12_GaInPAs = (v*(1-v)*(w*C12_GaInP + (1-w)*C12_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*C12_GaAsP + (1-v)*C12_InAsP))/(v*(1-v) + w*(1-w)) ;
C44_GaInPAs = (v*(1-v)*(w*C44_GaInP + (1-w)*C44_GaInAs))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*C44_GaAsP + (1-v)*C44_InAsP))/(v*(1-v) + w*(1-w)) ;

%% Ga_a_In_a-1_P_b_As_1-b_ Well Parameters
a = 0.19 ; %Taken from 'Laser Structure' PDF
b= 0.76 ; 

%% Ga_a_In_a-1_P_b_As_1-b_ Well
ECBGaInPAs_Well = (a*(1-a)*(b*ECBGaInP + (1-b)*ECBGaInAs))/(a*(1-a) + b*(1-b)) + (b*(1-b)*(a*ECBGaAsP + (1-a)*ECBInAsP))/(a*(1-a) + b*(1-b)) ;
% ^ Conduction band edge
EVBGaInPAs_Well = (a*(1-a)*(b*EVBGaInP + (1-b)*EVBGaInAs))/(a*(1-a) + b*(1-b)) + (b*(1-b)*(a*EVBGaAsP + (1-a)*EVBInAsP))/(a*(1-a) + b*(1-b)) ; 
% ^ Valence band edge
ESOGaInPAs_Well = (a*(1-a)*(b*ESOGaInP + (1-b)*ESOGaInAs))/(a*(1-a) + b*(1-b)) + (b*(1-b)*(a*ESOGaAsP + (1-a)*ESOInAsP))/(a*(1-a) + b*(1-b)) ;
% ^ Spin-orbit split-off band edge
%% Masses etc copying from original model, not sure whether to use 0.5 as in original or keep parameter
m_GaInAs_c = (1-y)*m_GaAs_c + y*m_InAs_c ;
m_GaInP_c = (1-t)*m_GaP_c + t*m_InP_c ;
m_GaAsP_c = (1-k)*m_GaAs_c + k*m_GaP_c ; 
m_InAsP_c = p*m_InAs_c + (1-p)*m_InP_c; 

m_GaInAs_hh001 = (1-y)*m_GaAs_hh001 + y*m_InAs_hh001 ;
m_GaInP_hh001 = (1-t)*m_GaP_hh001 + t*m_InP_hh001 ;
m_GaAsP_hh001 = (1-k)*m_GaAs_hh001 + k*m_GaP_hh001 ; 
m_InAsP_hh001 = p*m_InAs_hh001 + (1-p)*m_InP_hh001;

m_GaInAs_SO = (1-y)*m_GaAs_SO + y*m_InAs_SO ;
m_GaInP_SO = (1-t)*m_GaP_SO + t*m_InP_SO ;
m_GaAsP_SO = (1-k)*m_GaAs_SO + k*m_GaP_SO ; 
m_InAsP_SO = p*m_InAs_SO + (1-p)*m_InP_SO; 

%now interpolate for masses of quaternary

m_GaInPAs_c = (v*(1-v)*(w*m_GaInP_c + (1-w)*m_GaInAs_c))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*m_GaAsP_c + (1-v)*m_InAsP_c))/(v*(1-v) + w*(1-w)) ;
m_GaInPAs_h = (v*(1-v)*(w*m_GaInP_hh001 + (1-w)*m_GaInAs_hh001))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*m_GaAsP_hh001 + (1-v)*m_InAsP_hh001))/(v*(1-v) + w*(1-w)) ;
m_GaInPAs_s = (v*(1-v)*(w*m_GaInP_SO + (1-w)*m_GaInAs_SO))/(v*(1-v) + w*(1-w)) + (w*(1-w)*(v*m_GaAsP_SO + (1-v)*m_InAsP_SO))/(v*(1-v) + w*(1-w)) ;
%% Lattice parameter
%SC - swapped GaInAsSb for GaInAsP and deleted Bismuth related formulae
av = av_GaInPAs ; % Hydrostatic deformation potential of valence band.
ac = ac_GaInPAs ; % Hydrostatic deformation potential of conduction band.
a_0 = a_InP ; % change to match substrate lattice constant. SC - Substrate in InP
C11 = C11_GaInPAs ;
C12 = C12_GaInPAs ;
C44 = C44_GaInPAs ;


D_001 = 2*(C12/C11);

D_111 = 2*((C11+2*C12-2*C44)/(C11+2*C12+4*C44));

a_parallel = a_0; %a_0 is the lattice constant of the substrate 


%Deleted E+a parallel, E+a perpendicular and DeltaE's as Bismuth related


%% Schroedinger for well energy levels 
%Reaplaced GaInAsSbBi with GaInPAs
format long
hbar = 6.582119514*10^-16 ; % reduced Planc constant in eVs
L = 6*10.^-9; % Quantum well width in metres
mC = m_GaInPAs_c*(5.686462*10^-12) ; % Effective mass in conduction band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
mH = m_GaInPAs_h*(5.686462*10^-12) ; % Effective mass in valence band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
mS = m_GaInPAs_s*(5.686462*10^-12) ; % Effective mass in split-off band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
UC = ECBGaInPAs - ECBGaInPAs_Well;  % Well depth for Conduction Band. SC - NOT SURE... ORIGINAL(UC = ECBGaInAsVEC - ECBGaInPAs)
UH = abs(EVBGaInPAs - EVBGaInPAs_Well); % Well depth for Valence Band. SC - NOT SURE... ORIGINAL(UH = EVBGaInPAs - EVBGaInAsVEC)
US = abs(ESOGaInPAs - ESOGaInPAs_Well);  % Well depth for Spin-Orbit Split-off Band. SC - NOT SURE... ORIGINAL(UH = EVBGaInPAs - EVBGaInAsVEC)

syms EC
EC1 = vpasolve( tan(sqrt((2*mC*EC*L^2)/hbar^2)) == (2*sqrt((UC-EC)*EC))/(2*EC-UC), EC, [0.0000001, UC]); 
% ^ Offset of first energy level in conduction band. May have to adjust
% factor under UC until returns an error, use a value smaller than that to
% get lowest energy level.
syms EV
EV1 = vpasolve( tan(sqrt((2*mH*EV*L^2)/hbar^2))  == (2*sqrt((UH-EV)*EV))/(2*EV-UH), EV, [0.0000001, UH]); 
% ^ Offset of first energy level in valence band.May have to adjust
% factor under UH until returns an error, use a value smaller than that to
% get lowest energy level.
syms ES
ES1 = vpasolve( tan(sqrt((2*mS*ES*L^2)/hbar^2)) == (2*sqrt((US-ES)*ES))/(2*ES-US), ES, [0.0000001, US]); 
% ^ Offset of first energy level in split-off band.May have to adjust
% factor under US until returns an error, use a value smaller than that to
% get lowest energy level.

EnC = double(EC1+ECBGaInPAs);  %SC - Swapped GaInAsSbBi for GaInPAs
EnV = double(EVBGaInPAs-EV1);  %SC - Swapped GaInAsSbBi for GaInPAs
EnS = double(ESOGaInPAs-ES1);  %SC - Swapped GaInAsSbBi for GaInPAs

Band_Gap = (EnC)-(EnV)

Spin_Orbit_Splitting = (EnV)-(EnS) ;

Wavelength = ((4.135667516*10^-15)*299792458)/Band_Gap

[Star Strider]

Both ‘L’ and ‘Band_Gap’ are single points:

L =
     6.000000000000001e-09

Band_Gap =
   0.808866986571398

I see only constants in your code, (when I look at the Workspace window), no vectors. If you want to plot a range of ‘Band_Gap’ values as a function of a range of ‘L’ values, you have to create them all as vectors. Single points will only display if you plot them with a marker:

figure
plot(L, Band_Gap, 'p')
grid

However, this would not be very informative with only one point.

[SKC1997]

Sorry, yes that is my original code with L as a constant when I was modelling for that well width. I now want to see how Band_Gap would vary with L, so L would be a variable. How would I go about creating them as vectors? I'm basically MatLab illiterate. Cheers.

[Star Strider]

Choose a range of values for ‘L’, then create a vector for it.

For example:

L = linspace(1, 10, 150)*1E-9;                               % <—— CREATE VECTOR ‘L’

I would then create ‘EC1’ and everything else that uses ‘L’ to be numeric funcitons of ‘L’. (I tried to do this myself, but your code is too difficult for me to follow.) This may require that you use numeric solvers such as fsolve instead of the symbolic vpasolve or solve, and use various values of ‘L’ to calculate ‘EV1’ and the others in a loop. That would be inefficient and time-consuming, so it is likely best that you devise a different approach.

First, eliminate the syms call, and all the symbolic solve calls, so all the variables are now numeric.

The last section of your code is now:

%% Schroedinger for well energy levels 
%Reaplaced GaInAsSbBi with GaInPAs
format long
hbar = 6.582119514E-16 ; % reduced Planc constant in eVs
Lv = linspace(1, 10, 150)*1E-9;                               % <—— CREATE VECTOR ‘L’
mC = m_GaInPAs_c*(5.686462E-12) ; % Effective mass in conduction band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
mH = m_GaInPAs_h*(5.686462E-12) ; % Effective mass in valence band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
mS = m_GaInPAs_s*(5.686462E-12) ; % Effective mass in split-off band. Numeric term to convert from multiples of electron mass to eVm^-2s^2 so equation has consistent units.
UC = ECBGaInPAs - ECBGaInPAs_Well;  % Well depth for Conduction Band. SC - NOT SURE... ORIGINAL(UC = ECBGaInAsVEC - ECBGaInPAs)
UH = abs(EVBGaInPAs - EVBGaInPAs_Well); % Well depth for Valence Band. SC - NOT SURE... ORIGINAL(UH = EVBGaInPAs - EVBGaInAsVEC)
US = abs(ESOGaInPAs - ESOGaInPAs_Well);  % Well depth for Spin-Orbit Split-off Band. SC - NOT SURE... ORIGINAL(UH = EVBGaInPAs - EVBGaInAsVEC)

for k1 = 1:numel(Lv)
    L = Lv(k1);
    EC1(k1) = fsolve(@(EC) tan(sqrt((2*mC*EC*L^2)/hbar^2)) - (2*sqrt((UC-EC)*EC))/(2*EC-UC), 1);
    EV1(k1) = fsolve(@(EV) tan(sqrt((2*mH*EV*L^2)/hbar^2)) - (2*sqrt((UH-EV)*EV))/(2*EV-UH), 1);
    ES1(k1) = fsolve(@(ES) tan(sqrt((2*mS*ES*L^2)/hbar^2)) - (2*sqrt((US-ES)*ES))/(2*ES-US), 1);
end
% ^ Offset of first energy level in conduction band. May have to adjust
% factor under UC until returns an error, use a value smaller than that to
% get lowest energy level.

% ^ Offset of first energy level in valence band.May have to adjust
% factor under UH until returns an error, use a value smaller than that to
% get lowest energy level.

% ^ Offset of first energy level in split-off band.May have to adjust
% factor under US until returns an error, use a value smaller than that to
% get lowest energy level.
EnC = double(EC1+ECBGaInPAs);  %SC - Swapped GaInAsSbBi for GaInPAs
EnV = double(EVBGaInPAs-EV1);  %SC - Swapped GaInAsSbBi for GaInPAs
EnS = double(ESOGaInPAs-ES1);  %SC - Swapped GaInAsSbBi for GaInPAs
Band_Gap = (EnC)-(EnV)
Spin_Orbit_Splitting = (EnV)-(EnS) ;
Wavelength = ((4.135667516E-15)*299792458)./Band_Gap

figure
plot(Lv, real(Band_Gap))
grid

This runs without error, and is actually faster than I thought it would be.

Experiment to get the result you want.

[SKC1997]

Thanks for the help, I'll have a fiddle around with it.

[Star Strider]

My pleasure.

If my Answer helped you solve your problem, please Accept it!