function [C_fit,phi_fit] = original_shape_fit_Curvature(x,C_phi) % at points uniformly spaced by mean curvature
global ini_Val 

if nargin==1
    C_phi=zeros(ini_Val.ncontrol,1);
end



%Initialization
nn=ini_Val.ncontrol;
N=(nn-1)*40; % 120
gaussp=linspace(0,1,N);
%n_sample_fit=(nn-1)*20; % 100
Invariant_sample_bis = precal_generic(gaussp); 


[C,C1,C2,C3] = interpolate(Invariant_sample_bis,x);
%au = zeros(ngp_total,1);
%calculating the constants and surface properties at gauss point

au = (C1(:,1).^2 + C1(:,2).^2);
au = sqrt(au);
bu = -(C2(:,1).*C1(:,2))+(C2(:,2).*C1(:,1));
dau = ((C1(:,1).*C2(:,1))+(C1(:,2).*C2(:,2)))./au;
Hu = ((bu./(au.^2))+(C1(:,2)./C(:,1)))./(au);


tmp1 =  C1(:,1).*C2(:,1) + C1(:,2).*C2(:,2);
tmp2 = -C3(:,1).*C1(:,2) + C1(:,1).*C3(:,2);
tmp3 =  C2(:,2).* C(:,1) - C1(:,2).*C1(:,1);
Hpu = (tmp1.*(-3*bu./au.^2 - C1(:,2)./C(:,1)) ...
    + tmp2)./(au.^3) + tmp3 ./C(:,1).^2./au;



n_sample_fit = 0;

Hu_l = zeros(size(gaussp));
for k=3:length(au)-2
    Hu_l(k) = (abs(Hu(k)) + abs(Hu(k-1)))/2*Invariant_sample_bis.gweight;
end
Hu_avg = sum(Hu_l);

%Hu_avg=15;
n_sample_fit_l = zeros(length(gaussp)-1,1);
u_Hu_equi(1) = gaussp(1);
ctr =1;
divide = 2;
for k=2:length(au)
    Check_Hu = abs(Hu(k)) + abs(Hu(k-1))/2;
    if Check_Hu > Hu_avg & k>10 & k<length(au)-10
        %divide = max(2,round(Check_Hpu/Hpu_avg)); 
        n_sample_fit_l(k) = divide;
        for j=1:divide
            ctr=ctr+1;
            u_Hu_equi(ctr) = gaussp(k-1) + j/divide*(gaussp(k) - gaussp(k-1));
        end
    else
        n_sample_fit_l(k) = 1;
        ctr=ctr+1;
        u_Hu_equi(ctr) = gaussp(k);
    end
    n_sample_fit = n_sample_fit + n_sample_fit_l(k);
end





% 
% 
% u_max = 10;
% Del_u = zeros(n_sample_fit-1,1);
% for k=1:n_sample_fit-1
%     if abs(Hpu)> 1 & k > 3 & k < n_sample_fit-3       
%         Del_u(k) = 2/(abs(Hpu(k))+abs(Hpu(k+1)));
%     else
%         Del_u(k)=u_max;
%     end
% end
% 
% alpha = 1/sum(Del_u)
% pause(1);
% 
% 
% 
% u_Hpu_equi(1) = gaussp(1);
% for i=2:n_sample_fit-1
%     u_Hpu_equi(i) = u_Hpu_equi(i-1) + Del_u(i-1)*alpha;
% end
% u_Hpu_equi(n_sample_fit)=1;

figure(8)
plot(Hu(2:end),'ro')

figure(7)
plot(u_Hu_equi(2:end),'bx')


Invariant_u_Hu_equi = precal_generic(u_Hu_equi); 
[C_fit,C1,C2] = interpolate(Invariant_u_Hu_equi,x);
[phi_fit d_phi] = interpolate_func_2(Invariant_u_Hu_equi,C_phi);


% arclength=zeros(size(gaussp));
% arclength(1)=0;
% for i=2:length(au)
%     arclength(i)=arclength(i-1) + (au(i-1)+au(i))/2*Invariant_sample_bis.gweight;
% end
% 
% arcl_equi=linspace(0,arclength(end),n_sample_fit);
% 
% u_arcl_equi(1)=0;
% for i=2:n_sample_fit-1
%     [val, index]=min(abs(arclength-arcl_equi(i)));
%     if arclength(index)>arcl_equi(i)
%         index=index-1;
%     end
%     u_arcl_equi(i)=(arcl_equi(i)*(gaussp(index+1)-gaussp(index)) + ...
%         gaussp(index)*arclength(index+1) - gaussp(index+1)*arclength(index)) ...
%         / (arclength(index+1)-arclength(index));
%     
% end
% 
% %Original shape
% u_arcl_equi(n_sample_fit)=1;
% 
% Invariant_u_arcl_equi = precal_generic(u_arcl_equi); 
% [C_fit,C1,C2] = interpolate(Invariant_u_arcl_equi,x);
% [phi_fit d_phi] = interpolate_func_2(Invariant_u_arcl_equi,C_phi);
% 
% 
% figure
% plot(gaussp,au,'ko')
% figure
% plot(gaussp,arclength,'ro-')
% hold on
% for i=1:n_sample_fit
%  plot([u_arcl_equi(i) u_arcl_equi(i) ],[0 arcl_equi(i)],'gs-')
%  plot([0 u_arcl_equi(i) ],[arcl_equi(i) arcl_equi(i)],'gs-')
% end   