forked from I-STAR/SPEKTR
-
Notifications
You must be signed in to change notification settings - Fork 0
/
spektrBeers.m
140 lines (121 loc) · 5.22 KB
/
spektrBeers.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
function f = spektrBeers(q,elem_filters)
%%**************************************************************************
%% System name: SPEKTR
%% Module name: spektrBeers.m
%% Version number: 3
%% Revision number: 00
%% Revision date: 30-Jun-2015
%%
%% 2016 (C) Copyright by Jeffrey H. Siewerdsen.
%% I-STAR Lab
%% Johns Hopkins University
%%
%% Usage: f = spektrBeers(q, elem_filters)
%%
%% Inputs:
%% 'q' - X-Ray Energy Spectrum (is a [150 x 1] matrix), generated
%% from the function spektrSpectrum(...,..)
%% 'elem_filters' - [N x 2] matrix containing filters which are listed as
%% atomic number(column 1) & numeric filter thickness (column 2).
%% OR
%% {N x 2} cell with atomic number or symbol (column 1) &
%% numeric filter thickness (column 2)
%%
%% ie. f = spektrBeers(q,[13 1; 92 4]);
%% ie. f = spektrBeers(q, {'Al' 1; 'U' 4});
%% ie. f = spektrBeers(q, {13 1; 92 4});
%% Filter: 1mm of Aluminum & 4mm of Uranium
%%
%% Outputs:
%% Filtered X-ray Energy Spectrum, which is a 150x1 matrix, each matrix
%% element representing the # of photons / mAs / mm^2 at 100 cm from the source
%% for a 1 keV energy bin (energy bins range from 1-150 keV)
%%
%% Description:
%% This function will filter the energy spectrum, q, given a number of
%% filters & filter thicknesses
%%
%% Notes:
%%
%%*************************************************************************
%% References:
%%
%%*************************************************************************
%% Revision History
%% 0.000 2003 05 01 AW Initial code
%% 1.000 2004 03 15 DJM Initial released version
%% 2.000 2006 04 19 DJM Removed XLSread and replaced with .MAT
%% 3.000 2015 06 30 JGP Accept both atomic symbol or atomic number
%% as arguments
%%*************************************************************************
%%
% Date: May 2003
% Revisions History: June 17, 2003 - addition of " q=filtered; " at the
% end of the main for loop.
% By: Aaron Waese
% Where elem_filters will be a matrix [Nx2] which will indicate which filters will
% be used and their respective thicknesses. ie. [13 1;92 4]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% PARAMETERS
% Energy Vector
EnergyVector = 1:150;
% Column Identifiers for XL database sheets
%%% ... in xls database of mu/rho for elements
%%% AND for the xls database for mu/rho for compounds
XLColumn_LinearAttenuation = 9;
%%% ... in xls database for density of elements and compounds
XLColumn_Density = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize the matrix of the filtered spectra
filtered = zeros(150,1);
% Read in density of filter from density.xls
% This file consists of a 92x2 matrix with the following structure:(elements,densities)
% column 3 contains the censities in g/cm^3
load('spektrDensityElements.mat');
load('spektrMuRhoElements.mat'); % Let A be the temp matrix which has dim. 150x10
% Filter the unfiltered spectrum q, as many times as their are filters
for i=1:1:size(elem_filters,1),
% Read in attenuation coefficient data for the specified filter from
% murho_elements.mat
% In file 'murho_elements' there exists 92 sheets, each sheet contains elemental
% attenuation coefficients and attenuation coefficients (energy based).
% The structure of the data in the file is the following
% NIST DATA;
% A)column 1(energy) [MeV]
% B)column 2(attentuation/density)
% C)column 3(attentuation/density)en
% D)column 4(energy) [keV]
% E)column 5(attentuation/density)
% F)column 6(attentuation/density)en
% INTERPOLATED DATA;
% G)column 7 empty
% H)column 8(energy) [keV], 1 keV intervals, from 1-150 keV
% I)column 9(attentuation/density)
% J)column 10(attentuation/density)en
if iscell(elem_filters)
%If the row contains a string, I will turn it into the atomic number
%Else just store the atomic number in atomicNum.
if ischar(elem_filters{i,1})
atomicNum = spektrElement2Z(elem_filters{i,1});
else
atomicNum = elem_filters{i,1};
end
thickness = elem_filters{i,2};
else
atomicNum = elem_filters(i,1);
thickness = elem_filters(i,2);
end
% Read in attenuation coefficient data for the specified filter from
% murho_elements.xls
attenuation = A{atomicNum}(EnergyVector,XLColumn_LinearAttenuation); % 150x1 matrix (contains mass attenuation coefficients)
% take the 9th column of the spreadsheet A
% Convert interpolated data to attenuation via the following expression:
% attenuation = (attenuation/density)* density => [cm^2/g]*[g/cm^3]=[1/cm]
% The 0.1 accnts for converting the attenuation matrix to units of [1/mm]
attenuation = density(atomicNum,XLColumn_Density)*attenuation*(0.1);
% Apply Beer's Law to the spectrum, given a filter type/thickness
filtered = q.*exp(-attenuation*thickness);
q = filtered;
end
% Return filtered spectrum [150x1] MATRIX
f = filtered;