Temperature determination difference between codes that should do the same thing!

So, I do research of laser heating metals. We have this code (main code) that separates each hotspot into different qudrants to analyze. Occacsionally, one quadrant nees to be excluded from the calculations (2nd code). So the main code, uses the LT (left top) as the focus for analyzing. The 2nd code, needs to do the analysis with LD (Left down) as the focus instead adn exclude everything that has to do with LT.
So the issue is: the tempeartures calculated from the 2nd code are 300 to 400 K too low. (pics of ffure with each code.). If you have any isight I would greatly appreciate it!
Here is the normal code:
clear;close all;
[filename, pathname, filterindex] = uigetfile(‘*.txt’, ‘Enter input data file..’);
fullname = fullfile(pathname, filename);
mydata = importdata(fullname);
hotspot_ob=double(mydata);
scale= 120; % size of the box 90*90
subpixel=4;
hscale=scale/2;
[k,l]=size(hotspot_ob);
SR=ones(size(hotspot_ob));
box_Num = 1; % box 1 or 2
ND_Num = 0; % 0, .5, 1, 2
switch ND_Num
case 0
ND_col = 1;
case .5
ND_col = 2;
case 1
ND_col = 3;
case 2
ND_col = 4;
end

%LATEST from 08/09/2023
box(1).rs = [ 1.0 1.0 1.0 1.0
1.05 1.101 1.14 1.14
0.803 1.12 1.54 4.17
0.645 0.788 0.929 2.82];
box(2).rs = [1.0 1.0 1.0 1
0.296 0.312 0.325 0.424
0.331 0.457 0.632 1.735
0.245 0.294 0.348 1.077];
SR(1:fix(k/2),1:fix(l/2)) = box(box_Num).rs(1,ND_col);
SR(1:fix(k/2),fix(l/2)+1:l) = box(box_Num).rs(2,ND_col);
SR(fix(k/2)+1:k,1:fix(l/2)) = box(box_Num).rs(3,ND_col);
SR(fix(k/2)+1:k,fix(l/2)+1:l) = box(box_Num).rs(4,ND_col);
hotspot=hotspot_ob./SR;
[V_lt,H_lt]=find_max(hotspot_ob,1,fix(k/2),1,fix(l/2));
%V_lt=140;H_lt=166;
[V_rt,H_rt]=find_max(hotspot_ob,1,fix(k/2),fix(l/2)+1,l);
[V_ld,H_ld]=find_max(hotspot_ob,fix(k/2)+1,k,1,fix(l/2));
[V_rd,H_rd]=find_max(hotspot_ob,fix(k/2)+1,k,fix(l/2)+1,l);
%%find max for hotspot for later T correction Jie add on Mar2216, two ways:
%%neigther seem to work properly
%1st find the strongest point and make it as threshold–find corresponding
%intensities in other quadrants
hotspot_lt = hotspot(1:fix(k/2),1:fix(l/2)); [k_lt, l_lt] = size(hotspot_lt);
hotspot_rt = hotspot(1:fix(k/2),(fix(l/2)+1):l); [k_rt, l_rt] = size(hotspot_rt);
hotspot_ld = hotspot((fix(k/2)+1):k,1:fix(l/2)); [k_ld, l_ld] = size(hotspot_ld);
hotspot_rd = hotspot((fix(k/2)+1):k,(fix(l/2)+1):l); [k_rd, l_rd] = size(hotspot_rd);
[m_lt,n_lt]=find_max(hotspot_lt,1,k_lt,1,l_lt);
[m_rt,n_rt]=find_max(hotspot_rt,1,k_rt,1,l_rt);
[m_ld,n_ld]=find_max(hotspot_ld,1,k_ld,1,l_ld);
[m_rd,n_rd]=find_max(hotspot_rd,1,k_rd,1,l_rd);
misfit = zeros(4,1);
misfit(1,1) = (hotspot_rt(m_rt,n_rt)-hotspot_rt(m_lt,n_lt)).^2+(hotspot_ld(m_ld,n_ld)-hotspot_ld(m_lt,n_lt)).^2 …
+(hotspot_rd(m_rd,n_rd)-hotspot_rd(m_lt,n_lt)).^2;
misfit(2,1) = (hotspot_lt(m_rt,n_rt)-hotspot_lt(m_rt,n_rt)).^2+(hotspot_ld(m_rt,n_rt)-hotspot_ld(m_rt,n_rt)).^2 …
+(hotspot_rd(m_rt,n_rt)-hotspot_rd(m_rt,n_rt)).^2;
misfit(3,1) = (hotspot_lt(m_ld,n_ld)-hotspot_lt(m_ld,n_ld)).^2+(hotspot_rt(m_ld,n_ld)-hotspot_rt(m_ld,n_ld)).^2 …
+(hotspot_rd(m_ld,n_ld)-hotspot_rd(m_ld,n_ld)).^2;
misfit(4,1) = (hotspot_lt(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2+(hotspot_rt(m_rd,n_rd)-hotspot_rt(m_rd,n_rd)).^2 …
+(hotspot_ld(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2;
switch min(misfit)
case misfit(1,1)
mm = m_lt; nn = n_lt;
case misfit(2,1)
mm = m_rt; nn = n_rt;
case misfit(3,1)
mm = m_ld; nn = n_ld;
case misfit(4,1)
mm = m_rd; nn = n_rd;
end
Max_color = [hotspot_lt(mm,nn) hotspot_rt(mm,nn) hotspot_ld(mm,nn) hotspot_rd(mm,nn)]
m = min([k_lt;k_rt;k_ld;k_rd]); n = min([l_lt;l_rt;l_ld;l_rd]);
misfit2 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit2(countm,countn) = (hotspot_lt(countm,countn)-hotspot_lt(m_lt,n_lt)).^2 + (hotspot_rt(countm,countn)-hotspot_rt(m_rt,n_rt)).^2 …
+(hotspot_ld(countm,countn)-hotspot_ld(m_ld,n_ld)).^2+(hotspot_rd(countm,countn)-hotspot_rd(m_rd,n_rd)).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit2(countm,countn) == min(min(misfit2))
mmm = countm; nnn = countn;
end
end
end
Max_color2 = [hotspot_lt(mmm,nnn) hotspot_rt(mmm,nnn) hotspot_ld(mmm,nnn) hotspot_rd(mmm,nnn)]
%%
center_LT=[V_lt,H_lt];
center_LD=[V_ld,H_ld];
center_RT=[V_rt,H_rt];
center_RD=[V_rd,H_rd];
% Define the bounds of the matrix
k_max = size(hotspot, 1);
l_max = size(hotspot, 2);
% Ensure the indices are within valid bounds
V_lt_start = max(V_lt-hscale, 1);
V_lt_end = min(V_lt+hscale, k_max);
H_lt_start = max(H_lt-hscale, 1);
H_lt_end = min(H_lt+hscale, l_max);
V_ld_start = max(V_ld-hscale, 1);
V_ld_end = min(V_ld+hscale, k_max);
H_ld_start = max(H_ld-hscale, 1);
H_ld_end = min(H_ld+hscale, l_max);
V_rt_start = max(V_rt-hscale, 1);
V_rt_end = min(V_rt+hscale, k_max);
H_rt_start = max(H_rt-hscale, 1);
H_rt_end = min(H_rt+hscale, l_max);
V_rd_start = max(V_rd-hscale, 1);
V_rd_end = min(V_rd+hscale, k_max);
H_rd_start = max(H_rd-hscale, 1);
H_rd_end = min(H_rd+hscale, l_max);
% Now extract the subarrays with valid indices
LT = hotspot(V_lt_start:V_lt_end, H_lt_start:H_lt_end);
LD = hotspot(V_ld_start:V_ld_end, H_ld_start:H_ld_end);
RT = hotspot(V_rt_start:V_rt_end, H_rt_start:H_rt_end);
RD = hotspot(V_rd_start:V_rd_end, H_rd_start:H_rd_end);
% Find New centers and update them
[center_RT_new,Cor_RT] = imagematching (hotspot,LT,center_RT,scale); % use LT as target image
[center_LD_new,Cor_LD] = imagematching (hotspot,LT,center_LD,scale);
[center_RD_new,Cor_RD] = imagematching (hotspot,LT,center_RD,scale);
center_LD=center_LD_new;
center_RT=center_RT_new;
center_RD=center_RD_new;
[inter_center_RT,Cor_RT] = intermatching (hotspot,LT,center_RT_new,scale,subpixel);
[inter_center_LD,Cor_LD] = intermatching (hotspot,LT,center_LD_new,scale,subpixel);
[inter_center_RD,Cor_RD] = intermatching (hotspot,LT,center_RD_new,scale,subpixel);
center_LD=inter_center_LD;
center_RT=inter_center_RT;
center_RD=inter_center_RD;
Cor=[1 Cor_RT;Cor_LD Cor_RD];
V_ld=center_LD(1,1); H_ld=center_LD(1,2);
V_rt=center_RT(1,1); H_rt=center_RT(1,2);
V_rd=center_RD(1,1); H_rd=center_RD(1,2);
[XI_ld,YI_ld]=meshgrid(H_ld-hscale:H_ld+hscale,V_ld-hscale:V_ld+hscale);
[XI_rt,YI_rt]=meshgrid(H_rt-hscale:H_rt+hscale,V_rt-hscale:V_rt+hscale);
[XI_rd,YI_rd]=meshgrid(H_rd-hscale:H_rd+hscale,V_rd-hscale:V_rd+hscale);
LD=interp2(hotspot,XI_ld,YI_ld); % update images
RT=interp2(hotspot,XI_rt,YI_rt);
RD=interp2(hotspot,XI_rd,YI_rd);
lamda_LT=580*10^(-9);
lamda_RT=640*10^(-9);
lamda_LD=766*10^(-9);
lamda_RD=905*10^(-9);
planck=6.6*10^(-34);
Kb=1.38*10^(-23);
c=3*10^8;
Wien_LT=Kb/planck/c*log(2*pi*planck*c*c./LT/lamda_LT^5);
Wien_RT=Kb/planck/c*log(2*pi*planck*c*c./RT/lamda_RT^5);
Wien_LD=Kb/planck/c*log(2*pi*planck*c*c./LD/lamda_LD^5);
Wien_RD=Kb/planck/c*log(2*pi*planck*c*c./RD/lamda_RD^5);
[m,n]=size(LT);
for i=1:m
for j=1:n
X=[1/lamda_LT 1/lamda_RT 1/lamda_LD 1/lamda_RD]’;
% X=[1/lamda_650 1/lamda_800 1/lamda_900]’;
X_reg=[X ones(size(X))];
Y=[Wien_LT(i,j) Wien_RT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
% Y=[Wien_LT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
%R(i,j)=corrcoef(X,Y);
[T_inv,bint]=regress(Y,X_reg); %% bint is 95% confidence level of T
R(i,j)=corr(X,Y);
T(i,j)=1/T_inv(1,1); %% T_inv(1,1) is for T
Em(i,j)=exp(-T_inv(2,1)*planck*c/Kb); %% T_inv(2,1) is for emissivity
T_err(i,j)=(1/bint(1,1)-1/bint(1,2))/2; %% difference between upper limit of T – lower limit of T
T_low(i,j)=1/bint(1,2);
T_high(i,j)=1/bint(1,1);

end
end
%% max in four regimes Jie add on Dec 5th 2015
[max_LT,Indexing] = max(LT(:))
max_RD = max(max(RD));max_LD = max(max(LD));max_RT = max(max(RT));
Max_color4 = [LT(Indexing) RT(Indexing) LD(Indexing) RD(Indexing)]
%% max grid search Jie add on April 1st
[m,n]=size(RD);
misfit3 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit3(countm,countn) = (LT(countm,countn)-max_LT).^2 + (RT(countm,countn)-max_RT).^2 …
+(LD(countm,countn)-max_LD).^2+(RD(countm,countn)-max_RD).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit3(countm,countn) == min(min(misfit3))
m4 = countm; n4 = countn;
end
end
end
Max_color3 = [LT(m4,n4) RT(m4,n4) LD(m4,n4) RD(m4,n4)];
Max_Ih = [LT(m4,:); RT(m4,:); LD(m4,:); RD(m4,:)];
Max_Iv = [LT(:,n4)’; RT(:,n4)’; LD(:,n4)’; RD(:,n4)’];
%%
figure(‘Name’,’hotspot’);
imagesc(hotspot_ob);
figure(‘Name’,’hot spot in 4 regions’);
subplot(2,2,1)
imagesc(LT);
subplot(2,2,2)
imagesc(RT);
subplot(2,2,3)
imagesc(LD);
subplot(2,2,4)
imagesc(RD);
figure(‘Name’,’T profile (pixel)’);
imagesc(T);
[Max_T, H]=max(max(T)); % maximum Temperature, horizontal position
[Max_T, V]=max(max(T’)); % maximum Temperature, vertical position
grid=0.48; % each pixel is 0.5um
X_position=(-H+1)*grid: grid: (scale-H+1)*grid;
figure(‘Name’,’Hirizontal T profile throughh highest T point’)
plot (X_position, T(V, :));
title(‘horizotal plot’)
xlabel(‘horizontal positon/um’)
ylabel (‘temperature/K’);
figure(‘Name’,’Vertical T profile throughh highest T point’)
Y_position=(-V+1)*grid: grid: (scale-V+1)*grid;
plot(Y_position,T(:,H));
title(‘Vertical plot’)
xlabel(‘Vertical positon/um’)
ylabel (‘temperature/K’);
[Y_p,X_p]=meshgrid(X_position, Y_position);
figure(‘Name’,’3D T profile’)
surf(X_p,Y_p,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
zlabel(‘Temperature(K)’)
shading interp;
figure(‘Name’,’T profile (pixel)’);
imagesc(X_position,Y_position,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
colormap(‘turbo’)
colorbar;
figure(‘Name’,’T profile (pixel) contour’);
imagesc(X_position,Y_position,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
colormap(‘turbo’)
colorbar;
hold on; % Keep the image displayed while adding contour
contour(X_position, Y_position, T, 10, ‘LineWidth’, 1, ‘LineColor’, ‘w’); % 40 contours, white lines
hold off; % Release the hold on the current figure
figure(‘Name’,’T error’);
imagesc(X_position,Y_position,T_err);
figure(‘Name’,’Emissivity Profile’)
imagesc(X_position,Y_position,log10(Em));
figure(‘Name’,’3D Emissivity profile’)
surf(X_p,Y_p,log10(Em));
H=scale/2; V=scale/2;% give the point needed to plot T cross section
figure(‘Name’,’Vertical Em/Emax vs T trough highest T point’)
plot(T(:,H),Em(:,H)/max(Em(:,H)),’–rs’);
title(‘Vertical plot’)
xlabel(‘Temperature’)
ylabel (‘Em/Em_max’);
figure(‘Name’,’horizontal Em/Emax vs T trough highest T point’)
plot(T(V,:),Em(V,:)/max(Em(V,:)),’–rs’);
title(‘horizontal Em/Emax vs T’)
xlabel(‘Temperature’)
ylabel (‘Em/Em_max’);
% figure(12); % 3D image for intensity in 4 regions
% subplot(2,2,1)
% surf(X_p,Y_p,LT);
% subplot(2,2,2)
% surf(X_p,Y_p,RT);
% subplot(2,2,3)
% surf(X_p,Y_p,LD);
% subplot(2,2,4)
% surf(X_p,Y_p,RD)
figure (‘Name’,’Emissivity vs T over whole space’); plot(T,Em);

Image result:

2nd code:
clear; close all
[filename, pathname, filterindex]=uigetfile(‘*.txt’, ‘Enter input data file..’);
fullname=fullfile(pathname, filename);
mydata=importdata(fullname);
hotspot_ob=double(mydata);
scale=120;
subpixel=4;
hscale=scale/2;
[k,l]=size(hotspot_ob);
SR=ones(size(hotspot_ob));
box_Num = 1; % box 1 or 2
ND_Num = 0; % 0, .5, 1, 2
switch ND_Num
case 0
ND_col = 1;
case .5
ND_col = 2;
case 1
ND_col = 3;
case 2
ND_col = 4;
end

%LATEST from 08/09/2023
box(1).rs = [ 1.0 1.0 1.0 1.0
1.05 1.101 1.14 1.14
0.803 1.12 1.54 4.17
0.645 0.788 0.929 2.82];
box(2).rs = [1.0 1.0 1.0 1
0.296 0.312 0.325 0.424
0.331 0.457 0.632 1.735
0.245 0.294 0.348 1.077];
SR(1:fix(k/2),1:fix(l/2)) = box(box_Num).rs(1,ND_col);
SR(1:fix(k/2),fix(l/2)+1:l) = box(box_Num).rs(2,ND_col);
SR(fix(k/2)+1:k,1:fix(l/2)) = box(box_Num).rs(3,ND_col);
SR(fix(k/2)+1:k,fix(l/2)+1:l) = box(box_Num).rs(4,ND_col);
hotspot=hotspot_ob./SR;
[V_lt,H_lt]=find_max(hotspot_ob,1,fix(k/2),1,fix(l/2));
[V_rt,H_rt]=find_max(hotspot_ob,1,fix(k/2),fix(l/2)+1,l);
[V_ld,H_ld]=find_max(hotspot_ob,fix(k/2)+1,k,1,fix(l/2));
[V_rd,H_rd]=find_max(hotspot_ob,fix(k/2)+1,k,fix(l/2)+1,l);
hotspot_lt = hotspot(1:fix(k/2),1:fix(l/2)); [k_lt, l_lt] = size(hotspot_lt);
hotspot_rt = hotspot(1:fix(k/2),(fix(l/2)+1):l); [k_rt, l_rt] = size(hotspot_rt);
hotspot_ld = hotspot((fix(k/2)+1):k,1:fix(l/2)); [k_ld, l_ld] = size(hotspot_ld);
hotspot_rd = hotspot((fix(k/2)+1):k,(fix(l/2)+1):l); [k_rd, l_rd] = size(hotspot_rd);
%V_lt=98; H_lt=150;
[m_lt,n_lt]=find_max(hotspot_lt,1,k_lt,1,l_lt);
[m_rt,n_rt]=find_max(hotspot_rt,1,k_rt,1,l_rt);
[m_ld,n_ld]=find_max(hotspot_ld,1,k_ld,1,l_ld);
[m_rd,n_rd]=find_max(hotspot_rd,1,k_rd,1,l_rd);
misfit = zeros(4,1);
misfit(1,1) = (hotspot_rt(m_rt,n_rt)-hotspot_rt(m_lt,n_lt)).^2+(hotspot_ld(m_ld,n_ld)-hotspot_ld(m_lt,n_lt)).^2 …
+(hotspot_rd(m_rd,n_rd)-hotspot_rd(m_lt,n_lt)).^2;
misfit(2,1) = (hotspot_lt(m_rt,n_rt)-hotspot_lt(m_rt,n_rt)).^2+(hotspot_ld(m_rt,n_rt)-hotspot_ld(m_rt,n_rt)).^2 …
+(hotspot_rd(m_rt,n_rt)-hotspot_rd(m_rt,n_rt)).^2;
misfit(3,1) = (hotspot_lt(m_ld,n_ld)-hotspot_lt(m_ld,n_ld)).^2+(hotspot_rt(m_ld,n_ld)-hotspot_rt(m_ld,n_ld)).^2 …
+(hotspot_rd(m_ld,n_ld)-hotspot_rd(m_ld,n_ld)).^2;
misfit(4,1) = (hotspot_lt(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2+(hotspot_rt(m_rd,n_rd)-hotspot_rt(m_rd,n_rd)).^2 …
+(hotspot_ld(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2;
switch min(misfit)
case misfit(1,1)
mm = m_lt; nn = n_lt;
case misfit(2,1)
mm = m_rt; nn = n_rt;
case misfit(3,1)
mm = m_ld; nn = n_ld;
case misfit(4,1)
mm = m_rd; nn = n_rd;
end
Max_color = [ hotspot_rt(mm,nn) hotspot_ld(mm,nn) hotspot_rd(mm,nn)]
m = min([k_lt;k_rt;k_ld;k_rd]); n = min([l_lt;l_rt;l_ld;l_rd]);
misfit2 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit2(countm,countn) = (hotspot_rt(countm,countn)-hotspot_rt(m_rt,n_rt)).^2 …
+(hotspot_ld(countm,countn)-hotspot_ld(m_ld,n_ld)).^2+(hotspot_rd(countm,countn)-hotspot_rd(m_rd,n_rd)).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit2(countm,countn) == min(min(misfit2))
mmm = countm; nnn = countn;
end
end
end
Max_color2 = [hotspot_rt(mmm,nnn) hotspot_ld(mmm,nnn) hotspot_rd(mmm,nnn)]
%V_lt=126; H_lt=176;
%V_rt=125; H_rt=565;
%V_ld=398; H_ld=179;
%V_rd=388; H_rd=580;
center_LT=[V_lt,H_lt];
center_LD=[V_ld, H_ld];
center_RT=[V_rt,H_rt];
center_RD=[V_rd,H_rd];
LT=hotspot(V_lt-hscale:V_lt+hscale,H_lt-hscale:H_lt+hscale);
LD=hotspot(V_ld-hscale:V_ld+hscale,H_ld-hscale:H_ld+hscale);
RT=hotspot(V_rt-hscale:V_rt+hscale,H_rt-hscale:H_rt+hscale);
RD=hotspot(V_rd-hscale:V_rd+hscale,H_rd-hscale:H_rd+hscale);
% Find New centers and update them
[center_LT_new,Cor_LT] = imagematching (hotspot,RD,center_RT,scale); % use RT as target image
[center_LD_new,Cor_LD] = imagematching (hotspot,RD,center_LD,scale);
[center_RT_new,Cor_RT] = imagematching (hotspot,RD,center_RD,scale);
center_LD=center_LD_new;
center_LT=center_LT_new;
center_RT=center_RT_new;
[inter_center_LT,Cor_LT] = intermatching (hotspot,RD,center_LT_new,scale,subpixel);
[inter_center_LD,Cor_LD] = intermatching (hotspot,RD,center_LD_new,scale,subpixel);
[inter_center_RT,Cor_RT] = intermatching (hotspot,RD,center_RT_new,scale,subpixel);
center_LD=inter_center_LD;
center_LT=inter_center_LT;
center_RT=inter_center_RT;
Cor=[1 Cor_LT;Cor_LD Cor_RT];
V_ld=center_LD(1,1); H_ld=center_LD(1,2);
V_rt=center_LT(1,1); H_lt=center_LT(1,2);
V_rd=center_RD(1,1); H_rd=center_RD(1,2);
[XI_ld,YI_ld]=meshgrid(H_ld-hscale:H_ld+hscale,V_ld-hscale:V_ld+hscale);
[XI_lt,YI_lt]=meshgrid(H_lt-hscale:H_lt+hscale,V_lt-hscale:V_lt+hscale);
[XI_rd,YI_rd]=meshgrid(H_rd-hscale:H_rd+hscale,V_rd-hscale:V_rd+hscale);
LD=interp2(hotspot,XI_ld,YI_ld,’spline’); % update images
LT=interp2(hotspot,XI_lt,YI_lt,’spline’);
RD=interp2(hotspot,XI_rd,YI_rd,’spline’);
% Assuming you want to exclude 580 nm
lamda_LT=580*10^(-9);
lamda_RT=640*10^(-9);
lamda_LD=766*10^(-9);
lamda_RD=905*10^(-9);
planck=6.6*10^(-34);
Kb=1.38*10^(-23);
c=3*10^8;
%Wien_LT=Kb/planck/c*log(2*pi*planck*c*c./LT/lamda_LT^5);
Wien_RT=Kb/planck/c*log(2*pi*planck*c*c./RT/lamda_RT^5);
Wien_LD=Kb/planck/c*log(2*pi*planck*c*c./LD/lamda_LD^5);
Wien_RD=Kb/planck/c*log(2*pi*planck*c*c./RD/lamda_RD^5);
[m,n]=size(LD);
for i=1:m
for j=1:n
X=[1/lamda_RT 1/lamda_LD 1/lamda_RD]’;
%X=[1/lamda_650 1/lamda_800 1/lamda_900]’;
X_reg=[X ones(size(X))];
%Y=[Wien_LT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
Y=[Wien_RT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
%R(i,j)=corrcoef(X,Y);
[T_inv,bint]=regress(Y,X_reg);
R(i,j)=corr(X,Y);
T(i,j)=1/T_inv(1,1);
Em(i,j)=exp(-T_inv(2,1)*planck*c/Kb);
T_err(i,j)=(1/bint(1,1)-1/bint(1,2))/2;
T_low(i,j)=1/bint(1,2);
T_high(i,j)=1/bint(1,1);

end
end
%% max in four regimes Jie add on Dec 5th 2015
[max_RD,Indexing] = max(RD(:))
max_LD = max(max(LD));max_RT = max(max(RT));
Max_color4 = [RT(Indexing) LD(Indexing) RD(Indexing)]
%% max grid search Jie add on April 1st
[m,n]=size(RD);
misfit3 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit3(countm,countn) = (RT(countm,countn)-max_RT).^2 …
+(LD(countm,countn)-max_LD).^2+(RD(countm,countn)-max_RD).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit3(countm,countn) == min(min(misfit3))
m4 = countm; n4 = countn;
end
end
end
Max_color3 = [ RT(m4,n4) LD(m4,n4) RD(m4,n4)];
Max_Ih = [ RT(m4,:); LD(m4,:); RD(m4,:)];
Max_Iv = [ RT(:,n4)’; LD(:,n4)’; RD(:,n4)’];
figure(‘Name’,’hotspot’);
imagesc(hotspot_ob);
figure(‘Name’,’hot spot in 4 regions’);
subplot(2,2,2)
imagesc(RT);
subplot(2,2,3)
imagesc(LD);
subplot(2,2,4)
imagesc(RD);
figure(‘Name’,’T profile (pixel)’);
imagesc(T);
[Max_T, H]=max(max(T)); % maximum Temperature, horizontal position
[Max_T, V]=max(max(T’)); % maximum Temperature, vertical position
grid=0.48; % each pixel is 0.5um
X_position=(-H+1)*grid: grid: (scale-H+1)*grid;
figure(‘Name’,’Hirizontal T profile throughh highest T point’)
plot (X_position, T(V, :));
title(‘horizotal plot’)
xlabel(‘horizontal positon/um’)
ylabel (‘temperature/K’);
figure(‘Name’,’Vertical T profile throughh highest T point’)
Y_position=(-V+1)*grid: grid: (scale-V+1)*grid;
plot(Y_position,T(:,H));
title(‘Vertical plot’)
xlabel(‘Vertical positon/um’)
ylabel (‘temperature/K’);
[Y_p,X_p]=meshgrid(X_position, Y_position);
figure(‘Name’,’3D T profile’)
surf(X_p,Y_p,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
zlabel(‘Temperature(K)’)
shading interp;
figure(‘Name’,’T profile (pixel)’);
imagesc(X_position,Y_position,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
colormap(‘turbo’)
colorbar;
figure(‘Name’,’T profile (pixel) contour’);
imagesc(X_position,Y_position,T);
xlabel(‘Horizontal distance(um)’)
ylabel (‘Vetical distance(um)’);
colormap(‘turbo’)
colorbar;
hold on; % Keep the image displayed while adding contour
contour(X_position, Y_position, T, 10, ‘LineWidth’, 1, ‘LineColor’, ‘w’); % 40 contours, white lines
hold off; % Release the hold on the current figure
figure(‘Name’,’T error’);
imagesc(X_position,Y_position,T_err);
figure(‘Name’,’Emissivity Profile’)
imagesc(X_position,Y_position,log10(Em));
figure(‘Name’,’3D Emissivity profile’)
surf(X_p,Y_p,log10(Em));
H=scale/2; V=scale/2;% give the point needed to plot T cross section
figure(‘Name’,’Vertical Em/Emax vs T trough highest T point’)
plot(T(:,H),Em(:,H)/max(Em(:,H)),’–rs’);
title(‘Vertical plot’)
xlabel(‘Temperature’)
ylabel (‘Em/Em_max’);
figure(‘Name’,’horizontal Em/Emax vs T trough highest T point’)
plot(T(V,:),Em(V,:)/max(Em(V,:)),’–rs’);
title(‘horizontal Em/Emax vs T’)
xlabel(‘Temperature’)
ylabel (‘Em/Em_max’);
% figure(12); % 3D image for intensity in 4 regions
% subplot(2,2,1)
% surf(X_p,Y_p,LT);
% subplot(2,2,2)
% surf(X_p,Y_p,RT);
% subplot(2,2,3)
% surf(X_p,Y_p,LD);
% subplot(2,2,4)
% surf(X_p,Y_p,RD)
figure (‘Name’,’Emissivity vs T over whole space’); plot(T,Em);

image result:So, I do research of laser heating metals. We have this code (main code) that separates each hotspot into different qudrants to analyze. Occacsionally, one quadrant nees to be excluded from the calculations (2nd code). So the main code, uses the LT (left top) as the focus for analyzing. The 2nd code, needs to do the analysis with LD (Left down) as the focus instead adn exclude everything that has to do with LT.
So the issue is: the tempeartures calculated from the 2nd code are 300 to 400 K too low. (pics of ffure with each code.). If you have any isight I would greatly appreciate it!
Here is the normal code:
clear;close all;
[filename, pathname, filterindex] = uigetfile(‘*.txt’, ‘Enter input data file..’);
fullname = fullfile(pathname, filename);
mydata = importdata(fullname);
hotspot_ob=double(mydata);
scale= 120; % size of the box 90*90
subpixel=4;
hscale=scale/2;
[k,l]=size(hotspot_ob);
SR=ones(size(hotspot_ob));
box_Num = 1; % box 1 or 2
ND_Num = 0; % 0, .5, 1, 2
switch ND_Num
case 0
ND_col = 1;
case .5
ND_col = 2;
case 1
ND_col = 3;
case 2
ND_col = 4;
end

Image result:

%LATEST from 08/09/2023
box(1).rs = [ 1.0 1.0 1.0 1.0
1.05 1.101 1.14 1.14
0.803 1.12 1.54 4.17
0.645 0.788 0.929 2.82];
box(2).rs = [1.0 1.0 1.0 1
0.296 0.312 0.325 0.424
0.331 0.457 0.632 1.735
0.245 0.294 0.348 1.077];
SR(1:fix(k/2),1:fix(l/2)) = box(box_Num).rs(1,ND_col);
SR(1:fix(k/2),fix(l/2)+1:l) = box(box_Num).rs(2,ND_col);
SR(fix(k/2)+1:k,1:fix(l/2)) = box(box_Num).rs(3,ND_col);
SR(fix(k/2)+1:k,fix(l/2)+1:l) = box(box_Num).rs(4,ND_col);
hotspot=hotspot_ob./SR;
[V_lt,H_lt]=find_max(hotspot_ob,1,fix(k/2),1,fix(l/2));
[V_rt,H_rt]=find_max(hotspot_ob,1,fix(k/2),fix(l/2)+1,l);
[V_ld,H_ld]=find_max(hotspot_ob,fix(k/2)+1,k,1,fix(l/2));
[V_rd,H_rd]=find_max(hotspot_ob,fix(k/2)+1,k,fix(l/2)+1,l);
hotspot_lt = hotspot(1:fix(k/2),1:fix(l/2)); [k_lt, l_lt] = size(hotspot_lt);
hotspot_rt = hotspot(1:fix(k/2),(fix(l/2)+1):l); [k_rt, l_rt] = size(hotspot_rt);
hotspot_ld = hotspot((fix(k/2)+1):k,1:fix(l/2)); [k_ld, l_ld] = size(hotspot_ld);
hotspot_rd = hotspot((fix(k/2)+1):k,(fix(l/2)+1):l); [k_rd, l_rd] = size(hotspot_rd);
%V_lt=98; H_lt=150;
[m_lt,n_lt]=find_max(hotspot_lt,1,k_lt,1,l_lt);
[m_rt,n_rt]=find_max(hotspot_rt,1,k_rt,1,l_rt);
[m_ld,n_ld]=find_max(hotspot_ld,1,k_ld,1,l_ld);
[m_rd,n_rd]=find_max(hotspot_rd,1,k_rd,1,l_rd);
misfit = zeros(4,1);
misfit(1,1) = (hotspot_rt(m_rt,n_rt)-hotspot_rt(m_lt,n_lt)).^2+(hotspot_ld(m_ld,n_ld)-hotspot_ld(m_lt,n_lt)).^2 …
+(hotspot_rd(m_rd,n_rd)-hotspot_rd(m_lt,n_lt)).^2;
misfit(2,1) = (hotspot_lt(m_rt,n_rt)-hotspot_lt(m_rt,n_rt)).^2+(hotspot_ld(m_rt,n_rt)-hotspot_ld(m_rt,n_rt)).^2 …
+(hotspot_rd(m_rt,n_rt)-hotspot_rd(m_rt,n_rt)).^2;
misfit(3,1) = (hotspot_lt(m_ld,n_ld)-hotspot_lt(m_ld,n_ld)).^2+(hotspot_rt(m_ld,n_ld)-hotspot_rt(m_ld,n_ld)).^2 …
+(hotspot_rd(m_ld,n_ld)-hotspot_rd(m_ld,n_ld)).^2;
misfit(4,1) = (hotspot_lt(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2+(hotspot_rt(m_rd,n_rd)-hotspot_rt(m_rd,n_rd)).^2 …
+(hotspot_ld(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2;
switch min(misfit)
case misfit(1,1)
mm = m_lt; nn = n_lt;
case misfit(2,1)
mm = m_rt; nn = n_rt;
case misfit(3,1)
mm = m_ld; nn = n_ld;
case misfit(4,1)
mm = m_rd; nn = n_rd;
end
Max_color = [ hotspot_rt(mm,nn) hotspot_ld(mm,nn) hotspot_rd(mm,nn)]
m = min([k_lt;k_rt;k_ld;k_rd]); n = min([l_lt;l_rt;l_ld;l_rd]);
misfit2 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit2(countm,countn) = (hotspot_rt(countm,countn)-hotspot_rt(m_rt,n_rt)).^2 …
+(hotspot_ld(countm,countn)-hotspot_ld(m_ld,n_ld)).^2+(hotspot_rd(countm,countn)-hotspot_rd(m_rd,n_rd)).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit2(countm,countn) == min(min(misfit2))
mmm = countm; nnn = countn;
end
end
end
Max_color2 = [hotspot_rt(mmm,nnn) hotspot_ld(mmm,nnn) hotspot_rd(mmm,nnn)]
%V_lt=126; H_lt=176;
%V_rt=125; H_rt=565;
%V_ld=398; H_ld=179;
%V_rd=388; H_rd=580;
center_LT=[V_lt,H_lt];
center_LD=[V_ld, H_ld];
center_RT=[V_rt,H_rt];
center_RD=[V_rd,H_rd];
LT=hotspot(V_lt-hscale:V_lt+hscale,H_lt-hscale:H_lt+hscale);
LD=hotspot(V_ld-hscale:V_ld+hscale,H_ld-hscale:H_ld+hscale);
RT=hotspot(V_rt-hscale:V_rt+hscale,H_rt-hscale:H_rt+hscale);
RD=hotspot(V_rd-hscale:V_rd+hscale,H_rd-hscale:H_rd+hscale);
% Find New centers and update them
[center_LT_new,Cor_LT] = imagematching (hotspot,RD,center_RT,scale); % use RT as target image
[center_LD_new,Cor_LD] = imagematching (hotspot,RD,center_LD,scale);
[center_RT_new,Cor_RT] = imagematching (hotspot,RD,center_RD,scale);
center_LD=center_LD_new;
center_LT=center_LT_new;
center_RT=center_RT_new;
[inter_center_LT,Cor_LT] = intermatching (hotspot,RD,center_LT_new,scale,subpixel);
[inter_center_LD,Cor_LD] = intermatching (hotspot,RD,center_LD_new,scale,subpixel);
[inter_center_RT,Cor_RT] = intermatching (hotspot,RD,center_RT_new,scale,subpixel);
center_LD=inter_center_LD;
center_LT=inter_center_LT;
center_RT=inter_center_RT;
Cor=[1 Cor_LT;Cor_LD Cor_RT];
V_ld=center_LD(1,1); H_ld=center_LD(1,2);
V_rt=center_LT(1,1); H_lt=center_LT(1,2);
V_rd=center_RD(1,1); H_rd=center_RD(1,2);
[XI_ld,YI_ld]=meshgrid(H_ld-hscale:H_ld+hscale,V_ld-hscale:V_ld+hscale);
[XI_lt,YI_lt]=meshgrid(H_lt-hscale:H_lt+hscale,V_lt-hscale:V_lt+hscale);
[XI_rd,YI_rd]=meshgrid(H_rd-hscale:H_rd+hscale,V_rd-hscale:V_rd+hscale);
LD=interp2(hotspot,XI_ld,YI_ld,’spline’); % update images
LT=interp2(hotspot,XI_lt,YI_lt,’spline’);
RD=interp2(hotspot,XI_rd,YI_rd,’spline’);
% Assuming you want to exclude 580 nm
lamda_LT=580*10^(-9);
lamda_RT=640*10^(-9);
lamda_LD=766*10^(-9);
lamda_RD=905*10^(-9);
planck=6.6*10^(-34);
Kb=1.38*10^(-23);
c=3*10^8;
%Wien_LT=Kb/planck/c*log(2*pi*planck*c*c./LT/lamda_LT^5);
Wien_RT=Kb/planck/c*log(2*pi*planck*c*c./RT/lamda_RT^5);
Wien_LD=Kb/planck/c*log(2*pi*planck*c*c./LD/lamda_LD^5);
Wien_RD=Kb/planck/c*log(2*pi*planck*c*c./RD/lamda_RD^5);
[m,n]=size(LD);
for i=1:m
for j=1:n
X=[1/lamda_RT 1/lamda_LD 1/lamda_RD]’;
%X=[1/lamda_650 1/lamda_800 1/lamda_900]’;
X_reg=[X ones(size(X))];
%Y=[Wien_LT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
Y=[Wien_RT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
%R(i,j)=corrcoef(X,Y);
[T_inv,bint]=regress(Y,X_reg);
R(i,j)=corr(X,Y);
T(i,j)=1/T_inv(1,1);
Em(i,j)=exp(-T_inv(2,1)*planck*c/Kb);
T_err(i,j)=(1/bint(1,1)-1/bint(1,2))/2;
T_low(i,j)=1/bint(1,2);
T_high(i,j)=1/bint(1,1);

image result: So, I do research of laser heating metals. We have this code (main code) that separates each hotspot into different qudrants to analyze. Occacsionally, one quadrant nees to be excluded from the calculations (2nd code). So the main code, uses the LT (left top) as the focus for analyzing. The 2nd code, needs to do the analysis with LD (Left down) as the focus instead adn exclude everything that has to do with LT.
So the issue is: the tempeartures calculated from the 2nd code are 300 to 400 K too low. (pics of ffure with each code.). If you have any isight I would greatly appreciate it!
Here is the normal code:
clear;close all;
[filename, pathname, filterindex] = uigetfile(‘*.txt’, ‘Enter input data file..’);
fullname = fullfile(pathname, filename);
mydata = importdata(fullname);
hotspot_ob=double(mydata);
scale= 120; % size of the box 90*90
subpixel=4;
hscale=scale/2;
[k,l]=size(hotspot_ob);
SR=ones(size(hotspot_ob));
box_Num = 1; % box 1 or 2
ND_Num = 0; % 0, .5, 1, 2
switch ND_Num
case 0
ND_col = 1;
case .5
ND_col = 2;
case 1
ND_col = 3;
case 2
ND_col = 4;
end

Image result:

%LATEST from 08/09/2023
box(1).rs = [ 1.0 1.0 1.0 1.0
1.05 1.101 1.14 1.14
0.803 1.12 1.54 4.17
0.645 0.788 0.929 2.82];
box(2).rs = [1.0 1.0 1.0 1
0.296 0.312 0.325 0.424
0.331 0.457 0.632 1.735
0.245 0.294 0.348 1.077];
SR(1:fix(k/2),1:fix(l/2)) = box(box_Num).rs(1,ND_col);
SR(1:fix(k/2),fix(l/2)+1:l) = box(box_Num).rs(2,ND_col);
SR(fix(k/2)+1:k,1:fix(l/2)) = box(box_Num).rs(3,ND_col);
SR(fix(k/2)+1:k,fix(l/2)+1:l) = box(box_Num).rs(4,ND_col);
hotspot=hotspot_ob./SR;
[V_lt,H_lt]=find_max(hotspot_ob,1,fix(k/2),1,fix(l/2));
[V_rt,H_rt]=find_max(hotspot_ob,1,fix(k/2),fix(l/2)+1,l);
[V_ld,H_ld]=find_max(hotspot_ob,fix(k/2)+1,k,1,fix(l/2));
[V_rd,H_rd]=find_max(hotspot_ob,fix(k/2)+1,k,fix(l/2)+1,l);
hotspot_lt = hotspot(1:fix(k/2),1:fix(l/2)); [k_lt, l_lt] = size(hotspot_lt);
hotspot_rt = hotspot(1:fix(k/2),(fix(l/2)+1):l); [k_rt, l_rt] = size(hotspot_rt);
hotspot_ld = hotspot((fix(k/2)+1):k,1:fix(l/2)); [k_ld, l_ld] = size(hotspot_ld);
hotspot_rd = hotspot((fix(k/2)+1):k,(fix(l/2)+1):l); [k_rd, l_rd] = size(hotspot_rd);
%V_lt=98; H_lt=150;
[m_lt,n_lt]=find_max(hotspot_lt,1,k_lt,1,l_lt);
[m_rt,n_rt]=find_max(hotspot_rt,1,k_rt,1,l_rt);
[m_ld,n_ld]=find_max(hotspot_ld,1,k_ld,1,l_ld);
[m_rd,n_rd]=find_max(hotspot_rd,1,k_rd,1,l_rd);
misfit = zeros(4,1);
misfit(1,1) = (hotspot_rt(m_rt,n_rt)-hotspot_rt(m_lt,n_lt)).^2+(hotspot_ld(m_ld,n_ld)-hotspot_ld(m_lt,n_lt)).^2 …
+(hotspot_rd(m_rd,n_rd)-hotspot_rd(m_lt,n_lt)).^2;
misfit(2,1) = (hotspot_lt(m_rt,n_rt)-hotspot_lt(m_rt,n_rt)).^2+(hotspot_ld(m_rt,n_rt)-hotspot_ld(m_rt,n_rt)).^2 …
+(hotspot_rd(m_rt,n_rt)-hotspot_rd(m_rt,n_rt)).^2;
misfit(3,1) = (hotspot_lt(m_ld,n_ld)-hotspot_lt(m_ld,n_ld)).^2+(hotspot_rt(m_ld,n_ld)-hotspot_rt(m_ld,n_ld)).^2 …
+(hotspot_rd(m_ld,n_ld)-hotspot_rd(m_ld,n_ld)).^2;
misfit(4,1) = (hotspot_lt(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2+(hotspot_rt(m_rd,n_rd)-hotspot_rt(m_rd,n_rd)).^2 …
+(hotspot_ld(m_rd,n_rd)-hotspot_ld(m_rd,n_rd)).^2;
switch min(misfit)
case misfit(1,1)
mm = m_lt; nn = n_lt;
case misfit(2,1)
mm = m_rt; nn = n_rt;
case misfit(3,1)
mm = m_ld; nn = n_ld;
case misfit(4,1)
mm = m_rd; nn = n_rd;
end
Max_color = [ hotspot_rt(mm,nn) hotspot_ld(mm,nn) hotspot_rd(mm,nn)]
m = min([k_lt;k_rt;k_ld;k_rd]); n = min([l_lt;l_rt;l_ld;l_rd]);
misfit2 = zeros(m,n);
for countm = 1:m
for countn = 1:n
misfit2(countm,countn) = (hotspot_rt(countm,countn)-hotspot_rt(m_rt,n_rt)).^2 …
+(hotspot_ld(countm,countn)-hotspot_ld(m_ld,n_ld)).^2+(hotspot_rd(countm,countn)-hotspot_rd(m_rd,n_rd)).^2;
end
end
for countm = 1:m
for countn = 1:n
if misfit2(countm,countn) == min(min(misfit2))
mmm = countm; nnn = countn;
end
end
end
Max_color2 = [hotspot_rt(mmm,nnn) hotspot_ld(mmm,nnn) hotspot_rd(mmm,nnn)]
%V_lt=126; H_lt=176;
%V_rt=125; H_rt=565;
%V_ld=398; H_ld=179;
%V_rd=388; H_rd=580;
center_LT=[V_lt,H_lt];
center_LD=[V_ld, H_ld];
center_RT=[V_rt,H_rt];
center_RD=[V_rd,H_rd];
LT=hotspot(V_lt-hscale:V_lt+hscale,H_lt-hscale:H_lt+hscale);
LD=hotspot(V_ld-hscale:V_ld+hscale,H_ld-hscale:H_ld+hscale);
RT=hotspot(V_rt-hscale:V_rt+hscale,H_rt-hscale:H_rt+hscale);
RD=hotspot(V_rd-hscale:V_rd+hscale,H_rd-hscale:H_rd+hscale);
% Find New centers and update them
[center_LT_new,Cor_LT] = imagematching (hotspot,RD,center_RT,scale); % use RT as target image
[center_LD_new,Cor_LD] = imagematching (hotspot,RD,center_LD,scale);
[center_RT_new,Cor_RT] = imagematching (hotspot,RD,center_RD,scale);
center_LD=center_LD_new;
center_LT=center_LT_new;
center_RT=center_RT_new;
[inter_center_LT,Cor_LT] = intermatching (hotspot,RD,center_LT_new,scale,subpixel);
[inter_center_LD,Cor_LD] = intermatching (hotspot,RD,center_LD_new,scale,subpixel);
[inter_center_RT,Cor_RT] = intermatching (hotspot,RD,center_RT_new,scale,subpixel);
center_LD=inter_center_LD;
center_LT=inter_center_LT;
center_RT=inter_center_RT;
Cor=[1 Cor_LT;Cor_LD Cor_RT];
V_ld=center_LD(1,1); H_ld=center_LD(1,2);
V_rt=center_LT(1,1); H_lt=center_LT(1,2);
V_rd=center_RD(1,1); H_rd=center_RD(1,2);
[XI_ld,YI_ld]=meshgrid(H_ld-hscale:H_ld+hscale,V_ld-hscale:V_ld+hscale);
[XI_lt,YI_lt]=meshgrid(H_lt-hscale:H_lt+hscale,V_lt-hscale:V_lt+hscale);
[XI_rd,YI_rd]=meshgrid(H_rd-hscale:H_rd+hscale,V_rd-hscale:V_rd+hscale);
LD=interp2(hotspot,XI_ld,YI_ld,’spline’); % update images
LT=interp2(hotspot,XI_lt,YI_lt,’spline’);
RD=interp2(hotspot,XI_rd,YI_rd,’spline’);
% Assuming you want to exclude 580 nm
lamda_LT=580*10^(-9);
lamda_RT=640*10^(-9);
lamda_LD=766*10^(-9);
lamda_RD=905*10^(-9);
planck=6.6*10^(-34);
Kb=1.38*10^(-23);
c=3*10^8;
%Wien_LT=Kb/planck/c*log(2*pi*planck*c*c./LT/lamda_LT^5);
Wien_RT=Kb/planck/c*log(2*pi*planck*c*c./RT/lamda_RT^5);
Wien_LD=Kb/planck/c*log(2*pi*planck*c*c./LD/lamda_LD^5);
Wien_RD=Kb/planck/c*log(2*pi*planck*c*c./RD/lamda_RD^5);
[m,n]=size(LD);
for i=1:m
for j=1:n
X=[1/lamda_RT 1/lamda_LD 1/lamda_RD]’;
%X=[1/lamda_650 1/lamda_800 1/lamda_900]’;
X_reg=[X ones(size(X))];
%Y=[Wien_LT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
Y=[Wien_RT(i,j) Wien_LD(i,j) Wien_RD(i,j)]’;
%R(i,j)=corrcoef(X,Y);
[T_inv,bint]=regress(Y,X_reg);
R(i,j)=corr(X,Y);
T(i,j)=1/T_inv(1,1);
Em(i,j)=exp(-T_inv(2,1)*planck*c/Kb);
T_err(i,j)=(1/bint(1,1)-1/bint(1,2))/2;
T_low(i,j)=1/bint(1,2);
T_high(i,j)=1/bint(1,1);

image result: colormap, contour, image analysis MATLAB Answers — New Questions

Cart

Cart