C^JFS. /J/O •' A/£*?S /// .<*? °?% NOAA Technical Memorandum NESS 111 EARTH LOCATING IMAGE DATA OF SPIN-STABILIZED GEOSYNCHRONOUS SATELLITES Washington, D.C August 1980 U.S. DEPARTMENT OF COMMERCE / National Oceanic and Atmospheric Administration / National Environmental Satellite Service NOAA TECHNICAL MEMORANDUMS National Environmental Satellite Service Series The National Environmental Satellite Service (NESS) is responsible for ation of NOAA's environmental satellite systems. the establishment and oper- NOAA Technical Memorandums facilitate rapid distribution of material that may be preliminary in nature and so may be published formally elsewhere at a later date. Publications 1 to 20 and 22 to 25 are in the earlier ESSA National Environmental Satellite Center Technical Memorandum (NESCTM) series. The current NOAA Technical Memorandum NESS series includes 21, 26, and subsequent issuances. Publications listed below are available (also in microfiche form) from the National Technical Informa- tion Service, U.S. Department of Commerce, Sills Bldg., 5285 Port Royal Road, Springfield, VA 22161. Prices on request. Order by accession number (given in parentheses). Information on memorandums not listed below can be obtained from Environmental Data and Information Service (D822), 6009 Executive Boulevard, Rockville, MD 20852. M. H. 1966-75. Lushine, Satellite NESS 66 A Summary of the Radiometric Technology Model of the Ocean Surface in the Microwave Region. John C. Alishouse, March 1975, 24 pp. (COM-75-10849/AS) NESS 67 Data Collection System Geostationary Operational Environmental Satellite: Preliminary Report. Merle L. Nelson, March 1975, 48 pp. (COM-75-10679/AS) NESS 68 Atlantic Tropical Cyclone Classifications for 1974. Donald C. Gaby, Donald R. Cochran, James B. Lushine, Samuel C. Pearce, Arthur C. Pike, and Kenneth 0. Poteat, April 1975, 6 pp. (COM-75- 10676/AS) NESS 69 Publications and Final Reports on Contracts and Grants, NESS-1974. April 1975, 7 pp. (COM- 75-10850/AS) NESS 70 Dependence of VTPR Transmittance Profiles and Observed Radiances on Spectral Line Shape Parame- ters. Charles Braun, July 1975, 17 pp. (COM-75-11234/AS) NESS 71 Nimbus-5 Sounder Data Processing System, Part II: Results. W. L. Smith, H. M. Woolf , C. Hayden, and W. C. Shen, July 1975, 102 pp. (COM-75-11334/AS) NESS 72 Radiation Budget Data From the Meteorological Satellites, ITOS 1 and NOAA 1. Donald Flanders and William L. Smith, August 1975, 20 pp. (PB-246-877/AS) NESS 73 Operational Processing of Solar Proton Monitor Data (Revision of NOAA TM NESS 49). Stanley R. Brown, September 1975, 15 pp. (COM-73-11647) NESS 74 Monthly Winter Snowline Variation in the Northern Hemisphere From Satellite Records, Donald R. Wiesnet and Michael Matson, November 1975, 21 pp. (PB-248-437/6ST) NESS 75 Atlantic Tropical and Subtropical Cyclone Classifications for 1975. D. C. Gaby, J. B. B. M. Mayfield, S. C. Pearce, and K.0. Poteat, March 1976, 14 pp. (PB-253-968/AS) NESS 76 The Use of the Radiosonde in Deriving Temperature Soundings From the Nimbus and NOAA Data. Christopher M. Hayden, April 1976, 19 pp. (PB-256-755/AS) NESS 77 Algorithm for Correcting the VHRR Imagery for Geometric Distortions Due to the Earth Curvature, Earth Rotation, and Spacecraft Roll Attitude Errors. Richard Legeckis and John Pritchard, April 1976, 31 pp. (PB-258-027/AS) NESS 78 Satellite Derived Sea-Surface Temperatures From NOAA Spacecraft. Robert L. Brower, Hilda S. Gohrband, William G. Pichel, T. L. Signore, and Charles C. Walton, June 1976, 74 pp. (PB-258-026/AS) NESS 79 Publications and Final Reports on Satellite Service, June 1976, 10 pp. NESS 80 Satellite Images of Lake Erie Ice: June 1976, 15 pp. (PB-258-458/AS) NESS 81 Estimation of Daily Precipitation Over China and the USSR Using Satellite Imagery. Follansbee, September 1976, 30 pp. (PB-261-970/AS) NESS 82 The GOES Data Collection System Platform Address Code. Wilfred E. Mazur, Jr., October 1976, 26 pp. (PB-261-968/AS) NESS 83 River Basin Snow Mapping at the National Environmental Satellite Service. Stanley R. Schneider, Donald R. Wiesnet, and Michael C. McMillan, November, 1976, 19 pp. (PB-263-816/AS) NESS 84 Winter Snow-Cover Maps of North America and Eurasia From Satellite Records, 1966-1976. Michael Matson, March 1977, 28 pp. (PB-267-393/AS) NESS 85 A Relationship Between Weakening of Tropical Cyclone Cloud Patterns and Lessening of Wind Speed. James B. Lushine, March 1977, 12 pp. (PB-267-392/AS) NESS 86 A Scheme for Estimating Convective Rainfall From Satelliue Imagery. Roderick A. Scofield and Vincent J. Oliver, April 1977, 47 pp. (PB-270-762/AS) Contracts and Grants, NESS-1975. National Environmental (PB-258-450/AS) January-March 1975. Michael C. McMillan and David Forsyth, Walton A. (Continued on inside back cover) NOAA Technical Memorandum NESS 111 EARTH LOCATING IMAGE DATA OF SPIN-STABILIZED GEOSYNCHRONOUS SATELLITES Larry N. Hambrick Office of Systems Integration Dennis R. Phillips Scientific Programming and Applied Mathematics, Inc Contracts: DOC 7-35229, DOC 01-8-M01-5347 , NA79-KAC00026, and NA79-SAC00744 Washington, D.C. August 1980 -"gg**- UNITED STATES DEPARTMENT OF COMMERCE Philip M. Klutznick, Secretary NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION Richard A. Frank, Administrator National Environmental Satellite Service David S Johnson, Director ■\ ^^xr^ PREFACE The National Environmental Satellite Service (NESS) operates a central facility for determining accurate earth-location in- formation for the images derived from its Geostationary Opera- tional Environmental Satellite (GOES) . The facility is called the VISSR Image Registration and Gridding System (VIRGS) . The main purpose of this document is to set forth a mathe- matical formulation of the entire earth-location process for a spin-stabilized geosynchronous satellite. The implementation of this process on the VIRGS is described. Preceding the main topics is a discussion of the relevant aspects of the GOES sys- tem. An error analysis and some actual results are presented in the final section. The original objective of the VIRGS was that it provide a stand-alone capability to operationally produce highly accurate earth-location parameters. The "stand-alone" feature is impor- tant because it eliminates the expense of operating ranging stations and performing separate data processing tasks. "High accuracy" is as important; without it the users of GOES imagery are faced with the undesirable alternatives of either using poorly located data or determining accurate earth-location parameters on their own. The objectives of the VIRGS have not yet been fully real- ized. This is due in large to the lag-time between the devel- opment of a new technique and its full operational implementa- tion. A technique that uses the positions of stars in the VISSR image to determine some of the earth-location parameters has not yet been implemented. The star position technique streamlines the earth-location process and provides the needed accuracy and stand-alone features. Also remaining to be fully implemented is the closed-loop capability of the VIRGS--a capa- bility for complete quality control of the earth-location process. The authors' intent is to show that the objectives of the VIRGS can be achieved; and it is their hope that the demonstra- tion given herein will aid in that achievement. Some acknowledgements are deserved. Hank Schmidt began casually reading an early version of the document and ended up making significant improvements to the presentation. Jim Ellickson and Matt Jurotich reviewed an early version. Ron Gird is thanked for his encouragement and energetic participa- tion in the VIRGS implementation. Vince Oliver and others, inside and outside of NESS, who as users actively supported improvements to the operational earth-location of VISSR images, were instrumental in bringing about VIRGS. The McIDAS team at the University of Wisconsin (SSEC) planted the seed and then provided their know-how in delivering the basic VIRGS. Eric iii Smith, J.T. Young and others were responsible for the develop- ment and demonstration of the early navigation software on McIDAS. VIRGS is operated and maintained daily, in accordance with procedures and guidelines passed through management, by a group consisting of satellite meteorologists (oddly enough) , software analysts, computer operators, and electronic technicians. Members of this group are recognized for their daily efforts. IV Contents Preface i Abstract 1 I. Introduction 1 II. General review 3 A. The satellite 3 B. The VISSR 4 C. Functions of the SDB 5 Do The earth-location capability 7 III. Mathematical algorithms and models for the earth location process <,.... 11 A. Preview 11 B. Geometric transformations between image and earth coordinates 11 C. Finding the misalignment parameters and attitude. 21 D. Orbit determination algorithms 25 IV. Description of implementation 33 A. Orbit/attitude determination software 33 B. Standard geometry routine for users 35 C. VIRGS (the VISSR Image Registration and Gridding System) 36 D. Rationale for centralized service 40 E. Error budget and actual results 41 V. Summary 48 References 49 Appendices; A. Converting a pointing vector to line and element numbers B. Determination of attitude and orbit plane orientation C e Accuracy of approximate orbit equation Do Gaussian elimination E. Fortran listings of orbit and attitude determination routines; FINDAT and UPGORB F Fortran listings of standard geometry transformation routines; EARLOC and IMGLOC VI EARTH LOCATING IMAGE DATA OF SPIN-STABILIZED GEOSYNCHRONOUS SATELLITES Larry N. Hambrick Office of Systems Integration Dennis R. Phillips Scientific Programming and Applied Mathematics, Inc. ABSTRACT. Aspects of the GOES system rele- vant to earth-location are reviewed. The essential features of an earth-location capa- bility that satisfies the needs of users of the data are specified. A mathematical for- mulation of the geometry and orbit/attitude determinations, including the actual algo- rithms, is presented. The implementation of these methods on VIRGS is described. Final- ly, the rationale for a centralized, once- for-all users earth-location service and an error budget with some actual results are presented. I. INTRODUCTION The Geostationary Operational Environmental Satellites (GOES) are examples of spin-stabilized geosynchronous satellites. The GOES provide continuous viewing of the weather features over the Americas and their adjacent waters. The images, or views, are obtained every half-hour, and sometimes more frequently, from each of two satellites. The imaging instrument of the GOES is called a Visible and Infra-red Spin-Scan Radiometer (VISSR) . About a dozen users receive VISSR data directly from the GOES and over a hundred others are served by a Central Data Distribution Facility. The GOES program began in 1974 and is expected to remain in operation throughout the 1980's. During the next ten years the number of uses and users of GOES data is expected to increase. In order for the VISSR imagery to be useful, it must be earth-located. "Earth-location" means different things to different users. It can mean (1) the ability to determine 1 precisely the earth-location (latitude and longitude) of any feature in an image; (2) the presence of fixed surface features of earth in time-lapsed (animated) sequences of images; or (3) having geographic and political boundaries accurately laid over the images. Today many users of VISSR images use manual techniques to shift and align images. These methods are cumbersome, time consuming, and inaccurate. They provide, at best, partial solutions to the earth-location problem. A general solution --one that enables a user to map any picture element (VISSR f ield-of-view) into precise earth coordinates--must be based on a complete mathematical formulation of the problem. There are two reasons for this: (1) the basic geometry that describes a mapping or transformation from satellite/instrument (GOES/ VISSR) coordinates to earth coordinates is complex; and (2) the dynamics of the satellite's orbit and attitude must be modelled With the information needed to perform accurate earth- location inserted into the VISSR data stream, all direct VISSR users would benefit by not having to determine it on their own. II. GENERAL REVIEW A. The Satellite A satellite in an ideal geosynchronous orbit remains motion- less with respect to an earth-based observer. For this to occur, the satellite must be in the equatorial plane and the centrifugal force of its orbital motion must, at all times, equal the gravitational force between it and the earth. How- ever, there are perturbating forces that make it highly imprac- tical to maintain such an ideal orbit. The GOES orbit is not perfectly circular and it is slightly inclined to the equatorial plane. Furthermore, the perturba- tive forces, caused by anomolies in the earth's gravitational force field, the gravitational affect of the moon and the sun, and the pressure of the solar wind continually vary the orbit. Departure from an ideal geosynchronous orbit causes the satel- lite subpoint to trace out a "distorted figure eight" during its one-day orbital period. The VISSR is rigidly mounted in an upright position to an axis about which the satellite spins at 100 rpm. Ideally the attitude of the spin axis would be normal to the orbit plane. An ideal orbit and spin axis attitude would produce the same "perspective" in all VISSR images. In practice external torques from solar radiation and impulses from VISSR scan stepping cause the satellite spin axis to wobble around a vector normal to the orbit plane in a manner governed by rigid body dynamics. The practical manifestation is a nearly linear precession of the spin axis over a period of several days accompanied by a small nutation with a cycle period of several seconds. A set of orbit parameters determines the satellite's position and a set of attitude parameters determines the orientation of the satellite's spin axis relative to the earth, as functions of time. The satellite's position and orientation are control- led so as to stay within certain bounds and hence avoid exces- sive apparent earth motion in the VISSR images. The method of control is to fire the on-board rockets in a manner that appro- 3 priately alters the orbit and attitude. These are called "maneuvers" and such maneuvers are vital in maintaining the quality of the GOES service. When these maneuvers occur, sig- nificant discontinuities in the orbit and attitude result. The maneuvers are usually scheduled several days in advance and can be made as frequently as once a week. The orbit and attitude discontinuities caused by satellite manuevers interrupt the process of predicting the orbit and attitude state with a mathematical model. Thus the earth location capability is also interrupted. A navigational system, such as VIRGS, requires accurate orbit and attitude states within several hours fol- lowing a maneuver. B. The VISSR The imaging instrument on GOES is the VISSR. The VISSR is rigidly mounted to the satellite's body. An image is formed by the VISSR Field of View (FOV) scanning the earth, west to east, as the satellite spins. After each spin, a scan mirror is moved (stepped down one notch) so that the next scan sweep is slighly south of the previous sweep. Electronic sampling generates equal angle (or time) spacing between successive image elements on each scan line. Throughout each scan line the scan mirror is at a fixed position. The scan mirror stepping increment is 192 ^radians and the sampling interval is 84/\radians for the infrared channel(s) and 21^ radians for the visible channel. There is a single infra-red detector. It has a 192 by 84 //cradian field of view and provides 7 by 3 km subpoint resolution. There are 8 vertically aligned visible detectors which are sampled simulataneously giving a 21 / 4iradian square field of view, for 0.8 km resolution at the subpoint. The normal operation of the VISSR is to acquire a full earth disc image frame (1821 successive scans) once every 30 minutes. Sometimes the VISSR is commanded to acquire less than a full earth disc image by scanning over any interval of scan positions within the range of 1 to 1821. When this is done, images can be obtained for example during 3, 7, and 15 minute 4 intervals, depending on the size of the images. The earth- location capability described herein is compatible with these various imaginq modes. The mechanical alignment of the VISSR to the satellite spin axis and the sun pulse detector (see next paragraph) is imper- fect. This creates significant values for three misalignment angles that are analogous to the roll, yaw, and pitch angles of a spacecraft. The VISSR' s mechanical frame could be tilted forward (pitch) , tilted to the side (yaw) as well as being rotated around the vertical (roll). The determination of these misalignment angles, which are constant over periods of weeks, is a vital step in achieving accurate earth location. A sun sensor and an earth-edge detector are used to initiate the sampling along each scan line. During each spin of the satellite the delay between sensing the sun and detecting the earth's west edge is measured. The delay that is measured on one spin n is used to initiate sampling on spin n+1. This process merely ensures that the earth view is contained in the data acquisition interval. The task of precisely referencing the west border of the image frame on successive scan lines is performed on the ground in real-time by the Synchronous Data Buffer. C. Functions of the Synchronous Data Buffer (SDB) The VISSR data are transmitted in real-time at 28 Mbps* to a Command and Data Acquisition (CDA) facility on the ground. At the CDA an SDB is used to "stretch" the data stream for re- transmission through GOES at 1.74 Mbps. The stretching opera- tion uses the 340° portion of each satellite revolution, when the VISSR is not pointing at the earth, to rebroadcast the VISSR data at 1/16 the raw data rate. Users of the data are *Megabits per second (Mbps) therefore interested in the form of the SDB's VISSR output rather than the satellite's raw output. The SDB accomplishes other critical tasks in the process of stretching the data stream. From the point of view of earth- locating the data, one of the SDB's most critical functions is to center each VISSR scan line on the earth. The SDB, while maintaining synchronization with the satellite's spin rate, precisely times the start (west border) of each stretched scan line. To do this the time of the sun pulse and a "beta" angle are used. The "beta" angle is the angle subtended at the satellite by lines projected from the sun and earth to the satellite, in the satellite's spin plane. The beta angle is computed by a sun position model using the satellite's orbit and attitude and the VISSR misalignment parameters. This computation is readily available (as a by-product) from an earth-location capability such as VIRGS. The value of beta used on each scan line references each sample of that scan line with respect to the sun's position. The well-defined relation- ship between beta and every sample of the scan line provides an absolute reference. Hence, even if the "beta" provided to the SDB is inaccurate, the user of the data can still reference the imagery data to the sun's position. The "beta" used by the SDB is documented in the stretched VISSR data stream (SDB output) . The image frame as seen by the user consists of up to 1821 scan lines, each with 3822 infrared elements and 15288 visible elements across the line. Each infrared element is referenced to a coincident array of 8 by 4 visible picture elements (pixels) . The referencing along each scan line is accomplished by the SDB. The referencing across scan lines is governed by the VISSR step scan characteristics. For each scan line the actual VISSR scan line number and the value of "Beta" are contained in a "documentation block" that is prefixed to each line of stretched VISSR data. The documen- tation block also contains values of the parameters that speci- fy the orbit, the attitude, and the misalignment angles. This 6 block is described in the appendix of reference (6). These parameters are computed at the central facility (VIRGS) and transferred to the SDB. They are intended to furnish the users of VISSR data with the information needed to perform predictive earth-location. Clearly any significant errors in these para- meters force the user to either work with mislocated 'image data, or alternatively, to do the navigation independently. If the orbit-attitude documentation block contains precise, or the best achievable, values on a highly reliable basis, the user needs only a standardized geometry routine to which the docu- mented parameters are passed as inputs. Of interest to some types of users is the NESS standard full disc grids. They are embedded by the SDB as the 9th bit of each 8 bit infrared (IR) sample. The locations of the NESS standard grid points are generated by VIRGS and transferred to the SDB. For accuracies within the resolution of an IR pixel these grids can be used. Thus a user can avoid the large com- putational task of applying the geometry transform to some 30,000 grid points that outline political and geographic bound- aries. For users whose applications require full visible res- olution grid point location, the documented orbit and attitude parameters should be used. Figure 1 illustrates the overall data flow. D. The Earth-Location Capability The foundation of a VISSR earth-location capability is the geometrical transformation and its inverse. Given the earth- location parameters (the satellite's position and attitude and the VISSR 1 s misalignment angles) the transformation maps the J sample on the I scan line to the latitude and longi- tude of the point on the earth to wnich that sample corresponds The most convenient and accurate method for determining the values of the earth-location parameters uses earth landmarks and stars that are discernable in the VISSR image. The image coordinates (I, J) of known stars, over an interval of time, specify the values of the satellite's attitude and the VISSR' s misalignment angles. These values and the coordinates (I, J) GOES SYNCHRONOUS DATA BUFFER (SDB) STRETCHED 1.74MBPS 1 DIRECT READOUT USERS (WITH STANDARD GEOMETRY ROUTINES 1,'BETA' 2. ORBlT/ATTiTUDE PARAMETERS 3. STANDARD GRID CENTRAL EARTH-LOCATIOM FACILITY (VI RGS*) I NESS DBRECT READOUT * VISSR IMAGE REGISTRATION AND GRIDDING SYSTEM Figure 1. Overall data flow for earth locating VISSR dat of one or more recognizable landmarks (measured during a time interval close to the the time interval of the star measure- ments) provide a solution for the satellite's position as a function of time, i.e., for its orbit. Using star measurements in the navigational process separates the satellite's attitude and misalignment determination from the orbit determination. This disjointness greatly simplifies, and enhances the efficiency of, the orbit/attitude determination process. These advantages are maintained even when more sophisticated mathematical models are employed. This method also simplifies the structure of the software and streamlines the process so that the cpu time required for the orbit/attitude computation is negligible. Finally, the use of stars eliminates the need for the expensive operation of ground stations that provide ranging data which are currently being used to determine the orientation of the orbital plane. Implicit in the technique for determining attitude from star positions and orbit from landmark positions, are mathematical models of the dynamic behavior (or motion) of the spin axis and position. Such models propagate the motion from initial (epoch) values for the parameters. For GOES both the orbit and attitude models are relatively simple owing to the well-behaved motion of the satellite and its spin axis. Having relatively simple models suggests that if the orbit/attitude determination methods were cleverly designed, they could be implemented as a small, efficient software package. The less sophisticated models, i.e., those that account for fewer effects, can only be used to accurately predict satellite motion over short periods. Real-time data users need accurate parameters only a few minutes in advance, say in the prior image frame. The operators of the central facility, to avoid a continual real-time activity of supporting the SDB , would pre- fer a once-per-day cycle. This requires at least a 24 hour predictive capability. Being able to predict over several days or a week has some potential advantages in procedural efficien- cy. It would support such things as maneuver planning and 9 batch processing of the grid point locations during off hours. In brief, the prediction capability should be a minimum of one day, and preferably several days. Since maneuvers of the satellite occur frequently, and inevi- tably cause some degree of discontinuity in the earth-location process, an extremely important aspect of the earth-location capability is quick recovery of full accuracy following a ma- neuver. Predictions of the post-maneuver orbit and attitude are done but cannot be expected to be highly accurate. The modeling of maneuvers is not precise and many maneuvers are not conducted as planned. Using the "stars and landmarks" naviga- tion technique immediately following a maneuver to make succes- sively improved estimates of the orbit and attitude allows re- covery of full accuracy within a few hours. The amount of time required depends on the time-of-day that the maneuvers were performed, relative to the time when the stars and landmarks are visible. Orbit/attitude determination methods should en- able the best possible estimates of parameters with a minimal set of observations. Finally, the central facility must generate and transfer the necessary parameters to the SDB. These parameters, when their values are precise, allow the SDB to center the scan lines; and by using these parameters, users can perform accurate geometric transforms with standard geometry routines. 10 III. MATHEMATICAL ALGORITHMS AND MODELS FOR THE EARTH LOCATION PROCESS A. Preview The algorithms for navigating geosynchronous satellite images will now be described. These algorithms use the positions of recognizable stars and earth-based landmarks measured in VISSR image frames to determine the attitude, the misalignment angles and the orbit parameters. The set of algorithms is divided into three parts: (1) the transformation of image coordinates to earth coordinates, and the corresponding inverse transfor- mation, using values for the attitude, misalignment, and orbit parameters, (2) the determination of attitude and misalignment parameters from measured correspondences for stars between inertial pointing vectors and image frame positions, and (3) the determination of orbit parameters from the measured image positions of identifiable earth-based landmarks and attitude and misalignment parameters. B. Geometric Transformations Between Image and Earth Coordinates The approach for transforming image coordinates to earth coordinates is now outlined. The details of the mathematical formulation of the transformation follow shortly. The first step is to use the attitude and misalignment parameters to find the inertial pointing direction of the spin scan camera (e.g., the VISSR). The inertial pointing must be found as a function of line number, element number and image start time. Next, the position of the satellite as a function of time is found. For given line and element numbers the satellites orbital position and the spin scan camera's inertial pointing direction are known, and it is a straight forward procedure to determine the point on the earth's surface being observed by the satellite's camera. This is done by projecting a ray in the inertial pointing direction from the satellite position and determining the point where it would intersect an oblate spheriod (fig. 2) . 11 $ INERTIAL POIIMTING DIRECTION SATELLITE'S POSITION Figure 2. The inverse problem, i.e., determining the line and element number of the sample which observes a particular point on the earth's surface in a given image, is more complicated. The complications arise from not knowing the time that the point in question is being observed. Knowing the time is essential since the inertial pointing direction of the vector from the satellite to a particular earth location varies as a function of time due to orbital motion of the satellite and the rotation of the earth with respect to inertial space. Since this time is not known, a guess must be made. Using this "first guess" one normalizes the inertial vector from the satellite to the point on the earth's surface into an inertial pointing vector and then computes the line and element number at which the pointing direction of the spin scan camera was coincident with this pointing vector. In general, the time at which this par- ticular sample was taken is different from the first guess. The new time is taken and this procedure is repeated iterative- ly until two consecutive times agree to within the time between two successive scan line sweeps. The discussion shows that three basic computations must be made in order to transform image coordinates to earth coordi- nates and vice versa. These are: (1) the determination of the inertial pointing direction of the spin scan camera as a func- tion of line number, element number, picture time, beta count, and misalignment and attitude parameters, (2) the determination of the intersection between a ray and the surface of an oblate sphere, and (3) the determination of the line and element num- 12 ber at which the axis of view of the spin scan camera is paral- lel to a specific inertial pointing vector for a given set of beta counts and misalignment and attitude parameters. The reader should be aware that "spin scan camera" is substituted for "VISSR" to connote the generality of the formulation.. Also the "pointing direction" of a spin scan camera is the axis of the optical field of view of the camera or, for brevity, the "axis of view." Figure 3 illustrates the elements of the basic geometry. 1. Determination of Pointing Direction of the Spin Scan Camera The pointing direction of the spin scan camera is found as a function of time, beta counts, and line and element numbers. To determine this one must have an accurate attitude state, ice., the attitude itself and the precession of the attitude. One must also know the parameters that determine the state of misalignment between the satellite's body axis and its actual spin axis. Those parameters establish the stepping path of the spin scan camera as a function of line number. Nominally the spin axis and the body axis of the satellite are expected to coincide. However, in general, some small misalignments are observed. The misalignment causing a bias in the line positions where identifiable stars and earth-based landmarks appear in the image is designated as pitch and param- eterized with the angle C . The misalignment that causes a bias in the element direction is designated as roll and parame- terized with the angle p and, finally, the misalignment causing skew in the element direction as a function of line number is designated as yaw and parameterized with the angle r| . The pointing position at the start of each scan line is ref- erenced to the sun pulse detection for that scan line by the VISSR to sun pulse detector sweep angle, T , the current beta count, p , and whatever element bias and skew that are present. It is natural to study the position of the spin scan camera as a function of scan line number at those instances when the sun pulse is detected. Two satellite centered coordinate systems are created for this study. The first coordinate system has 13 >• I N c o •H 4J sin( i) , -sin( i) « cos(ft) , cos(i) ) (4.2) where fi is the ascending node and i is the inclination of the orbit. Since the inclination is always taken to be positive, the ascending node and the inclination can be determined as follows (in Fortran) : Q = ATAN2(a, -b) and (4.3) i = ASIN(SQRT(a 2 +b 2 ) ) . This completes the discussion of determining the orbit plane for the satellite. In order to study the angular motion of the satellite within its orbit plane, a planar coordinate system is set up in the orbital plane. The x-axis points from the center of the earth to the ascending node and the y-axis, also lying in the orbital plane, is 90° beyond the x-axis in the direction of satellite motion. Refer to figure 12. Angles in this planar coordinate system are measured as going from the x-axis towards the y-axis, i.e., in the direction of satellite motion. Specify angular positions occuring at time t. as a . . The argument of perigee is specified by oj,time is specified by t, and the time when the satellite is at perigee is specified by t . Within this framework the motion XT of a satellite in a Keplerian orbit is well represented by the following equation: a -oj-2 e . sin(ff) . cos(uj) + 2e . cos {a) . sin (w) =>\-it-t ), (4.4) where n is the mean motion constant. This representation is -3 accurate to within .12 km for eccentricities less than 10 Since only the form of this equation is of primary interest, the derivation verifying the accuracy is presented in Appendix C. Normally eq. (4.4) is written in the following equivalent form 28 EARTH Figure 12 t. - c, + c~*a. + c,sin( a.) + c.cos(a.) (4.5) 1 1 2 l 3 l 4 i where the index i indicates the times and angular positions corresponding to the landmarks and c,=tp-u/n, c~ = \/x\, c 3 =-2e*cos(co) /n, and c. = 2 £• sin(oj) /n. It is vital that the constants |c| be related directly to Keplerian orbit param- eters so that the reason for determining these constants becomes clear. Observe the following relations: n = 1/c 2 ' (4.6) where n is the mean motion constant; .2/3/ 2/3 .. -, a = r* k /n ' , (4.7) where a is the semimajor axis, r is the radius of the earth, and k is the gravitational constant of the earth; 2 2 Vs e = (c 3 + c 4 ) / V(2c 2 ) and (4.8) w = ATAN2(c 4 ,-c 3 ) , where e and u> are the eccentricity and argument of perigee, respectively. The inverse relationship for c 2 ,c, and c A are given later in (4.13). The inverse relationship between the mean motion constant n 29 and c~ is due to n being the multiplicative constant relating time to angular position. The relationship for determining the semimajor axis comes from Escabol's Methods of Orbit Determina - tion, henceforth referred to as EscaboL The argument of perigee is determined in the above manner because (c., -c-J points parallel to the vector (sin(w) , cos(w) ) due to the eccentricity and mean motion constant both being positive. The formula for eccentricity comes from the observation that (c 3 2 +c 4 2 )/c 2 =(-2€Cos(co) ) 2 + (2esin(w)) 2 = 4e 2 . (4.9) The computation for finding a mean anamoly for a given epoch is more involved. First one finds the time of perifocal passage by substituting w in (4.5) which is c, + c ? oo + c^sin(w) + c 4 cos(oj)=t . Since the vector (c 3 ,c 4 ) is parallel to (-cos(w) ,sin(w)) and hence perpendicular to (sin(oo) , cos(w)) , the last two terms cancel, giving t ^t+Cow. P J. £• Now, by definition, the mean anamoly at an epoch time t equals n»(t -t ) . Hence, adopting the programming symbol e p MEANOM for the mean anamoly, MEANOM = n(t - C, - c 2 «w) = (t - C^) /c 2 -co. (4.10) This completes the process of transforming the (cf into orbital elements. Now consider the problem of finding the jcj- from a set of measurements (a-, t .) . For this purpose linear regression is applied. The following expression is minimized: S= 2 ( c i +c 2 cy i + c 3 sin(a i ) + c 4 cos(o i )-t.) . (4.11) i=l To solve for the |c[, set the following four partial derivatives to zero: "8c, (c l' C 2' c 3 f C 4 } = °' 9c 2 (c l' c 2' c 3' c 4 } = °' — (c lf c 2 , c 3 , c 4 ) = and 30 9S . \ _ n 8c" 4 (C l' c 2' C 3' C 4 ) U ' This results in a matrix equation for the c.'s which is now presented with the following abuse of notation: SM stands for sum over i for 1=1,.=, N, A for the angles a., SN for sin(tt.) , CS for cos(ff-) , T for the time t., and 2 means ii i that the immediately proceeding term is squared. N SMA SMSN SMCS \ / c -\\ / SMT \ SMA SMA2 SMASN SMACS \ / C2 \_ / SMAT \ SMSN SMASN SMSN2 SMSNCS I* I c 3 / I SMSNT I SMCS SMACS SMSNCS SMCS2 / \c^/ \SMCST/ The details of the specific mathematical technique used to in- vert this matrix are given in Appendix D. Notice here that the inverse of the matrix is ill-conditioned over short time inter- vals. This results from the inability to resolve the differ- ence between a straight line and part of a sine wave over a short arc. For this reason, to achieve results comparable in accuracy to one's measurements, one should make landmark meas- urements distributed over a time interval at least 18 hours long. A convenient way to circumvent making such a long stretch of landmark measurements is to make 8-hour stretches of landmark measurements 24 hours apart. In daily operations it is often the case that one uses the semimajor axis computed from yesterday's and today's data but only uses today's data to compute the three along-track orbit parameters: mean anamoly? eccentricity, and argument of peri- gee. This is a two-step procedure. First compute the constant c~ by inverting the method used to compute the semi-major axis from the constant c~. Then compute the remaining c.'s by applying linear regression in the same manner used to find all four c.'s. In this case the resultant three-by-three matrix is inverted by applying Cramer's Rule. The along-track orbit parameters are then found in the same manner as before. This particular approach is potentially valuable for the time period immediately following orbit maneuvers if one can obtain 31 an estimate of the orbit's semimajor axis from energy consider- ations. The estimate for the orbit's mean anamoly can be refined when there are good estimates of the other orbit parameters by requesting that only one further parameter be generated. In this case c~, c^ and c. are generated by the following formulae: c 2 = k.(r/a) 3/2 c^ = -2 • e •cacos(w) and (4.13) c 4 = 2« e »cisin (w) . Finally, if only two along-track orbital elements are wanted, linear regression is used to solve for the parameters c, and c 2 r and these parameters are transformed into a semimajor axis and mean anamoly while at the same time making the eccen- tricity and argument of perigee equal to zero. There are also ways of computing c, and c 2 while retaining nonzero values for the eccentricity and argument of perigee. These are cur- rently not implemented. This ends the discussion of orbit determination. 32 IV. DESCRIPTION OF IMPLEMENTATION A. The Orbit and Attitude Determination Software The functions of the software implementation of the star nav- igation can be split into three categories: (1) the scheduled ingest of images containing either stars or earth-based land- marks and the measurements of the image locations of these fea- tures, (2) the maintenance of navigational files containing the measurements and parameters necessary for navigation and (3) the execution of software that determines the necessary orbit, attitude and misalignment parameters. Further details of the software execution for tasks (1) and (2) are omitted because they are not germane to the purpose of this paper. Instead the mathematical aspects of measurement errors in task (2) affec- ting the scheduling in task (1) are discussed. These aspects affect the scheduling of stars for the determination of the attitude state and misalignment parameters and the scheduling of earth based landmarks for the determination of orbital state. Star measurements should be distributed with respect to time and the right ascension of the stars 1 coordinates to provide an adequate data base for generating the attitude parameters, the precession parameters and the misalignment parameters. Simi- larly landmark measurements should be distributed appropriately with respect to time to provide an adequate data base for gen- erating orbit parameters. Such distributions of measurements are essential to prevent an ill-conditioned inverse of a matrix from blowing up the impact of measurement errors. It is impor- tant to know what distribution of particular measurements are necessary for the accurate determination of specific subsets of these parameters. Such knowledge enables one to establish a schedule for ingesting stars and earth based landmarks. For attitude determination alone, assuming the precession rate is small, at least two star measurements are necessary with the right ascension separation, mod (180°), being 15° or more. Mathematically the ideal case is to have the right as- censions separated by 90°. For attitude and pitch alignment determination, three star measurements are necessary, two with 33 right ascension separation, mod (180°) , of 15° or more and a different pair (i.e., a member of the first pair matched up with the third star measurement) with right ascension measure- ment separation, mod (360°) , of 15° or more. For determining the attitude and precession parameters, four star measurements consisting of two disjoint pairs with each pair of stars having right ascension separation, mod(180°) , of 15° or more and with time separation from the first pair to the second pair of at least 18 hours. In this case the same stars can occur in each pair. For determining the roll and yaw misalignment parame- ters, two star measurements are needed with line separation of 5,000 to 10,000 lines. Questions have been raised about the viability of routinely distinguishing stars from background noise. Such doubts need to be examined because of the central role star measurements play in a stand-alone navigational system. Since star measure- ments are not currently being used to determine the satellite's attitude, an auxiliary system of ground stations and a second computer has to be maintained. It is vitally important to consider each star measurement as a single point in a consistent data set. Regarding each star measurement as a separate entity makes the measurement task more difficult and less fruitful. Simply stated, at a single time the spin axis vector points in only one direction. Hence, if a set of star measurements yields large residuals after a spin axis determination, at least one measurement is wrong. The measurement ( s) whose residuals are noticely distinct from the residuals of other measuremenets should be deleted. This idea can also be used predictively since the change of attitude over a one day period is small. Hence the measured position of a star should be close to its predicted position. Use of this idea in a procedual manner should greatly reduce problems of distinguishing stars from noise. In addition, a star often occurs more than once on a given day. By adding the image frames containing a recurrent star, with an appropriate element lag displacement, the signal to 34 noise ratio of the star's occurrence will be significantly enhanced. Also, the full precision (6 bit) visible image should be used in scanning for stars. For the determination of the orbit's inclination and ascend- ing node, at least two landmarks are necessary with time sepa- ration of 4 to 6 hours. For the determination of the orbit's mean anamoly once the other five orbit parameters are known, only one landmark measurement is needed. For the determination of the orbit's mean anamoly, eccentricity, and argument of per- igee once the other three orbital parameters are known, at least three landmark measurements distributed over a six to eight hour interval are necessary. In order to determine all the orbital parameters at once, at least four landmark measure- ments distributed over a sixteen to eighteen hour interval are necessary . The computer memory space necessary for both the attitude and orbit determination programs is less than 16,000 words and their execution times are below 1 second. Both programs are easily transferable to other computers with Fortran compilers. To achieve such a transfer one has to provide files for a navi- gational data base, routines to monitor and maintain this data base and interface programs between the navigational data base and the orbit and attitude determination software. Fortran listings of the orbit and attitude determination rou- tines, UPGORB and FINDAT, are provided in Appendix E. B. Standard Geometry Routine for Users Direct readout users have access to the orbit, attitude, and misalignment parameters. These parameters as determined by VIRGS are updated in the stretched VISSR orbit/attitude block. A description of the format and definitions of this block is given in reference (6). A geometric transformation that uses these parameters as in- put is available from NESS in a standard form as two Fortran routines. "IMGLOC" transforms from earth latitude and longi- tude to image frame line and element. "EARLOC" is the inverse routine, transforming from VISSR line and element to earth lat- 35 itude and longitude. Listing of these two routines are pro- vided in Appendix F. These routines do not use "Beta", the angle subtended at the satellite by the sun and earth, even though it is provided in the documentation block. The SDB uses a "Beta" function that is precisely derived from the documented orbit, attitude, and misalignment parameters. Therefore, centering of scan lines on the earth disc has to be assumed by the user. This simplifica- tion's benefit for the user is that fewer parameters have to be stored. Modification of the routines to account for earth disc offsets is a small, straightforward task. An assembly language version of "IMGLOC" was implemented on the control computer of the Central Data Distribution Facility operated by NOAA/NESS for the GOES-TAP service. Using the earth-location parameters in the VISSR data stream, the routine locates the image sectors in near real-time. The routine has been operating since June, 1979. C. The VISSR Image Registration and Gridding System (VIRGS) VIRGS is the central earth location facility. VIRGS, in terms of hardware and basic software, is a copy of the McIDAS II (Man-computer Interative Data Access System) . McIDAS II was designed and built by the Space Sciences and Engineering Center (SSEC) of the University of Wisconsin as an upgraded version of the original McIDAS (3,4). The VIRGS was installed by SSEC at the World Weather Building in Camp Springs, Md in July 1978. The operational use of VIRGS began in May, 1979. The orbit/attitude determination techniques and software package described in this paper were developed by SPAAM, Incor- porated under contract to NESS during the past 2 years as a follow-on to work began by Dennis Phillips while at SSEC. The VIRGS (or Mcldas II) consists of a small computer (the Harris/6 with 64,000 words of memory), an on-line 40 Mbyte disc, and an interactive display console for the analyst. The important characteristic of Mcldas is the elaborate set of "key-ins" the analyst has available for processing and dis- play. The demonstrations of accurate VISSR image location were 36 first performed on the Mcldas in 1974 using the interactive key-ins to iteratively adjust orbit-attitude parameters to achieve accurate locations of visible landmarks in the images (2,5) . The VIRGS performs the central earth location process on a daily basis for each GOES in four basic steps: i. measurement of star and landmark positions in the vis- ible imagery (and sometimes in the infra-red imagery) . ii. determination of the misalignment angles and attitude from the star observations and the orbit from the landmark observations iii. production and transfer (to the SDB) of the "Beta" function, earth-location parameters for the documenta- tion block, and the locations of points in the NESS standard grid table, all for the upcoming 24-hour period, iv. monitoring of accuracy of the earth-location through- out the day, with updating if necessary. The primary roles of the human operators-analysts in performing this daily operation are: 1. Scheduling the ingest of image sectors containing stars and then the measurement, using a joystick con- trolled cursor, of the bright spot corresponding to each star. 2. Selecting cloud-free areas on earth that contain fa- vored landmarks, scheduling the ingest of the appro- priate sectors, and then measuring with the cursor the image positions of the recognizeable features. 3. Judging a sufficient set of star and landmark observa- tions and entering, in proper sequence, the several key-ins that determine misalignment parameters, atti- tude, and orbit. 4. Verifying that error residuals as output by various routines are suitably small and that the various files have been updated. 5. Entering key-ins that execute routines that produce files for the SDB, quality checking these files, and 37 transferring the files to the SDB. 6. Routinely checking that the parameters and embedded grids that are currently on line (i.e., in the data stream) are suitably accurate, generating updates for the SDB as required. 7. Coordinating with the satellite control center, on the scheduling of maneuvers and special VISSR imaging sequences, and appropriately changing the standard operational schedule to best deal with such events. Clearly the human plays a critical role in the VIRGS opera- tion. While the computational tasks are fully automated as Fortran routines executable by simple key-ins, the human is required to supervise their execution, as a safeguard against erroneous measurements under normal circumstances and as a re- serve of judgement against abnormal circumstances. Given the number of contingencies that exist daily, the VIRGS implementa- tion appears to have an effective mix of human and computer talents. An important feature of the VIRGS implementation is the clos- ed loop formed when the VISSR data, appended with the VIRGS produced earth-location information (parameters and embedded grids), are ingested and displayed at VIRGS. Figure 13 por- trays the loop. When ingesting VISSR data, VIRGS is emulating any other direct readout user and hence has the capability to monitor the performance of the system, end to end. The VIRGS operator should be the first data user to become aware of deg- radations in the earth-location accuracy. In most cases, the operator will be able to correct them quickly through coordina- tion with SDB operators or satellite control center. The issue of optimal versus minimal sets of star and landmark position measurements is discussed in section IV-A as a part of the orbit/attitude determination software descriptions. The screening of observations, based on the confidence factor re- sulting from any sort of difficulty in making the measurements, is an intuitive manual task for the analyst. Although most 38 r i D O Q < D O DC > cc C3 Q < h- CO ° £ i u oc £ cc "J — — CO > Q 3 5 1- co 111 2 UJ s < UJ cc »"■ •/ DISPLAY / OF STARS LANDMARKS / GRIDS CO cc UJ O < o 5 UJ — CO - 2 §8 to Q UJ uj u. 2 CO — I CO CC UJ i- UJ S < cc < a. UJ Q D - — u — D CC 1- C5 F Q 1- cc < < ■H 39 situations are clear-cut, some require the collection of addi- tional measurements for a reliable solution. D. Rationale for Centralized Service Since the earth-location parameters generated by VIRGS at the central facility are appended to the stretched VISSR data, other direct readout users do not need to exercise their own capability to generate these parameters, _if the centrally gen- erated parameters are sufficiently accurate. A centralized service is cost effective since all data user can avoid the expense of developing and operating their own earth-location capability. The operators of VIRGS strive for the greatest possible accuracy and use the built-in safeguards as required to guarantee a high reliability. The users, knowing that the earth-location parameters are readily available in the data stream, need only to employ a standard geometery transformation routine. The time consuming processing to derive orbit, attitude, and misalignment parame- ters is performed for users by VIRGS. Doubt that 1 pixel accuracy is achievable in an operational environment exists because the demonstration that such accuracy is achievable has only been done on the McIDAS System at the University of Wisconsin and hence has occured only in a univer- sity environment. The skills to achieve this accuracy are em- bodied in one person, Mr. John T. Young, and these skills in their totality have not been transferred to any other in- stallation. However, since the mathematical algorithms cur- rently available on the VIRGS System for determining the satel- lite's orbit/attitude states and misalignment parameters are a refined version of the orbit/attitude model of the McIDAS used by Mr. Young, one would expect that, after software and proced- ual deficiencies are recognized and removed, 1 pixel accuracy will be achieved. The actual results discussed in Section IV-E, obtained on VIRGS by Larry Hambrick, lend strength to this expectation. One advantage of having these algorithms in the form of computer programs (as they are now available) is that an operator can objectively generate parameters describing 40 satellite's orbit/attitude state without having a deep under- standing of orbit/attitude relationships. An accuracy of 1 pixel, as used herein, is defined as the root-mean-square of the earth location errors, line and element errors taken separately, distributed over the earth disc throughout a 24 hour period. One pixel is being used here as the angular distance between successive pixels in a full reso- lution picture. With this definition one is not guaranteed that any given point on any given frame will be located to within 1 pixel. A specific point at a given time might be in error by several pixels. However, loosening this requirement further, say to two or more pixels rms, significantly increases the possibility of larger errors. In terms of ground distance, a one pixel error spans one kilometer at the subpoint but spans two kilometers 50° from the subpoint. Finally, it is noted that, although 1 pixel accuracy rms is probably acceptable to all current users, even sharper accuracies are possible with the precision alignment capabilities inherent in the data stream. The software currently available on the VIRGS System has been verified to approach 1 pixel accuracy. Even though routine operators are not now expected to attain this level of accura- cy, it appears likely that information transfer and training will improve the routine operation. The goal is to make the routinely achieved accuracy a characteristic of the software package alone. This can be accomplished by improving the soft- ware packages to reduce the reliance on an operators' s judge- ment and by eliminating the impact of an operator's bias through appropriate training. E. An Error Budget and Some Actual Results 1. Error Budget The accuracy goal, stated here for overall system perform- ance, is one full resolution visible pixel, which is a 21 jara- dian angular interval at the satellite that translates into 0.8 km at the subsatellite point. The natural unit of error meas- ure for the software, satellite, and SDB is the full resolution 41 angular interval or pixel. All elements of the earth location process involving computations and data handling have been de- signed and implemented to achieve accuracies of a small frac- tion of a pixel. However, the end to end process has several potential error sources that must be dealt with and which limit the ultimately achieveable accuracy. The potential error sources are exhaustively listed below. 1. Tolerances in the orbit and attitude modeling. 2. Tolerances in the geometry transformation. 3. Measurements of landmark and star positions. 4. Misassignment of the earth coordinates of the landmarks. 5. Tolerances on the satellite's sun detector. 6. Biases in the SDB scan-line interpolation and sun position referencing. 7. Timing jitter in the satellite and SDB. 8. Tolerances in the solution for the orbit, attitude and misalignment parameters. The timing jitter of the satellite and SDB and the satel- lite's sun detection error might be taken as the hard limits on accuracy. There are, however, 22 time-counts per full resolu- tion pixel. The jitter amounts to no more than a few counts (i.e., less than a quarter of a pixel) . The sun detection er- ror is the combination of random errors in the leading-edge- detection, which should be less than two or three counts, and a slowly changing offset due to changes in the sun's relative position in the sun sensor's slot-shaped field of view. The slowly changing offset can be determined frequently enough to be accounted for in the "roll" parameter to within a tolerance of a fraction of a pixel. Errors in the SDB ' s scan line inter- polation and sun position reference combined, appear to be less than a half pixel. The accuracy of the geometry formulation is known to be well within 1 pixel and is intended to be within 0.1 of a pixel but cannot be truly verified to that level. The orbit model and the attitude model are intended to be within 1 pixel. The operation of these models has been verified at SSEC to be 42 within a tolerance approaching 1 pixel. An effort has been made to keep these models commensurate with 1 pixel accuracy. Nutation of the satellite's spin axis can be modelled and de- termined, but has not as yet been accounted for in the earth- location parameters. Nutation, if not accounted for, can give errors of two or three visible lines (peak) , cycling over about 8 scan lines. The potential accuracy of star and landmark position measure- ments (i.e., observed image frame coordinates) is within 1 pixel under normal conditions. Human error and variations in human judgement on features, can combine to reduce the routine- ly achieveable accuracy to about 1 pixel rms. The earth coor- dinates of landmarks can be verified to within a 1 pixel toler- ance. The star coordinates are available from standard tables to an accuracy of better than 0.2 pixel (i.e., 0.5 arc seconds in earth centered celestial coordinates) . The solution for orbit, attitude and misalignment parameters is susceptible only to biases in the measurements since it incorporates least- square-error regression. When the errors discussed above are included, the resulting accuracy for a data user with a standard geometry transform routine that uses the documented parameters should be 1 pixel rms. This figure represents errors over the entire earth disc and throughout the day. The effect of spin axis nutation con- tributes in the worst of cases to more than a 1 visible pixel rms error in line. It's contribution to overall accuracy is however greatly reduced because of its high frequency cycle. Also the nutation dynamically deforms earth features but a human or computerized recognition process filters out much of its effect. 2. Actual Results Figure 14 presents the results of a determination of orbit, attitude, and misalignment parameters from star and landmark measurements. The results shown are not the best obtainable and in fact are on the verge of being unsatisfactory. This set was selected for discussion here because the nature of the 43 c sD » to to rH r i sO SO CN s0 O r- r> -r ,-. rH rH T (N c i O to Hf H- U") LO) n to to LO o rH rH D o o O r j CN CN f\| r j C 1 cj CN rH i rj < 1 " ) c i c i N CM CN CN M c I C 1 CN c i ■ i CI LO LO o, IO IO LO in LO LO 10 LO L0 LO LO o ■J n tO to to o to IO to M to to ro to to /) rH 1 1 "' I 1— 1 1 rH 1 »i i rH 1 r-1 1 7 rH i rH 1 rH 1 rH i rH 1 H O o ^c 00 CO » <* to to O CT *t Tf <* O O H -J C 1 CM |H ,-H to to LO o LO SO .0 r-~- r^ -O, S i 1 1 1 1 1 1 1 i i 1 P [-«) X, 00 o LO >£> «tf to CT> CN to LO rH o r- \D CTi M «r> ct. LO a. r-> LO to LO to o i i n so to -^r »* a w rH rH i CI 1 to ,— rH rH r-i j w J to r-- CM r i CM r- to oc) to rH i---, LO) (N LO 30 t -[ SO c CN CM C i i <-: to rH LO <* SO LO LO LO Q ^ 1 i 1 1 I 1 i 1 i— i hJ PJ o i— i rH rH to IO o o rH LO o rH O ,-H o Q c LO LO LO LO o LO O o O r > to to to to to to o to U ►J r/j o O o o o o o o o o o CT> o o CT, CO o o o o o o o o o o o CT> o o D £ LO LO LO LO to IO 1/-. LO LO LO LO LO LO LO Xj s r— 1 rH rH rH 1— 1 — 1 rH rH rH "<* «tf ** rH rH H X 1— 1 CM CN to LO SO \D sC 1-^ - o o ■ l X rH rH rH .H rH H CN CN ' ' Q a CO oo 00 CO CO 00 oo CO 00 00 SO 0C CO Xj Q sD vO vD SO ^D vO SO v£> sO vO so -r sO vD D Q CM r i r i C 1 r | CN r j r i ri CN C 1 - 1 r j CM ■ 1 >< CT) CT. cr. ct. cn CT. CT. © CTi CT, CT, ct. CT, CTi T, >H r- r- t -» r- r - r~ t- r r- r- r> > r-- r-. DO "sf •H- -et *tf r)- LO LO to O O to CN tO to" LO to CN CM to to o o © o © CT, CJ CN Cl CN CN CM OJ rH rH CN (N CN n CN rH CJ OJ CI CN ri CJ c J CN CN U0 LO LTi LO 10 LO, LO LO, LO LO LO) LO IO LO LO to to 1 1 to 1 l to 1 LO I to 1 to 1 to 1 to 1 to i to 1 to 1 to to 1 1 to to rH L0 LO CT, rH IO to M3 \D CN CN LO LO LO LO LO LO CN to c J CN •-1 .— . O CT © © rH rH 1 1 >-H rH •^r i SO l v0 I 1 i i l LO LO 1 1 00 LO CT. LO ri ^r ^D oo o r-- o r- t-~ © to 00 o rH o to CT, rH •^f t-- Cl ^c ^H CT, ^ rH CM 1 1 CO CN oo CT, t-^ r^ to CN rH ^H •—< •vT 00 00 vO LO LO !"^ rH oc to vO rH rH rH T— 1 r-< Cl LO OJ to ^T CN tO CN O rH rH i— 1 to a-^ CT, CT, CT. CT, CT, CT. CT. CT, CT, © © r-~ r~- r> l~v r> r- r-> r~- r- f^ r- r^ r-. t~- r-^ •* **• -t Tf i- «* •* Hf Tt ^t «* Cm u a en ■H ic". o H F- P -J o 10, O © UJ OS a. 44 various error contributions is illustrated. The measured star and landmark positions are listed and iden- tified by satellite number (SS, West GOES is #24) , year (YY) and Julian (DDD) in the first column, SSYYDDD. The GMT of the start of the image frames in which the particular measurement where made is given as hours-minutes-seconds in the second column, HHMMSS. The LCODE given in the third column identifies the measurements as; a landmark, in Baja, CA 301 a landmark, in islands of the South Pacific 351 visible star #1 (|3-Orion, "Rigel") 353 visible star #3 (o'-Canis Minor A, "Procyon") 354 visible star #4 (a-Aquilae, "Altair") 355 visible star #5 (Y-Orion, "Bellatr ix" ) The star a-Orion ( "Betelgeuse" ) is also usually visible in the images. The fourth and fifth columns are respectively the line and element residuals (differences) between the observed image positions and the positions specified by the solution. The units for the residuals are the full resolution visible pixel. The last two columns give the latitude and longitude, in degrees-minutes-seconds (DDDMMSS) , of the satellite's subpoint. The solution was obtained from the measurements on day 268 and predicted to day 269. The prediction is a linear extrapo- lation in time of the attitude and an evaluation of the orbit model with epoch parameters (for GMT on day 268) at the par- ticular times on day 269. The solution residuals are roughly 0.5 rms in line and 1.5 rms in element. The prediction residu- als are roughly 2.5 rms in line and 1.5 rms in element. The measurements are spread over nearly 24 hours and geographically dispersed with landmark being far above and to the right of the subpoint and landmark 301 slightly below and to the left of the subpoint. The different stars vary in right ascension by about 220° and in declination between ±8°. The largest component of the prediction error, the "bulge" in line residuals during the period 17 to 22 hours GMT, suggests 45 an imprecise solution for the precesion rate of the satellite's attitude. The combination of all the other error sources dis- cussed above is seen as the 1- to 2-pixel variation that is somewhat random. 3. A Chronology of Results with VIRGS When the VIRGS became operational in May, 1979, the orbit/at- titude determination software was at an intermediate stage of its development. The level of accuracy achieved in the opera- tion with this version of the software during the first several months of operation was sporadic, ranging from 3 pixels rms on some days to 10 pixels on other days. On an average of once per week during this "break-in" period, severely degraded per- formance occurred as a result of procedural and training defi- ciences in dealing with non-routine matters such as satellite maneuvers. The intermediate version of the software, while still being used in the operation, was unable to adequately distinguish the orbit plane from the attitude of the spin axis using landmark and earth edge measurements. By August 1, 1979 the software had been revised by Dennis Phillips for the "star technique" described in this document. Preliminary tests of the revised software verified that the attitude and misalignment parameters could be accurately determined from the star positions alone, leaving the landmark positions for an independent orbit determination. During the months of August, September, and October 1979, the star technique's operational viability was tested by Larry Ham- brick. While having to work around the regular operation on VIRGS, image sectors containing stars were ingested from both the East and the West GOES. The star positions were measured and entered on the navigation files. The attitude and misa- lignment parameters were determined from the star measure- ments. With these parameters and the landmark measurements (i.e., from the ongoing operation), the orbit was determined. The accuracy of the predictions using these solutions were then checked. About fifty daily navigations were performed in this 46 manner. A formal demonstration was held in September, 1979. The following conclusions are drawn from the results of those tests: i. A sufficient set of stars, for the star technique, are routinely available in the images, ii. Using the star technique, VIRGS has a stand-alone orbit/attitude determination capability, iii. Several deficiencies remain in the revised software (see next paragraph); however, accuracies better than 3 pixels rms are routinely achievable on the predictions, using the star technique. At the time of these tests a discrepancy in "beta", apparent- ly between VIRGS and the SDB, existed at the level of about 2 pixels. Also there was a discrepancy, internal to VIRGS, be- tween the attitude and orbit routines as manifested in their residual outputs. The computation of the attitude's precession rates was done manually, a cumbersome task to do precisely. Work by Dennis Phillips subsequent to the tests is expected to have remedied all of these problems. A separate type of deficiency that had to be dealt with dur- ing these tests was a discrepancy in the earth coordinates of landmarks. Landmarks not on the American Continents appeared to have inconsistent earth coordinates with respect to American landmarks. The coordinates for all land features had been ob- tained from official detailed maps. An extensive test was con- ducted by Larry Hambrick on VIRGS (using the star technique) to establish the size of the discrepancy for about a dozen non- American land features. For some islands in the South Pacific, the discrepancy was as high as five image pixels. A discrepancy of at least two pixels was found for most of the non-American features. Based on this test it was concluded that a calibration of the coordinates of non-American features with respect to American land features was needed and would further improve the accuracy of VIRGS. 47 V. SUMMARY The accurate earth location of GOES/VISSR image data requires a complete mathematical formulation that is founded on the basic geometry. The use of stars together with landmarks as the observables in the images affords a small, efficient soft- ware package that employs relatively simple algorithms and models for determining the true misalignment angles, attitude, and orbit. The capability described in this document has been implemen- ted on VTRGS, the central VISSR earth-location facility oper- ated by NOAA/NESS and on the Mcldas at the University of Wis- consin. Direct readout users have convenient access to the VIRGS generated earth-location parameters, and by using a stan- dard geometry transform routine, can locate the VISSR data automatically. The NESS standard grids, also produced by VIRGS, are embedded in the VISSR data at the infrared pixel resolution. The VIRGS should be operated at the best achievable level of accuracy. A commitment to an accurate, reliable centralized service would save individual users a sizable task. The users face the choice of exercising their own VIRGS-type operation or simply using the earth-location parameters of the centralized service. 48 REFERENCES Mottershead, T. and Phillips, D. , 1976: Image Navigation for Geosynchronous Meteorological Satellites, presented at 7th Conference on Aerospace and Aeronautical Meteorology and Symposium on Remote Sensing from Satellites, Florida Institute of Technology, Melbourne, FL. Phillips, D., and E. A. Smith, 1974: Geosynchronous Satellite Navigation Model, in Studies of the Atmosphere Using Aerospace Probes 1974, Space Science and Engineering Center, University of Wisconsin Press, Madison, WI . 3. Proceedings of Interactive Video Displays for Atmospheric Studies, Sponsored by NSF, held at the University of Wisconsin June, 1977. 4. Smith, E. A., 1975: The Mcldas system, IEEE Trans. Geosci, Electronics , pp. 123-136. 5. Smith, E. A., and Phillips, D. R. , July, 1972, Automated Cloud Tracking Using Precisely Aligned Digital ATS Pictures, IEEE Trans, on Computers , pp. 715-729. 6. Westinghouse Electric Corporation, Report Entitled "Earth Location Equations," Contract NAS 5-23582, July 20, 1977. 49 APPENDIX A CONVERTING A POINTING VECTOR TO LINE AND ELEMENT NUMBERS The pointing vector in the first satellite centered coordi- nate system is designated as k. The projection of k onto the z-axis of the first satellite coordinate system is simply k»( - S), where S is the spin axis vector. In terms of scan line number and the misalignment parameters, the projection of the pointing direction of the spin scan camera onto the z-axis of the first satellite coordinate system is simply the z com- ponent of (2.5) or m^ , ocos ( (c* - H) . r, ) + m 3 3 #sin( (Cn-{) t») . (A.O) To determine the line number at which the spin scan camera is pointing in the direction k, simply equate (A.O) to k»(-S) and solve for 12 . Solving this equation requires a trigometric manipulation. 2 2 Vi Divide (A.O) by (m 3 1 +ir. 3 3 ) , set 9=ATAN2(m 3 x ,m 3 3 ) and 2 2/2 substitute sin(0) for itu ,/(m 3 . +m r . 3 ) and 2 ' 2 !/2 cos (9) for m 3 3 /(m 3 , +m-> 3 ) . The result is sin (9) o cos( (c^, -St ) «r« ) + cos (9) • sin( (c» - H) • r* ) =r- (-"S)/dn 3fl +m 3 ^ 3 ) 1 /2 Setting 4> = ASIN (k. (-"§) / (m^ x +m 3 3 ) 2 ) , using the sum angle sine angle formula on the left-hand side and applying the arcsine function to both sides gives: (Cg-2) .r^ + = .) ) (B.l) i=l 2 2 2 subject to the constraint a +b +c =1 and where for atti- tude determination cj> . equals TT /2+r -. ( c, -1 . ) , see figure 8, and for the orientation of the orbital plane all the , ' s equal 1T /2. The constraint is absorbed by setting c = 2 2 Vo -(l-(a +b ) ) n . Starting with a and b equal to zero sets the initial spin axis vector guess and the orbital plane per- pendicular guess pointing opposite to the spin axis of the earth. This is a good operational approximation and, conse- quently, reduces the number of iterations required during computation. The change of variables modifies the problem to minimizing the following: S(a,b) =2 (ax.+b yi -(l-(a 2 +b 2 ) ) /2 .z i -cos(4> i ) ) 2 . (B.2) i=l The first partials of S with respect to a and b are set equal to zero resulting in the following two equations: k 2 (x^a^ aJ • (ax i +by i +cz i -cos( 4^) ) =0 1 ^ 1 (B.3) 2 (y.-bz./ ) . (ax . +by . +cz . -cost §• ) ) =0 ._, J 1 l/c 1 2 1 1 1 2 2 ] A Note the use of the symbol c here for -(l-(a +b )) . Ma- nipulation yields the following matrix equation for the indi- vidual terms of the summation: 2 2 x. -z. +z .cos( <)>. ) /c x.y. 1111 i J l 2 2 , cos (.)-cx.z.+a x.z. /c+aby . z . /c l l 11 11 2 i l 2 y.cos(4 , -) - cy.z. +abx . z . /c+b y . z . /c jr 1 t 1 / : L 1 x ■}_' J i i' B-l The matrix on the left is inverted and the resulting system of equations is solved iteratively. Note that all the terms on the right-hand side outside of the constant terms are second order. The advantage of such a separation is that the linear equation is solved for first and then the results are updated with the second order terms. This approach is a variation of Newton's Method. The summation signs have been omitted to com- pact the presentation. Convergence is achieved within three or four iterations. For attitude determination with an unknown pitch misalignment angle, C, the following expression of three variables has to be minimized: k 2 S(a, b, c, £ )= s (ax +by +cz -cos(4> -£)) (B.5) ,111 i i=l 2 2 2 subject to the constraint a +b +c =1. The expression (B.5) is not easily amenable to the variation of Newton's Method used earlier to find the attitude and hence the expression is reformulated to an equivalent problem. The *2 2 2 constraint a +b +c =1 is absorbed by setting 2 2 x h 2 !/2 c= -(l-(a +b )) and u is set to sin(C) and -(1-u ) is set to cos(C) • Once u is solved for, C can be found as the 2 V2 arctangent of (-u/(l-u ) ) , or in Fortran as 2 Vi ATAN2 (u,- (1-u )')• With this formulation, the following is minimized: S(a,b,u)= Z (ax.+by.-(l-(a 2 +b 2 ) z.+ucos(4>.) + (l-uV 2 sin(4>. ) ) 2 . 1=1 (B-6) To solve this minimization problem the first partials of S with respect to a, b, and u are set equal to zero and the resulting expressions are rearranged to have a matrix times the vector (a, b, u) on the left-hand side and constant terms plus second order terms on the right-hand side", as in (B-4) . The inverse of the matrix is then applied to both sides and the resulting equations are solved by iteration using Newton's Method. Con- vergence is again achieved within three or four iterations. B-2 APPENDIX C ACCURACY OF APPROXIMATE ORBIT EQUATION (Refer to Section III-D, Equation 4.4) In this paper the equation n(t-T)= Y - 2esinY (C.l) has been used to approximate the motion of a satellite in a Keplerian orbit. Equation C.l is deriveable from (4.4) with the substitution Y = a - co. In eq. (C.l) , n is the mean motion constant, t is time, T is the time of a perifocus passage, Y is the true anamoly and e the eccentricity. The accuracy of this equation for eccentricities less than or equal to 0.001 will now be established. Consider the following three equations taken from Escabol: n(t-T) = E - € sin E, (C.2) = A-- 2 sin E = vi-€ sinY (C.3) 1+fcosY and cos E = cos Y + £ (C.4) 1+ecos Y where E is the eccentric anamoly. Multiplying (C.3) by cos Y and (C.4) by sin Y and subtracting yields; cos Y sin E - sin Y cos E = A- 2 e sinY cos Y t sinYcos Y - e » sinY (C.5) 1 + c • cos Y Using the sine angle difference formula changes that left hand side to sin (E-Y) . The next two inequalities are well know for |E-Y| within the range of consideration: Isin (E-Y)|>|E-Y| (l- (E-Y) 2 ) 3! and (C.6) Isin (E-Y) - (E-Y) 1<1 E-Y) 3) 3! (C.7) Using (C.5) and restricting e to be less than or equal to 0.001 yields Isin(E-Y) |<0. 001002. Combining this with (C.6) yields |E-Y|<. 001003 C-l Finally, (C.7) yields I sin (E-Y) - ( Now (C5) is rewritten as Isin (E-Y) - (E-Y) |< 1. 682xl0~ 10 (C.8) r~ 2 2 vl-e sinYcosY-sinYcosY+ € sinYcosY sin(E-Y)=-«sinY+ =— r ^7— — — 1 + f»cos Y Combining (C.9) with the inequality (C.8) and appropriately bounding the last term in (C.9) yields |E-Y + e sinYl < 1.6 x Next, (C.3) is rewritten as + e sinYl < 1.6 x 10" 6 (C„ 10) I 2 v l- e sin V- sin Y - e » cosYsinY /r . , -, v „ • v \ C • 1 1 ) sin E - sin « = r — tt~~ -—-—_ — _____ 1 + e . cos Y This yields the inequality Isin E - sinYl < .001002 (C.12) Combining (1.2) with the inequalities (CIO) and (C.12) yields In(t-T) - Y + 2e sinYl < 2.7 x 10~ 6 (C.13) This yields a positional accuracy to within 0.12 kilometers. Thus the accuracy of the approximation of using (C.l) is established. C-2 APPENDIX D I GAUSSIAN ELIMINATION The process of Guassian elimination is now recalled for the reader. One begins with a matrix pair (A:I) , where A is the matrix whose inverse is sought and I is the identity matrix, and applies a sequence of non-singular matrices B. to both members of the matrix pair from the left until B B ,...B,A is the identity matrix. At that point the n n-1 1 ■* product accumulated on the right, B B n ...B, , is the n n-1 1 inverse of the matrix A since its product with A yields the identify matrix. Usually the B.'s are selected so that the matrix multiplication affects only one row of the matrix pair Because the matrix (4.12) is symmetric, positive definite, i.e., a. .=a. ., and all the eigenvalues are positive, the matrix operands B. can be selected to act on more than one row in one operation. In this case the B.'s act on two rows. The matrix (4.12) is treated as if it were partitioned into four parts, each part consisting of a two by two matrix. Gaussian elimination is then applied as if each two by two submatrix was a single element. Start with the matrix pair /A, , A, \ /I 0\ and apply Next apply/ I -A, 2 \ with the result Vo I A l,l" A l,2 A 2,2 A 2,l °\ I 1 ~ A 1,2 A 2,2 A 2,2 A 2,l V \° A 2,2 D-l Then set S = (A, , - A 1 2 A 2 „ A^ ) and apply /S O _ 1 \o I A 1,2 A 2,2 "I A 2,2 Finally > apply/ I O\to both matrices with the result -1 " A 2,2 A 2,l S " SA 1,2 A 2,2 -1 -1 -1 -1 ~ A 2,2 A 2,1 S A 2,2 +A 2,2 A 2 , 1 SA 1, 2 A 2 , 2 The inverse of (4.12) is the right-hand term. Actually, com- puting the lower left-hand term is unnecessary since it is the transpose of the upper right term by symmetry. D-2 Appendix E Fortran Listings of Orbit and Attitude Determination Routines: FINDAT and UPGORB 16.050575 PAGE 1 SUBROUTINE FINDAT( NOPCLN .RADLIN .DECLIN .RASCEN , PCLN ) DOUBLE PRECISION UX.UY.UZ.PSI COMMON/MINCOM/UX(60) ,UY (60 ) , UZ (60 ) , PS I ( 80) , NP DATA RDPDG/1.745329252E-2/ C WRITTEN NOVEMBER 13,1979 BY DENNIS PHILLIPS TO SIMPLIFY ATTITUDE C COMPUTATION. PROGRAM RIGHTS BELONG TO SCIENTIFIC PROGRAMMING C AND APPLIED MATHEMATICS , INC. C THIS PROGRAM ELIMINATES THE USE OF STEEPEST DESCENT METHODS IN C THE ATTITUDE COMPUTATIONS AND INSTEAD REPLACES THESE METHODS WITH C LINEAR REGRESSION METHODS COMBINED WITH ITERATION STEPS. C TWO PROBLEMS ARE WORKED OUT IN THIS CODE. ONE PROBLEM IS TO C MINIMIZE THE SUM OVER I OF (A*X ( I ) + B*Y ( I )+C*Z ( I )-COS (PSI (I ) ) )**2. C WITH THE CONSTRAINT A**2+B**2+C**2=l . C THE OTHER PROBLEM IS TO MINIMIZE THE SUM OVER I OF C (A*X(I.)+B*Y(I) + C*Z(I)-C0S(PSI(I)-PSIREF))**2 OR FQUIVALENTLY THE C SUM OVER I OF ( A*X( I )+B*Y ( I ) +C*Z ( I )-V*COS ( PSI ( I ) ) -U*SIN ( PSI ( I ) ) )**2 C WITH THE CONSTRAINTS A**2+B**2+C**2 AND U**2+V**2=l. C *** C FOR PICTURE CENTER LINE SOLUTION GO TO 30 Q ### IF(NOPCLN.EQ.l)GO TO 30 C RETURN AND SEND ERROR MESSAGE FOR INSUFFICIENT MEASUREMENTS IF(NP.LT.2)G0 TO 75 C COMPUTE ATTITUDE C *#* C COLLECT REGRESSION COEFFICIENTS SMXSQ=0.0 SMXY=0.0 SMXZ=0.0 SMXCS=0.0 SMYSQ=0.0 SMYZ=0.0 SMYCS=0.0 SMZSQ=0.0 SMZCS=0.0 DO 10 1=1, NP X=UX(I) Y=UY(I) Z=UZ(I) ANG=PSI(I) CS-COS(ANG) SMXSQ=SMXSQ+X*X SMXY=SMXY+X*Y SMXZ=SMXZ+X*Z SMXCS=SMXCS+X*CS SMYSQ=SMYSQ+Y*Y SMYZ=SMYZ+Y*Z SMYCS=SMYCS+Y*CS SMZSQ=SMZSq+Z*Z 10 SMZCS=SMZCS+Z*CS C SET INITIAL VALUES AND UPDATE VALUES A=0.0 16.050575 PAGE B = 0„0 C— 1.0 ANEW=0.0 BNEW=0.0 SET ITERATION COUNT ITER=0 COMPUTE MATRIX COEFFICIENTS FOR MATRIX EQUATION (All (A21 A12) A22) = Fl = F2 20 A11=SMXSQ+SMZCS/C-SMZSQ A12-SMXY A21=SMXY A22=SMYSQ+SMZCS/C-SMZSQ C INVERT THE MATRIX A DET=1.0/(A11*A22-A12*A21 ) B11=A22*DET B12=-A12*DET B21=B12 B22=A11*DET COMPUTE CONSTANT PLUS SECOND ORDER TERMS OF THE VECTOR F F1=SMXCS-C*SMXZ+ANEW**2*SMXZ/C+ANSW*BNEW*SMYZ/C F2=SMYCS-C*SMYZ+BNEW**2*SMYZ/C+ANEW*BNEW*SMXZ/C C COMPUTE NEW A AND B VALUES ANEW=B11*F1+B12*F2 BNEW=B21*F1+B22*F2 C=-SQRT(1.0-(ANEW**2+BNEW**2 ) ) C COUNT ITERATION STEPS ITER=ITER+1 C RETURN AND SEND ERROR MESSAGE IF ITERATION LIMIT IS EXCEEDED IF(ITER.GE.10)GO TO 80 C CHECK CONVERGENCE CRITERIA IF(ABS (A-ANEW)+ABS(B-BNEW) .LT.1.0E-8)GO TO 25 C UPDATE VALUES A=ANEW 3=BNEW C ITERATE AGAIN GO TO 20 25 DECLIN=-90.0+ASIN(SQRT(A**2+B**2))/RDPDG RASCEN=0.0 IF(ABS(A)+ABS(B) .GT . 1 . 0E-8 )RASCEN=ATAN2 (B , A) /RDPDG RETURN C *#* C COMPUTE ATTITUDE AND PICTURE CENTER LINE OFFSET C *** 30 IF(NP.LT.3)G0 TO 85 C COLLECT REGRESSION COEFFICIENTS SMXSQ-0.0 SMXY=0.0 SMXZ=0.0 SMXCS-0.0 SMXSN=0.0 SMYSQ=0.0 SMYZ=0.0 SMYCS-0.0 18.050575 PAGE 3 SMYSN=0.0 SMZSQ-0.0 SMZCS=0.0 SMZSN=0.0 SMCSSQ=0.0 SMCSSN=0.0 SMSNSQ=0.0 DO 35 1=1, NP X = UX(I)' Y=UY(I) Z=UZ(I) ANG=PSI (I) SN=SIN(ANG) CS-COS(ANG) SMXSQ=SMXSQ+X*X SMXY=SMXY+X*Y SMXZ=SMXZ+X*Z SMXCS=SMXCS+X*CS SMXSN=SMXSN+X*SN SMYSQ=SMYSQ+Y*Y SMYZ=SMYZ+Y*Z SMYCS=SMYCS+Y*CS SMYSN-SMYSN+Y*SN SMZ3Q=SMZSQ+Z*Z SMZCS=SMZCS+Z*CS SMZSN=SMZSN+Z*SN SMCSSQ=SMCSSQ+CS*CS SMCSSN=SMCSSN+CS*SN 35 SMSNSQ=SMSNSQ+SN*SN C SET INITIAL AND UPDATE VALUES A = 0.0 B=0.0 C=-1.0 ANEW=0.0 BNEW=0.0 U=0.0 V = 1.0 UNEW=0.0 C SET ITERATION COUNT ITER=0 C COMPUTE COEFFICIENTS FOR MATRIX EQUATION (All A12 A13) A = Fl C (A21 A22 A23) B = F2 C (A31 A32 A33) U = F3 40 A11=SMXSQ-SMZSQ+V*SMZCS/C A12=SMXY A13=-SMXSN A21-A12 A22=SMYSQ-SMZSQ+V*SMZCS/C t A23=-SMYSN A31-A13 A32=A23 A33=C*SMZCS/V-SMCSSQ+SMSNSQ r #** 18.050575 PAGE 4 C THE INVERSION OF THE MATRIX A IS DONE BY PARTITIONING THE MATRIX C INTO FOUR PARTS: (All), (A12 A13) , (A21) AND (A22 A23) C (A31) (A32 A33) C THEN EACH MEMBER OF THE PARTITION IS TREATED AS A SINGLE ELEMENT C AND GAUSSIAN ELIMINATION IS CARRIED OUT STARTING IN THE LOWER LEFT C HAND CORNER Q *** C FIND THE INVERSE OF THE LOWER LEFT HAND CORNEP DET=1.0/(A22*A33-A23*A32) CU=A33*DET C22=A22*DET C12=-A23*DET C21=C12 C MULTIPLY- THIS INVERSE BY THE LAST TWO ENTRIES OF THE FIRST ROW B1=-A12*C11-A13*C21 B2=-A12*C12-A13*C22 C FIND FIRST ELEMENT OF INVERSE MATRIX CDIVID=1.0/(A11+B1*A21+B2*A31 ) B11=CDIVID C FIND THE SECOND AND THIRD ELEMENTS OF FIRST ROW OF MATRIX B12=B1*CDIVID B13=B2*CDIVID C SET TWO MORE ENTRIES BY SYMMETRY B21=B12 B31 = B13 C COMPUTE THE REMAINING TWO BY TWO PART B22=C11+CDIVID*B1*B1 B23=C12+CDIVID*B1*B2 C BY SYMMETRY B32=B23 B33=C22+CDIVID*B2*B2 C COMPUTE RIGHT-HAND SIDE OF LINEAR EQUATION F1=V*SMXCS-C*SMXZ+A**2*SMXZ/C+A*B*SMYZ/C-A*U*SMZSN/C F2=V*SMYCS-C*SMTZ+A*B*SMXZ/C+B**2*SMYZ/C-B*U*SMZSN/C F3=-V*SMCSSN+C*SMZSN-A*U*SMXCS/V-B*U*SMYCS/V+U**2*SMCSSN/V C COMPUTE NEW A, B AND U VALUES ANEW=B11*F1+B12*F2+B13*F3 BNEW=B21*F1+B22*F2+B23*F3 UNEW=B31*F1+B32*F2+B33*F3 ITER=ITER+1 C RETURN AND SEND ERROR MESSAGE IF ITERATION LIMIT IS EXCEEDED IF(ITER.GE.15)G0 TO 90 C CHECK CONVERGENCE CRITERIA IF(ABS(A-ANEW)+ABS(B-BNEW)+ABS(U-UNEW) .LT . 1.0E-8)GO TO 45 C UPDATE VALUES A=ANEW B=BNEW C=-SQRT(1.0-(A*A+B*B)) U=UNEW V=SQRT(1.0-U**2) C ITERATE AGAIN GO TO 40 C COMPUTE DECLINATION, RIGHT ASCENSION AND PICTURE CENTER LINE OFFST 18.050575 PAGE 5 45 DECLIN=-90.0+ASIN(SQRT(A*A+B*B) )/RDPDG RASCEN=0.0 IF(ABS(A)+ABS(B).GT.1.0E-8)RASCEN=ATAN2(B,A)/RDPDG PCLN=PCLN-ATAN2(U,V)/RADLIN RETURN C SEND ERROR MESSAGES 75 WRITE(6,95) 95 FORMAT (5X, 'INSUFFICIENT MEASUREMENTS FOR ATTITUDE DETERMINATION') RETURN 80 WRITE(6,96) 96 F0RMAT(5X, 'ITERATION LIMIT EXCEEDED IN ATTITUDE DETERMINATION') RETURN 85 WRITE(6,97) 97 F0RMAT(5X, 'INSUFFICIENT MEASUREMENTS FOR ATTITUDE AND PCLN DET . ' ) RETURN 90 WRITE(6,98) 98 F0RMAT(5X, 'ITERATION LIMIT EXCEEDED IN ATTITUDE AND PCLN DET.') RETURN END SIZE 945 01661 16.050575 PAGE 1 C UPGORB PHILLI 067S NAVL COMPUTE ORBIT FROM LANDMARKS , BETAS , AND ATTITU0004 C C THIS PROGRAM HAS TWO MODES: ONE MODE TO COMPUTE ALL THE ORBIT C PARAMETERS AND THE OTHER MODE TO COMPUTE JUST THE ALONG-TRACK C ORBIT PARAMETERS ASSUMING THE ORBIT'S INCLINATION AND C ASCENDING NODE HAVE BEEN PROVIDED. MIN(3) CONTROLS THE C NUMBER OF ALONG-TRACK ORBIT PARAMETERS TO BE COMPUTED. C MIN(6) IS THE FLAG CONTROLING THE ORBITAL PLANE CALCULATION. C C MIN(l) = SSYYDDD C MIN(2) = HHHH t THE FIRST TWO DIGITS ARE THE FIRST HOUR ON THE YYDD C TO START. THE LAST TWO DIGITS ARE THE NUMBER OF C HOURS FROM THE FIRST HOUR TO CONSIDER. C MIN(3) = NO. OF PARAMETERS TO COMPUTE. C MIN(4) = EPOCH OF THE CALCULATED ORBIT PARAMETERS. 00 C MIN(5) = PRINT AND STORE OPTION. C MIN(6) = ORBITAL PLANE CALCULATION FLAG. OP MEANS ON. SUBROUTINE MAIN 0012 REAL MEANOM,PTIMLM(100) ,XLINLM(100) f XELELM(100) 0013 REAL XLATLM(100) ,XLONLM(100) 0014 DOUBLE PRECISION BETA( 100 ) ,V ( 100) ,T (100 ) 0015 INTEGER HHHH 0016 INTEGER MIN(S) ,LCODE(100) ,ILNDBT(100) 0017 INTEGER LNDPNT(100) 0018 COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0019 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN0020 2, PRERAT.PREDIR, PITCH, YAW, ROLL, SKEW 0021 DATA MIN/'UPGORB',6*0/ 0022 DATA MAXLND/100/,NUMDAY/1/ 0023 CALL IQ(MIN) 0024 ISYD=MIN(1) 0025 HHHH=MIN(2) 0026 NUMDAYMHHHH/100+MOD (HHHH, 100) )/24+l 0027 C THE CALL TO SETUPN INITIALIZES THE PARAMETERS IN THE NAVCOM AND C NAVIN COMMON BLOCKS CALL SETUPN(ISYD,0) 0028 C FETCH THE LANDMARKS FOR THE TIME PERIOD CALL GTLM(ISYD,NUMDAY,MAXLND,NUMLND, 0029 * PTIMLM,LCODE,XLINLM,XELELM,XLATLM,XLONLM) 0030 C THIS CALL TO GTLMBT MATCHES BETA VALUES TO THE AVAILABLE LANDMARKS C BY INTERPOLATING OR EXTRAPOLATING THE BETA VALUES THAT ARE C AVAILABLE IN THE SAME IMAGE WITH THE LANDMARKS. WHEN SUCH BETAS C ARE AVAILABLE FOP A SPECIFIC LANDMARK, THE CORRESPONDING INDEX C POSITION OF THE LANDMARK IN THE ARRAY ILNDBT IS SET TO ONE? C OTHERWISE THIS VALUE IS SET TO A MINUS ONE. CALL GTLMBT ( I SYD,NUMDAY, HHHH, NUMLND ,PTIMLM .XLINLM ,XELELM , LCODE, * ILNDBT, BETA, T) C COUNT LANDMARKS HAVING ASSOCIATED BETA VALUES IMATCH=0 0033 DO 2 1=1, NUMLND 0034 IF( ILNDBT (I ) .EQ.1)IMATCH=IMATCH+1 0035 2 CONTINUE 0036 C TRANSMIT ERROR MESSAGE AND RETURN WHEN NOT ENOUGH LANDMARKS HAVE 18.050575 PAGE ASSOCIATED BETA COUNTS IF(IMATCH.GE.MIN(3))G0 TO 3 0037 CALL TQMESCNOT ENOUGH MATCHES BETWEEN LANDMARKS AND BETAS $', IMATC003S *H) 0039 RETURN 0040 CONTINUE 0041 IF(MIN(3).GE.1.AND.MIN(3) .LE.4)GO TO 4 0042 CALL TQMES( 'WRONG NUMBER OF ORBIT PARAMETERS TO BE COMPUTED??', 0043 *MIN(3)) 0044 RETURN 0045 ISEMI=MIN(3) 0046 EPTIME=FTIME(MIN(4)) 0047 IHNS=MIN(4) 0048 IORBPL=0 0049 IF(MIN(6),EQ.2H0P)IQRBPL=1 0050 THE ITERATION LOOP (DO 6 1=1,3) CORRECTS THE HEIGHT OF THE SATELLITE AS A FUNCTION OF TIME USING THE MOST RECENTLY COMPUTED VALUES FOP THE ORBIT PARAMETERS. WHEN NO ORBIT PARAMETERS ARE AVAILABLE, THE NOMINAL GEOSTATIONARY HEIGHT OF 42165.0 IS USED. THIS HEIGm. CORRECTION IS DONE IN ANGORB IETIMH=MIN(4) IETIMY=MOD( I SYD, 100000) ITER=5 MEANOM=0.0 DO 6 1=1, ITER THE CALL TO ANGORB SERVES TWO FUNCTIONS: (1) IF IORBPL=l , THE INCLINATION AND ASCENDING NODE OF THE ORBITAL PLANE ARE DETERMINED AND (2) THE TRUE ANAMOLY POSITIONS + OFFSET OF THE SATELLITE ARE COMPUTED. THE TRUE ANAMOLIES+OFFSET ARE STORED IN THE ARRAY V(I). CALL ANGORB ( I SYD, NUMLND ,MEANOM ,PTIMLM ,XLATLM,XLONLM , *XLINLM,XELELM,ILNDBT,BETA,V,T,NV,LNDPNT,IORBPL) 0059 COMPUTE THE ALONG-TRACK ORBIT PARAMETERS. THE NUMBER OF PARAMETERS COMPUTED DEPENDS ON THE INDEX ISEMI . HENCE I SEMI EQUALS ONE, TWO, THREE OR FOUR. 0054 0055 0052 CALL ORBPR(V,T,NV,MEANOM, ISEMI, EPTIME) 6 CONTINUE IOUT=l COMPUTE RESIDUALS FOR LANDMARK MEASUREMENTS USING THE NEWLY COMPUTED ORBITAL STATE AND OUTPUT RESIDUALS TO TERMINAL. CALL ORBRES(ISYD,IOUT,NV,T,XLINLM,XELELM,BETA,MEANOM, * LCODE,XLATLM,XLONLM,LNDPNT) I0UT=2 OUTPUT RESIDUALS TO PRINTER IF PRINTER OPTION IS ON. 0060 0061 0062 0064 0065 18.050575 PAGE 3 IF(MIN(5).FQ.2HPS.0R.MIN(5) .EQ . 2HSP .OR . MI N ( 5 ) .EQ. 1HP ) CALL ORBRES 0066 *(ISYD,IOUT,NV,T,XLINLM,XELELM,BETA,MEANOM,LCODE,XLATLM.XLONLM,LNDP *NT) C C STORE NAVIGATIONAL PARAMETERS IN NAVFILES IF STORAGE OPTION IS ON. C- IF(MIN(5).EQ.2HPS.0R.MIN(5).EQ.2HSP.0R. MI N( 5) .EQ . 1HS ) CALL BQ 0068 *(ISYD,IHMS,SEMIMA,OECCEN,ORBINC ,MEANOM , PERHEL , ASNODE ) 0069 IOUT=l 0070 C C OUTPUT RESULTANT ORBIT PARAMETERS TO THE TERMINAL CRT. C CALL ORBOUT(IOUT,ISYD,IHMS,SEMIMA,OECCEN,ORBINC,M£a!\0;i,?L'RHEL, 0071 * ASNODE) 0072 I0UT=2 0073 C C IF PRINTER OPTION IS ON, ALSO OUTPUT ORBIT PARAMETERS TO PRINTER. C IF(MIN(5).E0.2HPS.OR.MIN(5) . EQ. 2HSP .OR . MI N( 5) .EQ . IflP )CALL O'RBOUT 0074 * (IOUT, IS YD, IHMS, SEMI MA, OECCEN .ORBIN C ,MEANOM , PERHEL , ASNODE) 0075 RETURN 0076 END 0077 SIZE 2221 04255 16.050575 PAGE SUBROUTINE OEBOUT ( IOUT , IS YD , IHMS .SEMI.'IA ,OECCEN ,ORBINC ,MEANOM, * PERHEL,ASNODE) THE SUBROUTINE ORBOUT TRANSMITS THE NEWLY COMPUTED ORBIT PARAMETERS TO THE TERMINAL CRT OR TO THE PRINTER DEPENDING ON THE VALUE OF IOUT, IOUT=l MEANS TERMINAL. IOUT =2 MEANS PRINTER. INPUTS: ALL PARAMETERS OUTPUTS: NO NEW PARAMETERS OR VALUES ARE RETURNED TO CALLING PROGRAM. THE ORBITAL PARAMETERS ARE TRANSMITED TO THE OUTPUT DEVICE SPECIFIED BY IOUT. REAL MEANOM INTEGER M0UT(132) ,IARRAY(8) CALL YDDMY(ISYD,ID,IM,IY) IARRAY(l)=(IY-1900)*10000+lM*100+ID IARRAY(2)=IHMS IARRAY(3)=IROUND(SEMIMA*100.0) IARRAY(4)=IROUND(OECCEN*1000000.0) IARRAY(5)=IROUND(ORBINC*1000.0) IARRAY(6)-IRQUND(MEANOM*1000.0) I ARRAY ( 7 )=IROUND( PERHEL*1000 .0 ) IAPRAY(8)=IROUND(ASNODE*1000.0) TP TRANSMITS A HOLLERETH FIELD TO THE OUTPUT DEVICE, CRT OR PRINTER 0078 0079 0060 0081 0082 0063 0084 0085 0086 0087 0088 0089 0090 CALL TP(I0UT,132H UPGORB OUTPUT RESULTS CALL TP(I0UT,132H KEPLERIAN ORBIT PARAMETERS MVCHAR MOVEW A HOLLERETH FIELD INTO THE ARRAY MOUT FIRST BLANK OUT FIELD CALL MVCHAR ( 'BLA' , 'BLA' , MOUT ,1,132) WRITE FIELD LABELS CALL MVCHAR( *' ETIMY ETIMH SEMIMA ECCEN ORBINC MEANA PERIGEE *,1, MOUT, 1,64) CALL TP(IOUT,MOUT) BLANK FIELD AGAIN CALL MVCHAR ( 'BLA' , 'BLA ', MOUT , 1 , 132 ) CALL TP(IOUT,MOUT) DO 10 1=1,8 ENCODE ORBIT PARAMETERS IN INTEGER HOLLERETH FIELD CALL MVCHAR(IARRAY(I ) , 'INT ', MOUT , 1*8, 1 ) 10 CONTINUE TRANSMIT HOLLERETH FIELD TO OUTPUT DEVICE CALL TP(IOUT,MOUT) RETURN 0091 0092 )0093 0094 0095 ) 0096 0097 0098 ASNODE'0099 0100 0101 0102 0103 0104 0105 0106 0107 0108 18.050575 PAGE 2 END 0109 SIZE 400 00620 18.050575 PAGE 1 SUBROUTINE DQ( ISYD ,IHMS , SEMI MA ,OECCEN , ORB I NC , MEANOM .PERHEL ,ASNODE ) 0110 C C DO TRANSMITS THE ORBITAL PARAMETERS TO THE NAVIGATIONAL FILES C C INPUTS: ORBITAL PARAMETERS C OUTPUTS: THE ORBIT PARAMETERS ARE TRANSMITED TO THE NAVIGATIONAL C FILES BY NAVFIL. C REAL MEANOM 0111 INTEGER MESONE(10) ,MESTWO(10) 0112 DATA MESONE/'NAVFIL', 4,7*0/ 0113 DATA MESTWO/'NAVFIL',5,7*0/ 0114 CALL YDDMY(ISYD,ID,IM,IY) 0115 I ETYMD=(IY-1900)* 10000+1 M*100+ID 0116 ISEMI=IROUND(SEMIMA*100.0) 0117 I OECC=IROUND(OECC EN* 1000000.0) 0118 IASND=IROUND(ASNODE*1000.0) 0119 I 01 NC = I ROUND (ORBING* 1000.0) 0120 IPERH=IROUND(PERHEL*1000.0) 0121 I MEAN- IROUND( MEANOM* 1000.0) 0122 MES0NE(4)=ISYD 0123 MES0NE(6)=1 0124 MES0NE(7)=IETYMD 0125 MES0NE(8)=IHMS 0126 MES0NE(9)=ISEMI 0127 MES0NE(5)=MES0NE(9) 0128 MESONE(10)=IOECC 0129 CALL SQ(MESONE) 0130 MESTW0(4)=ISYD 0131 MESTW0(6)=1 0132 MESTW0(7)=I0INC 0133 MESTW0(8)=IMEAN 0134 MESTW0(9)=IPERH 0135 MESTW0(5)=MESTW0(9) 0136 MESTWO(10)=IASND 0137 CALL SQ(MESTWO) 0138 RETURN 0139 END 0140 SIZE 139 00213 18.050575 PAGE 1 SUBROUTINE GTLMBT ( ISYD ,NUMDAY ,HHHH , NUMLND , PTI MLN , XLINLM, XELELM f * LCODE,ILNDBT,BETA,T) C C GTLMBT GETS BETA COUNTS WITH CALLS TO GETBET AND CHECKS TO SEE C IF THE BETA COUNT PAIRS OCCUR AT A PICTURE TIME OF ANY OF THE C LANDMARKS. IF SO, A BETA VALUE FOR THE LANDMARK IS FOUND BY C INTERPOLATION OR EXTRAPOLATION. ALSO, THE LANDMARK CODE IS C EXAMINED TO DISCARD STARS. C INPUTS: ISYD. SATELLITE ID. YEAR DAY C NUMDAY. THE NUMBER OF DAYS OF LANDMARK MEASUREMENTS TO C CHECK. C HHHH. THE BEGINNING HOUR AND ENDING HOUR TO CONSIDER C LANDMARK MEASUREMENTS. C NUMLND: THE NUMBER OF LANDMARKS THAT CAN BE HANDLED. C OUTPUTS: PTIMLM. PICTURE START TIME OF LANDMARK C XLINLM. LINE NUMBER OF LANDMARK C XELELM. ELEMENT NUMBER OF LANDMARK C LCODE. LANDMARK CODE C ILNDBT. LANDMARK-BETA ASSOCIATION FLAG. EQUALS 1 WITH C BETA. EQUALS -1 WITHOUT BETA. C T. LANDMARK TIME. C IMPLICIT DOUBLE PRECISION (Q) DOUBLE PRECISION BETA( 1 ) ,RDPBT , BETDIF, T( 1 ) ,BETCMP REAL PTIMLN(l) ,XLI NLM( 1 ), XELELM ( 1 ) INTEGER HHHH 0147 INTEGER LCODE(l) ,ILNDBT(1) 0148 INTEGER NOMTIM(100) 0149 INTEGER ISCAN1 ( 100 ) , ISTIM1 ( 100 ) , ISMIL1 ( 100) ,IBET1 (100) 0150 INTEGER ISCAN2(100),ISTIM2(100) ,ISMIL2(100) ,IBET2(100) 0151 COMMON/NAVCOM/NAVDAY,LlNTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,0ECCEN,0RBINC,PERHEL,ASN0DE,N0PCLN,DECLIN,RASCEN,PICLIN 2, PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW COMMON/NAVINI/ 0152 1 EMEGA,AB,ASQ,BSQ,R,RSQ, 0153 2 RDPDG, 0154 3 NUMSEN,TOTLIN,RADLIN, 0155 4 TOTELE,RADELE,PICELE, 0156 5 CPITCH,CYAW,CROLL, 0157 6 PSKEW, 0158 7 RFACT,ROASIN,TMPSCL, 0159 8 B11,B12,B13,B21,B22,B23,B31,B32,B33, 0160 9 GAMMA, GAMDOT, 0161 A R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33, 0162 B PICTIM,XREF 0163 DATA MAXNUM/100/,NEGBET/3l44960/,RDPBT/9.989292881D-7/ 0164 DATA BETCMP/3144960.0D0/ 0165 DATA TWPI/6. 28318530/ ISAT=ISYD/100000 0166 IYR=MOD(ISYD, 100000) /1000 0167 IDAY=MOD(ISYD,1000) 0168 C COMPUTE BEGINNING AND END TIME FOR TIME INTERVAL CHECK BTIME=HHHH/100 0169 18.050575 PAGE 2 ETIME=BTIME+MOD(HHHH,100) 0170 DO 10 I=1,NUMLND 0171 10 ILNDBT(I)=-1 0172 C THE BETAS ARE FETCHED ONE DAY AT A TIME. DO 50 I=1,NUMDAY 0173 JSYD=ISAT*100000+IYR*1000+IDAY 0174 CALL GETBET(JSYD,MAXNUM,NUM,NOMTIM, 0175 * ISCAN1,ISTIM1,ISMIL1,IBET1, 0176 * ISCAN2,ISTIM2,ISMIL2,IBET2) 0177 IF(NUM.EQ.0) GO TO 40 0178 DO 30 J=1,NUM 0179 C COMPUTE BETA TIME USING THE I INDEX TO COUNT THE HOURS OF EACH DAY PTIME=FTIME(NOMTIM(J))+(I-1)*24.0 0180 C IF TIME LIES OUTSIDE SET INTERVAL FOR ORBIT COMPUTATION, EXIT. IF(PTIME.LT.BTIME.OR.PTIME.GT.ETIME)GO TO 30 0181 C THE BETA PAIRS ARE TAKEN ONE AT A TI^E USING THE INDEX J AND ALL C THE LANDMARKS ARE LOOKED THROUGH USING THE INDEX K. C DO 25 K=1,NUMLND 0182 C CHECK EACH LANDMARK OR STAR FOR MATCH. C IF IT IS A STAR, DISCARD. IF(MOD(LCODE(K),100) .GE.5)G0 TO 25 0183 C IF THE TIMES DON'T MATCH, DISCARD. IF(PTIME.NE.PTIMLN(K))GO TO 25 0164 C IF LANDMARK HAD PRIOR MATCH, DISCARD. IF(ILNDBT(K).EQ.l)GO TO 25 0185 C OTHERWISE LINEARLY INTERPOLATE OR EXTRAPOLATE MATCHING BETA C AND TIME VALUES AND SET MATCH ARRAY ILNDBT EQUAL TO ONE AT INDEX K IF(IBET1(J).LT.0) IBET1(J)=(IBET1(J).AND. '37777777) +NEGBET 0186 IF(IBET2(J).LT.0) IBET2 (J )=( IBET2 (J ) .AND . '37777777 ) +NEGBET 018? BETDIF=IBET2(J)-IBET1(J) 0188 IF(BETDIF.GT.BETCMP)BETDIF=BETDIF-2.0D0*BETCMP 0189 IF(BETDIF.LT.-BETCMP)BETDIF=BETDIF+2.0D0*BETCMP 0190 ISCNLN=(IR0UND(XLINLM(K))+3)/8-ISCANl(J)+l QTIME1=FTIME(ISTIM1(J) ) +( 1-1 )*24. 0D0 QTIME2=FTIME(ISTIM2(J) )+( 1-1 )*24 .0D0 QDTIME=QTIME2-QTIME1 QTM1ML=ISMIL1(J)/360000.0D0 QTM2ML=ISMIL2(J)/360000.0D0 QDTIME=(QTM2ML-QTM1ML)+QDTIME QTIME1=QTIME1+QTM1ML QTMPSL=QDTIME/(ISCAN2(J)-ISCAN1(J)) QBTPSL=BETDIF/(ISCAN2(J)-ISCAN1(J)) QPCTLN=(XELELM(K)-1.0)/(TOTELE-1.0)*(DEGELE/360.0) T(K)=QTIME1+QTMPSL*(ISCNLN+QPCTLN) BETA(K)=(IBET1(J)+QBTPSL*(ISCNLN))*RDPBT T(K)=T(K)+QTMPSL*BETA(K)/TWPI C LANDMARK POINTER ARRAY WHICH IS USED FOR RESIDUAL COMPUTATION. ILNDBT (K)=l 0197 25 CONTINUE 0198 30 CONTINUE 0199 C INCREMENT DAY WHILE ACCOUNTING FOR LEAPYEAR 16.050575 PAGE IDAY=IDAY+1 IF(IDAY.LE.365)G0 TO 40 IF(IDAY.EQ.366.AND.MOD(IYR,4).EQ.0)GO TO 40 IDAY=1 IYR=IYR+1 IYR=MOD(IYR,100) 40 CONTINUE 50 CONTINUE RETURN END SIZE 1317 02445 18.050575 PAGE SUBROUTINE ANGORB( ISYD .NUMLND .MEANOM ,PTIME , XLAT ,XLON , * XLIN.XELE.ILNDBT.BETA.V.T.NV.LNDPNT.IORBPL) 0205 WRITTEN 04/17/78 BY DENNIS PHILLIPS TO SET UP GEOMETRY TO FIND ORBITAL POSITION ANGLE, I. E. TRUE ANAMOLY+OFFSET POSITON, IN THE ORBITAL PLANE OF THE SATELLITE. MODIFICATION BY DENNIS PHILLIPS ON MAY 2, 1979 TO ACCOUNT FOR THE 0209 VARIATION OF THE SATELLITE'S HEIGHT. 0210 MODIFICATION 07/17/78 BY DENNIS PHILLIPS TO DETERMINE ORIENTATION OF THE ORBIT PLANE FROM LANDMARK MEASUREMENTS IN 0212 ADDITION TO THE ALONG-PLANE ORBIT PARAMETER DETERMINATION. 0213 INPUTS: ISYD. SATELLITE ID. YEAR DAY NUMLND. NUMBER OF LANDMARKS MEANOM MEAN ANAMOLY OF SATELLITE POSITION AT EPIC PTIME. PICTURE START TIME OF LANDMARK XLAT. LATITUDE OF LANDMARK XLON. LONGITUDE OF LANDMARK XLIN. LINE NUMBER OF LANDMARK XELE. ELEMENT NUMBER OF LANDMARK ILNDBT. LANDMARK-BETA CORRESPONDENCE FLAG BETA. BETA ANGULAR SWEEP TO PICTURE START IN RADIANS T. TIME OF LANDMARK MEASUREMENT IORBPL. ORBIT PLANE CALCULATION FLAG OUTPUTS: V. THE TRUE ANAMOLIES+OFFSET NV. NUMBER OF LANDMARKS WITH CORRESPONDING BETA VALUES LNDPNT. INDEX POINTING TO LANDMARKS WITH MATCHING BETAS ORBINC. ORBIT INCLINATION IN NAVCOM COMMON ASNODE. ASCENDING NODE IN NAVCOM COMMON. INTEGER LNDPNT(l) 0214 INTEGER ILNDBT(l) 0215 REAL MEANOM REAL PTIME(l),XLAT(l),XLON(l),XLlN(l),XELE(l) 0216 DIMENSION ESTST1(100),ESTST2(100),ESTST3(100) 0217 DOUBLE PRECISION XCOR.YCOR 0218 DOUBLE PRECISION BETA( 1 ) ,V (1 ) ,T (l ) 0219 DOUBLE PRECISION COSB.SINB.BETAR.TEMP , YCSNST .XCSNST .XNORM 0220 DOUBLE PRECISION SAMTIM 0221 DOUBLE PRECISION UX.UY.UZ.PS I DOUBLE PRECISION XSN.YSN.ZSN 0223 COMMON/MINCOM/UX(80) ,UY(80) ,UZ(80) ,PSI (80) ,NP 0225 COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0226 lIMH,SEMIMA,OECCEN t ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN0227 2 .PRERAT .PREDIR , PITCH , YAW .ROLL .SKEW 0228 COMMON/NAVINI/ 0229 1 EMEGA,AB,ASQ,BSQ,R,RSQ, 0230 2 RDPDG 0231 3 NUMSEN.TOTLIN.RADLIN, 0232 4 TOTELE.RADELE.PICELE, 0233 5 CPITCH.CYAW.CROLL, 0234 6 PSKEW, 0235 7 RFACT.ROASIN.TMPSCL, 0236 18.050575 PAGE 2 8 D11,D12,D13,D21,D22,D23,D31,D32,D33, 0237 9 GAMMA ,GAMDOT, 0238 A R0TMll t R0TM13,R0TM21,R0TM23 t R0TM31,R0TM33, 0239 B PICTIM.XREF 0240 COMMON/BETCOM/IAJUST r IBTCON,NEGBET,ISEANG 0241 DATA EH/42165.0/ 0242 DATA PI/3.14159265/ 0243 DATA IENT/1/ C SET UP PARAMETERS BEFORE LANDMARK LOOP IYD=MOD(ISYD, 100000) 0244 NV=0 0245 C INIALIZE ORBIT PARAMETERS' IE THE SEMI MAJOR AXIS IS AVAILABLE. NAVDAY=NAVDAT-1 IF(SEMIMA.GT.42000.0)CALL STVEC (0 .0 ,MEANOM ,XC ,YC , ZC ) NAVDAY=NAVDAY+1 SENANG=FLALO(ISEANG)*RDPDG 0247 C THIS IS THE LANDMARK LOOP. FOR EACH LANDMARK WITH MATCHING BETA C VALUES AN ANGULAR POSITION OF THE SATELLITE IN ITS ORBITAL PLANE C IS COMPUTED. DO 2 I=1,NUMLND 0263 IF(ILNDBT(I) .EQ.-DGO TO 2 0264 C COUNT MATCHES BETWEEN LANDMARKS AND BETAS NV=NV+1 0265 C TIME VALUES CORRESPONDING TO ORBITAL ANGULAR POSITIONS ARE ONLY C SORTED ON THE FIRST ENTRY. OTHERWISE THESE VALUES GET SHUFFLED C AND LOST. IF(IENT.EQ.1)T(NV)=T(I) TIME=T(NV) SAMTIM=T(NV) C LANDMARK POINTER ARRAY WHICH IS USED FOR RESIDULA COMPUTATION. LNDPNT(NV)=I 0269 C PRECESS ATTITUDE CALL ATPREC(TIME,DEC,RAS) C SET UP SATELLITE CENTERED COORDINATE SYSTEM DETERMINED BY C DECLINATION AND RIGHT ASCENSION DEC=DEC*RDPDG RAS=RAS*RDPDG SR=SIN(RAS) 0250 CR=COS(RAS) 0251 SD=SIN(DEC) 0252 CD=COS(DEC) 0253 X1=-SR 0254 Y1=CR 0255 Z1=0.0 0256 X2=-SD*CR 0257 Y2=-SD*SR 0258 Z2=CD 0259 X3=CD*CR 0260 Y3=CD*SR 0261 Z3=SD 0262 C FIND VECTOR FROM EARTH'S CENTER TO LANDMARK IN INERTIAL COODINATES C TO LANDMARK. JTIME=ITIME(TIME) 0271 18.050575 PAGE RA-RAERAC(IYD,JTIME,XL0N(I))*RDPDG YLAT=XLAT(I)*RDPDG YLAT=GE0LAT(YLAT,1) TANLAT=TAN(YLAT)**2 RR=SQRT((l.0+TANLAT)/(BSO+ASQ*TANLAT) )*AB X=COS(RA)*COS(YLAT)*RR Y=SIN(RA)*COS(YLAT)*RR Z=SIN(YLAT)*RR C GET VECTOR FROM EARTH TO SUN CALL SUNVEC ( IYD ,SAMT IM , XS N ,YSN ,ZSN ) C SET SATELLITE VECTOR TO ZERO FOR FIRST ITERATION XSAT=0.0 YSAT=0.0 ZSAT=0.0 H=HH IF(SEMIMA.LT.1.0)GO TO 1 C OTHERWISE COMPUTE THE SATELLITE VECTOR IN INERTIAL COORDINATES C AND USE HEIGHT TO IMPROVE POSITION ESTIMATE. CALL STVEC(TIME,MEANOM,XSAT,YSAT,ZSAT) H = SORT (XSAT**2+YSAT**2+ZSAT**2 ) 1 CONTINUE C COMPUTE X, Y COORDINATES OF THE VECTOR FROM THE SATELLITE TO THE C SUN IN OUR SATELLITE CENTERED COORDINATE SYSTEM. XCSNST=X1*(-XSN-XSAT)+Y1*(-YSN-YSAT) YCSNST=X2*(-XSN-XSAT)+Y2*(-YSN-YSAT)+Z2*(-ZSN-ZSAT) THE NEXT SIX STATEMENTS YIELD THE POINTING DIRECTION OF THE SPIN SCAN CAMERA IN A COORDINATE SYSTEM WITH THE Z-AXIS COINCIDIN( WITH THE VECTOR OPPOSITE THE SPIN AXIS AND THE X-AXIS PERPENDICULAR TO THE Z-AXIS AT AN ANGLE EQUAL TO SENANG FROM THE SUN PULSE DETECTOR. YLIN=(PICLIN-XLIN(I))*RADLIN CLIN=COS(YLIN) SLIN=SIN(YLIN) U=R0TM11*CLIN+R0TM13*SLIN V¥=R0TM21*CLIN+R0TM23*SLIN W=R0TM31*CLIN+R0TM33#SLIN THE ANGULAR SWEEP DISTANCE OF THE SPIN SCAN CAMERA FROM THE SUN AT THE TIME WHEN THIS PARTICULAR LANDMARK IS BEING VIEWED. THE (VV.U) TERM ACCOUNTS FOR THE ROLL AND YAW MISALIGNMENT EFFECTS. BETAR=BETA(I)-SENANG+XELE(I)*RADELE-ATAN2(VV,U) THE NEXT 8 STATEMENTS COMPUTE A UNIT VECTOR POINTING AT THE LANDMARK IN THE SATELLITE CENTERED COORDINATE SYSTEM. COSB=DCOS(BETAR) SINB=DSIN(BETAR) TEMP=COSB*XCSNST+SINB*YCSNST YCSNST=-SINB*XCSNST+COSB*YCSNST XCSNST=TEMP XNORM=DSQRT( (1 .0D0-W**2) /(XCSNST**2+YCSNST**2) ) XCSNST=XCSNST*XNORM YCSNST=YCSNST*XNORM C FINALLY, WE HAVE THE POINTING DIRECTION OF THE SPIN SCAN CAMERA C IN INERTIAL COORDINATES. XLNDST=-XCSNST*X1-YCSNST*X2-W*X3 0272 0273 0274 0275 0276 0277 0278 0279 0280 0281 0282 0283 0284 0285 0287 0292 0293 0294 0295 0296 0297 0298 0299 0300 0301 0302 0303 0304 0305 0306 0307 0308 0309 0310 18.050575 PAGE 4 YLNDST=-XCSNST*Y1-YCSNST*Y2-W*Y3 0311 ZLNDST=-XCSNST*Z1-YCSN*T*Z2-W*Z3 0312 C COMPUTE SLANT RANGE (w-xitNG) TO LANDMARK FROM SATELLITE. COSANG=XLNDST*X+YLNDST*Y+ZINDST*Z 0313 SLTRNG=-C0SANG+SQRT(H**2-RR**2+C0SANG**2) 0314 C THE INTERSECTION OF THE RAY EXTENDING FROM THE LANDMARK WITH A C SPHERE CENTERED AT THE EARTH'S CENTER AND RADIUS EQUAL TO H. ESTST1(NV)=X+XLNDST*SLTRNG 0315 ESTST2(NV)=Y+YLNDST*SLTRNG 0316 ESTST3(NV )=Z+ZLNDST*SLTRNG 0317 IF(IORBPL.NE.l)GO TO 2 0318 C THIS VECTOR IS NORMALIZED AND THE ANGLES PSI(I) ARE SET WHEN THE C ORIENTATION OF THE ORBITAL PLANE IS TO BE CALCULATED. YNORM=1.0/SQRT(ESTST1(NV)**2+ESTST2(NV)**2+ESTST3(NV)**2) 0319 UX(NV)=ESTSTl(NV)*YNORM 0320 UY(NV)=ESTST2(NV)*YN0RM 0321 UZ(NV)=ESTST3(NV)*YN0RM 0322 PSI(NV)=PI/2.0 2 CONTINUE 0324 IENT=0 C BRANCH IF THE INCLINATION AND ASCENDING NODE ARE NOT TO BE C CALCULATED. IF(IORBPL.NE.l)GO TO 3 0325 NP=NV 0326 NOPCLN=0 C CHANGE ON MARCH 22, 1980 BY DENNIS PHILLIPS TO REPLACE CALL TO C MINMIZ, THE OLD WAY TO COMPUTE ATTITUDE. INSTEAD FINDAT, A C SUBROUTINE USING A VARIATION OF NEWTON'S METHOD, IS CALLED TO C COMPUTE THE ORIENTATION OF THE ORBITAL PLANE. CALL FINDAT(NOPCLN ,RADLIN ,ORBNC ,ASNDE , PCLN ) ORBINC=90.0+ORBNC ASNODE=ASNDE+270.0 IF(ASNODE.GE.360.0)ASNODE=ASNODE-360.0 3 CONTINUE 0338 ASN=ASNODE*RDPDG 0339 CASN=COS(ASN) 0340 SASN=SIN(ASN) 0341 XINC=ORBINC*RDPDG 0342 CINC=COS(XINC) 0343 SINC=SIN(XINC) 0344 U1=CASN 0345 Vl-SASN 0346 W 1=0.0 034? U2=-SASN*CINC 0348 V2=CASN*CINC 0349 W2=SINC 0350 C C COMPUTE THE ORBITAL ANGLE POSITIONS OF THE SATELLITE AROUND C THE ORBIT PERPENDICULAR. CARE IS TAKEN TO AVOID A WRONG MULTIPLE C OF 2*PI. C DO 10 1=1, NV 0351 XC0R=U1*ESTST1 ( I ) +V1*ESTST2 ( I )+Wl*ESTST3( I ) 0352 18.050575 PAGE 7 10 YC0R=U2*ESTST1(I)+V2*ESTST2(I )+W2*ESTST3( I ) V(I)=DATAN2(YC0R,XC0R) IF(I.NE.1)G0 TO 5 ANG2^V(l) TIME=T(1) GO TO 7 CONTINUE K=IROUND((T(I)-T(l))/24.0+(ANG2-V(I ))/(2.0*PI)) V(I)=V(I)+2.0*PI*K CONTINUE CONTINUE RETURN END 0353 0354 0355 0356 0357 0358 0359 0360 0361 0362 0363 0364 0365 SIZE 1456 02660 18.050575 PAGE 1 SUBROUTINE STVEC ( SAMTIM .MEANOM ,X , Y , Z ) C COPIED 0460 BY DENNIS PHILLIPS FROM SATVEC. C STVEC COMPUTES THE INERTIAL POSITION (X, Y, Z) AT TIME SAMTIM. C INPUTS: SAMTIM ZULU TIME IN HOURS FOR WHICH POSITION IS REQUIRED C MEANOM MEAN ANAMOLY OF SATELLITE POSITION AT EPIC C OUTPUTS: X, Y, Z THE INERTIAL POSITION OF THE SATELLITE AT ZULU C TIME SAMTIM. DOUBLE PRECISION RDPDG ,RE,GRACON .DIFTI M ,EACAN1 .ECANOM , XMANOM DOUBLE PRECISION DABS ,DSQRT ,DS IN ,DCOS REAL MEANOM COMMON /NAVCOM/NAVDAY, LI NTOT.DEGLIN, IELTOT, DEGELE, SP INRA, IETIMY, I ET 1IMH,SEMIMA,0ECCEN,0RBINC,PERIGE,ASN0DE,N0PCLN,DECLIN,RASCEN,PICLIN 2, PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW DATA NAVSAV/0/ DATA GRACON,RE/.07436574D0,6378.388/ DATA RDPDG/. 0174532925D0/ C STVEC WILL RE-INITIALIZE ORBIT PARAMETERS IF NAVDAY IS CHANGED. IF(NAVDAY.EQ.NAVSAV) GO TO 10 0=RDPDG*ORBINC P=RDPDG*PERIGE A=RDPDG*ASNODE SO-SIN(O) CO=COS(0) SP=SIN(P)*SEMIMA CP=COS(P)*SEMIMA SA=SIN(A) CA=COS(A) PX=CP*CA-SP*SA*CO PY=CP*SA+SP*CA*CO PZ=SP*SO QX=-SP*CA-CP*SA*CO QY=-SP*SA+CP*CA*CO QZ=CP*SO SROME2=SQRT(1.0-OECCEN)*SQRT(1.0+OECCEN) XMMC=GRAC0N*RE*DSQRT(RE/SEMIMA)/SEMIMA EPSILN=1.0E-8 C COMPUTE THE DIFFERENCE IN TIMES OF ( NAVDAY , )-( IETIMY , IETIMH) . C THIS DIFFERENCE IS COMPUTED IN MINUTES. IEY=MOD( IETIMY /10 00, 100) INY=MOD(NAVDAY/1000,100) IYRMN=0 IF(IEY.EQ.INY)GO TO 5 IF(INY.LT.IEY)GO TO 3 IE=INY-1 DO 2 I=IEY,IE IYRMN=IYRMN+365*1440 IF(MOD(I,4).EQ.0)IYRMN=IYRMN+1440 2 CONTINUE GO TO 5 3 CONTINUE IE=IEY-1 DO 4 I=INY,IE IYRMN=IYRMN-365*1440 18.050575 PAGE IF(MOD(I,4).EQ.0)IYRMN=IYRMN-1440 4 CONTINUE 5 CONTINUE IED=MOD(IETIMY,1000) IND=MOD(NAVDAY,1000) IYRMN=IYRMN+(IND-IED)*1440 TE=FLOAT(IYRMN)-FTIME(IETIMH)*60.0 10 DIFTIM=SAMTIM*60.0+TE XMANOM=XMMC*DIFTIM+MEANOM*RDPDG ECANMl=XMANOM DO 20 1=1,20 ECANOM=XMANOM+OECCEN*DSIN(ECANMl) IF(DABS(ECANQM-ECANM1) .LT .EPSILN )GO TO 30 20 ECANMl=ECANOM 30 XQMEGA=DCQS(ECANOM)-OECCEN Y0MEGA=SR0ME2*DSIN(ECAN0M) X=XOMEGA*PX+YOMEGA*QX Y=XOMEGA*PY+YOMEGA*QY Z=XOMEGA*PZ+YOMEGA*QZ RETURN END SIZE 345 00531 18.050575 PAGE 1 SUBROUTINE ORBPR( V ,T ,N ,MEANOM , ISEMI ,EPTLME) C WRITEEN BY DENNIS PHILLIPS OF SCIENTIFIC PROGRAMMING AND APPLIED C MATHEMATICS INC. C ORBPR CONVERTS ANGULAR MEASUREMENTS (TRUE ANAMOLY+OFFSET) IN THE C ORBPR CONVERTS ANGULAR MEASUREMENTS (TRUE ANAMOLY+OFFSET) IN THE C ORBITAL PLANE WHICH HAVE ASSOCIATED TIME TAGS INTO KEPLERIAN C ORBITAL ELEMENTS. C INPUTS: V. TRUE ANAMOLY+OFFSET ANGLES. C T. TIME TAGS OF ANGLES C N. NUMBER OF ANGLES. C ISEMI. NUMBER OF KEPLERIAN ORBITAL PARAMETERS TO 3E C COMPUTED. C EPTIME. EPIC TIME FOR KEPLERIAN ORBITAL ELEMENTS. C OUTPUTS: MEANOM. MEAN ANAMOLY OF ORBIT AT EPIC TIME. C SEMIMA. ORBIT'S SEMIMAJOR AXIS IN NAVCOM COMMON BLOCK. C OECCEN. ORBIT'S ECCENTRICITY IN NAVCOM COMMON BLOCK C PERIGEE. ORBIT'S PERIGEE IN NAVCOM COMMON BLOCK. DOUBLE PRECISION T(1),V(1) REAL MEANOM COMMON/NAVCOM/NAVDAY.LINTOT ,DEGLIN , IELTOT .DEGELE , SPINRA , IETIMY , I ET 1IMH,SEMIMA,OECCEN,ORBINC,PERIGE,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 2,PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW DATA GRACON ,RE/ .07436574D0 ,6378.388/ DATA PI, RDPDG/3. 14159265,. 01745329252/ C COMPLETELY REWRITTEN SEPTEMBER 25, 1979 BY DENNIS PHILLIPS OF SPAA C M TO CLARIFY THE CODING METHODS AND TO CONFORM TO A STRUCTURED C PROGRAMMING PHILOSOPHY. ALSO, THE METHOD FOR COMPUTINGE THE C SEMIMAJOR AXIS IS MUCH IMPROVED. IF(ISEMI.NE.4)GO TO 20 C COMPUTE ALL FOUR ALONG-TRACK PARAMETERS BY LINEAR REGRESSION. CALL F0URCF(C1,C2,C3,C4,V,T,N) SEMIMA=(GRACON*C2*60.0)**(2.0/3.0)*PE CALL ORBTHR ( CI, C2,C3,C4, EPTIME, MEANOM) RETURN 20 IF(ISEMI.NE.3)G0 TO 30 C COMPUTE ORBIT'S MEAN ANAMOLY, ECCENTRICITY, AND ARGUMENT OR PERIGE C2=(SEMIMA/RE)**(3.0/2.0)/(60.0*GRACON) CALL THRC0F(C1,C2,C3,C4,V,T,N) CALL ORBTHR ( CI, C2,C3,C4, EPTIME, MEANOM) RETURN 30 IF(ISEMI.NE.2)G0 TO 40 C COMPUTE ORBIT'S SEMIMAJOR AXIS AND MEAN ANAMOLY WHILE ADJUSTING C FOR THE PERTURBATIVE EFFECTS OF A NONZERO ECCENTRICITY. C3=0.0 C4=0.0 DO 36 J=l,4 SME=0.0 SMESQ=0.0 SMT=0.0 SMET=0.0 XN=FLOAT(N) DO 35 1=1, N ANG=V(I) 18.050575 PAGE TT=T( I )-C3*SIN ( ANG)-C4*C0S (ANG) SME=SME+ANG SMESQ=SMESQ+ANG**2 SMT=SMT+TT 35 SMET=SMET+ANG*TT DET=1 .0/(XN*SMESQ-SME*SME ) C2=(SMET*XN-SME*SMT)*DET C3=-2.0*OECCEN*C2*COS(PERIGE*RDPDG) C4=2.0*C2*OECCEN*SIN(PERIGE*RDPDG) 36 CONTINUE C1=(SMT*SMESQ-SMET*SME)*DET SEMIMA=(GRACON*C2*60.0)**(2.0/3.0)*RE MEANOM= (EPTIME-C1 ) /( C2*RDPDG )-PERIGE MEANOM=AMOD(MEANOM, 360.0) IF(MEANOM.LT.0.0)MEANOM=MEANOM+360.0 RETURN 40 IF(ISEMI.NE.1)RETURN COMPUTE ORBIT'S MEAN ANAMOLY. C2=(SEMIMA/RE)**(3.0/2.0)/(60.0*GRACON) C3=-2.0*OECCEN*C2*COS(PERIGE*RDPDG) C4=2.0*C2*OECCEN*SIN(PERIGE*RDPDG) CALL 0NEPAR(C2,C3,C4,V ,T , N ,EPT 1MB, PERIGE.MEANOM) RETURN END SIZE 321 00501 18.050575 PAGE SUBROUTINE FOURCF ( CI ,C2 ,C3 ,C4 ,V ,T ,N ) C WRITTEN BY DENNIS PHILLIPS OF SCIENTIFIC PROGRAMMING AND APPLIED C MATHEMATICS, INC. C FOURCF USES LINEAR REGRESSION TO MINIMIZE THE SUM OVER I OF C (CI + C2*V(I) + C3*C0S(V(I)) + C4*SIN(V(I)) - T(I))**2. C INPUTS: TRUE ANAMOLY+OFFSET ANGULAR POSITIONS. C T. ANGULAR TIME TAGS. C N. NUMBER OF ANGULAR POSITIONS. C OUTPUTS: CI, C2, C3 AND C4 LINEAR PLUS SINE WAVE CONSTANTS OF C ORBITAL MOTION. DOUBLE PRECISION A (4 ,4 ) ,B (4 ,4 ) ,C (4 ) ,D (4 ) , V ( 1 ) ,T (l ) DO 5 1=1,4 JS = I D(I)=0.0 DO 5 J=JS,4 5 A(I,J)=0.0 C SUM REGRESSION TERMS INTO THE APPROPIATE ELEMENTS OF THE MATRIX A XK>N A(1,1)=XK DO 10 1=1, N ANG=V(I) TT=T(I) CANG=COS SANG=SIN A(1,2)=A A(2,2)=A A(1,3)=A A(1,4)=A A(2,3)=A A(2,4)=A A(3,3)=A A(3,4)=A A(4,4)=A D(1)=D(1 D(2)=D(2 D(3)=D(3 D(4)=D(4 10 CONTINUE C INVERT A CALL INVE(A,B) DO 20 1=1,4 C(I)=0.0 DO 20 K=l,4 20 C(I)=C(I)+B(I,K)*D(K) C1=C(1) C2=C(2) C3=C(3) C4=C(4) RETURN END (ANG) (ANG) (1,2)+ANG (2,2)+ANG**2 (1,3)+SANG (1,4)+CANG (2,3)+ANG*SANG (2,4)+ANG*CANG (3,3)+SANG**2 (3,4)+SANG*CANG (4,4)+CANG**2 ) + TT )+TT*ANG )+SANG*TT )+TT*CANG SIZE 266 00412 18.050575 PAGE SUBEOUTINE INVE(A,B) C WRITTEN BY DENNIS PHILLIPS C MATHEMATICS, INC. C INVERTS FOUR BY FOUR REGRES C INPUTS: FOUR BY FOUR MATRIX C OUTPUTS: B, THE INVERSE OF C A IS PARTITIONED INTO FOUR C RESULTANT TWO BY TWO MATRIX C BY USING GAUSSIAN ELIMINATI DOUBLE PRECISION A(4,4),B(4 DOUBLE PRECISION DET ,TEMP DET=1.0D0/(A(1,1)*A(2,2)-A( B(l,l B(l,2 B(2,l B(2 t 2 C(l,l C(2,l C(l,2 C(2,2 B(3,l B(4,l B(3,2 B(4,2 D(l,l D(l,2 D(2,l D(2,2 DET = 1 B(3,3 B(3,4 B(4,3 B(4,4 TEMP=B(3,1) B(3,l B(4,l TEMP=B(3 f 2) B(3,2 B(4,2 B(l,3 B(l,4 B(2,3 B(2,4 B(l,l B(l,2 B(2,l B(2,2 RETURN END =A(2,2)*DET =-A(l,2)*DET =B(l,2) =A(1,1)*DET =B(1,1)*A(1,3)+B(1,2) =B(2,1)*A(1,3)+B(2,2) =B(1,1)*A(1,4)+B(1,2) =B(2,l)*A(l,4)+B(2,2) =-A(l,3)*B(l,l)-A(2,3 =-A(l,4)*B(l,l)-A(2,4 =-A(l,3)*B(l,2)-A(2,3 =-A(l,4)*B(l,2)-A(2,4 =A(3,3)-A(1,3)*C(1,1) OF SCIENTIFIC PROGRAMMING AND APPLIED SION MATRIX. A. A. TWO BY TWO SUBMATRICES AND THEN THE CONSISTING OF SUBMATRICES IS INVERTED ON. ,4),C(2,2),D(2,2) 1.2)*A(1,2)) =A(3,4 =A(3,4 =A(4,4 0D0/(D =D(2,2 =-D(l,2)*DET =B(3,4! =D(1,1 =B(3,3 =B(4,3 =B(3,3 =B(4,3 =B(3,1 =B(4,1 =B(3,2 =B(4 t 2 =B(1,1 =B(1,2 =B(1,2 =B(2,2 ~A(l,3)*C(l,2) -A(1,4)*C(1,1) -A(1,4)*C(1,2) 1,1)*D(2,2)-D( *DET *A(2,3) *A(2,3) *A(2,4) *A(2,4) )*B(2,1) )*B(2,1) )*B(2-,2) )*B(2,2) -A(2,3)*C(2,1) -A(2,3)*C(2,2) -A(2,4)*C(2,1) -A(2,4)*C(2,2) 1,2)*D(2,1)) *DET *B(3,1)+B(3,4) *TEMP+B(4,4)*B *B(3,2)+B(3,4) *TEMP+B(4,4)*B *B(4,1) (4,1) *B(4,2) (4,2) -C(1,1)*B(3,1) -C(1,1)*B(3,2) -C(l,2)*B(4,i; -C(l,2)*B(4 t 2; -C(2,1)*B(3,2)-C(2,2)*B(4,2) SIZE 271 00417 Fl 18.050575 PAGE 1 SUBROUTINE ORBTHR(Cl ,C2 ,C3 ,C4 ,EPTIME, MEANOM) C PULLED OUT OF ORBPR BY DENNIS PHILLIPS ON MAY 2, 1979 0607 C CHANGED BY DENNIS PHILLIPS ON MAY 11, 1979 0608 C ANAMOLY AND EQUALS ONE WHEN THE COEFFICIENTS ARE COMPUTED FROM THE0610 C ECCENTRIC ANAMOLY. 0611 C COMPUTE ECCENTRICITY, ARGUMENT OF PERIGEE AND MEAN ANAMOLY. 0612 C COMPUTES ECCENTRICITY, ARGUMENT OF PERIGEE AND MEAN ANAMOLY FROM C CI, C2, C3 AND C4 COEFFICIENTS. C INPUTS: CI, C2, C3 AND C4 COEFFICIENTS. C EPTIME. THE EPIC TIME OF ORBIT COMPUTATION. C OUTPUTS: OECCEN . ORBIT ECCENTRICITY. C PERIGE.' ORBIT'S PEPIGEE C MEANOM. ORBIT'S MEAN ANAMOLY. C REAL MEANOM 0613 COMMON /NAVCOM/NAVDAY ,LI NTOT .DEGLIN , IELTOT ,DEGELE , SPINRA , IETIMY , I ET0614 1IMH,SEM IMA, OECCEN, ORBINC , PERIGE , ASNODE , NOPCLN .DECLIN .RASCEN ,PI CLIN0615 2,PRERAT,PREDIR,PITCH,YAW,RCLL,SKEW 0616 COMMON/NAVINI/ 0617 1 EMEGA,AB,ASC,BS0,R,RSO, 0616 2 RDPDG, 0619 3 NUMSEN,TOTLIN,RADLIN, 0620 4 TOTELE,RADELE,PICELE, 0621 5 CPITCH,CYAW,CROLL, 0622 6 PSKEW, 0623 7 RFACT,ROASIN,TMPSCL, 0624 8 B11,B12,B13,B21,B22,B23,B31 ,B32,B33, 0625 9 GAMMA, GAMDOT, 0626 A R0TM11,R0TM13,R0TM21 ,R0TM23 ,R0TM31 .R0TM33 , 0627 B PICTIM,XREF 0628 OECCEN=SGRT(C3**2+C4**2)/(2.0*C2^ PERIGE=0.0 PERIGE=ATAN2(C4,-C3) /RDPDG PERIGE=AMOD( PERIGE ,360.0) IF(PERIGE.LT.0.0)PERIGE=PERIGE+360.0 40 CONTINUE 0635 MEANOM=(EPTIME-Cl )/( C2*RDPDG ) -PERIGE MEANOM=AMOD(MEANOM, 360.0) 0637 IF (MEANOM. LT.0.0)MEANOM=MEANOM+360.0 0638 RETURN 0639 END 0640 SIZE 78 00116 18.050575 PAGE SUBROUTINE THRCOF( CI ,C2,C3 ,C4,E ,T ,N ) 0566 WRITTEN MAY 2,1980 BY DENNIS PHILLIPS FOR NESS. THRCOF FITS A SINE WAVE PLUS CONSTANT TO THE CURVE T (I )-C2*E(I ) . 0568 EFFECTIVELY A LEAST SQUARES FIT IS PERFORMED TO FIND Cl 9 C3 AND C40569 IN THE EXPRESSION C1+C2*E( I ) +C3*SIN (E( I ) ) +C4*COS(E( I ) )=T ( I ) . 0570 INPUTS: E. THE TRUE ANAMOLY + OFFSET SATELLITE POSITION ARRAY. T. THE TIME OF THE ESTIMATED SATELLITE POSITIONS. N. THE NUMBER OF ESTIMATED SATELLITE POSITIONS. C2. THE COEFFICIENT DETERMINING THE SATELLITE'S MEAN OUTPUTS: CI, C3 AND C4. THE CONSTANT PLUS SINE WAVE COEFFICIENTS DOUBLE PRECISION E(1),T(1) 0571 SMSN=0.0 0572 SMCS=0.0 0573 SMSNSQ=0.0 0574 SMSNCS=0.0 0575 SMCSSQ=0.0 0576 SMCN=0.0 0577 SMSNCN=0.0 0578 SMCSCN=0.0 0579 COLLECT REGRESSION COEFFICIENTS DO 10 1=1, N 0580 ANG=E(I) 0581 CN=T(I)-C2*ANG 0582 SN=SIN(ANG) 0583 CS=COS(ANG) 0584 SMSN=SMSN+SN 0585 SMCS=SMCS+CS 0586 SMSNSQ=SMSNSQ+SN*SN 0587 SMSNCS=SMSNCS+SN*CS 0588 SMCSSQ=SMCSSQ+CS*CS 0589 SMCN=SMCN+CN 0590 SMSNCN=SMSNCN+SN*CN 0591 SMCSCN=SMCSCN+CS*CN 0592 10 CONTINUE 0593 XN=N 0594 INVERT THRREE BY THREE REGRESSION MATRIX USING CRAMER'S RULE. DET=XN*(SMSNSQ*SMCSSQ-SMSNCS**2)-SMSN*(SMSN*SMCSSQ-SMCS*SMSNCS) 0595 *+SMCS*(SMSN*SMSNCS-SMCS*SMSNSQ) 0596 DET=1.0/DET 0597 C1=(SMCN*(SMSNSQ*SMCSSQ-SMSNCS**2)-SMSNCN*(SMSN*SMQSSQ-SMCS*SMSNCS0598 *)+SMCSCN*(SMSN*SMSNCS-SMCS*SMSNSQ))*DET 0599 C3=(XN*(SMSNCN*SMCSSQ-SMSNCS*SMCSCN )-SMSN* (SMCN*SMCSSQ-SMCS*SMCSCN06 00 *)+SMCS*(SMCN*SMSNCS-SMCS*SMSNCN) )*DET 0601 C4=(XN*(SMSNSQ*SMCSCN-SMSNCN*SMSNCS)-SMSN*(SMSN*SMCSCN-SMCN*SMSNCS0602 *)+SMCS*(SMSN*SMSNCN-SMCN*SMSNSQ) )*DET 0603 RETURN 0604 END 0605 SIZE 250 00372 18.050575 PAGE 1 SUBROUTINE ONEPAR ( C2 ,C3 ,C4 f V ,T , N t EPTIME .PERIGE ,MEANOM ) C WRITTEN BY DENNIS PHILLIPS OF SCIENTIFIC PROGRAMMING AND APPLIED C MATHEMATICS, INC. C COMPUTE THE ORBIT'S MEAN ANAMOLY BY REGRESSION. C INPUTS: C2, C3 AND C4 LINEAR AND SINE WAVE CONSTANTS OF ORBITAL C MOTION. C V. TRUE ANAMOLY + OFFSET ANGULAR POSITIONS. C T. ANGULAR TIME TAGS. C N. NUMBER OF ANGULAR POSITIONS. C EPTIME. EPIC TIME OF DESIRED MEAN ANAMOLY POSITION. C PERIGE. ARGUMENT OF PERIGEE. C OUTPUTS: MEANOM: MEAN ANAMOLY. DOUBLE PRECISION V(1),T(1) REAL MEANOM DATA RDPDG/. 01745329252/ SUM=0.0 XN-N DO 10 1=1, N ANG=V(I) 10 SUM=SUM+T(I)-C2*ANG-C3*SIN(ANG)-C4*C0S(ANG) C1=SUM/XN MEANOM=(EPTIME-Cl )/(C2*RDPDG )-PERIGE MEANOM=AMOD (MEANOM ,360 .0 ) IF ( MEANOM. LT.0.0)MEANOM=MEANOM+360.0 RETURN END SIZE 100 00144 18.050575 PAGE 1 SUBROUTINE ORBRES ( ISYD , IOUT , NUMLND, T,XLI N ,XELE , BETA , MEANOM 9 * LCODE,XLAT,XLON,LNDPNT) 036? C C ADAPTED 10/05/78 BY DENNIS PHILLIPS FROM SUBROUTINE RESIDU 0368 C OUTPUTS LINE AND ELEMENT RESIUDES FROM UPGORB 0369 C C INPUTS: ISYD. SATELLITE ID. YEAR DAY C IOUT. OUTPUT DEVICE. ONE MEANS CRT. TWO MEANS PRINTER. C NUMLND. THE NUMBER OF LANDMARKS TO COMPUTE .RESIDUALS FOR C T. THE TIME OF EACH LANDMARK MEASUREMENT. C XLIN. THE LINE NUMBER OF EACH LANDMARK MEASUREMENT. C XELE. THE ELEMENT NUMBER OF EACH LANDMARK MEASUEEMNT. C BETA. THE SWEEP ANGLE UP TO THE BEGINNING OF THE SCAN C LINE ON WHICH THE LANDMARK WAS MEASURED. C MEANOM MEAN ANAMOLY OF SATELLITE POSITION AT EPIC C LCODE. THE LANDMARK CODE. C XLAT. THE LATITUDE OF THE MEASURED LANDMARK. C XLON. THE LONGITUDE OF THE MEASURED LANDMARK. C LNDPNT. AN INDEX ARRAY POINTING AT LANDMARKS WITH C ASSOCIATED BETA COUNTS. C INTEGER LCODE(l),LNDPNT(l) 0370 REAL MEANOM REAL XLAT(l) f XLON(l),XLIN(l>,XELE(l) 0371 DOUBLE PRECISION T(1),BETA(1) 0372 INTEGER LINE(44) 0373 C0MM0N/NAVC0M/NAVDAY,LINT0T,DEGLIN,IELT0T,DEGELE,SPINRA,IETIMY 9 IET 1IMH,SEMIMA,0ECCEN,0RBINC,PERIGE,ASN0DE,N0PCLN,DECLIN,RASCEN 9 PICLIN 2, PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW COMMON/NAVINI/ 1 EMEGA,AB,ASQ,BSQ.R,RSQ, 2 RDPDG, 3 NUMSEN.TOTLIN.RADLIN, 4 TOTELE,RADELE,PICELE, 5 CPITCH,CYAW,CROLL, 6 PSKEW, 7 RFACT,ROASIN,TMPSCL, 8 B11,B12,B13,B21,B22,B23,B31,B32,B33, 9 GAMMA, GAMDOT, A R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33, B PICTIM,XREF IF(NUMLND.EQ.0) GO TO 91 3374 C RE-INITIALIZE STVEC WITH NEW ORBIT PARAMETERS NAVDAY-NAVDAY-1 CALL STVEC(0.0, MEANOM, X,Y,Z) NAVDAY=NAVDAY+1 CALL STVEC(0.0, MEANOM, X,Y,Z) C C TP SENDS A HOLLERETH FIELD TO THE OUTPUT DEVICE IOUT. C CALL TP(I0UT,132H LANDMARK RESIDUALS FROM UPGORB 0375 * 0376 * 0377 18.050575 PAGE 2 C SEND SECOND HOLLERETH FIELD CALL TP(I0UT,132H N SSYYDDD HHMMSS LCODE LINDIF ELEDIF LANDL0376 *AT LANDLON 0379 * )0380 INAV=1 0362 IYD=MOD(ISYD, 100000) 0383 DO 10 I=1,NUMLND 0384 C C MVCHAR MOVES A HOLLERETH FIELD INTO THE ARRAY LINE. C C FIRST BLANK OUT THE FIELD CALL MVCHARCBLA', 'BLA', LINE, 1,132) 0385 C PLACE INTEGER HOLLERETH FIELD IN ARRAY CALL MVCHAR(I, 'INT', LINE, 3,1) 0386 C PLACE INTEGER HOLLERETH FIELD IN ARRAY CALL MVCHAR(ISYD, 'INT', LINE, 11,1) 0387 TIME=T(I) 0388 JTIME-ITIME(TIME) 0389 C PLACE INTEGER HOLLERETH FIELD IN ARRAY CALL MVCHAR (JTIME, 'INT', LINE, 19,1) 0390 K=LNDPNT(I) 0391 C PLACE INTEGER HOLLERETH FIELD IN ARRAY CALL MVCHAR(LCODE(K) , 'INT', LINE, 27,1) 0392 ILINE=(XLIN(K)-1.0)/8.0 + 1.0 PTIME = T( I )-TMPSCL*FLOAT( I LINE) C TRANSFORM LANDMARK'S LATITUDE AND LONGITUDE TO IMAGE COORDINATES C USING THE NEW ORBIT PARAMETERS. CALL CALRES(IYD,T(I) .BETA (K ) ,MEANOM ,XLAT(K) ,XLON(K) ,YLIN,YELE) C CALCULATE THE DISTANCE BETWEEN THE CALCULATED AND MEASURED C IMAGE POSITIONS. DLIN=XLIN(K)-YLIN 0396 DELE=XELE(K)-YELE 0397 C TRANSFER RESIDUALS INTO A FLOATING HOLLERETH FIELD. CALL FF(8,2,DLIN t LINE,27) 0400 CALL FF(8 f 2, DELE, LINE, 35) 0401 C TRANSFER THE LATITUDE AND LONGITUDE OF LANDMARK INTO AN INTEGER CALL MVCHAR(ILALO(XLAT(K) ), 'INT', LINE, 51,1) 0402 C HOLLERETH FILE. CALL MVCHAR(ILALO(XLON(K) ), 'INT', LINE, 60,1) 0403 C OUTPUT GENERATED HOLLERETH LINE. CALL TP(IOUT,LINE) 0404 10 CONTINUE 0405 RETURN 0406 91 CALL TQMESCNO LANDMARKS FOR SSYYDDD: $', ISYD ) 0407 RETURN 0408 END 0409 SIZE 419 00643 18.050575 PAGE SUBROUTINE CALRES ( IYD.T ,BETA , MEANOM ,XLAT ,XLCN ,XLI N ,XELE) C WRITTEN BY DENNIS PHILLIPS OF SCIENTIFIC PROGRAMMING AND APPLIED C MATHEMATICS, INC. CALRES COMPUTES THE LINE AND ELEMENT POSITION C FOR A GIVEN (XLAT , XLON) USING THE BETA SWEEP ANGLE. C INPUTS: IYD. YEAR DAY C T. TIME OF LANDMARK SCAN LINE. C BETA. SWEEP ANGLE ASSOCIATED WITH LANDMARK SCAN LINE. C MEANOM MEAN ANAMOLY OF SATELLITE POSITION AT EPIC C XLAT. THE LATITUDE OF THE EARTH LANDMARK. C XLON. THE LONGITUDE OF THE EARTH LANDMARK. C OUTPUTS: XLIN. THE COMPUTED LINE POSITION. C XELE. THE COMPUTED ELEMENT POSITION. REAL MEANOM DOUBLE PRECISION T ,XSN , YSN , ZSN ,XSNST , YSNST ,XST , YST ,ANG ,TWPI , BETA 0411 COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0412 1IMH,SEMIMA.,0ECCEN,0RBINC , PERHEL .ASNODE , NOPCLN .DEC LIN ,RASCEN ,PICLI N0413 2, PRERAT.PREDIR, PITCH, YAW, ROLL, SKEW 0414 COMMON/NAVINI/ 0415 1 EMEGA,AB,ASQ,BSQ,R,RSQ, 0416 2 RDPDG, 0417 3 NUMSEN,TOTLIN,RADLIN, 0418 4 TOTELE,RADELE,PICELE, 0419 5 CPITCH.CYAW.CROLL, 0420 6 PSKEW, 0421 7 RFACT,ROASIN,TMPSCL, 0422 8 D11,D12,D13,D21,D22,D23,D31,D32,D33, 0423 9 GAMA,GAMDOT, 0424 A R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33, 0425 B PICTIM,XREF 0426 COMMON/BETCOM/IAJUST,IBTCON,NEGBET,ISEANG 0427 DATA TWPI/6.28318531D0/ 0428 TIME=T 0430 JTIME=ITIME(TIME) 0431 C GET VECTOR TO THE SUN IN INERTIAL COORDINATES CALL SUNVEC(IYD,T,XSN,YSN,ZSN) 0432 C GET VECTOR TO THE SATELLITE IN INERTIAL COORDINATES CALL STVEC( TIME, MEANOM, XS AT, YSAT.ZS AT) C COMPUTE THE X, Y, Z VECTOR POINTING TO THE LANDMARK IN INERTIAL C COORDINATES FROM THE EARTH'S CENTER. JTIME-ITIME(TIME) YLON=RAERAC(IYD,JTIME,XLON)*RDPDG TDIF=TIME-FTIME(JTIME) YLON=YLON+TDIF*TWP 1/24.0 YLAT=XLAT*RDPDG YLAT=GEOLAT(YLAT,l) TANLAT=TAN(YLAT)**2 RR=SQRT( (1.0+TANLAT)/(BSQ+ASQ*TANLAT) )*AB XE=COS(YLON)*COS(YLAT)*RR YE=SIN(YLON)*COS(YLAT)*RR ZE=SIN(YLAT)*RR SENANG=FLALO(ISEANG)*RDPDG 0439 C FETCH CURRENT ATTITUDE STATE AND SETUP SATELLITE CENTERED C COORDINATE SYSTEM WITH Z-AXIS OPPOSITE THE SPIN AXIS. 16.050575 PAGE 2 CALL ATPREC(TIME,DEC ,RAS^ DEC=DEC*RDPDG RAS=RAS*RDPDG SR-SIN(RAS) 0442 CR-COS(RAS) 0443 SD=SIN(DEC) 0444 CD=COS(DEC) 0445 X1=-SR 0446 Y1=CR 0447 Z1=0.0 0448 X2=-SD*CR 0449 Y2=-SD*SR 0450 Z2=CD 0451 X3=CD*CR Y3=CD*SR Z3=SD C FIND LINE NUMBER C0SANG=X3*(XE-XSAT)+Y3*(YE-YSAT)+Z3*(ZE-ZSAT) COSANG=COSANG/SQRT( (XE-XSAT )**2+ ( YE-YSAT)**2+ ( ZE-ZSAT )**2 ) THETA-ATAN2 (R0TM31 ,R0TM33 ) PHI=ASIN(C0SANG/SQRT(R0TM13**2+R0TM33**2) ) YLIN=PHI-THETA XLIN=PICLIN-YLIN/RADLIN C FIND ELEMENT NUMBER C FIND COORDINATES OF THE VECTOR FROM THE SATELLITE TO THE SUN C IN THE SATELLITE COORDINATE SYSTEM XSNST=X1*(-XSN-XSAT)+Y1*(-YSN-YSAT) 0452 YSNST=X2*(-XSN-XSAT)+Y2*(-YSN-YSAT)+Z2*(-ZSN-ZSAT) 0453 C FIND THE COORDINATES OF THE VECTOR FROM THE SATELLITE TO THE EARTH C LANDMARK IN THE SATELLITE COORDINATE SYSTEM. XST=X1* (XE-XSAT )+Y 1* ( YE-YSAT ) YST=X2*(XE-XSAT)+Y2*(YE-YSAT)+Z2*(ZE-ZSAT) C CORRECT FOR ROLL AND SKEW CLIN=COS(YLIN) SLIN=SIN(YLIN) U=R0TM11*CLIN+R0TM13*SLIN V=R0TM21*CLIN+R0TM23*SLIN ANG=DATAN2(YSNST,XSNST)-DATAN2(YST,XST)-BETA+SENANG+ATAN2(V,U) C NORMALIZE THIS ANGLE TO LIE BETWEEN -PI/2.0 AND PI/2.0 ANG-DMOD(ANG,TWPI ) IFUNG.GT. (TWPI/2.0D0) ) ANG=ANG-TWPI IF(ANG.LT. (-TWPI/2.0D0) ) ANG=ANG+TW'PI XELE=ANG/RADELE RETURN 0465 END 0466 SIZE 436 00664 0641 0642 0643 0644 0645 0646 0647 Appendix F Fortran Listings of Standard Geometry Transformation Routines: EARLOC and IMGLOC 18.050575 PAGE 1 SUEROUTINF EARLOC ( ISC , FRTIME, DATF , SLINF ,FLEM , ELAT , ELONG ) C C C PURPOSE: C FARLOC COMPUTES THE GEODETIC LATITUDE (ELAT) AND LONGITUDE C (FLOMG) CORRESPONDING TO THE GIVEN LINF (SLINF) AND ELEMENT C (ELEM) COORDINATES OF THE VISSR IMAGE FRAME STARTING AT FRTIME, C EARLOC PFRFORMS THE INVERSE TRANSFORM OF IMGLOC. C C C PROGRAMMER: LARRY HAM BRICK, NESS /OS I C C ORIGINATION DATE: FEBRUARY 18,1976 C C C ARGUMENTS: C C THE UNITS OF ELAT AND ELONG ARE DECIMAL DEGREES WITH ELAT C POSITIVE NORTH AND FLONG POSITIVE EAST. C C THE UNITS OF SLINF AND ELEM ARE SCAN LINES AND IR SAMPLES C RESPECTIVELY. THEIR FRACTIONAL PARTS SPECIFY HIGHER C RESOLUTIONS. C C NOTE: THE RANGE FOR ELEM IS 0.5 TO 3822.5. THE IR SAMPLE C NUMBER IS OBTAINED BY ROUNDING-OFF ELEM. THE VISIBLE C SAMPLE NUMBER IS OBTAINED BY ROUNDING-OFF ((ELEM-0.5) C *4 + 0.5). C C THE RANGF FOR SLINE IS 0.5 TO 1821.5. THE IR (SCAN) C LINE NUMBER IS OBTAINED BY ROUNDING-OFF SLINE. THE C VISIBLE LINE NUMBER IS OBTAINED BY ROUNDING-OFF C ( (SLINE-0.5)*8 + 0.5). C C ISC DESIGNATES THE SATFLLITE: 1 FOR EAST GOES, 2 FOR WEST. C C FRTIME IS THE GMT IN SECONDS OF THE START OF THE IMAGE FRAME. C C DATE SPECIFIES THE YEAR AND JULIAN DAY AS YYDDD. C C C C REFERENCES: C C THE MATHEMATICAL BASIS FOR THIS ROUTINE IS THE PAPER BY C DENNIS PHILLIPS AND C . T .MOTTERSHEAD . THE REPORT BY WESTING- C HOUSF (CONTRACT NAS 5-23582 .DATED JULY, 1977) IS MORE DIRECTLY C LINKED TO THE FORTRAN CODE AND GIVES DETAILED DEFINITIONS C AND ILLUSTRATIONS OF THE ORBIT/ATTITUDE PARAMETERS AVAILABLE C IN THE VISSR DOCUMENTATION BLOCK. C C EARLOC IS NOT IN THE MOST EFFICIENT COMPUTATIONAL FORM; C EMPASIS HERE IS ON CLARITY OF THE GEOMETRY. 18.050575 PAGE 2 C C DIMENSION PHI(6),S(3,3) COMMON / GRNICH / GRA1 ,GRA2 COMMON / TIM / TIMF1,DATE1,D COMMON / SATATT / SPRA1 , SPRA2, SPDC1 , SPDC2 COMMON / SATPOS / CX ( 11 ) ,C Y ( 11 ) ,CZ ( 1 1) COMMON / ALNANG / ZETA,RHO,ETA COMMON / SPNRAT / SPPER C C GRA1 GREENICH ANGLE IN DEGREES AT TIME TIME1 C GRA2 GREENWICH ANGLE IN DEGREES AT TIMF TIME1 + D C TIME1 GMT IN SECONDS OF EPOCH FOR NAVIGATION PARAMETERS C D PERIOD IN SECONDS OVER WHICH TIME IS NORMALIZED C DATF1 DATE AS YYDDD,YEAR AND JULIAN DAY, OF EPOCH FOR NAV . PAR C SPRA-,SPDC- RIGHT ASCENSION AND DECLINATION OF POSITIVE SPIN C AXIS IN DEGREES AT TIMES: (l)TIMEl AND (2) TIME1 C PLUS D C CX(I ),CY(I) ,CZ(I ),I=1 TO 11 CHEBYCHEV COEFFICIENTS REPRESENT C -ING THE SATELLITE POSITION IN KM BEGINNING AT EPOCH C (TIMED IN EARTH CFNTFRED INERTIAL COORDINATES C NORMALIZED WRT TIME OVER D C RHO IS ROLL OR ELEMFNT BIAS OF THE VISSR IN DFGRFFS C ZETA IS PITCH OR LINE BIAS OF THE VISSR IN DEGREES C ETA IS YAW OR SKEW OF THE VISSR IN DEGREES C SPPFR SPIN PERIOD OF SATFLLITF IN SECONDS C DATA A S B,ACQANG, TOTSMP, SCENCA,SCENCS, TOTSCL / 6378.144,6356.759, 1 9.1875,3822.0,45.0,4096.0,1821.0 / C C A: EQATORIAL RADIUS OF OBLATE ELLIPSOID EARTH IN KM C B: POLAR RADIUS OF OBLATE ILLIPSOID FARTH IN KM C ACQANG: DATA ACQUISTION ANGLE OF VISSR IN DEGREES C TOTSMP: TOTAL IR SAMPLES ACQUIRED IN SCAM LINI C SCENCA: VISSR SCAN ENCODER CHARACTERISTIC ANGLE IN DEGREES Z SCFNCS: NUMBER OF VISSR SCAN ENCODER POSITIONS C TOTSCL: TOTAL NUMBER OF SCAN LINES C PI=3. 1415926535 C C COMPUTE SOME CONSTANT PARAMETERS C E=(A**2 - B**2)/B**2 BSAS= B**2/A**2 AMUF= (2.0*ACQANG/TOTSMP )*PI/180.0 AMUL= (SCENCA/SCENCS)*PI/180.0 CE= (TOTSMP + 1.0)/2.0 CL= (TOTSCL + 1.0)/2.0 C C COMPUTE VISSR ALIGNMENT ANGLES C PITCH = ZETA*PI/180.0 YAW= ETA*PI/180.0 18.050575 PAGE 3 ROLL = RHO*PI/180.0 C C CHECK THAT ELEM AND SLIN? ARE IN LIMITS C IF( (ELFM.GT.3822.5hOR. ( ELEM .LT .0 . 5) ) GO TO 10 IF( (SLINE.GT.1821.5).0R. (SLINE .LT .0 .5) ) GO TO 10 GO TO 15 10 PRINT 30 30 F0RMAT(41H LINE OR ELEMENT NUMBER IS OUT OF LIMITS ) ELAT=999.9 FL0NG=999.9 GO TO 20 15 CONTINUE C C COMPUTE POLAR COORDINATES OF THE UNIT VIEW VECTOR IN C THF VISSR COORDINATE SYSTEM C AZM=AMUE*(ELEM-CE) ELV « AMUL*(CL-SLINE^ C C TRANSFORM THE UNIT VFCTOR TO SATELLITE INFRTIAL COORDINATES C VSl=COS (AZM+ROLL )*COS (PI TCH+ELV )-SI N (AZM+ROLL )*SI N ( YAW )*SI N (PI TCH 1 +ELV) VS2=+SI M( AZM+ROLL ) *COS ( P ITCH+ELV ) +COS ( AZM+ROLL) -SIN ( YAW )*SIN( PITCH 1 +ELV) VS3=-C0S(YAW)*SIN(PITCH+ELV} C C DETERMINE TIME OF THE GIVFN SAMPLE (APPROX. TO THE SCAN C LINE TIME) C T = FRTIME + AINT(SLINE +0.5)*SPPER C C TRANSFORM THF UNIT VIEW VECTOR TO EARTH CENTERED INERTIAL C COORDINATES C CALL SATCRD(T,S ) ITC1= S(l,l)*VSl + S(2,1)*VS2 + S(3,l v:s VS3 VC2 = S(1,2)*VS1 + S(2,2)*VS2 + S(3,2)*VS3 VC3 = S(1,3)*VS1 + S(2,3)*VS2 + S(3,3)*VS3 C C EXTEND THE UNIT VECTOR TO INTERSECTION tflTH THE EARTH C ELLIPSOID C CALL SATVFC(T,PCX,PCY,PCZ) F=(PCX*VC1+PCY*VC2+PCZ*VC3 + F*PCZ*VC3 ) /( 1 .0+E*VC3**2 ) G=(PCX**2+PCY**2+PCZ**2-A**2 + F*PCZ**2) /(l .0 + E*VC3**2) Q = -F - SQRT(F**2-G) C C COMPUTE THE VECTOR FROM THF EARTH CENTER TO THE SAMPLE C POINT IN EARTH CENTERED INFRTIAL COORDINATES C RC1 = PCX + Q*VC1 18.050575 PAGE 4 RC2 - PCY + Q*VC2 RC3 = PCZ + Q*VC3 C C COMPUTE THE EARTH GEODETIC LATITUDE AND LONGITUDE OE SAMPLE C ELAT = ATAN( (1.0/BSAS )*RC3/SQRT(RC1**2+RC2**2)) *180.0/PI C C THE LONGITUDE IS COMPUTED UNDER THF CONVENTION -180 TO +180 C WITH EAST OF GREENWICH + C ANG = ATAN2(RC2,RC1 ) C C COMPUTF THE GREENWICH ANGLE C CALL GRANG(T.W) C ELONG = ANG - W IF(ELONG.LE.-PI ) EL0NG=2 .0*PI+ELONG IF(ELONG.GT.PI ) FL0NG=EL0NG-2 . 0*PI ELONG = ELONG*180.0/PI 20 CONTINUE RETURN END SIZE 499 00763 18.050575 PAGE 1 SUBROUTINE IMGLOC ( ISC ,FRTI ME, DATE , ELAT , ELONG , SLINE , ELEM ) C C C PURPOSE: C IMGLOC COMPUTES THE LINE (SLINE) AND ELEMENT (ELEM) COOR- C DINATFS.IN THE VISSR IMAGE FRAME STARTING AT TIME (FRTIME), C CORRESPONDING TO THE GIVEN EARTH GEODETIC LATITUDE (ELAT) AND C LONGITUDE (ELONG). C IMGLCC PERFORMS THE INVERSE TRANSFORM OF EARLOC. C C C PROGRAMMER: LARRY HAMBRICK, NESS/OSI C C ORIGINATION DATE: FEBRUARY 18,1978 C C C ARGUMFNTS: C C THT UNITS OF FLAT AND ELONG ARE DECIMAL DEGREES WITH ELAT C POSITIVE NORTH AND ELONG POSITIVE EAST. C C THE UNITS OF SLINE AND ELEM ARE SCAN LINES AND IR SAMPLES C RESPECTIVELY. THEIR FRACTIONAL PARTS SPECIFY HIGHER C RESOLUTIONS. C C NOTE: THE RANGE FOR ELEM IS 0.5 TO 3822.5. THE IR SAMPLE C NUMBER IS OBTAINED BY ROUNDING-OFF ELEM. THE VISIBLE C SAMPLE NUMBER IS OBTAINED BY ROUNDING-OFF ((ELEM-0.5) C *4 + 0.5). C C THE RANGE FOR SLINE IS 0.5 TO 1821.5. THE IR (SCAN) C LINE NUMBER IS OBTAINED BY ROUNDING--OFF SLINE. THE C VISIBLE LINE NUMBER IS OBTAINED BY ROUNDING-OFF C ((SLINE-0.5)*8 + 0.5) . C C ISC DESIGNATES THE SATELLITE: 1 FOR EAST GOES, 2 FOR JEST. C C FRTIME IS THE GMT IN SECONDS OF THE START OF THE IMAGE FRAME. C C DATE SPECIFIES THE YEAR AND JULIAN DAY AS YYDDD. C C C REFERENCES: C C THE MATHEMATICAL BASIS FOR THIS ROUTINE IS THE PAPER BY C DENNIS PHILLIPS AND C .T .MOTTERSHEAD . THE REPORT BY tfESTING- C HOUSE (CONTRACT NAS 5-23582 , DATED JULY,19?7) IS MORE DIRECTLY C LINKED TO THE FORTRAN CODE AND GIVES DETAILED DEFINITIONS C AND ILLUSTRATIONS OF THE ORBIT/ATTITUDE PARAMETERS AVAILABLE C IN THE VISSR DOCUMENTATION BLOCK. C C IMGLOC IS NOT IN THE MOST EFFICIENT COMPUTATIONAL FORM; C EMPASIS HERE IS ON CLARITY OF THE GEOMETRY. 18.050575 PAGE 2 C C DIMENSION PHI (6),T(6),S(3,3) COMMON / GRNICH / GRA1,GRA2 COMMON / TIM / TIMEl.DATEl ,D COMMON / SATATT / SPRA1 , SPRA2, SPDC1 , SPDC2 COMMON / SATPOS / CX( 11 ) , C Y ( 11) ,CZ ( 1 1 ) COMMON / ALNANG / ZETA ,RHO , ETA. COMMON / SPNRAT / SPPER C C GRA1 GREENICH ANGLE IN DEGREES AT TIME TIME1 C GRA2 GREENWICH ANGLE IN DEGREES AT TIME TIME1 + D C TIME1 GMT IN SECONDS OF EPOCH FOR NAVIGATION PARAMETERS C D PERIOD IN SECONDS OVER WHICH TIME IS NORMALIZED C DATF1 DATE AS YYDDD.YEAR AND JULIAN DAY, OF EPOCH FOR NAV . PAR C SPRA-.SPDC- RIGHT ASCENSION AND DECLINATION OF POSITIVE SPIN C AXIS IN DEGREES AT TIMES: (l)TIMEl AND (2) TIME1 C PLUS D C CX(I ),CY(I),CZ(I) ,1=1 TO 11 CHEBYCHEV COEFFICIENTS REPRESENT C -ING THE SATELLITE POSITION IN KM BEGINNING AT EPOCH C (TIMED IN EARTH CENTERFD INERTIAL COORDINATLS C NORMALIZED WRT TIME OVER D C RHO IS ROLL OR ELEMENT BIAS OF THE VISSR IN DEGREFS C ZETA IS PITCH OR LINE BIAS OF THE VISSR IN DEGREES C ETA IS YAW OR SKEW OF THE VISSR IN DEGREES C SPPER SPIN PERIOD OF SATELLITE IN SECONDS C DATA STDFC, A, B.AC3ANG, TOTSMP, SCENCA.SCENCS, TOTSCL / 1 6.611,6378.144,6356.759,9.1375,3822.0,45.0,4096.0,1821.0 / C C STDEC:RATIO OF NOMINAL SATELLITE RANGE AND EARTH RADIUS C A: I0ATORIAL RADIUS OF OBLATF ELLIPSOID EARTH IN KM C B: POLAR RADIUS OF OBLATE ELLIPSOID EARTH IN KM DATA ACQUISTION ANGLE OF VISSR IN DEGREES TOTAL IR SAMPLES ACQUIRED IN SCAN LINE VISSR SCAN ENCODER CHARACTERISTIC ANGLE IN DEGREES NUMBER OF VISSR SCAN ENCODER POSITIONS TOTAL NUMBER OF SCAN LINES C ACQAMG C TOTSMP C SCENCA C SCENCS C TOTSCL C PI- 3.1415926535 C C COMPUTE SOME CONSTANT PARAMETERS C E=(A**2 - B**2)/B**2 BSAS= B**2/A*#2 AMUE= (2.0*ACOANG/TCTSMP )*PI/180.0 AMUL= (SCENCA/SCENCS )*PI/180.0 CE= (TOTSMP + 1.0V2.0 CL= (TOTSCL + 1.0)/2.0 C C COMPUTE VISSR ALIGNMENT ANGLES C PITCH = ZETA*PI/180.0 18.050575 PAGE 3 YAtf = ETA*PI/180.0 ROLL = RHO*PI/180.0 C C CONVERT GEODETIC LATITUDE AND LONGITUDE TO RADIANS C PLAMDA= ELAT*PI/180.0 GLCNG= ELONG*PI/180.0 C C COMPUTE GEOCENTRIC LATITUDE C GLAMDA=ATAN(BSAS*TAN(PLAMDA)) C C ESTIMATE THE TIME AT WHICH THE POINT WAS SCANNED C PHI(1)= ATAN( 1.0*SIN(PLAMDA)/(STDEC-COS(PLAMDA)))-PITCH T(1)=FRTIME-(AINT(PHI(1)/AMUL + .5 )-CL )*SPPER C C THIS LOOP ITERATES TO REFINE THE TIME ESTIMATE C DO 40 1=1,5 TI=T(I ) C C COMPUTE THE GREENWICH ANGLE C CALL GRANG(TI ,W) C C DETERMINE THE EARTH POINT VFCTOR IN EARTH CENTERED INERTIAL C COORDINATES C XLl=COS(GLAMDA)*COS(GLONG+W) XL2=C0S(GLAMDA)*SIN(GL0NG+W) XL3=SIN(GLAMDA) R=A/SQRT(1.0 + E*SIN(GLAMDA)**2) RC1=R*XL1 RC2=R*XL2 RC3=R*XL3 C C DETERMINE THE SATELLITE VECTOR C CALL SATVEC(TI ,PX,PY,PZ) C C CHECK WHETHER EARTH POINT IS IN VIEW C XLDP=XL1*PX + XL2*PY + XL3*PZ IF(XLDP.LT.A) GO TO 50 C C COMPUTE VIEW VFCTOR IN EARTH CENTERED INERTIAL COORDINATES C VC1=RC1-PX VC2=RC2-PY VC3=RC3-PZ i C TRANSFORM THE VIEW VECTOR TO SATELLITE INERTIAL COORDINATES 18.050575 PAGE 4 C CALL SATCRD(TI ,S) VS1=S(1,1)*VC1 + S(l,2)*VC2 + S(1,3)*VC3 VS2=S(2,ll*VCl + S(2,2)*VC2 + S(2,3)*VC3 VS3=S(3,1)*VC1 + S(3,2)*VC2 + S(3,3)*VC3 C C DFTFRMINE AZIMUTH AND ELFVATION OF EARTH POINT IN SAT. COORD. C SIGMA=ATAN(VS2/VS1) TANYAW=TAN(YAW) COSYAW=COS(YAW) XI=ATAN(VS3*TANYAW/SQRT(VS1**2+VS2**2-VS3**2*TANYa«/**2)) PHI (1+1 )=ATAN(-VS3/(C0SYAW*C0S(XI )*SQRT( VS1**2+VS2**2) ) ) -PITCH C C CHECK FOR CONVERGENCE OF THE TIME ESTIMATE C DSLINE=(AINT(PKI(I+l)/AMUL + . 5 )-AI NT ( PHI ( I )/AMUL + 0.5)) K = I+1 IF(DSLINE.GE.1.0) GO TO 15 GO TO 25 C C CORRICT THE ESTIMATE OF TIME C 15 T(I+1)=T(I)-DSLINE*SPPER 40 CONTINUE C C COMPUTE THE LINE AND ELEMENT C 25 THETA=XI+SIGMA-ROLL ELEM= CE + THETA/AMUE SLINE=CL-(PHI (K)/AMUL) C C ROUGH CHECKS ON REASONABLENESS OF RESULTS IF(FLEM.LT.0.5) GO TO 35 IF(ELEM.LE.3822.5) GO TO 30 35 PRINT 39 39 FORMAT(17H ELFMENT ERROR ) 30 IF(SLINE.LT.0.5) GO TO 45 IF(SLINE.LE.1821.5) GO TO 60 45 PRINT 49 49 F0RMAT(14H LINE ERROR ) GO TO 60 50 PRINT 55 55 FORMAT(18H BEHIND THE EARTH) ELEM=0.0 SLINF=0.0 60 CONTINUE RETURN END SIZE 556 01066 18.050575 PASE 1 SUBROUTINJ SATCRD(T.S) C COMPUTES THF BASF VECTORS ( THE MATRIX S) OF THE SATELLITE C COORDINATE SYSTEM IN TERMS OF THF EARTH CENTERED INFRTIAL C COORDINATE SYSTEM AT TIME (T) BASED ON THE SATELLITE C POSITION AND ATTITUDE COMMON / SATPOS / CX ( 11) ,C Y ( 1 1) ,CZ ( 11 ) COMMON / SATATT / SPRA1 , SPRA2 , SPDC1 , SPDC2 COMMON / TIM / TIMF1.DATE.D DIMFNSION S(3,3) CALL SPNATT(T,RASC,DECL) S(3,l) = COS(DICL^*COS(RASC) S(3,2)=C0S(DECL)*SIN(RASC) S(3,3)=SIN(DECL) CALL SATVEC(T,PX,PY,PZ) PDOTS=PX*S (3,1) +PY*S ( 3 ,2 )+PZ*S (3 ,3 ) PSQ=PX**2+PY**2+PZ**2 S1DEN=SQRT(PSQ~PD0TS**2) S(1,1)=(-PX+PD0TS*S(3,1) )/SlDEN S(1,2)=(-PY+PD0TS*S(3,2) )/SlDEN S(1,3)=(-PZ+PD0TS*S(3,3))/S1DEN S(2,l)=S(3,2)*S(l,3)-S(3,3)*S(l,2) S(2,2)=S(3,3)*S(1,1)-S(3,1)*S(1,3) S(2,3)=S(3,l)*S(l,2)-S(3,2)*S(l,l) RETURN END SIZE 145 00221 18.050575 PAG-F 1 SUBROUTINE SATVEC(T,PCX,PCY,PCZ) C COMPUTES THF COMPONENTS (PCX,PCY,AND PCZ) OF THF VECTOR TO THE C SATELLITE IN FARTH CENTERED INERTIAL COORDINATES AT C TIMF (T) BASED ON THE CHEBYCHEV COEFFICIENTS CX(I ) , CY ( I ) ,CZ( I ) PIMFNSION BX(13),BY(13),BZ(13) COMMON / SATPOS / CX ( 11 ) ,CY ( 11 ) , C Z ( 1 1 ) COMMON / TIM / TIME1,DATE,D CALL NTIM(T,U) DO 5 1=12,13 BX(I)=0.0 BY(I)=0.0 BZ(I)=0.0 5 CONTINUE DO 15 J=l,ll BX(12-J)=CX(12-J)+2.0*U*BX(l2-J+l)-BX(12-J+2) BY(12-J ^ = CY(12-jU2.0*U*BY(12-J + l)-BY(12-J + 2) BZ ( 12- J )=CZ ( 12-J ) +2 .0*U*BZ ( 12-J+l )-BZ( 12-J+2 ) 15 CONTINUE PCX=(BX(1)-BX(3) )/2.0 PCY=(BY(l)-BY(3))/2.0 PCZ=(BZ(l)-BZ(3))/2.0 RETURN END SIZE 156 00234 18.050575 PAGE 1 SUBROUTINF SPNATT (T , SPRASC , SPDECL) C COMPUTES THE RIGHT AS CENS ION ( SPRASC ) AND DECLINATION ( SPDFCD C OF THE SPIN AXIS AT TIME (T) WITH A LINEAR C INTERPOLATION BASED ON VALUFS OF SPRA1.2 AND SPDC1.2 C C THE INTERPOLATION TECHNIQUE USED HERE IS A SIMPLE APPROXIM- C ATION BUT IS SUFFICIENTLY" PRECISE FOR THF SMALL PRECESSION C RATES ENCOUNTERED. C C SPRA1,2 HAS RANGE -180 TO +180 DEGREES C SPDC1,2 HAS RANGE -90 TO +90 DEGREES C SPRASC HAS RANGE -PI TO +PI RADIANS C SPDFCL HAS RANGE -PI/2 TO +PI/2 RADIANS C NOTE IF THE ABSOLUTE VALUE OF THE DIFFERENCE BETWEEN C SPRA1 AND SPRA2 IS GREATER THAN 90 DEGREES, THE C INTERPOLATION OF DECLINATION WILL BE INCORRECT COMMON/ SATATT / SPRA1 ,SPRA2,SPDC1 ,SPDC2 COMMON / TIM / TIME1,DATE,D PI=3. 1415926535 CALL NTIM(T,U) RASC1=SPRA1*PI/180.0 RASC2=SPRA2*PI/180.0 DECL1=SPDC1*PI/180.0 DFCL2=SPDC2*PI/180.0 IF(RASC1.LT.0.0) RASC1=2.0*PI+RASC1 IF (RASC2.LT. 0.0) RASC2=2.0*PI+RASC2 SPRASC=0.5*(RASC2+RASC1+(RASC2-RASC1 )*U) SPRASC=AMOD(SPRASC,2.0*PI) IF(SPRASC .GT.PI ) SPRASC=SPRASC-2.0*PI SPDECL=0.5*(DECL2+DECL1+(DFCL2-DECL1)*U) RETURN END SIZE 101 00145 18.050575 PAGE SUBROUTINE NTIM(T.U) C COMPUTFS NORMALIZED TIME (U) FROM ACTUAL TIME ( T 1 BASED C ON TIME1 AND D. C TIME1 IS IN SECONDS C D IS 13 HOURS IN SFCONDS COMMON / TIM / TIME1,DATF,D IF (TIMF1 .GT.T) GO TO 5 U=2.0*(T-TIME1 )/D-l .0 GO TO 10 5 U=2.0*(T + 24. 0*60. 0*60. 0-TIMED/D-1.0 10 RETURN END SIZE 33 00041 18.050575 PAGE SUBROUTINE GRANG(T.W) COMPUTE GREENWICH ANGLEU) IN RADIANS FROM TIME(T) BASED ON VALUES OF GRA1 AND GRA2 ( RANGF IS -180 TO W HAS RANGE -PI TO +PI IN RADIANS COMMON / GRNICH / GRA1 ,GRA2 COMMON / TIM / TIME1,DATE,D PI=3. 1415926535 CALL NTIM(T,0) W1=GRA1*PI/180.0 H/2=GRA2*PI/180.0 IF(W1.LT.0.0) W1=2.0*PI + Wl IFU2.LT. 0.0) W2=2.0*PI + W2 IF(W2.LT.tfl) tf2=W2+2.0#PI W=0 . 5* ( W2+W 1+ ( W2-W1 ) *U ) W=AMOD(W,2.0*PI) IF(W.GT.PI) W=W-2.0*PI RETURN END IN DEGREES) +180 DEGREES SIZE 84 00124 •U.S. GOVERNMENT PRINTING OFFICE : I960 0-311-01,6/223 (Continued from inside front cover) NESS 87 NESS 88 NESS 89 NESS 90 NESS 91 NESS 92 NESS 93 NESS 94 NESS 95 NESS 96 NESS 97 NESS 98 NESS 99 NESS NESS NESS NESS NESS NESS NESS NESS NESS NESS NESS Jenifer H. Wartha, August 1977, 68 pp. (PB-276-386/AS) Standby Satellites. Bruce Sharts and Chris Dunker, September Catherine M. Frain (Compiler), Pope, September 1977, Atlantic Tropical and Subtropical Cyclone Classifications for 1976. D. C. Gaby, J. B. Lushine, B. M. Mayfield, S. C. Pearce, K.O. Poteat, and F. E. Torres, April 1977, 13 pp. (PB-269-674/AS) National Environmental Satellite Service Catalog of Products. Dennis C. Dismachek (Editor), June 1977, 102 pp. (PB-271-315/AS) A Laser Method of Observing Surface Pressure and Pressure-Altitude and Temperature Profiles of the Troposphere From Satellites. William L. Smith and C. M. R. Piatt, July 1977, 38 pp. (PB- 272-660/AS) Lake Erie Ice: Winter 1975-76. In-Orbit Storage of NOAA-NESS 1977, 3 pp. (PB-283-078/AS) Publications and Final Reports on Contracts and Grants, 1976. August 1977, 11 pp. (PB-273-169/AS) Computations of Solar Insolation at Boulder, Colorado. Joseph H. 13 pp. (PB-273-679/AS) A Report on the Chesapeake Bay Region Nowcasting Experiment. Roderick A. Scofield and Carl E. Weiss, December 1977, 52 pp. (PB-277-102/AS) The TIROS-N/NOAA A-G Satellite Series. Arthur Schwalb, March 1978, 75 pp. (PB-283-859/AS) Satellite Data Set for Solar Incoming Radiation Studies. J. Dan Tarpley, Stanley R. Schneider, J. Emmett Bragg, and Marshall P. Waters, III, May 1978, 36 pp. (PB-284-740/AS) Publications and Final Reports on Contracts and Grants, 1977. Catherine M. Frain (Compiler), August 1978, 13 pp. (PB-287-855/AS) Quantitative Measurements of Sea Surface Temperature at Several Locations Using the N0AA-3 Very High Resolution Radiometer. Laurence Breaker, Jack Klein, and Michael Pitts, September 1978, 28 pp. (PB-288-488/AS) An Empirical Model for Atmospheric Transmittance Functions and Its Application to the NIMBUS-6 HIRS Experiment. P.G. Abel and W.L. Smith, NESS, and A. Arking, NASA, September 1978, 29 pp. (PB-288-487/AS) 100 Characteristics and Environmental Properties of Satellite-Observed Cloud Rows. Samuel K. Beckman (in consultation). 101 A Comparison of Satellite Observed Middle Cloud Motion With GATE Rawinsonde Data. D. Herman, January 1979, 13 pp. (PB-292-341/AS) 102 Computer Tracking of Temperature-Selected Cloud Patterns. Lester F. 15 pp. (PB-292-159/AS) 103 Objective Use of Satellite Data To Forecast Changes in Intensity of Tropical Disturbances. Carl 0. Erickson, April 1979, 44 pp. (PB-298-915) 104 Publications and Final Reports on Contracts and Grants. Catherine M. September 1979. (PB80 122385) 105 Optical Measurements of Crude Oil Samples Under Simulated Conditions. John S. Knoll, October 1979, 20 pp. (PB80 120603) 106 An Improved Model for the Calculation of Longwave Flux at 1 1 urn. P. G. October 1979, 24 pp. (PB80 119431) 107 Data Extraction and Calibration of TIROS-N/NOAA Radiometers. Levin Lauritson, Gary J. and Frank W. Porto, November 1979. (PB80 150824) 108 Publications and Final Reports on Contracts and Grants. Catherine M. Frain, (Compiler), Aug- ust 1980. 109 Catalog of Products, Third Edition. Dennis C. Dismachek, Arthur L. Booth, and John A. Leese, April 1980, 130 pp. 110 GOES Data Collection Program. Merle Nelson, August 1980. Leroy January 1979, Frain, (Compiler), Warren A. Hovis and Abel and A. Gruber, Nelson, PENN STATE UNIVERSITY LIBRARIES AD0DQ72D4117fi NOAA SCIENTIFIC AND TECHNICAL PUBLICATIONS The National Oceanic and Atmospheric Administration was established as part of the Department of Commerce on October 3, 1970. The mission responsibilities of NOAA are to assess the socioeconomic impact of natural and technological changes in the environment and to monitor and predict the state of the solid Earth, the oceans and their living resources, the atmosphere, and the space environment of the Earth. The major components of NOAA regularly produce various types of scientific and technical informa- tion in the following kinds of publications: PROFESSIONAL PAPERS — Important definitive research results, major techniques, and special inves- tigations. CONTRACT AND GRANT REPORTS — Reports prepared by contractors or grantees under NOAA sponsorship. ATLAS — Presentation of analyzed data generally in the form of maps showing distribution of rainfall, chemical and physical conditions of oceans and at- mosphere, distribution of fishes and marine mam- mals, ionospheric conditions, etc. TECHNICAL SERVICE PUBLICATIONS — Re- ports containing data, observations, instructions, etc. A partial listing includes data serials; prediction and outlook periodicals; technical manuals, training pa- pers, planning reports, and information serials; and miscellaneous technical publications. TECHNICAL REPORTS — Journal quality with extensive details, mathematical developments, or data listings. TECHNICAL MEMORANDUMS — Reports of preliminary, partial, or negative research or technol- ogy results, interim instructions, and the like. Information on avaitability of NOAA publications can be obtained from: ENVIRONMENTAL SCIENCE INFORMATION CENTER (D822) ENVIRONMENTAL DATA AND INFORMATION SERVICE NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION U.S. DEPARTMENT OF COMMERCE 6009 Executive Boulevard Rockville, MD 20852 N0AA--S/T 80-183