WO2003048689A2 - Methods and apparatuses for reconstructing angle information - Google Patents

Abstract

Methods and apparatuses for determining angular information from a target. A plurality of sensors on an image plane are used to obtain measured amplitude data. An azimuth and elevation may be estimated using this measured amplitude data. A residual error between the measured amplitude data and calibrated amplitude data may then be determined. The angular information corresponding to the azimuth and elevation at which the residual error is minimized may then be outputted as representing or corresponding to the angular information of the target.

Description

DESCRIPTION

METHODS AND APPARATUSES FOR RECONSTRUCTING ANGLE INFORMATION

BACKGROUND OF THE INVENTION

This application claims priority to U.S. Provisional Patent Application Serial No. 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, A_meas means the measured amplitude from detector A, while A_eSt 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 AngleRecorfdof ************************************************************************ *

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. Q* *********************************************************************** *

USE MSIMSL USE MSFLIB C USE PORTLIB

IMPLICIT DOUBLE PRECISION (A-H,θ-Z)

Parameter (NDET = 7) Dimension Signal (NDET) Dimension AZ ( 1000) , EL(IOOO) Double Precision SNR, SNREST C

Q* ************* * **********************************************************

C All of these files are not necessarily used every time. These are C All of these files are output files.

0PEN(UNIT=7, FILE=' Amplitude_Noise.xls' ) Write (7,*) 'TargetRange AZtruth ELtruth SNR AZnoise ELnoise' OPEN(UNIT=8, FILE=' CalMatrix_20X20.txt ' )

0PEN(UNIT=9, FILE=' Errormap.xls' )

OPEN(UNIT=10, FILE=' Anglemeas.xls' )

OPEN(UNIT=ll_/ FILE=' ngleres.xls' )

Q* _****_***_***_*******************************************************_*****_{* *} c Nuber of Monte Carlo noise cases for shotgun plot and noise statistics

Nmontenoise = 100 (2******** *************************************** **************************

Do 6000 Iseekercase=l,3, 1

If (Iseekercase.EQ.1) Then Aztruth = 0.0 Eltruth = 0.0

0PEN(UNIT=12, FILE= ' AnglestatszO-0.xls ' ) Write(12,*) 'TargetRange, AZtruth, ELtruth, Nmontenoise,AZmean,

&ELmean, AZrmsOO/0 , ELrmsOO/0 '

Elself (Iseekercase .EQ.2) Then Aztruth = 5.00 Eltruth = 0.00

0PEN(UNIT=12, FILE= 'Anglestatsz5-0.xls ' )

Write (12, *) 'TargetRange, AZtruth, ELtruth, Nmontenoise, AZmean, &ELmean,AZrms@5/0, ELrms@5/0 ' Elself (Iseekercase.EQ.3) Then Aztruth = 10.0 Eltruth = 0.0 0PEN(UNIT=12, FILE= ' AnglestatszlO-0.xls ' )

Write ( 12 , * ) ' TargetRange , AZtruth, ELtruth, Nmontenoise , AZmean, &ELmean,AZrms@10/0,ELrms@10/0 '

Elself (Iseekercase . EQ .4 ) Then Aztruth = 4.00 Eltruth = 7.00 Elself (Iseekercase . EQ .5 ) Then Aztruth = 8.50 Eltruth = -10.75 Elself (Iseekercase . EQ .6 ) Then Aztruth = 0.50 Eltruth = -13.00 Elself (Iseekercase .EQ.7) Then Aztruth -12.75 Eltruth 12.75 Elself (Iseekercase .EQ.8) Then Aztruth = 16.9 Eltruth = 10.25 Elself (Iseekercase . EQ .9 ) Then Aztruth = -2.25 Eltruth = -5.00 Elself (Iseekercase .EQ.IO) Then Aztruth = -10.75 Eltruth = -3.00 Elself (Iseekercase .EQ.ll) Then Aztruth = 13.25 Eltruth = 11.75 Endif c AZtruth=float (Iseekercase) /10.0 c ELtruth=0.0

Q********************************************************************** _c* ****************** This is the MonteCarlo Noise Loop **************** Do 1000 lrange=100,6000, 100

Do 2000 Inoise=l, 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(lnoise) =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 ****************************

Q* ****************************************************************** * close (12) 6000 Continue Stop END _C******************_*******_****_***_******_********************************_***

C END of trash code used to excercize the subroutine for seeker noise.

Q*********************************************_**_*****_********************* Q* ************** * ************************************************** *******

Subroutine AngleStats (N,AZ, EL, AZmean, ELmean, AZrms, ELrms)

IMPLICIT DOUBLE PRECISION (A-H.O-Z) Dimension AZ(IOOO), 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-l) ) ELrms=SQRT(ELres/Dfloat (N-l) ) return end

_C************************************************************************* **************** *********************************************** ********** **************** *************************************** ******************

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 althought 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

Q* *********************************************************************** *

IMPLICIT DOUBLE PRECISION (A-H,θ-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, DET)

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.0d0/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 0., .001166,, 00., .001111,, 00.. .000055,, 00.. .000000,, 00., .000000, 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 This is the only required input file. c

If (Firstcall) then Firstcall= . alse . C Point of each of 7 lenses

Q ******Lerιs Pattern************************ C 4 3

Q******** ***************************** ******

Anglens(l,l) = 0.0 Anglens(2,l) = 0.0

Anglens(l,2) = 14.00 Anglens(2,2) = 0.0 Anglens(l,3) = 7.0

Anglens(2,3) = 12.12

Anglens(l,4) = -7.00

Anglens(2,4) = 12.12

Anglens(l,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

Q******** ******************* ********************************************** Q*** ********************************************************** ************

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.5d0-20.5d0) EL = (dfloa (IEL)*0.5d0-20.5d0)

Aztbl (IAZ) = AZ Eltbl(IEL) = EL C Determine total angle in radians to each lens Do 3000 Idet = 1 , 7

TotAng ( Idet ) =sqrt ( (Anglens ( 1 , 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 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. θpt_eff= 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.0e6 BWnoise = BW*PI Tau=l/ (2*pi*BW) ; BWfactor= (l-exp(-PlWid/Tau) )

C Power of the laser (Watts) POWER_LASER = 5.0e6 C Atmospheric extintion coefficient (km) C EXT_COEF = -0.09135 Visibility = 23.5d0 EXT_COEF=- ( 10.0d0**(-0.136d0+(1.16d0*logl0( 3.912d0/Visibility) ) ) ) C Target Reflectivity REFLTAR = 0.1

C Reflectivity of background (i.e. terrain) REFLBACK = 0.1 C Number of detectors per aperture Num_det = 1

C Number of apertures in the reciever NUM_APERT_RECV = 3 NUM_APERT_BACK = 4

C INCIDENCT SPECTRAL FLUX WATTS/ (M^*2 * MICRON) 1.06 microns SPEC_FLUX = 500.

C OPTICAL FILTER BANDWIDTH IN MICRONS FIL BW = .010

C Range in meters C Range = 6000

Endif C End of calculations done for first subroutine call only C c

Q********** ******************************************************** *******

C Note reference to optical gain at 0 radians offset.

C Note assumes designator fixed at 6km while target range is variable. EXTINCT = EXP( EXT_COEF* (TargetRange+DesignatorRange)/1000)

POWER = POWER_LASER*EXTINCT*REFLTAR*APERTURE* & FLOAT (NUM_APERT_RECV) *0PT_EFF*0PGAIN3 (0.0)/ & (3.14159*TargetRange**2*FLOAT(NUM_DET) ) BACKPOW = (TLFOV)**2*SPEC_FLUX*FIL_BW*APERTURE*NUM_APERT_BACK*

& REFLBACK*OPT EFF*OPGAIN3 (0.0) / (NUM DET*8.)

TotalPower = POWER + BACKPOW Call DETECT (TotalPower, GAIN, TEMP, BWnoise, ZNEP)

SNRlimit=1000.0 zneplim=sqrt (znep**2+ (l/SNRlimit*BWFactor*POWER) **2) SNR = power*BWFactor/zneplim Q* _{********* ***************************************************************} 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(l)=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 θpt_gain(l) = OPGain3 (TotAng (1) ) /OPGain3 (0.0) θpt_gain(2) = OPGain3 (TotAng (2) ) /OPGain3 (0.0) θpt_gain(3) = 0PGain3 (TotAng (3) ) /OPGain3 (0.0) θρt_gain (4 ) = OPGain3 (TotAng(4) )/OPGain3 (0.0) Opt_gain(5) = OPGain3 (TotAng (5) )/OPGain3 (0.0) θpt_gain(6) = OPGain3 (TotAng (6) )/OPGain3 (0.0) θpt_gain(7) = 0PGain3 (TotAng (7) )/0PGain3 (0.0)

C Add noise to signal

Signal (1) = power*OPT_GAIN(l) + znep*dnorin zrandO )

Signal (2) = power*0PT_GAIN(2) + znep*dnorin zrand( ) )

Signal (3) = power*0PT_GAIN(3) + znep*dnorin zrand( ) )

Signal (4) = power*OPT_GAIN(4) + zneρ*dnorin( zrand ( ) )

Signal (5) = power*OPT_GAIN(5) + znep*dnorin zrand ( ) )

Signal (6) = power*0PT_GAIN(6) + znep*dnorin ( zrand ( ) )

Signal (7) = power*OPT_GAIN(7) + znep*dnorin zrand () )

Normalize all measured amplitudes to maximum

Signalmax = MAX (Signal (1) , 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+(l/SNRlimit*BWFactor*POWER) ^►2) C SNR = power*BWFactor/znep

znepmeas=sqrt (znep**2+ (l/SNRlimit*BWFactor*Signalmax) **2) SNREST = (Signalmax * BWFactor) / znepmeas

Return End

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

Q* * ***********************************************************************

Q* ***************************************************************** *******

Subroutine SeekerNoise (Signal , AZ, EL)

Q** ************************************************************* ********** 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. _Q*********** *************_*************_***************************** *******

IMPLICIT DOUBLE PRECISION (A-H,θ-Z)

C Parameters C

Parameter (NAZ = 81) Parameter (NEL = 81) Parameter (NDET = 7) Dimension Signal (NDET) Dimension OPT_GAI (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)

Q* ******************************************************************* * Q* ******************************************************************** * Q* ******************************************************************** *

C Angle reconstruction from amplitudes for a given pulse begins here. C Do coarse grid seach to determine starting point for gradient search. Q**** ******************************************************** ********** Q* ******************************************************************** * 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.0d0) 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 .It. -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

Q* ******************************************************************** * Q******* ******** ****************************************** *************

C Gradient Search Algorithm starts here.

(2*** **************************************************************** *** c_{* ************************************************** *******************}

Residual = 9999.0 DelAZ = 0.5 DelEL = 0.5 Do While (DelAZ. GT.0.01.OR. RESest .LT. esidual) 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=l 0 ELfactor=l 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+l

If (RESest GT Residual or Itry gt 1000) Then DelAZ=DelAZ/5 0

DelEL=DelEL/5 0 ltry=0

Else

AZ=AZest EL=ELest

Endif _c

C TheEnd of loop converging on angle estimate C Reduce search deltas and continue

C This is an intermediate iteration result

C Write (10, 1500) AZtruth, ELtruth, AZ, EL, RESest 1500 Format (9F12 6)

Enddo ( * ******************************************************************** * Q******** ********************************************* *****************

C TheEnd of Angle Reconstruction Algorithm

(2********************************************************************** Q* ******* ******************************************** ****************** C C TheEnd of loop for particular AZ, EL Truth Data Case Study _c c c

5000 Return

End cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

(2** * * * * * * **************************************************** ************ Subroutine DetResidual (AZest, ELest , Signal , Res) IMPLICIT DOUBLE PRECISION (A-H,θ-Z) _c

C Parameters

C Parameter (NAZ = 81)

Parameter (NEL = 81) Parameter (NDET = 7)

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 Retrieve amplitudes for search grid Angles from Calibration DataBase

ISEARCH = .TRUE.

Call Intrp2 (Aztbl, Eltbl, Tgain(l, 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, Tgaind, 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, AZ , EL, & AZest, ELest, Sigest (5) , ISEARCH, IWARN)

Call Intrp2 (Aztbl, Eltbl, Tgaind, 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 = DMAX1 (Sigest (1) , 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) /1 .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) -Siges (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 + & F0Vres**2

Return END CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXXXXXXXXXXXXXX c

SUBROUTINE DETECT (BACKPOW, GAIN, TEMP, BW, ZNEP) IMPLICIT DOUBLE PRECISION (A-H,θ-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 OPl = AMP CURRENT NOISE, AMPS/RHZ

C ENOPl = 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 ZINOPl = AMPLIFIER CURRENT NOISE, AMPS

C ZINENOPl = 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

**OPEN FILES****************

C 0PEN(UNIT=12, FILE= ' OUTPUT.DAT ' )

C **SENSOR CHARACTERISTICS****************

ZIDS = 150. E-9 ZIDB = 2.0E-12 RES = 20.E+3 OPl = 1833.0E-15 ENOP1 = 1.3e-9

CT = 15.0E-12 ZKEFF = 0.02 WAVEL = 1.064 QE = 0.36 _c

C COMPUTE THE NOISE AND THE THRESHOLD

PI=3.141592654 QELECT = 1.6E-19

RNOISE = SQRT(4.*1.38E-23*TEMP*BW/RES)

ZINOP1 = OPl*SQRT(BW)

ZINENOPl = (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 + ZINOPl**2 + ZINENOPl**2 + ZINDET**2)

ZNEP = TOTALN/RESP

RETURN

END

(2**************************** * ******************************************** ******* ********************************************** ********************

FUNCTION ZRAND ()

IMPLICIT DOUBLE PRECISION (A-H,θ-Z) ZRAND = RANDO

IF (ZRAND .EQ. 1) THEN

ZRAND = 0.99999999999D0 ELSE IF (ZRAND .EQ. 0) THEN

ZRAND = l.E-20 ENDIF END FUNCTION n***** ***************************************************** *************** ₍2******** ***************************************************** ********

Double Precision Function 0Pgain3 (Theta) IMPLICIT DOUBLE PRECISION (A-H,θ-Z)

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. 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** ****************************************************************** C*********************************************************************

SUBROUTINE SKIPL(IUNIT, NLINES) IMPLICIT DOUBLE PRECISION (A-H.O-Z)

THIS ROUTINE IS USED TO SKIP LINES IN INPUT FILES SKIPL

DO 10 1=1, NLINES READ UNIT,*) 10 CONTINUE RETURN

END

C********************************************************************** C**********************************************************************

SUBROUTINE INTRP1 (XA, FA, MAXX, NX, X, F, ISEARCH, IWARN)

IMPLICIT DOUBLE PRECISION (A-H,θ-Z)

C****************************************************************************** C** C** Single interpolation subroutine Q**

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** * *********************************************************************** ******

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 , KX , FACTX, IWARN) ENDIF

F = FA(KXl) + (FA(KX2) - FA(KXl) ) *FACTX

RETURN END

(2** ****** ************************************************************

*******************************************************************************************

SUBROUTINE INTRP2 (XA, YA, FA, MAXX, MAXY, NX, NY, 1 X, Y, F, ISEARCH, IWARN)

IMPLICIT DOUBLE PRECISION (A-H,θ-Z)

C** **************************************************************************** C**

Q** Double interpolation subroutine C**

Q** ARGUMENT DEFINITIONS:

C** Q** 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) ** 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)

(2** ISEARCH - Logical flag indicating whether it is necessary to search C** the XA and YA arrays for the appropriate interpolation

Q** 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, KYI, KY2 , FACTX, FACTY DATA KX1, KX2, KYI, KY2 /4*l/

IF (ISEARCH) THEN

CALL SEARCH (XA, MAXX, NX, X, KX1 , KX2 , FACTX, IWARN) CALL SEARCH (YA, MAXY, NY, Y, KYI, KY2 , FACTY, IWARN)

ENDIF

TEMPX1 = FA(KX1,KY1) + (FA(KX2,KY1) FA (KX1 , KYI ) ) *FACTX TEMPX2 = FA(KX1,KY2) + (FA(KX2,KY2) FA (KX1 , KY2 ) ) *FACTX

F = TEMPX1 + (TEMPX2-TEMPX1)*FACTY

RETURN END

Q* **************************************************************************** * Q*** ***************************************************************** **********

SUBROUTINE SEARCH(XA, MAXX, NX, X, Kl , K2 , FACTOR, IWARN) IMPLICIT DOUBLE PRECISION (A-H,0-Z) Q*** ******************************************************** *******************

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**

Q* **************************************************************************** * 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(KD) 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** Kl - 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** _Q***** ***********************************_*************************_********_*****

LOGICAL IWARN DIMENSION XA(MAXX) C WRITE (*,*)' ENTER K2 = ' , K2

C

C-- Check for X outside of range of XA

IF (X.LT.XA(l)) THEN Kl=l

K2=l

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 _c

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( (Kl + 1) , NX) 20 CONTINUE

C WRITE (*,*)' K2 = ' , K2

IF ( X .GT. XA(K2) ) THEN K2 = K2 + 1 ICOUNT = ICOUNT + 1 GOTO 20

ENDIF

Kl = MAX0 ( (K2 - 1) , 1) FACTOR = (X - XA(K1) )/ (XA(K2) - XA(K1))

900 CONTINUE

C WRITE (*,*)' RETURN K2 = ' , K2 RETURN

END

Q* ************************************************************ *

(2* ************************************************************ *

SUBROUTINE UNIT(A,AU)

IMPLICIT DOUBLE PRECISION (A-H,0-Z) (2** **************************************************** *************

C C FUNCTION: COMPUTE A UNIT VECTOR C

Q* ******************************************** **********************

DIMENSION A ( 3) DIMENSION AU ( 3) c*

AMAG=VECMAG(A) DO 10 1=1,3

AU ( I ) =A ( I ) /AMAG 10 CONTINUE RETURN

END

Q* ************************************************ **************

Q* ************************************************************ *

SUBROUTINE VECSUB (A, B , R) IMPLICIT DOUBLE PRECISION (A-H,θ-Z) (2****** ******************************************** ************

C C FUNCTION: SUBTRACT TWO VECTORS

C

(2*** ************************************************* **********

DIMENSION A ( 3)

DIMENSION B ( 3) DIMENSION R ( 3) DO 10 1=1,3

R(I)=A(I)-B(I) 10 CONTINUE RETURN

END

Q* ************************************************************************* *

SUBROUTINE SCAVEC (S, V, SV,N)

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

Q* ************************************************************ * c

C SCAVEC MULTIPLIES A SCALAR TIMES A VECTOR c

(2****** ******************************************************** DIMENSION V(N),SV(N)

DO 10 1=1, N SV(I)=S*V(I)

10 CONTINUE RETURN END (2***** ********************************************* ************ _Q*********** **********************_{** *}*************************_*

DOUBLE PRECISION FUNCTION VECMAG(A) IMPLICIT DOUBLE PRECISION (A-H,0-Z) c_{*** ************************** *****************************} c

C FUNCTION: COMPUTE A VECTOR MAGNITUDE

C

Q* ******************************************************** * DIMENSION A(3)

R=O.0D0

DO 10 1=1,3

R=R+ A(I)*A(I) 10 CONTINUE

VECMAG = DSQRT(R) RETURN END Q*** ** ******************************************************************* c_{*** ***************************************************** **************}*_* * *********************************************************************** c_{*** ********************************************************* ************}

(2*** ********************************************************* ************ (2* ********************************************************************** *

Claims

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.

1 1. 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.