% drugexperiment.m
% an experimental script for drug data
% generation and analysis.
%
% coordinates of event window
A = 0;
B = 109;
C = 0;
E = 161;
% coordinates of center
X = 36; %X = 43;
Y = 65; %Y = 17;
% model parameters .. user must set Nsoln to correct function
M = 46000;
D = 100;
a = 0.0;
% discretization parameters
dr = 20;
T = 5;
dt = 1;
% event detection probability
P = 0.8;
% profile time
t0 = 3;
% event placement
RadialDist = 1;
% repetitions for average model
count = 5;
% for predictability...unrem the next line
rand('seed',2);
% generate the events
[events,spacenode] = RadSim(A,B,C,E,X,Y,M,D,a,dr,T,dt,P,RadialDist);
spacenode = spacenode + 1;
disp('Generated events: ')
size(events)
spacenodes = linspace(0,dr*spacenode,spacenode);
spacepoints = (spacenodes(1:spacenode-1)+spacenodes(2:spacenode))/2;
timenodes = linspace(0,dt*T,T+1);
% do radial plot of events with fronts
plotevents(A,B,C,E,X,Y,D,a,dr,T,dt,events);
% plot radial profiles of data and exact solution at select times
plotprofiles(A,B,C,E,X,Y,D,a,dr,T,dt,events);
% now get the profiles
profile = RadialProfile(events,X,Y,spacenodes,t0-dt/2,t0+dt/2);
profile = profile(1:spacenode-1); % trim profiles
figure
grid,hold on
colors = 'rgb';
title(strcat('January ',num2str(t0+1998),' profiles, P = ', num2str(P), ...
  ', M = ',num2str(M),', D = ',num2str(D),', a = ',num2str(a)));
xlabel('Radius');
ylabel('Function values')
plot(spacepoints,profile,strcat('-',colors(1)))
% drug_offenses:
dataevents = dlmread('1998_2002_DF_matlab.csv',',');
disp('Actual events: ')
size(dataevents)
dataevents(:,3) = dataevents(:,3)-1998;
dataprofile = RadialProfile(dataevents,X,Y,spacepoints,t0-dt/2,t0+dt/2);
plot(spacepoints,dataprofile,strcat('-',colors(2)));
legend(strcat('Computed profile of N(r,',num2str(t0),')'), ... 
  strcat('True profile of N(r,',num2str(t0),')'));
% find all data profiles and pad, if necessary
dataprofiles = RadialProfiles(X,Y,dr,timenodes,dataevents);
dataprofiles = [dataprofiles,zeros(T,spacenode-1-size(dataprofiles,2))];
errors = dataerrornorms(A,B,C,E,X,Y,M,D,a,dr,T,dt,RadialDist,count,dataprofiles,P)