Micro-miniature beacon transmit-only geo-location emergency system for personal security

[0002] This application is a continuation-in-part of prior application Ser. No. 09/052,066 filed Apr. 6, 1998, now U.S. Pat. No. 6,028,551, which was a continuation-in-part of prior application Ser. No. 08/904,377 filed Aug. 1, 1997 and subsequently abandoned, which was a continuation of the original application Ser. No. 08/355,901 filed Dec. 13, 1994 and subsequently abandoned.

[0003] The radio frequency spectrum has been used in the past for emergency rescue operations, most notably in shipping and aircraft rescues. More recently, geo-location devices have been used to track endangered species. As mobile telephony has developed, the ability to signal for help in emergency situations has also developed. However, in all of the above cases, the geo-location systems either require two-way communications, require relatively long time intervals to acquire the data needed for geo-location, or lack continuous coverage capability. In addition, the emergency transmit/receive units required are relatively large and not convenient to “wear” or conceal, and the battery systems and modes of operation require frequent battery recharging. Also inhibiting continuous use for emergency situations are the problems associated with incapacitation, such as the dialing sequence of operations necessary on current mobile phone systems.

[0004] There are a number of geo-location techniques that have been described in the past that appear in this specification. Prior art in emergency rescue at sea and for tracking of endangered species includes a number of existing satellite and ground systems dating back to World War II. These include position determination systems used in navigation, such as LORAN (a form of triangulation), Direction Finding (DF), Transit (Doppler shifts) and GPS (another form of triangulation), as well as emergency beacon location systems such as SARSAT (Doppler shifts) and Argos (Doppler shifts). A brief description of these techniques and the weaknesses relative to the proposed invention is provided in this section, as prior art background.

[0005] One of the earliest techniques, direction-finding (DF), was developed during World War II and evolved into systems in use today, such as LORAN and associated techniques using time/phase/angular difference of arrival, to determine angular direction and errors. DF techniques are suitable for use with ground-based cellular receiver networks, such as the current mobile telephone grids, but will not provide sufficient accuracy when used in conjunction with global capability satellite networks that will soon be available. This is fundamentally due to the great distances the satellites are from the transmitters (typically >700 km), and the impact created by this multiplicative factor that transforms the angular uncertainty of the measurement into a position error. However, DF may be used in conjunction with this invention for certain applications where rapid proximity to the transmitter is possible.

[0006] Another approach, developed over the last few decades, utilizes Doppler shifts in the frequency of received signals to determine the location of the transmitter. It is well known in physics text books that the shift in frequency depends on the factor (v/c) cos(θ), where the cosine term is the angle between the velocity vector of the satellite and the satellite-beacon line-of-site vector, and c is the speed of light. NASA has employed this technique for it's SARSAT program (search and rescue satellite). One of the drawbacks of that system is that it takes several satellite passes to determine location (it determines the zero-Doppler shift direction which provides a direction perpendicular to the satellite orbit), which means that accurate locations cannot be determined for the order of hours, the time it takes for multiple satellite passes. Newer approaches (such as the Orbcomm system, which will be fully operational in two years) can determine position to the order of a few hundred meters, but takes a complete pass (of the order of 10-15 minutes) of a single satellite to acquire sufficient data. In contrast, the proposed system will only need a simultaneous single data reception by a few visible satellites to determine location. This will minimize the energy required to provide a sufficient operational lifetime for the invention.

[0007] Time difference of arrival (TDOA) and time of arrival (TOA) techniques can be considered a subset of triangulation and direction-finding, but require a minimum of three simultaneous measurements to determine position unambiguously (by the intersection of the transit-time “circles” (TOA) or the intersection of hyperbolas from pairs of TDOA relative measurements). However, TDOA/TOA in combination with Doppler can reduce this constraint, provided that clock and/or relative times of arrival uncertainties can be reduced to provide the needed accuracy. This synergism can be used in the proposed invention to reduce the number of satellites that must receive the beacon signal.

[0008] The Global Positioning System (GPS), which was developed and launched by DoD and contractor Rockwell, depends on receiver units being able to lock-on to signals from three or more GPS satellites to calculate a position. As opposed to the system proposed herein, GPS transmitters are in space, and the receiver is carried by the user; thus it is a “receive-only” system. The position errors, depending on the code the receiver is authorized to receive, are extremely small (of the order of less than a meter for the military codes). However, the units are relatively expensive ($200-$500), and in certain instances may not be able to lock-on to these signals due to terrain blockages and satellite radiated power limitations.

[0009] Other geo-location techniques include use of the amplitudes of received signals to determine location (analogous to TOA) and inertial guidance. Amplitude techniques suffer from amplitude variations induced by weather, absorption and fading due to phase interference. Advances in inertial guidance systems and miniaturization of associated components designed to measure accelerations, made principally by the Department of Defense (DoD) funded programs, will enable accurate geo-location to be made. However the systems are currently extremely expensive ($20K-$80K for avionics-quality miniaturized gyros/accelerometers) and would require frequent reference updates (of the order of hours to days) depending on the accuracy errors tolerable.

[0010] The following references are provided as examples of prior art (all refer to systems that are at least twenty years old since conception and description in the open literature), and no claim to novelty is made for these techniques, individually, as they are described in the literature:

[0011] Transit—Johns Hopkins APL TECHNICAL DIGEST January-March 1981, Volume 2, Number 1.

[0012] LORAN—“LORAN Long Range Navigation”, J. A. Pierce, A. A. McKenzie and R. H. Woodward, McGraw-Hill Book Company, New York, 1948.

[0013] GPS—“Global Positioning System”, Vol #1 NAVIGATION, The Institute of Navigation, 815 Fifteenth Street, Suite 832, Washington, D.C. 20005, 1980:

[0014] “Principle of Operation of NAVSTAR and System Characteristics”, R. J. Milliken and C. J. Zoller

[0015] Argos—“Proceedings For The Eleventh Annual Gulf Of Mexico Information Transfer Meeting”, New Orleans, La., pg. 230-232, November 1990, U.S. Dept. of Interior Minerals and Management Service, Published October 1991, Contract #1435-0001-30499, MMS-91-0040, 524 pages.

[0016] SARSAT—“Ambiguity Resolution For Satellite Doppler Positioning System”, P. Argentiero and J. W. Marini, IEEE Transactions On Aerospace And Electronics,pg 439, May 1979.

[0017] DF—“Principles Of High-Resolution Radar”, August W. Rihaczek, Peninsula Publishing, P.O. Box 867, Los Altos, Calif. 94023,1985.

[0018] It should be noted that the geo-location determination in all of these systems is conceptually simple, relying on the intersection of solid geometry surfaces, a technique which can be found in many high school and college freshman-level text books.

[0019] The invention described herein is different in several respects from prior art described above. First, the geo-location techniques utilized are novel and different from those currently in use. All present systems, with the exception of single-pass Doppler, are “receive-only”, versus “transmit-only” for this invention. The device takes advantage of the new geo-location techniques through the use of micro-miniaturization and utilizes no power until activated by one of several sequences or methods, and thus requires little or no power until used, making miniature battery power sources feasible. Since it broadcasts infrequently when activated, the energy drain on batteries allows operation over sufficient periods of time to allow rescues to be made. Further, miniaturization is made possible by compact, omni-directional antenna designs, or easily configurable antenna placements, depending on the operational frequency of the unit. The principles used in these novel techniques have also been applied to the use of beacons which emit visible or infra-red (IR) radiation.

[0020] Finally, the use of microprocessor technologies allows the sequence of operations to be activated, allowing micro-miniaturization for concealment to be possible, as well as affordable. The use of this device with satellite constellations allows global coverage, if desirable. Other units can function with cellular or wireless networks for metropolitan area coverage.

[0021] The major components of the micro-miniature emergency geo-location beacon transmitter generated in accordance with the present invention consists of the following functional components and/or subsystems. A critical component is a miniaturized antenna system capable of providing a significant signal strength to satellite, cellular, or wireless radio receiver antennas that can be from 10 to 500 kilometers distant, depending upon the particular receiving system chosen. Several antenna configurations are possible, including multi-turn loops, center- or end-fed dipoles, or dielectric patch antennas. Each type of antenna has a somewhat different angular distribution of radiated power, which must be adjusted to optimize signal in the receiver direction (overhead for satellites and near-horizon for terrestrial receivers). The radiation patterns have been analyzed in textbooks (e.g., “Antenna Theory and Design,” by R. S. Elliott, Prentice-Hall, Inc., 1981) and by judicious selection of size and wavelength, the patterns can be selected for optimum use. For loop antennas with large wavelengths (compared to the antenna radius), the radiated power can be written as

[0022] P=(ka)²ηI²J₁²(kasin θ)n²/(2r²)

[0023] where η is the impedance of free space, n=# of turns and θ=0 points to the zenith. For a=2.5λ, this pattern is lobed, with a strong zenith radiation power. If the wavelength is small compared to the antenna radius, then the radiated power is

[0024] P=(ka)⁴I²η(sin²θ)n²/(32r²)

[0025] and for a=0.05 (λ), this pattern has a null in the zenith direction and maximum strength in the horizontal direction (θ=+90 /−90). For a traveling wave loop (ka=1), the modes “reverse,” with a null in the horizon direction and a maximum in the vertical direction. For linear dipoles, the power can be written $P=2\eta \frac{I^2}{{\left(4\pi r\right)}^2}\left[{\cos }^2\left(\frac{\pi }{2}\right)\cos \theta /{\sin }^2\theta \right]$

[0026] for λ=4 L, and for large wavelengths (2 L<<λ) as $P={\left(\mathrm{kL}\right)}^2\frac{\eta I^2}{\left(4\pi r^2\right)}{\sin }^2\theta$

[0027] In all cases, the radiation resistance R_radcan reveal the total radiated power by the formula $P_{r\mathrm{ad}}=\frac{I^2}{2}{R}_{r\mathrm{ad}}$

[0028] The electric field at a receiving antenna can be determined from the formula $E\left(v/m\right)=\frac{1}{r}\sqrt{P{\mu }_{o}c/\left(2\pi \right)}P=\mathrm{power}\mathrm{in}\mathrm{watts}$

[0029] This allows the current that flows in the receiving antenna to be calculated. The absorbed power also traditionally can be determined from a reciprocity calculation which yields

[0030] P_abs=P_radD D=directivity

[0031] A second key component is the communications processor. This unit provides a local oscillator to provide the carrier frequency for the communications system. Frequencies range from 100 MHz to 20 GHz, depending on the particular system (satellite, cellular, wireless, etc.). In addition, electronic circuitry is present to convert position or location (latitude, longitude) information (if available from GPS) into the proper carrier frequency modulation to encode the data. In addition, communications signal protocols are impressed on the carrier so that the receiving systems can open a phone link to the central monitoring unit. Identification codes also can appear in the carrier for phone link systems. For geo-location systems dependent on determining bearing to the beacon signal and/or Doppler information, only identification data needs to be encoded in the carrier. Also, included are amplifiers and impedance matchers for the antenna interface. This technology is well established; the unique feature in this invention is the micro-miniaturization to reduce the unit to the order of 1-2 centimeters in thickness. This can be accomplished via a full scale integrated circuitry concept, combining all key sub components on one board/chip. Since the unit will only broadcast intermittently, passive heat sinking can be utilized to reduce components to minimum size consistent with low duty cycle heating.

[0032] A third component of the device is a central processor unit (CPU) which controls the operation of the device. The major functions of the CPU include: device status (battery, test sequence, etc.), transfer of data to the communications processor, device activation (mechanical, voice code, timer, etc.), and control of the above device states.

[0033] Typical device status features include battery voltage and current capability (under load) testing at periodic intervals to ensure sufficient energy is available for beacon signals to exceed detection thresholds. In addition, proper operation of key sub components needs to be monitored. This can include presence of carrier signal at the local oscillator, modulation checks, and the presence of appropriate RF power at the antenna. To conserve power, test sequences would be activated manually prior to use.

[0034] The CPU will also transfer data from memory (identification codes and transmission protocols) to the communications processor. Data could include position information if GPS receivers are incorporated. The overall control of the device, via ROM programs resident in the CPU, allows the selection of key operating modes, test/status, and transmit.

[0035] There are several possible modes to activate the device. The simplest can be mechanical pressure (i.e., pushing a button), but care needs to be taken to prevent inadvertent activation by small children. This could include requiring multiple touches within a short time period for the beacon to turn on, which would greatly reduce the “false alarm” rate. Voice activation is another option, in case the holder of the device is bound and unable to reach the device with his hands. Again, a voice sequence that is not susceptible to accidental activation needs to be programmed into the device CPU. A third option is to have the device activate automatically after a fixed time intervat (that could be adjusted by the user). This would be the most reliable method, but would increase the likelihood of heavy beacon traffic in large population areas, and would also require more frequent recharging. Finally, the device could be activated by a transmitted signal to turn on the device for a single up-link broadcast.

[0036] Use of visible or IR beacons makes the implementation of the micro- miniaturization easier, since these radiation sources can be more compact, the radiation patterns are more easily tailored with gradient index lenses, as opposed to relatively large antennas for the radio-frequency spectrum, and the radiation can be easily modulated, in either amplitude or frequency, to encode the needed data and signal patterns.

[0037] FIG. 1. is a diagram showing LEO communication satellites for receiving beacon signals and higher orbit global positioning satellites (GPS).

[0038] FIG. 2. is a diagram showing terrestrial communications networks such as cellular telephone receiver modes or wireless radio data transmission system antennas.

[0039] FIG. 3. shows several different beacon antenna configurations and radiation patterns.

[0040] FIG. 4. shows a mechanical layout of a representative beacon design.

[0041] FIG. 5. shows the logic flow for the CPU program controlling the activation of the communications control hardware and the associated events taking place in the communications processor.

[0042] FIG. 6 shows a composite illustration of prior art in geo-location for triangulation and Doppler techniques found in numerous textbooks and journal references.

[0043] FIG. 7 shows the basic geometry used in implementing a least squares solution for each of the novel geo-location algorithms described and listed in code form in the Appendix, and a notional emit signal protocol.

[0044] FIG. 8 contains basic error determination equations used in estimating geo-location errors for the Doppler technique, along with the estimated error magnitudes for an Iridium-like satellite constellation.

[0045] FIG. 9 illustrates the basic components of the geo-location technique which utilizes visible or IR radiation emitting beacons, and the basic time correlation measurement technique used for all timing signals.

[0046] FIG. 10 illustrates an application of the geo-location technique for use in autonomous docking of two spacecraft requiring multiple geo-location of reference points to provide both relative position and orientation.

[0047] The three major configurations or modes of operation for this micro-miniature emergency geo-location system are shown in FIGS. 1, 2. The modes allowing global coverage are shown in FIG. 1. The less complex mode utilizes a low earth orbit (LEO) satellite constellation, in which the signal emitted by the emergency beacon 2 is detected by several (>2) satellites 4, and the location of the beacon is determined from the Doppler shifts at the receiving satellites and the satellite positions (which are determined from on-board GPS receivers).

[0048] An optional configuration, which could reduce the number of LEO satellites, would utilize a constellation of GPS 6 satellites with a GPS receiver in the emergency beacon, which would then provide location information in the emitted beacon signal and, therefore, would only require one LEO receiving satellite to function. However, this mode requires a larger beacon unit to accommodate the GPS receiver unit. It also requires 30-60 seconds to determine a latitude/longitude data pair and good receipt of GPS signals from at least 3 GPS satellites. In high population areas, a third mode of operation is possible, which utilizes a network of local receivers or nodes 8, as shown in FIG. 2.

[0049] These nodes can be cellular telephone modes, or can be wireless antennas operating in a slightly different bandwidth. The cellular nodes are typically a few miles apart, whereas the wireless RF/radio links can receive signals up to 30 kilometers from the transmitting unit. Both systems rely on direction finding techniques to determine the location of the emergency beacon signal.

[0050] An operations center, or central ground monitoring station 3, which is comprised of a satellite or network communications link, a monitoring and processing computer (for geo-location calculation, display of location on computer-stored GIS maps, and telephone lists of with user identification codes and emergency phone numbers), and a conventional modem/telephone link with automated dialing equipment to commercial telephone networks, is shown in FIGS. 1 and 2. The operations center is designed to receive information from all satellites 4, or cellular nodes 8, that pick up a signal from the emergency beacon 2, and signals from all satellites or nodes that were able to receive the beacon signal are either down-linked directly or relayed via cross-link(s) to a satellite in position to down-link to the central monitoring station. Data received by this central ground station with the appropriate identification code is automatically collected by monitoring computer 5, which then processes all data available, using Doppler, triangulation, or other techniques described herein (see Appendix with computer programs) or, if necessary, well known in the art, to provide a location for the emergency beacon. The location latitude/longitude are forwarded by the monitoring computer using a conventional telephone/modem system 7, to pre-determined private security personnel or to the nearest local police via a pre-determined/stored message.

[0051] A system that utilizes visible or IR signals for geo-location is shown in FIG. 9. Satellites with optical/IR detection telescopes 56, or high flying surveillance aircraft with similar detection payloads 58 detect the emissions from missiles 54 or aircraft beacons 52. The time of detection is recorded and the signals are relayed to the operations center or monitoring station 3 where the monitoring computer 5 processes the signals and calculates the location of the emitters.

[0052] The link budget margins for reception of the signal emitted by the beacon will determine the design of the beacon antenna and power supply. There are three antenna configurations that can be utilized, depending on the receiver network configurations shown in FIGS. 1, 2. The three options shown in FIG. 3 are dipole, loop, and dielectric patch antennas. The dipole 10 will be more easily adapted to systems above a few GHz, since the half-wave size will be compatible with very small size (few centimeters). They can also be used at lower frequencies, or be sewn onto clothing. Below a few GHz, multi-turn loop antennas 12 can provide sufficient power at sizes small compared to wavelength. Patch antennas 14 provide a broader radiation pattern, but lower radiated power. An estimate of power levels at receivers can be made using a half-wave dipole, using the previously provided formulas. For a 1 GHz carrier, the power density for a 45 degree angle from zenith is $P\left(45°\right)=P_{rad}/r^2\left[\frac{{\cos }^2\left(\cos 45\right)\left(\pi /2\right)}{{\sin }^{2}45}\right]\approx \frac{0.65P_{rad}}{r^2}$

[0053] For a slant path to an orbiting satellite of r=10⁶m, the power density per watt of radiated power will be

[0054] P=0.65E-12 w/m²per watt

[0055] Assuming the receiver has radiation resistance R_rad≈(0.707) 73, then antenna voltages of 2.5×10⁻⁶v/m and currents of the order of 5×10⁻⁸amps would result per watt. For the 5 watt FCC limits, this would result in 0.25 microamps, which should be detectable by the satellite pre-amplifiers. Allowances must be made for obstructions and scattering/absorption, which will reduce received signal strengths. For local networks, such as cellular and RF/wireless receivers, the signal power has already been demonstrated to function adequately (since signals are about 10³times higher than the satellite case). A quarter wave monopole provides reduced attenuation. As is standard in the art, the particular antenna configuration selected will be impedance-matched to the carrier local oscillator, to ensure maximum transfer of oscillator power output to the antenna.

[0056] Battery technology will allow use of miniature Hg/Zn batteries, which approach energy densities of 50 watt-hr/lb, sufficient for intermittent operation over long (hours) periods of time, for weights of several ounces. They also have long storage life (>4 years) when the device is not being used.

[0057] The central processor unit (CPU) and communications processor will utilize miniaturized versions of current cellular telephone technologies. However, since receive capability is not needed (unless command activated), the circuitry can be reduced in complexity, size, and power. When GPS accuracy is needed, an optional commercial GPS unit will be embedded, with simplifications based on elimination of display characteristics and potentially reduced numbers of satellite signals necessary for geo-location.

[0058] A mechanical layout of a representative beacon design is shown in FIG. 4. The beacon can be placed in pockets, strapped to legs or wrists, or be Velcro-attached or sewn into clothing. The device can be disguised as a watch, if worn in an open/visible location. The batteries 16 power the CPU 18 and communications processor 20 (with local oscillator). The antenna 22 is shielded from the electronics by a metal ground plane 24 which enhances radiated power. A GPS unit 26 is optional. A test sequence button 28 allows the unit to be tested for readiness, and also serves to provide activation with multiple pushes.

[0059] The CPU 18 handles control of the device, voice processing for activation (if included), and data set-up for the communications processor (ID., phone protocol, and location data (optional, if applicable)). FIG. 5 contains a top-level flow diagram for the CPU functions that control the communications aspects of the beacon. The CPU provides the signal protocol, which requires an intermittent single burst coded signal pattern. This notional pattern, shown at the bottom of FIG. 7, consists of a sub-pattern consisting of a time period of pure carrier signal (for phase loop lock-in reception purposes), a period of modulated carrier with the detection protocol required by the communications satellite, another period of modulated carrier with the user identification code embedded, a period of pseudo-random pulse modulation (for correlation and transit time extraction), and a final period of pure carrier signal. This pattern is repeated n times in a burst (where n is determined by a high probability that one sub-pattern signal will fall into the acquisition gate of the satellite) to insure that the signal is received by the satellite and forwarded as a normal message. The pattern of n repetitions can be intermittently broadcast at predetermined intervals, or activated by either the user carrying the beacon or via an activation command transmitted to the beacon (in which case the emit burst can be synchronized with the satellite receive window).

[0060] A computer program logic flow to control the activation, operation and test for the emergency beacon is shown in FIG. 5. The CPU is tied to the communications processor 20, which in turn is connected to the antenna unit 22 via a multiplexor contained within the communications processor to control the receive mode (only used to detect the presence of an optional remote activation signal), or the transmit mode. The CPU also performs test functions when activated.

[0061] The novel geo-location techniques claimed in this invention can be segmented into two categories, which relate to the hardware implementation; radiofrequency (RF) and electro-optical (EO/IR). However, the mathematical techniques are almost identical. All of the techniques described below rely on implementation of least square methods, which can be found in many college-level mathematics textbook or mathematics handbooks. For GPS-type and other signal transit time systems (e.g., TOA and TDOA) which utilize a matrix solution approach, the least squares entity that is subsequently operated on by partial derivatives (this differentiated expression is then set to zero for a minimization solution) to extract the equation for the best fit geo-location parameters can be written as $T=\sum _{i=1}^4{\left[\left(\stackrel{\_}{D_i}\right)-{\rho }_i\right]}^2$

[0062] where ρ_irepresents the measured parameter (pseudo-range, transit time, etc.) and {overscore (D_i+L )}={overscore (S_i+L )}−{overscore (X_i+L )} is the vector from the transmit location to the receiver satellite. The geometry is shown in FIG. 7 for both techniques. For Doppler techniques, the least square entity can be written as $\sum _{i=1}^4{\left[\left(\stackrel{\_}{D_i}\right)·\left(\stackrel{\_}{V_{T_i}}\right)-c\left|D_i\right|\left(\beta f_{m_i}-1\right)\right]}^2$

[0063] where βf_mis the ratio of the measured frequency divided by the nominal transmit frequency, c is the speed of light, and V_Tis the total velocity vector (satellite velocity plus the transmitter velocity due to the Earth rotation). Taking the partial derivatives of the above function with respect to the x, y and z coordinates yields the matrix equations for determining the best fit for the location of the transmitter. The matrix equations are shown coded in MathCad in the code listing section that follows. The above techniques do not require knowledge of the Earth terrain to determine a solution, provided that four measured independent parameters are available. This can be accomplished with as few as two satellites.

[0064] A more straight-forward approach can be used if the Earth terrain is known in the area of the transmit beacon. This approach utilizes a least square minimization of the difference between calculated and measured parameters (e.g., distance or time interval), determined over a grid which contains terrain altitude (above the Earth geoid) at each grid point. Due to Doppler effects from the rotating Earth, it should be possible to geo-locate with Doppler and TOA signals from one satellite. Alternatively, if accuracy can be relaxed, simple quadratic solid geometry equations (representing spheres, paraboloids and hyperboloids that can be generated from the measured parameters) can be used to solve for the transmit location with three independent parameters, as is frequently depicted in high school and college mathematics books showing how to solve for the intersection of simple quadratic surfaces. FIG. 6 illustrates the solid geometry approach as it relates to prior art references. The program for this grid-search minimization method is also shown in an Appendix below.

[0065] For geo-location determination of optical sources, the above techniques can be used after the basic parameter (time interval) is extracted from the measurement. In the case of non-cooperative sources (ones whose content cannot be programmed), a traditional radio-astronomy technique, called auto-correlation, is utilized to determine the time difference of signals detected at three or more satellites. A simple description of the technique is shown in the bottom of FIG. 9. The correlation integral is equivalent to “sliding” the transmitted wave form over the received wave form in small steps and summing the product of the two. amplitudes. This sum peaks when the two sets of signals perfectly overlap the amount of time shift necessary gives the true time interval between transmission and reception of the signal. For the GPS/TDOA/TOA systems described above, the correlation integral can be done digitally (with a DSP signal processing chip) off-line at the monitoring center. The auto-correlation determines the time difference between signal arrivals to much greater accuracy than the clock sampling time intervals recorded with the signals by the satellite optical detectors. Most optical emissions, for example from a missile exhaust plume, are sufficiently random (i.e., the intensity as a function of time exhibits random fluctuations) so as to have an effectively infinite “repeat interval”, so that any time ambiguity from a pattern “alias” (such as the relatively short “C/A code” in a GPS system) should not occur.

[0066] Systems such as those described above can be used for cooperative high flying aircraft (e.g., for an FAA global passive tracking system), provided that the optical beacon has a sufficiently long pulse repetition interval. Detection can done with standard solid state detector focal plane arrays (FPAs) using visible or IR semiconductor materials (e.g., silicon compounds for visible, HgCdTe for IR). Sufficient number of pixels are necessary (e.g., >1024×1024) to resolve closely spaced transmitting sources, which should be feasible for LEO satellites. The detailed design principles for such a detection system have been described by this inventor in U.S. Pat. No. 5,604,595, and thus will not be repeated in this specification. Since the geo-location calculation does not depend on precise direction measurements, the detection telescope alignment is not critical, which should lower satellite costs. In addition, narrow-band sources and detectors can be utilized, so that background signals from the ground or atmospheric phenomenon can be effectively suppressed. A notional system is shown in FIG. 9.

[0067] It should be noted that the very precise geo-location achievable with current GPS and DGPS systems will not likely be necessary for emergency rescue applications. Geo-location errors no less than of the order of 100 meters are required. Thus tolerances on components may be relaxed to achieve micro-miniaturization, and satellite overlap coverage constraints may also be relaxed, allowing for less costly space systems. Precise error calculations have been programmed for the various techniques described, including errors from atmospheric effects, satellite state vector errors and instrument errors such as frequency and clock precision, and appear in a commercially available program (“GLM”, copyright May 1998). Initial error calculations for a Doppler system using an Iridium-like satellite constellation are shown in FIG. 8.

[0068] It should be noted that for use of both the RF or visible/IR spectrum implementation of this invention, an alternate embodiment consistent with the principles of operation described herein and which can reduce the size and complexity of the beacon transmitter unit carried by a user, can operate as follows. For an RF system, the local oscillator or radiation source contained in the beacon unit can be removed, and the device can operate as an amplifier and re-transmission device, or transponder, which receives a signal from an off-user beacon unit at a distant location. This would allow a system at a central location(s) to provide the appropriate user ID encoded in the requisite pseudo-random signal to automatically trigger a geo-location determination, as opposed to in effect sending a request for geo-location, which then has to be processed by the beacon unit, and a beacon signal generated. This signal could be provided by one or more of the communications receiving nodes, which have to be capable of communications with the beacon in the preferred embodiment (to activate and/or receive a geo-location signal). In addition, the beacon unit no longer has to have the power and signal processing capability to generate a geo-location signal; it just has to amplify a received signal (if it has the correct ID) and re-transmit it, thus reducing power, weight and size constraints on a device which for some important applications, requires small, stealthy beacon units.

[0069] For a visible/IR system, the optical signal does not have to originate with the object to be located (e.g., missile plume or aircraft beacon). An alternate embodiment could consist of an optical source (e.g., a laser with telescope or “optical antenna”) which illuminates the object to be geo-located with the requisite pseudo-random optical signal, which is then reflected off of the object and in effect operates as a transponder (without amplification however). This embodiment could function for non-cooperative targets when no appropriate pseudo-random emissions are available for geo-location. The reflection off the object to be located can be from the natural reflectivity of the object surface, or for cooperative targets, can be enhanced by using retro-reflectors. Such a system can be applied to autonomous docking systems to resupply/repair satellites on orbit. For example, the satellites would have multiple retro-reflectors (each with slightly different frequencies “notches” to prevent ambiguity in the source of the received signals) at known locations relative to the docking mechanism, and the docking vehicle would have multiple laser transmitter(s) and receivers at known locations. With three or four sources and retro-reflectors the relative position and orientations of the two spacecraft could be determined and used to control thrusters to accomplish autonomous docking.

[0070] Another alternative embodiment of the invention can be utilized for geo-location in dense urban environments, in which beacon signal multi-pathing effects (e.g., reflections off of building surfaces) can lead to larger than desired geo-location errors). This technique would utilize an initial geo-location using the preferred embodiment, followed by a refined geo-location process utilizing the grid minimization technique shown in the Computer Code Appendix, which would include calculation and correction of multi-pathing effects. The topology of the buildings could be obtained from a computer file based on building blueprints, or created by photogrammetry of satellite or aerial images of the urban region.

[0071] MELS Geoloca Program

[0072] Case #1: Satellites With Doppler And Triangula (TOA)

[0073] Date: Sep. 28, 1996

[0074] Summary: This program is a prototype for the general satellite case and calculates position location from Doppler and TOA data from only one satellite. The read-in topography file is not yet implemented.

[0075] Program constants: $\mathrm{RadiusEarth}:=6378$ ${\mathrm{AltGrid}}_{1,1}:=0.100$ $\mathrm{latitude}:=38.8$ $\mathrm{longitude}:=-77.0$ $\mathrm{omegaEarth}:=7.29·{10}^{-5}$ $\mathrm{timesec}:=0$ $\mathrm{timemin}:=0$ $\mathrm{timehour}:=0$ $\mathrm{time}:=\mathrm{timehour}·3600+\mathrm{timemin}·60+\mathrm{timesec}$ $\mathrm{SatPosVect}:=\left(\begin{array}{c}0\\ 0\\ 7028\end{array}\right)$ $\mathrm{SatVelVect}:=\left(\begin{array}{c}0\\ 7.2\\ 0\end{array}\right)$ $c:=2.997930·{10}^5$ $\mathrm{freq}:=240.50·{10}^6$ $\mathrm{vorb}:=\left|\mathrm{SatVelVect}\right|\mathrm{deltafreq}:=4804.885115$ $\mathrm{toadif}:=19420.602·{10}^{-6}$ $\mathrm{xvec}:=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$ $\mathrm{yvec}:=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$ $\mathrm{zvec}:=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$ $\mathrm{rho}:=\mathrm{RadiusEarth}+{\mathrm{AltGrid}}_{1,1}$ $\mathrm{phi}:=\mathrm{longitude}·\frac{\pi }{180}+\mathrm{time}·\mathrm{omegaEarth}\mathrm{time}·\mathrm{omegaEarth}=0$ $\mathrm{theta}:=\frac{\left(90.000-\mathrm{latitude}\right)}{90.000}·\frac{\pi }{2}\mathrm{phi}=-1.343903524$ $\mathrm{RotMat}:=\left(\begin{array}{ccc}\mathrm{rho}·\sin \left(\mathrm{theta}\right)·\cos \left(\mathrm{phi}\right)& 0& 0\\ 0& \mathrm{rho}·\sin \left(\mathrm{theta}\right)·\sin \left(\mathrm{phi}\right)& 0\\ 0& 0& \mathrm{rho}·\cos \left(\mathrm{theta}\right)\end{array}\right)$ $X:=\mathrm{RotMat}·\mathrm{xvec}$ $Y:=\mathrm{RotMat}·\mathrm{yvec}$ $Z:=\mathrm{RotMat}·\mathrm{zvec}$ $X_0=1.118163188·{10}^3$ $Y_1=-4.843296871·{10}^3$ $Z_2=3.996541769·{10}^3$ $\begin{array}{ccc}\begin{array}{c}\mathrm{distrav}:=c·\mathrm{toadif}\\ \mathrm{costhDop}:=[\left(\frac{\mathrm{deltafreq}}{\mathrm{freq}}\right)·\frac{c}{\mathrm{vorb}}\\ \mathrm{electden}:=1·{10}^6\\ \mathrm{plasmaf}:=8.98·{10}^3·\sqrt{\mathrm{electden}}\\ \mathrm{plasmaf}=8.98·{10}^6\\ \mathrm{costhinc}:=\frac{\pi }{100}\end{array}& \mathrm{edenvect}:=\left[\begin{array}{c}0.2·{10}^6\\ 1.5·{10}^6\\ 1.3·{10}^6\\ 0.6·{10}^6\\ 0.4·{10}^6\\ 0.2·{10}^6\\ 0.1·{10}^6\\ 0.05·{10}^6\\ 0.03·{10}^6\\ 0.02·{10}^6\end{array}\right]& \mathrm{dh}:=100·{10}^3\end{array}$ $\begin{array}{cc}\mathrm{indexref}:=\sqrt{1-{\left(\frac{\mathrm{plasmaf}}{\mathrm{freq}}\right)}^2}& \mathrm{TEC}:=\left(\sum \mathrm{edenvect}\right)·\frac{\mathrm{dh}·{10}^6}{\cos \left(\mathrm{costhinc}\right)}\\ \mathrm{indexref}=0.999302661& \mathrm{TEC}=4.402172206·{10}^{17}\\ \mathrm{freq2}:=2.4& \mathrm{tiondelay}:=\frac{40.3}{c·{10}^3·{\mathrm{freq}}^2}·\mathrm{TEC}\\ \mathrm{freq1}:=\frac{\mathrm{freq}}{{10}^9}& \mathrm{tiondelay}=1.023105511·{10}^{-6}\\ \mathrm{toadif2}:=19400·{10}^{-6}& \mathrm{tatmdelay}:={\left(\frac{1}{\left|\frac{1}{c}-\frac{1.003}{c}\right|·8}\right)}^{-1}\\ \mathrm{kslope}:={\left[\left|\frac{\left(\frac{1}{{\mathrm{freq2}}^2}-\frac{1}{{\mathrm{freq1}}^2}\right)}{\left(\mathrm{toadif2}-\mathrm{toadif}\right)}\right|\right]}^{-1}& \mathrm{tatmdelay}=8.005523811·{10}^{-8}\end{array}$ $\mathrm{kslope}=1.203712159·{10}^{-6}$ $\mathrm{delta0}:=\frac{\mathrm{kslope}}{2}·\frac{\left|\frac{1}{{\mathrm{freq}}^2}+\frac{1}{{\mathrm{freq2}}^2}\right|}{\mathrm{toadif}+\mathrm{toadif2}}$ $\mathrm{delta0}=2.69158378·{10}^{-6}$ $\mathrm{grpvel}:=\frac{c}{\left[\mathrm{indexref}+\frac{1}{\mathrm{indexref}}·{\left(\frac{\mathrm{plasmaf}}{\mathrm{freq}}\right)}^2\right]}$ $\mathrm{grpvel}=2.995839426·{10}^5$ $\mathrm{GeoVect}:=\left(\begin{array}{c}{X}_0\\ Y_1\\ Z_2\end{array}\right)$ $\mathrm{DifVec}:=\mathrm{GeoVect}-\mathrm{SatPosVect}$ $\mathrm{cosang}:=\mathrm{GeoVect}·\frac{\mathrm{DifVec}}{\left(\left|\mathrm{GeoVect}\right|·\left|\mathrm{DifVec}\right|\right)}\mathrm{Ionosphere}\mathrm{angle}$ $\mathrm{ang1}:=\mathrm{acos}\left(\mathrm{cosang}\right)$ $\mathrm{ang}:=180-\frac{180}{\pi }·\mathrm{ang1}\mathrm{toaest}:=\frac{\left|\mathrm{DifVec}\right|}{c}$ $\mathrm{ang}=109.822404269\mathrm{toaest}=0.019420602$ $\mathrm{cosang2}:=\mathrm{DifVec}·\frac{\mathrm{SatVelVect}}{\left(\left|\mathrm{DifVec}\right|·\left|\mathrm{SatVelVect}\right|\right)}\mathrm{Doppler}\mathrm{angle}$ $\mathrm{ang2}:=\mathrm{acos}\left(\mathrm{cosang2}\right)·\frac{180}{\pi }\mathrm{ang2}=146.291602191\mathrm{cosang2}=-0.8318729$ $\mathrm{iondepth}:=400$ $\mathrm{delta0p}:={[\frac{{\left(\frac{\mathrm{iondepth}}{\cos \left(\mathrm{ang}·\frac{\pi }{180}\right)}\right)}^{-1}}{\left(\frac{1}{c}-\frac{1}{\mathrm{grpvel}}\right)}]}^{-1}$ $\mathrm{delta0p}=2.745681785·{10}^{-6}$

[0076] This section calculates the differences between measured and actual values for the TOA and Doppler angle parameters: the sum of the two differences is calculated for each point on the earth surface grid, and the minimum sum determines the location coordinates. An altitude profile file must be input for the general area around the estimated location.

[0077] latstart:=38.75

[0078] longstart:=−77.15

[0079] initheta:=90−latstart

[0080] initphi:=longstat

[0081] deltaphi:=0.01

[0082] delta:=−0.01

[0083] oldmin:=10⁹

[0084] thetasol:=0

[0085] phisol:=0 $\mathrm{coscalc}:=\left(\frac{\mathrm{deltafreq}}{\mathrm{freq}}\right)·\frac{c}{\mathrm{vorb}}$

[0086] coscalc=0.83187279

[0087] cosDop:=coscalc

[0088] cosDop=0.83187279

[0089] Dopang:=acos(cosDo)·180/π

[0090] Dopage=33.708397799

[0091] vorb=7.2

[0092] velmag:={square root over (SatVelVecta·SatVelVect)}

[0093] valmag=7.2

[0094] alt:=AltGrid_1.1

[0095] Start of loops

[0096] i:=0.10

[0097] j:=0.35

[0098] thetagrid_i,j:=((−i)·deltatheta+inittheta)·180/π

[0099] phigrid_i,j:=((j)·deltaphi+initphi)·180/π+time·omegaEarth

[0100] xg_i,j:=RadiusEarth·sin(thetagrid_i,j)·cos(phigrid_i,j)

[0101] yg_i,j:=RadiusEarth·sin(thetagrid_i,j)·sin(phigrid_i,j)

[0102] zg_i,j:=radiusEarth·cos(thetagrid_i,j)

[0103] Take into account the velocity of the Earth's rotation ${\mathrm{losx}}_{i,j}:={\mathrm{SatPosVect}}_0-{\mathrm{xg}}_{i,j}{\mathrm{losy}}_{i,j}:={\mathrm{SatPosVect}}_1-{\mathrm{yg}}_{i,j}{\mathrm{losz}}_{i,j}:={\mathrm{SatPosVect}}_2-{\mathrm{zg}}_{i,j}$ ${\mathrm{omegx}}_{i,j}:={\mathrm{yg}}_{i,j}{\mathrm{omegy}}_{i,j}:={\mathrm{xg}}_{i,j}$ ${\mathrm{omegv}}_{i,j}:=\left(\mathrm{RadiusEarth}+\mathrm{alt}\right)·\mathrm{omegaEarth}·\frac{\left({\mathrm{omegx}}_{\mathrm{ij}}·{\mathrm{losx}}_{i,j}+{\mathrm{omegy}}_{i,j}·{\mathrm{losy}}_{i,j}\right)}{\sqrt{{\left({\mathrm{losx}}_{i,j}\right)}^2+{\left({\mathrm{losy}}_{i,j}\right)}^2+{\left({\mathrm{losz}}_{i,j}\right)}^2}·\sqrt{{\left({\mathrm{omegx}}_{i,j}\right)}^2+{\left({\mathrm{omegy}}_{i,j}\right)}^2}}$ ${\mathrm{satvel}}_{i,j}:=\frac{{\mathrm{SatVelVect}}_0·{\mathrm{losx}}_{i,j}+{\mathrm{SatVelVect}}_1·{\mathrm{losy}}_{i,j}+{\mathrm{SatVelVect}}_2·{\mathrm{losz}}_{i,j}}{\sqrt{{\left({\mathrm{losx}}_{i,j}\right)}^2+{\left({\mathrm{losy}}_{i,j}\right)}^2+{\left({\mathrm{losz}}_{i,j}\right)}^2}}$

[0104] vrel_i,j:=satvel_i,j+omegv_i,j

[0105] ${\mathrm{cosDop}}_{i,j}:=\left(\frac{\mathrm{deltafreq}}{\mathrm{freq}}\right)·\frac{c}{{\mathrm{vrel}}_{i,j}}$ ${\mathrm{min1}}_{i,j}:=\sqrt{{\left[\left({\mathrm{SatPosVect}}_0\right)-\left({\mathrm{xg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{xg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]}^2+{\left[\left({\mathrm{SatPosVect}}_1\right)-\left({\mathrm{yg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{yg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]}^2+{\left[\left({\mathrm{SatPosVect}}_2\right)-\left({\mathrm{zg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{zg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]}^2}-c\mathrm{toadif}$ ${\mathrm{min2}}_{i,j}:=\sqrt{{\left[\left[\left({\mathrm{SatPosVect}}_0\right)-\left({\mathrm{xg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{xg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]·{\mathrm{SatVelVect}}_0\right]}^2+{\left[\left[\left({\mathrm{SatPosVect}}_1\right)-\left({\mathrm{yg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{yg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]·{\mathrm{SatVelVec}}_1\right]}^2+{\left[\left[\left({\mathrm{SatPosVect}}_2\right)-\left({\mathrm{zg}}_{i,j}+\mathrm{alt}·\frac{{\mathrm{zg}}_{i,j}}{\mathrm{RadiusEarth}}\right)\right]·{\mathrm{SatVelVect}}_2\right]}^2}-|{\mathrm{min1}}_{i,j}+c\mathrm{toadif}|·\left|{\mathrm{vrel}}_{i,j}\right|·\left|{\mathrm{cosDop}}_{i,j}\right|{\mathrm{mintot}}_{i,j}:=\sqrt{{\left({\mathrm{min1}}_{i,j}\right)}^2}+\sqrt{{\left({\mathrm{min2}}_{i,j}\right)}^2}$

[0106] M_i,j:=mintot_i,j

[0107] Geo:=min(M)

[0108] Geo=7.281331364·10⁻⁵

[0109] k:=0 . . . 10

[0110] l:=0 . . . 35

[0111] Valx_k,l:=if(M_k,l=Geo, k, 0)

[0112] mthval:=max(Valx)

[0113] Valy_k,l:=if(M_k,l=Geo, 1, 0)

[0114] mphval:=max(Valy)

[0115] mthval=5

[0116] mphval =15

[0117] min1_5,15=−6.858075085·10⁻⁵

[0118] min2_5,15=−4.232562787·10⁻⁶

[0119] M_5,15=7.281331364·10⁻⁵

[0120] geolat:=mthval·deltatheta+latstart

[0121] geolon:=mphval·deltaphi+longstart

[0122] geolat=38.8

[0123] geolon=−77

[0124] This program calulates geolocation from four satellites using the A Z measured Dopple shifts at the satellites and an initial estimate of the location. It is directly analogous to the GPS approach. This mod replaces measured values in a matrix with calculated values and is identical to the original DPS program (Aug. 31, 1997). $\mathrm{freq0}:=240.50·{10}^6$ $\mathrm{lat}:=38.8$ $\mathrm{lon}:=-77.0$ $\mathrm{omegaEarth}:=0.0\mathrm{oE}:=7.29·{10}^{-5}$ $\mathrm{timesec}:=0$ $\mathrm{timemin}:=0$ $\mathrm{timehour}:=0$ $\mathrm{time}:=\mathrm{timehour}·3600+\mathrm{timemin}·60+\mathrm{timesec}$ $\begin{array}{cccc}\mathrm{SatPosVect1}:=\left(\begin{array}{c}0\\ 0\\ 7028\end{array}\right)& \mathrm{SatPosVect2}:=\left(\begin{array}{c}0.0\\ 4969.5465\\ 4969.5465\end{array}\right)& \mathrm{SatPosVect3}:=\left(\begin{array}{c}0\mathrm{.0}\\ -4969\mathrm{.5465}\\ 4969\mathrm{.5465}\end{array}\right)& \mathrm{SatPosVect4}:=\left(\begin{array}{c}4969.5465\\ 0.0\\ 4969.5465\end{array}\right)\\ \mathrm{SatVelVect1}:=\left(\begin{array}{c}0\\ 7.2\\ 0\end{array}\right)& \mathrm{SatVelVect2}:=\left(\begin{array}{c}0\\ 5.0911688\\ -5.0911688\end{array}\right)& \mathrm{SatVelVect3}:=\left(\begin{array}{c}0\\ 5.0911688\\ 5.0911688\end{array}\right)& \mathrm{SatVelVect4}:=\left(\begin{array}{c}5\mathrm{.0911688}\\ 0.0\\ -5\mathrm{.0911688}\end{array}\right)\end{array}$ $c:=2.997930·{10}^5$ $\mathrm{RadiusEarth}:=6378$ $\mathrm{alt}:=0.1$ $\mathrm{xvec}:=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$ $\mathrm{yvec}:=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$ $\mathrm{zvec}:=\left(\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 1\end{array}\right)$ $\mathrm{rho}:=\mathrm{RadiusEarth}+\mathrm{alt}$ $\mathrm{phi}:=\mathrm{lon}·\frac{\pi }{180}+\mathrm{time}·\mathrm{omegaEarth}\mathrm{time}·\mathrm{omegaEarth}=0$ $\mathrm{theta}:=\frac{\left(90.000-\mathrm{lat}\right)}{90.000}·\frac{\pi }{2}\mathrm{phi}=-1.344$ $\mathrm{RotMat}:=\left(\begin{array}{ccc}\mathrm{rho}·\sin \left(\mathrm{theta}\right)·\cos \left(\mathrm{phi}\right)& 0& 0\\ 0& \mathrm{rho}·\sin \left(\mathrm{theta}\right)·\sin \left(\mathrm{phi}\right)& 0\\ 0& 0& \mathrm{rho}·\cos \left(\mathrm{theta}\right)\end{array}\right)$ $X:=\mathrm{RotMat}·\mathrm{xvec}$ $Y:=\mathrm{RotMat}·\mathrm{yvec}$ $Z:=\mathrm{RotMat}·\mathrm{zvec}$ $X=\left(\begin{array}{c}1.118163188·{10}^3\\ 0\\ 0\end{array}\right)$ $Y=\left(\begin{array}{c}0\\ -4.843296871·{10}^3\\ 0\end{array}\right)$ $Z=\left(\begin{array}{c}0\\ 0\\ 3.996541769·{10}^3\end{array}\right)$ $\mathrm{XP}:=\left(\begin{array}{c}{X}_0\\ Y_1\\ Z_2\end{array}\right)$ $\begin{array}{cc}\mathrm{Ve}:=\mathrm{XP}\times \mathrm{zvec}·\left(\mathrm{omegaEarth}\right)& \mathrm{fm1}:=2.40504804885·{10}^8\\ \mathrm{Vtl}:=\mathrm{Ve}+\mathrm{SatVelVect1}& \mathrm{fm2}:=2.40503637992·{10}^8\\ \mathrm{Vt2}:=\mathrm{Ve}+\mathrm{SatVelVect2}& \mathrm{fm3}:=2.405023247783·{10}^8\\ \mathrm{Vt3}:=\mathrm{Ve}+\mathrm{SatVelVect3}& \mathrm{fm4}:=2.405018767603·{10}^8\\ \mathrm{Vt4}:=\mathrm{Ve}+\mathrm{SatVelVect4}& \mathrm{freqdelm1}:=\frac{\mathrm{fm1}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm2}:=\frac{\mathrm{fm2}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm3}:=\frac{\mathrm{fm3}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm4}:=\frac{\mathrm{fm4}-\mathrm{freq0}}{\mathrm{freq0}}\end{array}$ $\mathrm{D1}:=\sqrt{{\left({\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D2}:=\sqrt{{\left({\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D3}:=\sqrt{{\left({\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D4}:=\sqrt{{\left({\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\right)}^2}$ $\begin{array}{cc}\mathrm{D1v}:=\left(\begin{array}{c}{\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\end{array}\right)& \\ \mathrm{D2v}:=\left(\begin{array}{c}{\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs1}:=\frac{\mathrm{D1v}·\mathrm{Vt1}}{\left|\mathrm{Vt1}\right|·\left|\mathrm{D1v}\right|}\\ \mathrm{D3v}:=\left(\begin{array}{c}{\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs2}:=\frac{\mathrm{D2v}·\mathrm{Vt2}}{\left|\mathrm{Vt2}\right|·\left|\mathrm{D2v}\right|}\\ \mathrm{D4v}:=\left(\begin{array}{c}{\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs3}:=\frac{\mathrm{D3v}·\mathrm{Vt3}}{\left|\mathrm{Vt3}\right|·\left|\mathrm{D3v}\right|}\\ & \mathrm{cs4}:=\frac{\mathrm{D4v}·\mathrm{Vt4}}{\left|\mathrm{Vt4}\right|·\left|\mathrm{D4v}\right|}\end{array}$ $A:=\left(\begin{array}{cccc}{\mathrm{Vt1}}_0-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_0}{\mathrm{D1}}& {\mathrm{Vt1}}_1-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_1}{\mathrm{D1}}& {\mathrm{Vt1}}_2-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_2}{\mathrm{D1}}& c·\mathrm{D1}·\left(\frac{\left|\mathrm{Vt1}\right|·\frac{\mathrm{cs1}}{c}+1}{\mathrm{freq0}}\right)\\ {\mathrm{Vt2}}_0-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_0}{\mathrm{D2}}& {\mathrm{Vt2}}_1-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_1}{\mathrm{D2}}& {\mathrm{Vt2}}_2-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_2}{\mathrm{D2}}& c·\mathrm{D2}·\left(\frac{\left|\mathrm{Vt2}\right|·\frac{\mathrm{cs2}}{c}+1}{\mathrm{freq0}}\right)\\ {\mathrm{Vt3}}_0-\left|\mathrm{Vt3}\right|·\mathrm{cs3}·\frac{{\mathrm{D3v}}_0}{\mathrm{D3}}& {\mathrm{Vt3}}_1-\left|\mathrm{Vt3}\right|·\mathrm{cs3}·\frac{{\mathrm{D3v}}_1}{\mathrm{D3}}& {\mathrm{Vt3}}_2-\left|\mathrm{Vt3}\right|·\mathrm{cs3}·\frac{{\mathrm{D3v}}_2}{\mathrm{D3}}& c·\mathrm{D3}·\left(\frac{\left|\mathrm{Vt3}\right|·\frac{\mathrm{cs3}}{c}+1}{\mathrm{freq0}}\right)\\ {\mathrm{Vt4}}_0-\left|\mathrm{Vt4}\right|·\mathrm{cs4}·\frac{{\mathrm{D4v}}_0}{\mathrm{D4}}& {\mathrm{Vt4}}_1-\left|\mathrm{Vt4}\right|·\mathrm{cs4}·\frac{{\mathrm{D4v}}_1}{\mathrm{D4}}& {\mathrm{Vt4}}_2-\left|\mathrm{Vt4}\right|·\mathrm{cs4}·\frac{{\mathrm{D4v}}_2}{\mathrm{D4}}& c·\mathrm{D4}·\left(\frac{\left|\mathrm{Vt4}\right|·\frac{\mathrm{cs4}}{c}+1}{\mathrm{freq0}}\right)\end{array}\right)$ $\begin{array}{cc}\mathrm{GEOM}:={\left(A^T·A\right)}^{-1}·A^T& A=\left[\begin{array}{cccc}1.15& 2.218& -3.119& 7.258\\ 0.511& 0.607& -5.536& 12.371\\ 2.178& 5.337& 3.196& 1.854\\ 3.653& -1.809& -5.455& 7.808\end{array}\right]\\ L:=\left[\begin{array}{c}\mathrm{D1v}·\mathrm{Vt1}-c·\mathrm{D1}·\mathrm{freqdelm1}\\ \mathrm{D2v}·\mathrm{Vt2}-c·\mathrm{D2}·\mathrm{freqdelm2}\\ \mathrm{D3v}·\mathrm{Vt3}-c·\mathrm{D3}·\mathrm{freqdelm3}\\ \mathrm{D4v}·\mathrm{Vt4}-c·\mathrm{D4}·\mathrm{freqdelm4}\end{array}\right]& \mathrm{GEOM}=\left[\begin{array}{cccc}-0.051& -0\mathrm{.138}& 0.118& 0.239\\ 0.614& -0\mathrm{.28}& -0.074& -0.111\\ -0.774& 0\mathrm{.399}& 0.283& 0.019\\ -0.374& 0\mathrm{.279}& 0.125& 0.004\end{array}\right]\end{array}$ $\mathrm{DR}:=\mathrm{GEOM}·L$ $L=\left[\begin{array}{c}8.305·{10}^{-4}\\ -5.465·{10}^{-4}\\ -1.141·{10}^{-4}\\ -1.745·{10}^{-4}\end{array}\right]$ $\mathrm{DR}=\left[\begin{array}{c}-2.201·{10}^{-5}\\ 6.907·{10}^{-4}\\ -8.966·{10}^{-4}\\ -4.784·{10}^{-4}\end{array}\right]$

[0125] This program calulates geolocation from two satellites using the measured transit times and Doppler shifts at the satellites and an initial estimate of the location. It is directly analogous to the GPS approach. This mod replaces measured values in a matrix with calculated values and is identical to the original DPS/TDOA program (Sep. 5, 1997). $\mathrm{freq0}:=240.50·{10}^6$ $\mathrm{lat}:=38.9$ $\mathrm{lon}:=-77.0$ $\mathrm{omegaEarth}:=0.0\mathrm{oE}:=7.29·{10}^{-5}$ $\mathrm{timesec}:=0$ $\mathrm{timemin}:=0$ $\mathrm{timehour}:=0$ $\mathrm{time}:=\mathrm{timehour}·3600+\mathrm{timemin}·60+\mathrm{timesec}$ $\begin{array}{cccc}\mathrm{SatPosVect1}:=\left(\begin{array}{c}0\\ 0\\ 7028\end{array}\right)& \mathrm{SatPosVect2}:=\left(\begin{array}{c}0.0\\ 4969.5465\\ 4969.5465\end{array}\right)& \mathrm{SatPosVect3}:=\left(\begin{array}{c}0\mathrm{.0}\\ -4969\mathrm{.5465}\\ 4969\mathrm{.5465}\end{array}\right)& \mathrm{SatPosVect4}:=\left(\begin{array}{c}4969.5465\\ 0.0\\ 4969.5465\end{array}\right)\\ \mathrm{SatVelVect1}:=\left(\begin{array}{c}0\\ 7.2\\ 0\end{array}\right)& \mathrm{SatVelVect2}:=\left(\begin{array}{c}0\\ 5.0911688\\ -5.0911688\end{array}\right)& \mathrm{SatVelVect3}:=\left(\begin{array}{c}0\\ 5.0911688\\ 5.0911688\end{array}\right)& \mathrm{SatVelVect4}:=\left(\begin{array}{c}5\mathrm{.0911688}\\ 0.0\\ -5\mathrm{.0911688}\end{array}\right)\end{array}$ $c:=2.997930·{10}^5$ $\mathrm{RadiusEarth}:=6378$ $\mathrm{alt}:=0.1$ $\mathrm{xvec}:=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$ $\mathrm{yvec}:=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$ $\mathrm{zvec}:=\left(\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 1\end{array}\right)$ $\mathrm{rho}:=\mathrm{RadiusEarth}+\mathrm{alt}$ $\mathrm{phi}:=\mathrm{lon}·\frac{\pi }{180}+\mathrm{time}·\mathrm{omegaEarth}\mathrm{time}·\mathrm{omegaEarth}=0$ $\mathrm{theta}:=\frac{\left(90.000-\mathrm{lat}\right)}{90.000}·\frac{\pi }{2}\mathrm{phi}=-1.344$ $\mathrm{RotMat}:=\left(\begin{array}{ccc}\mathrm{rho}·\sin \left(\mathrm{theta}\right)·\cos \left(\mathrm{phi}\right)& 0& 0\\ 0& \mathrm{rho}·\sin \left(\mathrm{theta}\right)·\sin \left(\mathrm{phi}\right)& 0\\ 0& 0& \mathrm{rho}·\cos \left(\mathrm{theta}\right)\end{array}\right)$ $X:=\mathrm{RotMat}·\mathrm{xvec}$ $Y:=\mathrm{RotMat}·\mathrm{yvec}$ $Z:=\mathrm{RotMat}·\mathrm{zvec}$ $X=\left(\begin{array}{c}1.116592389·{10}^3\\ 0\\ 0\end{array}\right)$ $Y=\left(\begin{array}{c}0\\ -4.836492993·{10}^3\\ 0\end{array}\right)$ $Z=\left(\begin{array}{c}0\\ 0\\ 4.005211178·{10}^3\end{array}\right)$ $\mathrm{XP}:=\left(\begin{array}{c}{X}_0\\ Y_1\\ Z_2\end{array}\right)$ $\begin{array}{cc}\mathrm{Ve}:=\mathrm{XP}\times \mathrm{zvec}·\left(\mathrm{omegaEarth}\right)& \mathrm{fm1}:=2.40504804885·{10}^8\\ \mathrm{Vtl}:=\mathrm{Ve}+\mathrm{SatVelVect1}& \mathrm{fm2}:=2.40503637992·{10}^8\\ \mathrm{Vt2}:=\mathrm{Ve}+\mathrm{SatVelVect2}& \mathrm{fm3}:=2.405023247783·{10}^8\\ \mathrm{Vt3}:=\mathrm{Ve}+\mathrm{SatVelVect3}& \mathrm{fm4}:=2.405018767603·{10}^8\\ \mathrm{Vt4}:=\mathrm{Ve}+\mathrm{SatVelVect4}& \mathrm{freqdelm1}:=\frac{\mathrm{fm1}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm2}:=\frac{\mathrm{fm2}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm3}:=\frac{\mathrm{fm3}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm4}:=\frac{\mathrm{fm4}-\mathrm{freq0}}{\mathrm{freq0}}\end{array}$ $\mathrm{D1}:=\sqrt{{\left({\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D2}:=\sqrt{{\left({\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D3}:=\sqrt{{\left({\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D4}:=\sqrt{{\left({\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\right)}^2}$ $\begin{array}{ccc}\mathrm{D1v}:=\left(\begin{array}{c}{\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\end{array}\right)& & \\ \mathrm{D2v}:=\left(\begin{array}{c}{\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs1}:=\frac{\mathrm{D1v}·\mathrm{Vt1}}{\left|\mathrm{Vt1}\right|·\left|\mathrm{D1v}\right|}& \mathrm{t1}:=0.019420602\\ \mathrm{D3v}:=\left(\begin{array}{c}{\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs2}:=\frac{\mathrm{D2v}·\mathrm{Vt2}}{\left|\mathrm{Vt2}\right|·\left|\mathrm{D2v}\right|}& \mathrm{t2}:=0.03310337\\ \mathrm{D4v}:=\left(\begin{array}{c}{\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\end{array}\right)& \mathrm{cs3}:=\frac{\mathrm{D3v}·\mathrm{Vt3}}{\left|\mathrm{Vt3}\right|·\left|\mathrm{D3v}\right|}& \mathrm{cs1}=0.832\\ & \mathrm{cs4}:=\frac{\mathrm{D4v}·\mathrm{Vt4}}{\left|\mathrm{Vt4}\right|·\left|\mathrm{D4v}\right|}& \mathrm{cs2}=0.63\end{array}$ $A:=\left(\begin{array}{cccc}\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& 1·\frac{c}{{\mathrm{freq0}}^2}\\ \frac{{\mathrm{D2v}}_0}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D2v}}_1}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D2v}}_2}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& -1·\frac{c}{{\mathrm{freq0}}^2}\\ {\mathrm{Vt1}}_0-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_0}{\mathrm{D1}}& {\mathrm{Vt1}}_1-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_1}{\mathrm{D1}}& {\mathrm{Vt1}}_2-\left|\mathrm{Vt1}\right|·\mathrm{cs1}·\frac{{\mathrm{D1v}}_2}{\mathrm{D1}}& c·\mathrm{D1}·\left(\frac{\left|\mathrm{Vt1}\right|·\frac{\mathrm{cs1}}{c}+1}{\mathrm{freq0}}\right)\\ {\mathrm{Vt2}}_0-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_0}{\mathrm{D2}}& {\mathrm{Vt2}}_1-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_1}{\mathrm{D2}}& {\mathrm{Vt2}}_2-\left|\mathrm{Vt2}\right|·\mathrm{cs2}·\frac{{\mathrm{D2v}}_2}{\mathrm{D2}}& c·\mathrm{D2}·\left(\frac{\left|\mathrm{Vt2}\right|·\frac{\mathrm{cs2}}{c}+1}{\mathrm{freq0}}\right)\end{array}\right)$ $A=\left[\begin{array}{cccc}-0.192& 0.832& 0.52& 5.183·{10}^{-12}\\ 0.08& 0.157& -0.423& -5.183·{10}^{-12}\\ 1.151& 2.214& -3.116& 7.245\\ 0.511& 0.602& -5.533& 12.361\end{array}\right]$ $\mathrm{GEOM}:={\left(A^T·A\right)}^{-1}·A^T$ $L:=\left[\begin{array}{c}\left|\mathrm{D1v}\right|-\mathrm{t1}·c\\ \left|\mathrm{D2v}\right|-\left|\mathrm{D1v}\right|-\left(\mathrm{t2}-\mathrm{t1}\right)·c\\ \mathrm{D1v}·\mathrm{Vt1}-c·\mathrm{D1}·\mathrm{freqdelm1}\\ \mathrm{D2v}·\mathrm{Vt2}-c·\mathrm{D2}·\mathrm{freqdelm2}\end{array}\right]$ $\mathrm{GEOM}=\left[\begin{array}{cccc}-1.771& -1.895& 0.952& -0.558\\ 0.813& 1.026& 0.088& -0.051\\ -0.032& -2.341& 0.211& -0.124\\ 0.019& -1.019& 0.051& 0.051\end{array}\right]$ $\mathrm{DR}:=\mathrm{GEOM}·L$ $L=\left[\begin{array}{c}-10.474\\ 2.723\\ 13.749\\ 44.649\end{array}\right]$ $\mathrm{DR}=\left[\begin{array}{c}1.573\\ -6.807\\ -8.666\\ 1.231·{10}^{-4}\end{array}\right]$

[0126] This program calulates geolocation from four satellites using the time difference of arrival (TDOA) at the satellites and an initial estimate of the location. It is necessary to correct for atmospheric effects in the same manner as GPS. It assumes an “inaccurate” clock signal is transmitted by the beacon, which is corrected via the algorithm technique. It is directly analogous to the GPS approach (Sep. 1, 1997). $\mathrm{freq0}:=240.50·{10}^6$ $\mathrm{lat}:=38.8$ $\mathrm{lon}:=-77.0$ $\mathrm{omegaEarth}:=0.0\mathrm{oE}:=7.29·{10}^{-5}$ $\mathrm{timesec}:=0$ $\mathrm{timemin}:=0$ $\mathrm{timehour}:=0$ $\mathrm{time}:=\mathrm{timehour}·3600+\mathrm{timemin}·60+\mathrm{timesec}$ $\begin{array}{cccc}\mathrm{SatPosVect1}:=\left(\begin{array}{c}0\\ 0\\ 7028\end{array}\right)& \mathrm{SatPosVect2}:=\left(\begin{array}{c}0.0\\ 4969.5465\\ 4969.5465\end{array}\right)& \mathrm{SatPosVect3}:=\left(\begin{array}{c}0\mathrm{.0}\\ -4969\mathrm{.5465}\\ 4969\mathrm{.5465}\end{array}\right)& \mathrm{SatPosVect4}:=\left(\begin{array}{c}4969.5465\\ 0.0\\ 4969.5465\end{array}\right)\\ \mathrm{SatVelVect1}:=\left(\begin{array}{c}0\\ 7.2\\ 0\end{array}\right)& \mathrm{SatVelVect2}:=\left(\begin{array}{c}0\\ 5.0911688\\ -5.0911688\end{array}\right)& \mathrm{SatVelVect3}:=\left(\begin{array}{c}0\\ 5.0911688\\ 5.0911688\end{array}\right)& \mathrm{SatVelVect4}:=\left(\begin{array}{c}5\mathrm{.0911688}\\ 0.0\\ -5\mathrm{.0911688}\end{array}\right)\end{array}$ $c:=2.997930·{10}^5$ $\mathrm{RadiusEarth}:=6378$ $\mathrm{alt}:=0.1$ $\mathrm{xvec}:=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$ $\mathrm{yvec}:=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$ $\mathrm{zvec}:=\left(\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 1\end{array}\right)$ $\mathrm{rho}:=\mathrm{RadiusEarth}+\mathrm{alt}$ $\mathrm{phi}:=\mathrm{lon}·\frac{\pi }{180}+\mathrm{time}·\mathrm{omegaEarth}\mathrm{time}·\mathrm{omegaEarth}=0$ $\mathrm{theta}:=\frac{\left(90.000-\mathrm{lat}\right)}{90.000}·\frac{\pi }{2}\mathrm{phi}=-1.344$ $\mathrm{RotMat}:=\left(\begin{array}{ccc}\mathrm{rho}·\sin \left(\mathrm{theta}\right)·\cos \left(\mathrm{phi}\right)& 0& 0\\ 0& \mathrm{rho}·\sin \left(\mathrm{theta}\right)·\sin \left(\mathrm{phi}\right)& 0\\ 0& 0& \mathrm{rho}·\cos \left(\mathrm{theta}\right)\end{array}\right)$ $X:=\mathrm{RotMat}·\mathrm{xvec}$ $Y:=\mathrm{RotMat}·\mathrm{yvec}$ $Z:=\mathrm{RotMat}·\mathrm{zvec}$ $X=\left(\begin{array}{c}1.118163188·{10}^3\\ 0\\ 0\end{array}\right)$ $Y=\left(\begin{array}{c}0\\ -4.843296871·{10}^3\\ 0\end{array}\right)$ $Z=\left(\begin{array}{c}0\\ 0\\ 3.996541769·{10}^3\end{array}\right)$ $\mathrm{XP}:=\left(\begin{array}{c}{X}_0\\ Y_1\\ Z_2\end{array}\right)$ $\begin{array}{cc}\mathrm{Ve}:=\mathrm{XP}\times \mathrm{zvec}·\left(\mathrm{omegaEarth}\right)& \mathrm{fm1}:=2.40504804885·{10}^8\\ \mathrm{Vtl}:=\mathrm{Ve}+\mathrm{SatVelVect1}& \mathrm{fm2}:=2.40503637992·{10}^8\\ \mathrm{Vt2}:=\mathrm{Ve}+\mathrm{SatVelVect2}& \mathrm{fm3}:=2.405023247783·{10}^8\\ \mathrm{Vt3}:=\mathrm{Ve}+\mathrm{SatVelVect3}& \mathrm{fm4}:=2.405018767603·{10}^8\\ \mathrm{Vt4}:=\mathrm{Ve}+\mathrm{SatVelVect4}& \mathrm{freqdelm1}:=\frac{\mathrm{fm1}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm2}:=\frac{\mathrm{fm2}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm3}:=\frac{\mathrm{fm3}-\mathrm{freq0}}{\mathrm{freq0}}\\ & \mathrm{freqdelm4}:=\frac{\mathrm{fm4}-\mathrm{freq0}}{\mathrm{freq0}}\end{array}$ $\mathrm{D1}:=\sqrt{{\left({\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D2}:=\sqrt{{\left({\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D3}:=\sqrt{{\left({\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D4}:=\sqrt{{\left({\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\right)}^2+{\left({\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\right)}^2+{\left({\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\right)}^2}$ $\mathrm{D1v}:=\left(\begin{array}{c}{\mathrm{SatPosVect1}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect1}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect1}}_2-{\mathrm{XP}}_2\end{array}\right)$ $\mathrm{D2v}:=\left(\begin{array}{c}{\mathrm{SatPosVect2}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect2}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect2}}_2-{\mathrm{XP}}_2\end{array}\right)$ $\mathrm{D3v}:=\left(\begin{array}{c}{\mathrm{SatPosVect3}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect3}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect3}}_2-{\mathrm{XP}}_2\end{array}\right)$ $\mathrm{D4v}:=\left(\begin{array}{c}{\mathrm{SatPosVect4}}_0-{\mathrm{XP}}_0\\ {\mathrm{SatPosVect4}}_1-{\mathrm{XP}}_1\\ {\mathrm{SatPosVect4}}_2-{\mathrm{XP}}_2\end{array}\right)$ $\begin{array}{cccc}\mathrm{t11}:=\frac{\left|\mathrm{D1v}\right|}{c}& \mathrm{t22}:=\frac{\left|\mathrm{D2v}\right|}{c}& \mathrm{t33}:=\frac{\left|\mathrm{D3v}\right|}{c}& \mathrm{t44}:=\frac{\left|\mathrm{D4v}\right|}{c}\\ \mathrm{t1}:=0.019420602& \mathrm{t2}:=0.033103370& \mathrm{t3}:=0.004962104& \mathrm{t4}:=0.020894343\end{array}$ $A:=\left(\begin{array}{cccc}\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& 1·c\\ \frac{{\mathrm{D2v}}_0}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D2v}}_1}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D2v}}_2}{\left|\mathrm{D2v}\right|}-\frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& -1·c\\ \frac{{\mathrm{D3v}}_0}{\left|\mathrm{D3v}\right|}-\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D3v}}_1}{\left|\mathrm{D3v}\right|}-\frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D3v}}_2}{\left|\mathrm{D3v}\right|}-\frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& -1·c\\ \frac{{\mathrm{D4v}}_0}{\left|\mathrm{D4v}\right|}-\frac{{\mathrm{D1v}}_0}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D4v}}_1}{\left|\mathrm{D4v}\right|}-\frac{{\mathrm{D1v}}_1}{\left|\mathrm{D1v}\right|}& \frac{{\mathrm{D4v}}_2}{\left|\mathrm{D4v}\right|}-\frac{{\mathrm{D1v}}_2}{\left|\mathrm{D1v}\right|}& -1·c\end{array}\right)$ $\mathrm{GEOM}:={\left(A^T·A\right)}^{-1}·A^T$ $A=\left[\begin{array}{cccc}-0.192& 0.832& 0.521& 2.998·{10}^5\\ 0.079& 0.157& -0.423& -2.998·{10}^5\\ -0.56& -0.917& 0.133& -2.998·{10}^5\\ 0.807& -0.059& -0.365& -2.998·{10}^5\end{array}\right]$ $L:=\left[\begin{array}{c}\left|\mathrm{D1v}\right|-c·\mathrm{t1}\\ \left|\mathrm{D2v}\right|-\left|\mathrm{D1v}\right|-\left(\mathrm{t2}-\mathrm{t1}\right)·c\\ \left|\mathrm{D3v}\right|-|\mathrm{D1v}|-\left(\mathrm{t3}-\mathrm{t1}\right)·c\\ \left|\mathrm{D4v}\right|-|\mathrm{D1v}|-\left(\mathrm{t4}-\mathrm{t1}\right)·c\end{array}\right]$ $\mathrm{GEOM}=\left[\begin{array}{cccc}0.112& -0.998& -0\mathrm{.15}& 1\mathrm{.26}\\ 0.849& 0.998& -0\mathrm{.15}& 0\\ 1.768& -1.018& 1\mathrm{.337}& 1\mathrm{.449}\\ -2.019·{10}^{-6}& -1.642·{10}^{-6}& -2\mathrm{.003}·{10}^{-6}& -1\mathrm{.708}·{10}^{-6}\end{array}\right]$ $\mathrm{DR}:=\mathrm{GEOM}·L$ $L=\left[\begin{array}{c}-6.858·{10}^{-5}\\ -2.304·{10}^{-5}\\ 1.681·{10}^{-4}\\ -3.593·{10}^{-5}\end{array}\right]$ $\mathrm{DR}=\left[\begin{array}{c}-5.515·{10}^{-5}\\ -1.064·{10}^{-4}\\ 7.497·{10}^{-5}\\ -9.911·{10}^{-11}\end{array}\right]$