%% Background countrate and doserate - detector is NaI(Tl) - cylindrical 
% (D=5.1 cm, H=5.1 cm)

% activity 
Ba133=2.63E6;
Cs137=3.97E6;
% Actual number of elements is numOfElements^2
numOfElements=20;
% D ,H (Model 44-62 Gamma detector)
D=5.1;% cm
H=5.1;% cm
% definition of detector
detectorAttenuation=[5.7263,0.1657,0.0765,0.0792];% cm^-1 linear attenuation coefficient
detectorSize=[D;H]; % cm
detectorHeight=50;% cm
detectorPosition=[50;50;detectorHeight];% [cm;cm;cm] center of base of the detector
density=3.667E-3; % kg/cm^3
detector=translate2struct(detectorPosition,detectorAttenuation,detectorSize,density);
%
photonEnergy=[0.08,0.30,0.36,0.662]*10^6; % eV
% definition of background
backgroundSourceActivity=[Ba133*0.34;Ba133*0.329;Ba133*0.62;Cs137*0.851];% Bq
backgroundApproxLimit=5; 
x=linspace(0,100,backgroundApproxLimit);
sourcePosition=zeros(3,backgroundApproxLimit^2);
for n=1:backgroundApproxLimit
    for m=1:backgroundApproxLimit
        sourcePosition(:,backgroundApproxLimit*(m-1)+n)=[x(n);x(m);0];
    end
end
sourceActivity=repmat(backgroundSourceActivity/(size(sourcePosition,2)),1,size(sourcePosition,2));
% calculation
[countRateOutput,countRateOtoutputPhotonEnergy,countRateRawData,doseRateRateOutput,doseRateOtoutputPhotonEnergy,doseRateRawData]=runModel( sourcePosition,sourceActivity,photonEnergy,detector,numOfElements);
% count rate figure
figure(1);hold on;
bar([countRateOtoutputPhotonEnergy(:,:,1),countRateOtoutputPhotonEnergy(:,:,2),countRateOtoutputPhotonEnergy(:,:,3),countRateOtoutputPhotonEnergy(:,:,4)])
ax=gca;
ax.XTick=1:4;
names={'0.08','0.3','0.36','0.662'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Counts per second');
set(ax,'ygrid','on') 
% dose rate figure
figure(2);hold on;
bar([doseRateOtoutputPhotonEnergy(:,:,1),doseRateOtoutputPhotonEnergy(:,:,2),doseRateOtoutputPhotonEnergy(:,:,3),doseRateOtoutputPhotonEnergy(:,:,4)])
ax=gca;
ax.XTick=1:4;
names={'0.08','0.3','0.36','0.662'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Dose rate [Gy/s]');
set(ax,'ygrid','on') 
%% Countrate and doserate - NaI(Tl)- cylindrical (D=5.1 cm, H=5.1 cm),
% Single spring laying in the corner of the tray

% activity 
Ba133=2.63E6;
Cs137=3.97E6;
Co60=44.2E6;
% Actual number of elements is numOfElements^2
numOfElements=20;
% D ,H (Model 44-62 Gamma detector)
D=5.1;% cm
H=5.1;% cm
% definition of detector
%  gamma   [MeV]     0.08,  0.30,  0.36,  0.662, 1.17,  1.33
detectorAttenuation=[5.7263,0.1657,0.0765,0.0792,0.0535,0.05];% cm^-1 linear attenuation coefficient
detectorSize=[D;H]; %cm
detectorHeight=50;% cm
detectorPosition=[50;50;detectorHeight];% [cm;cm;cm] center of base of the detector
density=3.667E-3;% kg/cm^3
detector=translate2struct(detectorPosition,detectorAttenuation,detectorSize,density);
%
photonEnergy=[0.08;0.30;0.36;0.662;1.17; 1.33]*10^6;% eV
% definition of background
backgroundSourceActivity=[Ba133*0.34;Ba133*0.329;Ba133*0.62;Cs137*0.851;0;0];% eV
backgroundApproxLimit=5;
x=linspace(0,100,backgroundApproxLimit);
sourcePosition=zeros(3,backgroundApproxLimit^2+1);
for n=1:backgroundApproxLimit
    for m=1:backgroundApproxLimit
        sourcePosition(:,backgroundApproxLimit*(m-1)+n+1)=[x(n);x(m);0];
    end
end
sourceActivity=repmat(backgroundSourceActivity/(size(sourcePosition,2)-1),1,size(sourcePosition,2));
% definition of spring
sourcePosition(:,1)=[0;0;0]; % [cm;cm;cm]
sourceActivity(:,1)=[0;0;0;0;Co60*0.9988;Co60]; %Bq spring activity
% calculation
[countRateOutput,countRateOtoutputPhotonEnergy,countRateRawData,doseRateRateOutput,doseRateOtoutputPhotonEnergy,doseRateRawData]=runModel( sourcePosition,sourceActivity,photonEnergy,detector,numOfElements);
% count rate figure
figure(3);hold on;
bar([countRateOtoutputPhotonEnergy(:,:,1),countRateOtoutputPhotonEnergy(:,:,2),countRateOtoutputPhotonEnergy(:,:,3),countRateOtoutputPhotonEnergy(:,:,4),countRateOtoutputPhotonEnergy(:,:,5),countRateOtoutputPhotonEnergy(:,:,6)])
ax=gca;
ax.XTick=1:6;
names={'0.08','0.3','0.36','0.662','1.17','1.33'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Counts per second');
set(ax,'ygrid','on') 
figure(4);hold on;
% dose rate figure
bar([doseRateOtoutputPhotonEnergy(:,:,1),doseRateOtoutputPhotonEnergy(:,:,2),doseRateOtoutputPhotonEnergy(:,:,3),doseRateOtoutputPhotonEnergy(:,:,4),doseRateOtoutputPhotonEnergy(:,:,5),doseRateOtoutputPhotonEnergy(:,:,6)])
ax=gca;
ax.XTick=1:6;
names={'0.08','0.3','0.36','0.662','1.17','1.33'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Dose rate [Gy/s]');
set(ax,'ygrid','on') 
%% Countrate and doserate - NaI(Tl)- cylindrical (D=5.1 cm, H=5.1 cm),
% spring laying in the middle of the tray

% activity 
Ba133=2.63E6;%Bq
Cs137=3.97E6;%Bq
Co60=44.2E6; %Bq
% Actual number of elements is numOfElements^2%
numOfElements=20;
% D ,H (Model 44-62 Gamma detector)
D=5.1;%cm
H=5.1;%cm
% definition of detector
%  gamma   [MeV]     0.08,  0.30,  0.36,  0.662, 1.17,  1.33
detectorAttenuation=[5.7263,0.1657,0.0765,0.0792,0.0535,0.05];% cm^-1 linear attenuation coefficient
detectorSize=[D;H]; %cm
detectorHeight=50;%cm
detectorPosition=[50;50;detectorHeight];% [cm;cm;cm]  center of base of the detector
density=3.667E-3;%kg/cm^3
detector=translate2struct(detectorPosition,detectorAttenuation,detectorSize,density);
%
photonEnergy=[0.08;0.30;0.36;0.662;1.17; 1.33]*10^6;
% definition of background
backgroundSourceActivity=[Ba133*0.34;Ba133*0.329;Ba133*0.62;Cs137*0.851;0;0];% Bq
backgroundApproxLimit=5;
x=linspace(0,100,backgroundApproxLimit);
sourcePosition=zeros(3,backgroundApproxLimit^2+1);
for n=1:backgroundApproxLimit
    for m=1:backgroundApproxLimit
        sourcePosition(:,backgroundApproxLimit*(m-1)+n+1)=[x(n);x(m);0];
    end
end
sourceActivity=repmat(backgroundSourceActivity/(size(sourcePosition,2)-1),1,size(sourcePosition,2));
% definition of spring
sourcePosition(:,1)=[50;50;0]; % [cm;cm;cm]
sourceActivity(:,1)=[0;0;0;0;Co60*0.9988;Co60]; % Bq spring activity
% calculation
[countRateOutput,countRateOtoutputPhotonEnergy,countRateRawData,doseRateRateOutput,doseRateOtoutputPhotonEnergy,doseRateRawData]=runModel( sourcePosition,sourceActivity,photonEnergy,detector,numOfElements);
% count rate figure
figure(5);hold on;
bar([countRateOtoutputPhotonEnergy(:,:,1),countRateOtoutputPhotonEnergy(:,:,2),countRateOtoutputPhotonEnergy(:,:,3),countRateOtoutputPhotonEnergy(:,:,4),countRateOtoutputPhotonEnergy(:,:,5),countRateOtoutputPhotonEnergy(:,:,6)])
ax=gca;
ax.XTick=1:6;
names={'0.08','0.3','0.36','0.662','1.17','1.33'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Counts per second');
set(ax,'ygrid','on') 
figure(6);hold on;
% dose rate figure
bar([doseRateOtoutputPhotonEnergy(:,:,1),doseRateOtoutputPhotonEnergy(:,:,2),doseRateOtoutputPhotonEnergy(:,:,3),doseRateOtoutputPhotonEnergy(:,:,4),doseRateOtoutputPhotonEnergy(:,:,5),doseRateOtoutputPhotonEnergy(:,:,6)])
ax=gca;
ax.XTick=1:6;
names={'0.08','0.3','0.36','0.662','1.17','1.33'};
ax.XTickLabels=names;
xlabel('Initial energy of photons [MeV]');
ylabel('Dose rate [Gy/s]');
set(ax,'ygrid','on') 
