% read in patient info
M = csvread('patient_info_stoch.csv');
% rate of decrease of pain without drugs or disease
relax_rate = M(1);
% baseline pain, equilibrium for patient in absence of drugs
base_pain = M(2);
% initial pain
init_pain = M(3);
% pain level noise
eps = M(4);

% read in simulation info
M = csvread('sim_info.csv');
% initial time
t0 = M(1);
% end time
tmax = M(2);
% time resolution
dt = M(3);
% time vector - change to any vector you want - e.g. could be:
%             tvec = [0 3.4 5.1 7 9.4 15 19.8 24];
tvec = t0:dt:tmax;

% read in drug info
drug_info = csvread('drug_info.csv');
% % marginal effect of drug 1 on pain relaxation rate
% k1 = M(1);
% % rate of elimination of drug 1
% kD1 = M(2);
% % amplitude of response to drug 1
% a1 = M(3);
% % drug information
% drug_info = [k1, kD1, a1];
% read in drug dosage
dose_times = csvread('dose_times.csv');

% read in "real" patient data
% M = csvread('pain_report_real.csv');
% report_times_real = M(1,:);
% report_pain_real = M(2,:);

noise = 0.01;

% all the b values you want to get dynamic pain values for
test_b = [1 1.5 2 2.5 4 6 8];

% run forloop to test all the b values you want:
for ii = 1:length(test_b)
    % run the simulation for each b value
    [tvec,P,D] = pain_sim_drug1_stoch(tvec,init_pain,relax_rate,test_b(ii),...
                                      noise,drug_info,dose_times);
    % write each simulation to a csv file
    csvwrite(['sim_pain_',num2str(ii),'.csv'], [tvec; P']);
end