**METHODS AND APPARATUSES FOR RECONSTRUCTING ANGLE INFORMATION**

Schorr, David (3810 Frodo Cove, Austin, TX, 78739, US)

*; (IPC1-7): G01D/*

**G01S3/784**US5724047A | 1998-03-03 | |||

US4979221A | 1990-12-18 |

1. | A method for determining angular information from a target, comprising: (a) obtaining measured amplitude data from a plurality of sensors; (b) estimating an azimuth and elevation using the measured amplitude data; (c) obtaining calibrated amplitude data corresponding to the azimuth and elevation; (d) determining a residual error between the measured amplitude data and the calibrated amplitude data; (e) repeating steps (b) through (d) until the residual error is minimized; and (f) outputting the angular information as corresponding to the azimuth and elevation at which the residual error was minimized. |

2. | The method of claim 1, wherein obtaining measured amplitude data comprises obtaining measured amplitude data from seven sensors. |

3. | The method of claim 1, wherein one or more of the sensors comprise a CCD, thermal imaging array, or a quadrant detector. |

4. | The method of claim 1, wherein obtaining measured amplitude data comprises defocusing of the target. |

5. | The method of claim 1, further comprising determining a power centroid of the target using the measured amplitude data. |

6. | The method of claim 1, wherein determining the residual error comprises summing differences between the measured amplitude data and the calibrated amplitude data. |

7. | A computer program, comprising computer or machine readable program elements translatable for implementing the method of claim 1. |

8. | An apparatus for performing the method of claim 1. |

9. | A method for determining angular information from a target, comprising: (a) obtaining a plurality of sensors on an image plane; (b) defocusing the target or a point target on the image plane such that one or more of the sensors respond to a target pulse; (c) obtaining measured amplitude data from the plurality of sensors; (d) estimating an azimuth and elevation using the measured amplitude data; (e) obtaining calibrated amplitude data corresponding to the azimuth and elevation; (f) determining a residual error between the measured amplitude data and the calibrated amplitude data; (g) repeating steps (d) through (f) until the residual error is minimized; and (h) outputting the angular information as corresponding to the azimuth and elevation at which the residual error was minimized. |

10. | The method of claim 9, wherein obtaining measured amplitude data comprises obtaining measured amplitude data from seven sensors. |

11. | The method of claim 9, wherein one or more of the sensors comprise a CCD, thermal imaging array, or a quadrant detector. |

12. | The method of claim 9, further comprising determining a power centroid of the target using the measured amplitude data. |

13. | The method of claim 9, wherein determining the residual error comprises summing differences between the measured amplitude data and the calibrated amplitude data. |

14. | A computer program, comprising computer or machine readable program elements translatable for implementing the method of claim 9. |

15. | An apparatus for performing the method of claim 9. |

16. | A system for determining angular information from a target, comprising: a plurality of sensors on an image plane; and a computer configured to: acquire measured amplitude data from the plurality of sensors; estimate an azimuth and elevation using the measured amplitude data; access calibrated amplitude data corresponding to the azimuth and elevation; determine a residual error between the measured amplitude data and the calibrated amplitude data; and output the angular information as corresponding to the azimuth and elevation at which the residual error is minimized. |

60/334,386 entitled, "Methods And Apparatuses for Reconstructing Angle Information"by David Schorr. That application is incorporated herein by reference in its entirety.

1. Field of the Invention The present invention relates generally to the fields of position sensors and the analysis of positional data. More particularly, it concerns methods and apparatuses for efficiently reconstructing angle information by analyzing data from a set of independent sensors. Even more particularly, it concerns the reconstruction of angle information by comparing a data set from independent sensors with a calibration data set and minimizing an error function between those sets of data to reconstruct angle information.

2. Description of Related Art The ability to gather data from one or more sensors and translate that data into meaningful positional information is important in a wide variety of disciplines. For instance, tracking objects is essential for commercial aircraft controllers and an assortment of military applications.

Over the years, sundry techniques aimed at determining positional information, such as angular information, from sensor readings have been developed. Although each has shown at least a degree of success in its respective application, room for significant improvement remains.

U. S. Patent No. 5,724, 047, which is hereby incorporated by reference, involves a precision direction finding system for making precision angle of arrival estimates for a signal received through two antenna elements separated in space. Phase interferometry is used to determine a precise angle of arrival, with multiple ambiguities due to the periodic nature of the

phase difference related to geometric angle. The interferometric ambiguities are resolved using the time difference of arrival (TDOA) of the signal at the two antenna elements. TDOA is measured using leading edge envelope detection for simple pulsed signals, and predetection correlation for phase and frequency modulated signals. Although such phase interferometry determinations may exhibit certain advantages over other conventional methods, shortcomings remain due to, at least in part, the relatively complex nature of the computations required to implement the methodology. Further, reliance upon two antenna elements, rather than a plurality of independent sensors, may inhibit at least some flexibility of such a system.

U. S. Patent No. 6,104, 345, which is hereby incorporated by reference, involves a method for the direction of arrival (DOA) tracking of at least one source, along a single azimuth axis or along both azimuth and elevation axes. The method includes the steps of selecting all the high peaks from a DOA function as potential track points, converting the potential track points into a plurality of tracks and selecting a true track from the plurality of tracks. This methodology also suffers from some shortcomings; for example, it may be desirable to avoid the particular computational complexity involved with selecting a true track from the plurality of potential track points.

U. S. Patent No. 4,562, 439, which is hereby incorporated by reference, involves an imaging radar seeker for producing two-dimensional images of a target, which can be mounted on a missile or other moving body, such as an automobile. A computer directs the seeker to operate sequentially in searching, tracking, and imaging modes. In the searching mode, a combination of circumferential rotation of an antenna of the seeker and frequency scanning of electromagnetic energy fed to the antenna enables the seeker to search for its target over a conical field-of-view or a wider, peripheral belt field-of-view. In the imaging mode, circumferential rotation of an antenna is stopped, and the tilt angle of a linear array of the antenna is stepped or continuously moved to compensate for radial movement of a radiated beam caused by frequency stepping imparted by a frequency synthesizer. This keeps the beam fixed in space and centered on the target. Inverse synthetic aperture imaging is used to create a two- dimensional image of the target wherein the first dimension (range) is obtained by performing inverse Fourier transforms on the echo signals, and the second orthogonal dimension (cross- range or doppler frequency) is obtained by performing Fourier transforms. The array can be a

linear array of E-plane stacked linear waveguide antenna elements operating in either the traveling wave mode or the standing wave mode. Such methodology too exhibits shortcomings at least due to its complexity, lack of flexibility, and its inability to more efficiently perform angle reconstruction determinations.

U. S. Patent No. 5,818, 393, which is hereby incorporated by reference, involves a fixed body wide field-of-view conformal antenna array suitable for broadband precision direction finding on missile platforms. The array is configured as multiple sub-arrays of spiral antennas that cover particular regions within the desired field-of-view of the entire array. A lower cost, more reliable and more accurate direction finding solution for missile needs is provided, primarily by the elimination of conventional radomes and antenna gimbal structures. The array can be configured to include multi-mode sensors. Although useful in this elimination of structures, this methodology exhibits room for significant improvement given its inability to more efficiently reconstruct angular information.

U. S. Patent No. 6,313, 794, which is hereby incorporated by reference, uses feedback from RF carrier frequency measurements to disassociate the emitter angle-of-arrival component in the ambiguous phase measurement from the initially unknown phase measurement integer ambiguities. It then resolves the ambiguities and obtains the correct emitter angle-of-arrival (AOA). This is accomplished by converting the actual interferometer baselines on which the unassociated pulse phase measurements were made at different emitter frequencies to a baseline set for a single-frequency equivalent interferometer array. This methodology, as with the methodologies described above, suffers from issues relating to complexity, flexibility, and the inability to more quickly and efficiently reconstruct high-resolution angular information.

U. S. Patent No. 5,657, 251, which is hereby incorporated by reference, involves a computer-implemented process for processing incoming target data from a focal plane or scanning radar to accomplish multiple Target Tracking. Inputs are pixel plane coordinates and intensity of target blips. The Intelligent Target Tracking Processor (ITTP) employs an optimal target tracking algorithm. An optimal observation-to-track assignment exists when all target blips in a new frame of target data are matched up with nearby tracks, such that the sum of all the distances from each target blip to its assigned track is minimized. An expert system is used to

control overall processing flow and provide efficient allocation of computing resources. Target blips without near neighbors are allowed to go directly to a real track table of established tracks, if their coordinates match-up with projected tracking gates. Otherwise, target blips are tested sequentially against two-frame, three-frame, and four-or higher-frame discriminants, to reject blips not belonging to established tracks. The ITTP can partition the pixel plane into"bite size partitions", each with a manageable number of target blips, which it handles sequentially. The ITTP is designed to handle hundreds or thousands of targets as a stand-alone processor as is required in space object tracking or military scenarios. The expert system maintains an optimum balance between correlating on existing tracks and discriminating against impossible tracks. A total of 26 different metric and radiometric target tracking discriminants are employed. The ITTP is a dynamically and optimally configured set of general purpose parallel processors.

Although useful for trajectory tracking for applications such as air-traffic controlling, this methodology nevertheless includes shortcomings; for instance, it does not allow for the reconstruction of high quality, high resolution angle information using a plurality of independent sensors or detectors on an image plane.

Problems enumerated above are not intended to be exhaustive, but rather are among many that tend to impair the effectiveness of previously known techniques concerning positional tracking and the determination of angle information. Other noteworthy problems may also exist; however, those presented above should be sufficient to demonstrate that methodology appearing in the art have not been altogether satisfactory. In particular, existing techniques do not adequately allow for the determination of angle information with improved resolution from independent sensors or detectors on an image plane.

SUMMARY OF THE INVENTION Shortcomings listed above are reduced or eliminated by the techniques disclosed herein.

These techniques are applicable to a vast number of applications, including general targeting and tracking applications.

In one respect, the invention is a method for determining angular information from a target. Measured amplitude data is obtained from a plurality of sensors. An azimuth and

elevation are estimated using the measured amplitude data. Calibrated amplitude data corresponding to the azimuth and elevation is obtained. A residual error between the measured amplitude data and the calibrated amplitude data is determined. Aspects of this process are performed iteratively until the residual error is minimized. The angular information, which corresponds to the azimuth and elevation at which the residual error was minimized, is then output.

In other respects, the invention is a computer program, comprising computer or machine readable program elements translatable for implementing the method described above. Further the invention is an apparatus for performing such a method.

In another respect, the invention is a method for determining angular information from a target. A plurality of sensors on an image plane are obtained. The target or a point target is defocused on the image plane such that one or more of the sensors respond to a target pulse.

Measured amplitude data from the plurality of sensors is obtained. An azimuth and elevation are estimated using the measured amplitude data. Calibrated amplitude data corresponding to the azimuth and elevation is obtained. A residual error between the measured amplitude data and the calibrated amplitude data is determined. Aspects of this process are performed iteratively until the residual error is minimized. The angular information, which corresponds to the azimuth and elevation at which the residual error was minimized, is then output.

In other respects, the invention is a computer program, comprising computer or machine readable program elements translatable for implementing the method described above. Further the invention is an apparatus for performing such a method.

In another respect, the invention is a system for determining angular information from a target. The system includes a plurality of sensors and a computer. The computer is configured to: (a) acquire measured amplitude data from the plurality of sensors; (b) estimate an azimuth and elevation using the measured amplitude data; (c) access calibrated amplitude data corresponding to the azimuth and elevation; (d) determine a residual error between the measured amplitude data and the calibrated amplitude data; and (e) output the angular information as corresponding to the azimuth and elevation at which the residual error is minimized.

Other features and associated advantages will become apparent with reference to the following detailed description of specific embodiments in connection with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS The following drawings form part of the present specification and are included to further demonstrate certain aspects of the present invention. The invention may be better understood by reference to one or more of these drawings in combination with the detailed description of specific embodiments presented herein.

FIG. 1 is a flowchart showing techniques for reconstructing angle information in accordance with embodiments of the present disclosure.

FIGS. 2A and 2B show exemplary data involving angle information in accordance with embodiments of the present disclosure.

DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS The present disclosure describes apparatuses and methods for efficiently determining high-quality, high-resolution angle information from several independent sensors on an image plane. Applications for this technology are vast. For instance, one may use the techniques disclosed herein in a variety of military applications. Similarly, one may find applications in aircraft traffic controlling, general target seeking and/or tracking, and any other application where angular information is required. With the benefit of the present disclosure, those having skill in the art will comprehend that the techniques disclosed herein may be modified and applied to a number of additional, different applications. The present disclosure and claims attached hereto cover all such modifications that fall within the scope and spirit of this disclosure.

In one embodiment, the techniques for determining angle information may take the form of what the inventor has coined an"angle reconstruction algorithm. "The angle reconstruction

algorithm is one method for improving the resolution and quality of angle information from relative amplitude measurement information given a set of independent sensors or detectors on an image plane. These sensors may be part of a CCD, thermal imaging array, a quadrant detector, or any other device suitable for gathering data relevant to angular information. The angle reconstruction process may involve producing a"match"between predicted or calibrated relative amplitude responses and measured relative amplitude responses. In one embodiment, this match may involve minimizing a residual error function, although those having skill in the art will understand that other mathematical techniques may be employed to determine a correspondence between the predicted or calibrated amplitudes and the measured responses.

In one embodiment, a LASER seeker field-of-view (FOV) may be broken up into a plurality of distinct pixels. In one embodiment, seven pixels may be used, although more or fewer pixels may be utilized depending upon preference and/or the specific application in question. A point source may be defocused on the image plane such that a plurality of sensors respond to a given target pulse. In the embodiment using seven pixels, for instance, a point source may be defocused such that at least 3 of 7 independent detectors respond to a target pulse.

The resulting amplitudes from each of the detectors (7 in the embodiment discussed above) for a given received pulse may then be compared with expected amplitudes from either predictions or calibration data for each location in the FOV. The location in the FOV that produces a match or correspondence, which may take the form of exhibiting the least residual error between the predicted response magnitudes and the actual response amplitude measurements, is then determined to be the target's location in the field of view. Hence, angle (and position) information may be efficiently obtained.

As will be understood by those having skill in the art, one may determine an appropriate correspondence between measured amplitudes and prediction or calibration sets of data in many ways, including the determination of a minimum residual error. In one embodiment, a rather robust and computationally-intensive approach is to accomplish a grid search of a region of interest. The residual between the amplitude response prediction and the actual measured amplitude response may be compared for all grid points. The grid point resulting in the lowest

residual is assumed to be the best estimate of the targets location in the FOV. Refinement of the grid search interval may be accomplished until a particular desired level of accuracy is achieved.

Other computationally-efficient methods of finding the location in the FOV with the lowest residual may be readily implemented according to programming and mathematical methods known in the art. Different methods may, of course, amount to a significant reduction in the computational effort required to find a target location in the field of view.

Angle reconstruction techniques of this disclosure applied to an imaging array or other sensors provide for angle resolution much greater than would be derived based on the estimated angular extent of a pixel on the image plane. Intentional defocusing of the target (when angular extent is less than 1 pixel) or point target allows for relative amplitude data to be collected on multiple adjacent pixels. The relative amplitude data from each pixel can then be utilized to determine the power centroid of an extended target or a point source in the FOV.

These techniques allow very noisy non-linear raw sensor data to be interpreted into angle of arrival (AOA) with the aid of a database and an appropriate calibration set of data or data forming the basis of a prediction set of data. Noise corruption of the measured signal levels often prevents a perfect match to angle of arrival to be made; however, extremely good estimates of angle of arrival may be readily made with an angle reconstruction algorithm as described herein. No loss of bandwidth need be realized since previous data does not have to be utilized in determining the angle of arrival output.

Apparatuses suitable for implementing aspects of the methods described herein may include equipment known in the art such as sensors, detectors, seekers, etc. Additionally, the apparatuses may include one or more computers to analyze the sensor data, in accordance with the techniques described herein via corresponding software. The one or more computers may be networked and may be in contact with one or more local or remote databases or storage networks that store calibration or prediction data sets relating to the sensors or arrays.

The following examples are included to demonstrate specific embodiments of this disclosure. It should be appreciated by those of skill in the art that the techniques disclosed in

the example which follow represent techniques discovered by the inventor to function well in the practice of the invention, and thus can be considered to constitute specific modes for its practice.

However, those of skill in the art should, in light of the present disclosure, appreciate that many changes can be made in the specific embodiments which are disclosed and still obtain a like or similar result without departing from the spirit and scope of the invention.

EXAMPLE 1 Turning to FIG. 1, there is shown a flowchart illustrating one embodiment of a suitable angle reconstruction algorithm in which 7 detectors are used. As the flowchart shows, amplitudes from the seven detectors are first measured and then normalized, if necessary. This information is used to form an initial estimate of both azimuth (AZ) and elevation (El) information. Using a database or another appropriate storage device or network, one may then extract amplitude information corresponding to the initially-estimated AZ and El from calibration or prediction data. This amplitude information may be normalized, if necessary.

With the measured and calibration amplitude information, both corresponding to the initially- estimated AZ and El, one may then determine a residual error function (or determine a correspondence between the data sets using an alternative method).

In FIG. 1, the amplitude residual is determined by considering the difference between the measured amplitude and the amplitude from the calibration set, for each and every detector (detectors A through G). Specifically, one may sum the squares of the differences of those amplitudes, as shown, to arrive at the residual value. In FIG. 1, Ameas means the measured amplitude from detector A, while Ast means the amplitude at the estimated AZ and El from the calibration set.

If this residual value is not the minimum value, the initial estimate of AZ and El is updated, as shown. With this new AZ and El estimate, amplitude information may again be extracted from the calibration data set-this time, the calibration information will correspond to newly-estimated AZ and El. Again, normalization of the calibration amplitude information may be required. Once again, a residual error function may be obtained by considering the

differences between the measured and calibration amplitudes. And, once again, if the residual is not minimized, the estimate of AZ and El may be updated, thereby repeating the process.

However, if the residual is minimized, one may output the AZ and El that accountable for the minimization. It is this AZ and El that represents the angle information of the target or object being studied.

FIGS. 2A and 2B show exemplary amplitude data from seven detectors.

This example and the disclosure herein shows that, in general, target location and angular information may be triangulated by comparison to a look up table (which may take the form of the calibration data set or prediction data discussed above). The look up table may be several dimensions, including 3. The look up table may be built-up via calibration steps. Specifically, one may position a target sequentially at various, known locations throughout a grid. At each location, one may take sensor readings from a plurality of sensors to obtain relative light amplitudes. Taking multiple readings of this sort builds up a library of amplitudes as a function of target locations. This library may then be used to compare with measured amplitudes in the fashion described herein (for example, by determining a residual error function) to quickly determine the angular and positional information of a target based on measured amplitude data from those same sensors.

EXAMPLE 2 Shown in this example is source code, written in FORTRAN, that is suitable for carrying out steps described herein. This source code is exemplary only and does not limit the scope of the claims appended hereto. It simply represents one specific embodiment for carrying out aspects of this disclosure and is included for the convenience of the reader in this regard. Those having skill in the art, with the benefit of this disclosure, will recognize that a wide variety of computational techniques and different types of corresponding source code and software may be used in implementing the concepts described herein.

Program AngleRecorof C*********************************************************** ************** C This code synthesises ideal detector amplitudes from AZ EL seeker angle C and corrupts the amplitudes with predicted noise. C The corrupted detector amplitudes are then utilized to reconstruct the C original seeker angles through an angle reconstruction algorithm. C The effects of the noise on the detector amplitudes in studied. C Need to check optical gain profile in data statements and optical efficiency. C*********************************************************** ************** USE MSIMSL USE MSFLIB C USE PORTLIB IMPLICIT DOUBLE PRECISION (A-H, O-Z) Parameter (NDET = 7) Dimension Signal (NDET) Dimension AZ (1000), EL (1000) Double Precision SNR, SNREST C C*********************************************************** ************** C All of these files are not necessarily used every time. These are C All of these files are output files. OPEN (UNIT=7, FILE='AmplitudeNoise. xls') Write (7, *)'TargetRange AZtruth ELtruth SNR AZnoise ELnoise' OPEN (UNIT=8, FILE='CalMatrix20X20. txt') OPEN (UNIT=9, FILE='Errormap. xls') OPEN (UNIT=10, FILE='Anglemeas. xls') OPEN (UNIT=11, FILE='Angleres. xls') C C*********************************************************** ************** c Nuber of Monte Carlo noise cases for shotgun plot and noise statistics Nmontenoise = 100 C*********************************************************** ************** Do 6000 Iseekercase=1, 3, 1 If (Iseekercase. EQ. 1) Then Aztruth = 0. 0 Eltruth = 0. 0 OPEN (UNIT=12, FILE='AnglestatszO-0. xls') Write (12, *)'TargetRange, AZtruth, ELtruth, Nmontenoise, AZmean, &ELmean, AZrms@0/0, ELrms@0/0' ElseIf (Iseekercase. EQ. 2) Then Aztruth = 5. 00 Eltruth = 0. 00 OPEN (UNIT=12, FILE='Anglestatsz5-0. xls') Write (12, *)'TargetRange, AZtruth, ELtruth, Nmontenoise, AZmean, &ELmean, AZrms@5/0, ELrms@5/0'

ElseIf (Iseekercase. EQ. 3) Then Aztruth = 10. 0 Eltruth = 0. 0 OPEN (UNIT=12, FILE='AnglestatszlO-O. xls') Write (12, *) 'TargetRange, AZtruth, ELtruth, Nmontenoise, AZmean, &ELmean, AZrms@10/0, ELrms@10/0' ElseIf (Iseekercase. EQ. 4) Then Aztruth = 4. 00 Eltruth = 7. 00 ElseIf (Iseekercase. EQ. 5) Then Aztruth = 8. 50 Eltruth =-10. 75 ElseIf (Iseekercase. EQ. 6) Then Aztruth = 0. 50 Eltruth =-13. 00 ElseIf (Iseekercase. EQ. 7) Then Aztruth =-12. 75 Eltruth = 12. 75 ElseIf (Iseekercase. EQ. 8) Then Aztruth = 16. 9 Eltruth = 10. 25 ElseIf (Iseekercase. EQ. 9) Then Aztruth =-2. 25 Eltruth =-5. 00 ElseIf (Iseekercase. EQ. 10) Then Aztruth =-10. 75 Eltruth =-3. 00 ElseIf (Iseekercase. EQ. lu) Then Aztruth = 13. 25 Eltruth = 11. 75 Endif c AZtruth=float (Iseekercase)/10. 0 c ELtruth=0. 0 c*********************************************************** *********** c******************* This is the MonteCarlo Noise Loop **************** Do 1000 Irange=100, 6000, 100 Do 2000 Inoise=1, Nmontenoise TargetRange = dmaxl (Float (Irange), 50. 0) DesignatorRange = 6000. 0 Call DetectorAmps (TargetRange, DesignatorRange, &AZtruth, ELtruth, SNREST, Signal) Call SeekerNoise (Signal, AZnoise, ELnoise) c Write (*, 1500) TargetRange, AZtruth, ELtruth, SNREST, AZnoise, ELnoise c Write (12, 1500) TargetRange, AZtruth, ELtruth, SNREST, AZnoise, ELnoise Write (7, 1500) TargetRange, AZtruth, ELtruth, SNREST, AZnoise, ELnoise 1500 Format (9F10. 3) c Save results to array for post processing statistics AZ (Inoise) =AZnoise EL (Inoise) =ELnoise 2000 Continue c*************** End of Monte Carlo Noise Loop ********************** c****************Calculate Summary Statistics Data ******************

Call Anglestats (Inoise-l, AZ, EL, AZmean, ELmean, AZrms, ELrms) Write (*, 1500) TargetRange, AZtruth, ELtruth, &dFloat (Inoise-1), AZmean, ELmean, AZrms, ELrms Write (12, 1500) TargetRange, AZtruth, ELtruth, &dFloat (Inoise-1), AZmean, ELmean, AZrms, ELrms 1000 Continue c*************** End of AZ, EL case study **************************** C*********************************************************** ********* close (12) 6000 Continue Stop END C*********************************************************** ************** C END of trash code used to excercize the subroutine for seeker noise. C*********************************************************** ************** C*********************************************************** ************** Subroutine AngleStats (N, AZ, EL, AZmean, ELmean, AZrms, ELrms) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Dimension AZ (1000), EL (1000) c calculate mean of n data points AZsum=0. 0 ELsum=0. 0 AZres=0. 0 ELres=0. 0 Do 1000 i=l, N AZsum=AZsum + AZ (i) ELsum=ELsum + EL (i) 1000 continue AZmean=AZsum/n ELmean=ELsum/n Do 2000 i=l, N AZres=AZres + (AZ (i)-AZmean) **2 ELres=ELres + (EL (i)-ELmean) **2 2000 continue AZrms=SQRT (AZres/Dfloat (N-1)) ELrms=SQRT (ELres/Dfloat (N-1)) return end C*********************************************************** ************** return end Subroutine DetectorAmps (TargetRange, DesignatorRange, &AZtruth, ELtruth, SNREST, Signal) C*********************************************************** ************** C This code generates ideal detector signal levels, then corrupts those signals C with the anticipated noise levels from the APD and C assumed internal SNR limits. C C As a note of caution, no approximation to scintillation is included C although this is felt to be a dominat influence at C intermediate ranges (i. e. from 500 meters to 3 km) C C Inputs : C TargetRange = range from Seeket to Target in meters C DesignatorRange = range from Designator to target in meters

C AZtruth = true AZ line of sight C ELtruth = true EL line of sight C Outputs : C SNREST = predicted signal to noise ratio C Signal = predicted normalized detector signal levels C corrupted with noise C C*********************************************************** ************** IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Parameters C Parameter (NAZ = 81) Parameter (NEL = 81) Parameter (NDET = 7) C Dimension Signal (NDET) Dimension OPT GAIN (NDET) Double Precision TotAng (NDET) Double Precision Theta Double Precision whatimeisit Logical Firstcall Dimension Anglens (3, NDET) COMMON/Gainlens/Nangles, Angles (200), OptTrans (200) COMMON/Caldata/AZtbl (NAZ), ELtbl (NEL), Tgain (NAZ, NEL, NDET) C CONSTANTS Parameter (PI=3. 1415925) Parameter (RTD = 180. OdO/PI) Parameter (DTR = PI/180. OdO) C Optical Gain vs Total Angle for individual Lens (30 data points) Data Nangles/30/ Data Angles/0. 0, 0. 5, 1. 0, 1. 5, 2. 0, &2. 5, 3. 0, 3. 5, 4. 0, 4. 5, &5. 0, 5. 5, 6. 0, 6. 5, 7. 0, &7. 5, 8. 0, 8. 5, 9. 0, 9. 5, &10. 0, 10. 5, 11. 0, 11. 5, 12. 0, &12. 5, 13. 0, 13. 5, 14. 0, 180. 0, 170*0/ c This data is from Zemax Data OptTrans/0. 146, 0. 147, 0. 154, 0. 166, 0. 168, &0. 170, 0. 156, 0. 150, 0. 141, 0. 137, &0. 128, 0. 120, 0. 113, 0. 107, 0. 102, &0. 090, 0. 087, 0. 077, 0. 072, 0. 061, &0. 055, 0. 046, 0. 039, 0. 030, 0. 025, &0. 016, 0. 011, 0. 005, 0. 000, 0. 000, 170*0/ c This data is a slice from measure calibration used for outdoor testing at TA-3 c Data OptTrans/0. 150, 0. 146, 0. 141, 0. 140, 0. 133, c & 0. 123, 0. 115, 0. 107, 0. 100, 0. 090,

c & 0. 082, 0. 073, 0. 065, 0. 057, 0. 049, c & 0. 041, 0. 035, 0. 028, 0. 023, 0. 019, c & 0. 015, 0. 011, 0. 009, 0. 007, 0. 005, c & 0. 003, 0. 002, 0. 001, 0. 001, 0. 000, 170*0/ Data Firstcall/. True./ C C*********************************************************** ********** c c C READ OPTICAL gain versus angle data. C READ OPTICAL gain versus angle data. C This is the only required input file. C If (Firstcall) then If (Firstcall) then Firstcall=. False. C Point of each of 7 lenses C ******Lens Pattern************************ C 4 3 C 5 1 2 C 6 7 C******************************************* Anglens (1, 1) = 0. 0 Anglens (2, 1) = 0. 0 Anglens (1, 2) = 14. 00 Anglens (2, 2) = 0. 0 Anglens (1, 3) = 7. 0 Anglens (2, 3) = 12. 12 Anglens (1, 4) =-7. 00 Anglens (2, 4) = 12. 12 Anglens (1, 5) =-14. 00 Anglens (2, 5) = 0. 00 Anglens (1, 6) =-7. 00 Anglens (2, 6) =-12. 12 Anglens (1, 7) =7. 0 Anglens (2, 7) =-12. 12 C*********************************************************** ************** C*********************************************************** ************** C Generate Analytical Calibration tables from C Analytical Optical Gain Profiles on 0. 5 degree intervals C This will be replaced with measured calibration data when available. Do 1000 IEL= NEL, 1,-1 Do 2000 IAZ= NAZ, 1,-1 AZ = (dfloat (IAZ) *0. 5dO-20. 5dO) EL = (dfloat (IEL) *0. 5dO-20. 5dO) Aztbl (IAZ) = AZ Eltbl (IEL) = EL

C Determine total angle in radians to each lens Do 3000 Idet = 1, 7 TotAng (Idet) =sqrt ( (Anglens (l, Idet)-AZ) **2+ (Anglens (2, Idet)-EL) **2) Tgain (IAZ, IEL, Idet) = OPGain3 (TotAng (Idet))/OPGain3 (0. 0) 3000 Continue C Write out calibration file for MatLab Plot utilities. C Write (8, 1600) AZ, EL, 1. 00 * Tgain (IAZ, IEL, 1), Tgain (IAZ, IEL, 2), C & Tgain (IAZ, IEL, 3), Tgain (IAZ, IEL, 4), C & Tgain (IAZ, IEL, 5), Tgain (IAZ, IEL, 6), C & Tgain (IAZ, IEL, 7) 2000 Continue 1000 Continue 1600 Format (17 (F10. 4,','), F10. 4) Write (*, *)'Calibration Table Generation Complete.' C C C Aperture of individual optic in cm^2 C Aperture of individual optic in cm^2 APERTURE = 0. 78/10000. 0 C Total FOV of individual lens in radians C (assumed half power full width) TLFOV = 14. 0*DTR C Fraction of Photons going thru aperture than hit detector C Includes numerical aperture effects but excludes spot size vs ferrul size. C Optical Gain Profile updated to include overall efficiency. Opteff= 1. 00 C APD Gain Gain = 20 C Temperature (in Kelvin) Temp = 300 C LASER Pulse Width (seconds) PlWid = 20e-9 C BW of filter on amplifier (Hz) BW = 6. Oe6 BWnoise = BW*PI Tau=l/ (2*pi*BW) ; BWfactor= (1-exp (-PlWid/Tau)) C Power of the laser (Watts) POWER LASER = 5. Oe6 C Atmospheric extintion coefficient (km) C EXT COEF =-0. 09135 Visibility = 23. 5d0 EXT_COEF=-(10. OdO** (-0. 136dO+ (1. 16dO*loglO (3. 912dO/Visibility)))) C Target Reflectivity REFLTAR = 0. 1 C Reflectivity of background (i. e. terrain) REFLBACK = 0. 1

C Number of detectors per aperture Numdet = 1 C Number of apertures in the reciever NUMAPERTRECV = 3 NUMAPERTBACK = 4 C INCIDENCT SPECTRAL FLUX WATTS/ (M-2 MICRON) 1. 06 microns SPEC FLUX = 500. C OPTICAL FILTER BANDWIDTH IN MICRONS FILBW =. 010 C Range in meters C Range = 6000 Endif C End of calculations done for first subroutine call only C C C*********************************************************** ************** C Note reference to optical gain at 0 radians offset. C Note assumes designator fixed at 6km while target range is variable. EXTINCT = EXP (EXTCOEF* (TargetRange+DesignatorRange)/1000) POWER = POWER LASER*EXTINCT*REFLTAR*APERTURE* &FLOAT (NUMAPERTRECV) *OPTEFF*OPGAIN3 (0. 0)/ &(3. 14159*TargetRange**2*FLOAT (NUMDET)) BACKPOW = (TLFOV) **2*SPECFLUX*FILBW*APERTURE*NUMAPERTBACK* &REFLBACK*OPTEFF*OPGAIN3 (0. 0)/ (NUMDET*8.) Total Power = POWER + BACKPOW Call DETECT (TotalPower, GAIN, TEMP, BWnoise, ZNEP) SNRlimit=1000. 0 zneplim=sqrt (znep**2+ (1/SNRlimit*BWFactor*POWER) **2) SNR = power*BWFactor/zneplim C C C Get ideal"truth amplitudes"from optical model for truth Az, El location C before adding noise to ideal amplitudes. C Determine total angle in DEGREES to each lens. TotAng (1) =sqrt ((Anglens (1, 1)-AZtruth) **2 &+ (Anglens (2, 1)-ELtruth) **2) TotAng (2) =sqrt ( (Anglens (1, 2)-AZtruth) **2 &+ (Anglens (2, 2)-ELtruth) **2) TotAng (3) =sqrt ( (Anglens (1, 3)-AZtruth) **2 &+ (Anglens (2, 3)-ELtruth) **2) TotAng (4) =sqrt ( (Anglens (1, 4)-AZtruth) **2 &+ (Anglens (2, 4)-ELtruth) **2) TotAng (5) =sqrt ( (Anglens (1, 5)-AZtruth) **2 &+ (Anglens (2, 5)-ELtruth) **2) TotAng (6) =sqrt ( (Anglens (1, 6)-AZtruth) **2 &+ (Anglens (2, 6)-ELtruth) **2) TotAng (7) =sqrt ( (Anglens (1, 7)-AZtruth) **2 &+ (Anglens (2, 7)-ELtruth) **2) C Determine Optical Gain for each lens Opt_gain (1) = OPGain3 (TotAng (1))/OPGain3 (0. 0) Ope gain (2) = OPGain3 (TotAng (2))/OPGain3 (0. 0) Ope gain (3) = OPGain3 (TotAng (3))/OPGain3 (0. 0)

Opt gain (4) = OPGain3 (TotAng (4))/OPGain3 (0. 0) Opt_gain (5) = OPGain3 (TotAng (5))/OPGain3 (0. 0) Opt gain (6) = OPGain3 (TotAng (6))/OPGain3 (0. 0) Opt_gain (7) = OPGain3 (TotAng (7))/OPGain3 (0. 0) C Add noise to signal Signal (1) = power*OPT GAIN (1) + znep*dnorin (zrand ()) Signal (2) = power*OPT GAIN (2) + znep*dnorin (zrand ()) Signal (3) = power*OPT GAIN (3) + znep*dnorin (zrand ()) Signal (4) = power*OPT GAIN (4) + znep*dnorin (zrand ()) Signal (5) = power*OPT GAIN (5) + znep*dnorin (zrand ()) Signal (6) = power*OPT GAIN (6) + znep*dnorin (zrand ()) Signal (7) = power*OPTGAIN (7) + znep*dnorin (zrand ()) C Normalize all measured amplitudes to maximum Signalmax = MAX (Signal (l), Signal (2), Signal (3), &Signal (4), Signal (5), Signal (6), &Signal (7)) Signal (1) = Signal (1)/Signalmax Signal (2) = Signal (2)/Signalmax Signal (3) = Signal (3)/Signalmax Signal (4) = Signal (4)/Signalmax Signal (5) = Signal (5)/Signalmax Signal (6) = Signal (6)/Signalmax Signal (7) = Signal (7)/Signalmax C Estimated SNR from measured pulse amplitude C SNRlimit=100. 0 C znep=sqrt (znep**2+ (1/SNRlimit*BWFactor*POWER) **2) C SNR = power*BWFactor/znep znepmeas=sqrt (znep**2+ (1/SNRlimit*BWFactor*Signalmax) **2) SNREST = (Signalmax * BWFactor)/znepmeas Return End End End CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXXXXXXXXXXXXX C*********************************************************** ************** C*********************************************************** ************** Subroutine SeekerNoise (Signal, AZ, EL) C*********************************************************** ************** C The seeker angle is then recostructed with the corrupted signals to approximate the C anticipated seeker angle noise. C Seeker angles are reconstructed from detector amplitudes C utilizing a coarse grid search followed C by a gradient search technique to minimize the residual. C*********************************************************** ************** IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Parameters C Parameter (NAZ = 81) Parameter (NEL = 81) Parameter (NDET = 7)

c-------------- Dimension Signal (NDET) Dimension OPT GAIN (NDET) Double Precision TotAng (NDET) Double Precision Theta Double Precision whatimeisit Logical Firstcall Dimension Anglens (3, NDET) COMMON/Gainlens/Nangles, Angles (200), OptTrans (200) COMMON/Caldata/AZtbl (NAZ), ELtbl (NEL), Tgain (NAZ, NEL, NDET) C C*********************************************************** ********** C*********************************************************** *********** C*********************************************************** *********** C Angle reconstruction from amplitudes for a given pulse begins here. C Do coarse grid seach to determine starting point for gradient search. C*********************************************************** *********** C*********************************************************** *********** C Initial center of search pattern and pattern density BestResidual = 9999. 0 Azcenter = 0. 0 Elcenter = 0. 0 DelAz = 1. 0 DelEl = 1. 0 9000 Continue Do 7000 IEL=-15, 15, 2 Do 8000 IAZ=-15, 15, 2 AZest = (dfloat (IAZ) *DelAz+Azcenter) If (Azest. gt. 15. OdO) Azest = 15. 0 If (Azest. lt.-15. OdO) Azest =-15. 0 ELest = (dfloat (IEL) *DelEl+ElCenter) If (Elest. gt. 15. OdO) ELest = 15. 0 If (ELest. lt.-15. OdO) ELest =-15. 0 C Determine residual for point in grid. Call DetResidual (AZest, ELest, Signal, Residual) C Save Az and El for search grid point having lowest residual C for further grid refinement If (Residual. LT. BestResidual) Then BestResidual=Residual Az=Azest El=Elest endif 8000 Continue 7000 Continue El=Elest C*********************************************************** *********** 8000 Continue 7000 Continue C Gradient Search Algorithm starts here. C*********************************************************** *********** C*********************************************************** *********** Residual = 9999. 0 DelAZ = 0. 5 DelEL = 0. 5 Do While (DelAZ. GT. 0. 01. OR. RESest. LT. Residual)

C Determine Base AZ, EL Residual Call DetResidual (AZ, EL, Signal, Residual) C Determine Local Residual Slope with AZ AZPDAZ=AZ + DelAZ Call DetResidual (AZPDAZ, EL, Signal, ResDAZ) AZReslope= (ResDAZ-Residual)/DelAZ C Determine Local Residual Slope with EL ELPDEL=EL+DelEL Call DetResidual (AZ, ELPDEL, Signal, ResDEL) ELReslope= (ResDEL-Residual)/DelEL C Determine direction weighting Factors (sin and cos) Denom=max (SQRT (AZReslope**2+ELReslope**2), 1. Od-6) AZfactor =-AZReslope/Denom ELFactor =-ELReslope/Denom If (Denom. LE. ld-6) then AZfactor=1. 0 ELfactor=1. 0 Endif C Determine New AZ, EL estimate and calculate residual AZest=AZ+AZfactor*DelAZ ELest=EL+ELfactor*DelEL Call DetResidual (AZest, ELest, Signal, RESest) C Accept new AZ, EL estimate if Residual improved C Otherwise assume step size to big (stepped over best solution) Itry=Itry+1 If (RESest. GT. Residual. or. Itry. gt. 1000) Then DelAZ=DelAZ/5. 0 DelEL=DelEL/5. 0 Itry=0 Else AZ=AZest EL=ELest Endif Endif C TheEnd of loop converging on angle estimate C Reduce search deltas and continue. C This is an intermediate iteration result. C C Write (10, 1500) AZtruth, ELtruth, AZ, EL, RESest C Write (10, 1500) AZtruth, ELtruth, AZ, EL, RESest 1500 Format (9F12. 6) Enddo C*********************************************************** *********** Enddo C TheEnd of Angle Reconstruction Algorithm C*********************************************************** *********** C*********************************************************** *********** C C TheEnd of loop for particular AZ, EL Truth Data Case Study C C End 5000 Return End

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXXXXXXXXXXXXX cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxXXXXXXXXXXXXX XXXXXXXXXxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx Subroutine DetResidual (AZest, ELest, Signal, Res) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C Parameters C Parameter (NAZ = 81) Parameter (NEL = 81) Parameter (NDET = 7) C Logical ISEARCH, IWARN Logical ISEARCH, IWARN Dimension Signal (NDET), Sigest (NDET) COMMON/Gainlens/Nangles, Angles (200), OptTrans (200) COMMON/Caldata/AZtbl (NAZ), ELtbl (NEL), Tgain (NAZ, NEL, NDET) C C Determine residual for all points in search grid C Determine residual for all points in search grid C Retrieve amplitudes for search grid Angles from Calibration DataBase ISEARCH =. TRUE. Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 1), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (1), ISEARCH, IWARN) ISEARCH =. FALSE. Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 2), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (2), ISEARCH, IWARN) Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 3), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (3), ISEARCH, IWARN) Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 4), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (4), ISEARCH, IWARN) Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 5), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (5), ISEARCH, IWARN) Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 6), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (6), ISEARCH, IWARN) Call Intrp2 (Aztbl, Eltbl, Tgain (1, 1, 7), NAZ, NEL, NAZ, NEL, &AZest, ELest, Sigest (7), ISEARCH, IWARN) C Normalize retrieved amplitudes from caibration database Sigestmax = DMAXl (Sigest (l), Sigest (2), Sigest (3), Sigest (4), &Sigest (5), Sigest (6), Sigest (7)) Sigest (1) = Sigest (1)/Sigestmax Sigest (2) = Sigest (2)/Sigestmax Sigest (3) = Sigest (3)/Sigestmax Sigest (4) = Sigest (4)/Sigestmax Sigest (5) = Sigest (5)/Sigestmax Sigest (6) = Sigest (6)/Sigestmax Sigest (7) = Sigest (7)/Sigestmax C Calculate Residual Component for Edge of FOV FOVest=SQRT (AZest**2+ELest**2) FOVres=0. 0 IF (FOVest. GT. 14. 0) FOVres = (FOVest-14. 0)/14. 0

C Calculate residual for each point in search grid C Res = Signal (1) * (Signal (1)-Sigest (1)) **2 + C & Signal (2) * (Signal (2)-Sigest (2)) **2 + C & Signal (3) * (Signal (3)-Sigest (3)) **2 + C & Signal (4) * (Signal (4)-Sigest (4)) **2 + C & Signal (5) * (Signal (5)-Sigest (5)) **2 + C & Signal (6) * (Signal (6)-Sigest (6)) **2 + C & Signal (7) * (Signal (7)-Sigest (7)) **2 + C & FOVres**2 C Calculate residual for each point in search grid Res = (Signal (1)-Sigest (1)) **2 + &(Signal (2)-Sigest (2)) **2 + &(Signal (3)-Sigest (3)) **2 + &(Signal (4)-Sigest (4)) **2 + &(Signal (5)-Sigest (5)) **2 + &(Signal (6)-Sigest (6)) **2 + &(Signal (7)-Sigest (7)) **2 + &FOVres**2 Return END END xxxxxxxxxxxxxxxxxxxxxxxxxxxx c SUBROUTINE DETECT (BACKPOW, GAIN, TEMP, BW, ZNEP) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C C THIS SUBROUTINE COMPUTES THE NOISE EQUIVALENT POWER, C INCLUDING BACKGROUND/BACKSCATTER EFFECTS C C C INPUTS- C C BACKPOW = BACKSCATTER POWER, WATTS C GAIN = APD GAIN C TEMP = TEMPERATURE OF AMP AND DETECTOR, KELVIN C BW = Effective Noise BANDWIDTH OF FILTER ON AMP, HZ C C OUTPUTS- C C ZNEP = NOISE EQUIVALENT POWER, WATTS C C VARIABLE DEFINITIONS- C C PI = PI C ZIDS = DARK SURFACE CURRENT OF APD AT 25C, AMPS C ZIDB = DARK, UNITY GAIN, BULK CURRENT OF APD AT 25C, AMPS C RES = FEEDBACK RESISTOR, OHMS C OP1 = AMP CURRENT NOISE, AMPS/RHZ C ENOP1 = AMP VOLTAGE NOISE, VOLTS/RHZ C CT = TOTAL CAPACITANCE AT INPUT TO AMP (INCLUDES APD CAP), FARADS C ZKEFF = EFFECTIVE IONIZATION RATIO IN APD C WAVEL = LIGHT WAVELENGTH, MICRONS C QE = QUANTUM EFFICIENCY C QELECT = CHARGE ON AN ELECTRON, COULOMBS C RNOISE = FEEDBACK RESISTOR NOISE, AMPS C ZINOP1 = AMPLIFIER CURRENT NOISE, AMPS C ZINENOP1 = AMPLIFIER VOLTAGE NOISE CONVERTED TO CURRENT, AMPS C FSNR = EXCESS NOISE FACTOR, DIMENSIONLESS

C RESP = DETECTOR RESPSONSIVITY, AMPS/WATT C ZIDSTC = ZIDS COMPENSATED FOR TEMPERATURE, AMPS C ZIDBTC = ZIDB COMPENSATED FOR TEMPERATURE, AMPS C TOTALN = TOTAL NOISE, AMPS C THRESH = DETECTION THRESHOLD C C C C **OPEN FILES**************** C OPEN (UNIT=12, FILE='OUTPUT. DAT') C **SENSOR CHARACTERISTICS**************** ZIDS = 150. E-9 ZIDB = 2. OE-12 RES = 20. E+3 OP1 = 1833. 0E-15 ENOP1 = 1. 3e-9 CT = 15. OE-12 ZKEFF = 0. 02 WAVEL = 1. 064 QE = 0. 36 C C COMPUTE THE NOISE AND THE THRESHOLD C COMPUTE THE NOISE AND THE THRESHOLD C PI=3. 141592654 QELECT = 1. 6E-19 RNOISE = SQRT (4. *1. 38E-23*TEMP*BW/RES) ZINOP1 = OP1*SQRT (BW) ZINENOP1 = (ENOP1*2. *PI*CT) *SQRT (0. 3333*BW*BW*BW) FSNR = ZKEFF*GAIN + (2.-1./GAIN) * (1-ZKEFF) RESP = GAIN*QE*WAVEL/1. 24 ZIDSTC = ZIDS/ (2** ( (298-TEMP)/10)) ZIDBTC = ZIDB*GAIN/ (2** ( (298-TEMP)/10)) ZIBG = BACKPOW*RESP ZIBG NOISE = SQRT (BW*2. *QELECT*ZIBG*GAIN*FSNR) c ZINDET = SQRT (BW*2. *QELECT* ( (ZIBG + ZIDBTC) * c C GAIN*FSNR + ZIDSTC)) ZINDETBG = SQRT (BW*2. *QELECT*ZIBG*GAIN*FSNR) ZINDETDB = SQRT (BW*2. *QELECT*ZIDBTC*GAIN*FSNR) ZINDETDS = SQRT (BW*2. *QELECT*ZIDSTC) ZINDET = SQRT (ZINDETBG**2 + ZINDETDB**2 + ZINDETDS**2) TOTALN = SQRT (RNOISE**2 + ZINOP1**2 + ZINENOP1**2 + ZINDET**2) ZNEP = TOTALN/RESP RETURN END C*********************************************************** ************** C*********************************************************** ************** FUNCTION ZRAND () IMPLICIT DOUBLE PRECISION (A-H, O-Z) ZRAND = RAND () IF (ZRAND. EQ. 1) THEN ZRAND = 0. 99999999999DO ELSE IF (ZRAND. EQ. 0) THEN ZRAND = 1. E-20 ENDIF

END FUNCTION C*********************************************************** ************** C*********************************************************** ********** Double Precision Function OPgain3 (Theta) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C THIS SUBROUTINE RETURNS OPTICAL GAIN BASED ON INTERPOLATING IN TABULAR DATA C C DECLARATIONS LOGICAL ISEARCH, IWARN COMMON/Gainlens/Nangles, Angles (200), OptTrans (200) C C ISEARCH =. TRUE. ISEARCH =. TRUE. IWARN =. False. CALL INTRP1 (Angles, OPTTRANS, 200, NANGLES, THETA, Opgain3 &, ISEARCH, IWARN) IF (IWARN) THEN PRINT*,'SUBROUTINE OPTICAL DATA OUT OF BOUNDS !', THETA WRITE (7, *)'SUBROUTINE OPTICAL DATA OUT OF BOUNDS !', THETA ENDIF END FUNCTION C*********************************************************** ********* END FUNCTION SUBROUTINE SKIPL (IUNIT, NLINES) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C C THIS ROUTINE IS USED TO SKIP LINES IN INPUT FILES C THIS ROUTINE IS USED TO SKIP LINES IN INPUT FILES C SKIPL DO 10 I=1, NLINES READ (IUNIT, *) 10 CONTINUE RETURN END C*********************************************************** *********** END SUBROUTINE INTRP1 (XA, FA, MAXX, NX, X, F, ISEARCH, IWARN) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** ******************* C** C** Single interpolation subroutine C** C** ARGUMENT DEFINITIONS : C** C** XA-Array of first independent variable (input) C** FA-Array of dependent variable (input) C** X-Value of first independent variable (input) C** F-Value of dependent variable (output) C** MAXX-Dimension of XA and first dimension of FA (input) C** NX-Actual number of values in XA array (input) C** ISEARCH-Logical flag indicating whether it is necessary to search C** the XA array for the appropriate interpolation C** range. C** IWARN-Logical flag indicating whether to warn of X outside C** of range. C**

c*********************************************************** ******************* LOGICAL ISEARCH, IWARN DIMENSION XA (MAXX), FA (MAXX) SAVE KX1, KX2, FACTX DATA KX1, KX2/2*1/ IF (ISEARCH) THEN CALL SEARCH (XA, MAXX, NX, X, KX1, KX2, FACTX, IWARN) ENDIF F = FA (KX1) + (FA (KX2)-FA (KX1)) *FACTX RETURN END RETURN END SUBROUTINE INTRP2 (XA, YA, FA, MAXX, MAXY, NX, NY, 1 X, Y, F, ISEARCH, IWARN) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** ******************* C** C** Double interpolation subroutine C** C** ARGUMENT DEFINITIONS : C** C** XA-Array of first independent variable (input) C** YA-Array of second independent variable (input) C** FA-Array of dependent variable (input) C** X-Value of first independent variable (input) C** Y-Value of second independent variable (input) C** F-Value of dependent variable (output) C** MAXX-Dimension of XA and first dimension of FA (input) C** MAXY-Dimension of YA and second dimension of FA (input) C** NX-Actual number of values in XA array (input) C** NY-Actual number of values in YA array (input) C** ISEARCH-Logical flag indicating whether it is necessary to search C** the XA and YA arrays for the appropriate interpolation C** range. C** IWARN-Logical flag indicating whether to warn of X or Y outside C** of range. C** C*********************************************************** ******************* LOGICAL ISEARCH, IWARN DIMENSION XA (MAXX), YA (MAXY), FA (MAXX, MAXY) SAVE KX1, KX2, KY1, KY2, FACTX, FACTY DATA KX1, KX2, KY1, KY2/4*1/ IF (ISEARCH) THEN CALL SEARCH (XA, MAXX, NX, X, KX1, KX2, FACTX, IWARN) CALL SEARCH (YA, MAXY, NY, Y, KY1, KY2, FACTY, IWARN) ENDIF TEMPX1 = FA (KX1, KY1) + (FA (KX2, KY1)-FA (KX1, KY1)) *FACTX TEMPX2 = FA (KX1, KY2) + (FA (KX2, KY2)-FA (KX1, KY2)) *FACTX F = TEMPX1 + (TEMPX2-TEMPX1) *FACTY RETURN END RETURN END SUBROUTINE SEARCH (XA, MAXX, NX, X, Kl, K2, FACTOR, IWARN)

IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** ******************* C** C** This routine is called by INTRP1, INTRP2, or INTRP3 to search through C** a vector table, XA, to find the subscripts Kl and K2 such that C** C** XA (K1) <= X <= XA (K2). C** C** In general, C** C** K2 = Kl + 1 C** C** although the subscripts may be equal if X lies on or outside either end C** of the vector. C** Note that the XA vector must be arranged in increasing order. C** C*********************************************************** ******************* C** C** ARGUMENT DEFINITIONS : C** C** IWARN-Logical flag indicating whether to warn of X or Y outside C** of range. C** FACTOR-Ratio (X-XA (K1))/ (XA (K2)-XA (K1)) to be used for C** linear interpolation by the calling routine. If X < XA (1) C** or X > XA (NX), then FACTOR = 0 is returned. C** K1-Largest index of XA array having value less than or equal C** to X (output) C** K2-Smallest index of XA array having value greater than or equal C** to X (output) C** MAXX-Dimension of XA (input) C** NX-Actual number of values in XA array (input) C** XA-Vector array of independent variable (input) C** X-Value of independent variable (input) C** C*********************************************************** ******************* LOGICAL IWARN DIMENSION XA (MAXX) C WRITE (*, *)' ENTER K2 =', K2 C C--Check for X outside of range of XA C--Check for X outside of range of XA IF (X. LT. XA (1)) THEN K1=1 K2=1 FACTOR=0. IWARN =. TRUE. IF (IWARN) THEN WRITE (*, *)'X. LT. XA (1) IN SEARCH ; X=', X WRITE (7, *)'X. LT. XA (1) IN SEARCH ; X=', X Read (*, *) ENDIF GOTO 900 ELSE IF (X. GT. XA (NX)) THEN K1=NX K2=NX FACTOR=0. IWARN =. TRUE. IF (IWARN) THEN WRITE (*, *)'X. GT. XA (NX) IN SEARCH ; X=', X WRITE (7, *)'X. GT. XA (NX) IN SEARCH ; X=', X Read (*, *) ENDIF GOTO 900

ENDIF ENDIF C--Find indices of XA for which XA (K1) <= X <= XA (K2) ICOUNT = 0 10 CONTINUE IF (X. LT. XA (K1)) THEN Kl = Kl-1 ICOUNT = ICOUNT + 1 GOTO 10 ENDIF K2 = MIN0 ((K1 + 1), NX) 20 CONTINUE C WRITE (*, *)' K2 =', K2 IF (X. GT. XA (K2)) THEN K2 = K2 + 1 ICOUNT = ICOUNT + 1 GOTO 20 ENDIF K1 = MAXO ( (K2-1), 1) FACTOR = (X-XA (K1))/ (XA (K2)-XA (K1)) 900 CONTINUE C WRITE (*, *)' RETURN K2 =', K2 RETURN END RETURN END SUBROUTINE UNIT (A, AU) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** ******** C C FUNCTION : COMPUTE A UNIT VECTOR C C*********************************************************** ******** DIMENSION A (3) DIMENSION AU (3) C* AMAG=VECMAG (A) DO 10 I=1, 3 AU (I) =A (I)/AMAG 10 CONTINUE RETURN END RETURN END SUBROUTINE VECSUB (A, B, R) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** *** C C FUNCTION : SUBTRACT TWO VECTORS C C*********************************************************** *** DIMENSION A (3) DIMENSION B (3) DIMENSION R (3)

DO 10 I=1, 3 R (I) =A (I)-B (I) 10 CONTINUE RETURN END END SUBROUTINE SCAVEC (S, V, SV, N) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C*********************************************************** *** C C SCAVEC MULTIPLIES A SCALAR TIMES A VECTOR C C*********************************************************** *** DIMENSION V (N), SV (N) DO 10 I=1, N SV (I) =S*V (I) 10 CONTINUE RETURN END C*********************************************************** *** C*********************************************************** *** DOUBLE PRECISION FUNCTION VECMAG (A) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C********************************************************** C C FUNCTION : COMPUTE A VECTOR MAGNITUDE C C********************************************************** DIMENSION A (3) R=O. ODO DO 10 I=1, 3 R=R+ A (I) *A (I) 10 CONTINUE VECMAG = DSQRT (R) RETURN END C*********************************************************** ************* C*********************************************************** ************* C*********************************************************** ************* C*********************************************************** ************* RETURN END

**Previous Patent:**METHOD AND APPARATUS FOR MAKING HIGH-PRECISION MEASUREMENTS

**Next Patent: FLUID FLOW METER FOR GRAVITY FED INTRAVENOUS FLUID DELIVERY SYSTEMS**