function [A,x,y,xm,ym,IxLim,IyLim]=readASC(FName,pth,xLim,yLim)
% [A,x,y,xm,ym,IxLim,IyLim]=readASC(FName,pth,xLim,yLim [,header])
% Read ASC files van NHI model betweem given coordinates
% Structuur
% ncols         1200
% nrows         1300
% xllcorner     0                sometimes xllcenter
% yllcorner     300000           sometimes yllcenter
% cellsize      250
% NODATA_value  -9999

if nargin<2 || isempty(pth), pth=[]; end

fprintf('Reading file    ''%s''\n',FName);

if ~isempty(pth) && pth(end)~=filesep, pth=[pth filesep]; end
fid=fopen([pth FName],'r');

fseek(fid,0,-1);
keyw=cell(1,6);
val =NaN(1,6);
for i=1:6
    keyw{i}=fscanf(fid,'%s',1);
    val(i) =fscanf(fid,'%f',1);
end
NCOL=val(strmatchi('ncols',keyw));
NROW=val(strmatchi('nrows',keyw));
XLL=val(strmatchi( 'xllcorner',keyw));
YLL=val(strmatchi( 'yllcorner',keyw));
CELLSZ=val(strmatchi('cellsize',keyw));
NOVAL=val(strmatchi('NODATA_value',keyw));

x=[XLL, XLL+CELLSZ*(1:NCOL)];                xm=0.5*(x(1:end-1)+x(2:end));
y=[YLL, YLL+CELLSZ*(1:NROW)]; y=y(end:-1:1); ym=0.5*(y(1:end-1)+y(2:end));

A=fscanf(fid,'%f',[NCOL,NROW])';
A(A==NOVAL)=NaN;

if nargin<4, xLim=x([1 end]); yLim=sort(y([1 end])); end

Ix=find(xm>xLim(1) & xm<xLim(2));
Iy=find(ym>yLim(1) & ym<yLim(2));

IxLim=Ix([1,end]);
IyLim=Iy([1,end]);

x=x(x>=xLim(1) & x<=xLim(2));
y=y(y>=yLim(1) & y<=yLim(2));

xm=xm(Ix);
ym=ym(Iy);

A=A(Iy,Ix);

fclose(fid);
