function run_williams_curved(id,zone,cparam,filres)

check=cparam.check;
Step=cparam.steps;
Rmax=cparam.radius;% taille de la zone d'extraction
mask_radius=cparam.mask_radius;% taille de la zone a exclure au niveau de la pointe
mask_width=cparam.mask_width;% taille de la zone a exclure le long des faces de la fissure
km_indices=-3:7;
modes=[1,2];
mu=cparam.mu;% module de cisaillement
kappa=cparam.kappa;% kolossov constant
%mu=0.5*E/(1+nu);
%kappa=(3-4*nu);% kolossov constant en deformation plane (echantillon massif)
pix2m=cparam.pix2m;% correspondance pixel/metre
maxiter=500;
filk=sprintf('%s-crack-%02d-sif.res',strrep(filres,'.res',''),id);
load(filres,'-mat','U','xo','yo','param','Nnodes','Nelems','model','conn','elt');
etype=min(elt);
model.zone=zone;
U=U(:,Step);
U=U(1:prod(Nnodes),:)+1i*U(prod(Nnodes)+(1:prod(Nnodes)),:);
sizeim=[param.roi(2)-param.roi(1),param.roi(4)-param.roi(3)]+1;
xoo=xo;yoo=yo;
zo=xo+param.roi(1)-1+1i*(yo+param.roi(3)-1);
conno=conn;
elto=elt;
conn=cell(1,numel(Step));
elt=cell(1,numel(Step));
xo=cell(1,numel(Step));
yo=cell(1,numel(Step));
Ut=cell(numel(Step),1);
Uw=cell(numel(Step),1);
Urbt=zeros(1,numel(Step));
checka=zeros(numel(Step),1);

xyc=zone{2};
zcini=xyc*[1;1i];
tend=xyc(end,:)-xyc(end-1,:);
xyc(end,:)=xyc(end,:)+1e0*tend;
a=0;
scal=2*mu*sqrt(2*pi);
scalamp=scal*(pix2m.^(1-(km_indices)*.5));

Uo=U;
Iterationfin=ones(1,size(U,2));
xytips=ones(size(U,2),1);
Tfin=ones(1,size(U,2));
Bfin=ones(1,size(U,2));
K1fin=ones(1,size(U,2));
K2fin=zeros(1,size(U,2));
Afin=ones(1,size(U,2));
Widthfin=ones(1,size(U,2));

xc=xyc(1,1);yc=xyc(1,2);
xf=xyc(2,1);yf=xyc(2,2);
zci=(xc+1i*yc);
zf=(xf+1i*yf);
tt=zci-zf;
tt=tt/abs(tt);
nn=tt*exp(1i*pi/2);
thetap=angle(tt);
xcp=xc;
ycp=yc;
zc=xyc*[1;1i];

dcum=abs(diff(zc-zc(1)));
if check
    hf= figure;
end
decc=0;
for iim=size(U,2):-1:1
    disp(sprintf('Step %2d / %2d',Step(iim),Step(size(U,2))));
    iteration=0;
    decx=0;dda=Inf;
    while ((iteration<maxiter) && (((abs(decx)>1)&&(abs(dda)>0.1))||iteration<2))
        iteration=iteration+1;
%    if Step(iim)==14
%       keyboard
%    end
        if abs(decx)>max(sizeim)
            if iteration==1
                dda=0;
            else
                dda=-1.9*abs(dda)*sign(decx);
            end
        elseif abs(decx)>max(model.mesh_size)
            dda=2*mean(model.mesh_size)*sign(decx);
        else
            dda=decx;
        end
        dda=0.5/(1+(iteration>40))/(1+(iteration>50))/(1+(iteration>100))/(1+(iteration>200))/(1+(iteration>20))/(1+(iteration>30))*dda;
        decc=decc+dda;
        if (-(decc)>dcum(1))&(size(xyc,1)>2)
%             keyboard
            a=a-(decc-dda+dcum(1));
            xyc=xyc(2:end,:);
            zc=xyc*[1;1i];
            dcum=abs(diff(zc-zc(1)));
            decc=0;
        else
            a=a+dda;
            xyc(1,1)=xyc(1,1)+dda*real(tt);
            xyc(1,2)=xyc(1,2)+dda*imag(tt);
            zc=xyc*[1;1i];
        end
        
        xc=xyc(1,1);yc=xyc(1,2);
        xf=xyc(2,1);yf=xyc(2,2);
        zci=(xc+1i*yc);
        zf=(xf+1i*yf);
        tt=zci-zf;
        tt=tt/abs(tt);
        nn=tt*exp(1i*pi/2);
        thetap=angle(tt);
        
        Uef=U(:,iim)*exp(-1i*thetap);
        
        
        
        [crackc,~,~,~]=GetSignedDistanceToCrack(xyc,zo,0);
        front=real((zo-zci)*conj(tt));
        crack=real((zo-zci)*conj(nn));
        
%         figure
%         plot(zo,'bx')
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')
%         
%         figure
%         scatter(real(zo),imag(zo),[],front)
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')
% 
%                 figure
%         scatter(real(zo),imag(zo),[],crack)
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')
% 
%                 figure
%         scatter(real(zo),imag(zo),[],crackc)
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')
%                 figure
%         scatter(real(zo),imag(zo),[],sign(crackc.*crack))
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')

        
        %            front=front-0.5*decx;
        z=front+1i*crack;
        dist=max(abs(z),1);
        angl=angle(z);
%                     hf1=figure;
%                     scatter(real(zo),imag(zo),[],front)
%                     hold on
%                     plot(zcini,'r-')
%                     plot(zc(1:end-1),'ro')
%                     hf2=figure;
%                     scatter(real(zo),imag(zo),[],crack)
%                     hold on
%                     plot(zcini,'r-')
%                     plot(zc(1:end-1),'ro')
%                      pause
%                     close(hf1)
%                     close(hf2)
        mask=(~((dist>Rmax)|(dist<mask_radius)|((abs(crack)<mask_width)&(front<0))));
        maskc=(~((dist>Rmax)|(dist<mask_radius)|((abs(crackc)<mask_width)&(front<0))));
        mask=double(mask&maskc);
%                  figure
%         scatter(real(zo),imag(zo),[],double(mask))
%         hold on
%         axis equal
%         plot(zf,'k+')
%         plot(zci,'ro')
%         plot(xyc(:,1),xyc(:,2),'r-')
%         quiver(real(zci), imag(zci), real(nn), imag(nn),100,'Color','k')
%         quiver(real(zci), imag(zci), real(tt), imag(tt),100,'Color','r')
       
        mask=diag(sparse(mask(:)));
        
        icon=0;
        phic=zeros(length(dist),length(modes)*length(km_indices));
        for m=1:length(modes)
            for ii=1:length(km_indices)
                icon=icon+1;
                meq=modes(m);
                heq=km_indices(ii);
                [harm1,amp1]=KMPotFourrier(meq,heq,kappa);
                feq=0*dist(:);
                for kk=1:3
                    feq=feq+amp1(kk)*exp(harm1(kk)*angl(:)*1i/2);
                end
                feq=feq.*(dist(:).^(heq/2));
                %               feq=feq.*exp(1i*thetap);
                
                phic(:,icon)=feq(:);
            end
        end
        
        M=real(phic'*mask*phic);
        F=real(phic'*mask*Uef);
        [LT,UT]=lu(sparse(M));
        Ukm=UT\(LT\F);
        %        Ukm=M\F;
        found=find(km_indices==1);
        K1=scalamp(found)*Ukm(found);
        K2=scalamp(found)*Ukm(found+length(km_indices));
        found=find(km_indices==-1);
        SK1=scalamp(found)*Ukm(found);
        found=find(km_indices==-3);
        SK3=scalamp(found)*Ukm(found);
        decx=-2*SK1/K1/pix2m;
        width=(sqrt(-8*SK3/K1))/pix2m;
        found=find(km_indices==2);
        %        T=[scalamp(found)*Ks(found,:)];
        T=[scalamp(found)*Ukm(found)]*4/sqrt(2*pi);
        found=find(km_indices==3);
        B=[scalamp(found)*Ukm(found)];
        if dda==0&&iteration>1
            decx=0;
        end
        disp(sprintf('At iteration %2d decx %f applied dda %f',iteration,decx,dda));
        
        if check%&&~(((iteration<maxiter) && (abs(decx)>1||iteration==0)))
            Xi=real(zo*exp(-1i*thetap));
            Yi=imag(zo*exp(-1i*thetap));
            maskn=diag(mask);
            maskn(maskn==0)=NaN;
            set(0,'CurrentFigure',hf);
            subplot(2,2,1)
            scatter(Xi,Yi,[],imag(maskn.*Uef))
            colorbar
            axis equal
            subplot(2,2,2)
            scatter(Xi,Yi,[],imag(maskn.*(phic*Ukm)));
            colorbar
            axis equal
            subplot(2,2,3:4)
            scatter(Xi,Yi,[],abs(imag(maskn.*Uef)-imag(maskn.*(phic*Ukm))))
            colorbar
            axis equal
            pause%(1)
        end
        
        
    end
   
    decc=decc+(abs(decx)<1)*decx;
    a=a+(abs(decx)<1)*decx;
            xyc(1,1)=xyc(1,1)+(abs(decx)<1)*decx*real(tt);
            xyc(1,2)=xyc(1,2)+(abs(decx)<1)*decx*imag(tt);

    Tfin(iim)=T;
    Bfin(iim)=B;
    K1fin(iim)=K1;
    if any(modes==2)
        K2fin(iim)=K2;
    end
    Afin(iim)=a*pix2m;
    Iterationfin(iim)=iteration;
    Widthfin(iim)=width*pix2m;
    xytips(iim)=xyc(1,:)*[1;1i];
    
    
    checka(iim)=decx;
    nodes=diag(mask)>0;
    [eltn,~]=GetEltsFromNodes(conno,elto,nodes,0);
    [pind,~,j1]=unique(conno(eltn,1:etype));
    conn{iim}=reshape(j1,[numel(eltn),etype]);
    xo{iim}=xoo(pind);
    yo{iim}=yoo(pind);
    elt{iim}=elto(eltn);
    Ut{iim}=[real(Uo(pind,iim));imag(Uo(pind,iim))];
    wil=(phic*Ukm)*exp(1i*thetap);
    Uw{iim}=[real(wil(pind));imag(wil(pind))];
            found=find(km_indices==0);
            Urbt(iim)=exp(1i*thetap)*(phic(1,found)*Ukm(found)+phic(1,found+length(km_indices))*Ukm(found+length(km_indices)));
end
K1=K1fin;
K2=K2fin;
T=Tfin;
B=Bfin;
da=Afin;
rho=Widthfin;
try
if numel(Step)>1
    param.deformed_image=param.deformed_image(:,Step);
else
    if iscell(param.deformed_image)
        param.deformed_image=param.deformed_image{:,Step};
    end
end
catch
end
param.result_file=filk;
param.detect=1;
rint=0;
U=Uw;
if isfield(model,'visu')
    model=rmfield(model,'visu');
end
save(filk,'Urbt','K1','K2','T','B','checka','xytips','Widthfin','da','rho','xo','yo','conn','U','Ut','Nnodes','Nelems','rint','elt','model','param');
%         param.result_file=strrep(filk,'.res','-w.res');
%         U=Uw;
%  save(param.result_file,'K1','K2','T','B','xytips','Widthfin','da','xo','yo','conn','U','Nnodes','Nelems','rint','elt','model','param');
%         param.result_file=strrep(filk,'.res','-dfit.res');
%         for iim=1:size(U)
%             U{iim}=Ut{iim}-Uw{iim};
%         end
%  save(param.result_file,'K1','K2','T','B','xytips','Widthfin','da','xo','yo','conn','U','Nnodes','Nelems','rint','elt','model','param');

end