function time_series = make_timeseries

%========== Options ==========%
filename = '/esdata/env/dpn06xsu/job/test_data/ERA40_test.nc';
var = 'T2';

%========== Data info ==========%
time_levs = 1460;		% Number of time-levels (dim 4) in netCDF file
dt = 0.25;			% Interval between time-levels (units of choice)
lat_name = 'latitude';		% Name of latitude variable in NetCDF file
lon_name = 'longitude';		% Name of longitude variable in NetCDF file
x_size = 320;			% Total number of gridpoints in the x-direction
y_size = 160;			% Total number of gridpoints in the y-direction
grid_resolution = 125;		% Resolution of the dataset (km)

%======= Interpolation options =======%
interp_to_lon = 0;		% Longitude of timeseries
interp_to_lat = 56;		% Latitude of timeseries
interp_radius = 00;		% Radius of interpolation (km). Set to 0 to use nearest gridpoints


%========== CODE STARTS BELOW =============%

% Get longitude and latitude variables
lat = double(ncread(filename,lat_name));
lon = double(ncread(filename,lon_name));

% Initialise timeseries
time_series = zeros(1,time_levs);

% Guess at the range of the dataset we have to scan
n_points = ceil(sqrt(2)*interp_radius/grid_resolution);

% Find the nearest grid points to interp_to_[lon/lat]
lon_near = min(min(dnearest(lon,interp_to_lon)));
lat_near = min(min(dnearest(lat,interp_to_lat)));

% Make the range of grid points to scan over
lon_lim = [lon_near-n_points lon_near+n_points];
lat_lim = [lat_near-n_points lat_near+n_points];

% Meshgrid and transpose lon/lat arrays
[lon lat] = meshgrid(lon,lat);
% lon = lon;lat = lat;

% Cycle through time levs, interpolating to a position if needed
for(t=1:time_levs)

  % Initialise arrays for interpolation
  lts = [];
  lns = [];
  vals = [];

  % Load the current time level of data
  field = double(ncread(filename,var,[t 1 1 1]-1,[1 1 y_size x_size]));
  % Scan over an area around interp_to_[lon/lat], add any grid points suitable close
  % to interp_to_[lon/lat] to the interpolation arrays
  cnt = 1;
  for(i=lon_lim(1):lon_lim(2))
    for(j=lat_lim(1):lat_lim(2))
      if(i<=0) 
	i2 = i+x_size;
      else
	i2 = i;
      end
      if(max(gc_dist(interp_to_lon,interp_to_lat))<=interp_radius)
	lts(cnt) = lat(j,i2);
	lns(cnt) = lon(j,i2);
	vals(cnt) = field(j,i2);
	cnt = cnt + 1;
      end
    end
  end

  % If interp_radius is > 0, interpolate all values in interpolation arrays to
  % interp_to_[lon/lat]. Otherwise just add the nearest grid point to the timeseries.
  if(interp_radius > 0)
    time_series(t) = griddata(lns,lts,vals,interp_to_lon,interp_to_lat);
  else
    time_series(t) = vals;
  end

end

end 



%=============== FUNCTIONS ================%

function index = dnearest(array,value)
% index = nearest(array,value)

array = abs(array - value);

index = find(array == min(min(array)));

end 

function d = gc_dist(lon, lat)
%Returns the cumulative great circle distance from a series of ([lon] [lat])
%points. Units will be in units of R (currently km)

R = 6378.1;

if(size(lon) ~= size(lat))
    error('lon and lat arrays must be the same size');
end

lon = lon*pi/180;
lat = lat*pi/180;

d = zeros(size(lon));

for(i=2:length(d))
    dlon = lon(i) - lon(i-1);
    dlat = lat(i) - lat(i-1);
    a = (sin(dlat/2)).^2 + cos(lat(i-1)) * cos(lat(i)) * (sin(dlon/2))^2;
    c = 2 * asin(min(1,sqrt(a)));
    d(i) = d(i-1) + (R * c);
end
end