How to calculate land surface temperture from Landsat8 imagery in Matlab?

I have downloaded some Landsat 8 imagery and I want to calculate the LST (Land Surface Temperature) from these images in MATLAB.

However, something does not seem to go right. For example, when I checked the calculated temperatures over a glacier surface, the displayed values are more than 44 degrees Celcius. I don't think that this is possible, as a melting ice surface is of course not able to get positive higher than 0 degrees. Can anyone help me to check where things go wrong? I have not been able to find the error yet.

The code I have used can be found under the image. The formulas I used were found on: https://giscrack.com/how-to-calculate-land-surface-temperature-with-landsat-8-images/. I have uploaded the corresponding geotiff images at: https://drive.google.com/file/d/1CXxt7hHWNbXujt6rMQPgGw103Zo6vbIc/view?usp=sharing.

Thanks!

clear all
close all

B4 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_SR_B4.tif'));
B5 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_SR_B5.tif'));
B10 = double(imread('LC08_L2SP_194028_20210706_20210713_02_T1_ST_B10.tif'));

%%

% Search for metafile and extract parameters

directory = pwd;                       % Full path of the directory to be searched in
filesAndFolders = dir('*MTL.txt');     % Returns all the files and folders in the directory
metadata_file = filesAndFolders.name;

fid  = fopen(metadata_file,'r');
text = textscan(fid,'%s','Delimiter','');
text = text{1};

% Band-specific multiplicative rescaling factor ML for band 10

idx = find(~cellfun('isempty',strfind(text,'RADIANCE_MULT_BAND_10')));
ML = string(text(idx));
ML = extractAfter(ML,"= ");
ML = str2double(ML);

% Band-specific additive rescaling factor AL from the metadata

idx = find(~cellfun('isempty',strfind(text,'RADIANCE_ADD_BAND_10')));
AL = string(text(idx));
AL = extractAfter(AL,"= ");
AL = str2double(AL);

% Band-specific thermal conversion constant from the metadata K1 from the metadata

idx = find(~cellfun('isempty',strfind(text,'K1_CONSTANT_BAND_10')));
K1 = string(text(idx));
K1 = extractAfter(K1,"= ");
K1 = str2double(K1);

% Band-specific thermal conversion constant from the metadata K1 from the metadata

idx = find(~cellfun('isempty',strfind(text,'K2_CONSTANT_BAND_10')));
K2 = string(text(idx));
K2 = extractAfter(K2,"= ");
K2 = str2double(K2);

clearvars idx fid metadata_file

%% Calculate TOA spectral radiance 
% TOA = ML * Qcal + AL 
%
% where 
% ML = Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number).
% Qcal = corresponds to band 10
% AL = Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number).

TOA = ML * B10 + AL;

%% TOA to Brightness Temperature conversion
% BT = (K2 / (ln (K1 / TOA) + 1)) - 273.15
%
% where 
% K1 = Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the thermal band number).
% K2 = Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the thermal band number).
% TOA (Top of Atmospheric) spectral radiance.

BT =  real(1321.0789 ./ log(774.8853./TOA +1.0)) - 273.15;

%% Calculate the NDVI
% NDVI = (Band 5 – Band 4) / (Band 5 + Band 4)

NDVI = (B5-B4)./(B5+B4);

%% Calculate the proportion of vegetation Pv
% Pv = ((NDVI – NDVImin) / (NDVImax – NDVImin))^2
% where 
% NDVI is the Normalized Difference Vegetation Index

PV = ((NDVI - min(NDVI(:))) ./ (max(NDVI(:)) - min(NDVI(:)))).^2;

%% Calculate Emissivity E
% E = 0.004 * Pv + 0.986
% where 
% Pv is the proportion of vegetation

E = (0.004 .* PV) + 0.986;

%% Calculate final Land Surface Temperature (LST)
% LST = (BT / (1 + (0.00115 * BT / 1.4388) * Ln(E)))
% where
% BT = Brightness Temperature
% E = Emissivity

LST = (BT ./ (1 + (0.00115 .* BT / 1.4388) .* log(E)));
LST(isnan(LST))=0;

%% Plot result

figure,
imagesc(LST),colorbar
caxis([0 80])