function [coeff,prob,rlo05,rhi95,a,b,c,d,e,f,g,h,i,j,avy1,avy2,avr,mdy1,mdy2,mdr] = coeff_ls(y1,y2)

%function [coeff,prob,rlo05,rhi95,a,b,c,d,e,f,avy1,avy2,avr,mdy1,mdy2,mdr] = coeff_ls(y1,y2)
%function [coeff,prob,rlo05,rhi95,a,b,c,avy1,avy2,avr,mdy1,mdy2,mdr] = coeff_ls(y1,y2)
%function [coeff,prob,rlo05,rhi95,a,b,c] = coeff_ls(ax,y1,y2)

% This function calculates the cross-correlation coefficient and least square solution of y = ax and y=bx+c
% y2 and y1 must be array of one column, and they should have the same dimension
% Have the transpose(') of all the vectors so they are vector, not row.

% Print out vectors
%y1
%y2

% y1 and y2 are vectors, and get correlation coeff between self and 2
% R(1,1) = cc between y1 and y1 = 1.0
% R(1,2) = cc between y1 and y2
% R(2,1) = cc between y2 and y1
% R(2,2) = cc between y2 and y2 = 1.0
% NOTE:  R(1,2)=R(2,1)

%[R] = corrcoef(y1,y2);
%coeff = R(1,2);
%prob = 0.;
%rlo05 = coeff;
%rhi95 = coeff;

% On cedar-l, the other corrcoef calls fail in getparams, strmatch, iscellstr, cellfun
%  with the error:  ??? Attempt to execute SCRIPT cellfun as a function.
%  This is with the probability, since it is returned as a zero from Sun-OS.
%[R,P] = corrcoef(y1,y2);
%prob = P(1,2)
% WARNING: DIVIDED BY ZERO IN NH 27 in corrcoef!
 [R,P,RLO,RHI] = corrcoef(y1,y2);
 coeff = R(1,2);
 prob = P(1,2);
 rlo05 = RLO(1,2);
 rhi95 = RHI(1,2);

% Added 12/07 to see lags of -2,-1,0,1,2 between y1 and y2
% Vsw best for 2,1,0,-1,-2
%clag = xcorr(y1,y2,2,'coeff')
% Option='coeff' means normalized cc so autocorr at 0 lag=1.0, as result of R(1,1) of corrcoeff
% c0=xcorr(y1,y2,'coeff') = cc in a vector of length 2N-1 (N=3468 hourly values in 1996 of HpeSHVBz>0)
%c0 = xcorr(y1,y2,'coeff');
% What is relation of clag (normalized) to corrcoef R??  (clag(middle)~0.80-0.89 ???)

% y2 = a * y1
% Find least-square solution
a = y1\y2;
%p1=polyfit(y1,y2,0);
% DIFFERENT ANSWERS and always the same (17.6022 for Ap, Kp and Vsw ???)
% BUT THIS should be a polynomial order of 1 with zero constant, so is OK!
%p1(1),a;

% Do point-by-point division
% NOTE:  ry1y2 is NOT y1/y2, but is actually y2/y1!!
ry1y2 = y1.\y2;
% Find ave ratio (y2 = avr * y1)
avr = mean(ry1y2);
% Find median ratio (y2 = mdr * y1)
mdr = median(ry1y2);
% Find aves and medians of y1 and y2
avy1=mean(y1);
avy2=mean(y2);
mdy1=median(y1);
mdy2=median(y2);

% y2 = b * y1 + c
x2 = [y1,ones(length(y1),1)];
temp =  x2\y2;
b = temp(1);
c = temp(2);
p2=polyfit(y1,y2,1);
% Same answers! b=p2(1), c=p2(2)=const
%p2(1),p2(2),b,c;

% y2 = d + e*y1 + f*y1^2
%A=[ones(length(y1), y1, y1^2)];
%temp2 = A\y2;
%d=temp2(1);
%e=temp2(2);
%f=temp2(3);
p3=polyfit(y1,y2,2);
f=p3(1);
e=p3(2);
d=p3(3);

% y3 = g + h*y1 + i*y1^2 + j*y1^3
p4=polyfit(y1,y2,3);
j=p4(1);
i=p4(2);
h=p4(3);
g=p4(4);
