function [X,T,VAR]=oct_get_xt(hisfile,gridfile,lonsec,latsec,vname,tindices,vlevel);
% Copy of get_xt using oct_ versions of helper functions

if nargin < 1
  error('You must specify a file name')
end
if nargin < 2
  gridfile=hisfile;
end
if nargin < 3
  lonsec=[12 18];
end
if nargin < 4
  latsec=-34;
end
if nargin < 5
  vname='temp';
end
if nargin < 6
  tindices=(1:3);
end
if nargin < 7
  vlevel=-10;
end

% Find maximum grid angle size (dl)
[lat,lon,mask]=oct_read_latlonmask(gridfile,'r');
[M,L]=size(lon);
dl=1.3*max([max(max(abs(lon(2:M,:)-lon(1:M-1,:)))) ...
        max(max(abs(lon(:,2:L)-lon(:,1:L-1)))) ...
        max(max(abs(lat(2:M,:)-lat(1:M-1,:)))) ...
        max(max(abs(lat(:,2:L)-lat(:,1:L-1))))]);

% Read point positions
[type,vlevel]=get_type(hisfile,vname,vlevel);
[lat,lon,mask]=oct_read_latlonmask(gridfile,type);
[M,L]=size(lon);

% Find minimal subgrids limits
minlon=min(lonsec)-dl;
minlat=min(latsec)-dl;
maxlon=max(lonsec)+dl;
maxlat=max(latsec)+dl;
sub=lon>minlon & lon<maxlon & lat>minlat & lat<maxlat;
if (sum(sum(sub))==0)
  error('Section out of the domain')
end
ival=sum(sub,1);
jval=sum(sub,2);
imin=min(find(ival~=0));
imax=max(find(ival~=0));
jmin=min(find(jval~=0));
jmax=max(find(jval~=0));

% Get subgrids
lon=lon(jmin:jmax,imin:imax);
lat=lat(jmin:jmax,imin:imax);
sub=sub(jmin:jmax,imin:imax);
mask=mask(jmin:jmax,imin:imax);

% Put latitudes and longitudes of the section in the correct vector form
if (length(lonsec)==1)
  disp(['N-S section at longitude: ',num2str(lonsec)])
  if (length(latsec)==1)
    error('Need more points to do a section')
  elseif (length(latsec)==2)
    latsec=(latsec(1):dl:latsec(2));
  end
  lonsec=0.*latsec+lonsec;
elseif (length(latsec)==1)
  disp(['E-W section at latitude: ',num2str(latsec)])
  if (length(lonsec)==2)
    lonsec=(lonsec(1):dl:lonsec(2));
  end
  latsec=0.*lonsec+latsec;
elseif (length(lonsec)==2 & length(latsec)==2)
  Npts=ceil(max([abs(lonsec(2)-lonsec(1))/dl ...
                  abs(latsec(2)-latsec(1))/dl]));
  if lonsec(1)==lonsec(2)
    lonsec=lonsec(1)+zeros(1,Npts+1);
  else
    lonsec=(lonsec(1):(lonsec(2)-lonsec(1))/Npts:lonsec(2));
  end
  if latsec(1)==latsec(2)
    latsec=latsec(1)+zeros(1,Npts+1);
  else
    latsec=(latsec(1):(latsec(2)-latsec(1))/Npts:latsec(2));
  end
elseif (length(lonsec)~= length(latsec))
  error('Section latitudes and longitudes are not of the same length')
end
Npts=length(lonsec);

% Get the subgrid
sub=0*lon;
for i=1:Npts
  sub(lon>lonsec(i)-dl & lon<lonsec(i)+dl & ...
      lat>latsec(i)-dl & lat<latsec(i)+dl)=1;
end

% Get the variable for the first time step.
VAR=zeros(length(tindices),Npts);
tindex=tindices(1);
var=get_hslice(hisfile,gridfile,vname,tindex,vlevel,type);
var=var(jmin:jmax,imin:imax);
mask=isfinite(var.*mask);
var=var(sub==1 & mask);

% Get the coefficients of the objective analysis
londata=lon(sub==1 & mask);
latdata=lat(sub==1 & mask);
coef=oacoef(londata,latdata,lonsec,latsec,100e3);

% Get the variable for the first time step.
mvar=mean(var);
var=mvar+coef*(var-mvar);
VAR(tindex,:)=var';

% Loop on time
for i=2:length(tindices)
  tindex=tindices(i);
  var=get_hslice(hisfile,gridfile,vname,tindex,vlevel,type);
  var=var(jmin:jmax,imin:imax);
  var=var(sub==1 & mask);
  mvar=mean(var);
  var=mvar+coef*(var-mvar);
  VAR(tindex,:)=var';
end

% Get the distances
x=spheric_dist(latsec(1),latsec,lonsec(1),lonsec)/1e3;

% Get the time from history file (native oct_netcdf)
try
  ncid = oct_netcdf.open(hisfile,'NC_NOWRITE');
  try
    t_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'scrum_time'));
  catch
    try
      t_full = oct_netcdf.getVar(ncid,oct_netcdf.inqVarID(ncid,'ocean_time'));
    catch
      t_full = [];
    end
  end
  if ~isempty(t_full)
    t = double(t_full(tindices));
  else
    t = tindices;
  end
  oct_netcdf.close(ncid);
catch
  t = tindices;
end

[X,T]=meshgrid(x,t);
return