diff --git a/qem/FitLocalErrors/FitLocalErrors.m b/qem/FitLocalErrors/FitLocalErrors.m
index 712897a7b8b887578690ad821905d2b59847677b..f86943f5d4d15800551b15dcd38fe036387a6e4a 100644
--- a/qem/FitLocalErrors/FitLocalErrors.m
+++ b/qem/FitLocalErrors/FitLocalErrors.m
@@ -2,6 +2,7 @@ function [knl, ksl, serr] = FitLocalErrors(....
     ReferenceRM,...
     PerturbedRM,...
     errorlocation,...
+    errorlength,...
     neig)
 %FITLOCALERRORS fit the best error to reproduced the difference between 
 % Reference and Perturebed RM at errorlocation s positions
@@ -10,6 +11,7 @@ function [knl, ksl, serr] = FitLocalErrors(....
 % ReferenceRM   : folder where the RM data is loaded or RM itself
 % PerturbedRM   : folder where the RM data is loaded or RM itself
 % errorlocation : list of s positions of error sources (
+% errorlength   : length of the error source (default 1e-6)
 % neig          : # of eigenvectors for svd (default: length(errorlocation))
 %
 % OUTPUT:
@@ -23,11 +25,18 @@ function [knl, ksl, serr] = FitLocalErrors(....
 % WARNING: 
 % errorlocation should not be inside dipoles! an
 % 
-%see also: qemextract atsbreak
+%see also: qemextract atsbreak atinsertelems
 if nargin<4
+    errorlength = 1e-6 * ones(size(errorlocation)); 
+end
+if nargin<5
     neig=length(errorlocation); 
 end
 
+if any(errorlength<0)
+    error('Error source lengths must be positive');
+end
+
 % load reference (normal and skew)
 [qemb, qemres, semres] = qemextract(ReferenceRM);
 
@@ -44,18 +53,77 @@ if any(errorlocation>C) | any(errorlocation<0)
     error(['Some locations are below 0 or above length of lattice (' num2str(C) ')']);
 end
 
-[rerr,errpos]=atsbreak(r,errorlocation);
-Lthin = 1e-6;
-
-rerr = atsetfieldvalues(rerr,errpos,'PassMethod', 'StrMPoleSymplectic4Pass');%'ThinMPolePass';
-rerr = atsetfieldvalues(rerr,errpos,'Class', 'Quadrupole');% NECESSARY for FIT function
-rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,1}, 0);
-rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,1},0);
-rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,2}, 0);
-rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,2},0);
-rerr = atsetfieldvalues(rerr,errpos,'MaxOrder', 1);
-rerr = atsetfieldvalues(rerr,errpos,'Length', Lthin);
-rerr = atsetfieldvalues(rerr,errpos,'NumIntSteps', 1);
+Lthin = errorlength;
+if  ~any(Lthin ~= 1e-6) % if all thin, keep existing behaviour 
+    [rerr,errpos]=atsbreak(r,errorlocation);
+    
+    rerr = atsetfieldvalues(rerr,errpos,'PassMethod', 'StrMPoleSymplectic4Pass');%'ThinMPolePass';
+    rerr = atsetfieldvalues(rerr,errpos,'Class', 'Quadrupole');% NECESSARY for FIT function
+    rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,1}, 0);
+    rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,1},0);
+    rerr = atsetfieldvalues(rerr,errpos,'PolynomB',{1,2}, 0);
+    rerr = atsetfieldvalues(rerr,errpos,'PolynomA',{1,2},0);
+    rerr = atsetfieldvalues(rerr,errpos,'MaxOrder', 1);
+    rerr = atsetfieldvalues(rerr,errpos','Length', Lthin');
+    rerr = atsetfieldvalues(rerr,errpos,'NumIntSteps', 1);
+else % long errors
+    
+    rerr =r; 
+    
+    % loop error locations
+    for ierr = 1:length(errorlocation)
+        [~,ipos]=atsbreak(rerr,errorlocation(ierr)); % find element where to insert 
+    
+        elem=atmultipole(['Error' num2str(ierr)],...
+            errorlength(ierr),[0 0],[0 0],...
+            'StrMPoleSymplectic4Pass',...
+            'Class','Quadrupole',...
+            'MaxOrder', 1,...
+            'NumIntSteps', 20);
+        
+        %find fraction to insert
+        Lins = rerr{ipos-1}.Length; % length of element where to insert the new element.
+        if Lins < errorlength(ierr)
+            error(['Not enough space at location :' num2str(errorlength(ierr)) ' m']);
+        end
+        
+        PMins = rerr{ipos-1}.PassMethod; % Passmethod of element where to insert the new element.
+        if ~strcmp(PMins, 'DriftPass')
+            disp(rerr{ipos-1});
+            error(['Not a DriftPass at location :' num2str(errorlength(ierr)) ' m']);
+        end
+        
+        
+        sp = findspos(rerr,ipos-1);
+        frac = (errorlocation(ierr)-sp)./Lins;
+        rerr = atinsertelems(rerr,ipos-1,frac,elem);
+    end
+    
+    % get error location indexes
+    errpos= find(atgetcells(rerr,'FamName','Error\w*'));
+       
+    % check that no negative drift is created.
+    poscheck = sort([errpos;errpos-1;errpos+1]);
+    DriftLengths = atgetfieldvalues(rerr,poscheck,'Length');
+    sposerr = findspos(rerr,poscheck);
+    if any(DriftLengths<0)
+        for ip = 1:length(poscheck)
+            disp([rerr{poscheck(ip)}.FamName '@' num2str(sposerr(ip),'%3.2f')...
+                ' : ' num2str(DriftLengths(ip),'%4.3f') ' m']);
+        end
+        error('Found Negative Drifts');
+    end
+    
+    % check lattice length unchanged.
+    C0 = findspos(r,length(r)+1);
+    Cnew = findspos(rerr,length(rerr)+1);
+    if Cnew-C0 > 1e6*length(errorlocation) 
+        disp(C0);
+        disp('became: ');
+        disp(Cnew);
+        error('Lattice length changed');
+    end
+end
 
 % redefine fields of qemres and semres for fit of rm difference with local error
 qemres.qpidx = errpos; % index of location to fit