AREA Format Information Copyright 1996 Space Science and Engineering Center (SSEC) University of Wisconsin - Madison All Rights Reserved Permission is granted to make and distribute verbatim copies of this document, provided the copyright notice and this permission are preserved on all copies. SSEC makes no warranty of any kind with regard to the software or accompanying documentation, including but not limited to the implied warranties of merchantability and fitness for a particular purpose. SSEC does not indemnify any infringement of copyright, patent, or trademark through use or modification of this software. Mention of any commercial company or product in this document does not constitute an endorsement by SSEC. Many of the designations used by manufacturers and sellers to distinguish their products are claimed as trademarks. Where those designations appear in this document, and SSEC was aware of the trademark claim, the designations are printed in caps or initial caps. The information in this document is subject to change without notice. Considerable effort has been expended to make this document accurate and complete, but SSEC cannot assume responsibility for inaccuracies, omissions, manufacturers' claims or their representations. University of Wisconsin - Space Science and Engineering Center March, 1993 Satellite imagery within McIDAS occurs in several different forms depending on the original data source. Regardless of the source, an image is basically a sequence of 'lines' arranged one below the previous one, numbered from top to bottom with the top line being number 1. Each line consists of a sequence of 'elements' arranged across the line and numbered left to right, the leftmost element being number 1. An element may be an 8, 16 or 32 bit quantity, depending on the image source and format. This line/element numbering scheme determines a pair of coordinates for each element, called the 'image coordinates' of the element. This coordinate system is defined only by the satellite/camera combination and is independent of how the data is stored. ********************************************************************* AREA FILE DESCRIPTION ********************************************************************* Satellite imagery is stored on disk in data structures called digital areas or just 'areas.' Each area is stored as a binary file containing all the information necessary to display and navigate the image. These files are named AREAxxxx, where xxxx is a number between 0000 and 9999. All four digits must be present. This number is referred to as the area file number. Each file has the following format: Block Size Area Directory entry (ID) 256 bytes, portions SSS specific Navigation entry variable, SSS specific Calibration entry variable, SSS specific Digital data variable, SSS specific Comment records (audit trail) variable, 80 character records (ASCII) The first 256 bytes of every file contain the Area Directory entry for the image. Each value in the entry is a 4-byte word (64 words). The description of each entry is described below. The next entry contains the NAV block for the image. Again, each value in this block is a 4-byte quantity holding the navigation parameters necessary to earth locate the image pixels. If no navigation is available for the image, the block is zero filled. In addition, word 35 of the area directory will contain the byte offset of the NAV block in the file. If necessary, the next block of bytes may contain a CAL entry. This entry is present if the data must be calibrated before it can be displayed. Like the Area Directory and NAV blocks, each value in a CAL block is a 4-byte quantity. If present, word 63 of the Area Directory entry is a byte offset to the CAL block's position within the file. If word 63 is zero, the CAL block is not contained in the file. Next comes the digital data values for the image (line by line, no record separators). The length of the DATA block is computed using values within the area directory entry. Word 34 of the Area Directory is the byte offset to the start of the DATA block in the file. Finally, a number of comment records may be appended (see W64 of the Area Directory). These records are usually used to keep an audit trail of modifications made to this particular data file. Each comment record is always 80 ASCII characters. AREA DIRECTORY Each area extracted from a satellite image has an associated record called the 'Area Directory.' This directory contains information about the area and the image it was created from. The data in the directory is stored as 32-bit (4-byte) two's complement binary integers or as ASCII character data. The length of the directory is 256 bytes (64 Words). The directory contents are: WORD NAME AND DESCRIPTION ---- ---------------------- W1. Contains zeros if the record is valid W2. Area Format: always 4 (as of June, 1985). W3. SSS: Satellite Sensor Source; see Table 1 W4. YYDDD: Nominal Year and Julian day of the data W5. HHMMSS: Nominal Time of the data W6. UPPERLEFTLINE: Image line coordinate of area line 0, elem. 0 W7. UPPERLEFTELE: Image element coordinate of area line 0, elem.0 W8. Not used. W9. NLINES: Number of lines in this digital area W10. NELES: Number of elements in each line W11. ELESIZ: Number of bytes/element (1, 2 or 4) W12. LINERES: Line Resolution; spacing in image-lines between consecutive area lines W13. ELERES: Element Resolution; spacing in image-elements between consecutive area elements W14. NCHANS: Maximum number of bands/line of the area W15. PRESIZ: Length of the line prefix in bytes; indicated by the sum of W49, W50, W51 (+ 4 if the validity code is present; see W36) W16. PROJ: McIDAS user project number under which the area was created W17. CREATION DATE: Area creation day in YYDDD format W18. CREATION TIME: Area creation time in HHMMSS format W19. FILTER MAP: (for Multi-Channel Images) a 32 bit vector If a bit=1, there is data for that band in the area; the LSB (rightmost) is for band 1 W20-24. SSS specific W25-32. Memo: 32 ASCII characters available for comments W33. Area file number (xxxx from the file name) W34. Byte offset to the start of the data within an area file W35. Byte offset to the start of the navigation block within an area file W36. Validity Code: If these bytes are nonzero, they represent a code which must match the first 4 bytes of the line prefix; if the code and prefix bytes are not equal, the line does not contain valid data and applications must treat as missing. W37-48. SSS specific W49. Line prefix documentation section length in bytes W50. Line prefix calibration section length in bytes W51. Line prefix level map section length in bytes W52. Image source type: 'VISR', 'VAS ', 'AAA ', 'ERBE', 'AVHR', ... W53. Calibration Type: units in which the digital data is stored (e.g. 'RAW ','TEMP','BRIT',...) W54-62. Internal use only W63. Byte offset to the start of the calibration block within the area file W64. Number of card image comments (80 bytes each) appended to the end of the digital data. SSS Code Sensor Source 0 Non-Image Derived Data 2 Graphics 3 MDR Radar 4 METEOSAT Visible (obsolete) 5 METEOSAT Infrared (obsolete) 6 METEOSAT Water Vapor (obsolete) 7 RADAR 8 Miscellaneous Aircraft Data (MAMS) 9 Raw Meteosat 12 GMS Visible 13 GMS Infrared 14 ATS 6 Visible 15 ATS 6 Infrared 16 SMS-1 Visible 17 SMS-1 Infrared 18 SMS-2 Visible 19 SMS-2 Infrared 20 GOES-1 Visible 21 GOES-1 Infrared 22 GOES-2 Visible 23 GOES-2 Infrared 24 GOES-3 Visible 25 GOES-3 Infrared 26 GOES-4 Visible (VAS) 27 GOES-4 Infrared & Water Vapor (VAS) 28 GOES-5 Visible 29 GOES-5 Infrared & Water Vapor (VAS) 30 GOES-6 Visible 31 GOES-6 Infrared 32 GOES-7 Visible 33 GOES-7 Infrared 36-40 NOAA Series Satellites (1-5) 41 TIROS-N 42 NOAA-6 43 NOAA-7 44 NOAA-8 45 NOAA-9 46-49 MARINER X Spacecraft 50 Hubble Space Telescope 54 Meteosat-3 55 Meteosat-4 56 Meteosat-5 57 Meteosat-6 60 NOAA-10 61 NOAA-11 62 NOAA-12 63 NOAA-13 64 NOAA-14 70 GOES-I (Imager) 71 GOES-I (Sounder) 70 GOES-J (Imager) 71 GOES-J (Sounder) 70 GOES-K (Imager) 71 GOES-K (Sounder) 70 GOES-L (Imager) 71 GOES-L (Sounder) 70 GOES-M (Imager) 71 GOES-M (Sounder) 80 ERBE 87 DMSP F-8 88 DMSP F-9 89 DMSP F-10 90 DMSP F-11 91 DMSP F-12 95 FY-1b 96 FY-1c 97 FY-1d Satellite Sensor Source numbers TABLE 1. DIGITAL DATA The data stored in an area is arranged in a two-dimensional array of lines and elements, like the image from which it was created. Each element in an area has a line number, starting with zero at the top, and an element number starting with zero on the left side of the line. This line/element number pair defines a coordinate system for the elements in the area, called the 'area coordinates.' *** Note this is a 'geographical' description of the data, that is, if a line of data is missing the corresponding place in the data file MUST be either filled with zeros or flagged using (non-matching) Val Codes. If all the elements of a satellite image were contained in an area, there would be no point in distinguishing between image and area coordinates. This is not usually the case however, since images are usually too large to process efficiently in toto. For example, a GOES VISSR image in the visible-light band contains 14568 lines with 15288 elements per line. This image requires over 200 megabytes of storage. What is actually stored in an area is a rectangular subset of an image which is produced from the image by sampling (in the vertical direction) and averaging (in the horizontal direction). For example, a subset may be a geographical subset of an image, as in a real-time VIS area of the USA made from a global image. In the case of VAS images, it is also possible to only include in the disk area a subset of the measured spectral bands, so that each element contains fewer data values than are contained in the original image. In order to map an area back to the original image the following formulas are used: Image_Line = UpperLeftLine + (Area_Line * LineRes) Image_Element = UpperLeftEle + (Area_Element * EleRes) UpperLeftLine is the image line coordinate of the first area line. UpperLeftEle is the image element coordinate of the first area element. When LineRes = EleRes = 1, the area is said to be at 'Resolution 1', or 'Full Visible Resolution.' When, for example, LineRes = EleRes = 4, every fourth line and element of an image originally at resolution 1 are included in the area. This area is said to be at 'Resolution 4'. An area may be viewed as a continuous stream of bytes numbered from 0. Within this stream of bytes, the area data is contained line-by-line, with the lines in order, first to last. Each line is further divided into two parts, the 'Line Prefix' and the actual line data (image elements) See Figure 1. The line prefix contains documentation about the image and the particular line. line prefix 1 line data 1 line prefix 2 line data 2 etc. |______________|_____________|_______________|_____________|___ ... 0 byte numbers increase --> AREA DATA LAYOUT Fig. 1 The size and content of the line prefix depend heavily upon the area type, which in turn is determined by the data source. The area type is given in word W52 of the Area Directory. Regardless of the area type, each line in the area has the same length prefix. This length, in bytes, is given in word W15 of the Area Directory. It may be zero; that is, there may be no line prefix for the area. The line prefix may contain a validity code in the first four bytes of the prefix. If this code is present, a nonzero and non-null value is contained in W36 of the Area Directory. The validity code at the start of the line prefix must match exactly the value in W36 of the directory. If not, the line is invalid and should be ignored. The remaining portion of the prefix may contain up to 3 regions. The first region, DOC, is reserved for information specific to a particular satellite line (such as the VISSR IR documentation). The length (in bytes) of this region is given in word W49 of the directory. The next region, CAL, holds the calibration coefficients for the data. These coefficients repeat one time for each band/channel in the image. For example, VAS images normally contain 13 bands per line, each of which has 8 bytes of calibration. An additional 12 bytes of calibration precede the repetition groups in the calibration section yielding a total of 116 bytes in a VAS calibration region. The length of the calibration region is given in directory word W50. The last region is the level map, LEV. For each band/channel on a particular line, there is a one-byte entry in the level map region containing the number of that band/channel. These are only meaningful if they are greater than zero. For an image containing a single band, the level region is optional; otherwise, it must be present. The length of this region is contained in W51 of the Area Directory. The lengths of these regions, when present, are always multiples of four. val code documentation calibration level |______________|_____________|_______________|_____________| 0 byte numbers increase --> AREA LINE PREFIX EXAMPLE Fig. 2 Each line in an area has the same total length. This length, in bytes, is also always a multiple of four. ********************************************************************* ADDITIONAL STRUCTURE OF METEOSAT PDUS AREAS (SSS = 54 thru 57) ********************************************************************* METEOSAT PDUS images have been remapped and calibrated at the ground station before the stretched signal is disseminated. Thus the navigation and calibration data sections are simpler. This data is always 8 bits and the different bands are stored in separate areas. All Labels and Headers are documented in "Meteosat High Resolution Image Dissemination", 1992, from EUMETSAT. PDUS AREA DIRECTORY W14 = 1 (separate bands are always stored separately) W19. Filter map values are: 0 for visible image 128 (8th bit from right) for IR 512 (10th bit from right) for WV W22. MIEC absolute calibration band value (IR or WV as appropriate) from calibration section of header Stored as scaled integer (XXXXX, value is .XXXXX) W23. Space count corresponding to calibration value from calibration section of header (see DOC section below) Stored as scaled integer (XXX, value is XX.X) W24. The physical sensor number from the header (1 or 2). W37. Line offset of the southeast corner of the area in image coordinates. This is a 16 bit value picked up from the header, right justified, and has 1 added to it. W38. Element offset of the southeast corner of the area in image coordinates. This is a 16 bit value picked up from the header, right justified, and has 1 added to it. W39. Satellite Center Longitude (of rectification). This is a 16 bit value picked up from the header and right justified. W44. zero W49. Line prefix documentation length: 24 W50. Line prefix calibration length: 0 (zero) W51. Line prefix level map: 0 (zero) W52. Image source type: 'MSAT' W53. Calibration Type: 'RAW ' W54. zero if data was ingested as sent (full resolution) one if data was sampled down (every other pixel) W55. bitmap indicating what data was sent in image this data was extracted from. Bits number right to left (lsb to msb). bit 0: 1 if visible was included in transmission, 0 if not. bit 1: 1 if ir was included in transmission, 0 if not. bit 2: 1 if water vapor included in transmission, 0 if not. all other bits = 0. PDUS AREA DATA Prefix description: Val Code is optional, but recommended to flag missing data. Missing data can not be simply omitted, but must have placeholder data (either all zeros or indicated by Val Code that does not match value in Area Directory). Line Doc: 24 bytes long. This is a copy of the label which arrives with every subframe. Line Cal: 0 bytes long. Not used Level: 0 bytes long. Not used. Each band of Meteosat is stored in a separate area. Data: Each value is 8 bits (as transmitted) and is stored West-to-East (the inverse of how it is transmitted) and North-to-South in the area (again the inverse of how it is transmitted). PDUS NAVIGATION The format of the PDUS navigation block is: W1. Navigation type: 'MSAT' W2. YYDDD - Year and Julian day of this navigation W3. HHMMSS - hour, minutes, seconds UT of this navigation W4. reference position for the telescope (0 for PDUS) W5. line number corresponding to the telescpe reference position (0 for PDUS ) W6. Center scan line (1250 for PDUS) W7. DDMMSS Center Longitude of Rectification (West positive). W8. unused (0) W9. unused (0) W10. YYDDD - Year and Julian day of this navigation W11-256. zeros, unused. PDUS CALIBRATION No Calibration file section is needed for PDUS. ********************************************************************* ADDITIONAL STRUCTURE OF NOAA AVHRR AREAS (SSS = 60 thru 64) ********************************************************************* AVHRR data is data captured during flyover. Similar formats are used for the stored data (GAC and LAC). See "NOAA Technical Memorandum NESS 107", 1988 for a complete description of the line Doc fields. AVHRR AREA DIRECTORY W14 = 5 (all 5 bands are normally stored in one file) W19. Filter map value is: 31 (all 5 bands). W49. Line prefix documentation length: 192 W50. Line prefix calibration length: 40 W51. Line prefix level map: 8 W52. Image source type: 'TIRU' W53. Calibration Type: 'RAW ' AVHRR AREA DATA This data (including doc, telemetry, etc) is transmitted as 10 bits (we store in 16 bits) and 5 bands (all 5 bands are stored interleaved in one area). The lines of data are stored in the order received. The 10 bits are mapped to 16 as follows: 0dddddddddd00000 where d is a bit of data (that is, the 10 bit value is shifted left 5 bits and zero filled). Prefix description: Val Code is optional, but recommended to flag missing data. Missing data can not be simply omitted, but must have placeholder data (either all zeros or indicated by Val Code that does not match value in Area Directory). Line Doc: 192 bytes. This is the Doc section from the signal transmission. Line Cal: 40 bytes of zeros. This is filled in during post-processing within McIDAS. Level: 8 bytes. Values 1 through 5 (left to right) with three pad zeros in successive bytes. 1 through 5 indicate the order the bands appear in the subsequent data section. These bands are numbered the same as in NESS 107, for example, 1 and 2 are the visible channels. Data: The 16 bit values are stored one 16 bit value from each channel (those that cover the same geography), then repeating, until all 2048 samples are included for each channel. For example, 1,2,3,4,5,1,2,3,4,5,1,2,... Data is stored approximately west to east as sensed and transmitted by the spacecraft. Note that the data pixels are transmitted interleaved in this manner. AVHRR NAVIGATION These parameters are normally derived from TBUS messages from GTS or another conventional data circuit. The format of the AVHRR navigation block is: W1. 'TIRO' W2. SSYYDDD SSS, year, julian day of this navigation W3. HHMMSS hours, minutes, seconds UT of this navigation W4. orbit type (always 1) W5. YYMMDD Epoch date (ETIMY) W6. HHMMSS Epoch time (ETIMH) W7. semimajor axis (SEMIMA), km*100 W8. orbital eccentricity (ECCEN), *1000000 W9. orbital inclination (ORBINC), degrees*1000 W10. mean anomaly (MEANA), degrees*1000 W11. argument of perigee (PERIGEE), degrees*1000 W12. right ascension of ascending node (ASNODE), degrees*1000 W13. number of samples per line (2048) W14. angular increment between samples, degrees*1000000 (left-to-right is positive) W15. fraction of second in epoch time W16-45. (reserved) W46. direction of satellite pass (-1=descending, +1=ascending) W47. number of the earliest line to navigate (image coordinates) W48. time at the start of the earliest line, milliseconds from start of day W49. time interval between lines, milliseconds (see word 53) W50. flag for inverted image (0=normal, 1-flipped) W51. number of lines in the area (for flipped) W52. number of elements in the are (for flipped) W53. same as 40 except in microseconds (preferred over 49) W54. time between individual pixels *100000000 W55-120. (reserved) W121-128. memo entry; up to 32 characters of comments AVHRR CALIBRATION No cal section is needed for AVHRR, since the data is calibrated line-by- line in the post-processing. ********************************************************************* ADDITIONAL STRUCTURE OF NOAA TIP AREAS (TOVS) (SSS = 60 thru 64) ********************************************************************* TIP data is extracted from AVHRR, GAC, or LAC. See "NOAA Technical Memorandum NESS 107", 1988 for a complete description of the line Doc fields, including TIP. TIP contains the TOVS (MSU and HIRS) instrument data as well as SSU and other infrequently sensed data. TIP AREA DIRECTORY W14 = 1 W19. Filter map value is: 1 W49. Line prefix documentation length: 196 W50. Line prefix calibration length: 0 (zero) W51. Line prefix level map: 0 (zero) W52. Image source type: 'TIRU' W53. Calibration Type: 'RAW ' TIP AREA DATA This data (including doc, telemetry, etc) is transmitted as 10 bits (we store in 16 bits). The 10 bits are mapped to 16 as follows: 0dddddddddd00000 where d is a bit of data (that is, the 10 bit value is shifted left 5 bits and zero filled). Prefix description: Val Code is optional, but recommended to flag missing data. Missing data can not be simply omitted, but must have placeholder data (either all zeros or indicated by Val Code that does not match value in Area Directory). Line Doc: 196 bytes. This is the Doc section from the signal transmission (see NESS 107 for description; 192 bytes from the signal starting with the Satellite ID field), plus 4 bytes (parity information, now unused). Line Cal: 0 bytes. Not needed for TIP area. Level: 0 bytes. Not needed for TIP area. Data: The TIP portion of each sector is stored as it is received (with 10 into 16 bit shift done as described above). TIP NAVIGATION The navigation section of TIP files is filled with zeros. TIP CALIBRATION No calibration section is needed for TIP files. ********************************************************************* ADDITIONAL STRUCTURE OF GOES VISSR AREAS (SSS = 30 thru 33) ********************************************************************* The following description is for type VISR (W52 in the Area Directory). GOES AREA DIRECTORY The following directory words are specific to GOES : W37-44. PDL: Exists if the image was made in mode AA or AAA packed byte format (Program Data Load from signal Doc) W45. If image is a Mode AA DS ('STEP','DWEL') these bytes indicate the origin of Band 8; not used for mode AAA. W46. Actual image start YYDDD W47. Actual image start HHMMSS W48. Actual starting scan GOES GENERAL The GOES VISSR instrument produces images in two spectral bands, dubbed VISIBLE and INFRARED (IR). A particular area contains only one of these bands. Which band is in an area may be determined from the satellite code (SSS) in W3 of the Area Directory. Generally, even numbered SSS represents visible areas; odd numbered SSS represents IR areas. Every element in a VISSR area contains one 8-bit pixel, representing raw data from the instrument. If the area contains IR data, the observed temperature may be inferred from the pixel value, according to the following formulae: T = 418 - B (B>176 OR B=176) T = 330 - (B / 2) (B<176 OR B=176) Where: T is the brightness temperature (degrees K). B is the pixel value (0 to 255). Note that for IR data, the highest pixel values correspond to the coldest temperatures (space is white). The line prefix in a VISSR area may be absent entirely, or may contain just the 4-byte validity code. For IR areas, the line prefix may include 128 bytes of IR documentation transmitted by the satellite with the data. Note that most of the useful information in the IR documentation is processed by SSEC into a more useful form in the navigation records, which are included on the tape. The highest resolution (i.e. lowest values of LINERES and ELERES in the Area Directory) possible for a visible area is 1; for an IR area it's 4, as longer wavelengths inherently have less resolution. For a GOES satellite, resolution 1 means approximately 1 km resolution at the satellite subpoint. ********************************************************************* ADDITIONAL STRUCTURE OF GOES VAS AREAS ********************************************************************* There are 2 types of VAS data: Mode AA and Mode AAA. Most satellite data after 24 March 1987 (Julian day 87083) will be AAA. In the Area Directory, W52 contains the Source Type which will contain: VISR - for Visible imagery (1 byte data) VAS - for VAS multibanded imagery (2 byte data) AAA - for VAS and IR imagery after 24 March 1987 (2 byte data) The VAS produces observation data for a given spatial location in up to 12 different IR spectral bands and a visible band. All or some of the IR bands may be included in a single VAS type area. The visible, however, must be contained in a separate area of VISSR type. As a result, it may require two areas to contain the total information transmitted by the satellite during a given time period. The characteristics of the VAS bands are given in Table 2 below. ----.------------------.----------.---------.-------- BAND| SPECTRAL FILTERS|PURPOSE |MAIN |OTHER NUM.| CENTER WIDTH| FOR |ABSORBING|SIGNIFICANT | UM 1/CM 1/CM|SOUNDING |GAS |EFFECTS ----.------ ------ ----.----------.---------.-------- 1 |14.73 678.7 10.|TEMP |CO2 |O3 2 |14.48 690.6 16.|TEMP |CO2 |O3 3 |14.25 701.6 16.|TEMP |CO2 |O3 4 |14.01 713.6 20.|TEMP |CO2 |O3 5 |13.33 750. 20.|TEMP |CO2,H2O |O3 6 | 4.525 2210. 45.|TEMP+CLOUD|N20 |SUN 7 |12.66 790. 20.|MOISTURE |H2O |CO2+DUST 8 |11.17 895. 140.|SURFACE |H2O |CO2+DUST 9 | 7.261 1377.2 40.|MOISTURE |H2O |CO2 10 | 6.725 1487. 150.|MOISTURE |H2O |-- 11 | 4.444 2250. 40.|TEMP+CLOUD|N2O,CO2 |SUN 12 | 3.945 2535. 140.|SURFACE |H2O |SUN+DUST ----|------ ------ ----|----------|---------|----------- TABLE 2. The number of spectral bands present in a VAS area may be determined from W14 of the Area Directory. The filter map (W19 of the Area Directory) describes the particular bands associated with an area. A bit is set to 1 in this word for each band appearing in the area. The number of bands will agree with the value implied by NCHANS (W14 of the Area Directory). The line prefix is always in a VAS type area. It is required in order to calibrate the raw values into temperatures. This cannot be done with fixed formulae as is done for VISSR IR. If W52 is AAA, most of the calibration information is not in the line prefix, but in a separate block (AUX) on the tape. The line prefix consists of: 4-byte validity code 512 bytes of IR common documentation 116 bytes of VAS calibration information organized as follows: Day (YYDDD) 4-byte binary integer Time of scan (HHMMSS) 4-byte binary integer Scan number (satellite coordinate line number) 4-byte binary 13 8-byte groups (1 per possible band), each containing: Channel number 2-byte binary integer Number of spins 2-byte binary integer Rawdeltaf 2-byte binary integer (not mode AAA) Ysubz 2-byte binary integer (not mode AAA) 4, 8, 12 or 16 bytes of level map information; one byte for each band plus up to 3 bytes to round to the nearest whole word The structure of a VAS area is complicated by two facts. First, every line may not contain all the spectral bands indicated in the band map (W19 in the Area Directory). Second, the order of the bands need not be the same on every line. What does appear on a given line is indicated in the level map section which acts as an index for the line. Only the leftmost N bytes of the level map contain nonzero data, N being the actual number of bands contained in each element of the line. The Ith byte of the level map corresponds to the Ith 16-bit pixel in each element of the line. Unused level bytes are filled with binary zeros, while the data in unused pixel locations may not be zero but in any case should be ignored. GOES CALIBRATION INFORMATION Channel numbers range from 1 to 38. It encodes the type, size and location of the VAS detector used to make the corresponding observation. The meaning of the channel numbers is shown in Table 3. CHANNEL DETECTOR SIZE LOCATION SPECTRAL BAND ------- -------- ---- -------- ------------- 1 HGCDTE LARGE UPPER 1 2 HGCDTE LARGE UPPER 2 3 HGCDTE LARGE UPPER 3 4 HGCDTE LARGE UPPER 4 5 HGCDTE LARGE UPPER 5 6 INSB LARGE UPPER 6 7 HGCDTE LARGE UPPER 7 8 HGCDTE LARGE UPPER 8 9 HGCDTE LARGE UPPER 9 10 HGCDTE LARGE UPPER 10 11 INSB LARGE UPPER 11 12 INSB LARGE UPPER 12 13 HGCDTE LARGE LOWER 1 14 HGCDTE LARGE LOWER 2 15 HGCDTE LARGE LOWER 3 16 HGCDTE LARGE LOWER 4 17 HGCDTE LARGE LOWER 5 18 INSB LARGE LOWER 6 19 HGCDTE LARGE LOWER 7 20 HGCDTE LARGE LOWER 8 21 HGCDTE LARGE LOWER 9 22 HGCDTE LARGE LOWER 10 23 INSB LARGE LOWER 11 24 INSB LARGE LOWER 12 25 HGCDTE SMALL UPPER 3 26 HGCDTE SMALL UPPER 4 27 HGCDTE SMALL UPPER 5 28 HGCDTE SMALL UPPER 7 29 HGCDTE SMALL UPPER 8 30 HGCDTE SMALL UPPER 9 31 HGCDTE SMALL UPPER 10 32 HGCDTE SMALL LOWER 3 33 HGCDTE SMALL LOWER 4 34 HGCDTE SMALL LOWER 5 35 HGCDTE SMALL LOWER 7 36 HGCDTE SMALL LOWER 8 37 HGCDTE SMALL LOWER 9 38 HGCDTE SMALL LOWER 10 ----------------------------------------- TABLE 3 NOTE: Channel 39 exists but has never been put into service. For a given spectral band, only one detector size is employed in any given area. However, two channels representing different positions of the detector for a particular band may appear in a single area. The two channels may not appear on the same line. For example, channels 8 and 20 may appear in the same area, but not channels 8 and 36. MODE AA (W52 is VAS): Transformation of VAS raw IR values into brightness temperatures is accomplished via the intermediate computation of calibrated VAS radiometric values (radiances) using the following formula: R = (P - YSUBZ) * 2**(-N) * 2**(F + DELTAF) where: R = radiance in ergs / (sec*steradians*cm) N = number of bits used from the pixel values; the satellite transmits 10-bit values, but by averaging in the dwell-sound mode, the effective precision is raised to 15 bits. For maximum accuracy, n is set to 15. P = the N most significant bits, not including the sign bit, of the raw pixel value YSUBZ = the N most significant bits, not including the sign bit, of the space view value (from the calibration section of common doc). F = the Channel Scale Factor, a function of the spectral band Table 4 lists the F factors: BAND F factor ------------------------ 1 8 2 8 3 8 4 8 5 8 6 4 7 8 8 8 9 7 10 7 11 4 12 2 ---------------------- TABLE 4 DELTAF = Line Scale Factor. This is a value between -5 and 5 obtained from RAWDELTAF in the calibration section as follows: Set L= the 4 l.s.b. of RAWDELTAF (L=MOD(RAWDELTAF,16)) DELTAF is obtained from Table 5. L DELTAF --------------------- 0 0 1 -1 2 -2 3 -3 4 -4 5 -5 6 N/A 7 N/A 8 0 9 1 10 2 11 3 12 4 13 5 14 N/A 15 N/A -------------------- TABLE 5 MODE AAA (W52 is AAA): Transformation of VAS raw IR values into brightness temperatures is accomplished via the intermediate computation of calibrated VAS radiometric values (radiances) using the following procedure. The following 512 bytes contain the calibration data as follows: W1. Satellite number (SSS code) W2. Date in the form YYDDD W3. Time in the form HHMMSS W4-W79. Radiance Equation Coefficients, array IAB(2,38) W80-W117. Radiance Equation Coefficients Scale Factors, IFAB(38) W118-W128. Zeros The array IAB contains 2 coefficients for each of the 38 channels and IFAB contains 1 scale factor for each channel. If the channel is ICHAN, compute radiance for raw value P by: AB1 = IAB(1,ICHAN) AB2 = IAB(2,ICHAN) FAB = 2.**(15-IFAB(ICHAN)) R = (AB2 * P/32. - AB1) / FAB Note: The raw value P is divided by 32 because the data is stored as 15-bit numbers but the coefficients expect 10-bit numbers. Once the radiance value is calculated, the following converts it to temperature and an 8-bit gray scale value (MODE AA or AAA). Conversion of radiance to brightness temperature is accomplished using the inverse Planck function: T = (C2 * V) / LN((C1 * V**3 / B) + 1) Where: T = temperature (degrees K) V = wave number of the radiation (1/CM) This is the median of the response characteristic of the filter used for a particular spectral band. B = the radiance in (WATTS * CM) / (M**2 * STERADIANS) This is 1000 times the radiance calculated above due to the difference in units. C1 = 1.191066E-8 WATTS/(M**2 * STERADIANS * CM**_4) C1 = 2*H*C**2, where H is Planck's constant (6.63E-34 Joules*secs). C2 = 1.438833 CM * Degrees K C2 = H*C/K, where K is the Boltzmann constant (1.281E-23 Joules/Degree) Following are FORTRAN subroutines to perform these calculations: 1. convert raw data into radiances 2. convert radiances into brightness temperatures 3. convert brightness temperatures to gray scale values This routine is useful for comparing VAS IR to VISSR IR areas. VISSR IR is at approximately the same wavelength as VAS band 8. Each routine uses the output of the previous routine in its calculations. These routines are used on the SSEC IBM McIDAS. REAL FUNCTION RADENC (VALUE, KBAND, DELF, YSUBZ) C CONVERTS RAW DATA VALUES TO RADIANCES C VALUE IS 16 BIT SIGNED RAW DATA VALUE (SIGN BIT ALWAYS EQUAL 0) C I.E. THE TEN BITS OF SATELLITE DATA LEFT JUSTIFIED IN 15 BITS C AND ZERO FILLED TO THE RIGHT. SINCE WE AVERAGE SUCCESSIVE C SPINS IN DWELL SOUND (FOR SAME CHANNEL), WE 'MANUFACTURE' C EXTRA PRECISION AND THAT SHOWS UP IN THE 'EXTRA' LEFT BITS. C KBAND IS SPECTRAL BAND NUMBER (1 TO 12) C DELF IS DECODED DELTA-F TAKE 4 LSBS OF DELTA-F REPORTED C IN THE CALIBRATION SECTION AND DECODE AS FOLLOWS: C 0=>0 , 1=>-1 , 2=>-2 , 3=>-3 , 4=>-4 , 5=>-5 , 6,7=>ILLEGAL C 8=>0 , 9=>1 , 10=>2 , 11=>3 , 12=>4 , 13=>5 , 14,15=>ILLEGAL C YSUBZ IS ALL 16 BITS OF Y SUB Z AS REPORTED IN THE C CALIBRATION SECTION CONVERTED TO REAL FORMAT IMPLICIT INTEGER (A-Z) INTEGER*2 VALUE INTEGER KBAND, DELF REAL YSUBZ REAL AMAX1, FLOAT INTEGER FTABL(12) DATA FTABL /5*8, 4, 8, 8, 7, 7, 4, 2 / TMPVAL = VALUE RADENC = FLOAT(TMPVAL) - YSUBZ RADENC = AMAX1 (0.0, RADENC) RADENC = RADENC * (2. ** (FTABL(KBAND) - 15 + DELF)) RETURN END REAL FUNCTION VASTBB(RAD,KBAND) C CONVERTS RADIANCES TO BRIGHTNESS TEMPERATURES C RAD IS THE RADIANCE COMPUTED BY 'RADENC' C KBAND IS THE SPECTRAL BAND NUMBER (1...12) DIMENSION FK1(12),FK2(12),TC1(12),TC2(12) DATA FK1/.37407E+4,.39157E+4,.40871E+4,.43417E+4,.50296E+4, 1 .12819E+6,.58519E+4,.84911E+4,.30936E+5,.38873E+5,.13611E+6, 2 .19511E+6/ DATA FK2/.97802E+3,.99304E+3,.10073E+4,.10278E+4,.10795E+4, 1 .31767E+4,.11354E+4,.12853E+4,.19778E+4,.21343E+4,.32408E+4, 2 .36542E+4/ DATA TC1/-.89414E-3,.24306E-2,.24917E-2,.34902E-2,.39097E-2, 1 .66916E-1,.63888E-2,.34408,.70558E-1,.78113,.58431E-1,.43968/ DATA TC2/.99998,.99995,.99995,.99993,.99993,.99983,.99991, 1 .99722,.99973,.99717,.99986,.99903/ EXPN=FK1(KBAND)/RAD+1. TT=FK2(KBAND)/ALOG(EXPN) VASTBB=(TT-TC1(KBAND))/TC2(KBAND) RETURN END INTEGER FUNCTION GRYSCL (TEMPK) REAL TEMPK REAL TLIM INTEGER CON1, CON2 DATA CON1 / 418 / DATA CON2 / 660 / DATA TLIM / 242.0 / C CONVERT A BRIGHTNESS TEMPERATURE TO A GREY SCALE FOR DISPLAY C TEMPK---TEMPERATURE IN DEGREES KELVIN IF (TEMPK .GE. TLIM) GOTO 100 GRYSCL = MIN0 (CON1 - IFIX (TEMPK), 255) GOTO 200 100 CONTINUE GRYSCL = MAX0 (CON2 - IFIX (2 * TEMPK), 0) 200 CONTINUE RETURN END GOES NAVIGATION INFORMATION To navigate an image is to associate Earth-based coordinates, usually latitude and longitude, with the pixels of the image. Navigation information for a digital area, when present in the McIDAS system, are supplied in the NAV block of the area file. Navigation information is most useful if used with the McIDAS navigation software used to convert image coordinates (line, element) to geographic coordinates (latitude, longitude). Following is a description of the contents of the NAV block of a McIDAS area file. The byte offset to the block is located by word 53 of the Area Directory entry described above. If present, the individual entries in the NAV block are as follows: W1. Navigation type: 4-byte ASCII. Can be: 'GOES' for GOES satellites 'PS ' for polar stereographic projection 'LAMB' for Lambert Conformal projection ' ' (or binary 0) for no navigation If the type is 'GOES', the succeeding words are: W2. SSSYYDDD: Satellite ID, year, and Julian day W3. HHMMSS Nominal start time of the image Orbit parameters: W4. Orbit type (always 1 for GOES satellite) W5. Epoch date (ETIMY): format YYMMDD W6. Epoch time (ETIMH): HHMMSS W7. Semimajor axis (SEMIMA), Km * 100 W8. Orbital eccentricity (ECCEN), * 100000 W9. Orbital inclination (ORBINC), Deg * 1000 W10. Mean anomaly (MEANA), Deg * 1000 W11. Argument of perigee (PERIGEE), Deg * 1000 W12. Right ascension of the ascending node (ASNODE), Deg * 1000 Attitude parameters: W13. Declination of the satellite axis (DECLIN), DDDMMSS (+ = NORTH) W14. Right ascension of the satellite axis (RASCEN), DDDMMSS W15. Picture center line number (PICLIN) Spin: W16. Spin period (SPINP); either the period of the satellite on the given day, in microseconds, or the spin rate in revolutions/minute Frame geometry: W17. Total sweep angle, line direction (DEGLIN), DDDMMSS W18. Number of scan lines (LINTOT), format NNLLLLL NN is number of sensors; LLLLL is number of scans (total number of actual lines is NN * LLLLL) W19. Total sweep angle, element direction (DEGELE), DDDMMSS W20. Number of elements in a scan line (ELETOT) Camera geometry: W21. Forward-leaning (PITCH), DDDMMSS W22. Sideways-leaning (YAW), DDDMMSS W23. Rotation (ROLL), DDDMMSS Other parameters: W24. (Reserved) W25. East-West adjustment value (IAJUST), in visible elements(+ OR -) W26. Time that IAJUST computed from the first valid landmark of the day (IAJTIM), HHMMSS W27. (Reserved) W28. Angle between the VISSR and sun sensor (ISEANG), DDDMMSS W29. (Reserved for later implementation of *SKEW*; must be 0) W30. (Reserved) BETAS for this area: W31. Scan line of the first beta W32. Time of the above scan line (BEGINNING), HHMMSS W33. Time of the above scan line (CONTINUED), MILLISECONDS * 10 W34. Beta count 1 W35. Scan line of the second beta W36. Time of the above scan line (BEGINNING), HHMMSS W37. Time of the above scan line (CONTINUED), MILLISECONDS * 10 W38. Beta count 2 GAMMA for this area: W39. GAMMA, element offset * 100; this is the nominal offset at time zero of this day W40. GAMMA-DOT, element drift per hour * 100. W41-120. (Reserved) W121-128. Memo entry: up to 32 ASCII characters of comments ********************************************************************* NAVIGATION SUBROUTINE PACKAGE ********************************************************************* Following are listings of FORTRAN/77 subroutines and functions that perform the steps required to navigate a satellite area. If these routines are used, they should be installed on the computer as a package since they share some data through COMMON blocks. Each routine contains a brief description of its function and the required inputs. To navigate a satellite image or area, the parameters contained in the NAV block must be moved to an array called IARR. IARR(1) should contain GOES in ASCII. The call ISTAT= NVXINI(1,IARR) initializes the navigation package. If ISTAT is returned as 0, the initialization was successful. This must be done prior to using the routines in the package. To transform line and element from the satellite area to earth coordinates, use the call ISTAT= NVXSAE(XLIN,XELE,0.,XLAT,XLON,XDUMMY). Inputs are REAL*4 XLIN= satellite line number XELE= satellite element number Both inputs are based on full resolution. Outputs are REAL*4 XLAT= latitude of pixel XLON= longitude of pixel North and West are positive. The parameters, 0. and XDUMMY, are used in 3-dimensional navigation. If ISTAT= 0, the transformation was successful. If ISTAT= -1, the transformation was not possible (e.g. arguments out of range). To transform earth coordinates (latitude, longitude) to satellite coordinates, call ISTAT= NVXEAS(XLAT,XLON,0.,XLIN,XELE,XDUMMY). The arguments and return value are the same as the above call. Normally, transformations are made to and from latitude and longitude. If this is not desired, such as with points near the pole, earth-based coordinates may be changed to rectangular coordinates (X,Y,Z). These coordinates have the origin at the earth's center with the x-axis passing through the Equator at 0 degrees, the y-axis passing through the Equator at 90 degrees east (-90 deg.) and the z-axis passing through the North Pole. The call ISTAT= NVXINI(2,'XYZ') causes the routines NVXSAE and NVXEAS to perform the following functions: ISTAT= NVXSAE(XLIN,XELE,0.,X,Y,Z) where XLIN= satellite line number XELE= satellite element number X,Y,Z are the rectangular coordinates described above. ISTAT= NVXEAS(X,Y,Z,XLIN,XELE,XDUMMY) arguments as described above. It is possible to return to latitude/longitude coordinates with the call ISTAT= NVXINI(2,LL), where LL contains the characters 'LL '. FUNCTION NVXINI(IFUNC,IARR) C C THIS ROUTINE SETS UP COMMON BLOCKS NAVCOM AND NAVINI FOR USE BY THE C NAVIGATION TRANSFORMATION ROUTINES NVXSAE AND NVXEAS. C NVXINI SHOULD BE RECALLED EVERY TIME A TRANSFORMATION IS DESIRED C FOR A PICTURE WITH A DIFFERENT TIME THAN THE PREVIOUS CALL. C IFUNC IS 1 (INITIALIZE; SET UP COMMON BLOCKS) C 2 (ACCEPT/PRODUCE ALL EARTH COORDINATES IN LAT/LON C FORM IF IARR IS 'LL ' OR IN THE X,Y,Z COORDINATE FRAME C IF IARR IS 'XYZ '. C THIS AFFECTS ALL SUBSEQUENT NVXEAS OR NVXSAE CALLS.) C IARR IS AN INTEGER ARRAY (DIM 128) IF IFUNC=1, CONTAINING NAV C PARAMETERS C INTEGER IARR(*),GOES,LL,XYZ CHARACTER*2 CLLSW COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 2,PRERAT,PREDIR,PITCH,YAW,ROLL,SKEW COMMON /BETCOM/IAJUST,IBTCON,NEGBET,ISEANG COMMON /VASCOM/SCAN1,TIME1,SCAN2,TIME2 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 ROTM11,ROTM13,ROTM21,ROTM23,ROTM31,ROTM33, B PICTIM,XREF COMMON /NVUNIT/ LLSW DATA MISVAL/Z80808080/ DATA JINIT/0/,GOES/4HGOES/,LL/4HLL /,XYZ/4HXYZ / C IF (JINIT.EQ.0) THEN JINIT=1 LLSW=0 JDAYSV=-1 JTIMSV=-1 ENDIF IF (IFUNC.EQ.2) THEN IF (IARR(1).EQ.LL) LLSW=0 IF (IARR(1).EQ.XYZ) LLSW=1 NVXINI=0 RETURN ENDIF JTYPE=IARR(1) IF (JTYPE.NE.GOES) GOTO 90 JDAY=IARR(2) JTIME=IARR(3) IF(JDAY.EQ.JDAYSV.AND.JTIME.EQ.JTIMSV) GO TO 10 C C-----INITIALIZE NAVCOM NAVDAY=MOD(JDAY,100000) DO 20 N=7,12 IF(IARR(N).GT.0) GO TO 25 20 CONTINUE GO TO 90 25 IETIMY=ICON1(IARR(5)) IETIMH=100*(IARR(6)/100)+NINT(.6*MOD(IARR(6),100)) SEMIMA=FLOAT(IARR(7))/100.0 OECCEN=FLOAT(IARR(8))/1000000.0 ORBINC=FLOAT(IARR(9))/1000.0 XMEANA=FLOAT(IARR(10))/1000.0 PERHEL=FLOAT(IARR(11))/1000.0 ASNODE=FLOAT(IARR(12))/1000.0 CALL EPOCH(IETIMY,IETIMH,SEMIMA,OECCEN,XMEANA) IF (IARR(5).EQ.0.OR.IARR(9).EQ.0) GOTO 90 DECLIN=FLALO(IARR(13)) RASCEN=FLALO(IARR(14)) PICLIN=IARR(15) IF (IARR(15).GE.1000000) PICLIN=PICLIN/10000. IF (IARR(13).EQ.0.AND.IARR(14).EQ.0.AND.IARR(15).EQ.0) * GOTO 90 SPINRA=IARR(16)/1000.0 IF(IARR(16).NE.0.AND.SPINRA.LT.300.0) SPINRA=60000.0/SPINRA IF (IARR(16).EQ.0) GOTO 90 DEGLIN=FLALO(IARR(17)) LINTOT=IARR(18) DEGELE=FLALO(IARR(19)) IELTOT=IARR(20) PITCH=FLALO(IARR(21)) YAW=FLALO(IARR(22)) ROLL=FLALO(IARR(23)) SKEW=IARR(29)/100000.0 IF (IARR(29).EQ.MISVAL) SKEW=0. C C-----INITIALIZE BETCOM IAJUST=IARR(25) ISEANG=IARR(28) IBTCON=6289920 NEGBET=3144960 C C-----INITIALIZE NAVINI COMMON BLOCK EMEGA=.26251617 AB=40546851.22 ASQ=40683833.48 BSQ=40410330.18 R=6371.221 RSQ=R*R RDPDG=1.745329252E-02 NUMSEN=MOD(LINTOT/100000,100) IF(NUMSEN.LT.1)NUMSEN=1 TOTLIN=NUMSEN*MOD(LINTOT,100000) RADLIN=RDPDG*DEGLIN/(TOTLIN-1.0) TOTELE=IELTOT RADELE=RDPDG*DEGELE/(TOTELE-1.0) PICELE=(1.0+TOTELE)/2.0 CPITCH=RDPDG*PITCH CYAW=RDPDG*YAW CROLL=RDPDG*ROLL PSKEW=ATAN2(SKEW,RADLIN/RADELE) STP=SIN(CPITCH) CTP=COS(CPITCH) STY=SIN(CYAW-PSKEW) CTY =COS(CYAW-PSKEW) STR=SIN(CROLL) CTR=COS(CROLL) ROTM11=CTR*CTP ROTM13=STY*STR*CTP+CTY*STP ROTM21=-STR ROTM23=STY*CTR ROTM31=-CTR*STP ROTM33=CTY*CTP-STY*STR*STP RFACT=ROTM31**2+ROTM33**2 ROASIN=ATAN2(ROTM31,ROTM33) TMPSCL=SPINRA/3600000.0 DEC=DECLIN*RDPDG SINDEC=SIN(DEC) COSDEC=COS(DEC) RAS=RASCEN*RDPDG SINRAS=SIN(RAS) COSRAS=COS(RAS) B11=-SINRAS B12=COSRAS B13=0.0 B21=-SINDEC*COSRAS B22=-SINDEC*SINRAS B23=COSDEC B31=COSDEC*COSRAS B32=COSDEC*SINRAS B33=SINDEC XREF=RAERAC(NAVDAY,0,0.0)*RDPDG C C-----TIME-SPECIFIC PARAMETERS (INCL GAMMA) PICTIM=FLALO(JTIME) GAMMA=FLOAT(IARR(39))/100. GAMDOT=FLOAT(IARR(40))/100. C C-----INITIALIZE VASCOM IF (JDAY/100000.GT.25.AND.IARR(31).GT.0) THEN C THIS SECTION DOES VAS BIRDS C IT USES TIMES AND SCAN LINE FROM BETA RECORDS SCAN1=FLOAT(IARR(31)) TIME1=FLALO(IARR(32)) SCAN2=FLOAT(IARR(35)) TIME2=FLALO(IARR(36)) ELSE C THIS SECTION DOES THE OLD GOES BIRDS SCAN1=1. TIME1=FLALO(JTIME) SCAN2=FLOAT(MOD(LINTOT,100000)) TIME2=TIME1+SCAN2*TMPSCL ENDIF C C-----ALL DONE. EVERYTHING OK 10 CONTINUE JDAYSV=JDAY JTIMSV=JTIME NVXINI=0 RETURN 90 NVXINI=-1 RETURN END C ******************************************************************** C FUNCTION NVXSAE(XLIN,XELE,XDUM,XPAR,YPAR,ZPAR) C TRANSFORMS SAT COOR TO EARTH COOR. C ALL PARAMETERS REAL*4 C INPUTS: C XLIN,XELE ARE SATELLITE LINE AND ELEMENT (IMAGE COORDS.) C XDUM IS DUMMY (IGNORE) C OUTPUTS: C XPAR,YPAR,ZPAR REPRESENT EITHER LAT,LON,(DUMMY) OR X,Y,Z DEPENDING C ON THE OPTION SET IN PRIOR NVXINI CALL WITH IFUNC=2. C FUNC VAL IS 0 (OK) OR -1 (CAN'T; E.G. OFF OF EARTH) C COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,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 ROTM11,ROTM13,ROTM21,ROTM23,ROTM31,ROTM33, B PICTIM,XREF COMMON /NVUNIT/ LLSW DATA PI/3.14159265/ C ILIN=NINT(XLIN) PARLIN=(ILIN-1)/NUMSEN+1 FRAMET=TMPSCL*PARLIN SAMTIM=FRAMET+PICTIM CALL SATVEC(SAMTIM,XSAT,YSAT,ZSAT) YLIN=(XLIN-PICLIN)*RADLIN YELE=(XELE-PICELE+GAMMA+GAMDOT*SAMTIM)*RADELE XCOR=B11*XSAT+B12*YSAT+B13*ZSAT YCOR=B21*XSAT+B22*YSAT+B23*ZSAT ROT=ATAN2(YCOR,XCOR)+PI YELE=YELE-ROT COSLIN=COS(YLIN ) SINLIN=SIN(YLIN) SINELE=SIN(YELE) COSELE=COS(YELE) ELI=ROTM11*COSLIN-ROTM13*SINLIN EMI=ROTM21*COSLIN-ROTM23*SINLIN ENI=ROTM31*COSLIN-ROTM33*SINLIN TEMP=ELI ELI=COSELE*ELI+SINELE*EMI EMI=-SINELE*TEMP+COSELE*EMI ELO=B11*ELI+B21*EMI+B31*ENI EMO=B12*ELI+B22*EMI+B32*ENI ENO=B13*ELI+B23*EMI+B33*ENI BASQ=BSQ/ASQ ONEMSQ=1.0-BASQ AQ=BASQ+ONEMSQ*ENO**2 BQ=2.0*((ELO*XSAT+EMO*YSAT)*BASQ+ENO*ZSAT) CQ=(XSAT**2+YSAT**2)*BASQ+ZSAT**2-BSQ RAD=BQ**2-4.0*AQ*CQ IF(RAD.LT.1.0)GO TO 2 S=-(BQ+SQRT(RAD))/(2.0*AQ) X=XSAT+ELO*S Y=YSAT+EMO*S Z=ZSAT+ENO*S CT=COS(EMEGA*SAMTIM+XREF) ST=SIN(EMEGA*SAMTIM+XREF) X1=CT*X+ST*Y Y1=-ST*X+CT*Y IF (LLSW.EQ.0) THEN CALL NXYZLL(X1,Y1,Z,XPAR,YPAR) ZPAR=0. ELSE XPAR=X1 YPAR=Y1 ZPAR=Z ENDIF NVXSAE=0 RETURN 2 NVXSAE=-1 RETURN END C ******************************************************************** C FUNCTION NVXEAS(XPAR,YPAR,ZPAR,XLIN,XELE,XDUM) C 5/26/82; TRANSFORM EARTH TO SATELLITE COORDS C ALL PARAMETERS REAL*4 C INPUTS: C XPAR,YPAR,ZPAR REPRESENT EITHER LAT,LON,(DUMMY) OR X,Y,Z DEPENDING C ON THE OPTION SET IN PRIOR NVXINI CALL WITH IFUNC=2. C OUTPUTS: C XLIN,XELE ARE SATELLITE LINE AND ELEMENT (IMAGE COORDS.) C XDUM IS DUMMY (IGNORE) C FUNC VAL IS 0 (OK) OR -1 (CAN'T; E.G. BAD LAT/LON) C 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 ROTM11,ROTM13,ROTM21,ROTM23,ROTM31,ROTM33, B PICTIM,XREF COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 2,PRERAT,PREDIR,PITCH,YAW,ROLL,SKEW COMMON/VASCOM/SCAN1,TIME1,SCAN2,TIME2 COMMON /NVUNIT/ LLSW DATA OLDLIN/910./,ORBTIM/-99999./,MAXITR/5/ C NVXEAS=0 IF (LLSW.EQ.0) THEN IF (ABS(XPAR).GT.90.) THEN NVXEAS=-1 RETURN ENDIF CALL NLLXYZ(XPAR,YPAR,X1,Y1,Z) ELSE X1=XPAR Y1=YPAR Z=ZPAR ENDIF XDUM=0.0 SAMTIM=TIME1 DO 50 I=1,2 IF(ABS(SAMTIM-ORBTIM).LT.0.0005) GO TO 10 CALL SATVEC(SAMTIM,XSAT,YSAT,ZSAT) ORBTIM=SAMTIM XHT=SQRT(XSAT**2+YSAT**2+ZSAT**2) 10 CT=COS(EMEGA*SAMTIM+XREF) ST=SIN(EMEGA*SAMTIM+XREF) X=CT*X1-ST*Y1 Y= ST*X1+CT*Y1 VCSTE1=X-XSAT VCSTE2=Y-YSAT VCSTE3=Z-ZSAT VCSES3=B31*VCSTE1+B32*VCSTE2+B33*VCSTE3 ZNORM=SQRT(VCSTE1**2+VCSTE2**2+VCSTE3**2) X3=VCSES3/ZNORM UMV=ATAN2(X3,SQRT(RFACT-X3**2))-ROASIN XLIN=PICLIN-UMV/RADLIN PARLIN=IFIX(XLIN-1.0)/NUMSEN IF(I.EQ.2) GO TO 50 SAMTIM=TIME2 OLDLIN=XLIN 50 CONTINUE SCNNUM=(IFIX(OLDLIN+XLIN)/2.0-1.0)/8.0 SCNFRC=(SCNNUM-SCAN1)/(SCAN2-SCAN1) XLIN=OLDLIN+SCNFRC*(XLIN-OLDLIN) SAMTIM=TIME1+TMPSCL*(SCNNUM-SCAN1) CALL SATVEC(SAMTIM,XSAT,YSAT,ZSAT) COSA=X*XSAT+Y*YSAT+Z*ZSAT CTST=0.0001*R*XHT+RSQ IF(COSA.LT.CTST) NVXEAS=-1 XSATS1=B11*XSAT+B12*YSAT+B13*ZSAT YSATS2=B21*XSAT+B22*YSAT+B23*ZSAT CT=COS(EMEGA*SAMTIM+XREF) ST=SIN(EMEGA*SAMTIM+XREF) X=CT*X1-ST*Y1 Y= ST*X1+CT*Y1 VCSTE1=X-XSAT VCSTE2=Y-YSAT VCSTE3=Z-ZSAT VCSES1=B11*VCSTE1+B12*VCSTE2+B13*VCSTE3 VCSES2=B21*VCSTE1+B22*VCSTE2+B23*VCSTE3 VCSES3=B31*VCSTE1+B32*VCSTE2+B33*VCSTE3 XNORM=SQRT(ZNORM**2-VCSES3**2) YNORM=SQRT(XSATS1**2+YSATS2**2) ZNORM=SQRT(VCSTE1**2+VCSTE2**2+VCSTE3**2) X3=VCSES3/ZNORM UMV=ATAN2(X3,SQRT(RFACT-X3**2))-ROASIN SLIN=SIN(UMV) CLIN=COS(UMV) U=ROTM11*CLIN+ROTM13*SLIN V=ROTM21*CLIN+ROTM23*SLIN XELE=PICELE+ASIN((XSATS1*VCSES2-YSATS2*VCSES1)/(XNORM*YNORM))/RADE 1LE XELE=XELE+ATAN2(V,U)/RADELE XELE=XELE-GAMMA-GAMDOT*SAMTIM RETURN END C ******************************************************************** C FUNCTION NVXOPT(IFUNC,XIN,XOUT) C IFUNC= 'SPOS' SUBSATELLITE LAT/LON C XIN - NOT USED C XOUT - 1. SUB-SATELLITE LATITUDE C - 2. SUB-SATELLITE LONGITUDE INTEGER SPOS DATA SPOS/4HSPOS/ REAL*4 XIN(*),XOUT(2) CHARACTER*4 CLIT,CFUNC CHARACTER*12 CFF NVXOPT=0 IF(IFUNC.EQ.SPOS) THEN CALL SATPOS(0,ITIME(PICTIM),X,Y,Z) CALL NXYZLL(X,Y,Z,XOUT(1),XOUT(2)) ELSE NVXOPT=1 ENDIF RETURN END C ******************************************************************** C SUBROUTINE SATPOS(INORB,NTIME,X,Y,Z) C $ CALCULATE SATELLITE POSITION VECTOR FROM THE EARTH'S CENTER. C $ INORB = (I) INPUT INITIALIZATION FLAG - SHOULD = 0 ON FIRST CALL TO C $ SATPOS, 1 ON ALL SUBSEQUENT CALLS. C $ NTIME = (I) INPUT TIME (HOURS, MINUTES, SECONDS) IN HHMMSS FORMAT C $ X, Y, Z = (R) OUTPUT COORDINATES OF POSITION VECTOR C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 DEGLIN,DEGELE,SPINRA,SEMIMA,OECCEN,ORBINC,PERHEL REAL*4 ASNODE,DECLIN,RASCEN,PICLIN,PRERAT,PREDIR,PITCH,YAW REAL*4 ROLL,SKEW REAL*4 X,Y,Z COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 2,PRERAT,PREDIR,PITCH,YAW,ROLL,SKEW C IF(INORB.NE.0)GO TO 1 INORB=1 PI=3.14159265D0 RDPDG=PI/180.0 RE=6378.388 GRACON=.07436574D0 SOLSID=1.00273791D0 SHA=100.26467D0 SHA=RDPDG*SHA IRAYD=74001 IRAHMS=0 O=RDPDG*ORBINC P=RDPDG*PERHEL A=RDPDG*ASNODE SO=SIN(O) CO=COS(O) 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=GRACON*RE*DSQRT(RE/SEMIMA)/SEMIMA 1 DIFTIM=TIMDIF(IETIMY,IETIMH,NAVDAY,NTIME) XMANOM=XMMC*DIFTIM ECANM1=XMANOM EPSILN=1.0E-8 DO 2 I=1,20 ECANOM=XMANOM+OECCEN*DSIN(ECANM1) IF(DABS(ECANOM-ECANM1).LT.EPSILN)GO TO 3 2 ECANM1=ECANOM 3 XOMEGA=DCOS(ECANOM)-OECCEN YOMEGA=SROME2*DSIN(ECANOM) XS=XOMEGA*PX+YOMEGA*QX YS=XOMEGA*PY+YOMEGA*QY ZS=XOMEGA*PZ+YOMEGA*QZ DIFTIM=TIMDIF(IRAYD,IRAHMS,NAVDAY,NTIME) RA=DIFTIM*SOLSID*PI/720.0D0+SHA RAS=DMOD(RA,2.0*PI) CRA=COS(RAS) SRA=SIN(RAS) X=CRA*XS+SRA*YS Y=-SRA*XS+CRA*YS Z=ZS RETURN END C ******************************************************************** C From this point on, the routines included are the McIDAS C subroutines called by the navigation routines to perform specific low C level tasks. The FORTRAN routines can be compiled along with the C routines above. The assembler routine must either be assembled or C rewritten in FORTRAN. The rewriting is not recommended. C ******************************************************************** C FUNCTION ICON1(YYMMDD) C CONVERTS YYMMDD TO YYDDD IMPLICIT INTEGER(A-Z) DIMENSION NUM(12) DATA NUM/0,31,59,90,120,151,181,212,243,273,304,334/ C YEAR=MOD(YYMMDD/10000,100) MONTH=MOD(YYMMDD/100,100) DAY=MOD(YYMMDD,100) IF(MONTH.LT.0.OR.MONTH.GT.12)MONTH=1 JULDAY=DAY+NUM(MONTH) IF(MOD(YEAR,4).EQ.0.AND.MONTH.GT.2) JULDAY=JULDAY+1 ICON1=1000*YEAR+JULDAY RETURN END C ******************************************************************** C SUBROUTINE EPOCH(IETIMY,IETIMH,SEMIMA,OECCEN,XMEANA) C EPOCH PHILLI 0173 NAV: FINDS TIME OF PERIGEE FROM KEPLERIAN EPOCH PARAMETER (PI=3.14159265) PARAMETER (RDPDG=PI/180.0) PARAMETER (RE=6378.388) PARAMETER (GRACON=0.07436574) LEAPYR(IY)=366-(MOD(IY,4)+3)/4 C XMMC=GRACON*SQRT(RE/SEMIMA)**3 XMANOM=RDPDG*XMEANA TIME=(XMANOM-OECCEN*SIN(XMANOM))/(60.0*XMMC) TIME1=FLALO(IETIMH) TIME=TIME1-TIME IDAY=0 IF(TIME.GT.48.0)GO TO 8 IF(TIME.GT.24.0)GO TO 1 IF(TIME.LT.-24.0)GO TO 2 IF(TIME.LT.0.0)GO TO 3 GO TO 4 8 TIME=TIME-48.0 IDAY=2 GO TO 4 1 TIME=TIME-24.0 IDAY=1 GO TO 4 2 TIME=TIME+48.0 IDAY=-2 GO TO 4 3 TIME=TIME+24.0 IDAY=-1 4 IETIMH=ITIME(TIME) IF(IDAY.EQ.0)RETURN JYEAR=MOD(IETIMY/1000,100) JDAY=MOD(IETIMY,1000) JDAY=JDAY+IDAY IF(JDAY.LT.1)GO TO 5 JTOT=LEAPYR(JYEAR) IF(JDAY.GT.JTOT)GO TO 6 GO TO 7 5 JYEAR=JYEAR-1 JDAY=LEAPYR(JYEAR)+JDAY GO TO 7 6 JYEAR=JYEAR+1 JDAY=JDAY-JTOT 7 IETIMY=1000*JYEAR+JDAY RETURN END C ******************************************************************** C FUNCTION RAERAC(IYRDY,IHMS,RAE) C RAERAC PHILLI 0174 NAV: CONVERTS EARTH LON TO CELESTIAL LON C INPUT PARAMETERS C IYRDY = YEAR AND JULIAN DAY (INTEGER YYDDD) C IHMS = TIME (INTEGER HHMMSS) C RAE = EARTH LONGITUDE (REAL*4 DEGREES) C FUNCTION VALUE IS CELESTIAL LONGITUDE IN REAL*4 DEGREES C REF: ESCOBAL, "METHODS OF ORBIT DETERMINATION." WILEY & SONS, 1965 C GREENWICH SIDEREAL TIME := ANGLE BETWEEN PRIME MERID. & 0 DEG R.A. C JULIAN DATE := # DAYS ELAPSED SINCE 12 NOON ON JAN 1, 4713 B.C. C APPROXIMATE FORMULA FOR GREENWICH SIDEREAL TIME AT 0Z: C G.S.T (DEG) = S(0) = 99.6909833 + 36000.7689*C + 0.00038708*C*C ,WHERE C C = TIME IN CENTURIES = (J.D. - 2415020.0) / 36525 C FOR G.S.T. AT OTHER TIMES OF (SAME) DAY, USE: C G.S.T AT TIME T = S(T) = S(0) + (T * DS/DT) C DS/DT = 1 + (1 / 365.24219879) = 1.00273790927 REVOLUTIONS/DAY C = 0.250684477 DEGREES/MINUTE C FROM TABLE, J.D. AT 0Z ON JAN 1,1974 = 2442048.5 C THEN S(0) AT 0Z ON JAN 1,1974 = 100.2601800 C DOUBLE PRECISION TIMDIF,RAHA,SOLSID,SHA SHA=100.26467D0 IRAYD=74001 IRAHMS=0 SOLSID=1.00273791D0 RAHA=RAE+TIMDIF(IRAYD,IRAHMS,IYRDY,IHMS)*SOLSID/4.0D0+SHA RAC=DMOD(RAHA,360.0D0) IF(RAC.LT.0.0)RAC=RAC+360.0 RAERAC=RAC RETURN END C ******************************************************************** C FUNCTION RACRAE(IYRDY,IHMS,RAC) C RACRAE PHILLI 0174 NAV: CONVERTS CELESTIAL ONG TO EARTH LON C INPUT PARAMETERS-- C IYRDY = YEAR AND JULIAN DAY (INTEGER YYDDD) C IHMS = TIME (INTEGER HHMMSS) C RAC = CELESTIAL LONGITUDE (REAL*4 DEGREES) C FUNCTION VALUE IS EARTH LONGITUDE IN REAL*4 DEGREES C DOUBLE PRECISION TIMDIF,RAHA,SOLSID,SHA SHA=100.26467D0 IRAYD=74001 IRAHMS=0 SOLSID=1.00273791D0 RAHA=RAC-SHA+TIMDIF(IYRDY,IHMS,IRAYD,IRAHMS)*SOLSID/4.0D0 RAE=DMOD(RAHA,360.0D0) IF(RAE.LT.0.0)RAE=RAE+360.0 RACRAE=RAE RETURN END C ******************************************************************** C FUNCTION TIMDIF(IYRDA1,IHMS1,IYRDA2,IHMS2) C TIME DIFFERENCE IN MINUTES C INPUTS (ALL INTEGER) C IYRDA1 FIRST YEAR AND JULIAN DAY (YYDDD) C IHMS1 FIRST TIME (HHMMSS) C IYRDA2 SECOND YEAR AND DAY (YYDDD) C IHMS2 SECOND TIME (HHMMSS) C FUNC VAL (REAL*8) IS TIME DIFFERENCE IN MINUTES (POSITIVE IF C FIRST DAY/TIME IS THE EARLIER OF THE TWO). C DOUBLE PRECISION TIMDIF,D1,D2,T1,T2 IY1=MOD(IYRDA1/1000,100) ID1=MOD(IYRDA1,1000) IFAC1=(IY1-1)/4+1 D1=365*(IY1-1)+IFAC1+ID1-1 IY2=MOD(IYRDA2/1000,100) ID2=MOD(IYRDA2,1000) IFAC2=(IY2-1)/4+1 D2=365*(IY2-1)+IFAC2+ID2-1 T1=1440.D0*D1+60.D0*FLALO(IHMS1) T2=1440.D0*D2+60.D0*FLALO(IHMS2) TIMDIF=T2-T1 RETURN END C ******************************************************************** C C VECTOR EARTH-CENTER-TO-SAT (FUNC OF TIME) SUBROUTINE SATVEC(SAMTIM,X,Y,Z) DOUBLE PRECISION TWOPI,PI720,DE,TE,DRA,TRA,DNAV,TDIFRA,TDIFE DOUBLE PRECISION PI,RDPDG,RE,GRACON,SOLSID,SHA DOUBLE PRECISION DIFTIM,ECANM1,ECANOM,RA,XMANOM,TIMDIF DOUBLE PRECISION DABS,DMOD,DSQRT,DSIN,DCOS,DATAN2 COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 2,PRERAT,PREDIR,PITCH,YAW,ROLL,SKEW DATA NAVSAV/0/ IF(NAVDAY.EQ.NAVSAV) GO TO 10 NAVSAV=NAVDAY PI=3.14159265D0 TWOPI=2.0*PI PI720=PI/720. RDPDG=PI/180.0 RE=6378.388 GRACON=.07436574D0 SOLSID=1.00273791D0 SHA=100.26467D0 SHA=RDPDG*SHA IRAYD=74001 IRAHMS=0 O=RDPDG*ORBINC P=RDPDG*PERHEL A=RDPDG*ASNODE SO=SIN(O) CO=COS(O) 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=GRACON*RE*DSQRT(RE/SEMIMA)/SEMIMA IEY=MOD(IETIMY/1000,100) IED=MOD(IETIMY,1000) IEFAC=(IEY-1)/4+1 DE=365*(IEY-1)+IEFAC+IED-1 TE=1440.0*DE+60.0*FLALO(IETIMH) IRAY=IRAYD/1000 IRAD=MOD(IRAYD,1000) IRAFAC=(IRAY-1)/4+1 DRA=365*(IRAY-1)+IRAFAC+IRAD-1 TRA=1440.0*DRA+60.0*FLALO(IRAHMS) INAVY=MOD(NAVDAY/1000,100) INAVD=MOD(NAVDAY,1000) INFAC=(INAVY-1)/4+1 DNAV=365*(INAVY-1)+INFAC+INAVD-1 TDIFE=DNAV*1440.-TE TDIFRA=DNAV*1440.-TRA EPSILN=1.0E-8 10 TIMSAM=SAMTIM*60.0 DIFTIM=TDIFE+TIMSAM XMANOM=XMMC*DIFTIM ECANM1=XMANOM DO 2 I=1,20 ECANOM=XMANOM+OECCEN*DSIN(ECANM1) IF(DABS(ECANOM-ECANM1).LT.EPSILN)GO TO 3 2 ECANM1=ECANOM 3 XOMEGA=DCOS(ECANOM)-OECCEN YOMEGA=SROME2*DSIN(ECANOM) Z =XOMEGA*PZ+YOMEGA*QZ Y =XOMEGA*PY+YOMEGA*QY X =XOMEGA*PX+YOMEGA*QX RETURN END C ******************************************************************** C SUBROUTINE NLLXYZ(XLAT,XLON,X,Y,Z) C-----CONVERT LAT,LON TO EARTH CENTERED X,Y,Z C (DALY, 1978) C-----XLAT,XLON ARE IN DEGREES, WITH NORTH AND WEST POSITIVE C-----X,Y,Z ARE GIVEN IN KM. THEY ARE THE COORDINATES IN A RECTANGULAR C FRAME WITH ORIGIN AT THE EARTH CENTER, WHOSE POSITIVE C X-AXIS PIERCES THE EQUATOR AT LON 0 DEG, WHOSE POSITIVE Y-AXIS C PIERCES THE EQUATOR AT LON 90 DEG, AND WHOSE POSITIVE Z-AXIS C INTERSECTS THE NORTH POLE. DATA ASQ/40683833.48/,BSQ/40410330.18/,AB/40546851.22/ DATA RDPDG/1.7453292E-02/,PI/3.1415926/ YLAT=RDPDG*XLAT C-----CONVERT TO GEOCENTRIC (SPHERICAL) LATITUDE YLAT=ATAN2(BSQ*SIN(YLAT),ASQ*COS(YLAT)) YLON=-RDPDG*XLON SNLT=SIN(YLAT) CSLT=COS(YLAT) CSLN=COS(YLON) SNLN=SIN(YLON) TNLT=(SNLT/CSLT)**2 R=AB*SQRT((1.0+TNLT)/(BSQ+ASQ*TNLT)) X=R*CSLT*CSLN Y=R*CSLT*SNLN Z=R*SNLT RETURN END C ******************************************************************** C SUBROUTINE NXYZLL(X,Y,Z,XLAT,XLON) C-----CONVERT EARTH-CENTERED X,Y,Z TO LAT & LON C-----X,Y,Z ARE GIVEN IN KM. THEY ARE THE COORDINATES IN A RECTANGULAR C COORDINATE SYSTEM WITH ORIGIN AT THE EARTH CENTER, WHOSE POS. C X-AXIS PIERCES THE EQUATOR AT LON 0 DEG, WHOSE POSITIVE Y-AXIS C PIERCES THE EQUATOR AT LON 90 DEG, AND WHOSE POSITIVE Z-AXIS C INTERSECTS THE NORTH POLE. C-----XLAT,XLON ARE IN DEGREES, WITH NORTH AND WEST POSITIVE DATA RDPDG/1.745329252E-2/ DATA ASQ/40683833.48/,BSQ/40410330.18/,AB/40546851.22/ C XLAT=100.0 XLON=200.0 IF(X.EQ.0..AND.Y.EQ.0..AND.Z.EQ.0.) GO TO 90 A=ATAN(Z/SQRT(X*X+Y*Y)) C-----CONVERT TO GEODETIC LATITUDE XLAT=ATAN2(ASQ*SIN(A),BSQ*COS(A))/RDPDG XLON=-ATAN2(Y,X)/RDPDG 90 RETURN END C ******************************************************************** C FUNCTION FLALO(M) C INTEGER ANGLE (FORMAT DDDMMSS) OR TIME (FORMAT HHMMSS) TO REAL*4 C M=ANGLE OR TIME C N=IABS(M) 2 FLALO=FLOAT(N/10000)+FLOAT(MOD(N/100,100))/60.0+FLOAT(MOD(N,100))/ *3600.0 IF (M.LT.0) FLALO=-FLALO RETURN END C ******************************************************************** C FUNCTION ITIME(X) C $ FLOATING POINT TIME TO PACKED INTEGER ( SIGN HH MM SS ) C $ X = (R) INPUT FLOATING POINT TIME C IF(X.LT.0.0)GO TO 1 Y=X I=1 GO TO 2 1 Y=-X I=-1 2 J=3600.0*Y+0.5 ITIME=10000*(J/3600)+100*MOD(J/60,60)+MOD(J,60) ITIME=I*ITIME RETURN END C ******************************************************************** C Produced by: University of Wisconsin Space Science & Engineering Center