function [retval,spacenodes,cumprof] = RadialProfiles(X,Y,dr,TimeNodes,Events)
% usage: [profiles,spacenodes,cumprofiles] = RadialProfiles(X,Y,dr,TimeNodes,Events)
% description: returns summary radial profiles of all events
% between values of the vector TimeNodes. Each profile is a row
% in the matrix profiles, one row per time interval given by TimeNodes.
% The profile counts events in space intervals of width dr up
% to the limit determined by the farthest event.
% Optionally returns the spatial nodes used in spacenodes and returns the
% cumulative profiles resulting from summing over times up to a given time.
% Here are the meanings of the input parameters:
% X,Y: coordinates of the center of the diffusion process
% dr: size of concentric ring in which activity is placed
% TimeNodes: nodes which form endpoints of time intervals
% Events: array of events observed

% local variables:
% ntimes: length of TimeNodes
% event: Events filtered by TimeNodes
% nevent: number of filtered events
% radius: array of distances of events from center
% rmax: maximum radius in all events
% nspace: number of space nodes
% spacenodes: actual space nodes
% tstep: index of an event in TimeNodes
% rstep: index of an event in spacenodes
%
% Last modified: 09/05/05

% filter out irrelevant events
ntimes = length(TimeNodes);
event = TimeFilter(Events,TimeNodes(1),TimeNodes(ntimes));
nevent = size(event,1);
% find max radius and distances from center
radius = abs((event(:,1)-X)+1i*(event(:,2)-Y));
rmax = max(radius);
% number of space nodes needed
nspace = round(rmax/dr)+2;
spacenodes = linspace(0,(nspace-1)*dr,nspace);
%spacenodes
% returned value of profiles
retval = zeros(ntimes-1,nspace-1);
% loop for sorting events into appropriate time and space
for ii = 1:nevent
  % locate time interval containing the event
  tstep = find(event(ii,3) < TimeNodes);
  tstep = tstep(1)-1;
  % locate space interval containing the event
  rstep = find(radius(ii) < spacenodes);
  rstep = rstep(1)-1;
  % insert event occurrence into profiles
  retval(tstep,rstep) = retval(tstep,rstep) + 1;
end
% weight the values by appropriate coefficients
retval = retval/(pi*dr);
cumprof = retval;
for ii=2:ntimes-1
  cumprof(ii,:) = cumprof(ii-1,:) + cumprof(ii,:);
end
% calculate the profiles for N(r,t)
retval = diag(1./diff(TimeNodes))*retval* ...
             diag(1./(spacenodes(1:nspace-1)+spacenodes(2:nspace)));
% calculate the profiles for C(r,t)
cumprof = cumprof*diag(1./(spacenodes(1:nspace-1)+spacenodes(2:nspace)));
%cumprof = 0;



 