function errmap(fconf)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function:  errmap.m
% Package:   HoP
% Author:    Franco Buffa
% Purpose:   
% Date:      
%	
% GNU Octave Version: 4.0.2
%
% revised by 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% errmap.m  is a function that 
%
% Output:
% 
%
% Input:
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 if nargin>0
  run(fconf); 
 else
  run('conf.m');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 load unwr.dat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% u0=sfactor*u0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 dx=D/(mp-1)/k0;

 W=reshape(u0,mp*mp,1);
 X1=ones(mp,mp).*(1:mp); % colonna
 Y1=transpose(X1);       % riga
 
 x=dx*(X1-mp/2);
 y=dx*(Y1-mp/2);
 rho2=x.^2+y.^2;
 delta2=x.^2-y.^2;
 
 X=reshape(x,mp*mp,1);
 Y=reshape(y,mp*mp,1);
 RHO2=reshape(rho2,mp*mp,1);
 DELTA2=reshape(delta2,mp*mp,1);
 
 X(blk)=[];
 Y(blk)=[];
 RHO2(blk)=[];
 DELTA2(blk)=[];
 W(blk)=[];
 if (~exist('aber','var'))
  aber=0;
 end
 if(aber)
  switch aber
   case 1
    a=[ones(length(X),1) X Y RHO2 RHO2.*X RHO2.*Y];
    c=a\W;
    u0=u0-(c(1)+x*c(2)+y*c(3)+rho2*c(4)+rho2.*x*c(5)+rho2.*y*c(6));
   case 2 
    a=[ones(length(X),1) X Y RHO2 RHO2.*X RHO2.*Y DELTA2];
    c=a\W;
    u0=u0-(c(1)+x*c(2)+y*c(3)+rho2*c(4)+rho2.*x*c(5)+rho2.*y*c(6)+delta2*c(7));   
   case 3 
    a=[ones(length(X),1) X Y DELTA2];
    c=a\W;
    u0=u0-(c(1)+x*c(2)+y*c(3)+delta2*c(4));   
   otherwise
    error("invalid aber value");
  endswitch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if(aber<3)
   printf('********************************\n');
   xf=c(5)*2*F^3*lambda/pi;
   yf=c(6)*2*F^3*lambda/pi;
   zf=c(4)*2*F^2*lambda/pi;
%  phi=atan2((c(3)*lambda/(2*pi)-yf/F),(c(2)*lambda/(2*pi)-xf/F))*180/pi  
%  theta=((c(3)*lambda/(2*pi)-yf/F)/sind(phi))*180/pi
   printf('Feed displacements:\n');
   printf('x=%f (mm)\n',1E3*xf);
   printf('y=%f (mm)\n',1E3*yf);
   printf('z=%f (mm)\n',1E3*zf);
   printf('********************************\n');
 end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 else
  a=[ones(length(X),1) X Y];
  c=a\W;
  u0=u0-(c(1)+x*c(2)+y*c(3));
 end

 u0(blk)=0; 
 u1=lambda/4/pi*sqrt(1.0-(x.^2+y.^2)/4/(F^2)).*u0;
% u1=lambda/4/pi*u0;
 u1=u1*1E3;
%i3=find(abs(u1)>tol3);
%u1(i3)=0;
 
 u1=sfactor*u1;

 save errmap.dat u1
 
 switch(panels)
  case(1)
   fname=[datapath,'jacks.dat'];
   load(fname);
   R=D/2;
   jx=round((jx)*r1/R)+mp/2;
   jy=round((jy)*r1/R)+mp/2;
   jy=mp+1-jy;
   for k=1:length(jx)
    u1(jx(k),jy(k))=tol3;
   end  
  case(2) 
   fname=[datapath,'panel.dat'];
   load(fname);
   p1=map0.*(map0==1);
   ip1=find(p1==0);
   u1(ip1)=tol3;
   u1(blk)=0;
 end 
 
 figure(1);
 imagesc(dx*(-mp/2:mp/2),dx*(-mp/2:mp/2),u1);
 caxis([-tol3 tol3]);
%grid minor on;
 axis('square');
 colormap(cmap);
 colorbar;
 title('error map (mm)');
%print -dpng errmap.png
%   figure(9);
%   surf(u1);
%   colormap('jet');
%   view(145,70); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
end