%% load dose response and fading runs
load drcData.mat
load fadingData.mat

LabDoseRate=0.09397;%Gy/s

%% do Monte Carlo runs to propogate various error sources
%all disc positions for fading runs
allPosnsFade=[9,11,14,16,19,21,24,27,30];
%effective fading times (s) between irradiation and measurement of dose
%response points
eft50=[0;457;688;1152;2078;3933;7641;0;457];
eft225=[0;624;854;1318;2245;4100;7808;0;624];
labDosesSbeta=[0;625;1250;2500;5000;10000;20000;0;625];
geolDoseRateArray=[4.408 0.879;4.408 0.879;4.482 0.879;4.482 0.879;4.493 0.879;4.493 0.879;4.622 0.879;4.628 0.879;4.648 0.879];%entry per position

numMCruns=50;

for position_i=1:length(allPosnsFade)%% loop through positions
    for signalTemp=[50,225]% and stimulation temperature
        GeolDoseRate=geolDoseRateArray(position_i,1);
        GeolDoseRateErr=geolDoseRateArray(position_i,2);
        
        %load fade data
        eval(['fadeLiTi=p' num2str(allPosnsFade(position_i)) 'LT' num2str(signalTemp) ';'])
        eval(['fadeti=p' num2str(allPosnsFade(position_i)) 't' num2str(signalTemp) ';'])   
        
        %load dose response data
        eval(['drcLxTx=p' num2str(allPosnsFade(position_i)) '_' num2str(signalTemp) ';'])
        eval(['drcEft=eft' num2str(signalTemp) ';'])
        
        if position_i==2%p.11 has 1 dose point omitted due to sensitivity change after run break
            drc_i=[labDosesSbeta([1:5,7:end]),drcLxTx,drcEft([1:5,7:end])];
        else
            drc_i=[labDosesSbeta,drcLxTx,drcEft];
        end
        fade_i=[fadeLiTi,fadeti];
        if signalTemp==50
            titleString=['p.' num2str(allPosnsFade(position_i)) ', IR_{50}'];
        else
            titleString=['p.' num2str(allPosnsFade(position_i)) ', pIRIR_{225}'];
        end
        [age_MCarray,age2D0]=main_MC(drc_i,fade_i,LabDoseRate,GeolDoseRate,GeolDoseRateErr,titleString,numMCruns);
        
        eval(['p' num2str(allPosnsFade(position_i)) 't' num2str(signalTemp) '_ageMCarray=age_MCarray;'])  
        eval(['p' num2str(allPosnsFade(position_i)) 't' num2str(signalTemp) '_2D0=age2D0;'])  
    end
    save ageMCarrays
end