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
2
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
3
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
4
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
5
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
6
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.
7
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).
8
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
9
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
10
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.
11
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
12
*********************************************************************
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:
13
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
14
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
15
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.
16
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))
17
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
18
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
19
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
20
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)
21
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
22
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))
23
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
24
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.
25
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
26
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
27
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
28
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)
29
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:
30
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
31
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)
32
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)
33
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
34
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
