function [a,b1, b2, R] =two_var_lreg(t1,t2,y);

% t1 = m2c(TS);
% t2 = m2c(TV);
% y = m2c(OLR);



%linear regression for dependent variable y and independent variable t
% a is intercept and b slope R is the R^2 value

if min(size(t1)) > 1 |min(size(t2)) > 1 |  min(size(y)) > 1
    error(' This program is designed for 1 D data')
end

if (length(t1) ~= length(t2)) | (length(t1) ~= length(y));
    error('All inputs must be the same length')
end


XX=[ones(size(t1)) t1 t2];

A=XX \ y;
a=A(1,1) 
b1=A(2,1) 
b2=A(3,1)

%we should be able to verify the above results with theory
B1 = (std(y)/std(t1))*(corr(t1,y)-corr(t1,t2)*corr(t2,y))/(1-corr(t1,t2).^2);
B2 = (std(y)/std(t2))*(corr(t2,y)-corr(t1,t2)*corr(t1,y))/(1-corr(t1,t2).^2);

y_out = a+ b1*t1 + b2*t2; 
MSE = mean((y_out-y).^2); %mean squared error 
R = (var(y)-MSE)/var(y);

y_out_min=a+b1*min(t1) + b2*min(t2);
y_out_max=a+b1*max(t1) + b2*max(t2);

% figure; clf;
% plot3(t1,t2, y,'r*', 'MarkerSize', 12)
% hold on; plot3([min(t1) max(t1)],[min(t2) max(t2)],[y_out_min y_out_max] , 'b', 'LineWidth', 2);
% title(['R squared is ', num2str(R)], 'FontSize', 17)
% set(gca, 'FontSize', 18);
