Method and Apparatus for Analyzing a Porous Nonhomogeneous Cylindrical Object

This invention relates to a method for nondestructively analyzing a porous, inhomogeneous, generally cylindrical object having n (where n is an integer) known molecular components to determine the density and the mass fraction of each of the n molecular components thereof.

Lightguide fiber of the type used to carry communication signals is fabricated by heating and drawing a portion of a lightguide preform comprised of a refractive core surrounded by glass cladding. One process which has proven extremely useful for the fabrication of preforms is the modified chemical vapor deposition technique (MCVD) whereby reactant-containing precursor gases, such as silicon tetrachloride and germanium tetrachloride are passed through the center of a starter tube. As the starter tube is heated and rotated, reaction by-products, in the form of submicron-sized doped glass particles, are deposited on the inside surface of the tube to form the refractive core of the preform. For a further, more detailed description of the MCVD process, reference should be had to U.S. Pat. No. 4,217,027 issued to J. B. MacChesney et al. on Aug. 12, 1980 and assigned to Bell Telephone Laboratories.

Other chemical vapor deposition methods such as the Vapor Phase Axial Deposition (VAD) and the Vapor Phase Radial Deposition (VRD) Techniques are presently being investigated to determine their suitability as potential processes for fabricating lightguide preforms. The VAD and the VRD techniques are both characterized by the deposit of submicron-sized doped glass particles on a bait rod to form a boule. Since the submicron-sized doped glass particles deposited on the bait rod during VAD or VRD processes appear as tiny particles of soot, the term "soot deposition" is often used to describe these two techniques. For the same reason, the boule fabricated by the VAD or VRD process is often referred to as a soot boule.

To fabricate a preform from the soot boule, the boule is first consolidated by sintering. Once sintered, the boule is then cladded by insertion into a silica glass tube. In certain instances, the soot boule may be partially cladded by depositing pure silicon dioxide thereon prior to sintering thereof.

Successful implementation of either the VAD or VRD process on a widespread basis is dependent on the quality of the preforms produced thereby. One useful technique for verifying the quality of the preform is to determine the characteristics, particularly the density and the composition, of the soot boule from which the preform is produced. Since the soot particles comprising the unconsolidated boule are generally opaque to visible light, profiling techniques which rely on the use of visible light such as that disclosed in U.S. Pat. No. 4,227,806 issued to L. S. Watkins on Oct. 14, 1980 and assigned to the assignee of the present invention, are not applicable. At present, analysis of soot boules to determine their suitability for use in the fabrication of lightguide preforms is accomplished by slicing the boule and microscopically examining the cross section thereof. Such an analysis technique renders the boule unsuited for later use.

A recent paper "Method of Inspecting Optical Fiber Preforms Using X-ray Absorption Measurements" by H. Takahashi, I. Nakamura, T. Tadatani and T. Kurowa, published in the Japanese Journal Furukawa Denko Jiho 68, pages 143-149 suggests x-ray absorption techniques used in tomographic applications, can be applied to measure the germanium concentration of preforms. However, the measurement method described in the above-referred paper is believed to be flawed since the authors do not appear to account for the possibility of having an inhomogeneous preform because they fail to disclose any method for measuring the density thereof. Consequently, the accuracy of the germanium concentration measurement obtained by this method is doubtful.

Accordingly, there is a need for a method for accomplishing nondestructive analysis of a soot boule.

The difficulties associated with characterizing the density and the composition of a porous, generally cylindrical, inhomogeneous object, such as a soot boule, comprised of n (where n is an integer), known molecular components, are overcome by the method of the present invention which comprises the steps of: scanning the object with a beam of photons having an axis perpendicular to the axis of the boule, each photon having a distinct energy within one of n-1 energy groups; measuring the intensity of the photons in each energy group upon the passage of the beam through the object at each of a plurality of separate heights t.sub.i of the beam as measured from the axis thereof to the center of the object; establishing each of n-1 attenuation coefficients .mu..sub.1 (r), .mu..sub.2 (r) . . . .mu..sub.n-1 (r) for the object at each height t.sub.i of the beam in accordance with the measured intensity of the photons within each of the n-1 energy groups, respectively; and calculating the density .rho.(r) and the mass fractions of each of the n molecular components of the object at each height t.sub.i of the beam in accordance with each of the n-1 attenuation factors .mu..sub.1 (r), .mu..sub.2 (r) . . . .mu..sub.n-1 (r).

An advantage of the present method is that the object is not physically altered during analysis thereof. Thus, the present method is well suited for use as an on-line technique for measuring the density and the mass fractions of the molecular components of soot boules fabricated by the VAD or VRD processes. From a knowledge of the mass fraction of germania in a boule, the mole fraction of germania can easily be determined which can then be used to calculate the refractive index of the boule.

FIG. 1 is a plan schematic view of an apparatus for nondestructively analyzing a soot boule;

FIG. 2 is a perspective illustration of a slab bombarded by a beam of high energy photons;

FIG. 3 is a table generally depicting the steps of the measurement method of the present invention; and

FIGS. 4 and 5, taken together, illustrate a flow chart representation of a program for computing the attentuation coefficients of the soot boule of FIG. 1.

FIG. 1 illustrates an apparatus 10 for nondestructively analyzing a generally porous, nonhomogeneous cylindrical object 12 (shown in cross section) such as an unconsolidated soot boule, comprised of n (where n is an integer) known molecular components, to measure the density and each of the n mass fractions thereof. The apparatus 10 comprises a plate 14 upon which is mounted a radiation source 16 on one side of a cutout 17 in the plate. The cutout 17 is sized to receive the boule 12. The radiation source 16 produces a beam 18 comprised of high energy photons, typically in the range of 0 to 120 keV, which bombard the soot boule 12 so as to pass therethrough in a direction perpendicular to the longitudinal axis of the boule.

In practice, the radiation source 16 comprises an x-ray tube, typically a model MCN 161 tube made by Phillips, N.V. of Holland which is excited from a suitable power supply (not shown) which comprises part of a Phillips Model MG 160 X-ray System. As compared to radiation sources which produce photons of a single energy, the x-ray tube 16 possesses a distinct advantage. The beam 18 produced by the tube 16 contains photons each having distinct energy within a continuous spectrum E.sub.K bounded by an upper and lower energy limits, E.sub.Khigh and E.sub.Klow, respectively. The magnitude of the upper energy limit E.sub.Khigh depends on the excitation of the x-ray tube 16. Typically, the upper limit E.sub.Khigh of the photon energy spectrum is 120 keV. The lower limit E.sub.Klow is essentially 0 keV.

In practice, the x-ray tube 16 is housed within a lead shield 20 which is provided with a very small opening 21 facing the cutout 17 in the plate 14 so that the beam 18 produced by the x-ray tube is well collimated. The size of the opening 21 in the shield 20 is typically on the order of 75 microns and has been greatly exaggerated in FIG. 1 for the purposes of illustration. A detector 22 typically a model 905-4 NaI scintillation detector and the circuitry associated with it made by EG&G, Ortec Division, Oakridge, Tenn., is mounted on the plate 14 on the opposite side of the cutout 17 from the x-ray tube 16 so as to be in registration with the beam 18 produced thereby. The detector 22 is housed within a lead shield 24 which has an opening 26 therein directly opposite to, and sized approximately the same as, the opening 21 in the shield 20 surrounding the x-ray tube 16.

The plate 14 is slidably mounted for movement as indicated by the double-ended arrow 27. A motor 28 has its shaft (not shown) coupled to drive a lead screw 29 in threaded engagement with a lead nut (not shown) mounted to the plate 14. As the lead screw 29 is rotatably driven by the motor 28, the plate 14 is displaced in the direction indicated by the arrow 27, causing the beam 18 to be scanned across the boule 12.

A pulse height analyzer/scaler 30, which in practice takes the form of a model 7100 multichannel analyzer made by EG&G Company, Ortec Division, is coupled to the output of the detector 22. The pulse height analyzer determines the energy dependent intensity distribution of the photons in the beam 18. The output pulses of the detector 22 are classified by the analyzer 30 according to their amplitude. In practice, the analyzer 30 is programmed to provide either 512 or 1024 channels (not shown), each channel having a storage register (not shown) associated therewith which is incremented when the amplitude of the output pulse of the detector 22 has the approximate amplitude value established for one of the particular channels. Since the magnitude of each output pulse of the detector 22 is proportional to the energy of the photon of the beam 18 impinging thereon, the analyzer effectively determines the number of the photons in the beam 18 as a function of their energy so as to effectively measure the intensity of the beam. The scaler/analyzer 30 is connected at its output to a processor 31 which takes the form of a digital computer programmed in the manner described hereinafter.

In order to fully comprehend the manner by which the apparatus 10 accomplishes nondestructive analysis of the boule 12, a brief discussion of the way in which a beam of high energy photons is attenuated upon the passage thereof through matter may be helpful. Referring to FIG. 2, there is shown a homogeneous slab 34 of a monoatomic material comprised of an atomic species i. The slab 34, which is of a thickness L, is bombarded by a beam 36 comprised of photons, each having the same energy E. As the beam 36 propagates through the slab 34, the intensity thereof is attenuated exponentially. Mathematically, the attenuation of the intensity of the beam 36 propagating through the slab 34 is given by the relationship

where I(E,O) and I(E,L) are the intensities of the beam upon entering and exiting the slab, respectively, and .mu..sub.i (E) is the attenuation coefficient of the slab.

The attenuation coefficient .mu..sub.i (E) for the slab 34, when bombarded by the beam 36 comprised of photons of energy E, is given by ##EQU1## where N.sub.a is Avogadro's number, A.sub.i is the atomic weight of the atoms of the species i comprising the slab, .rho..sub.i is the density of the slab and .sigma..sub.j.sup.(i) (E) is the probability or atomic cross section of a photon of energy (E) being removed by the interaction process j with an atom of the species i.

The attenuation of the beam 36 is not only dependent on the photon energy E thereof, but also upon the atomic species i of the slab 34 through which the beam propagates, as indicated by the use of the subscript and superscript i in equation (2). Certain atomic species are more efficient absorbers of high energy photons than others. For example, germanium is a much better absorber of high energy photons than silicon which in turn is a much better absorber than air. It is this very phenomena which is exploited by the apparatus 10 of FIG. 1 to measure the density and the mass fraction of each of the n molecular components of the soot boule 12.

Not only is the attenuation of the beam 36 of FIG. 2 dependent upon the atomic species i comprising the slab 34, but also, the physical state of the slab affects the attenuation of the beam. The dependence of the attenuation of the beam 36 on the physical state of the slab 34 can be eliminated by dividing equation (2) by the density .rho..sub.i of the slab. The term .mu..sub.i /.rho..sub.i is referred to as the mass attenuation coefficient.

The mass attenuation coefficient for a composite material comprised of n atomic components each of a species i is given by the sum of each individual product of the mass attenuation coefficient for each atomic species .mu..sub.i (E)/.rho..sub.i and the mass fraction x.sub.i thereof. Mathematically, this can be expressed by the relationship ##EQU2## where the term .mu.(E)/.rho. on the left-hand side of equation (3) represents the mass attenuation coefficient of the composite material. The subscript i on the right-hand side of equation (3) need not denote a particular atomic species but can also refer to a molecular species, such as germanium dioxide, silicon dioxide or air which is a bulk absorber of photons. When the subscript i refers to a molecular species, then the term .mu..sub.i (E)/.rho..sub.i refers to the mass attenuation coefficient of the particular molecular species i.

Referring back to FIG. 1, the unconsolidated boule 12 is usually comprised of three molecular components--air, amorphous silicon dioxide (SiO.sub.2) and amorphous germanium dioxide (GeO.sub.2), in decreasing order of concentration. With this knowledge, equation (3) can be specifically written for the soot boule 12 as follows ##EQU3## the terms .rho..sub.SiO.sbsb.2, .rho..sub.GeO.sbsb.2 and .rho..sub.air represent the individual densities of silicon dioxide, germanium dioxide and air, respectively, which can be readily obtained from a suitable reference text. The terms .mu..sub.SiO.sbsb.2 (E), .mu..sub.GeO.sbsb.2 (E) and .mu..sub.air (E) represent the linear attenuation coefficients for silicon dioxide, germanium dioxide and air, respectively. These quantities can be determined by means of a suitable calibration procedure described later. Note that since the boule 12 is inhomogeneous, the mass attenuation coefficient thereof is radially dependent so that density and each of the molecular mass fractions of the boule obtained from equation (4) will also be radially dependent.

Equation (4) is only valid when the photons bombarding the composite material are of a single energy E. However, as indicated previously, the beam 18 produced by the x-ray source 16 has a continuous energy spectrum E.sub.K. To account for the polychromatic nature of the beam 18, equation (4) must be modified as follows ##EQU4## where the term .mu.(r,E.sub.K) is given by ##EQU5## where E.sub.Khigh and E.sub.Klow represent the upper and lower energy bounds of the spectrum E.sub.K. The terms .mu..sub.SiO.sbsb.2 (E.sub.K), .mu..sub.GeO.sbsb.2 (E.sub.K) and .mu..sub.air (E.sub.K) are given by equations (7), (8) and (9), respectively. ##EQU6##

Determining the attenuation coefficient of a round, inhomogeneous object, such as the soot boule 12 of FIG. 1, is not as straightforward as determining the attenuation coefficient for the homogeneous slab 34 of FIG. 2. The attenuation coefficient of the slab 34 is easily determined from equation (1) by measuring the intensity of the beam upon entering and exiting the slab and the thickness L thereof. However, the attenuation coefficient of the boule 12 is dependent on its radius as well as the energy of the incident beam and thus, a more involved procedure is associated with its calculation.

Referring back to FIG. 1, the attenuation coefficient of the boule 12 can be determined in accordance with the attenuation of the beam 18 propagating therethrough. Such a determination can be facilitated by recognizing that the detector 22 actually senses the intensity of the beam 18 upon the propagation thereof between an entrance plane 37 and an exit plane 38 which are each parallel to the other and are each tangent to the boule 12 on opposite sides thereof.

If one temporarily assumes that the beam 18 is monochromatic, that is, each photon thereof has an energy E, then the intensity thereof exiting the plane 38 is given by ##EQU7## for the condition where the height t of the beam, that is the distance of the axis thereof from the center c of the boule 12, is less than the boule radius R. When the height t of the beam 18 is greater than the radius R of the boule 12, the intensity of the beam is given by the expression

As indicated by equation (11), air is an absorber of the beam 18, albeit a weak one. Thus, the intensity of the beam 18 is attenuated during the propagation thereof through space outside of the region between the entrance and exit planes 37 and 38, respectively. However, the attenuation of the intensity of the beam 18 suffered upon the propagation thereof through space outside of the region between the entrance and exit planes 37 and 38, respectively, is of little concern. The detector 22 effectively only measures the increased attenuation of the intensity of the beam 18 upon the propagation thereof inside the boule 12.

The term L(t) in equation (10) represents the path over which the beam 18 suffers attenuation inside the boule 12. By changing the variable integration from x, which represents the linear path between the entrance and exit planes 37 and 38, respectively, to r, which represents the distance along the radius R of the boule 12, the term .sub.L(t) .mu.(r,E)dx can be expressed as ##EQU8##

The ratio of the intensity of the beam 18 at the exit plane 38 for the condition t<R to the intensity of the beam for the condition t>R can be obtained by dividing equation (10) by equation (11). Taking the logarithm of the resulting equation and integrating the right-hand side by parts, one obtains ##EQU9## where the term .mu.(r)=.mu..sub.air for the condition t>r.

By differentiating equation (13) with respect to t and subsequently inverting the resultant integral, the attenuation coefficient .mu.(r,E) can be expressed as ##EQU10## From equation (14), the linear attenuation coefficient .mu.(r,E) can be determined in accordance with the measured ratio of the beam intensities i(E,t)/i(E,t>R).

Equation (14) is valid only when the beam 18 passing through the boule 12 is assumed to be monochromatic. However, the beam 18 produced by the x-ray tube 12 is polychromatic so that equation (12) must be modified as follows ##EQU11## where the term f(t) is given by ##EQU12##

The ratio i(E,t)/i(E,t>R) in equation (16) is the ratio of the intensity of those photons of the beam 18 which have energies within an incremental range .DELTA.E about an energy E within the bounds of the energy spectrum E.sub.K. Such a ratio is given by the count of photons associated with each separate one of the channels of the analyzer 30. The quantity f(t) in equation 16, which is characteristic of all the photons in the beam, is obtained by integrating the logorithm of the intensity ratio i(E,t)/i(E,t>r) of photons of any particular energy E within the energy limits of the beam 18. When the beam 18 is polychromatic, then a statement indicating that the attenuation coefficient .mu.(r) for photons within a certain energy group is established in accordance with the ratio of the intensity of photons leaving the boule to those entering the boule implies that the coefficient is calculated in accordance with equations (15) and (16). For a monochromatic beam of photons such a statement implies that coefficient is calculated in accordance with equation (14).

Once the attenuation coefficient of the boule 12 has been determined from equations (15) and (16) for the polychromatic beam 18, then equation (5) can be used to provide an expression characterizing the boule parameters .rho.(r), X.sub.SiO.sbsb.2 (r), X.sub.GeO.sbsb.2 (r) and X.sub.air (r) in terms of known quantities. However, since there are four unknown parameters, four equations relating these unknowns are necessary in order to find a unique solution therefor.

The law of conservation of mass may be utilized to generate two additional equations relating the unknown parameters of the boule 12 to known quantities. For a small volume element .DELTA.B of the boule 12 located a distance r from the center c thereof, the sum of the mass fractions of the constituents X.sub.SiO.sbsb.2 (r), X.sub.GeO.sbsb.2 (r) and x.sub.air (r) of the boule must equal unity. Mathematically, this may be expressed as

The partial volume of the individual constituents in the volume element .DELTA.B may be summed to yield the total volume of the element .DELTA.B. Mathematically, this may be expressed by the relationship ##EQU13##

A fourth equation relating the individual components of the boule 12 may be obtained by measuring the linear attenuation coefficient .mu.(r,E) thereof for photons within a second energy group, thereby yielding an additional equation of the same form as equation (5) but with different coefficients. The four equations may be written in a matrix equation ##EQU14## The arguments E.sub.a and E.sub.b of the attenuation coefficients in the first and second lines of equation (19) denote that the coefficients are valid for photons within each of two different energy groups E.sub.a and E.sub.b, respectively, within the spectrum E.sub.K of the beam 18. Each of the energy groups E.sub.a and E.sub.b within the energy spectra E.sub.K has an upper and lower energy bound E.sub.high and E.sub.low, respectively. By inverting equation (19), the column vector comprising the four unknown parameters can be obtained, thus providing a solution to the problem.

From a knowledge of the individual mass fractions X.sub.SiO.sbsb.2 (r), X.sub.GeO.sbsb.2 (r) and X.sub.air (r), the mole fraction of germanium in the boule 12 can be computed. The mole fraction of germanium in the boule 12 is a useful quantity because the refractive index of the boule has been found to be proportional thereto.

FIG. 3 is a block diagram depicting the individual steps associated with determining the density .rho.(r) and each of the mass fractions X.sub.SiO.sbsb.2 (r), X.sub.GeO.sbsb.2 (r) and X.sub.air (r) of the constituents of the boule 12. Initially, the apparatus 10 of FIG. 1 is calibrated (step 36) so as to determine the actual values for .mu..sub.air (E), .mu..sub.SiO.sbsb.2 (E) and .mu..sub.GeO.sbsb.2 (E) for the photons within the two energy groups E.sub.a and E.sub.b. Determination of the actual values of the attenuation coefficient of air, SiO.sub.2 and GeO.sub.2 is accomplished by directing the beam 18 through samples (not shown) of air, SiO.sub.2 and GeO.sub.2, respectively. The intensity of the photons in each of the energy groups E.sub.a and E.sub.b within the spectrum E.sub.K of beam 18 entering and leaving each sample, as well as the width of each sample, are measured. From this data, the attenuation coefficient of air, SiO.sub.2 and GeO.sub.2 can be calculated from equation (1) and equations (7 ), (8) and (9), respectively.

It is only necessary to calibrate the apparatus 10 of FIG. 1 once as long as each boule 12 profiled thereby is comprised of air, amorphous silicon dioxide and amorphous germanium dioxide. However, if there are additional constituent compounds in the boule 12, then the apparatus 10 must be recalibrated so as to determine the linear attenuation coefficient of each additional constituent for photons in each energy group.

After the apparatus 10 has been calibrated, then the boule 12 is positioned in the cutout 17 (FIG. 1) of the plate 14 (FIG. 1) and is subsequently scanned by beam 18 (step 38 of FIG. 3). As described, scanning of the boule 12 by the beam 18 is accomplished by moving the plate 14 in the direction of the double-ended arrow 27. Alternatively, the plate 14 could be kept stationary while the boule 12 was moved across the beam 18.

As the boule 12 is scanned, the intensity spectrum of the beam 18 at the exit plane 38 in FIG. 1 is recorded at each of a plurality of discrete heights t.sub.i (step 40). In practice, the intensity spectrum of the beam 18 is determined at each of 400 separate heights (i=400), each 0.01 cm apart during movement of the plate 14 in the direction indicated by the double-ended arrow 27.

As discussed previously, in order to solve equation (19) to obtain values for the density .rho.(r) and the individual mass fractions X.sub.SiO.sbsb.2 (r), X.sub.GeO.sbsb.2 (r) and X.sub.air (r) of the boule 12, the attenuation of the intensity of the beam 18 at each height t.sub.i must be obtained by measuring the attenuation of the photons with energies within each of two different energy groups E.sub.a and E.sub.b. Each of these energy groups may comprise more than one of the channels of the analyzer 30. As discussed previously, the beam 18 is polychromatic, that is to say, it contains photons of varying energy levels between E.sub.Khigh and E.sub.Klow. By measuring the attenuation of photons in each of two distinct energy groups E.sub.a and E.sub.b within the overall energy spectrum E.sub.K of beam 18, the boule 12 need only be scanned once. The pulse height scaler/analyzer 30 can readily be adjusted to count the number of photons propagating through the boule 12 which have an energy level within one of the two energy groups E.sub.a and E.sub.b within the energy spectra E.sub.K of the beam 18. From the data provided by the pulse height analyzer 30, two separate attenuation coefficients for the boule 12 can be obtained in accordance with equations (15) and (16) by substituting the upper and lower energy limits E.sub.High and E.sub.Low of each of energy group E.sub.a and E.sub.b for the upper and lower limits E.sub.Khigh and E.sub.Klow of the energy spectra E.sub.K of the beam 18.

Once the intensity of the beam 18 is measured in the above-described manner at each height t.sub.i, then a data value f(t.sub.i), representing the ratio of attenuated intensity of the beam to the incident intensity thereof for the photons within each of the two energy groups E.sub.a and E.sub.b is calculated (step 41). Mathematically, the data value f(t.sub.i) is determined in accordance with equation (20). ##EQU15## The particular values of the upper and lower energy boundary limits E.sub.1 and E.sub.2 in equation (20) correspond to the energy limits E.sub.High and E.sub.Low for the energy groups E.sub.a and E.sub.b.

Following step 41, a pair of linear attenuation coefficients .mu.(r,E.sub.a), .mu.(r,E.sub.b) for the boule 12 are computed at each height t.sub.i of the beam 18 for the photons within each of the two energy groups E.sub.a and E.sub.b (step 42). Calculation of the attenuation coefficients is accomplished by the computer 31 of FIG. 1 upon executing a program which is further described with respect to FIG. 4 four times, twice for each half of the boule 12. Once the two attenuation coefficients of the boule 12 have been calculated, then the density .rho.(r) and the mass fractions x.sub.air (r), x.sub.SiO.sbsb.2 (r) and x.sub.GeO.sbsb.2 (r) can be determined by solving equation (19) (step 44).

FIGS. 4 and 5, taken together, depict a flow chart representation of the program executed by the computer 31 of FIG. 1 to numerically solve equation (15) to determine each of the attenuation coefficients .mu.(r,E.sub.a) and .mu.(r,E.sub.b) of the boule 12 for the photons in each of the energy groups E.sub.a and E.sub.b of the beam 18 at each height t.sub.i thereof. Referring to FIG. 4, upon commencement of the program, the start step (step 50) is executed by the computer 31 to clear scratch pad memory locations therein which may not have been cleared previously. Thereafter, a numerical value for each of a pair of variables N and DEL is read from a memory location in the computer 31 (step 52). The variable N represents the number of intervals at which the intensity of the beam 18 is measured. In other words, the variable N corresponds to the number of discrete heights t.sub.i of the beam 18 for one half of the boule 12. The variable DEL represents the spacing between successive heights t.sub.i of the beam 18. Typically, N and DEL are assigned the values 200 and 0.01, respectively.

Next, a variable I is set initially equal to 1 (step 54). Thereafter, a stored data value, representing the Ith component of a vector quantity F(I), whose components each correspond to a separate one of the i values of f(t.sub.i) determined during step 38 of FIG. 3 for photons of the beam 18 within one of the two energy groups E.sub.a and E.sub.b, is read from a memory location in computer 31 (step 56). Following step 56, the Ith component of a vector quantity T(I) is established in accordance with equation

(step 58).

Thereafter, the value of I is augmented by 1 (step 60). The newly augmented value of I is then compared to the value of the quantity N+1 (step 62). For so long as the value of I does not exceed the value of the quantity N+1, then program execution branches back to step 56 and those following it to read in each of the N components of the vector F(I) and to compute the N separate components of the vector quantity T(I).

Once the value of I equals the value of the quantity N+1 then step 64 is executed whereupon the value of I is reset to 1. Next, a variable J is set equal to the quantity I-1 (step 66). The value of I is then tested (step 68) to see whether or not the current value thereof equals unity. Should the value of I be equal to 1, as will occur upon the initial execution of step 68, then the program execution branches to step 70 whereupon a vector quantity G(I), whose individual components each represent a separate one of the i values of the term df(t.sub.i)/dt, has the first component thereof set equal to zero to account for the fact that value of ##EQU16##

If, upon execution of step 68, the value of I is found not to be equal to 1, then the Ith component of the vector of G(I) is calculated (step 72) in accordance with equation

Those skilled in the art will recognize that the value assigned to the Ith component of the G(I) during step 72 is a numerical approximation of the value of the term df(t.sub.i)/dt.

Following either steps 70 or 72, the value of I is incremented by 1 (step 74). Thereafter, the newly augmented value of I is compared to the value of the quantity N+1 (step 76). For so long as the value of I is less than the value of the quantity N+1, then program execution branches back to step 66 and those following it to calculate the value for each of the N components of the vector G(I).

Should the current value of the variable I equa1 the value of quantity N+1, indicating that the N components of the vector G(I) have been calculated, then the computer 31 proceeds to numerically calculate an integral of the form ##EQU17## which corresponds to the integral of equation (15), using the trapezoidal rule.

In order to better understand the steps executed by the computer 31 to accomplish numerical calculation of the integral, some background discussion may prove helpful. The trapezoidal rule provides that an integral can be approximated by the sum of the individual discrete area elements under the curve thereof. Applying the trapezoidal rule to the integral on the left-hand side of equation (23), one obtains ##EQU18## where t.sub.i represents the Ith value of the height of the beam 18. Note that the first partial integral in equation (23) cannot be approximated by a trapezoidal area element since the term ##EQU19## becomes infinite for j=i.

In order to numerically calculate the integral on the right-hand side of equation (24), it is necessary to approximate the term g(t) by ##EQU20## for t.sub.i <t<t.sub.i+1. The constants A and B in equation (25) are given by ##EQU21## Substituting g(t) as given by equation (25) into the integral on the right-hand side of equation (24), one obtains an integral given as ##EQU22## which can be evaluated analytically with the solution being ##EQU23## Substituting this solution for the integral on the right-hand side of equation (24) one arrives at ##EQU24## where A and B are given by equations (26) and (27), respectively. Thus, equation 28 allows the integral in equation (23) to be calculated from the data previously computed during step 72.

In order to solve for the integral ##EQU25## in the manner set forth above, after step 76, the variable I is assigned the value of 2 (step 78) and the variable J is assigned the value of the quantity I+1 (step 80). Following step 80, step 82 is executed whereupon values for each of a pair of variables A and B are determined in accordance with equations (31) and (32)

As may be appreciated, the variables A and B in equations (31) and (32) represent numerical values of the arguments A and B, respectively, in equations (24) and (25), respectively, described earlier.

Also during step 82, a vector quantity DELMU(I) has its Ith component assigned a value in accordance with the equation (32)

As may be appreciated upon a comparison of equation (33) to equation (28), the Ith component of the vector DELMU(I) is a numeric approximation of the value of the first two terms of the left-hand side of equation (28) for the Ith value of the height t.sub.i of the beam 18.

Following step 82 of FIG. 4, the value of the variable I is then tested to determine whether or not the magnitude thereof is greater than the value of the quantity N-1 during step 84 shown in FIG. 5. Should the value of I be less than the value of the quantity N-1 indicating that the computer 31 has yet to complete calculation of the N separate components of the vector DELMU(I), then the variable J is assigned the value of the quantity I+1 (step 86) and thereafter, a variable D is assigned the value 1.0 (step 87).

Next, the value of the variable J is tested (step 88) by comparison thereof to the value of both the quantities I+1 and N. When the value of J is equal to the value of either the quantities I+1 or N, then, the variable D is reassigned the value 0.5 (step 90) before proceeding to step 92. Otherwise, when the value of J is unequal to the value of either of the quantities I+1 or N, program execution branches directly from step 88 to 92, omitting the execution of step 90.

Upon the execution of step 92, the value of the expression G(J)*(T(J)-T(I))*D/.sqroot.T(J)).sup.2 -(T(I)).sup.2 is added to the Ith component of the vector DELMU(I). Augmenting the Ith component of the vector DELMU(I) by this expression achieves a numerical approximation of the sum of the succeeding terms in equation (28) following the first two terms thereof. Since each component of the vector DELMU(I) is to represent the sum of the terms in equation (23) for the Ith value of the height t.sub.i of the beam 18 of FIG. 1, it is necessary to adjust the value of the variable D, as was accomplished during the step 90, to account for the factor of 1/2 found in that equation.

Following step 92, the value of the variable J is augmented by 1 (step 94) before comparing the current value of that variable to the value of the quantity N+1 (step 96). Should the value of J not equal the value of the quantity N+1, then program execution branches back to step 87 and those following it to accomplish a numerical approximation of the sum of all of the terms of equation (24).

Should the value of J be greater than the value quantity N+1 then, following step 96, program execution branches to step 98 whereupon the Ith component of the vector DELMU(I) is attenuated by 1/.pi. (step 98). Attenuating the value of the Ith component of the vector DELMU(I) during step 98 insures that the term 1/.pi. which is present in equation (15) is accounted for. The execution of step 98 also follows the execution of step 84 when the value of I equals the value of the quantity N-1 which occurs when all but the Nth component of the vector DELMU(I) has been calculated.

Following step 98, the Ith component of the vector DELMU(I) is then written into a memory location (not shown) within the computer 31 (step 100) Next, the value of the variable I is augmented by 1 (step 102). Thereafter, the value of the variable I is compared to the value of the quantity N-1 (step 104). If, during step 104, the value of I is unequal to the quantity N-1, then program execution branches back to step 80 and those following it so that all but the Nth component of the variable DELMU(I) are calculated.

Otherwise, when the value of I equals the quantity N-1, then after step 104, program execution branches to step 106 whereupon the Nth component of the vector DELMU(I) is set equal to 0 to account for the fact that, when the height t.sub.i of the beam 18 equals the radius R of the boule 12 of FIG. 1, the linear attenuation coefficient at that height is zero. Following step 106, the value of the Nth component of the vector DELMU(I) is written into a memory location in computer 31 (step 108) before program execution stops (step 110).

The above-described program is executed four times. The first two times, the program calculates the attenuation coefficient for each half of the boule 12 at each height t.sub.i for photons with energies within the group E.sub.a. Upon the third and fourth execution of the program, the attenuation coefficients for each half of the boule 12 at each height t.sub.i for photons with energies within the second group E.sub.b is determined.

The apparatus 10 of FIG. 1 can analyze inhomogeneous cylindrical objects such as soot boules having as many as n molecular components, provided that the scaler/analyzer 30 has at least n-1 channels and is capable of resolving the energy level of the photons of the beam 18 striking the detector 22 into one of n-1 discrete energy groups. From the data provided by the scaler/analyzer 30, the computer 31 can, by executing the program of FIG. 4 two times n-1 times, determine each of the n-1 attenuation coefficients of the object for the photons in each of the n-1 energy groups. Once the n-1 attenuation coefficients have been calculated, then the computer 31 can determine each of the mass fractions X.sub.1 (r), X.sub.2 (r), . . . X.sub.n (r) and the density .rho.(r) of the object using techniques similar to that used to solve equation (19).

It is to be understood that the various embodiments described herein are merely illustrative of the principles of the invention. Various modifications may be made thereto by persons skilled in the art which may embody the principles of the invention and fall within the spirit and scope thereof.