function [C_fit,phi_fit] = original_shape_fit_Curvature2(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)*20; % 120
gaussp=linspace(0,1,N);
n_sample_fit=(nn-1)*30; % 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;

% Hpu(1)=0;
% Hpu(2)=0;
% Del_u=gaussp(2)-gaussp(1);
% for i=3:length(au)-1
%     Hpu(i) = (Hu(i+1)-Hu(i-1))/Del_u/2;
% end
% Hpu(end)=0;
Sq_tmp = sqrt(1+Hpu.^2);

tmp_Hu = au.*abs(C1(:,2));

arclength_H=zeros(size(gaussp));
arclength_H(1)=0;
for i=2:length(au)
    arclength_H(i)=arclength_H(i-1) + (tmp_Hu(i)+tmp_Hu(i-1))/2*Invariant_sample_bis.gweight;
end
% arclength_H(end) = 2*arclength_H(end-1) - arclength_H(end-2);

arcl_equi_H=linspace(0,arclength_H(end),n_sample_fit);


u_arcl_equi_H(1)=0;


% for i=2:10
%     u_arcl_equi_H(i) = (i-1)/n_sample_fit;
% end


%u_min=0;
u_min=0.96/n_sample_fit;
u_max=1.04/n_sample_fit;
for i=2:n_sample_fit-1
    [val, index]=min(abs(arclength_H-arcl_equi_H(i)));
    if arclength_H(index)>arcl_equi_H(i)
        index=index-1;
    end
    u_arcl_equi_H(i)=(arcl_equi_H(i)*(gaussp(index+1)-gaussp(index)) + ...
        gaussp(index)*arclength_H(index+1) - gaussp(index+1)*arclength_H(index)) ...
        / (arclength_H(index+1)-arclength_H(index));
    if (u_arcl_equi_H(i) - u_arcl_equi_H(i-1))<u_min
        u_arcl_equi_H(i) = u_arcl_equi_H(i-1) + u_min;
    end
    if (u_arcl_equi_H(i) - u_arcl_equi_H(i-1))>u_max
        u_arcl_equi_H(i) = u_arcl_equi_H(i-1) + u_max;
    end
end
% 
% for i=n_sample_fit-10:n_sample_fit-1
%     u_arcl_equi_H(i) = (i)/n_sample_fit;
% end

%Original shape
u_arcl_equi_H(n_sample_fit)=1;



figure(7)
plot(diff(u_arcl_equi_H),'bx')
hold on



figure(8)
plot(arclength_H,'bx')
hold on


arclength_S=zeros(size(gaussp));
arclength_S(1)=0;
for i=2:length(au)
    arclength_S(i)=arclength_S(i-1) + (au(i-1)+au(i))/2*Invariant_sample_bis.gweight;
end

arcl_equi_S=linspace(0,arclength_S(end),n_sample_fit);

u_arcl_equi_S(1)=0;
for i=2:n_sample_fit-1
    [val, index]=min(abs(arclength_S-arcl_equi_S(i)));
    if arclength_S(index)>arcl_equi_S(i)
        index=index-1;
    end
    u_arcl_equi_S(i)=(arcl_equi_S(i)*(gaussp(index+1)-gaussp(index)) + ...
        gaussp(index)*arclength_S(index+1) - gaussp(index+1)*arclength_S(index)) ...
        / (arclength_S(index+1)-arclength_S(index));
    
end

%Original shape
u_arcl_equi_S(n_sample_fit)=1;

figure(7)
plot(diff(u_arcl_equi_S),'ro')
hold off

figure(8)
plot(arclength_S,'ro')
hold on

u_arcl_equi = (u_arcl_equi_H);
%u_arcl_equi = (u_arcl_equi_H)/2 + (u_arcl_equi_S)/2;


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   