function [r_ab,R,PLOT] =Circle_Para_LS(rou_j,theta_1_j,Plot_switch)
% --------------------------------------------------------------------------------------------------------------------------------------------------------
% A least-square fit for determining parameters of the circle equation (x?a)^2 + (y?b)^2 = R^2 using observations r_j:=[x_j, y_j]:=[?_j*cos(?_1_j), ?_j*sin(?_1_j)].
% --------------------------------------------------------------------------------------------------------------------------------------------------------
%-----------------------
% Syntax
%------------------------
% [d,AOOp,Center,Circ,R] =Circle_Para_LS(rou_j,theta_1_j) returns the downward FAC density
% [d,AOOp,Center,Circ,R] =Circle_Para_LS(...,'OptionalVariableName',OptionalVariableValue,...)
% MFACE_v1(...,'OptionalVariableName',OptionalVariableValue,...)
%OptionalVariableName is
%OptionalVariableName list:
%------------------------
% Inputs
%------------------------
%rou_j: ?_j(t), denoting the distance between O and the ring tip
%theta_1_j:  ?_1_j, when length(rou_j)==length(theta_1_j),
%                   ?_{j-1}_j, when length(rou_j)==length(theta_1_j)+1,
%                   (Here, ?_i_j denotes the cross angle between the  i-th and j-th cores)
%     Plot_switch: a key to contrle plot,  1 for plotting the LS results and 0 for not
%------------------------
% Outputs
%------------------------
%r_ab: complex coordinate of the LS-fitted circle center.  :
%(the orientation of  circle center:angle_r_ab=angle(r_ab)/pi*180;the distance between  r_ab and O: d=abs(r_ab):  )
%R: LS-fitted radius
% PLOT: a structure comprising the handles and data of the figure.
%------------------------
%Examples
%------------------------
%---% 1st example: LS fit to 3 points 
%rou_j=rand(3,1);theta_1_j=rand(3,1)*180;[r_ab,R,PLOT] =Circle_Para_LS(rou_j,theta_1_j,1)
%---% 2nd example: LS fit to 4 points, 
%rou_j=rand(4,1);theta_1_j=rand(3,1)*180;
%[r_ab,R,PLOT] =Circle_Para_LS(rou_j,theta_1_j,1)
%---% 3rd example
% rou_j=rand(3,1);theta_1_j=rand(3,1)*180;[r_ab,R,PLOT] =Circle_Para_LS(rou_j,theta_1_j,1);
% close all
% Inputs_complex=PLOT.outputCircData*5./(5+rand(size(PLOT.outputCircData)))+3;
% Inputs_complex=Inputs_complex'*Inputs_complex(1)
% [r_ab,R,PLOT2] =Circle_Para_LS(abs(Inputs_complex),angle(Inputs_complex(2:end))/pi*180,1);
% 
%%
% Copyright [2012] [Maosheng He and Keyan Fang]
% Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%     http://www.apache.org/licenses/LICENSE-2.0
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and limitations under the License.
%%
% Cite As
%Fang et al., (2021). Variations in tree-ring geometric center: implications for wind speed change and basal area increment ...
%------------------------
%Communication to:
%------------------------
% Maosheng He
% hmq512@gmail.com
% 20-04-2021
%%
if nargin<3,Plot_switch=0;end
rou_j=rou_j(:);
theta_1_j=theta_1_j(:);
theta_1_j=mod(theta_1_j,360);
% transfer ?_j-1_j to ?_1_j
if length(theta_1_j)< length(rou_j)
    theta_1_j=cumsum([0;theta_1_j(1:length(rou_j)-1)/360*2*pi]);
end
rou_j_comp=rou_j(:).*exp(1i*theta_1_j);%?_j_complex:=?_j*exp(1i*?_1_j);
rou_j_comp=rou_j_comp(:).';
%% LS fit to x^2+y^2+A*x+ B*y+C=0 derived  from (x?a)^2 + (y?b)^2 = R^2
xy1 = [real(rou_j_comp); imag(rou_j_comp); ones(size(rou_j_comp))];
ABC = [xy1']\abs(rou_j_comp)'.^2; % coefficents A, B, C and
r_ab=ABC(1)/2+ABC(2)/2*1i; % coordinates of the circle center
R = sqrt(abs(r_ab)^2+ABC(3)); % radius

%% plot
if Plot_switch,
    theta_i=linspace(0,2*pi,1000);
    PLOT.outputCircData=r_ab+[R*exp(1i*theta_i)];
    PLOT.gcf=figure;
    PLOT.gca=subplot(1,1,1);
    if length(rou_j_comp(:))<5
    PLOT.input=compass(rou_j_comp,'b');
    else
    PLOT.input=scatter(real(rou_j_comp),imag(rou_j_comp),'bx');
    end
    hold on;
    PLOT.output_circle=plot(real(PLOT.outputCircData),imag(PLOT.outputCircData),'r-');
    hold on; PLOT.output_center=scatter(real(r_ab),imag(r_ab),'r+');
    axis equal
else
    PLOT=nan;
end