Xray spectrum and attenuation

data from program spectrum

Xray = {8.000000, 0.002593, 0.002593, 0.000000, 9.000000, 0.006275, 0.006275, 0.000000, 10.000000, 0.011130, 0.011130, 0.000000, 11.000000, 0.016319, 0.016319, 0.000000, 12.000000, 0.021172, 0.021172, 0.000000, 13.000000, 0.025313, 0.025313, 0.000000, 14.000000, 0.028608, 0.028608, 0.000000, 15.000000, 0.031068, 0.031068, 0.000000, 16.000000, 0.032855, 0.032855, 0.000000, 17.000000, 0.033984, 0.033984, 0.000001, 18.000000, 0.034572, 0.034572, 0.000007, 19.000000, 0.034726, 0.034726, 0.000037, 20.000000, 0.034541, 0.034541, 0.000144, 21.000000, 0.034105, 0.034105, 0.000304, 22.000000, 0.033467, 0.033467, 0.000588, 23.000000, 0.032675, 0.032675, 0.001056, 24.000000, 0.031773, 0.031773, 0.001777, 25.000000, 0.030792, 0.030792, 0.002830, 26.000000, 0.029758, 0.029758, 0.004296, 27.000000, 0.028691, 0.028691, 0.006255, 28.000000, 0.027605, 0.027605, 0.008777, 29.000000, 0.026514, 0.026514, 0.011922, 30.000000, 0.025427, 0.025427, 0.015732, 31.000000, 0.024352, 0.024352, 0.017684, 32.000000, 0.023292, 0.023292, 0.019678, 33.000000, 0.022253, 0.022253, 0.021697, 34.000000, 0.021235, 0.021235, 0.023721, 35.000000, 0.020243, 0.020243, 0.025731, 36.000000, 0.019277, 0.019277, 0.027707, 37.000000, 0.018338, 0.018338, 0.029632, 38.000000, 0.017428, 0.017428, 0.031488, 39.000000, 0.016545, 0.016545, 0.033259, 40.000000, 0.015690, 0.015690, 0.034929, 41.000000, 0.014862, 0.014862, 0.035205, 42.000000, 0.014062, 0.014062, 0.035364, 43.000000, 0.013289, 0.013289, 0.035407, 44.000000, 0.012542, 0.012542, 0.035336, 45.000000, 0.011820, 0.011820, 0.035151, 46.000000, 0.011123, 0.011123, 0.034854, 47.000000, 0.010451, 0.010451, 0.034446, 48.000000, 0.009801, 0.009801, 0.033930, 49.000000, 0.009174, 0.009174, 0.033308, 50.000000, 0.008569, 0.008569, 0.032582, 51.000000, 0.007985, 0.007985, 0.031400, 52.000000, 0.007421, 0.007421, 0.030155, 53.000000, 0.006876, 0.006876, 0.028850, 54.000000, 0.006351, 0.006351, 0.027488, 55.000000, 0.005843, 0.005843, 0.026072, 56.000000, 0.005353, 0.005353, 0.024605, 57.000000, 0.004880, 0.004880, 0.023089, 58.000000, 0.004423, 0.004423, 0.021526, 59.000000, 0.003982, 0.003982, 0.019921, 60.000000, 0.003556, 0.003556, 0.018274, 61.000000, 0.003144, 0.003144, 0.016541, 62.000000, 0.002745, 0.002745, 0.014782, 63.000000, 0.002360, 0.002360, 0.013000, 64.000000, 0.001988, 0.001988, 0.011196, 65.000000, 0.001629, 0.001629, 0.009371, 66.000000, 0.001281, 0.001281, 0.007528, 67.000000, 0.000945, 0.000945, 0.005668, 68.000000, 0.000619, 0.000619, 0.003793, 69.000000, 0.000305, 0.000305, 0.001903, 70.000000, 0.000000, 0.000000, 0.000000, 4, 57.980, 34068396, 34068396, 341476, 59.320, 63089623, 63089623, 656110, 67.200, 20819575, 20819575, 258499, 69.100, 5047170, 5047170, 65109};

Dimensions[Xray]

XrayP  =  Partition[Xray, 4];

Dimensions[XrayP]

XrayPT = Transpose[XrayP];

plot of Xray spectrum

ListPlot[
       Table[{XrayPT[[1]][[n]], XrayPT[[2]][[n]]}, {n,
           1, 62}], {PlotRange -> All,  PlotJoined -> True,          PlotStyle -> Thickness[0.01]}]

attenuation functions

f[E_] :=
         0.597 * (10^-24)*           Exp[-0.0028 (E - 30)];

E is in keV

NA  =  6.02 * 10^23;

μcompton[E_, ρ_, Z_,            A_] :=  ρ NA (Z/A) f[E];

μphoto[E_,  ρ_,  Z_,            A_] :=  ρ NA (Z/A) 9.8 *           10^-24 (Z^3.8 /E^3.2);

ρ is density
Z is atomic number
A is atomic mass

function for spectral change

\(\(spectrumcomp[ρ_, \ Z_, \ A_, \ l_]\  :=          Table[{\(XrayPT[\([1]\)]\)[\([n]\)], 
		\(XrayPT[\([2]\)]\)[\([n]\)]\ Exp[\(-μcompton[\(XrayPT[\([1]\)]\)[\([n]\)], ρ, Z, A]\)\ *\ l]}, {n, 1, 62}];\)\)

\(\(spectrumphoto[ρ_, \ Z_, \ A_, \ l_]\  :=          Table[{\(XrayPT[\([1]\)]\)[\([n]\)], 
		\(XrayPT[\([2]\)]\)[\([n]\)]\ Exp[\(-μphoto[\(XrayPT[\([1]\)]\)[\([n]\)], ρ,
                       Z, A]\)\ *\ l]}, {n, 1, 62}];\)\)

\(\(spectrumtotal[ρ_, \ Z_, \ A_, \ l_]\  :=          Table[{\(XrayPT[\([1]\)]\)[\([n]\)], 
		\(XrayPT[\([2]\)]\)[\([n]\)]\ Exp[\(-\((μphoto[\(XrayPT[\([1]\)]\)[\([n]\)], ρ, Z, A] + μcompton[\(XrayPT[\([1]\)]\)[\([n]\)], ρ, Z,
                           A]\ )\)\)\ *\ l]}, {n, 1, 62}];\)\)

l is length of attenuator in cm

spectral change in water

Zwater  =  7.54;
Awater  =  16;

spectrumComptonWater[l_] :=
         spectrumcomp[1, Zwater, Awater, l];

spectrumPhotoWater[l_] :=
         spectrumphoto[1, Zwater, Awater, l];

spectrumTotalWater[l_] :=
         spectrumtotal[1, Zwater, Awater, l];

ListPlot[
       spectrumComptonWater[1], {PlotRange -> All,          PlotJoined -> True,          PlotStyle -> Thickness[0.01]}]

waterCompton =
       Show[Table[
           ListPlot[
             spectrumComptonWater[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

waterPhoto =
       Show[Table[
           ListPlot[
             spectrumPhotoWater[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

waterTotal =
       Show[Table[
           ListPlot[
             spectrumTotalWater[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

Show[waterCompton]

Show[waterPhoto]

Show[waterTotal]

spectral change in bone

Zbone  =  14.2;
Abone  =  21.4;

spectrumComptonBone[l_] :=
         spectrumcomp[2.93, Zbone, Abone, l];

spectrumPhotoBone[l_] :=
         spectrumphoto[2.93, Zbone, Abone, l];

spectrumTotalBone[l_] :=
         spectrumtotal[2.93, Zbone, Abone, l];

boneCompton =
       Show[Table[
           ListPlot[
             spectrumComptonBone[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

bonePhoto =
       Show[Table[
           ListPlot[
             spectrumPhotoBone[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

boneTotal =
       Show[Table[
           ListPlot[
             spectrumTotalBone[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

Show[boneCompton]

Show[bonePhoto]

Show[boneTotal]

spectral change in oil

Zoil  =  5.51;
Aoil  =  11;

spectrumComptonOil[l_] :=
         spectrumcomp[0.66, Zoil, Aoil, l];

spectrumPhotoOil[l_] :=
         spectrumphoto[0.66, Zoil, Aoil, l];

spectrumTotalOil[l_] :=
         spectrumtotal[0.66, Zoil, Aoil, l];

oilCompton =
       Show[Table[
           ListPlot[
             spectrumComptonOil[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

oilPhoto =
       Show[Table[
           ListPlot[
             spectrumPhotoOil[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

oilTotal =
       Show[Table[
           ListPlot[
             spectrumTotalOil[l], {PlotRange -> All,                PlotJoined -> True,                PlotStyle -> Thickness[0.01]}], {l, 0, 10}]]

Show[oilCompton]

Show[oilPhoto]

Show[oilTotal]

\(Show[
       GraphicsArray[{{oilCompton, oilPhoto,
             oilTotal}, 
		\ \ \ \ \ \ \ \ \ \ \ {waterCompton, \ waterPhoto,              waterTotal}, 
			{boneCompton,              bonePhoto, \ boneTotal}}]]\)

Attenuation due to a wedge shaped object

consider  a parallel beam geometry with a wedge shape, what is the beam intensity as a function of position.
Recall that intensity is the power in the spectrum, therefore the integral of the intensity times the energy.
(this neglects the variation in the stopping power as a function of energy in the detector).

water

Intensity [l_]  :=        Sum[spectrumTotalWater[
                 l][[n]][[1]] * spectrumTotalWater[
                 l][[n]][[2]], {n, 1, 62}]

Wedge = Table[Intensity[l], {l, 0, 20, 0.5}];

pWedgeWater =
       ListPlot[Wedge, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

Wedgel =
         Table[Log[Wedge[[n]]/Wedge[[1]]],  {n, 1, 
41}];

plWedgeWater =
       ListPlot[Wedgel, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

oil

IntensityOil [l_]  :=        Sum[spectrumTotalOil[
                 l][[n]][[1]] * spectrumTotalOil[
                 l][[n]][[2]], {n, 1, 62}]

WedgeOil = Table[IntensityOil[l], {l, 0, 20, 0.5}];

pWedgeOil =
       ListPlot[WedgeOil, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> {Thickness[0.01], RGBColor[1, 0, 
0]}}]

WedgelOil =
         Table[Log[WedgeOil[[n]]/WedgeOil[[1]]],  {n, 1,
             41}];

plWedgeOil =
       ListPlot[Wedgel, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> {Thickness[0.01],                RGBColor[1, 0, 0]}}]

bone

IntensityBone[l_]  :=        Sum[spectrumTotalBone[
                 l][[n]][[1]] * spectrumTotalBone[
                 l][[n]][[2]], {n, 1, 62}]

WedgeBone = Table[IntensityBone[l], {l, 0, 20, 0.5}];

pWedgeBone =
       ListPlot[WedgeBone, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> {Thickness[0.01], RGBColor[0, 0, 
1]}}]

WedgelBone =
         Table[Log[WedgeBone[[n]]/WedgeBone[[1]]],  {n, 1,
             41}];

plWedgeBone =
       ListPlot[WedgelBone, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> {Thickness[0.01], RGBColor[0, 0, 
1]}}]

comparison

Show[pWedgeWater, pWedgeOil,
       pWedgeBone, {PlotRange -> All}]

Show[plWedgeWater, plWedgeOil,
       plWedgeBone, {PlotRange -> All}]

Show[plWedgeWater, plWedgeOil,
       plWedgeBone, {PlotRange -> {0, -4}}]

Attenuation due to a cylender in a rectangular shaped object

consider  a parallel beam geometry
cylender is bone and rectangle is water
z is length along detector
l is length of water
w is width of water
r is radius of bone

\(\(spectrumTotalObject[lw1_, lb_, lw2_] =          Table[{\(XrayPT[\([1]\)]\)[\([n]\)], 
		\(XrayPT[\([2]\)]\)[\([n]\)]\ Exp[\(-\((μphoto[\(XrayPT[\([1]\)]\)[\([n]\)], 1,
                           Zwater,
                           Awater] + μcompton[\(XrayPT[\([1]\)]\)[\([n]\)], 1, Zwater,
                           Awater]\ )\)\)\ *\ \((lw1\  + \ lw2)\)]\ *               Exp[\(-\((μphoto[\(XrayPT[\([1]\)]\)[\([n]\)], 2.93, Zbone,
                           Abone] + μcompton[\(XrayPT[\([1]\)]\)[\([n]\)],
                           2.93, Zbone, Abone]\ )\)\)\ *\ lb]}, {n, 1,
             62}];\)\)

Intensity [l_, w_, r_, z_]  :=        If[Abs[z] >
           l, 
Sum[XrayPT[[1]][[n]] * XrayPT[[2]][[n]],  {n, 1, 62}], 
   If[
           Abs[z] >
             r, 
     Sum[spectrumTotalWater[
                     w][[n]][[1]] * spectrumTotalWater[
                     w][[n]][[2]], {n, 1,
               62}], 
          Sum[spectrumTotalObject[w/2 - Sqrt[r^2 - z^2], Sqrt[r^2 - z^2],                      w/2 - Sqrt[
                         r^2 - z^2]][[n]][[1]] *     
                  spectrumTotalObject[w/2 - Sqrt[r^2 - z^2], Sqrt[r^2 - z^2],                      w/2 - Sqrt[r^2 - z^2]][[n]][[2]], {n, 1,
               62}]]]

Object1 =
         Table[Intensity[5, 10, 2, z], {z, -10, 10,  .5}];

pObject1 =
       ListPlot[Object1, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

Object2 =
         Table[Intensity[5, 10, 1, z], {z, -5, 5,  .25}];

pObject2 =
       ListPlot[Object2, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

Object3 =
         Table[Intensity[5, 10, 0.5, z], {z, -2.5,
             2.5,  .125}];

pObject3  =        ListPlot[Object3, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

Object4 =
         Table[Intensity[5, 10, 0.25, z], {z, -1.25,
             1.25,  .0625}];

pObject4 =
       ListPlot[Object4, {PlotRange -> All, PlotJoined -> True,
           PlotStyle -> Thickness[0.01]}]

Comparison

Note that the horizontal axis is incorrect in this composite plot but the vertical (intensity) axis is right.

Show[pObject1, pObject2, pObject3, pObject4]

Show[pObject1, pObject2, pObject3,
       pObject4,  {PlotRange -> {0, 4}}]