What does this section of code do?

I'm adapting a pre-made set of code and was wodering what this section of it actually does.

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.

Full code shown below

 %% 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

Answer 1

[Walter Roberson]

it tries to solve nonlinear equations numerically over the range of values given in the last parameter .