United States Atmospheric Research and Exposure EPA/600/8-88/003
Environmental Protection Assessment Laboratory March, 1989
Agency Research Triangle Park, NC 27711
Research and Development
S.EPA
Project Report
The Complex Terrain
Dispersion Model
Terrain Preprocessor
System
User's Guide and
Program Descriptions
-------
United States Atmospheric Sciences
Environmental Protection Research Laboratory
Agency Basearch Triangle Park NC 27711
Research and Development December 1987
EPA/600/8-88/003
PROJECT REPORT
THE COMPLEX TERRAIN DISPERSION MODEL
TERRAIN PREPROCESSOR SYSTEM
USER'S GUIDE AND
PROGRAM DESCRIPTION
-------
United States Atmospheric Sciences
Environmental Protection Research Laboratory
Agency Research Triangle Park NC 27711
Reeearch and Development DeceBber 1987
f/EPA
PROJECT REPORT
THE COMPLEX TERRAIN DISPERSION MODEL
TERRAIN PREPROCESSOR SYSTEM
USER'S GUIDE AND
PROGRAM DESCRIPTION
-------
-------
THE COMPLEX TERRAIN DISPERSION MODEL (CTDM)
TERRAIN PREPROCESSOR SYSTEM - USER GUIDE
AND PROGRAM DESCRIPTION
by
Michael T. Mills1
Robert J. Paine1
Elizabeth M. Insley2
Bruce A. Egan1
^•ERT, Inc.
696 Virginia Road, Concord, Massachusetts 01742
2Sigina Research Corp.
394 Lowell Street. Suite 12
Lexington, Massachusetts 02173
Contract Ho. 68-02-3421
Project Officer
Peter L. Finkelstein
Meteorology Division
Atmospheric Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ATMOSPHERIC SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
The Information 1n this document has been funded by the United
States Environmental Protection Agency under Contract No. 68-02-3421
to ERT, Inc. It has been subjected to the Agency's peer and admini-
strative review, and it has been approved for publication as an EPA
document. Mention of trade names or commercial products does not
constitute endorsement or recommendation for use.
11
-------
ABSTRACT
This report describes the operation of a terrain preprocessor which
approximates actual terrain features with mathematical functions. The best-
fit parameters for these functions are used by the Complex Terrain Dispersion
Model (CTDH) in the calculation of lateral and vertical streamline
displacement, an important step in the calculation of concentrations at hill
receptor locations.
The CTDM Terrain Preprocessor is a series of 3 programs which process
digitized contour data to provide hill shape parameters in a format suitable
for direct input to CTDH. The first program, FITCOM, accepts as input a
user-defined hill in terms of its maximum elevation and the x,y coordinates
of the hill center. The elevation and point coordinates of individual
contours are then input from a master file. After evaluation and editing,
each contour is processed by numerical integration to determine the following
parameters for an equivalent ellipse: semi-major and semi-minor axis
lengths; contour centreid coordinates; and the orientation of the ellipse.
These parameters are input to the second preprocessor program, HCRIT, which
determines, for the portion of the hill above a given critical elevation, the
best-fit inverse polynomial profiles along the hill's major and minor axes.
The center coordinates of the fitted hill are calculated as the mean of the
ellipse center coordinates for those contours above a given critical
elevation. The orientation of the fitted hill is calculated as a vector
average of the ellipse orientations weighted by the ellipse eccentricity.
HCRIT provides an input file for CTDtt which contains the following
information for each critical elevation:
e Ellipse parameters corresponding to contours at user-specified
elevations
• Coordinates of the center of the fitted hill
e Orientation of the major axis of the fitted hill with respect to
north
e The length scale and exponents for the inverse polynominal fits
along the hill major and minor axes.
•
The third program, PLOTCOH, generates the following screen displays to aid in
the evaluation of the hill fitting process:
e Map of input contours
e Map of digitized contours which have been qualified and edited
e Map of the digitized contours and their associated fitted ellipses
111
-------
• For each critical cut-off elevation, a map showing the digitized
contours and the contours of the fitted hill at elevations
corresponding to the elevations of those digitized contours above
the critical elevation.
This report has been submitted in partial fulfillment of Contract
68-02-3421 by ERT under the sponsorship of the U.S. Environmental Protection
Agency.
1v
-------
CONTENTS
Abstract
Figures vi
Tables vi
Symbols and Abbreviations vii
1. Introduction 1
1.1 Development of the Complex Terrain Dispersion Model ... 1
1.2 Requirements for a CTDM Terrain Preprocessor Program. . . 1
1.3 Summary of Preprocessor Operation 2
1.4 Organization of the Manual 4
2. Rules for Terrain Fitting 5
2.1 Fitting of Ellipses to Digitized Contours 5
2.1.1 Rules for Contour Digitization 5
2.1.2 Contour'Qualification and Editing 6
2.1.3 Calculation of the Area and Centreid of
Each Digitized Contour 18
2.1.4 Determination of Contour Orientation 18
2.1.5 Calculation of the Semi-Major and Semi-Minor
Axis Lengths for the Contour Elliptical
Representation 20
2.2 Mathematical Representation of a Cut-Off Hill 21
2.2.1 'Best-Fit Ellipse at a Critical Elevation 21
2.2.2 The Inverse Polynomial Profile 24
3. System Operation 29
3.1 FITCON 29
3.2 HCRIT 33
3.3 PLOTCOH 35
4. Input Requirements . 40
4.1 FITCON 40
4.2 HCRIT 41
4.3 PLOTCOH 43
5. Output Description 45
5.1 FITCOM 45
5.2 HCRIT 47
5.3 PLOTCOH ^ 51
References 54
APPENDIXES
A. Selection of Terrain Features for CTDM 55
B. Derivation of Equations for the Area, Centreid
Coordinates, and Second Moment of a Digitized Contour . .
C. Program Test Case 68
D. Program Listings 113
-------
FIGURES
Number Page
1 Example of digitized contour 7
2 The closing of an incomplete contour as performed by
program FITCON 9
3 Illustration of point filtering and contour completion . . 11
4 Selection of an acceptance angle for contour completion . . 12
5 Contour completion with the hill center used as the
reflection point - acceptance angle equals 180 degrees . . 14
6 Contour completion with the hill center used as the
reflection point - acceptance angle equals 31.9 degrees . . 15
7 Contour completion with the incomplete contour centroid
used as the reflection point - acceptance angle equals
51.5 degrees 16
8 * Analysis of multiple contours at the same elevation .... 17
9 Axes for the calculation of second moments 19
10 Inverse polynomial hill profile for 5 different exponent
values 25
TABLES
Bomber Page
1 Format for the CONTOUR Master File 42
2 Format for the plot file generated by FITCON 48
3 Format for the file generated by FITCON for input to HCRIT. 50
4 Format for the plot file generated by HCRIT 52
5 Format for the file generated by HCRIT for Input to CTDM . 53
vi
-------
LIST OF SYMBOLS AND ABBREVIATIONS
SYMBOL
A
a
abs(A)
E
ECC
ECCj
h
H,.
HC
Area enclosed by a polygon formed by the straight lines
connecting digitized contour points
Area of trapezoidal element k
Length of an ellipse semi-major axis
Absolute value of A
Calculated semi-major axis length for contour j
Length of an ellipse semi-minor axis
Calculated semi-minor axis length for contour j
Quantity minimized to determine best fit values of L and P
Perpendicular distance from contour point k to an axis line,
which lies within the contour plane and passes through the
contour centroid
Elevation used in the interpolation of ellipse parameters
Ellipse eccentricity
Calculated eccentricity for the elliptical representation of
contour j
Elevation of a point on the hill surface
Critical dividing streamline height
Elevation of contour j
Critical cutoff elevation
Elevation of the uppermost contour on the hill
Elevation of hill top
ith critical height
Inverse polynominal length scale
The x-component of the first moment of trapezoidal element k
-------
M2 Second 'moment of an ellipse about its minor axis
M2fc Second- moment for trapezoidal segment k about a line passing
through the contour centroid
Nc Number of contours for a hill above the critical cutoff elevation
HHC Number of critical elevations
Np Number of digitized points for a contour
« 3.14159265 ...
P Inverse polynominal exponent
Pc Interpolated ellipse parameter (not including the orientation)
r Calculated radius of a contour determined to be circular
R Maximum value of R^ for m » 1,18
Radius of gyration of the digitized contour about a line passing
through the contour centroid and making an angle 6^ with
respect to the positive x-axis
Second moment of the digitized contour about a line passing
through the contour centroid and making an angle 6^ with
respect to the positive x-axis
oa Standard deviation parameter for the Gaussian terrain
distribution shape along the major axis
«b Standard deviation parameter for the Gaussian terrain
distribution snaps along the minor axis
6 Orientation angle for the minor axis of the fitted hill
6e Orientation angle obtained for a critical elevation by vector
interpolation.
6j Value of ^ associated with the maximum value of Rgg for
the contour j
OB Angle with respect to the x-axis of the m th line through the
centroid of the contour in the plane of the contour (m * 1,18)
Wfc ' Distance along an axis line, within the contour plane and passing
through the contour centroid, between the intersection of
perpendiculars from adjacent contour points (k and k+1)
x Distance from a specified origin toward the east
x* Distance measured along the hill major axis
Xg Calculated x-coordinate for the contour centroid
vi i i
-------
Xej x-coordinate of the centreid of digitized contour j
Xjl x-coordinate of the centreid of a fitted inverse polynomial hill
Xj Major or minor axis length for contour j
Xfc x-coordinate for contour point k
y Distance from a specified origin toward the north
y* Distance measured along the hill minor axis
Yc Calculated y-coordinate for the contour centroid
Ycj y-coordinate of the centroid of digitized contour j
YH y-coordinate of the centroid of a fitted inverse polynomial hill
yjt y-coordinate for contour point k
ABBREVIATION
ASRL Atmospheric Sciences Research Laboratory
CCB Cinder Cone Butte
CTMD Complex Terrain Model Development
CTDM Complex Terrain Dispersion Model
DOS Disk Operating System
EPA Environmental Protection Agency
EPRI Electric Power Research Institute
FSPS Full Scale Plume Study
HBR Hogback Ridge
SHIS Small Hill Impaction Study
TPP Tracy Power Plant
«
USGS United States Geological Survey
1x
-------
-------
SECTION 1
INTRODUCTION
1.1 Development of the Complex Terrain Dispersion Model (CTDM)
CTDM is a model designed to estimate ground level concentrations on
elevated terrain during periods in which the atmosphere is stably
stratified. The model provides concentration estimates for receptors on a
single isolated hill for a single averaging period. The model can accept
multiple terrain features; however, the flow is only influenced by one
hill at a time.
The central feature of CTDM is its use of a critical
dividing-streamline height (Hg) to separate the flow into two discrete
layers. This basic concept was suggested .by theoretical arguments of
Drazin (1961) and Sheppard (1956) and was demonstrated through laboratory
experiments by Riley et al. (1976), Brighton (1978), Hunt and Snyder
(1980), Snyder et al. (1980) and Snyder and Hunt (1984). The flow below
He is restricted to lie in a nearly horizontal plane, allowing little
motion in the vertical. Consequently, plume material below Hg travels
along and around the terrain, rather than up and over the terrain. The
flow above Hg is allowed to rise up and over the terrain. Two separate
components of CTDM compute ground-level concentrations resulting from
material in each of these flows.
An important step in the calculation of concentrations at receptors
above Hg is the determination of lateral and vertical streamline
displacements. The calculation of these displacements for a hill of
arbitrary shape would require the use of an elaborate numerical model and
significant computing resources, neither of which can be justified on the
basis of increased accuracy of the concentration predictions. The current
version of the model is designed to run on a microcomputer.
If one assumes that the portion of the hill above He can be fit to
a simple mathematical surface, then the lateral and vertical streamline
displacements can be estimated from analytical expressions which can be
rapidly evaluated. For CTDM to make use of this idealized terrain, the
model must have access to hill fit parameters for a range of He values.
1.2 Requirements for a CTDM Terrain Preprocessor
CTDM requires much more information about hills than other screening
models. CTDH needs a 3-dimensional representation of each hill.
Therefore, the Terrain Preprocessor produces an analytical description of
-------
the hill shape. Although CTDM will accept several distinct hills, the
Terrain Preprocessor will process only one hill at a time. Hence, the
.Terrain Preprocessor must be run for each hill and the resulting files may
be appended to one large terrain file for input to CTDM. One constraint
of CTDM is that in the calculations.only one isolated hill is considered
at a time. A discussion whose purpose is to aid the user in selecting
distinct terrain features is included in Appendix A.
Since CTDM is designed for regulatory applications, an objective
method is needed to characterize actual terrain in terms of a mathematical
shape. In the absence of such a method, two users analyzing the same hill
(with the same contours) could arrive at significantly different
representations for the fitted hill. The preprocessor provides a display
of the actual and fitted hill to enable the user to determine whether the
fit is reasonable from a physical standpoint, or whether a subfeature of
the digitized terrain should be isolated for analysis.
1.3 Summary of Preprocessor Operation
Two programs must be run to generate terrain input parameters to CTDM
for a given hill. A third program allows the user to display the contours
for the actual and fitted hills. The first program, FITCON, asks the user
to define a hill in terms of its name, identification number, maximum
elevation and x,y coordinates of the hill center. The user then specifies
the name of a master file of digitized contour data and a file to be used
for diagnostic output during the fitting process. In the master file, the
following data is provided for each contour:
• Contour identification number,
• Contour elevation,
• Humber of digitized points,
• A code indicating whether a contour is input as complete or
incomplete,
• x,y coordinates of the digitized contour points.
The user chooses one of the following 3 methods for selection of contours
from the master file: (1) all contours selected, (2) contours selected
based upon a range of user-specified contour identification numbers, or
(3) the specification of a file containing the contour identification
numbers for the hill in question. Before a contour is accepted fur
processing, it must pass a number of tests. •Incomplete contours are
closed by a reflection of points through the hill center or contour
centroid. The program provides special processing for those contours
which are found to be a series of multiple contours at the same
elevation. After qualification and editing (described in Section 2.1.2),
the area and centroid coordinates of each contour are determined by
numerical integration. Each contour is then fit to an ellipse by first
-------
finding the slope of the line through the centroid in the plane of the
contour.which gives the largest second moment for the area within the
contour. In the determination of this maximum second moment for a contour
to 10* resolution, eighteen lines having equal angular spacing are used.
The line associated with the maximum second moment is assumed to define
the orientation of the minor axis of the ellipse representing the
contour. The lengths of the semi-major and semi-minor axes for this
ellipse are calculated from the analytical expressions for the area and
second moment of an ellipse.
These fitted ellipse parameters for each contour are input to the
second preprocessor program, HCRIT, which determines, for the portion of
the hill above a given critical elevation, the best-fit inverse polynomial
profiles along the hill major and minor axes. The center coordinates of
the fitted hill are calculated as the mean of the ellipse center
coordinates for those contours above a given critical elevation, the
orientation of the fitted hill is calculated as a vector average of the
ellipse orientations, weighted by the. ellipse.eccentricity. The user can
specify the critical elevations to be used by HCRIT in two ways. The
first option is to have each contour elevation, with the exception of the
uppermost, serve as a critical elevation. Alternatively, the user can
specify a number of equally spaced critical elevations between the lowest
and uppermost contour. The lowest critical elevation must be at or below
the lowest stack or tower base elevation for model input. Sometimes it
may be necessary to extrapolate imaginary heights from the hill base down
to below the stack base elevation. In the inverse polynomial fit to the
hill profile, a critical elevation is treated as an effective hill base.
HCRIT provides an input file for CTDH which contains the following
information for each critical elevation:
• Ellipse parameters corresponding to the contour at the critical
elevation* (these parameters are interpolated in the case where
a critical elevation does not correspond to a contour elevation),
• Coordinates of the center of the fitted hill,
• Orientation of the major axis of the fitted hill with respect to
north,
• The length scales and exponents for the inverse polynomial fits
along the hill major and minor axes.
The third preprocessor program, PLOTCOV, uses plot files from FITCOM
and HCRIT to generate the following screen displays which aid in the
evaluation of the hill-fitting process:
* It should be noted that the term "critical elevation" is only used here
to indicate that the same elevations are used in the specification of
ellipse parameters and cutoff hills within CTDtf. CTOH uses these
parameters to determine the characteristics of the ellipse at plume
height if the plume is below the computed critical dividing-streamline
height for a given hour.
-------
• Map of digitized contours either as they were input or after they
have been qualified and edited,
• Map of the digitized contours and their associated fitted ellipses,
• For each critical cutoff elevation, a map showing the digitized
contours and the contours of the fitted hill at elevations
corresponding to the elevations of those digitized contours above the
critical elevation.
1.4 Organization of the Manual
This manual is designed for users requiring different levels of detail
regarding system operation. Users wishing to simply run the Terrain
Preprocessor System should consult Section 4, which gives detailed input
requirements for the system. For each of the inputs, references are made
to those portions of the manual giving more detailed information. These
cross references are also provided for the output items described in
Section 5. Hew users of the system should also read Section 3, which
covers the operation of the system. Finally, those users requiring a more
detailed discussion of the terrain-fitting process should consult
Section 2.
-------
SECTION 2
RULES FOR TERRAIN FITTING
This section describes in detail the rules followed by the Terrain
Preprocessor.System in the fitting of terrain features to mathematical
shapes. Also given in this section are the rules which must be followed by
the system user in the preparation of input data. Throughout this
discussion, a clear distinction is made between the rules which must be
followed by the user and the rules followed by the system during the
process of terrain data qualification,.editing and fitting. These
different types of rules are discussed in the same section because it is
important for the user to understand how decisions made by .the system
depend upon user inputs. The user rules discussed in this section are,
however, restated in the listing of input requirements given in Section 4.
2.1 Fitting of Ellipses to Digitized Contours
For an isolated hill, the selection of contours for the fitting of the
hill is relatively straightforward. For more complex terrain, the user
must decide which features are to be included in the description of a
"hill". .For example, if the contours for 2 adjacent peaks are input to the
program, it trill attempt to fit the peaks with a single inverse polynomial
hill. Although this may be appropriate for peaks very close together in
comparison to the distance to the source, the user may need to fit each of
these terrain features in the absence of the other for some receptors (see
discussion in Appendix A).
The user must first identify, from an examination of topographic maps,
those-features which represent individual hills to be used in the
modeling. Bach of these hills should be assigned a name, identification
number, and a maximum elevation. The user must also specify the x,y
coordinates of the hill center. This center does not have to coincide with
the location of the point of maximum elevation of the hill. Since it is
only used in the completion of incomplete contours (see Section 2.1.2.1),
the hill center should correspond roughly to the mean center of the hill
contours which have been input as complete. All contours in the study area
oust be assigned identification numbers and the correspondence between
hills and contours determined. For each hill, a file of contour
identification numbers must be prepared with one contour identification
number on each line of the file. As mentioned earlier, the same contour
may be assigned to more than one*hill.
2.1.1 Rules for Contour Digitization
Since a given contour may be assigned to more than one hill, the
contour parameters must be placed by the user in a master file, which is
read during each'run of the FITCON program for a given hill. Each hill is
then characterized by the user in terms of a set of contour identification
numbers, which is also input to FITCON.
-------
In the master file prepared by the user, each contour is'described
first by a record giving the contour identification number, the elevation
of the contour, the number of digitized points on the contour, and an
indicator (CFLAG) which specifies whether the contour is being input as
open (CFLAG»0) or closed (CFLAG=1). Following this record are a number of
records, each giving the x,y coordinates of the digitized contour points.
For the convenience of the user, all master file parameters are input in
free format. An example of a digitized contour is given in Figure 1. In
the subsequent fitting of this contour to an ellipse, the FITCOU program
assumes that the contour is a polygon in the horizontal (x,y) plane with
the sides of the polygon formed by straight lines connecting adjacent
points input by the user. A sufficient number of points should be selected
by the user to define the basic shape of the contour. An unnecessarily
large number of contour points trill slow down the process by which a
contour is fitted to an ellipse. A maximum of 1000 digitized points are
allowed for each contour. This number could be increased by changing the
value of the parameter NPCMAX and modifying the appropriate DIMENSION
STATEMENTS in the FITCON main program and associated subroutines.
The contour points for a given contour must be input by the user in a
consecutive order, either clockwise or counter-clockwise. All contour
elevations given in the master file must have the same units and origin as
the hill-top elevation (which is specified interactively by the user during
each run of FITCON) and the stack base elevation (which is input to CTDM).
All contour x,y point coordinates must have the same units and zero
reference point as the hill center x,y coordinate, which is also specified
interactively by the user. Obviously, the same scale should be used for
both x and y. The positive x-axis must point to the true east and the
positive y-axis must point to true north. The receptor locations, which
are input directly to CTDM separately, must be specified according to the
same x,y coordinate system as the digitized contour points. As long as the
consistency requirements mentioned above are met, the user is free to use
any coordinate system for specifying the contour elevation and digitized
point coordinates. For a CTDM run, the user will be required to furnish
the factors needed to convert elevations and distances to meters.
There may be cases in which it is not practical or desirable to
digitize the full length of a contour. In this case, the program FITCOU
will complete the contour according to a procedure described in Section
2.1.2.1. There is the requirement that the digitized incomplete contour
input by the user be continuous from the first to the last point. The
program is not designed to complete contours which have been digitized in a
piecewise fashion. The user may also find it necessary to specify multiple
contours at the same elevation as a single contour for the purpose of
determining an equivalent elliptical contour -at that elevation. The rules
for input of multiple contours at the same elevation are given in Section
2.1.2.2.
2.1.2 Contour Qualification and Editing
The program FITCOU will retrieve data from the master file for
contours specified by the user for the hill in question. Before a contour
is accepted for processing, it must be subjected to a number of tests. On
the basis of these tests, the contour is either accepted or rejected.
Before a contour is accepted, however, it may require editing by the
program, as in the case of an incomplete contour or a contour which
represents multiple contours at the same elevation.
6
-------
3000
2500
2000
> 1500
1000
500
500
1000 1500
XDtetwica
2000 2500
Figure 1. Example digitized contour.
-------
The program first determines whether the.contour elevation is greater
than the hill top elevation input by the user. If so, the contour is
rejected. Contours having less than 3 points or more points than the
maximum allowed are also rejected.
Due to errors in the digitization process, a contour may actually
cross itself. This problem can lead to computed values of the contour area
and second moment having opposite signs. If this is found to be the case,
then the contour is rejected. Even if the signs of the area and second
moment are the same, the problem of the contour crossing itself will be
revealed from the display of input contours generated by program PLOTCON.
If the input of an additional contour from the master file would cause
the maximum number of contours (200) to be exceeded, then the input of
contours is halted and a warning message is written to the diagnostic
output file. The ellipse fitting procedure is then carried out using the
200 contours input up to that point. The maximum number of allowed
contours can be increased by changing the value of the variable HCMAX and
the appropriate array dimensions in the FITCOH main program.
2.1.2.1 Contour Completion
A common situation which arises in the association of terrain contours
with individual hills is shown in Figure 2. Both HILL 1 and HILL 2 are
enclosed by a common contour which eventually closes at a rather large
distance from the centers of both hills. To obtain a mathematical
description of the surface of HILL 1, it is necessary to digitize the 3
component contours associated with the hill. If, however, the lowest of
these contours were digitized over its full extent (not shown in Figure 2),
it would strongly bias the computed parameters (center location,
orientation and shape) for the fitted hill, which would bear little
resemblance to HILL 1 as shown in Figure 2. In this case the user should
only digitize that portion of the contour which is associated with HILL 1,
allowing the program FITCOH to complete the partially digitized contour as
shown in Figure 2. The procedures followed by the program in this contour
completion are described below.
A contour completion code (CFLAG) must be specified in the first
record for each contour in the master file. A complete contour is one for
which the coordinates of the first and last digitized points are
identical. If it is determined that the first and last digitized contour
points have identical coordinates and the contour completion code has a
value of 1, then the contour is accepted for processing as it stands. If
the first and last points are identical and the completion code- has a value
other than 1, then a warning message is written to the diagnostic output
file and the contour is accepted for processing as it stands. If the first
and last points are not identical and the completion code is equal to 1,
then the program FITCOH assumes that the user intended to complete the
contour, but, in fact,, did not. In this case, a warning message is written
to the diagnostic output file and a point, having the same coordinates as
the initial point, is added to the contour. If the addition of this point
causes an exceedance of the allowed maximum number of digitized contour
points, then the contour is rejected. If not, then the contour is accepted
for processing. If the first and last contour points are not identical and
the contour completion code is not equal to 1, then the program FITCOH will
complete the contour by the selective removal and addition of contour
points.
8
-------
EXPLANATION
Contour ttrws M displayed on a map
Digitized points along a map contour
firit or last contour point digitized for the incomplete contour
^Points addtd by tn« FTTCON contour completion algorithm
CanMrofHtU.1
Figure 2. The closing of an incomplete contour as performed by program FITCON.
-------
To minimize the probability that the completed contour will cross
itself, the incomplete contour can be edited by FITCON (if this option is
selected by the user) prior to the addition of points to complete the
contour. The objective of this editing is to obtain a sequence of points
whose order is consistently clockwise or counter-clockwise as viewed from
the hill center location, whose x,y coordinates were specified
interactively by the user. The resulting series of points is a subset of
the original set of points. The parameter used in this editing process is
the filtering angle, which is input interactively by the user during
program FITCON execution. This filtering angle must be no smaller than 1
degree and no larger than 22.5 degrees. Until additional experience has
been gained in the application of FITCOH to a wider range of complex
terrain settings, this filtering angle should be set to 1 degree by the
user. This filtering angle is divided into 360 degrees and the result
rounded to the nearest integer to obtain the value for the total number of
angular sectors to be used in the filtering process. Moving in order from
point to point, the program calculates the heading and distance from the
hill center x,y position to the x,y position of the contour point. The
angular sector which contains this heading is then determined. If another
contour point was previously found to occupy this sector and its distance
to the hill center is smaller than the distance from the current point,
then the current point is discarded. If, in moving from one point to the
next, the change in sector number is more than 1, the program assigns a
"pseudopoinf to each of the intermediate sectors. The distance to each
pseudopoint is determined through simple interpolation by angular sector
between the distances to the current and proceeding, points. If the sector
associated with a pseudopoint corresponds to a sector occupied by a
previously evaluated (and subsequently retained) actual point (or
pseudopoint) located at-a greater distance than the current pseudopoint,
then the previous point (or pseudopoint) is discarded and the current
pseudopoint retained for future comparisons. If its distance to the hill
center is smaller than that of the pseudopoint, the previous point (or
pseudopoint) is retained and the current pseudopoint is discarded. The
pseudopoints are only used for making decisions -regarding the retention or
elimination of actual points during the filtering process. Once the
filtering process has been completed these pseudopoints are discarded. The
point filtering process is illustrated in Figure 3.
•
If two actual points are sufficiently close together, the filtration
process will cause an unwarranted removal of one of the points simply
because the points occupy the same sector. This problem is reduced by
choosing a relatively small filtration angle in comparison to the angular
separation between adjacent contour points as viewed from the hill center.
The selection of a very small filtration angle, however, can needlessly
•low down the filtering process. For this reason, the minimum filtration
angle in the program is currently set at 1 degree.
After the contour points have been filtered, the contour is completed
through the addition of points (see Figure 3).. The locations of these
points are determined by a reflection of existing points through the hill
center or incomplete contour centroid. The first step -in this process is
to determine the angle formed by the lines from the hill center to the
first and last digitized contour points. As shown in Figure 4, the
10
-------
A A A
A A A
EXPLANATION
^-Digitized contour point* retained after angular filtering
RnH or last digitized contour point
Digitized contour points removed by angular filtering
HIM center
»Points added by contour completion
Figure 3. Illustration of point filtering and contour completion.
11
-------
Case 1
Points Input Counter Clockwise
-------
computation of this positive acceptance angle depends upon the sense in
which the contour points have been input (clockwise or counter clockwise).
An incomplete contour is considered to have its points input in a clockwise
sense if the area of this contc .r is found to be positive after it has been
completed through the addition of a point corresponding to the hill
center. Otherwise, the order of input for the incomplete contour points is
assumed to be counter clockwise.
If the angle (shown as ^ F in Figure 4) is found to be less than 90
degrees, then the point reflection is performed using the incomplete
contour centroid (determined without the addition of the hill center
point), rather than the hill center. The rationale for this decision can
be understood by examining Figures 5 through 7. Jhe contour completion
shown in Figure 5 is reasonable from a physical standpoint. In this case
the hill center was used as the reflection point, since the angle formed by
joining the hill center to the first and last contour points is greater
than 90 degrees. In Figure 6, however, the hill center is shown to be
relatively close to one segment of the incomplete contour, giving an
acceptance angle of less than 90 degrees. If, in this case, the hill
center were used as the reflection point, then the completed contour would
have the unrealistic shape shown in Figure 6. The shape of the' completed
contour becomes more realistic when the centroid of the uncompleted contour
is used for reflection instead of the hill center (see Figure 7).
nevertheless, the user should exercise care in the choice of the hill
center (see Appendix A).
Ho matter which of the two points is used in the point reflection, the
rule for the addition of points is the same. Moving from the first
digitized point to the last, a new point is located along a line joining
the contour point and the reflection point (hill center or incomplete
contour centroid) at a distance from the contour point which is equal to
twice the distance from the point to the reflection point. If this new
point falls within the acceptance angle, then the point is retained.
Otherwise, the point is discarded. If the addition of a point will cause
the maximum number of allowable points to be reached, then this point is
sat equal to the initial point, thereby prematurely ending the contour
completion process.
2.1.2.2 Dealing with'Multiple Contours at the Same Elevation
In a complex terrain situation, there may be more than one contour at
the same elevation (see Figure 2). If each of these contours is associated
with a well-defined terrain feature, then contours should be input to the
program in conjunction with separate hills. If this is not the case, the
program will accept the input of more than one contour at the same
elevation. The user must input the points as if they lay on a single
contour. If two separate contours having the same elevation are input from
the master file, the second contour will be discarded. Consider the three
contours shown in Figure 8, which fall inside a single contour. The user
must input the coordinates of the beginning point of each contour twice
(i.e. the contour must be closed by the user). The acceptable input
sequence is shown in Figure 8(a). For the calculation of the area,
centroid coordinates, and second moments for this generalized contour
(actually containing 3 contours), the program renumbers the contour points
as shown in Figure 8(b). The two zero area segments connecting the three
contours effectively transform the three contours into a single contour.
13
-------
I
EXPLANATION
—• Digitized contour points
>-• Additional points generated by the FITCON contour completion
Firat or lut dJgitind contour point
Hill owner M input by the UMT
Figure 5. Contour completion with the hill center used as the reflection point
acceptance angle equals 180 degrees.
14
-------
EXPLANATION
SigKized contour points
-••Additional points generated by reflection through the hill center
first or last digitized contour point
Hill center as input by the user
Figure 6. Contour completion with the hill center used as the reflection point
acceptance angle equals 31.9 degrees'.
15
-------
EXPLANATION
»Digitized contour point*
e—•—• Additional points generated by reflection through the incomplete contour cemroid
S First or test Digitized contour point
HiH center ss input by the user
A Computed incomplete contour centroid
Figure 7. Contour completion with the incomplete contour centroid used as the
reflection point - acceptance angle equals 51.5 degrees.
16
-------
a) User input sequence of points.
b) Modified point sequence used by the program.
Figure 8. Analysis of multiple contours at the same elevation.
17
-------
2.1,3 Calculation of the Area and Cent ro id of Each Digitized Contour
The area, A, of the polygon formed by the Np digitized contour
points (after editing) is given by
N
See Section B.I for the derivation of Equation (1).
The value for A will be positive if the points are input in a
clockwise sense and negative if the points are input in a counter-clockwise
sense. The value reported in the diagnostic output file for the contour
area is the absolute value of A. The contour centroid coordinates X,.,
Yc are given by
N
V2 * (2a)
N
•> t \
(2b)
+yk)/2
The calculated values for Xc and Yc will be the same if the
contour points are input in the clockwise or counter-clockwise sense. See
Section B.2 for a derivation of Equations (2a) and (2b).
2.1.4 Determination of Contour Orientation
The next step in the contour fitting is, to assume that each digitized
contour can be approximated by an ellipse having the same centroid and area
aa the contour. The orientation of this ellipse is determined by first
taking the second moment of the contour area about each of eighteen axes
passing through the centroid of the contour in the plane of the contour
(see Figure 9). The second moment, SMg, of the digitized contour polygon
about a line passing through the contour centroid and making an angle
OB with respect to the positive x-axis is given by:
18
-------
3000
2500
2000
1500
1000
500
m = 8 (70°)
m = 7 (60°)
m = 6 (50°)
m = 5 (40°)
m = 13(120°)
14(130°)
ms16(150°)
^
17(160°)
2500
Figure 9. Axes for the calculation of second moments.
g
o
§
19
-------
where
) cos + (Y- ~ *> sin
- -(xk-Xc) sin 6,,! + (yfc-Yc) cos 6m
m '» 1,18 .
For contour points input in a clockwise sense, the value of SH^ will
be positive. It will be negative for points input in a counter clockwise
sense. In determining the axis having the greatest second moment, the
parameter actually used by the program is the -radius of gyration", Rgm,
given by
1/2
(4)
ffl
The value R-JJ will be positive no matter which direction the contour
points are input. The axis having the greatest value of Rgm, R». is
assumed to define the orientation, 6, of the minor axis of the ellipse.
See Section B.3 for a derivation of equations (3) and (4).
2.1.5 Calculation of the Semi-Major and Semi-Minor Axis Lengths for
the Contour Elliptical Representation
In determining the length, a, of the semi-major axis of the fitted
ellipse, the following property of an ellipse is used:
a • 2Rg . (5)
The length, b, of the semi-minor axis of the ellipse is calculated from the
formula for the area of an ellipse as follows:
. (6)
«a
where abs(A) » absolute value of A .
Once a and b are known, then the eccentricity of the ellipse, ECC, can be
calculated as follows:
2 .2 1/2
ECC - <* - b > . ' (7)
If the range of values of R^ for m * 1,18 is less than 1 percent of the
maximum value of R., then the contour is assumed to be circular with a
radius, r, given by
r- (**U)
20
-------
The same assumption of a circular contour is also made if, for some reason,
the computed value for b is found to be greater than a. In the tests of
the program to date, this condition (b > a) has not been observed.
2.2 Mathematical Representation of a Cut-Off Hill
The first preprocessor program, FITCOtI, determines an elliptical
representation for each digitized contour associated with a given hill.
Each contour can then be specified in terms of an elevation, centroid
coordinates, semi-major and semi-minor axis lengths, and the orientation of
the minor axis with respect to the positive x-axis. To calculate
concentrations at hill receptors under stable conditions, the CTOM requires
the following information for each of a series of critical cutoff
elevations:
• The characteristics of the ellipse corresponding to the cut-off
elevation. If this cutoff elevation coincides with a contour
elevation, then the ellipse parameters are those determined by
FITCON for the contour. Otherwise, the ellipse parameters are
obtained by interpolation using parameters for the contours above
and below the critical elevation.
• The center x,y coordinates and the orientation of the fitted hill
corresponding to the portion of the actual hill above the
critical elevation.
• The parameters which give the best inverse polynomial fit to the
cut-off hill profile along the cutoff hill major and minor axes.
These parameters are calculated by HCRIT, the second of the 3 preprocessor
programs, and are passed in a file to CTDtt. The assumptions used in the
calculation of these parameters are described in this section.
2.2.1 Best-Fit Ellipse at a Critical Elevation
The file written by program FITCOtl for input to program HCRIT contains
the following parameters for each digitized contour processed:
e Contour elevation
e x,y coordinates of the fitted ellipse centroid
e Lengths of the semi-major and semi-minor axes for the fitted
ellipse
*
e Eccentricity of the fitted ellipse
e Orientation of the fitted ellipse semi-minor axis with respect to
the positive x-axis
21
-------
The same assumption of a circular contour is also made if, for some reason,
the computed value for b is found to be greater than a. In the tests of
the program to date, this condition (b > a) has not been observed.
2.2 Mathematical Representation of a Cut-Off Hill
The first preprocessor program, FITCON, determines an elliptical
representation for each digitized contour associated with a given hill.
Each contour can then be specified in terms of an elevation, centroid
coordinates, semi-major and semi-minor axis lengths, and the orientation of
the minor axis with respect to the positive x-axis. To calculate
concentrations at hill receptors under stable conditions, the CTDH requires
the following information for each of a series of critical cutoff
elevations:
• The characteristics of the ellipse corresponding to the cut-off
elevation. If this cutoff elevation coincides with a contour
elevation, then the ellipse parameters are those determined by
FITCON for the contour. Otherwise, the ellipse parameters are
obtained by interpolation using parameters for the contours above
and below the critical elevation.
• The center x,y coordinates and the orientation of the fitted hill
corresponding to the portion of the actual hill above the
critical elevation.
• The parameters which give the best inverse polynomial fit to the
cut-off hill profile along the cutoff hill major and minor axes.
These parameters are calculated by HCRIT, the second of the 3 preprocessor
programs, and are passed in a file to CTDM. The assumptions used in the
calculation of these parameters are described in this section.
2.2.1 Best-Fit Ellipse at a Critical Elevation
The file written by program FITCOtr for input to program HCRIT contains
the following parameters for each digitized contour processed:
• Contour elevation
• x,y coordinates of the fitted ellipse centroid
• Lengths of the semi-major and semi-minor axes for the fitted
ellipse
*
• Eccentricity of the fitted ellipse
• Orientation of the fitted ellipse semi-minor axis with respect to
the positive x-axis
21
-------
If there is only one contour, then the values for the ellipse centroid
coordinates are set equal to the values for the single contour. The
semi-major and semi-minor axis lengths for the critical elevation contour
are extrapolated by assuming a zero-area contour at the hill-top
elevation. Also, if the extrapolated value for the major axis is less than
the extrapolated value for the minor axis, then both axes are set equal to
the square root of their product. Finally, if the interpolated or
extrapolated axis lengths are less than the corresponding axis lengths for
the first contour above this critical elevation, then the axis lengths are
set equal to the corresponding values for the first contour above this
critical elevation.
The orientation of the ellipse at the critical elevation, 6e, is
computed by vector interpolation of contour orientations weighted by the
contour eccentricities. The equations for the interpolation
(extrapolation) of the value for the ellipse orientation at the critical
elevation are given below:
6e » tan"1 (SUMY/SUMX) (10)
where
SUM* - ECCT cos 6T * "- (BCCU cos 9U - ECCT cos 9.) (ll(a))
u L V K-, — B. / tl H L> L
/•••L \
( _ j! j
SUMY - BCC, sin 0. * *- (ECC__ sin Ou - ECC. sin 9.) (1Kb))
U L \ B-. — K. / H It L L
eccentricity at EL
ECCH • eccentricity at Eg
&L • orientation at EL
OH " orientation at EH .
If the critical elevation, E, is lower than the lowest contour elevation,
then Equations (10) through (11) still apply except that EL and EH are
defined as follows:
EL » elevation of the lower digitized contour
EH * elevation of the contour immediately above EL
If there is only, one contour, then the value for the ellipse orientation at
the critical elevation is set equal to the orientation for the single
contour. The parameters for the ellipse at the critical elevation are
written by program HCRIT to a file which is input directly to CTDM. These
ellipse parameter records are written to the file in the order of
increasing critical elevation. Before they are written to this file, the
ellipse orientations are modified so that they are expressed in terms of
degrees clockwise from north for the contour major axis.
23
-------
2.2.2 The Inverse Polynomial Profile
In addition to the ellipse parameters at the critical elevation, the
CTDH requires that, for each critical elevation, the hill be fit along its
major and minor axes in the vertical plane with an inverse polynomial
profile. In this fit, the critical elevation is considered as the base of
the cut-off hill. The location, shape and orientation of this fitted hill
depend upon the ellipse parameters for the contours above the critical
elevation. In the following discussion, the procedure for fitting an
inverse polynomial profile to a cutoff hill will be described. First,
however, the mathematical properties of the inverse polynomial profile will
be investigated.
2.2.2.1 Mathematical Properties
A hill having an inverse polynomial profile along the x-direction can
be characterized by the following expression:
(12)
where
h » elevation of the hill surface at a horizontal distance x from
the hill center along the x-direction
hg * base elevation of the hill which would correspond to a
particular critical cut-off elevation
AH - ox - h0
h-j a hill top elevation
t» » profile exponent
L • half-height length scale
An example of this profile is shown in Figure 10 for a range of exponent
values.
To determine the best fit values of P and L for a series of actual
terrain points along the x-direction, it is necessary to express Equation
(12) in the following form:
.lnx-lnL-±ln (. . - 1 (13)
Consider a collection of Hc terrain points with coordinates (Xj. hj)
with h<, < hj < bj. and Xj > 0. The best-fit values for P and L can
be obtained by minimizing the following expression, x2t with respect to
1/P and In L:
24
-------
0
o
X
-s 1
-s
n
o
3
§
o
I
o
u
o
o
u
fr
_ O
O
D
+
o
a.
O
O
u
o
IO
o
s
O O
10 O
25
-------
2
J = 1
in x - in L - in
- 1
J o
(14)
(15a)
)
(15b)
Solving Equations 15 (a) and (b>, the following results are obtained for P
and L:
H *SUM3 - SUM1**2
B *SUM4 - SUM1*SUM2
c
L.
/SUM2*SUM3 - SUM1*SUMA
*SUH3 - SUM1**2
)
*fhere
* - multiplication
** • exponentiation
SUM1
SUM2 .£ In x.
J - 1 J
SUM3
SUM4 . I in x. In
J - 1
26
-------
2.2.2.2 Procedure for Fitting the Profile to a Cutoff Hill
The center coordinates of the fitted hill (XH, YH) are calculated
as the arithmetic mean of the centroid coordinates (Xcj , Ycj) of the
Hc individual contours above the critical cutoff elevation:
(18b)
Ho weighting with respect to contour area is performed since this
would minimize the influence of the higher elevation contours in the
determination of the fitted hill centroid. The orientation angle, 6, for
the minor axis of the fitted hill is assumed to be a vector average of the
individual contour minor axis orientations, weighted by the contour
eccentricities ;
•c
I ECC, sin 6
_ i J
p
•c
I ECC. cos 6
_
6 - tan - - | (19)
•c
where
6j * orientation angle for the ellipse representing contour j
ECCj - eccentricity for the ellipse representing contour j
The eccentricity is used as a weighting factor since the orientation of a
nearly circular contour (ECC-0) should not be considered when determining
the orientation of the fitted cutoff hill.
At this point, all the ellipses above the cutoff elevation have been
assigned a common center (X|{, Yg) and a common orientation (6) of the
semi-minor axis with respect to the positive. x-axis. What remains to be
determined are the best-fit inverse polynomial parameters (P and L) for the
hill major and minor axes. The best-fib parameters for the cutoff hill
major axis (Pa and La) are calculated using Equations (16) and (17)
with Xj • aj (the semi-major axis length for the ellipse representing
contour j). The number of terms in the summation is equal to He, the
number of contours above the critical cutoff elevation. The best-fit
parameters for the cutoff hill minor axis (P^ and LD) are calculated
using Equations (16) and (17) with Xj * bj (the semi-minor axis length
for the ellipse representing contour j).
27
-------
If the best-fit values of either Pa and PJ, ace found to be
negative (implying an increase of contour axis length with height), then
this negative value is set equal to its absolute value.
Once the fitted hill parameters have been determined, they are written
by program HCRIT to the CTDM input file with the parameter records in the
order of increasing hill cut-off elevation. The cut-off hill orientations,
written to the CTDM input file, are expressed in terms of degrees clockwise
from north for the fitted hill major axis.
It should be noted that CTDM actually fits the inverse polynomial
profiles along the major and minor axes of the cut-off hill to the
following Gaussian shapes:
(20a)
h « AH exp (=~) + hQ . (20b)
where
oa » standard deviation parameter for the Gaussian terrain .
distribution shape along the major axis,
«b - standard deviation parameter for the Gaussian terrain
distribution shape along the minor axis.
The assumption is then made within CTDM that the hill terrain shape in 2
dimensions is given by the following expression:
(21)
where
x « distance along the major axis,
y » distance along the minor axis.
This 2-dimensional Gaussian hill will have elliptical contours and will
have a Gaussian profile along an arbitrary direction, 6, with respect to
the major axis with the standard deviation parameter, «r, given by
22 22
(o * sin 6 +
-------
SECTION 3
SYSTEM OPERATION
The CTDM Terrain Preprocessor described in this manual is designed to
run on IBM-PC compatible microcomputers for which a FORTRAN compiler and
BASICA interpreter (or compiler) are available. It is recommended that
the microcomputer contain the appropriate math coprocessor for faster
program execution.
The Terrain Preprocessor actually consists of 3 programs which are
executed consecutively. The first program, written in FORTRAN, is named
FITCON. It determines the best-fit ellipse for each input digitized
contour. These ellipse parameters are used.by the second FORTRAN program,
HCRIT, to calculate, for each of a number of elevations, the ellipse
parameters at that elevation and the best fit inverse polynomial hill
shape for the portion of the hill above that elevation. Since both FITCON
and HCRIT are written in FORTRAN 77, their use is not restricted to
microcomputers. These programs can be run on any mini or mainframe
computer for which a FORTRAN 77 compiler is available. The computer
system must also allow for interactive input and output. The third
program, PLOTCON, displays the input digitized contours, the fitted
ellipse for each contour, and the contours of the fitted inverse
polynomial hills. The operational characteristics of FITCON, HCRIT, and
PLOTCON are described in the following 3 sections. Program listings are -
given in Appendix D of this manual.
3.1 FITCON
The program FITCON is written in the FORTRAN 77 language. The $LARGE
metacommand must precede the FITCON main program- if the Microsoft1"
compiler is used to compile the program. The comments section of the main
program is extensive. The first portion of this section is devoted to a
summary of program operation. This is followed by a detailed alphabetical
glossary of every variable used in the program. The same level of detail
of commenting is followed in the FITCON subroutines.
The main program comments section is followed by a number of TYPE,
DIMENSION and DATA statements. Following these statements, several
program constants are set. These include file unit numbers and maximum
array dimensions. • . .
The program first asks the user to specify-the name of the master
file containing the digitized contour information for the hill in
question. The format for this file, is given in Section 4 tjf this
manual. For each digitized contour in the master file, there is a record
giving the contour identification number, contour elevation, number of
digitized contour points, and a variable which indicates whether the
contour is to be considered open or closed. The contour elevations must
have the same units as the other elevations input to CTDM. This record is
29
-------
followed by the x,y coordinates of each digitized point on the contour.
These coordinates must have the same units as the source and receptor
coordinates used in CTDM. This master file of digitized contours raay be
generated by an automated contour digitization procedure or prepared by
use of the DOS editor EDLIN or some other text editing program.
The user is then asked to input the name of the file to be used for
program diagnostic output. If the name CON is specified for the output
file, then the output is routed directly to the terminal. If the name PRN
is specified for the output file, then the output will be routed directly
to the line printer. If the name NUL is used for the output file, then no
output file will be generated. By specifying a name for the file such as
OUTPUT, the output file can be examined following program execution by use
of the TYPE or PRINT commands or by-inputting the output file to a file
editing program. A detailed description of this output file is given in
Section 5 of this manual. The file contains the following information:
• An echo of inputs specified by the user;
• Listings of input contour data;
• Listings of modified contour point coordinates;
• Messages generated during the contour editing and qualification
process;
• The area, centroid coordinates, orientation and semi-axes
lengths for the elliptical representation of each input
digitized contour.
The user must then specify an identification number and a name for
the hill in question. These two parameters are used.by the plot program
to ensure that plot files generated by programs FITCON and HCRIT
correspond to the same hill. These identification parameters are also
used by the CTDM itself. The user is also asked for the hill-top
elevation and x.y coordinates of the hill center. The hill-top elevation
is passed to program HCRIT to be used in the calculation of best-fit
inverse polynomial parameters for cut-off hills. This elevation must have
the same units as other elevations input to CTDM. The hill center x,y
coordinates are used by FITCON only for the completion of incomplete
contours. They must have the same units as the contour point coordinates
input from the master file. Even if contours are input as complete,
however, these coordinates should have realistic values since the plot of
input digitized contours will be scaled to include this hill center
location. The user is then asked whether poi,nt filtering is to be used
prior to the completion of an incomplete contour. The point filtering
process was described in detail in Section 2.1.2.1. If the answer is yes,
the user is asked to input a filtering angle between 1 and 22.5 degrees.
A value of 1 degree is currently recommended.
There are 3 ways provided by FITCON for selection of hill contours
from the contour master file:
30
-------
• All contours selected from the file,
• All contours having identification numbers falling within a user
specified range,
• A file containing contour identification numbers for the hill in
question.
If the third of these options is selected by the user, then the file
of hill contour identification numbers is sorted using subroutine ISORT to
facilitate the retrieval of contour data from the master file.
The user is then asked Whether a plot containing the digitized and
fitted contours is to be generated. If so, then the user must specify a
name for the plot file. It must have a different name than the plot file
specified by the user during the subsequent execution of the program
HCRIT. After opening this plot file, FITCON also opens a scratch file
which will temporarily hold some of the plot parameters. In the first
record of the plot file, the character expression "FITCOM" is written to
indicate that the plot file is being written by program FITCON. The hill
identification number and hill name are then written to the second record
of the plot file. The hill center coordinates are written to the third
plot file record.
FITCON then begins reading the digitized contour parameters from the
master file and subjecting these contours to a qualification and editing
process. The first qualification step is to determine whether the contour
identification number matches the user-specified selection criteria. If
not, the x,y coordinates for the contour are bypassed using subroutine
SKIPCN and the next contour is input from the master file. If the contour
identification number does meet the selection criteria, the contour is
then subjected to the following additional tests:
• The contour elevation must be less than the hill-top elevation,
• The contour must not have an.elevation equal to that of a
contour which has been previously 'accepted from the master file,
•
• The contour must have at least 3 points,
• The contour must have no more than the maximum allowed number of
digitized points (currently set at 1000).
If the contour passes the tests listed above, the program FITCON
subroutine MULTC determines whether the contour is actually a set of
multiple contours at the same elevation. If so, subroutine MULTC carries
out the point reassignment procedure described in Section 2.1.2.2. The
contour is rejected during this process if the number of maximum points is
exceeded or if the last in.a series of multiple contours is found not be
be closed. If the contour is found to be a single contour (i.e. no
contour closure before the final contour point), then this contour is
subjected to the contour completion process described in Section 2.1.2.1.
31
-------
This procedure is performed within subroutine CONCOMP and its associated
subroutine VECTOR. Before the contour is finally accepted for processing,
however, it must pass the following final -*o tests:
• The absolute value of the contour area obtained by numerical
integration, using subroutine ARCH, must be greater than zero.
• The area and maximum second moment (calculated by subroutine
SMOMNT) for the contour must have the same sign so that a real
value for the radius of gyration may be calculated. A contour
can fail this test if it crosses itself.
Once the contour has-been accepted and if a plot has been requested,
then both the unedited and edited contour point x,y coordinates are
written to the plot scratch file. The plot boundaries for the display of
unedited and edited contours are then updated to reflect the extent of
this newly accepted contour. The centroid coordinates, orientation and
serai-axes lengths for the elliptical representation of this edited contour
are calculated, according to the methods described in Sections 2.1.3
through 2.1.5, and stored in the appropriate arrays. The next contour is
then read from the master file and the qualification, editing and fitting
process is repeated.
The input of digitized contours is halted when the end of the master
file is reached or when the input of an additional contour would exceed
the maximum number of contours allowed for a hill (currently set at 200).
If the third contour selection mode was used, the user is informed,
through a message written to the diagnostic output file, whether or not
all the contours- requested were found in the master file. If no contours
were selected, an error message is written both to the screen and the
diagnostic output file and program execution is terminated. In addition
to the messages written to the diagnostic output file during the contour
qualification and editing process, the user is kept informed by screen
messages as to the disposition of contours input from the master file.
Once all contours have been processed, the master file is closed and
the following information is written to the plot file:
• The number of contours processed,
e The identification numbers of the processed contours sorted in
ascending order,
• The plot boundary limits for both unedited and edited contours.
The scratch plot file containing the coordinates of the unedited and
edited contours is then rewound and the contour coordinates transferred to
the plot file. Once this is completed, the scratch plot file is closed
and deleted. The ellipse parameters for each processed contour are then
written to the plot file.
The user is then asked for the name of the file which will be input
to the program HCRIT. The following information is then written to this
file:
32
-------
• Hill identification number and name;
• Hill-top elevation;
• Processed contour identification numbers sorted in ascending
order;
• The elevation, centroid coordinates, semi-axes lengths,
eccentricity and orientation of the ellipse associated with each
processed contour.
The contour identification numbers are sorted and passed to HCRIT
only for subsequent consistency checks within the plot program PLOTCON.
These checks will be discussed in Section 3.3.
3.2 HCRIT
The HCRIT program is used to calculate (1) ellipse parameters for a
series of elevations and (2) the inverse polynomial fit parameters for the
portion of the hill above each of these elevations, which is assumed to
act effectively as the base of a cutoff hill. When HCRIT is executed, the
user is first asked to specify the name of the file generated by FITCOtI
for input to HCRIT. The program then asks the user for the name of the
file to be passed to CTOM. The user must then specify whether or not a
plot showing the fitted contours is to be generated for each critical
elevation. If a plot is required, then the user must specify the name of
the plot file, which must have a name different from the plot file
generated by FITCOM.
HCRIT then reads the following information from the file passed from
FITCOH:
• The hill identification number and name;
• Hill-top elevation;
• lumber of fitted contours;
• Sorted identification numbers for the fitted contours;
• Elevations, centroid coordinates, semi-axes lengths,
eccentricities and orientations of the fitted contours.
The following information is then transferred to the HCRIT plot file:
*
• The character expression "HCRIT" to identify the plot file as
one generated by the HCRIT program,
• The hill identification number and name,
•
• Humber of fitted contours,
33
-------
• Sorted identification numbers for fitted' contours,
• Hill-top elevation.
The FITCON file is then closed and the fitted contour parameters are
sorted in ascending order with respect to contour elevation. This pointer
sort is carried out by subroutine PSORTR. These sorted contour elevations
are written to the HCRIT plot file.
Two options are available for specifying critical elevations in the
HCRIT program. The user may specify that each contour level (with the
exception of the uppermost contour) is to be used as a critical
•levation. This mode can be used only if there is more than one contour.
The other option is to ask for up to 20 critical elevations to be evenly
spaced between a user-supplied lower elevation and the uppermost contour
elevation. If the second option is selected and H^c critical elevations
are requested with the first elevation beginning at HC^, the remaining
%C~1 critical elevations are determined as follows:
HCX
where
* elevation of the ith critical elevation
elevation of the uppermost contour of the N fitted contours
2,
The value of HC^ specified by the user must be at least one elevation
unit below hg.
Note: CTDM requires that the lowest critical elevation be at or below
stack base elevation. If there are no actual terrain contours near this
elevation, the HCRIT program calculates best-fit ellipses and cut-off hill
shapes at the required heights by extrapolation from the lower part of the
hill for which contour input is available.
Once the critical elevations have been specified, then the number of
these elevations is written to the HCRIT plot file. Also, the hill
identification number, number of critical elevations, hill-top elevation
and hill name are written to the CTDM input file. The program then
determines the best-fit ellipse at each critical elevation according to
the methods described in Section 2.2.1. HHC records are written to the
CTDM input file, each containing the critical elevation, ellipse centroid
coordinates, orientation of the ellipse major axis with respect to north,
and the lengths of the ellipse semi-major and semi-minor axes.
Using the methods described in Section 2.2.2.2., the program then
calculates the following parameters for the cut-off hill corresponding to
each critical elevation:
e Critical elevation,
34
-------
• x,y coordinates of the cut-off hill center,
• Orientation of the major axis of the fitted cut-off hill with
respect to north,
• Inverse polynomial best-fit exponents for the major and minor
axes of the fitted cut-off hill,
• Inverse polynomial best-fit length scales for the fitted cut-off
hill major and minor axes.
In fitting the inverse polynomial profile parameters for the cut-off hill,
only contour elevations greater than one elevation unit above the critical
elevations are used in the fit. These parameters are written to both the
HCRIT plot file and the CTOM input file.
3.3 PLOTCON
The PLOTCOH program is used to display the output from the FITCOH and
HCRIT programs. The user must have an Advanced BASIC(BASICA) interpreter
or compiler to run PLOTCOH. A special version of PLOTCON(HPLTCOH) has
been written in HBASIC for users of the Hercules1" Graphics Board with an
IBM PC, XT or AT. PLOTCOH will not run, however, on non-IBM personal
computers with a Hercules1" Graphics Board. The interpretive version of
PLOTCOH is run by typing the command BASICA PLOTCON.
The first statement of PLOTCOH erases the function key display from
the 25th line of the display, making that line available for program use.
The user is then asked to input the name of the plot file generated by
program FITCOH. This file is opened and the first record is read to
determine whether the expression "FITCOH" is present. If not, an error
massage is written to the screen and the user is asked again to input the
file name. If the file is accepted, the following information is input
from the plot file:
• Hill identification number,
• Hill name,
• Hill center coordinates,
• Humber of fitted contours,
• Contour identification numbers,
«
• Maximum and minimum x and y coordinates for the unedited and
edited contours.
The user is then asked whether the plot is to be low resolution (320
points horizontal by 200 points vertical) color or high resolution (640
points horizontal by 200 points vertical) black and white. For the low
resolution color display, the text and actual contours (unedited or
edited) are given in white, the fitted contours in magenta, and the
background in light blue. If the Hercules1" version of PLOTCOH
35
-------
(HPLTCOH) is used, only one type of display (720 points horizontal by 348
points vertical - black and white) is available. For each type of
display, the scale parameters are assigned values which ensure that each
plot, when output to a graphics printer, will have the same vertical and
horizontal scale and will occupy the maximum area when viewed on the
display terminal. For the low resolution (320 x 200) color option, the
vertical and horizontal scales for the plots on the display terminal are
virtually equal. To obtain the same vertical and horizontal scale for the
hardcopy in high resolution (640 x 200), however, the horizontal scale for
the terminal display had to be made slightly larger than the vertical
scale. Plots obtained with the Hercules1" Graphics Board have the same
horizontal and vertical resolution on both the screen display and
printout. Once the type of display has been selected, the user must
specify whether the actual hill contours are to be plotted as unedited or
edited. As discussed in Section 2.1.2.1, the edited contours represent
contours which were input as open, but have been closed, in some cases by
the point reflection process (either with or without angular filtering).
After these inputs have been completed, PLOTCOH sets a number of plot
boundaries, scale factors and colors. The actual contours are then
plotted with special logic designed to handle the case of multiple
contours at the same elevation. Since coordinates for both unedited and
edited contours are contained in the same file, logic is provided for
skipping the unwanted point sets (unedited or edited). Once the actual
contours have been plotted on the screen, the program plots a square (3
points by 3 points) in white corresponding to the user-supplied location
for the hill center. This plot of actual contours and the hill center is
stored in an array for later display. The program then displays the hill
name and the text "INPUT CONTOURS" at the top of the screen. The program
pauses until the user strikes any key. If the "ESC" (escape) key is
pressed, then execution of the program is terminated. During this and
subsequent pauses in program execution, the user has the option of
printing out the plot by holding down the "Shift" key and pressing the
"PrtSc" key*. Hercules1" Graphics Board users must also press the 0
(zero) key following the "Shift PrtSc" command.
PLOTCOH then displays the elliptical representation of each contour,
overlayed on the plot of actual contours. For each fitted contour, the
following parameters are input from the FITCOH plot file: x.y coordinates
of the ellipse centreid, lengths of the ellipse semi-major and semi-minor
axes, and the orientation of the ellipse minor axis with respect to the
positive x-axis. Starting at the end of the semi-major axis and moving
counter-clockwise from the ellipse centroid, ellipse points are generated
*The user must have one of the commonly available 9-pin graphics printers
to exercise this option. Also, for displays with a resolution of 320 x
200 (color) and 640 x 200 (black and white), the user must have run the
program GRAPHICS.COM prior to PLOTCOH execution. Hercules1" Graphics
Board users must run the program HGC.COM before HPLTCOH execution using
the command HGC PRINT.
36
-------
at 3 degree intervals, transformed to the x,y coordinate system, and then
connected with straight lines. After the fitted contours have been
plotted, the hill name and text "FITTED COHTOURS" are written to the top
line of the display. The program execution then pauses until the user
presses any key. If the ESC key is pressed, program execution is
terminated.
Once a key is pressed, the program returns to text mode and the user
is asked whether a display of fitted cutoff hill contours is to be
generated. If the answer is yes, the user is asked to input the name of
the plot file generated by program HCRIT. If the first record in this
file does not contain the character expression "HCRIT", then the user is
asked again to specify the plot file name. If the file is accepted,
PLOTCOH checks whether there is agreement between the FITCOH and HCRIT
plot files with respect to the hill identification number, hill name,
number of contours, and contour identification numbers. If any of these
parameters do not agree, then an error message is given and program
execution is terminated. If agreement is found, then the program returns
to graphics mode and the elevations of the hill top and individual
contours are input from the HCRIT plot file. The program also reads the
number of critical elevations and begins a loop over these elevations.
For each critical elevation, the following parameters are input from the
HCRIT plot file:
• Critical elevation,
• x.y coordinates of the fitted hill centroid,
• Orientation of.the fitted hill minor axis with respect to the
positive x-axis,
• Inverse polynomial parameters (exponent and length scale) for
the fitted hill major and minor axes.
The program then determines the orientation of the cutoff hill major axis
with respect to the positive x-axis and retrieves the background plot of
digitized actual contours, which was previously generated and stored in an
array. The program then plots the contours of the fitted cutoff hill at
those elevations which are at least one elevation unit above the critical
cutoff elevation. Por display of these contours, the following
generalized version of the inverse polynomial profile (Equation 12) is
used:
AH (23)
37
-------
'where
hj a contour elevation
ho a hill cutoff elevation
AH a hT - ho
h-r a hill top elevation
pa» pb * ' best-fit inverse polynomial exponent parameters for the
fitted hill major and minor axes obtained from Equation
(16)
La. LI, a best-fit inverse polynomial length scale parameters for
the fitted hill major and minor axes obtained from
Equation (17)
x*, y' a coordinates of a point on the inverse polynomial contour
at elevation hj with x' measured along the fitted hill
major axis and y' measured along the fitted hill minor
axis.
Equation (23) may be rewritten as follows:
(wit) ?a + (wit) ?b ' l (24)
where
AFIT - L
BFIT
If Pa and Pj, are equal to 2, then Equation (24) represents an
ellipse. For values of Pa and Pj, less than 2,, the contour shape at
the axis points becomes sharper with the contour- acquiring a diamond shape
for Pa » Pb - 1. For Pa and Pj, values less than 1, the contour
acquires a concave appearance. For values of Pa and PJ, larger than 2,
the contour acquires a rectangular shape with rounded corners. Due to the
wide variation in shapes possible with Equation (24), the fitted contour
is not plotted in a parametric fashion (r, 6) as was done in the case of
the fitted ellipses. Instead, 400 evenly spaced x' values along the major
axis from - AFIT to + AFIT are selected and the corresponding y' values
38
-------
determined. Also, 400 evenly spaced y' values along the minor axis from -
BFIT to -I- BFIT are selected and the corresponding x' values determined.
After transformation to the x.y coordinate system, these 800 contour
points are simply plotted and not connected.
After all the fitted hill contours for a given critical cut-off
elevation have been plotted, the program displays the calculated centroid
of the fitted hill as a 3 x 3 array of points. The hill name and critical
elevation value are then displayed on the top line of the screen. The
program execution then pauses until the user presses any key. Program
execution terminates if the ESC key is pressed. Once a key is pressed,
fitted hill contours are plotted for the next highest critical elevation.
This process is repeated until all critical elevations have been analyzed.
It is important to note that with respect to the cut-off hill, the
only parameters used by the CTDM "LIFT" calculation are the cut-off hill
center, orientation and inverse polynomial profile coefficients for the
major and minor axes. The fitted hill contours are only displayed to
determine whether their size and location correspond to the actual contour
locations, Which in many cases they will not. In this regard, the use of
Equation (24) to represent the shape of these contours.is arbitrary. In
fact, CTDM fits the inverse polynomial profile along the major and minor
cutoff hill axes to Gaussian shapes. The resulting Gaussian-shaped hill
trill have contours which are elliptical. In the CTDM "WRAP" calculation,
the individual fitted ellipses constructed by program FITCOM are used.
39
-------
SECTION 4
•
INPUT REQUIREMENTS
In this section, the requirements for interactive and batch inputs
will be described for the programs FITCON, HCRIT and PLOTCON. References
are provided to discussions in Sections 2 and 3. Some of the inputs
required by the programs HCRIT and PLOTCON are given in output files
generated by FITCON and HCRIT. These files will be described in
Section 5, which deals with program output.
4.1 FITCON
The first program of the Terrain Preprocessor System to be run for a
given hill is program FITCON. The following interactive inputs must be
specified by the user when running program FITCON for a given hill:
• Contour master file name (maximum of 14 characters)*.
(See Section 2.1.1, paragraph 1 and Section 3.1, paragraph 3.)
• Diagnostic output file name (maximum of 14 characters)*.
(See Section 3.1, paragraph 4.)
• Hill identification number (1-99).
(See Section 2.1; Section 3.1, paragraph 5; and Appendix A.)
• Hill name (1 to 15 characters).
(See Section 3.1, paragraph 5.)
• Hill-top elevation in user units (maximum of 10 digits including
the decimal point).
(See Section 3.1, paragraph 5 and Section 2.2.2.1.)
• Hill center x-coordinate in user units (maximum of 10 digits
including the decimal point).
(See Section 3.1, paragraph 5 and Section 2.1.2.1.)
• Hill center y-coordinate in user units (maximum of 10 digits
including the decimal point).
(See Section 3.1, paragraph 5 and Section 2.1.2.1.)
• A yes (Y) or no (H) answer as to whether angular filtering is to
be used in the contour completion process.
(See Section 3.1, paragraph 5 and Section 2.1.2.1.)
*A maximum of 2 characters for the drive specifier and 8 characters for
the file name followed by a period and a maximum of 3 characters for the
file name extension.
40
-------
• Angular filter size (1-22.5 degrees) to be used in the contour
completion analysis (maximum of 10 digits including the decimal
point). Required only if the answer to the previous question is
yes. Recommended value is 1. (See Section 3.1, paragraph 5 and
Section 2.1.2.1.)
• Mode for selection of contours from the master file (1-3)
- Mode 1 - all contours in master file selected for the hill.
Mode 2 - range of identification numbers used in the
selection of contour numbers for the hill.
Mode 3 - identification numbers for the hill contours
specified in a separate file.
(See Section 3.1, paragraphs 6 and 7.)
• Lower bound for the contour identification number range for
contour selection (1-9999). Required only if the contour
selection mode is 2. (See Section 3.1, paragraphs 6 and 7.)
• Upper bound for the contour identification range for contour
selection (1-9999). Required only if the contour selection mode
is 2~. (See Section 3.1, paragraphs 6 and 7.)
• Mama of the file holding the identification numbers of the
contours to be selected from the master file (maximum of 11
characters). Required only if the contour selection mode is 3.
(See Section 3.1, paragraphs 6 and 7.)
• A y«s (Y) or no (H) answer as to whether a plot of the actual
and fitted contours is to be generated. (See Section 3.1,
paragraph 8.)
• Bame of.the PITCOH plot file (maximum of 14 characters).
Required only if the answer to the previous question is yes.
(See Section 3.1, paragraph 8.)
• Hame of the file containing fitted contour output to be passed
to program KCRIT (maximum of 14 characters). (See last 2
paragraphs of Section 3.1.)
The first batch file required by FITCOtf is the master file containing
digitized contour parameters. The input format for this file is specified
in Table 1. For more information related to the preparation of this file,
sea Sections 2.1.1 and 2.1.2. A second batch input file is required only
if contour selection Mode 3 is specified. This file contains the contour
identification numbers for the selected contqurs for the hill in
question. These identification numbers (1-9999) must be given one integer
number per record in free format.
4.2 HCRIT
Following the successful completion of the FITCOtf run, program HCRIT
is then run for the hill in question. Aside from the file of fitted
contour parameters passed to it by program FITCOtf, program HCRIT requires
only the user-supplied interactive inputs given below:
41
-------
TABLE 1
FORMAT FOR THE CONTOUR MASTER FILE
Record Parameter
Croup Haute Format Description
1* IDC Free Contour identification number
(integer)
HCON Free Contour elevation (user units)
NPC • Free Number of input points for contour
(integer)
CFLAG Free Contour flag (integer)
0 » Open
1 - closed
2* XCON, Free x, y coordinates of the NPC contour
YCOM** points (user units)
*Record group repeated for each contour
**Values may span more than one record. A convenient format for file review
and editing is to specify x and y for a single contour point on one record.
-------
• Name of the file containing fitted contour parameters generated
by program FITCON (maximum of 14 characters) must be the same as
the last file name input to program FITCON. Also see Section
. 3.2, paragraphs 1 and 2.)
• Name of the output file to be passed to CTDH (maximum of 14
characters). (See Section 3.2, paragraphs' 1 and 6.)
• A yes (Y) or no (N) answer as to whether plots of fitted hill
contours are to be generated for each critical elevation.- (See
Section 3.2, paragraph 1.)
• HCRIT plot file name (maximum of 14 characters). Required only
if the answer to the previous question is yes. (See Section
3.2, paragraph 1.)
• Selection mode for critical elevations (1-2).
Mode 1 - critical elevations at all contour elevations
except uppermost
- Mode 2 - critical elevations at all evenly spaced levels
between a user-supplied elevation and the uppermost contour
elevation (see Section 3.2, paragraphs 4 through 6.)
• Number of critical elevations (1-200). Required only if
critical elevation selection mode 2 is used. (See Section 3.2,
paragraphs 4 through 6.)
• Lowest critical elevation (maximum of 10 digits including the
decimal point). Must be at least one elevation unit below the
highest contour elevation. CTDM requires this elevation to be
. at or below stack base elevation. (See Section 3.2, paragraphs
4 through 6.)
4.3 PLOTCON
With the exception of the plot files generated by FITCON and HCRIT.
the only inputs to PLOTCON are specified interactively by the user. These
input parameters are listed below:
• Name of the plot file from program FITCON (maximum of 14
characters) oust be the same file name as that specified
interactively during program FITCON execution. (See
Section 4.1.)
• Display selection option (1 or 2).
Option 1 - low resolution with color
Option 2 - high' resolution black and white. Not required
if the Hercules1" version of the program is used. (See
Section 3.3, paragraph 3.)
• Contour type selection option (1 or 2)
Option 1 - Unedited contours
Option 2 - Edited contours.
(See Section 3.3, paragraph 3.)
43
-------
A yes (Y) or no (N) answer as to whether a plot of fitted
contours is to be generated for each critical elevation.
(See Section 3.3, paragr ?h 6.)
Name of the plot file from program HCRIT (maximum of 14
characters). Required only if the answer to the previous
question is yes. Must be the same file name as that specified
interactively during program HCRIT execution. (See Section 4.2.)
-------
SECTION 5
OUTPUT DESCRIPTION
In this section, the output of the 3 terrain preprocessor programs
FITCOH, HCRIT and PLOTCON is described. There are 4 types of output
generated by the 3 programs which constitute the Terrain Preprocessor:
• Screen diagnostics,
• Diagnostic information written to an output file or printer,
• Files passed from one program to the next,
• Graphical output.
Program FITCON utilizes the first 3 types of output, while HCRIT only
passes files to PLOTCON and CTDM. Program PLOTCON generates both
graphical output and screen diagnostics.
5.1 FITCOH
The following screen diagnostic outputs are provided by program
FITCOH:
• A message that a particular contour (identified by its ID
number) has been rejected. User is directed to consult the
diagnostic output file after program completion. The diagnostic
output file is described below.
• A message that a particular contour (identified by its ID
number) has been accepted.
• A message that no contours were requested (i.e. the contour ID
file is empty). Applies only if contour selection mode 3 is
used.
• A message that no contours were selected from the master file.
The following information is written by program FITCOH to the
diagnostic output file:
e Hill identification number and name. (See Section 3.1,
paragraph 5.)
e Angular filter size specified by the user (only if angular
filtering has been requested by the user). (See Section
2.1.2.1.)
e Angular filter size after it has been modified to divide evenly
into 360 degrees (only if angular filtering has been requested
by the user). (See Section 2.1.2.1.)
45
-------
• Hill-top elevation and hill center coordinates. (See Section
3.1, paragraph 5.)
• A message that all contours in the master file are to be
selected (only if contour selection mode 1 is used). (See
Section 3.1, paragraphs 6 and 7.)
• Lower and upper identification number limits for selection of
contours from the master file (only if contour selection mode 2
is used). (See Section 3.1, paragraphs 6 and 7.)
• List of sorted identification numbers for contours to be
selected from the master file (only if contour selection mode 3
is used). (See Section 3.1, paragraphs 6 and 7.)
• Contour identification number and elevation. (See Section
2.1.1.)
• x.y coordinates input for a contour. (See Section 2.1.1.)
• Modified x.y coordinates for a contour. (See Section 2.1.2.)
• Contour area and centroid coordinates (See Section 2.1.3.)
• Semi-axis lengths, eccentricity, and orientation for the ellipse
representing a contour. (See Section 2.1.4.)
In addition to the items of information listed above, the following
messages are written by FITCOK to the diagnostic output file if errors are
found:
• Maximum number of contours reached. (See Section 2.1.2,
paragraph 4.).
• Contour elevation greater than the hilltop elevation - contour
rejected. (See Section 2.1.2, paragraph 2.)
e Previously accepted contour has the same elevation as the
current contour - current contour rejected. Multiple contours
at the same elevation must be input as a single contour. (See
Section 2.1.2.2.)
e Contour has fewer than 3 points - contour rejected. (See
Section 2.1.2, paragraph 2.)
«
e Contour has more than the maximum number of allowable points -
contour rejected. (See Section 2.1.1, paragraph 2.)
• Maximum number of points exceeded in the contour point
reassignment process for multiple contours at the same elevation
- contour rejected. (See Section 2.1.2.2.)
e The last in a series of multiple contours was found not to be
closed - contour rejected. (See Section 2.1.2.2.)
e Contour found to be a single contour (i.e. no contour closure
was found before the final contour point). (See Section
2.1.2.2.)
46
-------
• . Point reassignment for the multiple contour successfully
completed. (See Section 2.1.2.2.)
• Point reassignment for the multiple contour successfully
completed after the point input order of one or more component
contours was reversed to make the order of point input for each
component contour the same as the first component contour. (See
Section 2.1.2.2.)
• Contour completion halted due to exceedance of the maximum
number of allowable contour points. The final contour point
will have coordinates equal to those of the initial point. (See
Section 2.1.2.1*)
• Contour specified as closed was found to be open. Added final
point assumed to be the same as the initial point. (See Section
2.1.2.1.)
e Maximum number of contour points exceeded in the closing
operation. Final point is replaced by the initial point. (See
Section 2.1.2.1.)
e Contour specified as open was found to be closed. (See Section
2.1.2.1.)
e Contour area found to be zero - contour rejected. (See
Section 2.1.2 paragraph 2.)
e A real value for the radius of gyration could not be computed -
contour rejected. (See Section 2.1.2, paragraph 3.)
• The relative difference between the maximum and minimum radii of
gyration for the contour was found to be less than 1 percent.
Contour assumed to be circular. (See Section 2.1.5.)
e Ho contour identification numbers were found in the contour ID
file (only applies if contour selection mode 2 is used). (See
Section 3.1, paragraphs 6 and 7.)
e Ho contours were selected from the master file.
Program FITCOH constructs 2 output files which are subsequently read by
the HCRIT and PLOTCOH programs. The file passed to the PLOTCOH program
contains information required to plot actual and fitted contours for the
hill in question. The format for this file is given in Table 2. The
second file, which is passed to program HCRIT, contains the parameters for
the fitted ellipses, along with the elevations of. the hill top and
contours. The format for this file is given in Table 3.
5.2 HCRIT
The only output froa program HCRIT consists of two files. The first
is a plot file containing parameters required by PLOTCOH for the display
of contours on cut-off hills for a number of cut-off elevations. The
47
-------
TABLE 2
FORMAT FOR THE PLOT FILE GENERATED BY FITCON
Record
Group
1
2
2
3
3
4
5*
6
7
8**
8
*
10
10
11
12*
12
12
Parameter
Name
-
IDHILL
HNAME
XHTOP
YHTOP
HC
IDC(J),
J-l. HC
XMIN1, XMAX1,
YMIN1, YMAX1
XMIH2, XMAX2,
YMIN2, YMAX2
NPCSV
HCON(J)
XCONSV(K) ,
YCONSV(K) ,
K-l, HPCSV
MFC
HCON(J)
XCON(JC) ,
YCON(K) ,
K-l, HPC
XCM
YCM
A
Columns
1-6
1-2
4-18
1-15
16-30
1-10
1-10
1-60
1-60
1-10
11-25
1-30
1-10
11-25
1-30
1-15
16-30
31-45
Format
A6
12
A15
E15.4
E15.4
no
no
4E15.4
4E15.4
110
E15.4
2E15.4
110
E15.4
2E15.4
«
E15.4
ELS. 4
E15.4
Description
"FITCON"
Hill identification number
Hill name
Hill center x-coordinate
Hill center y-coordinate
Number of contours accepted
Sorted contour ID numbers
Boundary limits for unedited
contours
Boundary limits for edited
contours
Number of contour points for
an unedited coutour
Elevation of contour J
x,y coordinates for an
unedited contour
Number of contour points for
an edited contour
Elevation of contour J
x,y coordinates for an
edited contour
x-coordinate for a contour
eentroid
y-coordinate for a contour
eentroid
Semi-major axis length for
the elliptical representation
of the contour
48
-------
TABLE 2 (Continued)
Record Parameter
Group Maine Columns Format Description
12 B 46-60 E15.4 Semi-minor axis length for
the elliptical representation
of the contour
12 OREN 61-75 E15.4 Orientation (degrees) of the
semi-minor axis of the
ellipse with respect to the
positive x-axis (east)
*Repeat for each contour
**Groups 8 through 11 repeated for each contour
49
-------
TABLE 3
FORMAT FOR THE FILE GENERATED BY FITCON FOR INPUT TO HCRIT
Record
Group
1
1
2
3
4*
4
Parameter
Name
IDHILL
HNAME
HTOP
MC
HCON
XCH
Columns
1-2
4-18
IS
10
1-15
16-30
Format
12
A15
E15.4
110
E15.4
E15.4
Description
Hill identi;
Hill name
Hill-top eli
Number of a<
Contour ele<
x-coordinati
4
4
YCH
B
BCC
OREH
31-45
46-60
61-75
76.-90
91-105
centroid
E15.4 /-coordinate of contour
centroid
E15.4 Ellipse semi-major axis
length
E15.4 Ellipse semi-minor axis
length
E15.4 Eccentricity of the ellipse
E15.4 Orientation (degrees) of the
ellipse semi-minor axis with
respect to the positive
x-axis (east)
*Repeat for each contour
50
-------
format for this file is given in Table 4. The second is a file of
parameters required by CTDM for the calculation of plume transport
parameters for a number of critical elevations. The format for this file
is given in Table 5.
5.3 PLOTCOH
The program PLOTCOH generates output in the form of screen displays
and screen diagnostics. The program first displays actual contours in
either unedited or edited format. This display is then overlayed with a
display of fitted ellipses corresponding to the actual contours. Finally,
a set of inverse polynomial contours is displayed for the fitted hill
corresponding to each of a series of cut-off elevations. The following
screen error mesages are generated by PLOTCOH:
• FITCOH plot file not found,
• HCRIT plot file not found,
• FITCOH and HCRIT hill identification numbers do not match,
• FITCOH and HCRIT hill names do not match,
• FITCOH and HCRIT number of contours do not match,
• FITCOH and HCRIT contour identification numbers do not match.
51
-------
TABLE 4
FC 1AT FOR THE PLOT FILE GENERATED BY HCRIT
Record
Group
1
2
2
3
4*
5
6*
7
g**
8
8
8
8
Parameter
Name
-
IDHILL
HNAME
NC
IDC(J),
J»l. NC
HTOP
HCON
NCR
HC
ZHTOPF.
YHTOPF
ORENF
PA. PB
LA. LB
Columns
1-6
1-2
4-18
10
10
15
15
10
15
16-45
46-60
61-90
91-120
Format
A6
12
A15
110
110
E15.4
E15.4
110
E15.4
2E15.4
E15.4
2E15.4
2E15.4
Description
"HCRIT"
Hill identification number
Hill name
Number of contours
Sorted contour ID numbers
Hill-top elevation
Sorted contour elevations
Number of critical elevations
Critical elevation
x.y coordinates for the
fitted cut-off hill centroid
Orientation (degrees) of the
fitted cut-off hill minor
axis with respect to the
positive x-axis (east)
Inverse polynominal exponent
parameters for the major and
minor fitted hill axes.
Inverse polynominal length
scale parameters for the
major and minor fitted hill
axes
"Repeated for each contour
**Repeated for each critical elevation
52
-------
TABLE 5
FORMAT FOR THE FILE GENERATED BY HCRIT FOR INPUT TO CTDM
Record
Group
1
1
I
I
2*
2
2
2
3*
3
3
3
3
Parameter
Name
IDHILL
NCR
HTOP
HNAME
HC
XCM, YCM
ONOR
A, B
HC
ZHIOPF,
YHTOPF
OHOR**
PA, PB
LA. LB
Columns
6-7
9-10
21-30
31-45
1-10
11-30
31-40
41-60
1-10
10-30
31-40
41-60
61-80
Format
12
12
E10.4
A15
F10.3
2E10.4
F10.3
2F10.3
F10.3
2E10.4
F10.3
2F10.3
2F10.3
Description
Hill identification number
Number of critical elevations
Hill-top elevation
Hill name
Critical elevation
x,/-coordinates of the
ellipse centroid for the
critical elevation
Orientation (degrees) of the
ellipse major axis with
respect to north
Semi-major and semi-minor
axes lengths for the ellipse
at the critical elevation
'Critical elevation '
x,y coodinates for the
fitted cut-off hill centroid
Orientation of the fitted
cut-off hill major axis with
respect to north (degrees)
Inverse polynominal exponent
parameters for the major and
minor fitted hill axes
Inverse polynominal length
.scale parameters for the
major and minor fitted hill
axes
*• There are NCR group 2 records followed by NCR group 3 records
** Same as the variable name used for the third parameter of record group 2,
the value is, however, different.
53
-------
REFERENCES
Brighton, P.W.M. 1978. Strongly Stratified Flow Past Three-Dimensional
Obstacles. Quart. J. R. Meteor. Soe.. 104; 289-307.
Drazin, P.G. 1961. On the Steady Flow of a Fluid of Variable Density Past
an Obstacle. Tellus 13: 239-251.
Holzworth, G.C. 1980. The EPA Program for Dispersion Model Development
for Sources in Complex Terrain. Second Joint Conference on
Applications of Air Pollution Meteorology, New Orleans, LA. AMS,
Boston.
Hovind, E.L., M.W. Edelstein, and V.C. Sutherland, 1979. Workshop on
Atmospheric Dispersion Models in Complex Terrain. EPA-600/9-79-041.
U.S. EPA. Research Triangle Park, B.C.
Hunt, J.C.R. and W.H. Snyder 1980. Experiments on Stably and Neutrally
Stratified Flow Over a Model Three-Dimensional Hill. J. Fluid Meeh..
96: 671-704.
Ril*y. J.J., Liu, H.T. and Geller, E.W. 1976. A Numerical and
Experimental Study of Stably Stratified Flow Around Complex Terrain.
EPA Report No. EPA-600/4-76-021, Research Triangle Park., NC, 41p.
Sheppard, P.A. 1956. Airflow Over Mountains. Quart. J.R. Meteor. Soe..
82: 528-529.
Snyder. W.H., R.E. Britter and J.C.R. Hunt 1980. A Fluid Modeling Study
of the Flow Structure and Plume Impingement on a Three-Dimensional
Hill in Stably Stratified Flow. Proc. Fifth Int. Conf. on Wind Engr.
(J.E. Cermak, ed.), 1: 319-329, Pergamon Press, NY, NY.
•
Snyder. W.H. and J.C.R. Hunt 1984. Turbulent Diffusion from a Point
Source in Stratified and Neutral Flows Around a Three Dimensional
Hill; Part II - Laboratory Measurement of Surface Concentrations.
Atmos. Envir.. 18: 1969-2002.
Strimaitis, D.6.. R.J. Paine, B.A. Egan, and R.J. Yamartino, 1987. EPA
Complex Terrain Model Development. Final Report. EPA Contract No.
68-02-3421. Prepared for U.S. EPA, Research Triangle Park, NC.
54
-------
APPENDIX A
SELECTION OF TERRAIN FEATURES
FOR CTOH
55
-------
APPENDIX A
SELECTION OF TERRAIN FEATURES FOR CTDM
The CTDM user is required to assign a terrain element (hill number)
to each model receptor near a terrain feature. (For receptors away from
terrain features (hill number 0), a flat terrain calculation is performed.)
CTDM considers flow around isolated hills and is not designed to
account for complex interactions among several terrain features adjacent
to one another. Therefore, CTDM operates under the assumption that the
meteorological input represents local flows with the steering of very
large terrain features already accounted for. It is recommended that the
most local terrain element be defined as a hill for each receptor.
In areas where hills are not isolated, the user faces a decision as
to which portion of a large hill complex should be chosen for use in
CTDM. In such cases, an incomplete hill (with some contours not closed)
is often the most appropriate choice. Such incomplete hills are typically
portions of a more complex hill that are on the side facing the source in
question.
An example of how an incomplete hill is defined is shown in
Figure A-l. In this figure, the monitor locations labeled 1, 3, 4, and 6
reside on a portion of the much larger Piedmont Mountain complex that also
encompasses monitors 5, 7, 8, and 9. In general, it is- appropriate to
isolate a portion of a hill from a larger complex if the receptor points
in question are at or above a closed contour on the hill portion. This is
the case (although marginally so) for monitors 1, 3, 4, and 6. In this
example, the data points for contours, that are not closed on this hill are
•topped during the digitization process at the hill portion boundary as
defined by the user (e.g., see Figure A-l). The Terrain Preprocessor will
then complete these contours through the use of symmetry about a hill
center specified by the user. In this way, the isolated hills required by
CTDM are created from non-isolated hill input. The local hill shape for
monitor locations 1, 3, 4, and 6 is retained by using the smallest hill
portion possible. The part of the isolated hill that is "made up" from
the non-isolated hill for use in CTDM is generally to the lee of the
monitor locations, and is therefore not a significant factor in the model
calculations.
A sample terrain configuration featuring a complex of two local hills
as part of a larger hill system is shown in Figure A-2. For plumes from
stack A, receptors 1 and 2 can be considered to be on two separate local
hills (divided by the dashed line in the figure). For stack B, however,
the influence of Hill A as well as Hill B must be considered for receptor
2 (the Terrain Preprocessor would consider a multi-peak hill for this
56
-------
01
01
I
o
IM
•v4
a
0
4J
01
4J
w
O .C
4J
a 3
I
•0
i
57
-------
Stack
• Stack A
Figure A-2.
Example illustrating bow to select terrain features for
various CTDM receptor locations. Receptors 1 and 2 should
be associated with local hill* A and B, respectively, for
•tack A (these hills are separated by the dashed line).
For stack B, receptor 2 oust consider hill A and B
together. Receptor 3 can be associated with local hill A
for both stacks. For receptor 4, a combined terrain
feature including both hills oust be used for both stacks.
58
a
I
-------
instance). The complication for receptor 2 makes it necessary to consider
separate CTDH runs for stacks A and B. In order to use receptor 2 for
Stacks A and B in a single run, a user has to create two coincident
receptors for that position and assign different hill numbers for each,
depending on the stack. Afterward, the results for the inappropriate
receptor for each stack must be ignored.
The user must input a hill center for each CTDM hill. This center is
used only in the contour completion algorithm of the terrain preprocessor;
it does not affect closed contours. The hill center does not necessarily
have to be at the .physical center of the hill. There is no program input
restriction on where incomplete contours end, since the terrain
preprocessor will accept incomplete contours and complete them through a
reflection algorithm about the specified hill center. However, it is
advisable to define hill such that contours extend at least halfway around
the specified hill center. The center should be chosen to minimize
crossover of completed ("edited") contours onto any existing closed
contours. The graphical display of the edited contours will aid the user
in determining if the contour completion looks reasonable or if the hill
center position should be adjusted. Fortunately, the arbitrary
positioning of the hill center appears to have a minor effect upon the
maximum concentrations on the hill if an adequate receptor coverage is
provided (Strimaitis et al., 1987).
If a receptor does not appear to be associated with any hill (i.e.,
is in "simple" or flat terrain, with the plume height far above the
terrain features), the user may still model this receptor with CTDM by
specifying hill number '0' in the receptor input file.
59
-------
APPENDIX B
DERIVATION OF EQUATIONS FOR
THE AREA, CENTROID COORDINATES.
AND SECOND MOMENT OF A DIGITIZED CONTOUR
60
-------
APPENDIX B
DERIVATION OF EQUATIONS FOR
THE AREA, CENTROID COORDINATES,
AND SECOND MOMENT OF A DIGITIZED CONTOUR
B.I Derivation of Equation (1) for the Area of a Digitized Cjmtour
An example digitized contour is shown in Figure B-l. The Np
individual digitized points (x^, y^) have been connected by straight
lines. Consider the trapezoid formed by connecting the following points
in order (xk, y^) , (xk+1, yfc+i)t
-------
Yk+1
k+1
Vk
Vk-1
• Trapezoidal Element
Figure B-l. Example Digitized Contour.
62
-------
An evaluation of Equation (B-2) gives
. < , (x
It 2 x. - x 1+1
2 (xi+l*i * Vi+l} (xi+l + V * 3
A summation of the Mx^ terms over k will give the total x-component of
the first moment. Dividing this sum by the contour area (positive or
negative) gives the x-coordinate of the centroid (See Equation (2a)). The
y-coordinate of the centroid is obtained by interchanging x and y in the
formula for the centroid x-coordinate and multiplying by -1 to correct for
the fact that the contour area was computed by integration in the x,y
coordinate system.
B.3 Derivation of Equation (3) for the Second Moment of a Digitized
Contour About a Line Passing Through the Contour Centroid
The calculation of a contour second moment about a line passing
through the contour centroid is performed by calculating the second
moments of individual trapezoidal sections and summing up the
contributions. One such trapezoidal section is shown in Figure B-2. It ~
is constructed by dropping perpendiculars from the points (x^, y^> and
(xfc+i, yfc+i) to the line passing through the contour centroid and
making an angle 6 with respect to the x-axis . The- value of D^ can be
calculated by use of the identity that the magnitude of the cross product
of two vectors is equal to the product of the magnitudes of the vectors
times the sine of the angle between the vectors. Both of these vectors
lie in the plane of the contour. The first is the unit vector along the
line passing through the contour centroid. The second vector is formed by
connecting the centroid and the point (x^, y^) . The cross product of
these vectors is given by
(cos 6 + sin 6 ) • [(x - X + ( - y - (B-5)
where
i, j and Tic are unit vectors in the x,y and z directions.
Equation (B-5) can be rewritten as
- I,.) sin 6 + (jfc - Yc)cos 0 . (B-6)
•ote that Dk will be positive for points lying to the left of the line
passing through the centroid and negative for points lying to the right of
the line.
63
-------
k+1
Figure B-2. Digitized Contour Illustrating the Second Moment
Calculation.
-------
The final side of the trapezoidal section has a length W^ (s?-e
Figure B-2). It connects the foot of the perpendicular D^ to the . oot
of perpendicular D^^. The length of W^, which may be either positive
or negative, is calculated by use of the identity that the scalar product
of two vectors is equal to the product of the magnitude of the two vectors
times the cosine of the angle between the vectors. Both of these vectors
lie in the plane of the contour. The first is the unit vector along the
line passing through the contour centroid. The second vector is formed by
connecting the centroid and the point (x^, y^) . The dot product of
these vectors is the distance from the centroid to the foot of
perpendicular D^ and is given by
(cos 8 + sin 8 ) •
- (xk - X.) cos 8 + (yk - YC) sin 6 . (B-7)
The distance from the centroid to the foot of perpendicular Dfc+i is
given by
- *c)cos 0 + (y - Y)sin 6 . (B-8)
The distance, W^, from the foot of perpendicular Dfc to the foot of
perpendicular Dfc+i is calculated by subtracting Equation (B-7) from
Equation (B-8):
cos 6 -h (y - yk) sin 0 . (B-9)
The second moment, M2^, about the B-axis for the trapezoidal segment
shown in Figure B-3 is given by
°k Vl a-D ' 2. -
M2. - J W-ada + J \L (1-
o
' 5
-------
da
Figure B-3. Trapezoidal Segment for Second Moment Calculation.
• k-1
Contr ibufes a "Negative"
Second Moment
Figure B-4. Hon-Trapezoidal Segments for Contours Crossing the Axis.
66
-------
there can be cases in which D^ and Dfc+i have opposite signs. In this
case, there is no longer a trapezoidal segment, but two right similar
triangles with opposite sides equal to D^ and D^-i-l- If one uses
Equation (B-ll) to calculate the second moment contribution of this
segment, the result will correspond to the difference between the second
moment of the triangle with side D^ and the triangle with side D^^.
As shown in Figure B-4, this negative contribution of the triangle with
side Dfc+i is needed to properly offset a portion of the second moment
associated with the trapezoidal segment having sides D^^ and Dfc+2-
The computation of the total second moment of the contour (about the axis)
obtained through a summation of Equation (B-ll) over k is therefore exact.
B.4 Calculation of the Radius of Gyration of an Ellipse about its Minor
Axis
Consider an ellipse with a semi-minor axis of length, b, and a
semi-major axis of length, a. The equation of this ellipse, with its
major axis along the x-axis and its centroid at the origin, is given by
The second moment, M2, about the minor axis is then given by
M2 - 4j x y dx (B-13)
0
a
4bJ i
o
(B-14))
4
The radius of gyration, Eg, is.therefore,given by
Rg ' (area of^he ellipse) (B~15)
1/2
f • (B-16)
67
-------
APPENDIX C
PROGRAM TEST CASE
68
-------
APPENDIX C
PROGRAM TEST CASE
For the purpose of illustrating the operation of the Terrain
Preprocessor, a single "gullied hill" having 7 closed contours is used.
This example does not exercise the contour completion or multiple contour
(at a single elevation) processing capabilities of program FITCON. These
were tested, however, by the cases run for the CTDM model evaluation
studies.
The file names chosen for this test case (DATCOM, OUTPUT, PLOT1,
COHOUT, PLOT2 and CTDMIH) are arbitrary. Different names could be selected
by the user. To obtain higher resolution plots for this report, the
Hercules™ version of PLOTCON (HPLTCON) has been used.
69
-------
COHTOUR MASTER FILE (DATCON)
70
-------
1 .1000E+02 35
.OOOOE+00 .4000E+01
.OOOOE+00 .2400E+02
.1000E+01 .2800E+02
.2000E+01 .2900E+02
.4000E+01 .3000E+02
.5000E+01 .3000E+02
.8000E+01 .2800E+02
.1000E+02 .2400E+02
.1300E+02 .2300E+02
.1500E+02 .2400E+02
.1700E+02 .2700E+02
.1900E+02 .2900E+02
.2200E+02 .3000E+02
.2400E+02 .3000E+02
.2600E+02 .2900E+02
.2800E+02 .2700E+02
.3000E+02 .2300E+02
.2900E+02 .2200E+02
.2600E+02 .2100E+02
.2100E+02 .1900E+02
.2000E+02 .1700E+02
.2000E+02 .1600E+02
.2100E+02 .1400E-K02
.2200E+02 .1300E+02
.2600E+02 .1000E-K02
.2800E+02 .8000E-K01
.3000E+02 .4000E-K01
.3000E+02 .3000E+01
.2800E+02 .2000E+01
.2400E+02 .2000E+01
.1700E+02 .3000E-H01
.1100E+02 .2000E+01
.5000E+01 .OOOOE+00
.2000E+01 .1000E+01
.OOOOE+00 .4000E+01
2 .2000E+02 32
.1000E+01 .7000E+01
.1000E+01 .2100E+02
.2000E+01 .2500E+02
.3000E+01 .2800E+02
.5000E+01 .2900E+02
.7000E+01 .2700E+02
.9000E+01 .2300E+02
.1200E+02 .2100E+02 -
.1300E+02 .2100E+02
.1600E+02 .2300E+02
.2100E+02 .2800E+02
.2300E+02 .2900E+02
.2700E+02 .2600E+02 .
.2800E+02 .2500E+02
.2800E+02 .2400E+02
.2300E+02 .2400E+02
.2000E+O2 .2300E+02
.1900E+02 .2200E+02
.1800E+02 ' .1800E+02
.1800E+02 .1600E+02
.1900E+02 .1300E+02
.2300E+02 .9000E+01
71
-------
.2600E+02 .6000E+01
.2600E+02 .5000E+01
.2400E+02 .5000E+01
.2100E+02 .6000E+01
.1700E+02 .7000E+01
.1300E+02 .6000E+01
. 1000E-I-02 .5000E+01
.5000E+01 .4000E+01
.3000E+01 .4000E+01
.1000E+01 .7000E+01
3 .3000E+02 27
.2000E+01 .9000E+01
.2000E+01 .2100E+02
.3000E+01 .2500E+02
.5000E+01 .2700E+02
.7000E+01 .2500E+02
.9000E+01 .2100E+02
.1100E+02 .1900E+02
.1300E+02 .1900E+02
.1800E+02 .2400E+02
.2200E+02 .2700E+02
.2500E+02 .2700E+02
.2600E+02 .2500E+02
.2300E+02 .2500E+02
.2000E+02 .2400E+02
.1800E+02 .2200E+02
.1700E+02 .2000E+02
.1600E+02 .1700E-1-02
.1700E+02 .1300E+02
.2000E+02 .1000E+02
.2300E+02 .7000E+01
.2200E-t-02 .7000E-1-01
.1900E+02 .8000E+01
.1600E+02 .9000E+01
.1300E+02 .9000E+01
.7000E+01 .7000E+01
.3000E+01 .7000E+01
.2000E+01 .9000E+01
4 .4000E+02 18
.3000E+01 .1000E+02
.3000E+01 .2100E+02
.4000E+01 .2400E+02
.5000E+01 .2500E+02
.6000E+01 .2400E+02
.8000E+01 .2000E+02
.1000E+02 .1800E+02
.1300E+02 .1700E+02
.14 OOE+02 .1600E+02
.1500E+02 .1300E+02
.1700E+02 .1100E+02
.1800E+02 .1000E+02
.1700E+02 .1000E+02
.1400E+02 .1100E+02
.13 OOE+02 .1100E+02
.8000E+01 .1000E+02
.5000E+01 .9000E+01
.3000E-1-01 .1000E+02
5 .5000E+02 11
.40002+01 .1200E+02
.4000E+01 .2100E+02
72
-------
.5000E+01
.6000E+01
. 8000E+01
.1000E+02
. 1100E+02
.1000E+02
.8000E+01
.5000E+01
.4000E+01
6
.5000E+01
.5000E+01
.6000E+01
.6000E+01
.8000E+01
.9000E+01
.9000E+01
.8000E+01
.6000E+01
.5000E+01
7
.6000E+01
. 6000E+01
.8000E+01
.8000E+01
.7000E+01
.6000E+01
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1400E+02
.1200E+02
.1100E+02
.1100E+02
.1200E+02
.6000E+02
. .1300E+02
.2100E+02
.2100E+02
.1900E+02
.1600E+02
.1500E+02
.1300E+02
.1200E+02
.1200E+02
.1300E+02
.7000E+02
.1400E+02
.1700E+02
.1500E+02
.1400E+02
.1300E+02
.1400E+02
10
73
-------
FITCOH EXECUTION WITH INTERACTIVE INPUT
74
-------
r 1 I UUIN
ENTER CONTOUR MASTER FILE NAME -> DATCON
ENTER DIAGNOSTIC OUTPUT FILE NAME -> OUTPUT
ENTER HILL ID NUMBER(l-99) -> 2
ENTER HILL NAMEC1-15CHAR.) -> BULLIED HILL
INPUT HILL TOP ELEVATION -> 80
INPUT HILL CENTER X-COORDINATE -> 7
INPUT HILL CENTER Y-COORDINATE -> l N
SPECIFY CONTOUR SELECTION MODE
1.) ALL CONTOURS SELECTED
2.) SELECT RANGE OF CONTOUR IDs
3.) INPUT FILE WITH CONTOUR IDs
CHOICE?(1,2,OR 3) -> 1
PLOT REQUESTED?(Y/N) -> Y
ENTER PLOT FILE NAME -> PLOT1
Please wait...Contour data being processed
Contour ID 1 has been accepted
Contour ID S has been accepted
Contour ID 3 has been accepted
Contour ID 4 has been accepted
Contour ID 5 has been accepted
Contour ID 6 has been accepted
Contour ID 7 has been accepted
ENTER FILE NAME FOR FITTED CONTOUR OUTPUT -> CONOUT
Stop - Program terminated.
DJ\>
75
-------
FITCOH DIAGNOSTIC OUTPUT FILE (OUTPUT)
76
-------
HILL NUMBER 2 IS GULLIED HILL
HILL TOP ELEVATION- .8000E+02
HILL CENTER X-COORDINATE- .7000E+01
HILL CENTER Y-COORDINATE- .1450E+02
ALL CONTOURS IN FILE DATCON SELECTED FOR INPUT
CONTOUR ELEVATION FOR CONTOUR ID 1 - .1000E+02
X-Y COORDINATES INPUT FOR CONTOUR ID 1
.OOOOE+00
.OOOOE+00
.1000E+01
.2000E+01
.4000E+01
.5000E+01
.8000E+01
.1000E+02
.1300E+02
.1500E+02
.1700E+02
.1900E+02
.2200E+02
.2400E+02
.2600E+02
.2800E+02
.3000E+02
.2900E+02
.2600E+02
.2100E+02
.2000E+02
.2000E+02
.2100E+02
.2200E+02
.2600E+02
.2800E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.1700E+02
.1100E+02
.5000E+01
.2000E+01
.OOOOE+00
.4000E+01
.2400E+02
.2800E+02
.2900E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.2300E+02
.2400E+02
.2700E+02
.2900E+02
.3000E+02
.3000E+02
.2900E+02
.2700E+02
.2300E+02
.2200E+02
.2100E+02
.1900E+02
.1700E+02
.1600E+02
.1400E+02
.1300E+02
.1000E-t-02
.8000E+01
.4000E+01
.3000E+01
.2000E+01
.2000E+01
.3000E+01
.2000E+01
.OOOOE+00
.1000E+01
.4000E+01
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
X-Y COORDINATES(EDITED) FOR CONTOUR ID
1 - 35
.OOOOE+00
.OOOOE+00
.1000E+01
.2000E+01
.4000E+01
.2400E+02
•2800E+02
.2900E+02
77
-------
.4000E+01
.5000E+01
.8000E+01
.1000E+02
.1300E+02
.1500E+02
.1700E+02
.1900E+02
.2200E+02
.2400E+02
.2600E+02
.2800E+02
.3000E+02
.2900E+02
.2600E+02
.2100E+02
.2000E+02
.2000E+02
.2100E+02
.2200E+02
.2600E+02
.2800E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.1700E+02
.1100E+02
.5000E+01
.2000E+01
.OOOOE+00
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.2300E+02
.2400E+02
.2700E+02
.2900E+02
.3000E+02
.3000E+02
. 2900E+02
.2700E+02
.2300E+02
.2200E+02
.2100E+02
.1900E+02
.1700E+02
.1600E+02
.1400E+02
.1300E+02
.1000E+02
.8000E+01
.4000E+01
.3000E+01
.2000E+01
.2000E+01
.3000E+01
.2000E+01
.OOOOE+00
.1000E+01
,4000E+01
CONTOUR AREA- .6620E+03
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLZPSE PARAMETERS FOR CONTOUR ZO
.1315E+02
.1468E+02
SEMI-MAJOR AXIS LENGTH- .1632E+02
SEMI-MINOR AXIS LENGTH- .1291E+02
ELLIPSE ECCENTRICITY- .6113E+00
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS-120.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID 2 - .2000E+02
X-Y COORDINATES INPUT FOR CONTOUR ID 2
.1000E+01
.1000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01
.9000E+01
.1200E+02
.1300E+02
.1600E+02
.2100E+02
.2300E+02
.70002+01
.2100E+02
.2500E+02
.2800E+02
.2900E+02
.2700E+02
.2300E+02
.2100E+02
.2100E+02
.2300E+02
.2800E+02
.2900E+02
/8
-------
.2700E+02
.2800E+02
.2800E+02
.2300E+02
.2000E+02
.1900E+02
.1800E+02
.1800E+02
.1900E+02
.2300E+02
.2600E+02
.2600E+02
.2400E+02
.2100E+02
.1700E+02
.1300E+02
.1000E+02
.5000E+01
.3000E+01
.1000E+01
.2600E+02
.2500E+02
.2400E+02
.2400E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1300E+02
.9000E+01
.6000E+01
.5000E+01
.5000E+01
.6000E+01
.7000E+01
.6000E+01
.5000E+01
.4000E+01
.4000E+01
.7000E+01
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
X-Y COORDINATES(EDITED) FOR CONTOUR ID
2 -
32
.1000E+01
.1000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01.
.9000E+01
.1200E+02
.1300E+02
.1600E+02
.2100E+02
.2300E+02
.2700E+02
.2800E+02
.2800E+02
.2300E+02
.2000E+02
.1900E+02
.1800E+02
.1800E+02
.1900E+02
.2300E-K02
.2600E+02
.2600E+02
.2400E+02
.2100E+02
.1700E+02
.1300E+02
.1000E-I-02
.5000E+01
.3000E+01
.1000E+01
.7000E+01
.2100E+02
.2500E+02
.2800E+02
.2900E+02
.2700E+02
.2300E+02
.2100E+02
.2100E+02
.2300E+02
.2800E+02
.2900E+02
.2600E+02
.2500E+02
.2400E+02
.2400E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1300E+02
.9000E-I-01
.6000E+01
.5000E+01
.5000E+01
.6000E+01
.7000E+01
.6000E+01
.5000E+01
.4000E+01
.4000E+01
.7000E+01
79
-------
CONTOUR AREA- .3970E+03
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLIPSE PARAMETERS FOR CONTOUR ID
1128E+02
, 1543E+02
SEMI-MAJOR AXIS LENGTH- .1377E+02
SEMI-MINOR AXIS LENGTH- .9175E+01
ELLIPSE ECCENTRICITY- .7458E+00
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS-120.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID 3
X-Y COORDINATES INPUT FOR CONTOUR ID
.3000E+02
•2000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01
.9000E+01
.1100E+02
.1300E+02
.1800E+02
.2200E+02
.2500E+02
.2600E+02
.2300E+02
.2000E+02
.1800E+02
.1700E+02
.1600E+02
.1700E+02
.2000E+02
.2300E+02
•2200E+02
,1900E+02
.1600E+02
.1300E+02
.7000E+01
•3000E+01
.2000E+01
.9000E+01
.2100E+02
.2500E+02
.2700E+02
.2500E+02
.2100E+02
.1900E+02
.1900E+02
.2400E+02
.2700E+02
.2700E+02
.2500E+02
.2500E+02
.2400E+02
.2200E+02
.2000E+02
.1700E+02
.1300E+02
.1000E+02
.7000E+01
.7000E+01
.8000E+01
.9000E+01
.9000E+01
.7000E+01
.7000E+01
.9000E+01
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
X-Y COORDINATES(EDITED) FOR CONTOUR ID
3 - 27
.2000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01
.9000E+01
.1100E+02
.1300E+02
.9000E+01
.2100E+02
.2500E+02
.2700E+02
.2500E+02
.2100E+02
.1900E+02
.1900E+02
30
-------
. 1800E+02
.2200E+02
.2500E+02
.2600E+02
.2300E+02
.2000E+02
. 1800E+02
.1700E+02
. 1600E+02
.1700E+02
.2000E+02
.2300E+02
.2200E+02
.1900E+02
.1600E+02
.1300E+02
.7000E+01
.3000E+01
.2000E+01
. 2400E+02
. 2700E+02
. 2700E+02
. 2500E+02
. 2500E+02
. 2400E+02
. 2200E+02
. 2000E+02
. 1700E+02
. 1300E+02
. 1000E+02
.7000E+01
.7000E+01
.8000E+01
.9000E+01
.9000E+01
.7000E+01
.7000E+01
. 9000E+01
CONTOUR AREA- .2425E+03
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLIPSE PARAMETERS FOR CONTOUR ID
.1010E+02
.1569E+02
SEMI-MAJOR AXIS LENGTH- .1138E+02
SEMI-MINOR AXIS LENGTH- .6782E+01
ELLIPSE ECCENTRICITY- .8030E+00
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS-110.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID
4 - .4000E+02
X-Y COORDINATES INPUT FOR CONTOUR ID
.3000E+01
.3000E+01
.4000E+01
,5000£+01
.6000E+01
.8000E+01
,1000E+02
.1300E+02
.1400E+02
.1500E+02
.1700E+02
,1800E+02
,1700E+02
,1400E+02
.1300E+02
.8000E+01
.5000E+01
.3000E+01
.1000E+02
.2100E+02
.2400E+02
.2500E+02
.2400E+02
.2000E+02
. 1800E-I-02
.1700E+02
.1600E+02
.1300E+02
.1100E+02
.1000E+02
.1000E+02
.1100E+02
.1100E+02
.1000E+02
.9000E+01
.1000E+02
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(.I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
18
81
-------
X-Y COORDINATES(EDITED) FOR CONTOUR ID
.3000E+01
.3000E+01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1300E+02
.1400E+02
.1500E+02
.1700E+02
.1800E+02
.1700E+02
.1400E+02
.1300E+02
.8000E+01
.5000E+01
.3000E+01
.1000E+02
.2100E+02
.2400E+02
.2500E+02
.2400E+02
.2000E+02
.1800E+02
.1700E+02
.1600E+02
.1300E+02
.1100E+02
.1000E+02
.1000E+02
.1100E+02
.1100E+02
.1000E+02
.9000E+01
. 1000E-I-02
CONTOUR AREA- .1190E+03
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLIPSE PARAMETERS FOR CONTOUR ID
.7972E+01
.1532E+02
SEMI-MAJOR AXIS LENGTH- .8283E+01
SEMI-MINOR AXIS LENGTH- .4573E+01
ELLIPSE ECCENTRICITY- .8337E+00
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS- 40.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID 5
X-Y COORDINATES INPUT FOR CONTOUR ID
.5000E+02
.4000E+01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1100E+02
.1000E+02
.8000E+01
.5000E+01
.4000E+01
.1200E+02
.2100E+02
.2300E+02
.2200E+02
. 1800E-I-02
.1600E+02
.1400E+02
.1200E+02
.1100E+02
.1100E+02
.1200E+02
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
X-Y COORDINATES(EDITED) FOR CONTOUR ID
5 -
11
.4000E+01
.4000E+01
.5000E+01
.6000E+01
.1200E+02
.2100E+02
.2300E+02
.2200E+02
-------
8000E+01
•1000E+02
•1100E+02
.1000E+02
.8000E+01
.5000E+01
•4000E+01
.1800E+02
.1600E+02
.1400E+02
.1200E+02
.1100E+02
.1100E+02
.1200E+02
CONTOUR AREA- .5300E+02
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLIPSE PARAMETERS FOR CONTOUR ID
.6679E+01
.1574E+02
SEMI-MAJOR AXIS LENGTH- .5957E+01
SEMI-MINOR AXIS LENGTH- .2832E+01
ELLIPSE ECCENTRICITY- .8798E+00
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS- 20.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID
6 -
.6000E+02
X-Y COORDINATES INPUT FOR CONTOUR ID
.5000E+01
.5000E+01
.6000E+01
.6000E+01
.8000E+01
.9000E+01
.9000E+01
.8000E+01
.6000E+01
.5000E+01
.1300E+02
.2100E+02
.2100E+02
,1900E+02
.1600E+02
, 1500E-I-02
,1300E+02
.1200E+02
•1200E+02
.1300E+02
CONTOUR HAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID
X-Y COORDINATES(EDITED) FOR CONTOUR ID
.5000E+01 .1300E+02
.5000E+01 .2100E+02
.6000E+01 .2100E+02
.6000E+01 .1900E+02
.8000E+01 .1600E-I-02
.9000E+01 .1500E+02
.9000E+01 .13 OOE+02
.8000E+01 .1200E+02
.6000E+01 .1200E+02
.5000E+01 .13 OOE+02
CONTOUR AREA- .2250E+02
X-COORDINATE OF CONTOUR CENTROID-
Y-COORDINATE OF CONTOUR CENTROID-
ELLIPSE PARAMETERS FOR CONTOUR ID
SEMI-MAJOR AXIS LENGTH- .4562E+01
6 - 10
.6585E+01
.1544E+02
83
-------
SEMI-MINOR AXIS LENGTH- .1570E+01
ELLIPSE ECCENTRICITY- ' .9389E+oo
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS- 20.00 DEGREES
CONTOUR ELEVATION FOR CONTOUR ID 7 - .7000E+02
X-Y COORDINATES INPUT FOR CONTOUR ID 7
.6000E+01 .1400E+02
.6000E+01 .1700E+02
.8000E+01 .1500E+02
.8000E+01 .1400E+02
.7000E+01 .1300E+02
.6000E+01 .1400E+02
CONTOUR WAS FOUND TO BE A SINGLE CONTOUR.
(I.E. NO CONTOUR CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT.)
MODIFIED NUMBER OF POINTS FOR CONTOUR ID 7 - 6
X-Y COORDINATES(EDITED) FOR CONTOUR ID 7
.6000E+01 .1400E+02
.6000E+01 .1700E+02
.8000E+01 .1500E+02
.8000E+01 .1400E+02
.7000E+01 .1300E+02
.6000E+01 .1400E+02
CONTOUR AREA- .5000E+01
X-COORDINATE OF CONTOUR CENTROID- .6867E+01
Y-COORDINATE OF CONTOUR CENTROID- .1480E+02
ELLIPSE PARAMETERS FOR CONTOUR ID 7
SEMI-MAJOR AXIS LENGTH- .1764E+01
SEMI-MINOR AXIS LENGTH- .9025E+00
ELLIPSE ECCENTRICITY- .SSSIE+OO
ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE
POSITIVE X-AXIS- 20.00 DEGREES
-------
FITCOH PLOT FILE (PLOT 1)
85
-------
FITCON
2 GULLIED HILL
.7000E+01
7
1
2
3
4
5
6
7
.OOOOE+00
.OOOOE+00
35
.OOOOE+00
.OOOOE+OQ
.1000E+01
.2000E+01
.4000E-I-01
.5000E-I-01
.8000E+01
.1000E+02
.1300E+02
.1500E+02
.1700E+02
.1900E+02
.2200E+02
.2400E+02
.2600E+02
.2800E+02
.3000E-I-02
.2900E+02
.2600E+02
.2100E+02
.2000E+02
.2000E-I-02
.2100E-I-02
.2200E+02
.2600E-I-02
.2800E+02
.3000E-I-02
.3000E+02
.2800E-I-02
.2400E-I-02
. 1700E-I-02
.1100E+02
.5000E+01
.2000E+01
.OOOOE+00
35
.OOOOE+00
.OOOOE+00
.1000E+01
.2000E+01
.4000E+01
.5000E+01
.8000E+01
.1000E+02
.1300E+02
.1450E+02
.3000E+02
.3000E+02
.1000E+02
.4000E+01
.2400E+02
.2800E+02
.2900E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.2300E+02
.2400E+02
.2700E+02
.2900E+02
.3000E+02
.3000E+02
.2900E+02
.2700E+02
<2300E+02
.2200E+02
.2100E+02
.1900E+02
.1700E+02
.1600E+02
.1400E+02
.1300E+02
.1000E+02
.8000E+01
.4000E+01
.3000E+01
.2000E+01
.2000E+01
.3000E+01
.2000E+01
.OOOOE+00
.1000E+01
.4000E+01
.1000E+02
.4000E+01
.2400E+02
.2800E+02
.2900E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.2300E+02
, OOOOE+00
, OOOOE+00
.3000E+02
.3000E+02
86
-------
.1500E+02
.1700E+02
.1900E+02
.2200E+02
.2400E+02
.2600E+02
.2800E+02
.3000E+02
.2900E+02
.2600E+02
.2100E+02
.2000E+02
.2000E+02
.2100E+02
.2200E+02
.2600E+02
.2800E+02
.3000E+02
.3000E+02
.2800E+02
.2400E+02
.1700E+02
.1100E+02
.5000E+01
.2000E+01
.OOOOE+00
32
.1000E+01
.1000E+01
.2000E+01
.3000E+01
.5000E-I-01
.7000E+01
. 9000E-I-01
.1200E+02
.1300E+02
.1600E+02
.2100E+02
.2300E+02
.2700E+02
.2800E+02
.2800E+02
.2300E+02
.2000E+02
.1900E+02
.1800E+02
.1800E+02
.1900E+02
.2300E+02
.2600E+02
.2600E+02
.2400E-I-02
.2100E+02
.1700E+02
.1300E+02
.1000E+02
.5000E+01
.3000E+01
.1000E+01
32
.2400E-I-02
.2700E-I-02
.2900E+02
.3000E+02
.3000E-I-02
.2900E-I-02
.2700E+02
.2300E-I-02
.2200E-I-02
.2100E-I-02
. 1900E-I-02
.1700E+02
.1600E+02
.1400E+02
.1300E-I-02
.1000E+02
.8000E-t-01
.4000E-I-01
.3000E-I-01
.2000E+01
.2000E-I-01
.3000E-I-01
.2000E+01
.OOOOE+OO
.1000E-t-01
.4000E+01
.2000E+02
.7000E+01
.2100E+02
.2500E+02
.2800E+02
.2900E+02
.2700E+02
.2300E-I-02
.2100E+02
.2100E+02
.2300E+02
.2800E+02
.2900E+02
.2600E+02
.2500E+02
.2400E+02
.2400E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1300E+02
.9000E+01
.6000E+01
.5000E+01
.5000E-I-01
. 6000E-I-01
.7000E+01
.6000E+01
.5000E+01
.4000E+01
.4000E+01
.7000E-I-01
.2000E+02
87
-------
.1000E+01
.1000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01
.9000E+01
.1200E-I-02
.1300E+02
.1600E+02
.2100E+02
.2300E+02
.2700E+02
.2800E+02
.2800E+02
.2300E+02
.2000E+02
.1900E+02
.1800E-I-02
.1800E+02
.1900E+02
.2300E+02
.2600E+02
.2600E+02
.2400E+02
.2100E+02
.1700E+02
.1300E+02
.1000E+02
.5000E+01
.3000E+01
.1000E+01
27
.2000E+01
.2000E+01
.3000E+01
.5000E-I-01
.7000E+01
.9000E+01
.1100E+02
.1300E+02
.1800E+02
.2200E+02
.2500E+02
.2600E+02
.2300E-I-02
.2000E+02
. 1800E-I-02
.1700E+02
. 1600E+02
.1700E+02
.2000E-I-02
.2300E+02
.2200E+02
. 1900E+02
. 1600E+02
. 1300E+02
.7000E+01
.3000E+01
.2000E-H01
.7000E+01
.2100E+02
.2500E+02
.2800E+02
.2900E+02
.2700E+02
.2300E+02
.2100E+02
.2100E+02
.2300E+02
.2800E+02
.2900E+02
.2600E+02
.2500E+02
.2400E+02
.2400E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1300E+02
.9000E+01
.6000E-I-01
.5000E+01
.5000E+01
.6000E+01
.7000E+01
.6000E+01
.5000E-I-01
,4000E-I-01
.4000E+01
.7000E+01
.3000E-I-02
.9000E+01
.2100E+02
.2500E+02
.2700E+02
.2500E+02
.2100E+02
.1900E+02
.1900E+02
.2400E+02
.2700E+02
.2700E+02
.2500E+02
.2500E+02
.2400E+02
.2200E+02
.2000E+02
.1700E+02
.1300E+02
.1000E-I-02
.7000E+01
.7000E+01
. 8'OOOE+01
.9000E+01
.9000E+01
.7000E-I-01
.7000E+01
.9000E+01
38
-------
27
.2000E+01
.2000E+01
.3000E+01
.5000E+01
.7000E+01
.9000E+01
.1100E+02
.1300E+02
.1800E+02
.2200E+02
.2500E+02
.2600E+02
.2300E+02
.2000E+02
.1800E+02
.1700E+02
.1600E+02
.1700E+02
.2000E+02
.2300E+02
.2200E+02
.1900E+02
.1600E+02
.1300E+02
.7000E+01
.3000E+01
.2000E+01
19
.3000E+01
.3000E-I-01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1300E+02
.1400E+02
.1500E+02
.1700E+02
.1800E+02
.1700E+02
.1400E+02
.1300E+02
.8000E+01
.5000E+01
.3000E+01
18
.3000E+01
.3000E+01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1300E+02
.1400E+02
.1500E+02
.1700E+02
.1800E+02
.3000E+02
.9000E+01
.2100E+02
.2500E+02
.2700E+02
.2500E+02
.2100E+02
.1900E+02
.1900E+02
.2400E+02
.2700E+02
.2700E+02
.2500E+02
.2500E+02
.2400E+02
.2200E+02
.2000E+02
.1700E+02
.1300E+02
.1000E+02
.7000E+01
.7000E+01
.8000E-I-01
.9000E+01
.9000E+01
.7000E+01
.7000E+01
.9000E+01
.4000E+02
.1000E+02
.2100E+02
.2400E+02
.2500E+02
.2400E+02
.2000E+02
.1800E+02
.1700E+02
.1600E+02
.1300E+02
.1100E+02
.1000E+02
.1000E+02
.1100E-I-02
.1100E+02
.1000E+02
.9000E+01
.1000E+02
.4000E+02
.1000E+02
.2100E+02
.2400E+02
.2500E+02
.2400E+02
.2000E+02
.1800E+02
.1700E+02
.1600E+02
.1300E+02
.1100E+02
.1000E+02
89
-------
.1700E+02
.1400E+02
.1300E+02
.8000E+01
.5000E+01
.3000E+01
11
.4000E+01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1100E+02
.1000E+02
.8000E+01
.5000E+01
.4000E+01
11
.4000E+01
.4000E+01
.5000E+01
.6000E+01
.8000E+01
.1000E+02
.1100E+02
.1000E+02
.8000E+01
.5000E+01
.4000E+01
10
.5000E+01
.5000E+01
.6000E+01
.6000E+01
.8000E+01
.9000E+01
.9000E+01
.8000E+01
.6000E+01
.5000E+01
10
.5000E+01
.5000E+01
.6000E+01
.6000E+01
.8000E+01
.90002+01
.9000E+01
.8000E+01
.6000E+01
.5000E+01
6
.6000E+01
.6000E+01
.8000E+01
.8000E+01
.7000E-I-01
.6000E+01
6
.1000E+02
.1100E+02
.1100E+02
.1000E+02
.9000E+01
.1000E+02
.5000E+02
.1200E+02
.2100E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1400E+02
.1200E+02
.1100E+02
.1100E+02
.1200E+02
.5000E+02
.1200E+02
.2100E+02
.2300E+02
.2200E+02
.1800E+02
.1600E+02
.1400E+02
.1200E+02
.1100E+02
.1100E+02
.1200E+02
.6000E+02
.1300E+02
.2100E+02
.2100E+02
.1900E+02
.1600E+02
.1500E+02
.1300E+02
.1200E+02
.1200E+02
.1300E+02
.6000E+02
.1300E+02
.2100E+02
.2100E+02
.1900E+02
.1600E+02
.1500E+02
.1300E+02
.1200E+02
.1200E+02
.1300E+02
.7000E+02
.1400E+02
.1700E+02
.1500E+02
.1400E+02
. 1300E-I-02
.1400E+02
.7000E+02
90
-------
6000E+01
6000E+01
8000E+01
8000E+01
7000E+01
6000E+01
1315E+02
1128E+02
1010E+02
7972E+01
6679E+01
6585E+01
6867E+01
.1400E+02
.1700E+02
. 1500E+02
.1400E+02
.1300E+02
.1400E+02
. 1468E+02
. 1543E+02
.1569E+02
.1532E+02
. 157.4E+02
. 1544E+02
. 1480E+02
. 1632E+02
.1377E+02
. 1138E+02
.8283E+01
.5957E-H01
.4562E+01
.1764E+01
.1291E+02
.9175E-H01
. 6782E+01
.4573E+01
.2832E+01
.1570E+01
.9025E+00
.1200E+03
.1200E+03
.1100E+03
.4000E+02
.2000E+02
.2000E+02
.20003+02
-------
FITCOH OUTPUT FILE FOR HCRIT INPUT (CONOUT)
92
-------
B GULLIED HILL
.BOOOE+O8
7
1
2
3
5' . •
6
7
. 1000E+08 . 1319E+08 . 146BE+08 . 1638E+OS .1891E+Oe .6113E+00 . 120OE+03
.EOOOE+Oe .118BE+O8 . 13A3E+08 .1377E+O8 .917SE+O1 .7
-------
HCRIT EXECUTION WITH INTERACTIVE INPUT
-------
A:\T£RRAIN>HCRIT
ENTER INPUT PILE NAME(FROM FITCON) -> CONOUT
ENTER OUTPUT FILE NAME(FOR CTDM) ->TERRAIN
PLOT REQUESTED?(Y/N) -> Y
ENTER PLOT FILE NAME -> PLOT2
SPECIFY CRITICAL HEIGHT SELECTION MODE
1.) AT ALL CONTOUR ELEVATIONS EXCEPT UPPERMOST
2.) EVENLY SPACED BETWEEN A USER SUPPLIED ELEVATION
AND THE UPPERMOST CONTOUR ELEVATION
CHOICE?(1 OR 2) -> 1
A:\TERRAIN>
95
-------
HCRIT PLOT FILE (Plot 2)
96
-------
NCR IT
a GULLIED HILL
7
I
e
3
*
3
&
7
.BOOOE+OB
. IOOOE+OB
.BOOOE+OB
.3OOOE+OB
.4OOOE+OB
.3OOOE+OB
.&OOOE+OB
.7000E+OB
.
. 1OOOE+O8 .BB<»7E+OI . 13<»OE+OB .4B6BE+O8 . 1791E+O1 .I467E+O1 .AB93E+O1 .3333E+O1
.BOOOE+Oe .7641E+O1 . 13<»OE+OB .37<»3C+OS . I799E+O1 . 13<*OE+OI .93B3E+OI .S62^E+01
.3000E+03 .7086E+01 . 1332E+OB .B<»7IE+OB . 1B33E>O1 . 1676E>O1 .<»<»6<»E+O1 .BO7OE+O1
.OB .BOOOE+OB .1BO3E+O1 .19BIE+OI .3633E+O1 . 1589E+O1
.3000E+OB .67B6E+01 . 131BE>OB .BQOOE^OB .!<»99E-fOl .BSO'.E+Ol .BB37E+O1 . 119OE+O1
.6000E+OB .6B67E+01 . 1<»BOE+OB .BOOOE+OB .BOOOE+OI .BOOOE+O1 .1764E+01 .9OB5E+OO
-------
HCRIT OUTPUT FILE FOR INPUT TO CTDH (TERRAIN)
98
-------
2 6
10.000
20.000
30.000
40.000
50.000
60.000
10.000
20.000
30.000
40.000
50.000
.1315E+02
. 1128E-I-02
.1010E+02
.7972E+01
.6679E+01
.6585E+01
.8247E+01
.7641E+01
.7026E+01
.6710E+01
.6726E+01
. 8000E-I-02
.1468E+02
.1543E+02
.1569E+02
.1532E+02
.1574E+02
.1544E+02
.1540E+02
.1540E+02
.1532E+02
. 1533E+02
. 1512E+02
GULLIED HILL
60.
60.
70.
140.
160.
160.
131.
142.
155.
160.
160.
000
000
000
000
000
000
319
548
290
000
000
16.
13.
11.
s:
5.
4.
1.
1.
1.
1.
1.
320
770
380
283
957
562
791
799
853
805
459
12.
9.
6.
4.
2.
1.
1.
1.
1.
1.
2.
910
175
782
573
832
570
467
540
676
921
504
6
5
4
3
2
.295
.383
.464
.633
.837
3.233
2.624
2.070
1.589
1.190
60.000 .6867E+01 .1480E+02 160.000 2.000 2.000 1.764* .902
-------
PLOTCON (HPLTCON) EXECUTION WITH INTERACTIVE INPUT
100
-------
D:\>PLOTCON"
INPUT NAME OF PLOTFILE FROM PROGRAM FITCON—>? PLOT1
SELECT TYPE OF DISPLAY
1.) Low resolution with color
2.) High resolution black and whits
Choic«?(l or 2)—>? 2
SELECT THE CONTOUR TYPE FOR DISPLAY
1.) Un«dit»d Contours
2.) Edited Contours
Choice?U or 2)—>? 2
DISPLAY FITTED CUTOFF HILL CONTOURS?(Y/N)->? Y
INPUT NAME OF PLOTFILE FROM PROGRAM HCRIT? PLOT2
101
-------
PLOT OF ACTUAL INPUT CONTOURS
102
-------
103
-------
PLOT OF ACTUAL CONTOURS AND FITTED ELLIPSES
104
-------
105
-------
PLOT OF ACTUAL CONTOURS AND INVERSE POLYNOMIAL
CONTOURS AT ACTUAL CONTOUR ELEVATIONS
ABOVE THE CRITICAL CUTOFF ELEVATION
106
-------
o
a
CB
u
107
-------
0
CM
o
a
a>
i-i
u
a
u
u
108
-------
«0
I
*•*
cu
to
o
109
-------
o
(9
-------
o
in
o
cu
B
U
u
111
-------
o
vO
O
a
4)
!•*
cu
«5
U
(4
u
112
-------
APPENDIX D
PROGRAM LISTINGS
113
-------
FITCOH MAIN PROGRAM AND SUBROUTINES
114
-------
PROGRAM FITCON
C***PROGRAM TO FIT DIGITIZED CONTOURS TO ELLIPTICAL SHAPES. PROGRAM
C***C£NERATES A FILE OF ELLIPTICAL CONTOUR PARAMETERS TO BE USED BY
C***PROGRAM HCRIT TO PERFORM THE CRITICAL HEIGHT ANALYSIS FOR THE
C***HILL IN QUESTION. A PLOT FILE IS ALSO GENERATED FOR SUBSEQUENT
C***DISPLAY OF DIGITIZED AND FITTED CONTOURS.
C***
C GLOSSARY OF TERMS
C A(J)-CALCULATED SEMI-MAJOR AXIS LENGTH(USER COORDINATES) FOR THE
C ELLIPTICAL REPRESENTATION OF CONTOUR J
C AFIL-ANGULAR FILTER SIZE(1 TO 22.5 DEGREES) INPUT BY THE USER
C FOR THE CONTOUR COMPLETION ANALYSIS. MODIFIED AFTER INPUT
C SO THAT IT DIVIDES EVENLY INTO 360 DEGREES.
C ANGLE(M)-(M-1)*10.0 WHERE M-1,18
C ANS-CHARACTER*! VARIABLE HOLDING THE ANSWER TO A YES(Y) OR NO(N)
C QUESTION
C AR-VALUE OF THE AREA RETURNED BY A CALL TO SUBROUTINE ARCM.
C AR WILL BE POSITIVE IF THE CONTOUR POINTS ARE GIVEN IN A '
C CLOCKWISE FASHION AND NEGATIVE IF THE CONTOUR POINTS ARE GIVEN
C IN A COUNTER-CLOCKWISE FASHION
C .ARCH-SUBROUTINE TO CALCULATE THE CONTOUR AREA AND CENTROID
C COORDINATES
C AREA-AREA OF A GIVEN CONTOUR(-ABS(AR))
C B(J)-CALCULATED SEMI-MINOR AXIS LENGTH(USER COORDINATES) FOR THE
C ELLIPTICAL REPRESENTATION OF CONTOUR J
C CTLAG-CONTOUR CLOSURE INDICATOR
C -0(CONTOUR OPEN)
C -1(CONTOUR CLOSED)
C CH(M)-COS(PI*(M-l)*10.0/180.) WHERE M-1,18
C CONIN-UNIT NUMBER FOR FILE CONTAINING CONTOUR IDs FOR THE HILL
C IN QUESTION
C CONFILE-CHARACTER*15 VARIABLE GIVING THE NAME OF THE FILE
C CONTAINING CONTOUR IDs FOR THE HILL IN QUESTION
C CONCOMP-SUBROUTINE WHICH ADDS POINTS TO COMPLETE A CONTOUR
C COOT-UNIT NUMBER FOR OUTPUT FILE COUTFILE WHICH WILL BE INPUT TO
C THE CRITICAL HEIGHT ANALYSIS PROGRAM
C COUTFIL£-CHARACTER*15 VARIABLE GIVING THE NAME OF THE OUT&JT FILE
C CONTAINING THE FITTED HILL PARAMETERS WHICH WILL BE INPUT
C TO THE CRITICAL HEIGHT ANALYSIS PROGRAM(HCRIT)
C DFTOL-DISTANCE FROM FIRST TO LAST CONTOUR POINT(USER COORDINATES)
C DOUT-UNIT NUMBER FOR FILE CONTAINING DIAGNOSTIC OUTPUT
C DOUTFILE-CHARACTER*15 VARIABLE GIVING THE NAME OF THE FILE
C CONTAINING DIAGNOSTIC OUTPUT FOR THE HILL IN QUESTION
C ECC(J)-ECCENTRICITY OF THE ELLIPSE REPRESENTING CONTOUR J
C -SQRT(A(J)**2-B(J)**2)/A(J)
C HCON(J)-ELEVATION OF HILL CONTOUR J(USER COORDINATES)
C HCONT-VALUE OF HCON(J) FOR A PARTICULAR CONTOUR J
C HNAME-CHARACTER*15 VARIABLE GIVING THE RILL NAME
C HTOP-HILL TOP ELEVATION(USER COORDINATES)
C ICL-SMALLEST ID(1-9999) NUMBER FOR THE CONTOUR GROUP(INPUT ONLY FOR
C ICMODE-2)
C ICMODE-CONTOUR INPUT MODE FOR THE HILL IN QUESTION
C -1(ALL CONTOURS IN THE MASTER FILE SELECTED FOR INPUT)
C -2(CONTOUR ID RANGE SPECIFIED FOR INPUT)
C -3(INPUT FILE WITH CONTOUR IDs SPECIFIED)
C ICU-LARGEST ID NUMBER(1-9999) FOR THE CONTOUR GROUP(INPUT ONLY FOR
C ICMODE-2)
FIT00010
FIT00620
FIT00030
FIT00040
FITOOpSO
FIT00660
FIT90070
FITOOdSO
FITOOU90
FIT00100
FITOOllO
FITOO:i20
FTTOO:L30
FITOO:L40
FIT00150
FITOO!L60
FIT00170
FITOOPSO
FIT00190
FIT00200
FIT00210
FIT00220
FTT00230
FIT00240
FIT00250
Firoo^eo
FITOOJ70
FITOO280
FIT00290
FITOO)oo
FITOO?10
FTT00320
FTT00330
FTT00340
FITOO!} 50
FIT00360
FIT00370
FITOO JJ 80
FIT00390
FTT00400
FIT00410
FIT00420
FITOO430
FIT00440
FIT00450
FIT00460
FIT00470
FIT00480
FIT00490
FITOO
FITOO
FITOO
FITOO
FITOO
FITOO
00
10
20
30
40
50
FIT00560
FITOO?70
FITOOPSO
FIT00590
FIT00600
115
-------
C IDC(J)-ID NUMBER FOR CONTOUR J WHICH HAS BEEN SELECTED FROM THE FIT00610
C CONTOUR MASTER FILE FIT00620
C IDCPK(I)-ID NUMBER FOR THE Ith CONTOUR SPECIFIED IN FILE CONFILE FIT00630
C IDHILL-HILL ID NUMBER(l-999) SPECIFIED BY THE USER FIT00640
C IN-UNIT NUMBER FOR CONTOUR MASTER FILE FIT00650
C ISMFLG-COHPLETION CODE RETURNED BY SUBROUTINE SMOMNT FIT00660
C -0(RADIUS OF GYRATION WAS CALCULATED) FIT00670
C -1(RADIUS OF GYRATION COULD NOT BE CALCULATED) FIT00680
C J-CURRENT NUMBER OF CONTOURS INPUT FROM THE MASTER FILE FOR THE FIT00690
C HILL IN QUESTION(AFTER QUALIFICATION AND EDITING) FTT00700
C LTPR-WORKING ARRAY USED BY SUBROUTINE ISORT FIT00710
C MASTER-CHARACTER*15 VARIABLE GIVING THE NAME OF THE MASTER FILE FIT00720
C CONTAINING THE CONTOUR ELEVATIONS AND POINT COORDINATES FIT00730
t MCFLAG-MULTIPLE CONTOUR SUBROUTINE COMPLETION CODE RETURNED FROM FIT00740
C SUBROUTINE MULTC FIT00750
C -0(MAXIMUM NUMBER OF POINTS EXCEEDED IN THE CONTOUR POINT FIT00760
C REASSIGNMENT PROCESS—CONTOUR REJECTED) FIT00770
C . -1(THE LAST IN A SERIES OF MULTIPLE CONTOURS WAS FOUND NOT FIT00780
C TO BE CLOSED—CONTOUR REJECTED) FIT00790
t . -2(CONTOUR WAS FOUND TO BE A SINGLE CONTOUR(I.E. NO CONTOUR FIT00800
C CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT)) FIT00810
C -3(POINT REASSIGNMENT FOR THE MULTIPLE CONTOUR WAS FIT00820
C SUCCESSFULLY COMPLETED) FIT00830
C HOTOTAL NUMBER OF CONTOURS SELECTED FROM THE MASTER FILE FOR THE FIT00840
C HILL IN QUESTION FTT00850
t HCID-NUMBER OF REQUESTED CONTOUR IDs CONTAINED IN CONFILE FIT00860
C HCMAX-MAXIMUM NUMBER OF CONTOURS ALLOWED FIT00870
C MCT2-2*NC FIT00880
C MFIL-INT(3«0./AFIL) FIT00890
C MPC-NUMBER OF POINTS ON A CONTOUR . FIT00900
C VPCMAX-MAXIMUM NUMBER OF POINTS PER CONTOUR ALLOWED FIT00910
C FPCSV-NUMBER OF POINTS ON A CONTOUR PRIOR TO CONTOUR COMPLETION FIT00920
6 HSLOPE-HUMBER OF LINES USED IN THE DETERMINATION OF THE LINE, FIT00930
C PASSING THROUGH THE CONTOUR CENTROID, WHICH GIVES THE FIT00940
C MAXIMUM RADIUS OF GYRATION FOR THE DIGITIZED CONTOUR FIT00950
t OREN(J)-ANGLE CORRESPONDING TO THE ORIENTATION OF THE SEMI- FTT00960
C MINOR AXIS OF CONTOUR J. THE POSSIBLE ORIENTATIONS REPRESENTFIT00970
t THE FOLLOWING DIRECTIONS WITH RESPECT TO THE POSITIVE FIT00980
C' X-AXXS:0,10,20,30,40,30,60,70,80f90,100,110,120,130,140,150,FIT00990
C 160,AND 170 DEGREES FIT01000
t ORENT-CONTOUR MINOR AXIS ORIENTATION CORRESPONDING TO THE MAXIMUM FIT01010
C RADIUS OF GYRATION RETURNED BY SUBROUTINE SMOMNT. ORENT IS FIT01020
C SIMPLY A TEMPORARY HOLDING VARIABLE FOR OREM(J) FIT01030
C PI-3.14159265 FIT01040
C PFILE-CHARACPER*15 VARIABLE GIVING THE NAME OF THE PLOT FILE FIT01050
C PFLAG-PLOT GZNERATON INDICATOR FIT01060
C -0(NO PLOT GENERATED) • FIT01070
-1(PLOT GENERATED) FIT01080
RAD-RADIUS OF THE EQUIVALENT CIRCULAR CONTOUR(USER COORDINATES) FIT01090
RG-MAXIMUM RADIUS OF GYRATION CONSIDERING THE 18 ORIENTATIONS OF FIT01100
AXES PASSING THROUGH THE CONTOUR CENTROID IN THE PLANE OF FIT01110
THE CONTOUR(USER COORDINATES) * FIT01120
RGRAT-THE RATIO OF THE DIFFERENCE BETWEEN THE MAXIMUM AND MINIMUM FIT01130
RADII OF GYRATION(CONSIDERING THE 18 ORIENTATIONS OF AXES FIT01140
PASSING THROUGH THE CONTOUR CENTROID) TO THE MAXIMUM RADIUS FIT01150
C OF GYRATION. USED TO DETERMINE WHETHER AN INPUT CONTOUR SHOULDFIT01160
C BE REPRESENTED BY A CIRCLE FIT01170
C SKIPCN-SUBROUTINE TO SKIP OVER CONTOUR POINTS FOR CONTOURS WHICH AREFIT01180
C NOT PROCESSED . FIT01190
C SMOMNT-SUBROUTINE WHICH CALCULATES THE MAXIMUM RADIUS OF GYRATION FIT01200
116
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
•c
c
c
c
c
c
c
c
c
c
c
c
c
c
c***
c***
FOR AM INPUT CONTOUR BY CONSIDERING THE CALCULATED RADII OF FIT01210
GYRATION FOR 18 LINES OF EQUAL ANGULAR SPACING WHICH PASS FIT01220
THROUGH THE CONTOUR CENTROID IN THE PLANE OF THE CONTOUR FIT01230
SN(M)-SIN(PI*(M-l)*10.0/180.) WHERE M-1,18 FIT01^40
ISORT-SUBROUTINE TO CARRY OUT A SIMPLE PARAMETER SORT FIT01250
UPL-ONIT NUMBER FOR THE FILE PFILE FIT01260
UPSCR-UNIT NUMBER FOR THE SCRATCH FILE "PSCRAT" FITOIJTO
XOCONTOUR CENTROID X-COORDINATE RETURNED BY A CALL TO ARCM. FIT01280
XC IS SIMPLY A TEMPORARY HOLDING VARIABLE FOR XCM(J) FIT01290
XCM(J)-CALCULATED X-COORDINATE(USER COORDINATES) OF THE CENTER OF FIT01300
MASS OF CONTOUR J FIT01310
XCON(K)-X-COORDINATE
-------
C*** FIT01810
!C* "INITIALIZE THE ANGLE ARRAY TO BE USED FOR THE CONTOUR FIT01820
^•••ORIENTATION ANALYSIS. " ' FIT01830
DATA ANGLE/0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,110., FIT01840
&120.,130.,140.,ISO.,160.,170./ FIT01850
C***INITIALIZE SINE AND COSINE ARRAYS TO BE USED FOR THE CONTOUR FIT01860
C***ORIENTATION ANALYSIS. FIT01870
DATA SN/0.,0.1736,0.3420,0.5,0.6428,0.7660,0.8660,0.9397, FIT01880
40.9848,1.0,0.9848,0.9397,0.8660,0.7660,0.6428,0.5,0.3420,0.1736/ FIT01890
DATA CN/1.0,0.9848,0.9397,0.8660,0.7660,0.6428,0.5,0.3420,0.1736, FIT01900
&0.0,-0.1736,-0.3420,-0.5,-0.6428,-0.7660,-0.8660,-0.9397,-0.9848/ FIT01910
C***SPECIFY FILE UNIT NUMBERS. FIT01920
CONIN-14 FIT01930
IN-1S FIT01940
DOUT-1S FIT01950
COUT-17 FIT01960
UPL-18 FIT01970
UPSCR-19 FIT01980
'C***SPECIFY CONSTANTS. FIT01990
PI-3.14159265 FIT02000
NCMAX-200 . FIT02010
NPCMAX-1000 FIT02020
NSLOPE-18 FIT02030
C*** FIT02040
C*** . FIT02050
C INPUT FILE NAMES(MASTER FILE AND DIAGNOSTIC OUTPUT FILE) AND FIT02060
C HILL IDENTIFICATION INFORMATION. FIT02070
C*** FIT02080
C*** FIT02090
;C***INPUT NAMES FOR THE CONTOUR MASTER FILE AND THE DIAGNOSTIC FIT02100
C***OUTPUT FILE. . FIT02110
! 5 WRITE(*,10) • FIT02120
! 10 FORMAT(/,IX,'ENTER CONTOUR MASTER FILE NAME -> «\) FIT02130
READ(*,'(A)•) MASTER ' FIT02140
IF(MASTER.EQ.' •) GO TO 5 . FIT02150
1 15 WRITE(*,20) . FTT02160
20 FORMAT(/,IX,'ENTER DIAGNOSTIC OUTPUT FILE NAME -> *\) FIT02170
READ(*,'(A)') DOUTFILE ' FIT02180
!* IF(DOUTFILE.EQ.• •) GO TO 15 FIT02190
C***OPEN THE CONTOUR MASTER FILE AND THE DIAGNOSTIC OUTPUT FILE. FIT02200
OPEN(IH,FILE-MASTER,STATUS-'OLD') FIT02210
OPEH(DOUT,FILE-DOUTFILE,STATUS-'NEW') FIT02220
;C***INPUT HILL IDENTIFIER NUMBER AND HILL NAME. FIT02230
! 25 WRITE(*,30) FIT02240
30 FORMAT(/,IX,'ENTER HILL ID NUMBER(1-99) -> >\) FIT02250
READ(*,'(BN,I2)',ERR-25) IDHILL FIT02260
! IF(IDHILL.EQ.O) GO TO 25 - FIT02270
! 35 WRITE(*,40) FIT02280
40 FORMAT(/,IX,'ENTER HILL NAME(1-15CHAR.) -> «\) FIT02290
READ(*,'(A)') HNAME FIT02300
IF(HNAME.EQ.' •) GO TO 35 FIT02310
WRITE(DOUT,50) IDHILL,HNAMB • . FIT02320
50 FORHAT(/,1X,'HILL NUMBER',I4,1X,'IS',1X,A15) FIT02330
C*** FIT02340
C*** FIT023SO
C INPUT THE HILL TOP ELEVATION AND THE COORDINATES OF THE FIT02360
C HILL CENTER. FIT02370
C*** • FIT02380
C*** FIT02390
C***INPUT THE HILL TOP ELEVATION. FIT02400
118
-------
60 CONTINUE
WRITE(*,70)
70 FORMAT(/,IX,'INPUT HILL TOP ELEVATION -> '\)
READ(*,'(BN,F10.0)',ERR-60) HTOP
C***INPOT THE HILL CENTER X AND Y COORDINATES.
80 WRITE(*, HO)
110 FORMAT(/,IX,'INPUT HILL CENTER X-COORDINATE -> •
READ(*,'(BN,F10.0)',ERR-80) XHTOP
115 WRITE(*,120)
120 FORMAT(/,IX,'INPUT RILL CENTER Y-COORDINATE -> '
READ(*,'(BN,F10.0)»,ERR-115) YHTOP
C***DETERMINE WHETHER ANGULAR FILTERING IS TO BE USED.
WRITE(*,1201)
1201 FORMAT(/,IX,'ANGULAR FILTERING?(Y/N) -> -\)
READ(*,'(A)') ANS
IF(ANS.EQ.'Y'.OR.ANS.EQ.'y') GO TO 1202
NFIL-0
GO TO 1261
1202 CONTINUE
\)
\)
FIT02410
FIT024*0
FIT024JO
FIT02440
FIT02450
FIT02460
FTT02470
FIT024IO
FIT02490
FIT02500
FIT02510
FIT02520
FIT02530
FIT02540
FIT02S50
FIT02560
FIT02570
FIT02580
FIT02590
C***INPOT ANGULAR FILTER SIZE FOR USE IN THE CONTOUR COMPLETION ANALYSISFIT02600
121 WRITE(*,122) FIT02610
122 FORMAT(/,IX,'INPUT ANGULAR FILTER SIZE FOR CONTOUR COMPLETION(1-22FIT02620
«.5 DEC.) -> «\) FIT02630
READ(*,'(BN,F10.0)',ERR-121) AFIL FIT02640
IF(AFIL.GE.1.0.AND.AFIL.LE.22.5) GO TO 124 FIT02630
WRITE(*,123) FIT02660
123 FORMAT(/,IX,'***ERROR*** SPECIFIED FILTER SIZE OUT OF RANGE—TRY AFIT02670
*GAHt«) FIT02680
GO TO 121 FIT02690
C***WRTTE SPECIFIED ANGULAR FILTER SIZE TO THE DIAGNOSTIC OUTPUT FILE. FTE02700
124 WRITE(DOOT,123) AFIL FIT02710
125 FORMAT(/,IX,'SPECIFIED ANGULAR FILTER SIZE FOR CONTOUR COMPLETIOK-FIT02720
«',F10.3,IX,'DEGREES') FTT02730
C***MODIFY FILTER SIZE SO THAT IT DIVIDES EVENLY INTO 360 DEGREES. FIT02740
NFIL-INT(360./AFIL) FIT02750
C***WRITE MODIFIED ANGULAR FILTER SIZE TO THE DIAGNOSTIC OUTPUT FILE. FTT02760
WRITE(DOOT,126) AFIL FIT02770
126 FORMAT(/,IX,'MODIFIED ANGULAR FILTER SIZE-',F10.3,IX,'DEGREES') FIT02780
1261 CONTINUE . - FIT02790
C***KAKE SURE THAT MAP BOUNDARIES INCLUDE HILL CENTER COORDINATES. FIT02800
XMINl-XHTOP FIT02810
YMIN1-YHTOP FIT02820
XMAX1-XHTOP FIT02830
YMAX1-YHTOP FTT02840
XMIN2-XHTOP FIT02850
YMIN2-YHTOP FIT028CO
XMAX2-XHTOP . FIT02870
YMAX2-YHTOP FIT02880
C***WRTTE THE KILL TOP ELEVATION AND HILL CENTER COORDINATES TO THE ?IT02890
C***DIAGNOSTIC OUTPUT FILE. FIT02900
WRITE(DOUT,130) HTOP,XHTOP,YHTOP FIT02910
130 FORMAT(/,IX,'HILL TOP ELEVATION-',E12.4,/, ' FIT02920
«1X,'HILL CENTER X-COORDINATE-',E12.4,/, . FIT02930
«1X,'HILL CENTER Y-COORDINATE-',E12.4) . FIT02940
C*** • FIT02950
C*** FIT029<0
C SPECIFICATION OF CONTOURS TO BE INPUT FROM THE MASTER FILE USING' FIT02970
C ONE OF 3 METHODS FIT02980
C*** FIT02990
C*** FIT03000
119
-------
C***ASK THE USER TO SPECIFY THE MODE OF CONTOUR SELECTION FROM
C***THE CONTOUR MASTER FILE.
135 WRITE(*,140)
WRITE(*,142)
: 140 FORMAT(//,22X,'SPECIFY CONTOUR SELECTION MODE',/,
422X,'!.) ALL CONTOURS SELECTED',/,
&22X,'2.) SELECT RANGE OF CONTOUR IDs',/,
i22X,'3.) INPUT FILE WITH CONTOUR IDs')
142 FORMAT(/,26X,'CHOICE?(1,2,OR 3) -> *\)
READ(*,'(BN,I3)«,ERR-135) ICMODE
GO TO 150
GO TO 170
GO TO 210
IF(ICMODB.EQ.l)
IF(ICMODE.EQ.2)
IF(ICMODE.EQ.3)
WRITE(*,146)
146 FORMAT(/,IX,'***ERROR***
GO TO 135
FIT03010
FIT03020
FIT03030
FIT03040
FIT03050
FIT03060
FTT03070
FTT03080
FIT03090
FIT03100
FIT03110
FIT03120
FIT03130
FTT03140
MODE SELECTION OUT OF RANGE—TRY AGAIN')FIT03150
FIT03160
1)
C***USE ALL CONTOURS IN THE MASTER FILE(CONTOUR SELECTION MODE
150 WRITE(DOUT,160) MASTER
> 160 FORMAT(/,IX,'ALL CONTOURS IN FILE ',A15,IX,'SELECTED FOR INPUT')
GO TO 300
C***INPUT THE SMALLEST AND LARGEST ID NUMBERS FOR THE GROUP OF
C***CONTOURS(CONTOUR SELECTION MODE NUMBER 2).
I 170 WRITE(*,180)
i 180 FORMAT(/,IX,'INPUT SMALLEST ID NUMBER(1-9999) FOR CONTOUR GROUP ->FIT03240
! * »\) FIT03250
FIT03170
FIT03180
FIT03190
FIT03200
FIT03210
FIT03220
FIT03230
READ(*,'(BN,I4)',ERR~170) ICL
FIT03260
IF(ICL.EQ.O) GO TO 170 FIT03270
1 185 WRITE(*,190) FIT03280
I 190 FORMAT(/,IX,'INPUT LARGEST ID NUMBER(1-9999) FOR CONTOUR GROUP -> FIT03290
*'\) FIT03300
i READ(*,'(BN,I4)',ERR-185) ICU FIT03310
i IF(ICU.EQ.O) GO TO 185 FIT03320
! IF(ICU.GB.ICL) GO TO 195 •• FIT03330
i WRITE(*,191) FIT03340
1 191 FORMAT(/,IX,'***ERROR*** LOWER SERIAL NUMBER GREATER THAN UPPER—FIT03350
! ATRY AGAIN*) FIT03360
: GO TO 170
i 195 CONTINUE
C***WRITE ID RANGE FOR CONTOUR SELECTION TO THE DIAGNOSTIC OUTPUT FILE
I WRITE(DOUT,200) MASTER,ICL,ICU
I 200 FORMAT(/,IX,'CONTOURS SELECTED FROM MASTER FILE ',A15,/,
«1X,'HAVE ID NUMBERS BETWEEN',15,IX,'AND',15)
GO TO 300
C***IHPUT THE NAME OF THE FILE CONTAINING THE CONTOUR ID NUMBERS FOR
C***THE HILL IN QUESTION(CONTOUR SELECTION MODE NUMBER 3).
! 210 WRITE(*,220)
: 220 FORMAT(/,IX,'ENTER CONTOUR ID FILE NAME ->-'\)
I READ(*,'(A)') CON7ILE
I IF(CONFILE.EQ.' *) GO TO 210
C***OPEN CONTOUR ID FILE.
! OPEN(CONIN,FILE-CONFILE,STATUS-'OLD')
C***INPUT ID NUMBERS FROM THE CONTOUR ID FILE.
C***SET COUNTER FOR CONTOUR ID«.
I NdD-1
! 230 CONTINUE
C***READ THE NEXT ID NUMBER.
' READ(CONIN,*,END-270) IDCPK(NCID)
I IfCID-NCID+l
C***CHECK TO SEE IF THE NUMBER OF CONTOURS IS GREATER THAN THE MAXIMUM
C***AMOUNT.
FIT03370
FIT03380
FTT03390
FIT03400
FIT03410
FIT03420
FIT03430
FIT03440
FIT03450
FIT03460
FIT03470
FIT03480
FIT03490
FIT03500
FIT03510
FIT03520
FIT03530
FIT03540
FIT03550
FIT03560
FIT03570
FIT03S80
FIT03590
FIT03600
120
-------
IF(NCID.GT.NCMAX) GO TO 250 FIT03610
GO TO 230 FIT03620
250 WRITE(DOUT,260) NCMAX FIT03630
260 FORMAT(/,IX,•***WARNING***MAXIMUM NUMBER OF CONTOURS(',14,') REAC FIT03640
tHED') FIT03650
C***DETERMINE WHETHER ANY CONTOURS HAVE BEEN REQUESTED. IF NOT, WRITE FIT03660
C***AN ERROR MESSAGE TO BOTH THE DIAGNOSTIC OUTPUT FILE AND THE SCREEN FIT03670
C***AND THEN EXIT THE PROGRAM. FIT03680
270 NCID-NCID-1 FIT03690
IF(NCID.EQ.O) GO TO 1000 FIT03700
WRITE(DOUT,280) NCID,MASTER,IDHILL,HNAME FIT03710
280 FORMAT(/,IX,14,IX,'CONTOURS TO BE SELECTED FROM MASTER FILE ', FIT03720
*A15,/,1X,'FOR HILL',14,'(«,A15,')',//,IX,'IDs REQUESTED:') FIT03730
C***SORT LIST OF CONTOUR IDs IN ASCENDING ORDER. - FIT03740
CALL ISORT(IDCPKfNCID,LPTR) FIT03750
WRITE(DOUT,290) (IDCPK(I),I-1,NCID) FIT03760
290 FORMAT(IX,15) FIT03770
C***CLOSE THE CONTOUR ID FILE. FIT03780
CLOSE(COKIN,STATUS«'KEEP>) FIT03790
300 CONTINUE • FIT03800
C*** FIT03810
C*** FIT03820
C DETERMINE WHETHER A PLOT IS TO BE GENERATED, INPUT PLOT FILE NAME, FIT03830
C AND OPEN THE PLOT FILE. IF PLOT IS REQUESTED, ALSO OPEN A SCRATCH FIT03840
C FILE "PSCRAT". FIT03850
C*** FIT03860
C*** FIT03870
C***ASX WHETHER A PLOT IS TO BE GENERATED. FIRST, INITIALIZE THE PLOT FIT03880
C***FLAG INDICATOR TO CORRESPOND TO A "NO" ANSWER. FIT03890
PFLAG-0 FIT03900
WRITE(*,310) FZT03910
310 FORMAT(/,IX,'PLOT REQUESTED?(Y/N) -> (\) FIT03920
READ(*,'(A)') ANS FIT03930
IF(AHS.EQ.'Y'.OR.ANS.EQ.'y«) PFLAG-1 FIT03940
IF(PFLAG.EQ.O) GO TO 315 ' FIT03950
C***ASX USER TO INPUT THE NAME OF THE PLOT FILE. FIT03960
3101 WRITE(*,311) FIT03970
311 FORMAT(/,IX, 'ENTER PLOT FILE NAME -> «\) .FIT03980
READ(*,'(A)') PFILE • FIT03990
IF(PFILE.EQ.' ') GO TO 3101 FIT04000
C***OPEN THE PLOT FILE AND THE SCRATCH FILE. FIT04010
OPEN (UPL,FILE-PFILE.STATUS-'NEW') FIT04020
OPEN(UPSCR,FILE«'PSCRAT>,STATUS-'NEW') FIT04030
IF(PFLAG.NE.l) GO TO 315 FIT04040
C***WRITE "FITCON" TO THE FIRST RECORD OF THIS PLOT FILE TO INDICATE FIT04050
C***THAT THE PLOT FILE IS BEING GENERATED BY PROGRAM FITCON. FIT04060
WRITE(UPL,3111) . FIT04070
3111 FORMAT('FITCON') FIT04080
C***WRITE THE HILL ID NUMBER AND NAME TO THE PLOT FILE. FIT04090
WRITE (UPL,312) IDHILL,RNAME FIT04100
312 FORHAT(I2,1X,A15) FIT04110
C***WRITE THE HILL CENTER COORDINATES TO THE PLOT FILE. . FIT04120
WRITE(UPL,313) XHTOP,YHTOP FIT04130
313 FORMAT(2E15.4) FIT04140
315 CONTINUE FIT04150
WRZTE(*,316) FIT04160
316 FORMAT(/,IX,'Plaa«a wait...Contour data b«ing procasced',/) FIT04170
C*** FIT04180
C*** . FIT04190
C INPUT AND EDIT CONTOUR DATA. FIT04200
121
-------
C*** FIT04210
C*** FIT04220
C***SET CONTOUR COUNTER. FIT04230
J-l FIT04240
320 CONTINUE FIT04250
C***CHECK WHETHER THE MAXIMUM NUMBER OF CONTOURS HAVE BEEN INPUT. FIT04260
IF(J.OT.NCMAX) GO TO 670 FIT04270
C***INPUT THE ID NUMBER* ELEVATION, NUMBER OF POINTS, AND CONTOUR FIT04280
C***CLOSURE INDICATOR FOR THE NEXT CONTOUR. FIT04290
READ(IN,*,END-700) IDC(J),HCON(J),NPC,CFIAG FTT04300
IF(ICKODE.NE.2) GO TO 340 FIT04310
C***CONTOUR SELECTION MODE 2 FIT04320
C***DETERMINE WHETHER THE CONTOUR ID NUMBER FALLS WITHIN THE BOUNDS FIT04330
C***SPECIFIED BY THE USER. IF NOT, READ DATA FOR ANOTHER CONTOUR FROM FIT04340
C***THE MASTER FILE. FIT04350
IF(IDC(J).LT.ICL.OR.IDC(J).GT.ICU) GO TO 355 FIT04360
GO TO 360 FIT04370
340 IF(ICMODE.NE.3) GO TO 360 FIT04380
C***CONTOUR SELECTION MODE 3 FIT04390
C***DETERM1NE WHETHER THE ID NUMBER FOR THE CONTOUR INPUT FROM THE FIT04400
C***MAST£R FILE MATCHES ONE OF THE SORTED ID NUMBERS INPUT FROM COMPILE.FIT04410
C***IF NOT, READ DATA FOR ANOTHER CONTOUR FROM THE MASTER FILE. FIT04420
DO 350 I-1,NCID FIT04430
C***SINCE IDCPK ARRAY VALUES HAVE BEEN SORTED IN ASCENDING ORDER, THE FIT04440
C***CURRENT ID NUMBER FROM THE MASTER FILE CAN SOMETIMES BE ELIMINATED FIT04450
C***WITHOUT HAVING TO GO THROUGH THE ENTIRE LIST OF IDCPK ARRAY VALUES. FIT04460
IF(IDC(J).LT.IDCPK(I)) GO TO 355 FIT04470
IF(IDC(J).EQ.IDCPK(I)) GO TO 360 FIT04480
350 CONTINUE FIT04490
355 CALL SKIPCN(IH,NPC) FIT04500
GO TO 320 . FIT04510
360 CONTINUE FIT04520
C***CHECK WHJfi'HEK THE CONTOUR ELEVATION IS GREATER THAN THE KILL TOP FIT04530
!C***BLBVATION. IF SO, WRITE AN ERROR MESSAGE AND DISCONTINUE PROCESSING FIT04540
;C***THE CONTOUR. FIT04550
IF(HCON(J).LT.HTOP) GO TO 375 FIT04560
WRITE(*,365) IDC(J) FIT04570
365 FORMAT(/,IX,'Contour ID ',14,IX,'ha« baan rajactad',/,IX, FIT04580
ft'—Saa diagnostic output fila aftar program completion') FIT04590
WRTTE(DOUT,370) IDC(J),HCON(J),HTOP FTT04600
370 FORMAT(/,IX,'***ERROR*** CONTOUR ID',15 ,IX,'DOES NOT RAVE AN ELEFIT04610
tVACTON LESS THEN THE KILL TOP',/,IX,'CONTOUR ELEVATION-»,E12.4, FIT04620
*/,lX,'HILL TOP ELEVATION-',E12.4,/,IX,'CONTOUR WILL NOT BE PROCESSFIT04630
&ED',/)
CALL SKEPCN(IN,HPC)
GO TO 320
C***FIND WHETHER THE CONTOUR HAS AN ELEVATION WHICH IS THE SAME AS A
,C***CONTOUR WHICH HAS BEEN PREVIOUSLY ACCEPTED. IF SO, WRITE AN ERROR
C***MESSAGE AND DISCONTINUE PROCESSING THE CONTOUR. MULTIPLE CONTOURS
C***AT THE SAME ELEVATION MUST BE INPUT AS A SINGLE CONTOUR.
375 IF(J.EQ.l) GO TO 38O
DO 376 JJ-1,JM1
JJK-vTJ •
IT(ABS(HCON(J)-HCON(JJ)).LE.1.0E-15) GO TO 377
376 CONTINUE
GO TO 380
377 WRITE (DOUT, 378) IDC(JJK) ,HCON(J)
378 FORMAT (/, IX, '***ERROR*** PREVIOUSLY ACCEPTED CONTOUR ID1, 15, IX,
ft 'ALSO HAS',/, IX, 'AN ELEVATION OF' , £15. 4, IX, '—CONTOUR REJECTED',
FIT04640
FIT04650
FIT04660
FIT04670
FIT04680
FIT04690
FIT04700
FTT04710
FIT04720
FIT04730
FIT04740
FTT04750
FIT04760
FIT04770
FIT04780
FIT04790
FIT04800
122
-------
&/,lX,'MULTIPLE CONTOURS AT THE SAME ELEVATION MUST BE INPUT AS A SFIT04810
SINGLE CONTOUR') FIT04820
WRITE(*,365) IDC(J) FIT04830
CALL SKIPCN(IN,HPC) . FIT04840
GO TO 320 FIT04850
C***CHEOC WHETHER THE CONTOUR HAS FEWER THAN 3 POINTS. IF SO, WRITE AN FIT04860
C***ERROR MESSAGE AND DISCONTINUE PROCESSING THE CONTOUR. FIT04870
380 IF(NPC.GT.2) GO TO 400 FIT04880
WRITE(*,365) IDC(J) FIT04890
WRITE(DOUT,390) IDC(J),NPC FIT04900
390 FORMAT(/,IX,'***ERROR*** CONTOUR ID', IS,IX,'HAS FEWER THAN 3 POIFIT04910
! WrrS.',/,14X,'CONTOUR HILL NOT BE PROCESSED',/) FZT04920
CALL SKIPCN(IN,NPC) FTT04930
I GO TO 320 FIT04940
C***CHECK WHJgi'EEK THE MAXIMUM NUMBER OF CONTOUR POINTS HAS BEEN EXCEEDEDFIT04950
C***IF SO, WRITE AN ERROR MESSAGE AND DISCONTINUE PROCESSING THE CONTOURFIT04960
400 IF(NPC.LT.NPCMAX) GO TO 420 FTT04970
WRITE(*,36S) IDC(J) FIT04980
WRITE(DOUT,410) IDC(J),NPC,NPCMAX FTT04990
410 FORMAT(/,IX,«***£RROR*** CONTOUR ID',15,IX,'HAS',15,IX,'POINTS.•,FIT05000
*/,14X,'MAXIMUM ALLOWED IS',IS,'. CONTOUR WILL NOT BE PROCESSED.') FITOS010
CALL SKIPCN(IH,HPC)
GO TO 320
C***WRITB THE CONTOUR ELEVATION TO THE DIAGNOSTIC OUTPUT FILE.
420 WRITE(DOUT,425) IDC(J),HCON(J)
425 -FORMAT(/,IX,'CONTOUR ELEVATION FOR CONTOUR ID', 15,IX,'-•,E12.4)
C***INPUT X,Y COORDINATES OF CONTOUR POINTS.
READ(IN,*)
-------
XCONSV(K)-XCON(K) FIT05410
YCONSV(K)-YCON(K) FIT05420
460 CONTINUE FIT05430
IF(MCFLAG.GE.3) GO TO 530 FIT05440
C***P£RFORM EDIT CHECKING FOR A SINGLE CONTOUR. FIT05450
C***FIND THE DISTANCE(DFTOL) FROM THE FIRST TO THE LAST CONTOUR. FIT05460
DFTOL»SQRT((XCON(NPC)-XCOH(l))**2+(YCON(NPC)-YCON(l))**2) FIT05470
C***IF THIS DISTANCE IS EFFECTIVELY ZERO AND THE CONTOUR HAS BEEN FIT05480
C***SPECIFIED AS CLOSED, THEN CONTINUE PROCESSING THE CONTOUR. FIT05490
IF(DFTOL.LT.1.0E-15.AND.CFLAG.EQ.l) GO TO 530 FIT05500
C***IF THIS DISTANCE IS EFFECTIVELY ZERO AND THE CONTOUR HAS BEEN FIT05510
C***SPECIFIED AS OPEN, THEN WRITE A WARNING TO THE DIAGNOSTIC OUTPUT FIT05520
C***FILB AND CONTINUE PROCESSING THE CONTOUR AS IF IT WERE CLOSED. FIT05530
IF(DFTOL.LT.1.0E-15.AND.CFLAG.NE.l) GO TO 510 FTT05540
C***IF THIS DISTANCE IS SIGNIFICANTLY GREATER THAN ZERO AND THE CONTOUR FIT05550
C***HAS BEEN SPECIFIED AS CLOSED, THEN ADD TO THE CONTOUR A FINAL POINT FIT05560
C***WHICH HAS THE SAME COORDINATES AS THE FIRST POINT. IF THE ADDITION FIT05570
C***OF THIS POINT CAUSES THE NUMBER OF CONTOUR POINTS TO EXCEED THE FIT05580
C***HAXIMUM ALLOWABLE, THEN SUBSTITUTE THE FIRST CONTOUR POINT FOR THE FTT05590
C***LAST CONTOUR POINT AND CONTINUE PROCESSING THE CONTOUR AS IF IT WEREFIT05600
C***CLOSED. THE APPROPRIATE WARNINGS ARE WRITTEN TO THE DIAGNOSTIC FIT05610
C***OUTPUT FILE. FIT05620
IF(DFTOL.GB.1.0E-15.AND.CFLAG.EQ.l) GO TO 470 FIT05630
C***IF THIS DISTANCE IS SIGNIFICANTLY GREATER THAN ZERO AND THE CONTOUR FIT05640
C***HAS BEEN SPECIFIED AS OPEN, THEN CALL SUBROUTINE CONCOMP TO ADD FIT05650
C***POINTS TO COMPLETE THE CONTOUR. FIT05660
CALL CONCOMP(XCON,YCON,NPC,NPCMAX,XHTOP,YHTOP,AFIL,NFIL,DOUT) FIT05670
GO TO 330 FIT05680
470 IF(NPC.EQ.NPCHAX) GO TO 490 FTT05690
HPONPC+1 FTT05700
XCON(NPC)-XCON(1) • FTT05710
YCON(NPC)-YCON(1) FTT05720
WRITE(DOUT,480) ' . FTT05730
480 FORMAT(/,IX,'***WARNING***CONTOUR SPECIFIED AS CLOSED WAS FOUND TOFIT05740
t BE OPEN.',/,14X,'ADDED FINAL POINT IS ASSUMED TO BE THE SAME AS TFIT05750
tHE INITIAL POINT.') FIT05760
SO TO 530 FTT05770
490 XCON(NPC)-XCON(1) FTT05780
YCON(NPC)-YCON(1) FTT05790
WRITE(DOUT,500) FIT05800
; 500 FORMAT(/,IX,'***WARNING***CONTOUR SPECIFIED AS CLOSED WAS FOUND TOFIT05810
; t BE OPEN.*,/,14X, 'ADDED FINAL POINT IS ASSUMED TO BE THE SAME AS TFIT05820
; CHE INITIAL POINT',//,IX, '***WARNING***MAXIMUM NUMBER OF CONTOUR POFIT05830
I 4INTS EXCEEDED IN THE CLOSING OPERATION.',/,14X,'FINAL POINT IS REPFIT05840
I 4LACED BY THE INITIAL POINT.') FIT05850
! GO TO 530 FTT05860
| 510 WRITE (DOUT, 520) . FTT05870
| 520 FORMAT(/,IX,'***WARNING***CONTOUR SPECIFIED AS OPEN WAS FOUND TO BFIT05880
! tE CLOSED1) FTT05890
! 530 CONTINUE FIT05900
! C***WRITE THE EDITED NUMBER OF CONTOUR POINTS TO THE DIAGNOSTIC OUTPUT FIT05910
i C***FILE. . FIT05920
! WRITE(DOUT,531) IDC(J),NPC FIT05930
531 FORMAT(/,IX,'MODIFIED NUMBER OF POINTS FOR CONTOUR ID',15,IX,'-•, FTT05940
I *I5) FIT05950
! C***WRITE THE EDITED CONTOUR POINT COORDINATES TO THE DIAGNOSTIC OUTPUT FIT05960
C***FILE. . FIT05970
! WRITE(DOUT,532) IDC(J) FIT05980
i 532 FORMAT(/,IX,'X-Y-COORDINATES(EDITED) FOR CONTOUR ID',15,/) FIT05990
WRITE(DOUT,450) (XCON(K),YCON(K),X-1,NPC) FIT06000
124
-------
c***
c***
C CALCULATE THE AREA AND CENTER OF MASS FOR THE INPUT CONTOUR.
FIT06010
FIT06020
FIT06030
FIT06040
FIT06050
CALL ARCH(XCON,YCON,AR,XC,YC,NPC) FIT06060
AREA-ABS(AR) FIT06070
C***DETERMINE WHETHER THE CALCULATED AREA OF THE CONTOUR IS EFFECTIVELY FIT06080
C***Z£RO. IF SO, WRITE AN ERROR MESSAGE AND DISCONTINUE PROCESSING THE FIT06090
C***CONTOUR. FIT06100
IF(AREA.GT.1.0E-15) GO TO 550 FIT06110
WRITE(*,365) IDC(J) FIT06120
WRITE(DOUT,540) FIT06130
540 FORMAT(/,IX,'AREA FOUND TO BE EFFECTIVELY ZERO—CONTOUR REJECTED1)FITO6140
GO TO 320 FITO6150
550 CONTINUE FIT06160
C***CALCULATE THE MAXIMUM RADIUS OF GYRATION AND THE ASSOCIATED MINOR FIT06170
C***AXIS ORIENTATION FOR THE CONTOUR. FIT06180
CALL SMOMNT(XCON,YCON,AR,NSLOPE,SN,CN,ANGLE,NPC, FITO6190
iXC,YC,RG,RGRAT,ORENT,ISMFLG) FIT06200
C***DETERMINE WHETHER A REAL VALUE FOR THE RADIUS OF GYRATION HAS BEEN FIT06210
C***CALCULATED FOR THE CONTOUR. IF NOT, WRITE AN ERROR MESSAGE AND FIT06220
C***DISCONTINUE PROCESSING THE CONTOUR. FIT06230
IF(ISMFLG.EQ.O) GO TO 555 FIT06240
WRITE(*,365) IDC(J) FIT06250
WRITE(DOUT,551) FIT06260
551 FORMAT(/,IX,'CONTOUR REJECTED—A REAL VALUE FOR THE RADIUS OF GYRAFIT06270
4TION COULD HOT BE*,/,IX,'COMPUTED. THIS CAN OCCUR IF THE CONTOUR IFIT06280
63 VERY TORTUOUS AND TOO FEW POINTS WERE',/,IX,' USED IN ITS DIGITIFIT06290
4ZATION OR IF A VERY TORTUOUS CONTOUR HAS BEEN INPUT AS',/,lX,'INCOFIT06300
tMPLETE EVEN IF A SUFFICIENT NUMBER OF POINTS HAVE BEEN USED IN THEFIT06310
*',/,IX,'DIGITIZATION PROCESS.',//,IX,'SOLUTION—MANUALLY COMPLETE FIT06320
&AND/OR REDIGITIZE THE CONTOUR')
GO TO 320
555 CONTINUE
XCM(J)-XC
i YCM(J)-YC
C***WRITE THE CALCULATED CONTOUR AREA AND CENTROID COORDINATES TO THE
C***DIACNOSTIC OUTPUT FILE.
WRITE(DOUT,5«0) AREA,XCH(J),YCM(J)
; 5«0 FORMAT(/,IX,'CONTOUR AREA-',E12.4,/,
: *1X,'X-COORDINATE OF CONTOUR CENTROID-',E12.4,/,
i ilX,'Y-COORDINATE OF CONTOUR CENTROID-',E12.4)
C***EDIT CHECKS HAVE BEEN COMPLETED. CONTOUR HAS BEEN ACCEPTED FOR
C***PROCESSING.
17 A PLOT HAS BEEN REQUESTED, WRITE THE CONTOUR COORDINATES(BOTH
UNEDITED AND EDITED) TO THE SCRATCH FILE •PSCRAT" AND UPDATE THE
PLOT BOUNDARIES TO REFLECT THE BOUNDARIES OF THE NEWLY INPUT
CONTOUR.
1C***
c***
c***
IT(PFLAG.EQ.O) GO TO 575
WRITE(UPSCR,570) NPCSV,HCON(J)
570 FORMAT(110,£15.4)
DO 572 K-l,NPCSV
WRITE(UPSCR,571) XCONSV(X),YCONSV(K)
571 FORMAT(2E15.4)
IF(XCONSV(X).GT.XMAX1) XMAX1-XCONSV(K)
FIT06330
FIT06344
FIT06350
FIT0636e
FIT06370
FTT06380
FIT06390
FITO6400
FITO6410
FIT06420
FIT06430
FIT06440
FIT06450
FTT06460
FIT06470
FIT06480
FIT06490
FIT06500
FIT06510
FIT06520
FIT06530
FIT06540
FIT06550
FTT06560
FIT06570
FIT06580
FIT06590
FIT06600
125
-------
XMIN1-XCONSV(K)
YMAXl-YCONSV(K)
YMINl-YCONSV(K)
IF(XCONSV(K).LT.XMIN1)
IF(YCONSV(K).GT.YMAX1)
IF(YCONSV(K).LT.YMIN1)
572 CONTINUE
WRITE(UPSCR,570) NPC,HCON(J)
DO 574 K-1,NPC
WRITE(UPSCR,571) XCON(K),YCON(K)
IF(XCON(K).GT.XMAX2) XMAX2-XCON(K)
.LT.XHIN2)
.GT.YMAX2)
.LT.YMIN2)
IF(XCON(K)
IF(YCON(K)
IF(YCON(K)
CONTIHUE
XHIN2-XCON(K)
YMAX2-YCON(K)
YKtN2-YCON(K)
574
575 CONTIHUE
C***
c***
C COMPOTE THE PARAMETERS FOR THE ELLIPTICAL REPRESENTATION OF THE
C CONTOUR.
C***
c***
OREN(J)-ORENT
C***CALCULATE THE SEMI-MAJOR AXIS LENGTH FOR THE EQUIVALENT ELLIPSE
C***USING THE RELATIONSHIP, FOR AN ACTUAL ELLIPSE, BETWEEN THE
C***SEMI-MAJOR AXIS LENGTH AND THE RADIUS OF GYRATION ABOUT AN AXIS
C***WHICH COINCIDES WITH THE SEMI-MINOR AXIS OF THE ELLIPSE.
A(J)-2.*RG
C***CALCULATE THE SEMI-MINOR AXIS LENGTH FOR THE EQUIVALENT ELLIPSE
C***USING THE FORMULA FOR THE AREA OF AN ELLIPSE AND THE PREVIOUSLY
C***DETERMINED VALUE FOR THE SEMI-MAJOR AXIS LENGTH.
B(J)-AREA/(PI*A(J))
C***DETERMINE WHETHER THE CONTOUR SHOULD BE CONSIDERED CIRCULAR
C***FXRST TEST FOR CIRCULAR CONTOUR—CALCULATED SEMI-MINOR AXIS
C***LENGTH GREATER THAN OR EQUAL TO SEMI-MAJOR AXIS LENGTH.
IF(A(J).GT.B(J)) GO TO 390
WRITE(DOUT,580)
980 FORMAT(/,IX,'CALCULATED ELLIPSE SEMI-MINOR AXIS LENGTH WAS FOUND',
t' TO BE GREATER THAN*,/,IX,'OR EQUAL TO THE CALCULATED SEMI-MAJOR'
*,' AXIS LENGTH—CONTOUR ASSUMED TO BE CIRCULAR')
GO TO 610
C***SECOND TEST FOR CIRCULAR CONTOUR—DETERMINE WHETHER THE RELATIVE
C***DIFFERENCE BETWEEN THE MAXIMUM AND MINIMUM RADII OF GYRATION FOR
C***THE CONTOUR IS LESS THAN 1 PERCENT.
990 IF(RGRAT.GT.O.Ol) GO TO 620
WRITE(DOUT,600)
600 FORMAT(/,IX,'THE RELATIVE DIFFERENCE BETWEEN THE MAXIMUM AND',
*' MINIMUM RADII OF GYRATION',/,IX,'FOR THE CONTOUR IS LESS THAN',
t* 1 PERCENT—CONTOUR ASSUMED TO BE CIRCULAR')
C***SET BOTH THE SEMI-MAJOR AND SEMI-MINOR AXIS LENGTHS EQUAL TO THE
C***RADIUS OF A CIRCLE WITH AREA EQUAL TO AREA.
610 RAD-SQRT(AREA/PI)
C***THE ECCENTRICITY OF A CIRCLE IS ZERO.
ECC(J)-0.
A(J)-RAD
B(J)«RAD
GO TO 630
C***CALCULATE THE ECCENTRICITY OF THE ELLIPSE REPRESENTING THE CONTOUR.
620 ECC(J)-SQRT(A(J)**2-B(J)**2)/A(J)
C***WRITE ELLIPSE FIT PARAMETERS TO THE DIAGNOSTIC OUTPUT FILE.
630 WRITE(DOUT,640) IDC(J)
640 FORMAT(/,IX,'ELLIPSE PARAMETERS FOR CONTOUR ID',I5,/)
WRITE(DOUT,650) A(J),B(J),ECC(J),OREN(J)
FIT06610
FIT06620
FIT06630
FIT06640
FIT06650
FIT06660
FIT06670
FITO«680
FIT06690
FIT06700
FIT06710
FIT06720
FIT06730
FTT06740
FIT06750
FIT06760
FIT06770
FIT06780
FIT06790
FIT06800
FIT06810
FIT06820
FIT06830
FIT06840
FIT06850
FIT06860
FIT06870
FIT06880
FIT06890
FIT06900
FIT06910
FIT06920
FIT06930
FIT06940
FIT06950
FIT06960
FIT06970
FIT06980
FIT06990
FIT07000
FIT07010
FIT07020
FIT07030
FIT07040
FIT07050
FIT07060
FIT07070
FIT07080
FIT07090
FIT07100
FIT07110
FIT07120
FIT07130
FIT07140
FIT07150
FIT07160
FIT07170
FIT07180
FIT07190
FIT07200
126
-------
650 FORMAT(IX,'SEMI-MAJOR AXIS LENGTH-',E12.4,/, FIT07210
«1X,'SEMI-MINOR AXIS LENGTH-',E12.4,/, FIT07220
*1X,'ELLIPSE ECCENTRICITY-',E12.4,/, FIT07230
«1X,'ORIENTATION OF SEMI-MINOR AXIS WITH RESPECT TO THE1, FIT07240
«/,IX,'POSITIVE X-AXIS-',F«.2,IX,•DEGREES•) FIT07250
C***UPDATE THE CONTOUR COUNTER AND READ DATA FOR A NEW CONTOUR FROM THE FIT07260
C***MAST£R FILE. FIT07270
HRITE(*,660) IDC(J) FTT07280
660 FORMAT(/,IX,'Contour ID ',14,IX,'has been accepted') FIT07290
J-J+1 FIT07300
•• GO TO 320 'FIT07310
670 -WRITE(DOUT,260) FIT07320
C***END OF CONTOUR MASTER FILE REACHED • FIT07330
! 700 MOJ-1 FIT07340
C***CLOSE THE MASTER FILE. FIT07350
CLOSE(IN,STATUS-'KEEP') FIT07360
C***IF THE CONTOUR ID NUMBERS(FOR CONTOUR SELECTION FROM THE MASTER FIT07370
C***FILE) WERE INPUT FROM CONFILE, CHECK WHETHER THE NUMBER OF CONTOURS FIT07380
C***REQUESTED MATCHES THE NUMBER ACTUALLY SELECTED FROM THE MASTER FILE.FIT07390
C***IF NOT, WRITE A WARNING MESSAGE Td THE DIAGNOSTIC OUTPUT FILE. FTT07400
IF(ICMODE.EQ.3.AND.NC.NE.NCID) WRITE(DOUT,710) ' FIT07410
710 FORMAT(/,IX,•***WARNING***NUMBER OF CONTOURS SELECTED FROM THE1, FIT07420
: «' MASTER FILE DOES NOT',/,14X,'MATCH THE NUMBER REQUESTED') FIT07430
C***CHECK WHETHER ANY CONTOURS HAVE BEEN SELECTED FROM THE MASTER FILE. FIT07440
C***IF NOT, WRITE AN ERROR MESSAGE BOTH TO THE DIAGNOSTIC OUTPUT FILE FIT07450
C***AND THE SCREEN AND THEN EXIT THE PROGRAM. FIT07460
I IF(NC.EQ.O) GO TO 1010 FIT07470
C*** . FIT07480
C*** FTT07490
C***WRITE THE OUTPUT FILES FOR SUBSEQUENT PROCESSING BY THE PLOT PROGRAMFIT07500
C***AND THE CRITICAL HEIGHT ANALYSIS PROGRAM(HCRIT). . FIT07510
'C*** FIT07520
C*** . FIT07530
C***SORT THE ID NUMBERS FOR THE CONTOURS WHICH WERE FINALLY SELECTED. FIT07540
CALL ISORT(IDC,NC,LPTR) FIT07550
C***CHECK WHETHER PLOT HAS BEEN REQUESTED. IF SO, WRITE TO THE PLOT FIT07560
C***FILE THE INFORMATION NECESSARY TO SUBSEQUENTLY PLOT THE INPUT FIT07570
C***DIGITIZED CONTOURS. FIT07580
1 IF(PFLAG.EQ.O) GO TO 770 FIT07590
C***REWIND THE SCRATCH FILE. FIT07600
REWIND UPSCR FIT07610
C***WRITE THE NUMBER OF*CONTOURS TO THE PLOT FILE. FIT07620
WRITE(UPL,720} NC . FIT07630
720 FORMAT(110) FIT07640
C***WRITE THE SORTED CONTOUR ID NUMBERS TO THE PLOT FILE. FIT07650
C***NOTE: THE CONTOUR IDs ARE SORTED ONLY FOR SUBSEQUENT ID CHECKS IN FIT07660
C***THE PLOT PROGRAM. THE DIGITIZED CONTOURS INPUT TO THE PLOT PROGRAM FIT07670
C***DO NOT RAVE TO BE SORTED. FIT07680
WRITE(UPL, 730) (IDC(J),J-1,NC) FIT07690
730 FORMAT(IIO) FIT07700
C***WRITE THE PLOT BOUNDARY LIMITS FOR BOTH UNEDITED AND EDITED FIT07710
C***CONTOURS TO THE PLOT FILE. . FIT07720
WRITE(UPL,740) XMIN1,XMAX1,YMIN1,YMAX1 FIT07730
740 FORMAT(4E15.4) FTT07740
WRITE(UPL,740) XMIN2,XMAX2,YMIN2,YMAX2 FIT07750
C***TRANSFER THE DIGITIZED CONTOUR COORDINATES FROM THE SCRATCH FILE FIT07760
C***TO THE PLOT FILE. FIT07770
NCT2-2*NC FIT07780
DO 760 J-1,NCT2 FIT07790
READ(UPSCR,570) NPC.HCONT • FIT07800
127
-------
WRITE (UK,, 570) HPC,HCONT FIT07810
DO 750 K-1,NPC FIT07820
READ (UPSCR ,-571) XCOH(K),YCON(K) FIT07830
WRITE (UPL, 571) XCOK(K),YCOH(K) PTT07840
750 CONTINUE FIT07850
760 CONTINUE ' FIT07860
C***CLOSE AND DELETE THE SCRATCH- FILE.' nT07870
CLOSE(UPSCR,STATUS-•DELETE•) FIT07 880
C***OP£N THE OUTPUT FILE FOR THE CRITICAL HEIGHT ANALYSIS PROGRAM. FIT07890
770 CONTINUE FIT07900
775 WRITE(*,780) FIT07910
780 FORMAT(/,IX,'ENTER FILE NAME FOR FITTED CONTOUR OUTPUT -> f\) FIT07920
READ(*,f(A)1) COUTFILE FIT07930
IF(COUTFILB.EQ.' ') GO TO 775 FIT07940
OPEN(GOUT,FILE-CODTFILE,STATUS-'NEW') FIT07950
C***WRITE THE HILL ID NUMBER AND NAME TO THE FITTED CONTOUR OUTPUT FIT07960
C***FILE. FIT07970
WRITE(COUT,312) IDHILL,HNAME FIT07980
C***WRITE THE ACTUAL HILL TOP ELEVATION TO THE FITTED CONTOUR OUTPUT FIT07990
C***FILE. FIT08000
WRITE(COUT,790) HTOP FIT08010
790 FORMAT(£15.4) FIT08020
C***WRITE THE NUMBER OF CONTOURS TO THE FITTED CONTOUR OUTPUT FILE. FIT08030
WRITE(COUT,720) NC FIT08040
C***WRITE THE SORTED CONTOUR IDs TO THE FITTED CONTOUR OUTPUT FILE.. FIT08050
C***NOTB: CONTOUR IDs ARE SORTED ONLY FOR SUBSEQUENT ID CHECKS IN FIT08060
C***THE PLOT PROGRAM. FITTED CONTOUR PARAMETERS DO NOT HAVE TO BE FIT08070
C***SORTED FOR INPUT TO THE CRITICAL HEIGHT ANALYSIS PROGRAM. FIT08080
WRITE(COUT,730) (IDC(J),J-1,NC) FIT08090
C***WRITE THE CONTOUR FIT PARAMETERS TO THE FITTED CONTOUR OUTPUT FIT08100
C***FILE. FIT08110
DO 810 J-1,NC • FIT08120
WRITE(COUT,800) HCON(J) ,XCM(J) ,YCH(J) ,A(J) ,B(J) ,ECC(J) ,OREN(J) FIT08130
800 FORMAT(7E15.4) FIT08140
810 CONTINUE . . - FIT08150
IF(PFLAG.EQ.O) GO TO 840 FTT08160
C***WRITE THE CONTOUR FIT PARAMETERS TO THE PLOT FILE. . FIT08170
DO 830 J-1,NC FIT08180
WRITE (UPL,820) XCM(J),YCH(J)/A(J),B(J),OREN(J) FIT08190
820 FORMAT(5E15.4) FIT08200
830 CONTINUE FIT08210
840 CONTINUE FIT08220
C***ANALYSIS COMPLETED—EXIT PROGRAM. FIT08230
GO TO 2000 FIT08240
1000 WRITE(DOUT,1005) FIT08250
1005 FORMAT(/,IX,'***ERROR*** NO CONTOURS WERE REQUESTED—EXIT PROGRAMFIT08260
ft1) . FIT08270
WRITE(*,1005) FIT08280
GO TO 2000 FIT08290
1010 WRITE(DOUT,1015) FIT08300
1015 FORMAT(/,IX,•***ERROR*** NO CONTOURS SELECTED FROM MASTER FILE—EFIT08310
ftXIT PROGRAM') • . FIT08320
WRITE(*,1015)
C***DELETE SCRATCH FILE AND PLOT FILE.
CLOSE(UPSCR,STATUS*'DELETE')
CLOSE(UPL,STATUS-'DELETE')
2000 CONTINUE
STOP
END
FIT08330
FIT08340
FIT08350
FIT08360
FIT08370
FIT08380
FIT08390
128
-------
; SUBROUTINE AHCM(XCON,YCON,AR,XC,YC,NPC)
C***SUBROUTINE TO CALCULATE THE AREA AND CENTROID X-Y COORDINATES
C***FOR AN INPUT CONTOUR
DIMENSION- XCON(1),YCON(1)
AR-0.
XOO.
YC-0.
NPCMl-NPC-1
DO 100 K-1,NPCM1
AR-AR+(YCON(K+l)+YCON(K))*(XCON(K+l)-XCON(K))/2.
XOXC+0. 5* (XCON(K+l) *YCON(K) -XCON(K) *YCON(K+l)) *
«(XCON(K+1)+XCON(K) ) + (YCON (X+l)- YCON (K)) *
4(XCON(K+1)**2+XCON(K)*XCON(K+1)+XCON(K)**2)/3.
YC-YC+0.5* (YCON(K+1) *XCON(K) -YCON(K) *XCON(K+1)) *
«(YCON(K+1)+YCON(K))+(XCON(K+1)-XCON(K))*
«(YCON(K+1) **2+YCON(K) *YCON(K+1)+YCON(K) **2)/3.
XOO CONTINUE
C***CLOSE CONTOUR FOR PURPOSES OF CALCULATING THE AREA AND CENTROID.
C***THIS IS REQUIRED FOR USE OF THE SUBROUTINE BY SUBROUTINE CONCOMP,
C***THE CONTOUR COMPLETION SUBROUTINE.
AR-AR*(YCON(1)+YCON(NPC))*(XCON(1)-XCON(NPC))/2.
XOXC+0 . 5* (XCON (1) *YCON (NPC) -XCON (NPC) * YCON (1)) *
&(XCON(1)+XCON(NPC))+(YCON(1)-YCON(NPC))*
« (XCON(l) **2+XCON(NPC) *XCON(1) -t-XCON(NPC) **2)/3.
YC>YC+0.5*(YCON(1)*XCON(NPC)-YCON(NPC)*XCON(1))*
4(YCON(l)-l-YCON(NPC))-l-(XCON(l)-XCON(NPC))*
C**
«(YCON(1)**2+YCON(NPC)*YCON(1)+YCON(NPC)**2)/3.
'CHECK FOR ZERO AREA CONTOUR
IF(ABS(AR).LT.1.0E-15) RETURN
XOXC/AR
YC—YC/AR
RETURN
END
ARC00010
ARC00020
ARC00030
ARC00040
ARC00050
ARC00060
ARC00070
ARC00080
ARC00090
ARC00100
ARC00110
ARC00120
ARC00130
ARC00140
ARC00150
ARC00160
ARC00170
ARC00180
ARC00190
ARC00200
ARC00210
ARC00220
ARC00230
ARC00240
ARC00250
ARC00260
ARC00270
ARC00280
ARC00290
ARC00300
ARC00310
ARC00320
ARC00330
129
-------
SUBROUTINE CONCOMP(XCON,YCON,NPC,NPCMAX,XHTOP,YHTOP,AFIL,NFIL, CCP00010
6DOUT) . CCP00020
C***THIS SUBROUTINE COMPLETES A CONTOUR WHICH HAS BEEN INPUT FROM THE CCP00030
C***CONTOUR MASTER FILE AS AN INCOMPLETE CONTOUR. THE FIRST STEP IN THISCCP00040
C* "COMPLETION PROCESS IS THE ELIMINATION OF THOSE POINTS WHICH LIE IN CCP00050
C***THE SAME SECTOR AS (1)AN ACTUAL CONTOUR POINT WHICH IS CLOSER TO CCP00060
C***THE HILL CENTER OR (2)A SEGMENT OF A LINE CONNECTING ADJACENT POINTSCCP00070
C***WHICH IS CONTAINED WITHIN THE SECTOR AND WHOSE APPROXIMATE MIDPOINT CCP00080
C***IS CLOSER TO THE HILL CENTER THAN THE POINT IN QUESTION.
C***
ic***
C GLOSSARY OF TERMS
C***
c***
C
C
C
C
C
C
C
C
C
C
C
C
C
-------
DO 1 ISEO1,NTIL CCP00670
DXSTM(XSEC)-0. CCP00680
IR(ISEC)-0 - CCP00690
1 CONTINUE CCP00700
NCOUNT-0 CCP00710
N7ILM-FLOAT(NTIL)/2.+0.500001 CCP00720
XP-XCON(l) CCP00730
YP-YCON(l) CCP00740
C***CALCULATE THE SECTOR AND DISTANCE FROM THE HILL CENTER FOR THE FIRSTCCP00750
C***CONTOUR POINT. CCP00760
CALL VECTOR(XHTOP,YHTOP,XP,YP,ANGLE,DX,DY) CCP00770
ISOLD-ANGLE/AFIL CCP00780
IF(ISOLD.LT.l) ISOLD-1 CCP00790
IF(ISOLD.GT.NFIL) ISOLD-NFIL CCP00800
DOLD-SQRT((XP-XHTOP)**2+(YP-YHTOP)**2) CCP00810
C***CHOOSE THE CLOSEST POINT TO THE HILL CENTER LOCATION FOR EACH CCP00820
C***SECTOR OF ANGULAR WIDTH AFIL. CCP00830
DO 3 K-1,NPC CCP00840
XP-XCON(K) CCP00850
YP-YCON(K) CCP00860
DIST-SQRT((XP-XHTOP)**2+(YP-YHTOP)**2) CCP00870
CALL VECTOR(XHTOP,YHTOP,XP,YP,ANGLE,DX,DY) CCP00880
ISEC-ANGLE/AFIL CCP00890
IF(ISEC.LT.l) ISE01 CCP00900
IF(ISEC.GT.NFIL) ISEONFIL CCP00910
C***DETERMINE WHETHER A CONTOUR POINT OR SEGMENT HAS ALREADY APPEARED INCCP00920
C***THIS ANGULAR SECTOR. CCP00930
IF(IR(ISEC).NE.O) GO TO 2 CCP00940
NCOUNT-NCOUNT+1 . CCP00950
C***ACCEPT THE POINT(XCON(X),YCON(X)) AND INITIALIZE THE ARRAYS IR AND CCP00960
C***DISTM. CCP00970
XCONS(NCOUNT)-XP CCP00980
YCONS(NCOUNT)-YP CCP00990
IR(ISEC)-NCOONT CCP01000
DISTM(ISEC)-SQRT((XP-XHTOP)**2+(YP-YHTOP)**2) CCP01010
GO TO 1000 CCP01020
C***D£TERMINE WHETHER THE DISTANCE FROM THE HILL CENTER TO THE POINT IS CCP01030
C***LESS THAN THE CURRENT MINIMUM DISTANCE FOR THIS SECTOR. CCP01040
2 CONTINUE CCP01050
IF(DIST.GE.DISTM(ISEC)) GO TO 1000 CCP01060
C***REINITIALIZE DISTM(ISEC) TO CORRESPOND TO THE DISTANCE FROM THE CCP01070
C«**HILL CENTER FOR THE CONTOUR POINT IN QUESTION. CCP01080
DISTM(ISEC)-DIST CCP01090
C***TEMPORARILY SAVE THE PREVIOUS VALUE OF IR(SEC). CCP01100
IROLD-IR(ISEC) . CCP01110
IF(IROLD.EQ.9999) GO TO 201 CCP01120
C***FLAG THE X-COORDINATB OF THE PREVIOUSLY CLOSEST CONTOUR POINT CCP01130
C***(IROLD) FOR LATER ELIMINATION OF THE POINT. CCP01140
XCONS(IROLD)-1.0E+13 CCP01150
201 NCOUNT-NCOUNT+1 CCP01160
C***REINITIALIZE THE IR ARRAY TO CORRESPOND TO THE NUMBER OF THE CONTOURCCP01170
C***POINT IN QUESTION. • CCP01180
IR(ISEC)-NCOUNT CCP01190
C***ACCEPT THE CONTOUR POINT NCOUNT. CCP01200
XCONS(NCOUNT)-XP CCPO1210
YCONS(NCOUNT)-YP CCPO1220
C***HANDLE SECTORS BETWEEN THOSE SECTORS OCCUPIED BY THE CURRENT AND CCP01230
C***PREVIOUS CONTOUR POINTS. CCPO1240
C***DETERMINE HOW MANY SECTORS RAVE BEEN CROSSED. IF MORE THAN ONE HAS CCPO1250
C***BEEN CROSSED, THEN DEAL WITH "PSEUDOPOINTS" WHICH OCCUPY SECTORS CCP01260
C***BETWEEN THE CURRENT AND PREVIOUS CONTOUR POINTS. CCP01270
1000 ITEST-IABS(ISEC-ISOLD) CCPO1280
C***DETERMINE WHETHER MORE THAN ONE SECTOR HAS BEEN CROSSED. . CCP01290
IF(ITEST.LE.l.OR.ITEST.EQ.NFIL-l) GO TO 14QO CCPO1300
C***FOUR CASES: CCP01310
C***(l) ISEOISOLD; ISEC-ISOLO-MFILM CCP01320
131
-------
C***(2) ISEOISOUD; ISEC-ISOLDISEC; ISOLD-ISEO-NFILM CCP01340
C***(4) ISOLD>ISEC; ISOLD-ISEC
-------
IF(IR(I).EQ.9999) GO TO 1202 CCP01990
IROLD-IR(I) CCP02000
XCONS(IROLD)-1.0E+15 . CCP02010
1201 HI (I)-9999 CCP02020
1202 DISTM(I)-DTEST CCP02030
COUNT-COUNT+1. CCP02040
1210 CONTINUE CCP02050
GO TO 140O CCP02060
C***CASE II • CCP02070
1300 OIF-FLOAT(NFIL-(ISEC-ISOLD)) CCP02080
COUNT-1. • CCP02090
I C***LOOP MUST BE BROKEN INTO 2 PARTS. CCP02100
C***PART ONE CCP02110
DO 13SO I-ISOLD-1,1,-1 CCP02120
DTEST»DOLD+(COONT/DIF)*(DIST-DOLD) CCP02130
IP(DTEST.GT.DISTH(I).AND.IR(I).NE.O) GOTO 1350 CCP02140
IF(IR(I).EQ.O) GO TO 1311 CCP02150
IF(IR(I).EQ.9999) GO TO 1312 CCP02160
IROLD-IR(I) CCP02170
i XCONS(IROLD)-1.0E+1S CCP02180
1311 HI (I)-9999 CCP02190
i 1312 BISTM(I)-DTEST CCP02200
COUNT-COUNT+1. CCP02210
1350 CONTINUE CCP02220
:c***PART TWO CCP02230
DO 1360 I-NFIL,ISEC+1,-1 CCP02240
DTEST-DOLD+(COUNT/DIF)*(DIST-DOLD) CCP02250
I IF(DTEST.GT.DISTH(I).AND.IR(I).NE.O) GO TO 1360 CCP02260
IF(IR(I).EQ.O) GOTO 1351 CCP02270
IF(IR(I).EQ.9999) GO TO 1352 CCP02280
: IROLD-IR(I) CCP02290
XCONS(IROLD)-1.0E+1S CCP02300
j 1351 IR(I)-9999 CCF02310
! 1352 DISTH(I)-DTEST CCP02320
I COUNT-COUNT+1. CCP02330
' 1360 CONTINUE CCP02340
' C***SAVE THE SECTOR NUMBER AND DISTANCE FOR COMPARISON WITH THE NEXT CCP02350
C***POINT. CCP02360
1400 ISOLD-ISEC . CCP02370
DOLD-DIST CCP02380
3 CONTINUE * CCP02390
NPC-1 CCP02400
DO 4 K-1,NCOONT CCP02410
IF(XCONS(K).GT.1.0E+14) GO TO 4 ' CCP02420
XCON(NPC)-XCONS(K) CCP02430
TCON(NPC)-TCONS(K) CCP02440
NPC-HPC+l CCP02450
4 CONTINUE CCP02460
; NPOHPC-1 CCP02470
! 45 CONTINUE CCP02480
\ C***CKLL SUBROUTINE ARCM TO DETERMINE THE AREA OF THE INCOMPLETE CCP02490
: C***CONTOUR. IF THE AREA IS NEGATIVE, THEN THE CONTOUR POINTS HAVE CCP02500
C***BEEN INPUT IN A COUNTERCLOCKWISE SENSE. IF THE AREA IS POSITIVE, CCP02510
O**THEN THE CONTOUR POINTS HAVE BEEN INPUT IN A CLOCKWISE SENSE. CCP02520
C***THIS INFORMATION IS REQUIRED BY THE CONTOUR COMPLETION ALGORITHM. CCP02530
C***THE X AND Y COORDINATES OF THE INCOMPLETE CONTOUR CENTER OF CCP02540
I C***MASS(XCM,YCM) ARE NOT USED BY THE CONTOUR COMPLETION ALGORITHM. CCP02550
! C***FIRST ADD THE POINT (XHTOP,YHTOP) TO THE CONTOUR ONLY FOR THE CCP02560
C***PURPOSE OF THIS DIRECTION DETERMINATION. CCP02570
! METH-0 CCP02580
XREF-XHTOP CCP02590
YREF-YHTOP CCP02600
NPCPl-HPC-t-1 CCP02610
C***CHECK WHETHER CONTOUR COMPLETION WILL CAUSE THE NUMBER OF CONTOUR CCP02620
:O**POINTS TO EXCEED THE MAXIMUM. IF SO, SET THE COORDINATES OF THE CCP02630
! C***FINAL POINT EQUAL TO THOSE OF THE INITIAL POINT AND PRINT A WARNING CCP02640
133
-------
C***MESSAGE. CCP02650
IF(NPCPl.LT.NPCMAX) GO TO 5 CCP02660
XCON(NPC)-XCON(1) CCP02670
YCON(NPC)-YCON(1) • CCP02680
WRITE(DOUT,20) NPCMAX CCP02690
RETURN CCP02700
5 XCON(NPCP1)-XHTOP CCP02710
YCON(NPCP1)-YHTOP CCP02720
CALL ARCM(XCON,YCON,AREA,XCMfYCM,MPCPl) CCP02730
6 XP-XCON(l) CCP02740
YP-YCON(l) CCP02750
C***FIND THE HEADING AND X,Y COMPONENTS OF THE VECTOR FROM THE HILL CCP02760
C***TOP X,Y POINT TO THE FIRST CONTOUR POINT. CCP02770
CALL VECTOR(XREF,YREF,XP,YP,ANGLE,DX,DY) CCP02780
IF(AREA.LT.O.) ANG2-ANGLE CCP02790
IF(AREA.GE.O.) ANG1-ANGLE CCP02800
XP-XCON(NPC) CCP02810
YP-YCON(NPC) CCP02820
C***FIND THE HEADING AND X,Y COMPONENTS OF THE VECTOR FROM THE HILL CCP02830
C***TOP X,Y POINT TO THE LAST CONTOUR POINT. CCP02840
CALL VECTOR(XREF,YREF,XP,YP,ANGLE,DX,DY) CCP02850
IF(AREA.LT.O.) ANG1-ANGLE CCP02860
IF(AREA.GE.O.) ANG2-ANGLE CCP02870
IF(METH.EQ.l) GO TO 7 CCP02880
ADIF-ANG2-ANG1 CCP02890
IF(ADIF.LT.O.) ADIF-360.+ADIF CCP02900
C***IF THE ANGULAR DIFFERENCE BETWEEN THE VECTORS IS LESS THAN CCP02910
C***90 DEGREES, USE THE CENTROID OF THE CONTOUR FOR THE REFLECTION CCP02920
C***POINT INSTEAD OF THE HILL CENTER. CCP02930
IF(ADIF.GT.90.) GO TO 7 CCP02940
CALL ARCM(XCON,YCON,AREA,XCM,YCM,NPC) CCP02950
XREF-XCM CCP02960
YREF-YCM CCP02970
METH-1 CCP02980
GO TO 6 CCP02990
7 CONTINUE ' CCP03000
C***SAVE THE NUMBER OF ORIGINAL CONTOUR POINTS. CCP03010
NPCO-NPC CCP03020
C***» THE CASE OF POSITIVE(NEGATIVE) CONTOUR AREA, DETERMINE, FOR EACH CCP03030
C***CONTOUR POINT, WHETHER THE HEADING OF THE VECTOR FROM THE CONTOUR CCP03040
C***POINT TO THE HILL TOP X,Y POINT LIES BETWEEN THE HEADINGS OF THE CCP03050
C***VECTORS FROM THE HILL TOP X,Y POINT TO THE FIRST(LAST) CONTOUR POINTCCP03060
•C***AND FROM THE HILL TOP X,Y POINT TO THE LAST(FIRST) CONTOUR POINT. CCP03070
C***IF THIS IS SO, THEN ASSIGN AN ADDITIONAL CONTOUR POINT AT THE CCP03080
C**«IERMINATION OF THE VECTOR HAVING A HEADING EQUAL TO THAT OF THE CCP03090
C***VECTOR FROM THE ORIGINAL CONTOUR POINT TO THE HILL TOP X,Y POINT CCP03100
C***AND HAVING A LENGTH EQUAL TO TWICE THE LENGTH OF THIS VECTOR. IF CCP03110
.C***THE ADDITION OF A CONTOUR POINT WOULD CAUSE THE MAXIMUM NUMBER OF CCP03120
C***CONTOUR POINTS TO BE EQUALED, THEN THE COORDINATES OF THIS CONTOUR CCP03130
C***POINT ARE SET EQUAL TO THE COORDINATES OF THE FIRST CONTOUR POINT CCP03140
C***AND THE CONTOUR COMPLETION PROCESS IS HALTED. CCP03150
DO 100 K-1,NPCO ' CCP03160
XP-XCON(K) CCP03170
YP-YCON(K) CCP03180
CALL VECTOR(XP,YP,XREF,YREF,ANGLE,DX,DY) CCP03190
IF(ANG1.GT.ANG2) GO TO 40 CCP03200
IF(ANGLE.GT.ANG1.AND.ANGLE.LT.ANG2) GO TO 10 CCP03210
GO TO 100 CCP03220
10 NPC-NPC+1 CCP03230
IF(NPC.LT.NPCKAX) GO TO 30 CCP03240
XCON(NPC)-XCON(1J CCP03250
YCON(NPC)-YCON(1) CCP03260
WRITE(DOUT,20) NPCMAX CCP03270
20 FORMAT(/,IX,'***WARNING***CONTOUR COMPLETION HALTED DUE TO EXCEEDACCP03280
&NCE OF1,/,IX,'MAXIMUM NUMBERC,13,IX,') OF CONTOUR POINTS1,/, CCP03290
SIX,'THE FINAL CONTOUR POINT WILL HAVE COORDINATES EQUAL TO THOSE OCCP03300
134
-------
4F THE IKITIAL POINT1) CCP03310
RETURN CCP03320
30 XCON(NPC)-XREF+DX . CCP03330
I YCON(NPC)-YREF+DY CCP03340
I GO TO 100 CCP03350
: 40 IF(ANGLE.LT.ANG1.AND.ANGLE.GT.ANG2) GO TO 100 CCP03360
NPONPOH CCP03370
IF(NPC.LT.NPCKAX) GO TO 50 CCP03380
XCON(NPC)-XCON(1) CCP03390
YCON(NPC)-YCON(1) CCP03400
WRITE(DOOT,20) NPCMAX CCP03410
RETURN CCP03420
50 XCON(NPC)-XREF+DX CCP03430
YCON(NPC)-YREF+DY CCP03440
100 CONTINUE CCP03450
C***CLOSE OUT THE CONTOUR BY ADDING A FINAL POINT WITH COORDINATES CCP03460
C***EQUAL TO THOSE OF THE INITIAL POINT. CCP03470
NPC-NPC+1 CCP03480
XCON(NPC)-XCON(1) CCP03490
YCON(NPC)«YCON(1) CCP03500
RETURN CCP03510
END CCP03520
135
-------
i SUBROUTINE ISORT(LIST,NDL,LPTR)
C***MERGE EXCHANGE SORT
£***NUMBER OF COMPARISONS»N*LOG(N)/LOG(2)
C***LIST-ARRAY TO BE SORTED
C***NDL-NUMBER OF ARRAY ELEMENTS TO BE SORTED
C***LPTR-WORKING ARRAY
DIMENSION LIST(1),LPTR(1)
C***CHECK INITIAL ORDER
IF(NDL.LE.l) RETURN
DO 10 1-2, NDL
IP(LIST(I-1).GT.LIST(I)) GO TO 20
10 CONTINUE
RETURN
C***BEGIN SORT
20 L2I-1
DO 100 1-1,20
M-l
L2IH-L2I
, L2I-2*L2I
IF(L2IH.GT.NDL) GO TO 110
JUP-NDVL2I+1
DO 90 J-1,JUP
N-M+L2IH
CT(N.GT.NDL) GO TO 90
KLO-K
KUP-MINO(KLO+L2I-1,NDL)
MUP-KLOH-L2IH-1
DO 60 K-KLO,KUP
IF(M.GT.NDL) GO TO 30
IF(M.GT.KUP) GO TO 30
IT(M.GT.MUP) GO TO 40
IT(LIST(M).GT.LIST(N)) GO TO 40
30 NL-M
M-M+1
GO TO SO
40 ML-N
SO
60
70
80
90
100
110
LPTR(K)-LIST(NL)
CONTINUE
DO SO KHKLO,KUP
LIST(K)-LPTR(K)
CONTINUE
M-KLO+L2I
CONTINUE
CONTINUE
RETURN
END
IS000010
ISO00020
ISO00030
ISO00040
IS000050
ISO00060
ISO00070
IS000080
ISO00090
ISO00100
IS000110
IS000120
ISO00130
IS000140
IS000150
IS000160
ISO00170
ISO00180
IS000190
ISO00200
ISO00210
ISO00220
IS000230
ISO00240
ISO00250
ISO00260
IS000270
ISO00280
IS000290
ISO00300
IS000310
ISO00320
ISO00330
ISO00340
ISO00350
ISO00360
ISO00370
ISO00380
ISO00390
ISO00400
ISO00410
IS000420
ISO00430
IS000440
ISO00450
ISO00460
IS000470
136
-------
SUBROUTINE MULTC(XCON,YCON,NPC,NPCMAX,MCFLAG) MTC00010
C***THIS SUBROUTINE DETERMINES WHETHER A CONTOUR IS REALLY A SERIES MTC00020
C***OF MULTIPLE CONTOURS. IF THIS IS FOUND TO BE THE CASE, THEN THE MTC00030
C***CONTOUR POINT NUMBERING SCHEME IS MODIFIED TO SHOW A SERIES OF MTC00040
C***CONTOURS CONNECTED TO THE FIRST CONTOUR IN THE SERIES BY INFINITELY MTC00050
C***THIN STRIPS FOR THE PURPOSE OF CALCULATING THE AREA, COKTROID MTC00060
C***COORDINATES, AND SECOND MOMENTS OF THE COMPONENT CONTOURS TAKEN MTC00070
C***AS A GROUP. MTC00080
C*** MTC00090
C*** MTC00100
C GLOSSARY OF TERMS MTC00110
C*** MTC00120
C*** MTC00130
C ISO-SIGN(+ OR -) OF THE AREA OF THE FIRST COMPONENT CONTOUR MTC00140
C ISN-SIGN(+ OR -) OF THE AREA OF THE Nth COMPONENT CONTOUR(NOT MTC00150
C INCLUDING THE FIRST COMPONENT CONTOUR) MTC00160
C K-COUNTER FOR THE INPUT SET OF CONTOUR POINTS MTC00170
C XCOUNT-COUNTER FOR THE FINAL SET OF CONTOUR POINTS MTC00180
C KFIN(N)-ENDING VALUE OF KCOUNT FOR THE Nth COMPONENT CONTOUR MTC00190
C (NOT COUNTING THE FIRST COMPONENT CONTOUR) MTC00200
C KSTART(N)-STARTING VALUE OF KCOUNT FOR THE Nth COMPONENT CONTOUR MTC00210
C (NOT COUNTING THE FIRST COMPONENT CONTOUR) MTC00220
C MCFLAG-SUBROUTINE COMPLETION CODE MTC00230
C -0(MAXIMUM NUMBER OF POINTS EXCEEDED IN THE CONTOUR POINT MTC00240
C REASSIGNMENT PROCESS—CONTOUR REJECTED) MTC00250
C -1(THE LAST IN A SERIES OF MULTIPLE CONTOURS WAS FOUND MTC00260
C HOT TO BE CLOSED—CONTOUR REJECTED) MTC00270
C -2(CONTOUR HAS FOUND TO BE A SINGLE CONTOUR(I.E. NO CONTOUR MTC00280
C CLOSURE WAS FOUND BEFORE THE FINAL CONTOUR POINT)) MTC00290
C -3(POINT REASSIGNMENT FOR THE MULTIPLE CONTOUR WAS MTC00300
C SUCCESSFULLY COMPLETED) MTC00310
C -4(ALL COMPONENT CONTOURS NOT INPUT WITH POINTS IN THE SAME MTC00320
;C ORDER. THE ORDER OF POINT INPUT FOR THE COMPONENT CONTOURS MTC00330
C HAS BEEN MADE THE SAME AS THE FIRST COMPONENT CONTOUR. MTC00340
C FOLLOWING THIS ACTION, THE POINT REASSIGNMENT FOR THE MTC00350
C MULTIPLE CONTOUR WAS SUCCESSFULLY COMPLETED.) MTC00360
C NCON-NUMBER OF COMPONENT CONTOURS NOT INCLUDING THE FIRST MTC00370
C (INCREMENTED DURING THE COURSE OF THE ANALYSIS) MTC00380
C HN-POINT COUNTERd TO (KFIN(N)-KSTART(N))+l) WITHIN COMPONENT MTC00390
C CONTOUR X MTC00400
C ICON-ARRAY CONTAINING X COORDINATES OF THE INITIAL AND FINAL SET MTC00410
C OF CONTOUR POINTS MTC00420
C YCON-ARRAY CONTAINING Y COORDINATES OF THE INITIAL AND FINAL SET MTC00430
C OF CONTOUR POINTS MTC00440
C XCONS,YCONS-WORKING ARRAYS FOR CONTOUR POINT REASSIGNMENT MTC00450
C*** . MTC00460
C*** MTC00470
DIMENSION XCON(1000),YCON(1000),XCONS(1000),YCONS(1000) MTC00480
DIMENSION KSTART(500),KFIN(500) MTC00490
NCON-1 MTC00500
XCONS(1)-XCON(1) MTC00510
YCONS(1)-YCON(1) • MTC00520
DO 100 K-2.NPC MTC00530
KSAVE-K MTC00540
XCONS(K)-XCON(K) . MTC00550
YCONS(X)-YCON(K) MTC00560
C***DETERMINE WHETHER THE CONTOUR CLOSES BEFORE THE LAST CONTOUR POINT MTC00570
C***HAS BEEN REACHED. IF SO, ASSUME THE CONTOUR IS COMPOSED OF MULTIPLE MTC00580
C***CONTOURS AT THE SAME ELEVATION. CONTINUE WITH THE ANALYSIS. IF NOT, MTC00590
C***THEN RETURN TO THE MAIN PROGRAM WITH A COMPLETION CODE OF 2. MTC00600
XF(ABS(XCON(K)-XCON(1)).LT.1.0E-15.AND.ABS(YCON(K)-YCON(1)).LT. MTC00610
41.0E-1S.AND.X.NE.NPC) GO TO 110 MTC00620
100 CONTINUE ' MTC00630
MCFLAG-2 MTC00640
GO TO 400 MTC00650
110 CONTINUE MTC00660
137
-------
4***DETERMINE THE AREA OF THE FIRST COMPONENT CONTOUR AND ITS SIGN FOR MTC00670
C***LATER USE. MTC00680
CALL ARCM(XCON,YCON,AREA,XCM,YCM,KSAVE) MTC00690
ISO— 1 MTC00700
IF(AREA.LT.O.) ISO— 1 MTC00710
KSP1-KSAVE+1 MTC00720
KSTART(1)-KSP1 MTC00730
C***STORE THE COORDINATES OF THE FIRST POINT OF THE SECOND COMPONENT MTC00740
C***CONTOUR IN THE TEMPORARY STORAGE ARRAYS. MTC00750
i XCONS(KSP1)-XCON(KSP1) MTC00760
: YCONS(KSP1)-YCON(KSP1) MTC00770
I KSP2-KSAVE+2 _ MTC00780
C***IF ONLY ONE ADDITIONAL POINT HAS BEEN SPECIFIED AFTER THE FIRST KTC00790
<***CONTOUR CLOSURE, THEN RETURN TO THE MAIN PROGRAM WITH A COMPLETION MTC00800
C***CODE OF 1. MTC00810
IF(KSP2.LE.NPC) GO TO 150 MTC00820
MCFLAG-1 MTC00830
GO TO 400 MTC00840
i 150 CONTINUE MTC00850
C***SPECIFY THE BEGINNING POINT OF THE SECOND COMPONENT CONTOUR AS KTC00860
C***(XCOMP,YCOMP). MTC00870
i XCOMP-XCON(KSPl) MTC00880
; YCOMP-YCON(KSPl) MTC00890
KCOUNT-KSP2 . MTC00900
i K-KCOUNT MTC00910
C***UP TO THIS POINT THE NUMBER OF THE INPUT AND MODIFIED CONTOUR POINTSMTC00920
C***IS STILL THE SAME. NOW ENTER THE LOOP WHICH CARRIES OUT THE POINT MTC00930
C***REASSIGNMENT PROCESS. MTC00940
200 CONTINUE MTC00950
i XCONS (KCOUNT) -XCON (K) MTC00960
i YCONS (KCOUNT) -YCON(K) MTC00970
C***HAS THE NEXT CLOSURE BEEN REACHED? IF SO, RETURN TO THE POINT MTC00980
C***OF FIRST CLOSURE(XCON(1),YCON(1)) BEFORE CONTINUING. MTC00990
IF(ABS(XCON(K)-XCOMP) .GT.1.0E-15.OR.ABS (YCON(K) -YCOMP) .GT.1.0£-15)MTC01000
tGO TO 210. MTC01010
C***SPECIFY THE END POINT FOR COMPONENT CONTOUR NCON. MTC01020
KFIN(NCON) -KCOUNT MTC01030
C***nrCREMENT COUNTER FOR THE SET OF MODIFIED CONTOUR POINTS. MTC01040
: KCOUNT-KCOUNT-H . MTC01050
C***CRECK WHETHER THE MAXIMUM NUMBER OF CONTOUR POINTS HAS BEEN EXCEEDEDMTC01060
IF(KCOUNT.LE.NPCMAX) GO TO 205
I MCFLAG-O
GO TO 400
205 CONTINUE
C***RETURN TO CLOSURE POINT FOR FIRST COMPONENT CONTOUR.
XCONS(KCOUNT)-XCON(1)
YCONS(KCOUNT)-YCON(l)
C***IHCREMENT COUNTER FOR THE ORIGINAL SET OF POINTS.
K-K+1
-------
C***CLOSURE. MTC01330
XCOMP-XCON(K) MTC01340
YCOMP-YCON(K) MTC01350
C***INCREM£NT THE COUNTER FOR THE INPUT SET OP CONTOUR POINTS. MTC01360
210 K-K+1 MTC01370
C***DETERMINE WHETHER THE NUMBER OF INPUT CONTOUR POINTS HAS BEEN MTC01380
C***EXKAUSTED. MTC01390
IF(K.LE.NPC) GO TO 250 MTC01400
MCFLAG-1 MTC01410
GO TO 400 MTC01420
250 CONTINUE MTC01430
C***INCREMENT THE COUNTER FOR THE MODIFIED SET OF CONTOUR POINTS. MTC01440
KCOUNT-KCOUNT+1 MTC01450
C***CHECK WHETHER THE NUMBER OF CONTOUR POINTS HAS BEEN EXCEEDED. MTC01460
IF(KCOUNT.LE.NPCMAX) GO TO 200 MTC01470
MCFLAG-0 MTC01480
GO TO 400 MTC01490
300 CONTINUE MTC01500
NPOKCOUNT MTC01510
C***TRANSFER THE POINT COORDINATES FROM THE TEMPORARY HOLDING ARRAYS MTC01520
C***TO THE INITIAL POINT COORDINATE ARRAYS. . MTC01530
DO 350 K-1,KCOUNT MTC01540
XCON(X)-XCONS(K) MTC01550
YCON(K)-YCONS(K) MTC01560
350 CONTINUE MTC01570
KCFLAG-3 MTC01580
C***DETERMINE WHETHER ALL COMPONENT CONTOURS HAVE THEIR POINTS INPUT MTC01590
C***IN THE SAME SENSE(CLOCKWISE OR COUNTER-CLOCKWISE). IF NOT, MODIFY MTC01600
C***THB INPUT SEQUENCES TO REFLECT THE SEQUENCE USED FOR THE FIRST MTC01610
C***COMPONENT CONTOUR. MTC01620
DO 390 N-1,NCON MTC01630
NN-0 MTC01640
DO 370 K-KSTART(N),KFIN(N),1 . MTC01650
XN-NN+1 MTC01660
ICONS(NN)-XCON(K) MTC01670
YCONS(NN)»YCON(K) MTC01680
370 CONTINUE . MTC01690
C***FIND THE AREA OF THE COMPONENT CONTOUR AND ITS SIGN. IF THE SIGN MTC01700
C***OF THE AREA Iff DIFFERENT FROM THE SIGN OF THE AREA OF THE INITIAL MTC01710
C***COMPONENT CONTOUR, THEN REVERSE THE INPUT ORDER OF THE COMPONENT MTC01720
C***CONTOUR POINTS. MTC01730
CALL ARCM(XCOKS,YCONS,AREA,XCM,YCM,NN) MTC01740
ISN-1 • MTC01750
IF(AREA.LT.O.) ISN—1 MTC01760
IF(ISN.BQ.ISO) GO TO 390 MTC01770
MCFLAG-4 . MTC01780
DO 380 K-KSTART(N),KFIN(N),1 MTC01790
XCON(K)-XCONS(NN) MTC01800
YCON(K)-YCONS(NN) MTC01810
NN-NN-1 MTC01820
380 CONTINUE MTC01830
390 CONTINUE • MTC01840
400 CONTINUE MTC01850
RETURN MTC01860
END MTC01870
139
-------
SUBROUTINE SKIPCN(IN,NPC) SKP00010
C***SUBROUTINE TO SKIP CONTOUR POINTS FOR CONTOURS WHICH ARE NOT SKP00020
C***TO BE PROCESSED SKP00030
READ(IN,*) (XDUM,YDUM,K-1,NPC) SKP00040
RETURN SKPOOOSO
END SKP00060
140
-------
SUBROUTINE SMOHNT(XCON,YCON,AR,NSLOPE,SN,CN,ANGLE,NPC, SMO00010
«XC,YC,RG,RGRAT,ORENT,ISMFLG) SMO00020
C***SUBROUTINE TO CALCULATE THE SECOND MOMENTS AND RADII OF GYRATION SMO00030
C***FOR THE INPUT CONTOUR ABOUT AXES OF EQUAL ANGULAR SPACING AND SMO00040
C***WHICH PASS THROUGH 'THE CENTROID OF THE CONTOUR IN THE PLANE OF THE SMOOOOSO
C***CONTOUR . • SMO00060
C***INITIALIZE VALUES FOR THE MAXIMUM AND MINIMUM RADII OF GYRATION SMO00070
DIMENSION SN(1),CN(1),XCON(1),YCON(1),ANGLE(1) SMOOOOSO
ISMFLG-0 SMO00090
RGMAX-0. SM000100
RGMIN-l.OE+15 SMO00110
C***BEGIN LOOP OVBR AXIS ORIENTATION VALUES SM000120
DO 200 M-1,NSLOPE SMO00130
C***INITIALIZE THE SECOND MOMENT FOR THIS AXIS TO ZERO SMO00140
SMOK-0. SMO001SO
C***BEGIN LOOP OVER CONTOUR POINTS SMO00160
C***D1»PERPENDICULAR DISTANCE TO THE AXIS LINE FROM CONTOUR POINT K SMO00170
C***D2-PERPENDICULAR DISTANCE TO THE AXIS FROM CONTOUR POINT K+l SMO00180
C***W*DISTANCE ALONG THE AXIS LINE BETWEEN THE INTERSECTION OF SMO00190
C*** PERPENDICULARS FROM ADJACENT CONTOUR POINTS(K AND K+l) SMO00200
NPCMl-NPC-1 SM000210
Dl—(XCON(1)-XC)*SN(M) + (YCON(1)-YC)*CN(M) SMO00220
DO 100 K-1,NPCM1 - SMO00230
D2— (XCON(K+1)-XC)*SN(M)+(YCON(K+1)-YC)*CN(M) SMO00240
W-(XCON(K+1)-XCON(K) ) *CN(M)+(YCON(K+1)-YCON(K)) *SN(M) SMO00250
SMOM«SMOM+(W/12.)*(D2**3+D1**3+D1*D2**2+D2*D1**2) SMO00260
01-02 SM000270
100 CONTINUE SM000280
IF(SMOH/AR.LT.O.O) ISMFLG-1 SMO00290
IF(ISMFLG.EQ.l) RETURN SMO00300
C***CALCULATE THE RADIUS OF GYRATION OF THE CONTOUR ABOUT THIS AXIS SMO00310
RG-SQRT(SMOM/AR) SMO00320
C***UPDATE THE VALUES FOR THE MAXIMUM AND MINIMUM RADII OF GYRATION SMO00330
C***AND SAVE THE ORIENTATION INDEX FOR THE AXIS HAVING THE CURRENT SMO00340
C***LARGEST RADIUS OF GYRATION . SMO00350
IF(RG.LT.RGMAX) GO TO 150 SMO00360
' RGMAX-RG . . SM000370
MMAX-M . SMO00380
150 IF(RG.GT.RGMIN). GO TO 200 SMO00390
RGMIN-RG SMO00400
200 CONTINUE . SMO00410
RG-RGMAX SMO00420
RGRAT»TRGMAX-RGMIN)/RGMAX ' SMO00430
ORENT-ANGLE(MHAX) SMO00440
RETURN SM000450
END SMO00460
141
-------
SUBROUTINE VECTOR(XBEG,YBEG,XEND,YEND,ANGLE,DX,DY) . VEC00010
C***SUBROUTINE CALCULATES THE DIRECTION(DEGREES) AND X,Y COMPONENTS VEC00020
C***FOR A VECTOR FROM (XBEG.YBEG) TO (XEND,YEND). THE COMPUTED VEC00030
C***DIRECTIONS RANGE FROM 0 TO 360 DEGREES VEC00040
PI-3.14159265 VEC00050
DX-XEND-XBEG VEC00060
DY-YEND-YBEG VEC00070
I?(ABS(DX).LT.1.0E-15.AND.ABS(DY).LT.l.OE-15) GO TO 10 VEC00080
ANGLE-(180./PI)*ATAN2(DY,DX) VEC00090
GO TO 20 VEC00100
10 WGLE-0. VEC00110
20 I?(ANGLE.LT.O.) ANGLE-360.+ANGLE VEC00120
RETURN VEC00130
END VEC00140
142
-------
HCHIT MAIN PROGRAM AND SUBROUTINE PSORTR
143
-------
PROGRAM HCRIT HCT00010
C***PROGRAM TO FIT ELLIPTICAL CONTOURS TO AN INVERSE POLYNOMIAL HILL HCT0002a
C***SHAPE FOR A RANGE OF USER SPECIFIED CRITICAL CUTOFF ELEVATIONS. HCT00030
C***THE PROGRAM PROVIDES CRITICAL ELEVATION HILL PARAMETERS FOR HCT00040
C***INPUT TO THE COMPLEX TERRAIN DISPERSION MODEL(CTDM). HCT00050
C*** HCT00060
C*** HCT00070
C GLOSSARY OF TERMS HCT00080
HCT00090
HCT00100
A(J)-SEMI-MAJOR AXIS LENGTH FOR CONTOUR J(USER COORDINATES) HCT00110
AI-INTERPOLATED VALUE OF A(J) TO A GIVEN CRITICAL ELEVATION HCT00120
AS(J)-TEMPORARY A(J) STORAGE ARRAY USED IN SORTING HCT00130
ANS-CHARACTER*! VARIABLE HOLDING THE ANSWER TO A YES(Y) OR NO(N) HCT00140
QUESTION HCT00150
B(J)-SEMI-MINOR AXIS LENGTH FOR CONTOUR J(USER COORDINATES) HCT00160
BI-INTERPOLATED VALUE OF B(J) TO A GIVEN CRITICAL ELEVATION HCT00170
BS(J)-TEMPORARY B(J) STORAGE ARRAY USED IN SORTING HCT00180
ECC(J)-ECCENTRICITY OF CONTOUR J HCT00190
ECCS(J)-TEMPORARY ECC(J) STORAGE ARRAY USED IN SORTING HCT00200
FCONFILE-CHARACTER*15 VARIABLE CONTAINING THE INPUT FILE NAME FOR HCT00210
THE FITTED CONTOUR PARAMETERS GENERATED BY PROGRAM FITCON HCT00220
FEXT-EXTRAPOLATION FACTOR USED TO ASSIGN THE SEMI-MAJOR AND HCT00230
SEMI-MINOR AXIS LENGTHS FOR THE CASE OF ONE CONTOUR AND HCT00240
A CRITICAL ELEVATION BELOW THAT SINGLE CONTOUR HCT00250
FRACT-FRACnONAL DIFFERENCE OF THE CRITICAL ELEVATION BETWEEN HCT00260
ADJACENT CONTOUR ELEVATIONS HCT00270
HC(I)-ARRAY OF CRITICAL ELEVATIONS HCT00280
HCLOW-THE LOWEST CRITICAL ELEVATION(INPUT FOR CRITICAL ELEVATION HCT00290
SELECTION MODE 2) HCT00300
HCON(J)-ELEVATION OF CONTOUR J(USER COORDINATES) HCT00310
HCONM1-HCX)N(NC)-1. HCT00320
HCONS(J)-TEMPORARY HCON(J) STORAGE ARRAY USED IN SORTING HCT00330
HHILL-HEIGHT OF THE RILL TOP ABOVE A GIVEN CRITICAL ELEVATION HCT00340
HNAME-CHARACTER*15 VARIABLE GIVING THE HILL NAME HCT00350
HTOP-HILL TOP ELEVATION(USER COORDINATES) HCT00360
ICHMOD-CRmCAL ELEVATION INPUT MODE FOR THE RILL IN QUESTION HCT00370
-1(CRITICAL ELEVATIONS WILL BE AT CONTOUR ELEVATIONS WITH HCT00380
THE' EXCEPTION OF THE UPPERMOST CONTOUR) HCT00390
-2(CRITICAL ELEVATIONS EVENLY SPACED BETWEEN A USER SUPPLIED HCT00400
C***
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
• C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
LOWER ELEVATION AND THE ELEVATION OF THE UPPERMOST CONTOUR)HCT00410
IDC(J)-ID NUMBER FOR CONTOUR J
IDHILL-HILL ID NUMBER(1-999)
IN-UNIT NUMBER FOR THE FITTED HILL INPUT FILE FROM PROGRAM FITCON
LA-LENGTH PARAMETER FOR AN INVERSE POLYNOMIAL FIT ALONG THE HILL
MAJOR AXIS(USER COORDINATES)
LB-LENGTH PARAMETER FOR AN INVERSE POLYNOMIAL FIT ALONG THE HILL
MINOR AXIS(USER COORDINATES)
LPTR-WORKING ARRAY USED IN THE POINTER SORT(PSORTR)
MOOT-UNIT NUMBER FOR THE FILE (MOUTFILE) CONTAINING TERRAIN
PARAMETERS WHICH ARE PASSED TO CTDM
MOUTFILE-CHARACTER*15 VARIABLE CONTAINING THE OUTPUT FILE NAME FOR
THE PARAMETERS TO BE PASSED TO CTDM
HC-NUMBER OF FITTED CONTOURS INPUT FROM FCONFILE
NCHMAX-MAXIMUH NUMBER OF CRITICAL ELEVATIONS WHICH CAN BE ANALYZED
HCON-NUMBER OF CONTOURS USED IN FITTING A RILL FOR A GIVEN
ELEVATION
NCR-NUMBER OF CRITICAL ELEVATIONS USED
NPTR-ARRAY CONTAINING THE SORTED POINTERS RETURNED FROM SUBROUTINE
PSORTR
ONOR-MAJOR AXIS ORIENTATION IN DEGREES CLOCKWISE FROM NORTH(FOR A
CONTOUR OR A FITTED RILL)(BETWEEN 0 AND 180 DEGREES)
OREN(J)-ORIENTATION ANGLE OF THE CONTOUR J SEMI-MINOR AXIS WITH
RESPECT TO THE POSITIVE X-AXIS
ORENF-ORIENTATON OF THE MINOR AXIS OF A FITTED HILL AS MEASURED
COUNTER CLOCKWISE FROM THE POSITIVE X-AXIS(EAST)
HCT00420
HCT00430
HCT00440
HCT00450
HCT00460
HCT00470
HCT00480
HCT00490
HCT00500
HCT00510
HCT00520
HCT00530
HCT00540
HCT00550
CRITICAL HCT00560
HCT00570
HCT00580
HCT00590
HCT00600
HCT00610
HCT00620
HCT00630
HCT00640
HCT00650
HCT00660
144
-------
c
c
c
c
c
c
c
c
c
c
•c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c***
c***
ORENI-INTERPOLATED VALUE OF OREN(J) TO A GIVEN CRITICAL ELEVATION HCT00670
ORENS(J)-TEMPORARY OREN(J) STORAGE ARRAY USED IN SORTING HCT00680
PA-EXPONENT FOR AN INVERSE POLYNOMIAL FIT ALONG THE HILL MAJOR AXIS HCT00690
PB-EXPONENT FOR AN INVERSE POLYNOMIAL FIT ALONG THE HILL MINOR AXIS HCT00700
PFILE-CHARACTER*15 VARIABLE GIVING THE NAME OF THE CRITICAL HCT00710
ELEVATION PLOT FILE. THIS NAME MOST BE DIFFERENT THAN THE NAMEHCT00720
OF THE PLOT FILE GENERATED BY PROGRAM FITCON. BOTH PLOT FILES HCT00730
ARE EVENTUALLY BE INPUT TO THE PLOT PROGRAM.
PFLAG-PLOT GENERATION INDICATOR
-0(NO PLOT GENERATED)
-1(PLOT GENERATED)
PSORTR-SUBROUTINE FOR SORTING POINTERS(CALLED TO SORT CONTOUR
FIT PARAMETERS BY CONTOUR ELEVATION(ASCENDING ORDER))
SUM1,SUM2A,SUM2B,SUM3,SUM4A,SUM4B-SUMMATION VARIABLES USED IN THE
CALCULATION OF BEST FIT INVERSE POLYNOMIAL HILL PROFILES
SUMX,SUMY-INTERMEDIATE VARIABLES USED IN THE DETERMINATION OF THE
ORIENTATIONS OF INTERPOLATED CONTOURS AND FITTED HILLS
UPL-UNIT NUMBER FOR THE CRITICAL ELEVATION PLOT FILE
XCM(J)-X-COORDINATE OF THE CONTOUR J CENTROID(USER COORDINATES)
XCMI-INTERPOLATED VALUE OF XCM(J) TO A GIVEN CRITICAL ELEVATION
XCMS(J)-TEMPORARY XCM(J) STORAGE ARRAY USED IN SORTING
XHTOPF-AVERAGE OF THE X-COORDINATES OF CONTOUR CENTROIDS ABOVE
A GIVEN CRITICAL ELEVATION
YCM(J)-Y-COORDINATE OF THE CONTOUR J CENTROID(USER COORDINATES)
YCMI-INTERPOLATED VALUE OF YCM(J) TO A GIVEN CRITICAL ELEVATION
YCMS(J)-TEMPORARY YCM(J) STORAGE ARRAY USED IN SORTING
YHTOPF-AVERAGE OF THE Y-COORDINATES 'OF CONTOUR CENTROIDS ABOVE
A GIVEN CRITICAL ELEVATION
DIMENSION A(200),AS(200),B(200),BS(200),ECC(200),ECCS(200),
*HCON(200),HCONS(200),IDC(200),OREN(200),ORENS(200),XCM(200),
4XCMS(200),YCM(200),YCHS(200),LPTR(200),NPTR(200),HC(200)
REAL*4LA,LB
INTEGER UPL
CHARACTER*! ANS
CHARACTER*1S FCONFILE,MOUTFILE,PFILE,HNAME
C***
c***
C INITIALIZATION OF VARIABLES
C***
c***
C***SPECIFY FILE UNIT NUMBERS.
IN-14
MOOT-IS
UPL-16
C***SPECIFY CONSTANTS.
PI-3.14159263
NCHMAX-200
C***
C***
C INPUT FILE NAMES.
C***
C***
C***ENTER THE NAME OF THE INPUT FILE CONTAINING THE CONTOUR FIT
C***PARAMETERS GENERATED BY PROGRAM FITCON.
5 WRITE(*,10)
10 FORMAT(/,IX,'ENTER INPUT FILE NAME(FROM FITCON) -> *\)
READ(*,'(A)•) FCONFILE
IF(FCONFILE.EQ.« •) GO TO S
C***OPEN THE INPUT FILE.
. OPEN(IN,FILE-FCONFILE,STATUS-IOLDI)
C***ENTER THE NAME OF THE OUTPUT FILE WHICH WILL BE PASSED DIRECTLY
C***TO CTDM.
15 WRITE(*,20)
20 FORMAT(/,IX,'ENTER OUTPUT FILE NAME(FOR CTDM) ->'\)
HCT00740
HCT00750
HCT00760
HCT00770
HCT00780
HCT00790
HCT00800
HCT00810
HCT00820
HCT00830
HCT00840
HCT00850
HCT00860
HCT00870
HCT00880
HCT00890
HCT00900
HCT00910
HCT00920
HCT00930
HCT00940
HCT00950
HCT00960
HCT00970
HCT00980
HCT00990
HCT01000
HCT01010
HCT01020
HCT01030
HCT01040
HCT010SO
HCT01060
HCT01070
HCT01080
HCT01090
HCT01100
HCT01110
HCT01120
HCT01130
HCT01140
HCT011SO
HCT01160
HCT01170
HCT01180
HCT01190
HCT01200
HCT01210
HCT01220
HCT01230
HCT01240
HCT01250
HCT01260
HCT01270
HCT01280
HCT01290
HCT01300
HCT01310
HCT01320
145
-------
READ(*,'(A) •) MOOTFILE
IF(MOUTFILE.EQ.« •) GO TO 15
•C***OPEN THE OUTPUT FILE.
; OPEN (MOUT,FILE-MOUTFILE, STATUS- 'NEW')
C***
C***
C DETERMINE WHETHER A PLOT IS TO BE GENERATED.
C*** •
C***FIRST, INITIALIZE THE
-C***"NO" ANSWER.
PFLAG-0
WRITE(*,30)
PLOT FLAG INDICATOR TO CORRESPOND TO A
30 FORMAT (/, IX, 'PLOT REQUESTED? (Y/N) -> «\)
'
PFLAG-1
READ(*,(A)') ANS
IF(ANS.EQ. 'Y'.OR.ANS.EQ. 'y')
IF(PFLAG.EQ.O) GO TO 50
C***INPUT THE NAME OF THE PLOT FILE.
35 WRITE(*,40)
40 FORMAT (/, IX, 'ENTER PLOT FILE NAME -> «\)
READ(*,'(A)') PFILE
; IF ( PFILE. EQ. ' •) GO TO 35
C***OPEN THE PLOT FILE.
i OPEN (OTL.FILE-PFILE, STATUS- 'NEW')
C***WRITE "HCRIT" TO THE FIRST RECORD OF THIS PLOT FILE TO INDICATE
C***THAT THIS PLOT FILE IS GENERATED BY PROGRAM HCRIT.
WRITE (UPL, 45)
45 FORMAT (' HCRIT1)
; 50 CONTINUE
*C***
;C***
C READ DATA FROM THE FITTED CONTOUR FILE.
C***
i C***
I C***INPUT THE HILL ID NUMBER AND NAME.
I READ(IN,60) IDHILL,HNAME
i 60 FORMAT(I2,1X,A15) '
j C***INPUT THE HILL TOP ELEVATION. .
I RZAD(IN,70) HTOP . . .
; 70 FORMAT(E15.4)
jC***INPUT THE NUMBER OF FITTED CONTOURS.
: READ(IN,80) NC .
: 80 FORMAT (110)
! C***INPUT THE SORTED CONTOUR IDs. THESE IDs ARE WRITTEN TO THE
i C***PLOT FILE AND COMPARED WITH THE SORTED IDs WRITTEN TO THE
i C***PLOT FILE WRITTEN BY FITCON. THIS CHECK PREVENTS THE COMPARISON
C***OF AN ACTUAL AND A FITTED CONTOUR WHICH IN FACT REPRESENT DIFFERENT
C***CONTOURS.
I READ(IN,80) (IDC(J) ,J-1,NC)
; IF(PFLAG.EQ.O) GO TO 85
I C***WRITE TO THE PLOT FILE THE RILL ID NUMBER, KILL NAME, NUMBER
! C***OF FITTED CONTOURS, THE SORTED CONTOUR IDs, AND THE HILL TOP
| C***ELEVATION.
! WRITE(UPL,60) IDHILL,HNAME .
i WRITE (UPL, 80) NC .
! WRITE(UPL,80) (IDC(J) ,J»1,NC)
WRITE (UPL, 70) HTOP .
i 85 CONTINUE •
: C***INPUT THE CONTOUR FIT PARAMETERS FOR THE HILL IN QUESTION. .
I DO 100 J-1,NC
I READ (IN, 90) HCON(J),XCM(J),YCM(J),A(J),B(J),ECC(J),OREN(J)
! 90 FORMAT (7E15. 4)
! 100 CONTINUE
|C***CLOSE THE FITTED CONTOUR INPUT FILE.
CLOSE (IN, STATUS- 'KEEP')
HCT01330
HCT01340
HCT01350
HCT01360
HCT01370
HCT01380
HCT01390
HCT01400
HCT01410
HCT01420
HCT01430
HCT01440
HCT01450
HCT01460
HCT01470
HCT01480
HCT01490
HCT01500
HCT01510
HCT01520
HCT01530
HCT01540
HCT01550
HCT01560
HCT01570
HCT01580
HCT01590
HCT01600
HCT01610
HCT01620
HCT01630
HCT01640
HCT01650
HCT01660
HCT01670
HCT01680
HCT01690
HCT01700
HCT01710
HCT01720
HCT01730
HCT01740
HCT01750
HCT01760
HCT01770
HCT01780
HCT01790
HCT01800
HCT01810
HCT01820
HCT01830
HCT01840
HCT01850
HCT01860
HCT01870
HCT01880
HCT01890
HCT01900
HCT01910
HCT01920
HCT01930
HCT01940
HCT01950
HCT01960
HCT01970
HCT01980
146
-------
C*** HCT01990
C SORT THE CONTOUR PARAMETERS BY CONTOUR ELEVATION (IN ASCENDING ORDER) HCT02000
C BY USE OF A POINTER SORT.
C***
CALL PSORTR(HCON,NC,NPTR,LPTR)
C***REORDER THE CONTOUR FIT PARAMETERS BASED UPON THE RESULTS OF THE
C***POINTER SORT.
DO 110 J-1,NC
HCONS(J)-HCON(NPTR(J))
AS(J)-A(NPTR(J))
BS(J)-B(NPTR(J))
ECCS(J)-ECC(NPTR(J))
ORENS(J)-OREN(NPTR(J))
XCMS(J)-XCM(NPTR(J))
YCMS(J)-YCM(NPTR(J))
110 CONTINUE
DO 120 J-1,NC
i HCON(J)-HCONS(J)
A(J)-AS(J)
B(J)-BS(J)
ECC(J)-ECCS(J)
OREN(J)-ORENS(J)
XCM(J)-XCMS(J)
YCM (J)-YCMS (J)
120 CONTINUE
i IF(PFLAG.EQ.O) GO TO 123
iC***WRITE THE SORTED CONTOUR ELEVATIONS TO THE PLOT FILE.
DO 122 J-1,NC
WRITE(UPL,121) HCON(J)
121 FORMAT(E15.4)
122 CONTINUE
123 CONTINUE
iC*«*
! C***
C DETERMINE CRITICAL ELEVATIONS TO BE USED FOR FITTING CUTOFF
C HILLS.
C***
i C***
iC***TWO MODES ARE AVAILABLE FOR THE INPUT OF CRITICAL ELEVATIONS.
C***THE USER MAY SPECIFY THAT EACH CONTOUR LEVEL(WITH THE EXCEPTION
C***OF THE UPPERMOST CONTOUR) IS TO BE SPECIFIED AS A CRITICAL
C***ELEVATION, OR THE USER MAY ASK FOR UP TO A MAXIMUM OF NHCMAX
,C***CRTTICAL ELEVATIONS EVENLY SPACED BETWEEN A USER SPECIFIED LOWER
C***BLEVATION AND THE UPPER MOST CONTOUR OF THE HILL.
; 125 WRITE(*,130)
130 FORMAT(//,22X,'SPECIFY CRITICAL HEIGHT SELECTION MODE',/,
HCT02010
HCT02020
HCT02030
HCT02040
HCT02050
HCT02060
HCT02070
HCT02080
HCT02090
HCT02100
HCT02110
HCT02120
HCT02130
HCT02140
HCT02150
HCT02160
HCT02170
HCT02180
HCT02190
HCT02200
HCT02210
HCT02220
HCT02230
HCT02240
HCT022SO
HCT02260
HCT02270
HCT02280
HCT02290
HCT02300
HCT02310
HCT02320
HCT02330
HCT02340
HCT02350
HCT02360
HCT02370
HCT02380
HCT02390
HCT02400
HCT02410
HCT02420
HCT02430
HCT02440
HCT02450
HCT02460
HCT02470
HCT02480
HCT02490
HCT02500
HCT02510
HCT02520
HCT02530
C22X,'!.) AT ALL CONTOUR ELEVATIONS EXCEPT UPPERMOST',/,
C22X,*2.) EVENLY SPACED BETWEEN A USER SUPPLIED ELEVATION',/
S26X,'AND THE UPPERMOST CONTOUR ELEVATION',/,
S26X,'CHOICE?(1 OR 2) -> *\)
READ(*,'(BN,I1)«,ERR-125) ICHMOD
IF(ICHMOD.EQ.l) GO TO 150
IF(ICHMOD.EQ.2) GO TO 170
WRTTE(*,140)
140 FORMAT(/,IX,»***ERROR*** SELECTION MODE OUT OF RANGE—TRY AGAIN')HCT02S40
I GO TO 125 . HCT02550
|C***CRITICAL ELEVATION SELECTION MODE 1 HCT02560
|C***THE NUMBER OF CRITICAL ELEVATIONS WILL BE ONE LESS THAN THE NUMBER HCT02S70
iC***OF CONTOURS. HCT02580
i ISO NCR-NC-1 HCT02S90
IF(NCR.CT.O) GO TO 155 HCT02600
i WRTTE(*,1S1) HCT02610
151 FORMAT(/,IX,'SINCE THERE IS ONLY ONE CONTOUR, THE CONTOUR SELECTIOHCT02620
i CN MODE 1 CANNOT BE USED.',/,IX,'MODE 2 WILL BE USED INSTEAD.')- HCT02630
:C***RESET THE CRITICAL ELEVATION SELECTION MODE. HCT02640
147
-------
ICHMOD-2
GO TO 170
155 DO 160 I-1,HCR
HC(I)-HCON(I)
160 CONTINUE
GO TO 250
C***CRITICAL ELEVATION SELECTION MODE 2
C***READ IN NUMBER OF CRITICAL ELEVATIONS.
170 WRITE(*,180)
HCT02650
HCT02660
HCT02670
HCT02680
HCT02690
HCT02700
HCT02710
HCT02720
HCT02730
180 FORMAT(/,1X, • INPUT THE NUMBER OF CRITICAL ELEVATIONS (1-200) -> '\)HCT02740
. READ(*,'(BN,I3) «,ERR-170) NCR HCT02750
k' IF(NCR.LE.NCHMAX.AND.NCR.NE.O) GO TO 200 HCT02760
WRITE <*, 190) HCT02770
190 FORMAT (/, IX, '***ERROR*** NUMBER OF CRITICAL ELEVATIONS OUT OF RANHCT02780
*E—TRY AGAIN') HCT02790
GO TO 170 HCT02800
C***INPUT THE LOWEST CRITICAL ELEVATION. HCT02810
200 WRITE(*,210) HCT02820
210 FORMAT (/, IX, 'INPUT THE LOWEST CRITICAL ELEVATION -> *\) HCT02830
215 R£AD<*,«(BN,F10.0)) «,ERR-200) HCLOW HCT02840
C***CHECK WHETHER THE LOWEST CRITICAL ELEVATION IS OVER 1 ELEVATION UNITHCT02850
C***ABOVE THE HIGHEST CONTOUR ELEVATION. IF SO, ASK THE USER TO INPUT HCT02860
C***ANOTHER VALUE FOR THE LOWEST CRITICAL ELEVATON.
: HCONMl-HCON(NC)-!. •
IF(HCLOW.LT.HCONMl) GO TO 230
WRITE (*, 220) HCOHM1
220 FORMAT (/, IX, 'LOWEST CRITICAL ELEVATION MUST BE LESS THAN',E15.4
;' »,/, IX, 'TRY AGAIN') '
: GO TO 215
C***ASSIGN THE CRITICAL ELEVATIONS.
C***HCLOW WILL BE THE FIRST ELEVATION. THERE WILL BE NCR-1 ADDITIONAL
C***CRITICAL ELEVATIONS ABOVE HCLOW HAVING A SPACING EQUAL TO DELC,
C***WHERE DELO(HCON(NC)-HLOW)/NCR. THE HIGHEST CRITICAL ELEVATION
C***WILL BE A DISTANCE OF DELC BELOW THE UPPERMOST CONTOUR LEVEL.
230 DELC- (HCON(NC) -HCLOW) /FLOAT (NCR) .
DO 240 I-1,NCR
HC(I)-HCLOW+(I-1)*DELC
240 CONTINUE
i 250 CONTINUE
IF(PFLAG.EQ.O) GO TO 251
C***WRITB THE NUMBER OF CRITICAL ELEVATIONS TO THE PLOT FILE.
1 WRITE (UPL, 80) NCR
251 CONTINUE
t***ASSIGNMENT OF CRITICAL ELEVATIONS COMPLETED.
C***
t***
C WRITE THE HILL ID, THE NUMBER OF CRITICAL ELEVATIONS, THE RILL
C TOP ELEVATION, AND THE HILL NAME TO THE CTDM INPUT FILE.
1C***
C***
! WRITE (MOUT, 260) IDHILL,NCR,HTOP,HNAME
! 260 FORMAT(5X,I2,1X,I2,10X,E10.4,10X,A15) '
C***
I
FOR EACH CRITICAL ELEVATION, DETERMINE THE PARAMETERS WHICH BEST
DESCRIBE THE ELLIPTICAL TERRAIN CONTOUR AT THAT ELEVATION. THESE
HCT02870
HCT02880
HCT02890
HCT02900
HCT02910
HCT02920
HCT02930
HCT02940
HCT02950
HCT02960
HCT02970
HCT02980
HCT02990
HCT03000
HCT03010
HCT03020
HCT03030
HCT03040
HCT03050
HCT03060
HCT03070
HCT03080
HCT03090
HCT03100
HCT03110
HCT03120
HCT03130
HCT03140
HCT03150
HCT03160
HCT03170
HCT03180
HCT03190
HCT03200
PARAMETERS ARE WRITTEN TO THE CTDM INPUT FILE FOR USE IN THE "WRAP" HCT03210
C
e
e
e
c***
c***
PLUMB CALCULATION IN CTDM. 17 THE CRITICAL ELEVATION DOES NOT CO-
INCIDE WITH AN INPUT CONTOUR(I.E. ICHMOD-2), THEN THE PARAMETERS
MUST BE DETERMINED BY A SIMPLE INTERPOLATION OF FITTED CONTOUR
PARAMETERS BASED ON ELEVATION. THE INTERPOLATION OF THE OREN-
TATION VALUES IS A VECTOR INTERPOLATION WITH THE VECTORS WEIGHTED
WITH THE ECCENTRICITY OF THE CONTOUR.
I7(ICHMOD.EQ.2) GO TO 290
HCT03220
HCT03230
HCT03240
HCT03250
HCT03260
HCT03270
HCT03280
HCT03290
HCT03300
148
-------
C***CASE I—CRITICAL ELEVATIONS COINCIDE WITH CONTOUR ELEVATIONS. HCT03310
DO 280 J-1,NCR HCT03320
C***FIND THE ORIENTATION OF THE MAJOR AXIS MEASURED CLOCKWISE FROM HCT03330
C***NORTH(BETWEEN 0 AND 180 DEGREES). HCT03340
ONOR-ISO.-OREN(J) HCT03350
IF(ONOR.LT.O.) ONOR-360.+ONOR HCT03360
IF(ONOR.GT.180.) ONOR-ONOR-180. HCT03370
WRITE (MOOT, 270) HC (J)-, XCH (J) ,YCH(J) ,ONOR,A(J) ,B(J) HCT03380
270 FORMAT(F10.3,2E10.4,3F10.3) HCT03390
280 CONTINUE HCT03400
' GO TO 360 HCT03410
C***CASE 2—CRITICAL ELEVATIONS EVENLY SPACED BETWEEN HCLOW AND THE HCT03420
C***UPPERMOST CONTOUR HCT03430
290 DO 350 1-1,NCR HCT03440
DO 300 J-1,NC HCT03450
JK-J HCT03460
IF(HCON(J).GT.HC(I)) GO TO 310 HCT03470
300 CONTINUE HCT03480
310 IF(JK.GT.l) GO TO 320 HCT03490
C***IF THE CRITICAL ELEVATION IS BELOW THE LOWEST CONTOUR, THEN • HCT03500
C***EXTRAPOLATE THE VALUES FOR THE CONTOUR ORIENTATION, CENTROID HCT03510
C***COORDINATES, AND SEMI-MAJOR AND SEMI-MINOR AXIS LENGTHS USING THE HCT03520
C***VALUES OF THESE PARAMETERS FOR THE LOWEST TWO CONTOURS. IF THERE HCT03530
:C***IS ONLY ONE CONTOUR, THEN THE VALUES FOR THE ORIENTATION AND HCT03540
C***CENTROID COORDINATES OF THE CRITICAL ELEVATION CONTOUR ARE SET HCT03550
C***EQUAL TO THE CORRESPONDING VALUES FOR THE SINGLE CONTOUR. THE HCT03560
C***SEMI-MAJOR AND SEMI-MINOR AXIS LENGTHS FOR THE CRITICAL ELEVATION HCT03570
C***CONTOUR ARE EXTRAPOLATED BY ASSUMING A ZERO AREA CONTOUR AT THE HCT03580
C***HILL TOP ELEVATION. HCT03590
IF(MC.EQ.l) GO TO 315 HCT03600
JK-2 HCT03610
GO TO 320 HCT03620
315 XOd-XCM(l) HCT03630
- YCMI-YCM(l) HCT03640
ORENI-OREN(l) HCT03650
FEXT-(HTOP-HC(I))/(HTOP-HCON(1)) HCT03660
AI-A(1)*FEXT HCT03670
BI-B(1)*FEXT HCT03680
GO TO 340 HCT03690
C***INT£RPOLATE TO FIND CONTOUR PARAMETERS AT THE Ith CRITICAL HCT03700
C***ELEVATION. HCT03710
320 FRACT-(HC(I)-HCON(JK-l))/(HCON(JK)-HCON(JK-l)) HCT03720
XCMI-XCM(JK-l)+FRACT*(XCM(JK)-XCM(JX-l)) HCT03730
YCMI-YCM(JK-l)+FRACT*(YCM(JK)-YCM(JK-l)) HCT03740
AI-A(JK-l)+FRACT*(A(JK)-A(JK-l)) HCT03750
BI-B(JK-1)+FRACT*(B(JK)-B(JK-l)) HCT03760
C***DO NOT ALLOW AI AND BI TO DECREASE WITH ELEVATION. HCT03770
IF(AI.LT.A(JK-1)) AI-A(JK-l) HCT03780
IF(BI.LT.B(JK-1)) BI-B(JX-l) HCT03790
C***INTERPOLATB THE ORIENTATION VECTORIALLY WITH THE ELLIPSE HCT03800
C***ECCENTRICITY USED AS A WEIGHTING FACTOR. HCT03810
SUMX-ECC(JK-1)*COS(PI*OREN(JK-1)/180.)+FRACT*(ECC(JK)* HCT03820
4COS(PI*OREM(JR)/180.)-ECC(JK-l)*COS(PI*OREN(JK-1)/180.)) HCT03830
SUMY«ECC(JK-l)*SIH(PI*OREN(JK-l)/180.)-fFRACT*(ECC(JR)* HCT03840
«SIH(PI*OREN(JK)/180.)-ECC(JK-1)*SIN(PI*OREN(JK-1)/180.)) HCT03850
C***AVOID CALLING THE ATAN2 FUNCTION WITH BOTH ARGUMENTS BEING HCT03860
C***EFFECTIVELY ZERO. . HCT03870
IF(ABS(SUMX).LT.1.0E-8.AND.ABS(SUMY).LT.l.OE-8) GO TO 330 HCT03880
ORENI-(180./PI)*ATAN2(SUMY,SUMX) HCT03890
GO TO 340 HCT03900
330 ORENI-O. HCT03910
C***IF THE EXTRAPOLATION PROCESS GIVES AN ELLIPSE WITH A MINOR AXIS HCT03920
C***GREATER THAN A MAJOR AXIS, THEN ASSUME THAT THE AXES ARE EQUAL HCT03930
C***AND THAT THE ELLIPSE HAS THE SAME AREA. HCT03940
340 IF(AI.GE.BI) GO TO 345 HCT03950
AI-SQRT(AI*BI) HCT03960
149
-------
BI-AI
: 345 CONTINUE
C***FIND THE ORIENTATION OF THE INTERPOLATED CONTOUR MAJOR AXIS AS
C***MEASURED CLOCKWISE FROM NORTH ( BETWEEN 0 AND 180 DEGREES).
! OHOR-180.-ORENI
IF(ONOR.LT.O.) ONOR-360.+ONOR
IFfONOR.GT.180.) ONOR-ONOR-180.
WRITE (MOOT, 270) HC(I) ,XCMI,YCMI,ONOR,AI,BI
350 CONTINUE
! 360 CONTINUE
C***THE WRITING OF BEST FIT CONTOUR ELLIPSE PARAMETERS FOR CUTOFF
C***ELEVATIONS TO THE CTDM INPUT FILE HAS BEEN COMPLETED.
C***
'C DETERMINE THE FITTED HILL PARAMETERS FOR EACH CRITICAL CUTOFF
C ELEVATION AND WRITE THESE PARAMETERS TO BOTH THE PLOT FILE AND
C THE CTDM INPUT FILE.
c***
; DO 500 I-1,NCR
C***ZERO OUT SUMMATION VARIABLES.
SUM1-0.
; SUM2A-0.
SUM2B-0.
; SUM3-0.
SUM4A-0.
! SUM4B-0.
; SDMX-0.
: SUMY-0.
XHTOPF— 0 . »
\ YHTOPF-0.
NCOH-0
;C***CALCULATE THE HILL HEIGHT ABOVE THE CRITICAL HEIGHT.
; HHILL-HTOP-HC(I)
DO 400 J-l.NC
C***CONTOUR ELEVATIONS USED IN FITTING THE PORTION OF THE HILL ABOVE
C***THE CRITICAL ELEVATION MUST BE AT LEAST ONE UNIT ABOVE THE CRITICAL
C***ELEVATION.
< IF(HCON(J).LE.HC(I)+1.) GO TO 400
BCON-NCON+1
FJ-ALOG (HHILL/ (HCON ( J) -HC ( I) ) -1 . )
SUM1-SUM1+FJ
SUM3»SUM3+FJ**2
SUM2A-SUH2A+ALOG(A(J))
; SUM2B-SUM2B+ALO<3(B(J)
SUM4A-SUM4A+ALOG(A(J) ,
SUM4B-SUM4B+ALOG(B(J)
*FJ
*FJ
SUMX-StJMX+ECC(J) *COS(PI*OREN(J)/180.)
SUMY-SUMY+ECC(J)*SIN(PI*OREN(J)/180.)
XHTOPF-XHTOPF-I-XCM(J)
YHTOPF-YHTOPF+YCM(J)
400 CONTINUE
: IF(NCON.EQ.l) GO TO 410
LA-EXP((SUM2A*SUM3-SUM4A*SUM1)/(NCON*SUM3-SUM1**2))
! LB-EXP({SUM2B*SUM3-SUM4B*SUM1)/(NCON*SUM3-SUM1**2))
PA-(NCON*SUM3-SUM1**2)/(NCON*SDM4A-SUM1*SUM2A)
! PS- (NCON*SUM3-SUM1**2 ) / (NCON*SUM4B-SOMl*StJH2B)
C***NEGATTVE EXPONENTS NOT ALLOWED
! PA-ABS(PA)
I PB-ABS(PB)
GO TO 420
C***IF ONLY ONE CONTOUR IS USED IN THE RILL FIT, ONE MUST ASSUME
C***THAT THE EXPONENTS IN'THE INVERSE POLYNOMIAL FIT ARE BOTH 2.
410 CONTINUE
PA-2.
PB-2.
HCT03970
HCT03980
HCT03990
HCT04000
HCT04010
HCT04020
HCT04030
HCT040-TO
HCT04050
HCT04060
HCT04070
HCT04080
HCT04090
HCT04100
HCT04110
HCT04120
HCT04130
HCT04140
HCT04150
HCT04160
HCT04170
HCT04180
HCT04190
HCT04200
HCT04210
HCT04220
HCT04230
HCT04240
HCT04250
HCT04260
HCT04270
HCT04280
HCT04290
HCT04300
HCT04310
HCT04320
HCT04330
HCT04340
HCT04350
HCT04360
HCT04370
HCT04380
HCT04390
HCT04400
HCT04410
HCT04420
HCT04430
HCT04440
HCT04450
HCT04460
HCT04470
HCT04480
HCT04490
HCT04500
HCT04510
HCT04520
HCT04530
HCT04540
HCT04550
HCT04560
HCT04570
HCT04580
HCT04590
HCT04600
HCT04«10
HCT04620
ISO
-------
LA-A(NC)/(HHILL/(HCON(NC)-HC(I))-!.)**(!-/PA) HCT04630
LB-B(NC)/(HHILL/(HCON(NC)-HC(I))-1.)**(1./PB) HCT04640
C***AVOID CALLING THE ATAN2 FUNCTION WITH BOTH ARGUMENTS BEING. HCT04650
C***EFFECTIVELY ZERO. HCT04660
420 IF(ABS(SUHX).LT.l.OE-8.AND.ABS(SUMY).LT.l.OE-8) GO TO 430 HCT04670
ORENF-(180./PI)*ATAN2(SUMY,SUMX) . HCT04680
GO TO 440 . HCT04690
430 ORENF-0. HCT04700
C***FIND THE ORIENTATION OF THE MAJOR AXIS AS MEASURED CLOCKWISE FROM HCT04710
C***NORTH(BETWEEN O AND 180 DEGREES). HCT04720
440 ONOR-180.-ORENF . HCT04730
IF(ONOR.LT.O.) ONOR-360.+ONOR HCT04740
IF(ONOR.GT.180.) ONOR-ONOR-180. HCT04750
XHTOPF-XHTOPP/FLOAT(NCON) HCT04760
YHTOPF-YHTOPF/FLOAT(NCON) HCT04770
IF(PFLAG.EQ.O) GO TO 455 HCT04780
C***WRXT£ THE FITTED HILL PARAMETERS TO THE PLOT FILE. HCT04790
WRITE(UPL,450) HC(I),XHTOPF,raTOPF,ORENF,PA,PB,LA,LB HCT04800
450 FORMAT(8E15.4) HCT04810
455 CONTINUE HCT04820
C***WRTTE THE FITTED HILL PARAMETERS TO THE CTDM INPUT FILE. HCT04830
WRITE(MOOT,460) HC(I),XHTOPF,YHTOPF,ONOR,PA,PB,LA,LB HCT04840
460 FORMAT(F10.3,2E10.4,5F10.3) HCT04850
500 CONTINUE HCT04860
STOP HCT04870
END HCT04880
151
-------
SUBROUTINE PSORTR( ARRAY, NDL,NPTR,LPTR)
C***POINTER SORT USING THE MERGE EXCHANGE METHOD
C***NUMBER OF COMPARISONS-N*LOG(N)/LOG(2)
C***ARRAY-REAL ARRAY TO BE SORTED
C***NDL-NUMBER OF ELEMENTS OF ARRAY TO BE SORTED
C***NPTR-POINT£R ARRAY
C***LPTR-WORKING ARRAY
DIMENSION ARRAY(1),NPTR(1) ,LPTR(1)
C***CHECK INITIAL ORDER
Il-NPTR(l)
IF(MDL.LE.l.AND.Il.EQ.l) RETURN
IF(Il.LT.l.OR.Il.GT.NDL) GO TO 30
DO 20 1-2, NDL
I2-NPTR(I)
IF(I1.EQ.I2) GO TO 30
IF(I2.LT.1.OR.I2.GT^NDL) GO TO 30
rP(ARRAY(Il) .GT.ARRAY(I2)) GO TO 30
11-12
20 CONTINUE
RETURN
C***SET UP POINTER ARRAY
30 DO 40 I-1,NDL
MPTR(I)-I
1 40 CONTINUE
C***BEGIN THE SORT
IF(NDL.LE.l) RETURN
L2I-1
DO 120 1-1,20
M-l
L2IH-L2I
L2I-2*L2I
IF(L2IH.GT.NDL) GO TO 130
JUP-NDL/L2I+1
DO 110 J-1,JUP-
N-H+L2IH
IF(N.GT.NDL) .GO TO 110
KLO-M
XUP-iflNO (KLO+L2I-1,NDL)
MUP-XLO+L2IH-1
DO 80 K-KLO,KUP
IF(N.GT.NDL) GO TO 50
IF(N.GT.KUP) GO TO 50
IF(M.GT.MUP) GO TO 60
»( ARRAY (HPTR(M)).GT. ARRAY (NPTR(N) ) ) GO TO 60
50. NL-M
M-M+1
GO TO 70
60 WL-N
70 LPTR(K)-NI»TR(NL)
80 CONTINUE
90 DO 100 K-KLO,KUP
VPTR(X)-LPTR(K)
100 CONTINUE
M-KLO+L2I
110 CONTINUE
IF(L2I.GE.NDL) GO TO 130
120 CONTINUE
130 RETURN
END
PS000010
PSO00020
PSO00030
PS000040
PS000050
PS000060
PSO00070
PSO00080
PSO00090
PSO00100
PSO00110
PS000120
PS000130
PS000140
PS000150
PS000160
PS000170
PS000180
PSO00190
PS000200
PS000210
PS000220
PSO00230
PS000240
PSO00250
PS000260
PS000270
PS000280
PS000290
PS000300
PS000310
PSO00320
PSO00330
PS000340
PSO00350
PSO00360
PSO00370
PSO00380
PSO00390
PS000400
PSO00410
PSO00420
PSO00430
PSO00440
PSO00450
PS000460
PSO00470
PSO00480
PSO00490
PSO00500
PSO00510
PSO00520
PSO00530
PSO00540
PSO00550
PSO00560
PSO00570
PSO00580
PSO00590
PSO00600
152
-------
PLOTCOH
153
-------
10 'Program to plot contours for actual and fitted hills on a display
20 'terminal with 320(horizontal)x200(vertical) resolution in color
30 'or 640(horizontal)x200(vertical) resolution in black and white
40 'Clear the screen.
50 CTiP
60 'Disable the display of function keys to allow more space for
70 'plotting.
80 KEY OFF
90 DEFINT I-N
100 'Dimension the arrays for contour elevations, contour identification
110 'numbers(from both FITCON and HCRIT), and the array for storing the
120 'plot of digitized contours(unedited or edited).
130 DIM HCON(200),IDC1(200),IDC2(200),IAR(8002)
140 LOCATE 12,15
150 'Input the name of the plot file from program FITCON.
160 INPUT " INPUT NAME OF PLOTFILE FROM PROGRAM FITCON—>";PLOT1$
170 ON ERROR GOTO 3190
180 OPEN PLOT1$ FOR INPUT AS fl
190 ON ERROR GOTO 0
200 'Make sure that this plot file was generated by program FITCON.
210 INPUT!l, PF1$
220 IF PF1$-"FITCON" THEN GOTO 280
230 CLS
240 LOCATE 10,15
250 PRINT PLOT1S « IS NOT A FILE GENERATED BY PROGRAM FITCON-TRY AGAIN1*
260 CLOSE fl
270 GOTO 140
280 CLS
290 'Input the hill identification number, hill name, hill center
300 'coordinates, number of fitted contours, and the identification
310 'numbers for the fitted contours.
320 INPUT!1, IDH1,HNAME1$
330 INPUT*1, XHTOP,YHTOP
340 INPUT*1, NCI
350 FOR J-l TO NCI
360 INPUT*1,IDC1(J)
370 NEXT J
380 'Input the plot boundaries for the unedited contours.
390 INPUT*1, XMIN1,XMAX1,YMIN1,YMAX1
400 'Input the plot boundaries for the edited contours.
410 INPUT*1, XMIN2,XMAX2,YMIN2,YMAX2
420 LOCATE 10,22
430 'Select the type of display.
440 PRINT "SELECT TYPE OF DISPLAY*
450 PRINT
460 PRINT TAB(22) "1.) Low resolution with color"
470 PRINT TAB(22) "2.) High resolution black and white"
480 PRINT
490 INPUT » Choice? (l o^r 2)—>";RFLAG*
500 CLS
510 LOCATE 10,22
520 'Select the type of contours to be displayed.
530 PRINT "SELECT THE CONTOUR TYPE FOR DISPLAY"
540 PRINT
550 PRINT TAB(22) "1.) Unedited Contours"
560 PRINT TAB(22) "2.) Edited Contours"
570 PRINT
580 INPUT « Choice?(1 or 2)—>";DFLAG%
590 CLS
600 'Set plot boundaries, scale factors, and colors.
154
-------
610 SCRCX-320!:DSCRX-468!:SCRCY-104!:DSCRY-190!:RATIO-1.3574
620 IF RFLAG*-! THEN SCRCX-1601:DSCRX»2051:RATIO-1.5437
630 IF DFLAG%-2 THEN GOTO 690
640 XO(XMINl+XMAXl)/2!
650 YO(YMINH-YMAXl)/2!
660 DX-XMAX1-XMIN1
670 DY-YMAX1-YMIN1
680 GOTO 730
690 XO(XMIN2+XMAX2)/2!
700 YC-(YMXN2+YMAX2)/2!
710 DX-XMAX2-XMIN2
720 DY-YMAX2-YMIN2
730 IF DX/DY-2 AND ABS(X-Xl)<1E-15 AND ABS(Y-Yl)<1E-15 THEN GOTO 1310
1120 'Scale the point X,Y for plotting.
1130 XS-SCRCX+(X-XC)*DSCRXr>DD
1140 YS«SCRCY-(Y-YC)*DSCRYDDD
1150 IF DUPFLG%»0 GOTO 1270
1160 'One of the multiple contours has been closed. Move to the new point
1170 'without drawing a line. Substitute the current point for the
1180 'previous individual contour beginning point.
1190 XOLD-X
1200 YOLD-Y
155
-------
1210 DUPFLG*-0
1220 PSET(XS,YS),IC
1230 GOTO 1310
1240 'Determine whether one of the individual multiple contours has been
1250 'closed. If so, set the closure indicator DUPFLG* to 1 and
1260 'increment the contour closure counter IFR by 1.
1270 IP ABS(X-XOLD)<1E-15 AND ABS(Y-YOLD)<1E-15 THEN DUPFLG*-1:IFR-IFR+1
1280 'Draw a line from the previous point to the current point.
1290 LINE -(XS,YS),IC
1300 'End loop over contour points.
1310 NEXT K
1320 I? DFLAG%<> 1 THEN GOTO 1390
1330 'Skip over edited contours.
1340 INPOTfl, NPC
1350 FOR K-l TO NPC
1360 INPUTfl,DUMX,DUMY
1370 NEXT K
1380 'End loop over contours.
1390 NEXT J
1400 'Scale hill center coordinates.
1410 XSHOSCRCX+(XHTOP-XC)*DSCRXDDD
1420 YSHC»SCRCY-(YHTOP-YC)*DSCRYDDD
1430 XUL-XSHC-1
1440 XLR-XSHC+1
1450 YUL-YSHC-1
1460 YLR-YSHC+1
1470 'Plot a 3x3 box of points centered at the hill center.
1480 LINE(XUL,TOL)-(XLR,YLR),IC,BF
1490 I? RFLAG*-! THEN GXMX%»319 ELSE GXMX%»639
1500 'Store the plot of digitized contours in array IAR.
1510 GET (0,0)-(GXMX%,199),IAR
1520 PRINT HNAMEIS • INPUT CONTOURS'1
1530 'Pause until user presses any key. Program will terminate if the
1540 'user presses the ESC key.
1550 GOSUB 3410
1560 CLS
1570 'Change color to magenta for plotting fitted contours.
1580 IF RFLAG%-1 THEN IC-2
1590 'Restore the plot of digitized contours.
1600 POT (0,0),IAR,PSET
1610 'Begin loop over contours.
1620 FOR J-l TO NCI
1630 'Input ellipse parameters for each contour: ellipse centroid
1640 'coordinates, semi-axes lengths, and the orientation of the minor
1650 'axis with respect to the positive x-axis.
1660 INPUT!1, XCM,YCM,A,B,OR£N
1670 'Determine the orientation of the major axis with respect to the
1680 'positive x-axis
1690 OREN-OREN-901
1700 CSE-COS(.017453*OREN)
1710 SNE-SIN(.017453*OREN)
1720 XP-A
1730 XFIT-XCM+XP*CSE
1740 YFIT-YCM+XP*SNE
1750 XS-SCRCX+(XFIT-XC)*DSCRXDDD
1760 YS-SCRCY-(YFIT-YC)*DSCRYDDD
1770 'Move to a point at the end of the ellipse semi-major axis.
1780 PSET(XS,YS)
1790 A2-AA2
1800 B2-B*2
156
-------
1810 'Draw an ellipse with 120 points.
1820 FOR L-l TO 120
1830 THO-L*.05276
1840 R-SQR(1!/(COS(THC)*2/A2+SIN(THC)*2/B2))
1850 XP-R*COS(THC)
1860 YP-R*SIN(THC)
1870 XFIT-XCM+XP*CSE-YP*SNE
1880 YFIT-YCM+XP*SNE+YP*CSE
1890 XS-SCRCX+(XFIT-XC)*DSCRXDDD
1900 YS-SCRCY-(YFIT-YC)*DSCRYDDD
1910 LINE -(XS,YS),IC
1920 NEXT L
1930 'End loop over contours.
1940 NEXT J
1950 PRINT HNAME1$ " FITTED CONTOURS*
1960 'Pause until the user presses any key. If the user presses the ESC
1970 'key, then program execution will terminate.
1980 GOSUB 3410
1990 CLS
2000 'Begin plotting contours for fitted cutoff hills.
2010 'Go to text mode for user input.
2020 SCREEN 2:SCREEN 0
2030 LOCATE 12,19
2040 'Determine whether fitted hill contours are to be displayed.
2050 INPUT " DISPLAY FITTED CUTOFF HILL CONTOURS?(Y/N)->";ANS$
2060 IF ANS$-"N" THEN SYSTEM
2070 IF ANS$-"n" THEN SYSTEM
2080 CLOSE fl
2090 LOCATE 14,15
2100 'Input the name of the plot file from program RCRIT.
2110 INPUT " INPUT NAME OF PLOTFILE FROM PROGRAM HCRIT";PLOT2$
2120 ON ERROR GOTO 3220
2130 OPEN PLOT2$ FOR INPUT AS fl
2140 ON ERROR GOTO 0
2150 'Make sure the plot file was generated by program HCRIT.
2160 INPUTfl, PF2$
2170 IF PF2$-"HCRIT" THEN GOTO 2230
2180 CLS
2190 LOCATE 12,20
2200 PRINT PLOT2$ " IS NOT A FILE GENERATED BY PROGRAM HCRIT-TRY AGAIN"
2210 CLOSE fl
2220 GOTO 2090
2230 CLS
2240 'Check whether the hill identification number, hill name, number
2250 'of fitted contours, and contour identification numbers match
2260 'those from the FTTCON plot file.
2270 INPUTfl, IDH2,HNAME2$
2280 17 IDH20IDH1 THEN GOTO 3250
2290 IF HNAME1$<>HNAME2$ THEN GOTO 3280
2300 INPUTfl, NC2
2310 IF NC10NC2 THEN GOTO 3310
2320 FOR J-l TO NC2
2330 INPUTfl, IDC2(J)
2340 IF IDC1(J)OIDC2(J) THEN GOTO 3340
2350 NEXT J
2360 'Return to graphics mode.
2370 IF RFLAG4-1 THEN SCREEN 1:IO2 ELSE SCREEN 2:IC-1
2380 IF RFLAG*»1 THEN COLOR 9,1
2390 'Input hill top elevation and contour elevations.
2400 INPUTfl, HTOP
157
-------
2410 FOR J-l TO NC2
2420 INPUT*!, HCON(J)
2430 NEXT J
2440 'Input number of critical elevations.
2450 INPUT*1, NCR
2460 'Begin loop over critical elevations.
2470 FOR 1-1 TO NCR
2480 'For each critical elevation, input the critical elevation, cutoff
2490 'hill centroid coordinates, orientation of the hill minor axis
2500 'with respect to the positive x-axis, and the inverse polynomial
2510 'fit parameters for each hill axis.
2520 INPUTtl, HC,XHTOPF,YHTOPF,ORENF,PA,PB,RIA,RLB
2530 'Determine the orientation of the major axis with respect to the
2540 'positive x-axis.
2550 ORENF-ORENF-90!
2560 CSE-COS(.017453*OR£NF)
2570 SNE-SIN(.017453*ORENF)
2580 'Retrieve background plot of digitized contours(unedited or edited).
2590 PUT (0,0),IAR,PSET
2600 'Begin loop over contours.
2610 FOR J-l TO NC2
2620 'Contours must be at least one elevation unit higher than the
2630 'critical elevation if their elevations are to be used for the display
2640 'of contours on the cutoff hill.
2650 IF HCON(J)<-HC+1! THEN GOTO 3020
2660 FLOG-LOG((HTOP-HC)/(HCON(J)-HC)-1!)
2670 AFIT-RIA*EXP((1!/PA)*FLOG)
2680 BFIT-RLB*EXP((11/PB)*FLOG)
2690 'The equation for the inverse polynomial contour is
2700 • (XP/AFIT)**PA+(YP/BFIT)**PB-1
2710 'in the coordinate system in which the x and y primed axes
2720 'coincide with the major and minor axes of the hill respectively.
2730 'Begin loop to calculate 800 contour point coordinates.
2740 FOR L-l TO 200
2750 17 L>99 GOTO 2810
2760 'Let x primed be the independent variable.
2770 XPOL-L*.01*AFIT
2780 YPOL-BFIT*(ll-(XPOL/AFrr)APA)A(ll/PB)
2790 GOTO 2840
2800 'Let y primed be the independent variable.
2810 YPOL-(L-100)*.01*BFIT
2820 XPOL-AFIT*(ll-(YPOVBFIT)APB)*(l!/PA)
2830 'First quadrant(x primed-*,y primed-*)
2840 XP-XPOL
2850 YP-YPOL
2860 GOSUB 3460
2870 'Second quadrant (x primed-*, y primed—)—moving clockwise.
2880 XP-XPOL
2890 YP— YPOL
2900 GOSUB 3460
2910 'Third quadrant (x primed—,y primed—)
2920 XP—XPOL
2930 YP—YPOL
2940 GOSUB 3460
2950 'Fourth quadrant (x primed—,y primed-*)
2960 XP—XPOL
2970 YP-YPOL
2980 GOSUB 3460
2990 'End contour point loop.
3000 NEXT L
158
-------
3010 'End contour loop.
3020 NEXT J
3030 XSHCF-SCRCX+(XHTOPF-XC)*DSCRXDDD
3040 YSHCF-SCRCY-(YHTOPF-YC)*DSCRYDDD
3050 XUL-XSHCF-1
3060 XLR-XSHCF+1
3070 YUL-YSHCF-1
3080 YLR-YSHCF+1
3090 'Plot a 3x3 box of points centered about the cutoff hill centroid.
3100 LINE (XUL,YUL)-(XLR,YLR)»IC,BF
3110 PRINT HNAME2$ " ECRIT-" HC
3120 'Pause until user strikes a key. If the ESC key is pressed, then
3130 'execution of the program is terminated.
3140 GOSUB 3410
3150 CLS
3160 'End loop on critical elevations.
3170 NEXT I
3180 SYSTEM
3190 IP ERR-53 THEN PRINT "FITCON PLOT FILE NOT FOUND-Press any key"
3200 GOSUB 3410
3210 SYSTEM
3220 17 ERR-53 THEN PRINT "HCRIT PLOT FILE NOT FODND-Press any key"
3230 GOSUB 3410
3240 SYSTEM
3250 PRINT "FITCON AND HCRIT RILL IDs DO NOT MATCH-Press any key"
3260 GOSUB 3410
3270 SYSTEM
3280 PRINT "FITCON AND HCRIT HILL NAMES DO NOT MATCH-Press any key"
3290 GOSUB 3410
3300 SYSTEM
3310 PRINT "FITCON AND HCRIT NUMBER OF CONTOURS DO NOT MATCH-Press any key"
3320 GOSUB 3410
3330 SYSTEM
3340 PRINT "FITCON AND 'HCRIT CONTOUR IDs DO NOT MATCH-Press any key"
3350 GOSUB 3410
3360 SYSTEM
3370 END
3380 'Subroutine which causes program execution to pause until a key
3390 'is struck. If the ESC key is pressed, then program execution
3400 'will be terminated.
3410 A$-INKEY$: 17 A$-"« THEN 3410
3420 17 A$-CHR$(27) THEN SYSTEM
3430 RETURN
3440 'Subroutine to rotate points into the xfy coordinate system before
3450 'plotting
3460 X7IT-XHTOP7+XP*CSE-YP*SNE
3470 Y7IT-YHTOP7+XP*SNE+YP*CSE
3480 XS-SCRCX+(XFIT-XC)*DSCRXDDD
3490 YS-SCRCY-(Y7IT-YC)*DSCRYDDD
3500 PSET(XS^YS),IC
3510 RETURN
159
-------
HPLTCOH
160
-------
10 'Program to plot contours for actual and fitted hills on a display
20 'terminal with 720(horizontal)x348(vertical) resolution(black and
30 'white) driven by a Hercules Graphics Board
40 'Clear the screen.
50 CLS
60 'Disable the display of function keys to allow more space for
70 'plotting.
80 KEY OFF
90 DEFINT I-N
100 'Dimension the arrays for contour elevations, contour identification
110 *numbers(from both FITCON and HCR1T), and the array for storing the
120 'plot of digitized contours(unedited or edited).
130 DIM HCOH(200),IDC1(200),IDC2(200),IAR(15662)
140 LOCATE 12,15
150 'Input the name of the plotfile from program FITCON
160 INPUT " INPUT NAME OF PLOTFILE FROM PROGRAM FITCON—>";PLOT1$
170 ON ERROR GOTO 2940
180 OPEN PLOT1$ FOR INPUT AS 11
190 ON ERROR GOTO 0
200 INPUT!1, PF1$
210 IF PF1$-"FITCON" THEN GOTO 270
220 CLS
230 LOCATE 10,15
240 PRINT PLOT1$ " IS NOT A FILE GENERATED BY PROGRAM FITCON-TRY AGAIN"
250 CLOSE f-1
260 GOTO 140
270 CLS
280 'Input the hill identification number, hill name, hill center
290 'coordinates, number of fitted contours, and the identification
300 'for the fitted contours.
310 INPUT*1, IDH1,HNAME1$
320 INPUT*1, XHTOP,YHTOP
330 INPUT*1, NCI
340 FOR J-l TO NCI
350 INPUT*1,IDC1(J)
360 NEXT J
370 'Input the plot boundaries for the unedited contours.
380 INPUT*!, XMIN1,XMAX1,YMIN1,YMAX1
390 'Input the plot boundaries for the edited contours.
400 INPUT*!, XMIN2,XMAX2,YMIN2,YMAX2
410 'Set plot boundaries and scale factors.
420 SCRCX-3601:DSCRX-499!:SCRCY-180!:DSCRY-333I:RATIO-1.4286
430 CLS
440 LOCATE 10,22
450 'Select the type of contours to be displayed.
460 PRINT "SELECT THE CONTOUR TYPE FOR DISPLAY"
470 PRINT
480 PRINT TAB(22) "1.) Unedited Contours"
490 PRINT TAB(22) "2.) Edited Contours"
500 PRINT
510 INPUT " . Choice?(1 or 2)—>";DFLAG%
520 CLS
530 IF DFLAG*»2 THEN GOTO 590
540 XO(XMIN1+XMAX1)/21
550 YO(YMINl+YMAXl)/2!
560 DX-XMAX1-XMIN1
570 DY-YMAX1-YMIN1
580 GOTO 630
590 XO(XMIN2+XMAX2)/2!
600 YO(YMIN2+YMAX2)/2!
161
-------
610 DX-XMAX2-XMIN2
620 DY-YMAX2-YMIN2
630 IF DX/DY-2 AND ABS(X-X1)<1E-15 AND ABS(Y-Y1)<1E-15 THEN GOTO 1170
980 'Scale the point x,y for plotting.
990 XS-SCRCX+(X-XC)*DSCRXDDD
1000 YS-SCRCY-(Y-YC)*DSCRYDDD
1010 IF DUPFLG%-0 GOTO 1130
1020 'One of the multiple contours has been closed. Move to the new point
1030 'without drawing a line. Substitute the current point for the
1040 'previous individual contour beginning point.
1050 XOLD-X
1060 YOLD-Y
1070 DUPFLG%-0
1080 PSET(XS,YS)
1090 GOTO 1170
1100 'Determine whether one of the individual multiple contours has been
1110 'closed. If so, set the closure indicator DUPFLG% to 1 and
1120 'increment the contour closure counter IFR by 1.
1130 IF ABS(X-XOLD)<1E-15 AND ABS(Y-YOLD)<1E-15 THEN DUPFLG%-1:IFR«IFR+1
1140 'Draw a line from the previous point to the current point.
1150 LINE -(XS,YS)
1160 'End loop over contour points.
1170 NEXT K
1180 IF DFIAG*01 THEN GOTO 1250
1190 'Skip over edited contours.
1200 INPUTfl,NPC
162
-------
121C FOR K-l TO NPC
1220 INPUT!1 DUMX,DUMY
1230 NEXT K
1240 'End loop over contours.
1250 NEXT J
1260 'Scale hill center coordinates.
1270 XSHC-SCRCX+ (XHTOP-XC) *DSCRXODD
1280 YSHOSCRCY- (YHTOP-YC) *DSCRYDDD
1290 XUL-XSHC-1
1300 XLR-XSHC+1
1310 YUL-YSHC-1
1320 YLR-YSHC-P1
1330 'Plot a 3x3 box of points centered at the hill center.
1340 LZNE(XUL,YUL)-(XLR,¥LR),,BF
1350 'Store the plot of digitized contours in array IAR.
1360 GET (0,0)-(719,347),IAR
1370 PRINT HNAME1$ " INPUT CONTOURS"
1380 'Pause until user presses a key. Program will terminate if the
1390 'user presses the ESC key.
1400 GOSUB 3160
1410 CLS ' .
1420 'Restore the plot of digitized contours.
1430 PUT (0,0),IAR,PSET
1440 'Begin loop over contours.
1450 FOR J-l TO NCI
1460 'Input ellipse parameters for each contour: ellipse centreid
1470 'coordinates, semi-axes lengths, and the orientation of the minor
1480 'axis with respect to the positive x-axis.
1490 INPUT*1, XCM,YCM,A,B,OR£N
1500 'Determine the orientation.of the ellipse major axis with respect
1510 'to the positive x-axis.
1520 OREN-OREN-901
1530 CSE-COS(.017453*OREN)
1540 SNE-SIN<.Or7453*OR£N)
1550 XP-A
1560 XFIT-XCM+XP*CSZ
1570 YFIT-YCM+XP*SNE
1580 XS-SCRCX+(XFXT-XC)*DSCRXDDD
1590 YS«SCRCY-(YFIT-YC)*DSCRYDDD
1600 'Move to a point at the end of the ellipse semi-major axis.
1610 PSET(XS,YS)
1620 A2-AA2
1630 B2»BA2
1640 'Draw an ellipse with 120 points.
1650 FOR L-l TO 120
1660 THC—L*. 05276
1670 R-SQR(li/(COS(THC)A2/A2+SIN(THC)*2/B2))
1680 XP-R*COS(THC)
1690 YP-R*SIN(THC)
1700 XFIT-XCM+XF*CSE-YI>*SNE
1710 YFIT-YCK+XP*Sl!E+YP*CS£
1720 XS-SCRCX-l-(X7IT-XC)*DSCRXDDD
1730 YS-SCRCY-(YFIT-YC)*DSCRYDDD
1740 LINE -(XS,YS)
1750 NEXT L
1760 'End loop over contours.
1770 NEXT J
1780 PRINT HNAME1$ " FITTED CONTOURS"
1790 'Pause until the user presses a key. If the user presses the ESC
1800 'key, then program execution will terminate.
163
-------
1810 GOSUB 3160
1820 CLS
1830 'Begin plotting contours for fitted cutoff hills.
1840 LOCATE 12,19
1850 'Determine whether fitted hill contours are to be displayed..
1860 INPUT " DISPLAY PITTED CUTOFF HILL CONTOURS?(Y/N)->";ANS$
1870 IF ANS$-"N" THEN SYSTEM
1880 IF ANS$-"n" THEN SYSTEM •
1890 CLOSE fl
1900, LOCATE 14,20
1910 'Input the name of the plot file from program HCRIT.
1920 INPUT " INPUT NAME OF PLOTFILE FROM PROGRAM HCRIT";PLOT2$
1930 ON ERROR GOTO 2970
1940 OPEN PLOT2S FOR INPUT AS II
1950 ON ERROR GOTO 0
1960 'Make sure that the plot file was generated by program HCRIT.
1970 INPUTfl, PF2$
1980 IF PF2$-"HCRIT" THEN GOTO 2040
1990 CLS
2000 LOCATE 12,15
2010 PRINT PLOT2$ " IS NOT A FILE GENERATED BY PROGRAM HCRIT-TRY AGAIN"
2020 CLOSE fl
2030 GOTO 1900
2040 CLS
2050 'Check whether the hill identification number, hill name, number
2060 'of fitted contours, and contour identification numbers match.
2070 INPUTfl, IDH2,HNAME2$
2080 IF IDH20IDH1 THEN GOTO 3000
2090 IF HNAMB1$<>HNAME2$ THEN GOTO 3030
2100 INPUTfl, NC2
2110 IF NC10NC2 THEN GOTO 3060
2120 FOR J-l TO NC2
2130 INPUTfl, IDC2(J)
2140 IF IDC1(J)OIDC2(J) THEN GOTO 3090
2150 NEXT J
2160 'Input hill top elevation and contour elevations.
2170 INPUTfl, HTOP
2180 FOR J-l TO NC2
2190 INPUTfl, HCON(J)
2200 NEXT J
2210 'Input the number of critical elevations.
2220 INPUTfl, NCR
2230 'Begin loop for critical elevations.
2240 FOR 1-1 TO NCR
2250 'For each critical elevation, input the critical elevation, cutoff
2260 'hill centroid coordinates, orientation of the hill minor axis
2270 'with respect to the positive x-axis, and the inverse polynomial
2280 'fit parameters for each hill axis.
2290 INPUTfl, HC,XHTOPF,YHTOPF,ORENF,PA,PB>RLA,RLB
2300 'Determine the orientation of the major axis with respect to the
2310 'positive x-axis.
2320 ORENF-ORENF-90!
2330 CSE-COS(.017453*OR£NF)
2340 SNE-SIN(.017453*ORENF)
2350 'Retrieve the background plot of digitized contours(unedited or edited)
2360 PUT (0,0),IAR,PSET
2370 'Begin loop over contours.
2380 FOR J-l TO NC2
2390 'Contours must be at least one elevation unit higher than the
2400 'critical elevation if their elevations are to be used for the display
164
-------
2410 'of contours on the cutoff hill.
2420 IF HCON(J)<-HC-H! THEN GOTO 2770
2430 FLOG-LOG((HTOP-HC)/(HCON(J)-HC)-1!)
2440 AFIT-RIA*£XP((1!/PA)*FLOG)
2450 BFIT-RLB*£XP((1!/PB)*FLOG)
2460 'The aquation for the inverse polynomial contour is
2470 • (XP/AFIT)**PA+(YP/BFIT)**PB-1
2480 'in the coordinate system in which the x and y primed axes
2490 'coincide with the major and minor axes of the hill respectively.
2500 'Begin loop to calculate 800 contour point coordinates.
2510 FOR XXL TO 200
2520 IF L>99 GOTO 2580
2530 'Let x primed be the independent variable.
2540 XPOL-L*.01*AFIT
2550 YPOL-BFIT*(1!-(XPOI/AFIT)APA)A(1!/PB)
2560 GOTO 2610
2570 'Let y primed b« the independent variable.
2580 YPOL-(L-100)*.01*BFIT
2590 XPOL-AFIT*(11-(YPOL/BFIT)APB)A(1!/PA)
2600 'First quadrant(x primed-*,y primed—t-)
2610 XP-XPOL
2620 YP-YPOL
2630 GOSUB 3210
2640 'Second quadrant (x primed—t-,y primed—)--moving, clockwise
2650 XP-XPOL
2660 YP—YPOL
2670 GOSUB 3210
2680 'Third quadrant (x primed— ,y primed—)
2690 XP—XPOL
2700 YP—YPOL
2710 GOSUB 3210
2720 'Fourth quadrant (x primed— ,y primed—t-)
2730 XP—XPOL
2740 YP-YPOL
2750 GOSUB 3210
2760 NEXT L
2770 NEXT J
2780 XSHCF-SCRCX+(XHTOPF-XC)*DSCRXDDD
2790 YSHCF-SCRCY-(YHTOPF-YC)*DSCRYDDD
2800 XUL-XSHCF-1
2810 XLR-XSHCF+1
2820 YUL-YSHCF-1
2830 YLR-YSHCF+1
2840 'Plot a 3x3 box of points centered about the cutoff hill centroid.
2850 LINE (XUL,YUL)-(XLR,YLR),,BF
2860 PRINT HNAME2$ " ECRIT-" HC
2870 'Pause until the user strikes a key. If the ESC key is pressed, then
2880 'execution of the program is terminated.
2890 GOSUB 3160
2900 CLS '
2910 'End loop on critical elevations.
2920 NEXT I
2930 SYSTEM
2940 IF ERR-53 THEN PRINT "FITCON PLOT FILE NOT FOUNO-Press any key"
2950 GOSUB 3160
2960 SYSTEM
2970 IF ERR-53 THEN PRINT "HCRIT PLOT FILE NOT FOUND-Press any key"
2980 GOSUB 3160
2990 SYSTEM
3000 PRINT "FITCON AND HCRIT HILL IDs DO NOT MATCH-Press any key"
165
-------
3010 GOSUB 3160
3020 SYSTEM
3030 PRINT "FITCON AND HCRIT HILL NAMES DO NOT MATCH-Press any key"
3040 -GOSUB 3160
3050 SYSTEM
3060 PRINT "FITCON AND HCRIT NUMBER OF CONTOURS DO NOT MATCH-Press any key"
3070 GOSUB 3160
3080 SYSTEM
3090 PRINT "FITCON AND HCRIT CONTOUR IDs DO NOT MATCH-Press any key"
3100 GOSUB 3160
3110 SYSTEM
3120 END
3130 'Subroutine which causes program execution to pause until a key
3140 'is struck. If the ESC key is pressed, then program execution
3150 'will be terminated..
3160 A$-INKEYS: IF A$-"" THEN 3160
3170 IF AS-CHR$(27) THEN SYSTEM
3180 RETURN
3190 'Subroutine to rotate points into the x,y coordinate system before
3200 'plotting
3210 XFIT-XHTOPF+XP*CSE-YP*SNE
3220 YFIT-YHTOPF+XP*SNE+YP*CSE
3230 XS-SCROC+(XFIT-XC)*DSCRXDDD
3240 YS-SCRCY-(YFIT-YC)*DSCRYDDD
3250 PSET(XS,YS)
3260 RETURN
166
------- |