eddy/apogee_control
1
0
Fork 0
mirror of https://github.com/eddyem/apogee_control.git synced 2025-12-06 02:35:16 +03:00
Code Issues Packages Projects Releases Wiki Activity
apogee_control/airmass.c
Eddy 4776ef09d1 cloned
2015-02-10 00:08:55 +03:00

1158 lines
64 KiB
C
Raw Blame History

/*******************************************************************************
*
* airmass -- A program to compute the airmass for a given zenith angle.
*
* Copyright (C) 2000 Reed D. Meyer (rxm128@yahoo.com).
* You are permitted to copy and distribute this program as much as you want,
* as long as you leave the source code, including this comments section,
* intact and unmodified.
*
* To compile, type: gcc airmass.c -o airmass -s -O3 -lm
* (This command should work in any environment where the GNU C compiler,
* gcc, is installed, including Linux and most UNIX environments)
*
* To run, type:
*
* airmass [-a <alt>] [-d <day>] [-H <h>] [-l <lat>] [-L <wavelength>]
* [-p <P>] [-t <T>] [-v] <zenith>
*
* where:
* <zenith> is the APPARENT zenith angle (not the true angle in the
* absence of atmosphere) at which you want the airmass;
* -a <alt> sets the altitude of the observer to <alt> meters above mean
* sea level (the default is 0 meters);
* -d <day> specifies the date, in days since the beginning of the year
* (0 = midnight Jan. 1; the default is 80, corresponding roughly to
* the vernal equinox); don't worry about fractions of a day, since the
* seasonal effect on airmass is relatively small, as discussed below;
* -H <h>, if specified, tells the program to additionally compute an
* estimate of the water vapor column density and airmass, using a
* CRUDE model for the water vapor profile based on the local relative
* humidity at ground level, <h>. <h> is in percent (100 = saturated).
* The local humidity is NOT used in the computation of the standard
* (dry atmosphere) airmass. See EXTINCTION below for details.
* -l <lat> sets the observer's latitude to <lat> degrees (negative = south;
* default is 45 degrees);
* -L <wavelength> changes the observing wavelength in Angstroms (default
* 5500 A; don't worry too much about an exact figure here, as the
* wavelength dependence is very small, as discussed below);
* -p <P> specifies the local pressure in millibars, presumably read off
* a barometer at the time of observation (default is 1013.25 mb,
* the so-called standard atmosphere);
* -T <T> specifies the local temperature in Celsius, again at the time
* of observation (default is 15 C, again the standard atmosphere);
* -v is a flag that increase the "verbosity" of the output. In this
* context, the program computes the relative improvement it makes
* over some standard airmass formulas astronomers use (e.g. polynomial
* expansion in (sec(z) - 1) terms) -- see VERBOSITY OPTION below
* for details.
*
*
* VERSION HISTORY
* ===============
*
* 1.0 (2000 Aug. 13): first public release.
*
*
* BACKGROUND
* ==========
*
* This program, and the necessary underlying algorithm, was developed to
* find out to just what degree the standard approximate formulas for airmass
* are inaccurate. Many astronomers use the following polynomial in
* (sec(z) - 1) to calculate the airmass for a given zenith angle z:
*
* 2 3
* X = secz - 0.0018167(secz - 1) - 0.002875(secz - 1) - 0.0008083(secz - 1)
*
* (that is, those astronomers who do not use the simpler approximation
* X = sec(z), or the formula used in IRAF which is from Allen's _Astrophysical
* Quantities_.) The fact is, though, that some astronomers take this formula
* practically as "truth" and do not stop to think about its origin.
* This formula is a fit to a table of theoretical airmass as a function of z,
* computed at the turn of the century, when the structure of the atmosphere
* was less well-understood and equations were simplified to ease numerical
* integration since there wasn't the benefit of computers. To the best of
* my knowledge, the polynomial approximation above was first presented by
* R.H. Hardie in _Astronomical Techniques_ (W.A. Hiltner, ed.) (Chicago:
* U. Chicago Press), p. 180 (1962). This was a fit to computations by
* A. Bemporad (Mitt. d. Grossh. Sternwarte Heidelberg 4 (1904)), published
* again by E. Schoenberg in _Handbuch der Astrophysik_ (Berlin: J. Springer),
* 2, 268-273 (1929).*
* So there are two potential sources of inaccuracy in using this method:
* First, and perhaps the greatest source, is simply that the formula is
* a fitted function to a table of discrete values (and, moreover, was not
* designed to fit well the entire range of zenith angles). The second is
* that the underlying table of theoretical values itself is based on an
* old model atmosphere, and used simplifying assumptions which further
* degrade the accuracy of the model.
* What if, instead, we were to use the power of modern computers to
* obviate the need for most simplifying assumptions and, in addition,
* calculate the airmass directly for any given z rather than make use of
* interpolating or best-fit functions? This is what this program was
* designed to do. In carrying out such a "modernization" of the airmass
* function, we might as well incorporate more modern theories of the
* structure of the atmosphere, too.
*
* *Interestingly, the literature seems confused by Bemporad's work: some
* think the (sec(z) - 1) fit was to real DATA but in fact Bemporad just
* tabulated the results of his theoretical airmass formula. It is pleasing
* to note that even the old German references can be found in the Yale
* Astronomy Library.
*
*
* DERIVATION
* ==========
*
* We assume a spherically symmetric atmosphere of constant chemical
* composition. The equation for the column density, which can be derived
* strictly from Snell's Law with the above assumption, is (c.f.
* Schoenberg (1929)):
*
* /\ r_m
* | rho dr
* | -----------------------------
* N(z) = | 2 2 2 (1)
* | r_o mu_o sin z
* | sqrt(1 - -----------------)
* | 2 2
* \/ r_o r mu
*
* where the subscript o refers to values at the observer, and r is the
* distance from the center of the earth, mu is the air's index of refraction
* at r, rho is the air's density at r, r_m is the upper limit on the
* integration (taken as the top of the mesosphere; see below) and z
* is the apparent zenith angle for which we're calculating the airmass.
* Note that this is just the expression for the column density of air
* directly overhead, with a correction factor (sqrt(...) in the denominator)
* to account for the light's actual path through the atmosphere when the
* source is not actually at the zenith.
* To get the airmass as it's traditionally defined (i.e., normalized to
* 1 at the zenith), simply divide by the column density at z = 0:
*
* N(z)
* X(z) = ------ . (2)
* N(0)
*
* The Clausius-Mossotti equation (see, for example, B. Garfinkel, AJ 72, 235
* (1967)) gives the index of refraction in terms of the density to an
* extremely good approximation; rewriting the equation gives
*
* 2 3 + 4 c rho
* mu = ------------- (3)
* 3 - 2 c rho
*
* where c is a wavelength-varying constant depending only on the chemical
* composition of the gas. A formula for air's index of refraction at a
* given wavelength appearing in the CRC Handbook of Chemistry and Physics
* can be rewritten to yield c:
*
* 13.412 0.3777 -7 273.15 K
* c ~ [2875.66 + -------- + ----------] 10 R [-----------] (4)
* 2 -8 4 -16 101325 Pa
* L 10 L 10
*
* where L is the wavelength in Angstroms and R is the ideal gas constant.
* (c has units of volume/mass so R needs to be expressed in energy/mass/K).
*
* The main problem is finding an expression for rho(r). We do this by
* assuming a temperature profile and then combining the ideal gas law
* (which expresses density in terms of pressure and temperature) with the
* condition of hydrostatic equilibrium (which gives density in terms of
* pressure). This yields a unique density profile, given an assumed
* temperature profile. The problem is the temperature profile. This can
* only be determined uniquely by solving the equations of fluid dynamics
* for the atmosphere coupled with knowledge of the radiative processes in
* the atmosphere with appropriate boundary conditions, i.e., the incredibly
* complex problem which still makes accurate weather forecasting impossible.
* So we assume a temperature profile which is based on average measured
* values of temperatures at different heights in the atmosphere. As I will
* argue below under ASSUMPTIONS, I feel that the approximate temperature
* profile probably does not deviate strongly from the real profile at the
* observer under most conditions, and since this is the primary assumption
* made in calculating the airmass for dry air, said airmass probably doesn't
* differ strongly from the real value.
* For our temperature profile, we start by assuming the atmosphere is
* broken down into several layers within which the temperature changes
* linearly with altitude. In any given layer n, the temperature is
* T(Q) = T_n + beta_n (Q - Q_n), where beta_n is the slope of the temperature
* rise with altitude within the layer n, T_n is the temperature at the bottom
* of that layer, Q is the height (defined below), and Q_n is the height of
* the bottom of the layer. This is the methodology used in the _U.S. Standard
* Atmosphere, 1962_ (Washington, DC: U.S. Govt. Printing Office), wherein a
* "standard atmosphere" appropriate for the continental United States is
* defined (a "standard atmosphere" being useful for aircraft and so forth).
* We take their bottom eight temperature layers (corresponding roughly to
* the troposphere, stratosphere, and mesosphere; top of eighth layer is at
* 90 km above mean sea level). But in order to generalize their work to make
* it more appropriate for any latitude and any season, we change the values
* defining the lowest two layers (the troposphere and lower stratosphere).
* We let the boundary between those two layers (the tropopause; Q_1) vary as
* a function of time of year and latitude. The same thing goes for beta_0,
* the temperature falloff in the first layer. The bottom of the first layer
* (Q_0) is set to the height Q of the observer's location. The boundary
* condition "temperature at bottom of first layer" (T_0) is set to the
* actual local temperature at the observer. The condition T_1 equals the
* temperature at the top of the first layer, i.e., as determined from the
* above expression for T with n = 0 and Q = Q_1. Finally, beta_1 is set
* by forcing the temperature at the top of the second layer to match T_2
* (fixed at 216.65 K, the value in the U.S. Standard Atmosphere). Although
* the temperatures in the third layer and higher (altitudes > 20 km) *do*
* vary with season and latitude, the variations are less severe than in the
* troposphere and lower stratosphere, and fortunately for us the atmosphere
* has already greatly thinned out by that point so these layers contribute
* little to the airmass calculation anyway.
* To determine Q_1 and beta_0 (the two quantities allowed to vary with
* time and latitude, and which define the lower two layers along with T_0),
* I extracted a table of values from two plots appearing on pages 48-49 of
* _General Meteorology_, 3rd ed., by H.R. Byers (New York: McGraw-Hill)
* (1959). Byers provides a plot of the temperature distribution within the
* troposphere as a function of latitude for summer, and another for winter.
* Keep in mind that these are average values. (In particular, the location
* of the tropopause at any given time is far more complex than the smooth
* line drawn in Byers' plots, due to the movement of the jet streams.) I
* then fit a polynomial to the data to get a smooth function for beta_0 for
* any given latitude. I fit another polynomial to the data to get a function
* for "r1", the altitude of the troposphere (which is related to Q_1). The
* polynomials fit quite well and in any case are far more than adequate given
* the fact that Byers' plots are time averages. beta_0 needed only a
* 4th-order polynomial to fit the data but r1 required a 10th-order
* polynomial. You can find the values of the best-fit coefficients by
* looking at the declarations of betapoly[] and r1poly[] in the source code.
* (The units are K/m and km, respectively.) Then, to simulate the effect of
* seasonal changes, the odd coefficients in the polynomials sinusoidally vary
* with day of the year (period = 1 year; no change at day number 202,
* corresponding to the peak of summer).
* Now that we have values for the beta_n's, T_n's, and Q_n's (either
* from the _U.S. Standard Atmosphere_ or the method above for the bottom two
* layers), the simplicity of the temperature model leads to easily solvable
* equations which result in a closed-form expression for the density:
*
* /
* | T_n (1 + g_E /(R beta_n))
* | [------------------------]
* rho | T_n + beta_n (Q - Q_n) , beta_n <> 0
* ----- =_/ (5)
* rho_n \
* | g_E
* | exp[ ------- (Q_n - Q) ]
* | R T_n , beta_n = 0
* \
*
* where as before, subscripts n refer to the bottom of some temperature layer
* n (0 <= n < 8); g_E is the standard value of the acceleration of gravity
* at the earth's surface (9.80665 m/s^2); and all other symbols are as
* defined above.
* In order to use (5) in Equation (1), we need to write r in terms of Q.
* (The only other quantity in the integrand of (1) which varies with r is mu,
* and we already have an expression for that in terms of rho.) So, let us
* discuss the definition of Q. The _U.S. Standard Atmosphere_ defines the
* geopotential altitude (here referred to as "h") as the integral from 0 to r
* of (g/g_E) dr (compare to their equation I.2.5-(1)), where g is the gravity
* at some r. They then define their eight lowest temperature layers in terms
* of the geopotential altitude h. The advantage of changing variables to h
* is that it lets us forget the fact that the g field is varying, thus
* simplifying the calculations. To write the equations more simply, I
* further define the quantity Q by Q = h - (r_E^2 / r_msl), where r_E is the
* "standard" radius of the earth (i.e., the radius you would get from the
* formula g_E = G M_E / r_E^2, where M_E is the mass of the Earth, G is the
* gravitational constant, and g_E is the "standard" value of gravity as
* before); and r_msl is the actual radius of the earth at the latitude of
* the observer (i.e., the distance between the earth's center and a point
* at mean sea level at the latitude of the observer; "msl" means "mean sea
* level"). r_msl varies because the earth is an oblate spheroid; the reason
* for distinguishing r_E and r_msl is that Earth's gravity at the equator is
* slightly less than at the poles, because the earth's radius is larger at
* the equator. Note that Q is always negative.
* We ignore the fact that the gravity vector is not exactly equal to
* the gravity vector of a spherical, nonrotating Earth (see ASSUMPTIONS
* below). Combining our definitions for h and Q lets us write r in terms
* of Q:
*
* 2
* r_E
* r = - ----- (6)
* Q
*
* Similarly we can find a Q for any desired r; for example, Q_o (the lower
* limit on the integrand in Equation (8), and the Q of the observer) is
* related to the observer's altitude "a" above mean sea level by
* Q_o = -r_E^2 / (r_msl + a).
* We calculate r_msl from the analytic geometry of the oblate spheroid,
* and the standard definition of the geographic latitude. (Latitudes of
* actual places on the earth are always geographic latitudes.) The result is
*
* 4 4 4 2
* 2 b + (a - b ) cos phi
* r = ------------------------- (7)
* msl 2 2 2 2
* b + (a - b ) cos phi
*
* where a is the major axis of the spheroid (the Earth's equatorial radius),
* b is the spheroid's minor axis (the Earth's polar radius), and phi is the
* observer's latitude. (Incidentally, the geocentric latitude phi_c can be
* found from: tan phi_c = (b^2 / a^2) tan phi.)
* Plugging (6) into (1) gives
*
* /\ Q_m 2
* | rho r_E dQ
* | -------------------------------
* N(z) = | 2 2 2 (8)
* | 2 Q mu_o sin z
* | R sqrt(1 - -----------------)
* | 2 2
* \/ Q_o Q_o mu
*
* (Q_m being the Q at the top of the highest temperature layer; for this
* layer, h = 88.743 km. We take this as a suitable upper limit on the
* integration because the mass above this height is negligible; the air
* density at this altitude is 5 orders of magnitude less than at sea level.
* A practical reason for choosing this height for our upper limit is that
* at about this height (the top of the mesosphere), the molecular weight of
* the air begins to change appreciably. While we can safely assume that all
* the air below this altitude has a fixed proportion of all the components
* important for extinction at visible and near-infrared wavelengths (besides
* water vapor, ozone, and aerosols), this assumption is not a good one
* above this rough altitude.
* Plugging Equations (3), (4), and (5) into (8) and integrating gives
* N(z); then we set z = 0 to use the same formulas and integration routines
* to get N(0); Equation (2) then gives us the desired airmass at z.
*
*
* EXTINCTION
* ==========
*
* Of course, the whole point of calculating airmass, for astronomical
* purposes anyway, is to determine the atmospheric extinction. In addition
* to calculating the airmass, the program also computes the column density,
* because it is perhaps a more appropriate figure for determining extinction.
* After all, the airmass only tells you what the extinction is going to be
* at some zenith angle, RELATIVE to at the zenith; you still need to know
* the extinction at the zenith! The standard approximate formulas for
* airmass cannot provide this important bit of information, which is directly
* related to the column density.
* Several sources contribute to the atmospheric extinction at any given
* wavelength. For visible and infrared wavelengths, we can limit the
* discussion to four primary contributors; the others are small enough to
* be negligible. (Actually, there is an EXTREMELY important fifth component,
* namely suspended water droplets, commonly known as clouds; but in this
* discussion we'll assume that you're not trying to do high-precision
* photometric or spectroscopic work through clouds, so we'll ignore them.)
* The four components are: dry air, water vapor, ozone, and dust. "Dust"
* means any particulate matter suspended in the air (i.e., aerosols), and
* "dry air" signifies clean (dust-free) air with zero humidity, and having
* a constant chemical composition.*
* The total extinction, in magnitudes, is given approximately by:
*
* E = 2.5 log(e) (k_air N_air + k_O3 N_O3 + k_dust N_dust + k_H20 N_H20)
*
* where the various k's are the extinction coefficients, in cm^2/g, for each
* contribution; the N's are the column densities (in g/cm^2) of each
* contribution; and the subscripts denote the contributions: air = dry,
* clean air of constant (or negligible) ozone composition; O3 = ozone;
* dust = suspended particles of any sort; H20 = water vapor. The k's are
* wavelength-dependent, making E a function of wavelength. This program
* computes the column density N_air, and optionally (with the "-H" flag on
* the command line) a crude estimate for N_H20. (If your k's are in units of
* cm^2 rather than cm^2/g, convert the column densities to number densities
* by multiplying by A/m:
*
* N (number density, units 1/cm^2) = N (g/cm^2) * A/m
*
* where A is Avogadro's constant and m is the mean molecular weight of the
* constituent. For dry air, m = 28.9644 g/mol, so the number density equals
* 2.07915 * 10^22 * the column density in g/cm^2 for the dry air component.)
* k_air * N_air is approximately 0.2 at visual wavelengths, for dry air,
* observing at the zenith.
* Let us discuss the "-H" option briefly and how a crude estimate of
* the water vapor airmass and column density is calculated. It is impossible
* to calculate the water vapor values from a simple theory; it requires the
* same sort of full-blown analysis that makes the calculation of the
* temperature profile impossible. Thus, one needs either a direct
* measurement of the water vapor as a function of altitude (for example,
* through LIDAR), or one must make some approximation. Based on an
* admittedly brief look into real vapor profiles, I propose the following
* extremely simplistic model, which does have the desired feature of being
* constrained by the actually measured value of the humidity at the
* observatory:
*
* /
* | H_0 , Q_0 < Q < Q_1
* |
* | Q_2 - Q
* H = _/ H_0 [ ----------- ]
* \ Q_2 - Q_1 , Q_1 < Q < Q_2
* |
* |
* | 0 , otherwise
* \
*
* where the Q's are as defined in the section "DERIVATION", H is the relative
* humidity at some point Q, and the H_n's are the relative humidities at the
* Q_n's.
* The primary motivation behind the above linear model (constant H up
* to the troposphere; linear falloff through the lower stratosphere; zero
* humidity above that) is that, as I just said, it is constrained by the
* measured value of the humidity at ground level (H_0). In fact, it is
* quite strongly constrained, since H throughout the troposphere is assumed
* equal to this value. Believe it or not, even though this is a very lame
* model, I believe it does have some merit, for the following reason: looking
* at a couple typical vapor profiles, the relative humidity did not tend to
* vary THAT MUCH throughout the lower troposphere -- "THAT MUCH" meaning
* "deviations larger than a factor of two, over large regions of the
* troposphere, seemed rare". Furthermore, there is one great numerical
* property of the relative humidity which "limits the carnage" as far as
* keeping the REAL relative humidity from really deviating by orders of
* magnitude from the assumed value of H_0. Which is this: a relative humidity
* above 100% is forbidden. If the relative humidity exceeds 100%, it will
* not be by very much, and will probably not last for long, because the
* water vapor will try to condense as soon as it can (e.g., it spots a dust
* particle floating by and grabs onto it). So, for all intents and purposes,
* we are numerically limited to values between 0% and 100% everywhere in the
* atmosphere, which is still a huge range of possible humidities but at least
* it is not all the real numbers.
* I am interested in seeing how well the above crude formula for H
* holds up against real data (both real measured humidity profiles as well
* as water vapor extinction computed from accurate photometry at some
* wavelength where water vapor dominates). If you can supply any feedback
* on this, or if you know of a better humidity model which doesn't depend
* on complicated fluid dynamics equations, please do let me know.
* Anyway, once we have an assumed humidity profile, it is straightforward
* to calculate the water vapor column density and airmass. The definition
* of relative humidity gives the water vapor pressure, given the vapor
* pressure at saturation. The Clausius-Clapeyron equation gives the
* saturation vapor pressure as a function of temperature. For the Clausius-
* Clapeyron equation, see Equation 4.4 in A.S. Monin, _Weather Forecasting
* as a Problem in Physics_ (Cambridge, MA: MIT Press), p. 15 (1972); we
* rewrite Monin's formulation of the equation so that it gives us exactly
* what we want, the saturation vapor pressure as a function of temperature,
* relative to some saturation pressure at some standard temperature. These
* two conditions (the humidity profile and the Clausius-Clapeyron equation)
* let us solve for the vapor pressure. Then, we assume local thermodynamic
* equilibrium (i.e., that the local temperature of the water vapor equals the
* temperature of dry air at that point -- a good assumption) and that the
* water vapor pressure behaves independently of the other gases in the
* neighborhood (not completely true, but certainly a good enough assumption
* given the poorly-known humidity profile). Finally, we need only apply the
* ideal gas law (very close to valid for water vapor) and we have our density
* as a function of location in the atmosphere. We use the same algorithm
* as with dry air, from that point on, to compute the column density and
* hence the airmass of the water vapor at some zenith angle z.
*
* *In our airmass and column density calculations, ozone is included in the
* dry air component (with a constant fraction of the total air pressure
* assumed), showing up as a small contribution to the mean molecular weight.
* The mean molecular weight is the only way ozone has any impact in the
* dry air calculation. Because the ozone concentration does change
* significantly throughout the atmosphere (it peaks in the lower
* stratosphere, and anyone who follows the news knows that there are "ozone
* holes", i.e., ozone is also a function of latitude and longitude), its
* column density CANNOT be calculated with the current code. You should
* determine the ozone extinction contribution from a separate algorithm
* and then apply that contribution to the total extinction using the formula
* for E above. (Since ozone is a small contribution to the mean molecular
* weight, do not worry about its effect in the calculation of the dry
* air column density or airmass.)
*
*
* NOTE ON ASSUMPTIONS
* ===================
*
* Here is a rough analysis of the expected error in the calculation of the
* column density and airmass from using various assumptions. As far as I
* know, the following list contains all the assumptions this program makes
* use of. The list is sorted in order of their expected impact for most
* zenith angles (the worst assumption listed first).
* Note that this assumption list is only for the DRY AIR column density
* and airmass. It is not a list of the assumptions on the water vapor column
* density or airmass which the program optionally calculates (see EXTINCTION
* above for more on that). Nor does it deal with ozone or aerosol column
* densities or airmasses.
* Since the exact quantitative impacts of the assumptions depend on
* multiple factors, most particularly the zenith angle of interest, I only
* treat them in a qualitative sense, using adjectives like "moderate",
* "small", etc. to denote their expected relative importance. From most to
* least in terms of impact, the adjectives are "moderate", "small", "very
* small", "minimal", "insignificant", and "unimportant". Based on some
* VERY crude and quick calculations, "moderate" would correspond to a percent
* error (from using the approximation) of roughly 1% or less; "small" would
* correspond to something like 0.2% or less; "very small" smaller still;
* "minimal" meaning something like 0.01% if not considerably less;
* "insignificant" here meaning "probably falls into the class 'unimportant',
* but I'm not sure, and in any case is probably not a bigger assumption than
* something in the 'minimal' class"; and "unimportant" meaning around 0.001%
* if not considerably less.
* One thing to keep in mind is that these descriptions really only apply
* to LOW TO MODERATE ZENITH ANGLES (i.e., z less than ~ 60 degrees). Many
* of the assumptions tend to become worse with z; those that are believed to
* have a significant z dependence (so that the assumptions are likely to be
* notably worse at extreme zenith angles) have the remark "z dependence."
* For example, the temperature profile assumption might become quite severe
* (i.e., larger than the rough 1% upper limit corresponding to the so-called
* "moderate" class) for zenith angles near 90 degrees.
* Another thing to keep in mind is that under no circumstances are the
* computed column density or airmass of higher accuracy than about 1 part in
* 10^8. This is because we consider the numerical integration "converged"
* when the column density changes by less than 1 part in 10^8 between
* iterations. One could always change this tolerance limit, but there seems
* little point when certain assumptions like the temperature profile likely
* contribute a far larger fractional error than 1 in 10^8. (The convergence
* criterion is the motivation, by the way, for printing the results to 8
* significant digits. If all the following assumptions were negligible, in
* principle all 8 digits would be valid.)
*
* TEMPERATURE PROFILE: SMALL to MODERATE for column density, SMALL for
* airmass. z dependence (possibly strong). This is almost certainly
* the biggest assumption for the dry air calculations. The problem is
* that we cannot compute this profile theoretically, and it would be
* impractical for most observatories to launch radiosondes or whatever
* to probe the atmosphere overhead. Fortunately, Nature is kind to us
* in that the day-to-day fluctuations in the temperature at a given
* altitude are generally small compared to the absolute temperature;
* that is, the percent error of the fluctuations is relatively small
* when the temperature is expressed in Kelvins. According to
* Figure II.2.1(b) in the _U.S. Standard Atmosphere_, at any given
* altitude the temperature will tend to be within about 10 K of the
* average value, over much of the altitude range of interest. Thus
* the typical percent error is perhaps 5% at any given altitude and
* time under most conditions. The variance is greatest at ground level,
* but we correct for that to a great extent by *incorporating* the
* ground temperature explicitly into our calculations of the temperature
* profile for the bottom two layers. This should eliminate much of the
* variance at the lowest altitudes. Furthermore, systematic offsets
* in one direction at one range of altitudes get balanced to some
* extent by offsets at another range, so the overall effect is lessened.
* And, unlike the U.S. Standard Atmosphere, we incorporate seasonal and
* latitude trends into our calculation of the bottom two layers, which
* helps reduce the offset between the theoretical profile and the real
* one. Based on a crude experiment, I expect the error from not knowing
* the true profile to affect the column density by < 0.8%, and the
* airmass by < 0.3%. These are conservative estimates, and the real
* error might be considerably less.
* SPHERICAL SYMMETRY: SMALL under most conditions. z dependence (possibly
* strong). Except for the calculation of r_msl (see the discussion of
* the next assumption), complete spherical symmetry is assumed. This,
* of course, is never absolutely true, for otherwise there would not
* be, for example, high-pressure and low-pressure weather systems since
* all locations on earth would have the same weather conditions at a
* given altitude. The question really is to what extent these deviations
* from spherical symmetry affect the true airmass. These deviations
* do not affect the column density at z=0 (the zenith) but in principle
* would for any other z. Weather systems are essentially confined to
* the troposphere and lower stratosphere, i.e., roughly the lowest 20 km.
* Thus, except for extreme zenith angles, large weather patterns would
* not really affect our airmass that much because the horizontal distance
* starlight would travel through the weather system would not be great
* (less than 20 km for z = 45 degrees). So, only local deviations from
* symmetry are important. Small-scale effects such as turbulence would
* average out over the light's path, so only mesoscale deviations could
* be important, i.e., regions with sizes on the order of 20 km * sec(z)
* or less, but larger than ~ 1 km. Furthermore, the dominant impact of
* these varying mesoscale regions would be a varying temperature, but
* this should be on the same order of magnitude if not considerably
* smaller than the vertical temperature variations, and therefore the
* effect of this assumption would be about the same if not much smaller
* than the "temperature profile" assumption. One special case of
* deviation from spherical symmetry is discussed under the next
* assumption heading.
* GRAVITY VECTOR: SMALL. As the _U.S. Standard Atmosphere_ points out, the
* true gravity vector at any point does not generally point straight
* down to the center of the earth nor have the magnitude that you would
* calculate from Newton's Law of Gravitation for a sphere. There are
* two main reasons. The somewhat more important one, numerically, is the
* centrifugal force caused by the rotation of the earth. Additionally,
* the earth is really an oblate spheroid (ironically, the earth is
* oblate precisely BECAUSE of the centrifugal force). In principle, we
* could account for these effects in the computer code, since the theory
* is quite well known. However, there really is no reason to complicate
* the algorithm in this manner, because the effects are rather small.
* Compared to the formula for a gravitating sphere, the first term of
* the centrifugal correction is ~ 580 times smaller, and the first term
* of the oblate spheroid correction is ~ 920 times smaller. These are
* both on the order of the error of the "temperature profile" assumption
* above. Thus, it did not seem necessary to complicate the code unless
* we have better information on the temperature profile. The oblateness
* of the earth *is* accounted for in one way, however, since it is
* easy to implement it in the code. It is in the calculation of r_msl,
* the true distance between mean sea level on some point on the earth's
* surface and the earth's center. (See Equation (7) under DERIVATION
* above.)
* CONSTANCY OF CHEMICAL COMPOSITION: VERY SMALL. z dependence. This is the
* assumption that all the components of dry, clean air have the same
* partial pressures throughout the atmosphere region over which we
* integrate (i.e., throughout the lowest 90 km). For all the atmospheric
* gases which are important to extinction in visible and near-infrared
* wavelengths, except water vapor and ozone, this is a very good
* approximation. It is such a good approximation that the U.S. Standard
* Atmosphere assumes a constant molecular weight of 28.9644 throughout
* the lowest 90 km. Water vapor is considered separately (it is not
* a component of "dry, clean air", of course); ozone is a special case
* in that it is included into the dry air calculations through its
* molecular weight, but is considered constant throughout the atmosphere
* even though it's not. (Ozone's contribution to the total molecular
* weight is small in any case, no matter how we consider its contribution
* to dry air.) See also the "airmass independent of contaminants"
* assumption below.
* IDEAL GAS: MINIMAL. In many sources in the literature, dry air is said
* to behave very similarly to a ideal gas. The deviation from an ideal
* gas is without question smaller than some of the assumptions above.
* HYDROSTATIC EQUILIBRIUM: MINIMAL. For hydrostatic equilibrium to be
* seriously violated, a parcel of air would practically have to be in
* a state of free fall. Small-scale regions of turbulence would average
* out over a light's path, so a large region would have to be undergoing
* this type of motion extreme. This seems quite unlikely except maybe
* near a tornado. In which case, you wouldn't be observing. Thus,
* departure from hydrostatic equilibrium is expected to have a quite
* minor effect on the true airmass.
* AIRMASS INDEPENDENT OF CONTAMINANTS: INSIGNIFICANT. z dependence.
* "Contaminants" here means water vapor, ozone, or aerosols (where
* ozone abundance varies with height). As stated in the section
* "IMPORTANCE OF ATMOSPHERIC VARIABLES", water vapor is a very weak
* contribution to the total airMASS. The same goes for ozone and
* aerosols. We treat those components separately when considering
* EXTINCTION (q.v.), since water vapor, ozone, and aerosols do contribute
* heavily to the total extinction, even though they contribute so little
* to the mass of dry air.
* CLAUSIUS-MOSSOTTI EQUATION: INSIGNIFICANT. z dependence. As it is, the
* Clausius-Mossotti Equation describes pretty well the relation between
* index of refraction and density for many materials, including air (c.f.
* Garfinkel 1967). Furthermore, the index of refraction for air anywhere
* within the earth's atmosphere is going to be pretty close to 1, no
* matter what the functional form of mu(rho) is. Finally, mu has only
* a rather weak impact on the airmass calculation (mainly because mu
* does stay close to 1); it is more important in calculating refraction.
* All these facts combined together illustrate that any errors which
* arise from using the Clausius-Mossotti Equation to calculate mu are
* probably quite small.
* ALL MASS BELOW 90 KM: UNIMPORTANT. The density at 90 km is about 10^5
* times less than at sea level, and continues to fall off nearly
* exponentially above that (see e.g. Figure I.2.11 in the _U.S. Standard
* Atmosphere_). Thus, the contribution to the column density, and hence
* airmass, is about 10^5 times smaller than the contribution from the
* lower layers, so that this assumption contributes an error of around
* 0.001%. If you desire more precision, you could change the upper
* limit on the integration in the code, but to do it properly you would
* have to deal with the physics of the thermosphere, which does not
* have quite the same properties as the lower layers (for example, the
* assumption "constancy of chemical composition" breaks down).
* CRC FORMULA FOR MU(WAVELENGTH): UNIMPORTANT. This means the error in
* calculating the index of refraction mu for a given wavelength, because
* the formula in the CRC Handbook is not exact. The wavelength of
* observation already has almost no effect on the airmass calculation
* (see IMPORTANCE OF ATMOSPHERIC VARIABLES below); the error in using
* the CRC formula (which itself is not a bad estimator for the index of
* refraction, at least for reasonable temperatures and pressures) must
* therefore have even less of an impact. I feel pretty confident in
* saying that this is a good assumption. Note that the code does not
* include an additional term in the formula for mu which accounts for
* humidity; this is safe because of the smallness of the term and
* because airmass depends on mu so weakly.
*
*
* VERBOSITY OPTION
* ================
*
* With the -v ("verbose") flag, the program also spits out the airmass
* estimated by the polynomial approximation formula, and the percent error
* in this approximation (relative to the supposed "true" airmass computed
* by this program). It also computes the result and error of the simple
* formula X ~ sec(z), and of the airmass formula used within IRAF. The
* latter formula comes from J.A. Ball, _Algorithms for the HP-45 and HP-35_
* (1975 ed.) (Cambridge, MA: Center for Astrophysics), with Q assumed to be
* 750. (Ball's algorithm comes from C.W. Allen, _Astrophysical Quantities_
* (3rd ed.) (London: Athlone), p. 133 (1973).)
* For all but extreme zenith angles (i.e., for z < 86 degrees), the
* sec(z) - 1 polynomial is the best of the three approximations. The
* polynomial only goes bad at extreme zenith angles because it was not
* designed to fit these large angles. Above about 86 degrees, the IRAF
* formula is a better approximation. (The sec(z) - 1 polynomial goes negative
* above z ~ 88.35, which is of course physically impossible.) But at these
* extremes it is clearly best to use the airmass algorithm utilized in this
* program.
* For low to moderate zenith angles, i.e., z < 70 roughly, it seems
* perfectly safe to use the sec(z) - 1 polynomial to approximate the airmass
* unless very high accuracy is required or atmospheric conditions are highly
* abnormal. The error of that approximation is about 0.09% for typical
* atmospheric conditions. The sec(z) - 1 polynomial error exceeds 1% at
* about z = 85 degrees, again for typical conditions. The approximation
* rapidly deteriorates above this z.
* To judge the conditions over which the sec(z) - 1 polynomial is
* suitable, here is a table for a few large zenith angles. Again, this is
* for the default values of atmospheric variables. Also tabulated are
* percent errors for the polynomial when atmospheric conditions are at the
* most extreme one could expect to encounter at any reasonable point on the
* earth's surface. The sec(z) and IRAF formula errors are also included.
*
* PERCENT ERRORS RELATIVE TO THIS PROGRAM
* z (deg) sec(z)-1 (extreme) IRAF (extreme) sec(z) (extreme)
* ======= =================== ================ =================
* 60 0.0346 0.12 0.111 0.20 0.310 0.40
* 65 0.0537 0.18 0.169 0.29 0.474 0.60
* 70 0.0889 0.27 0.273 0.46 0.774 0.96
* 75 0.159 0.44 0.489 0.77 1.41 1.7
* 80 0.280 0.90 1.04 1.7 3.16 3.8
* 81 0.294 1.0 1.25 2.0 3.87 4.6
* 82 0.279 1.2 1.52 2.4 4.84 5.8
* 83 0.186 1.7 1.88 3.0 6.20 7.4
* 84 0.106 2.5 2.37 3.9 8.20 9.8
* 85 0.942 4.1 3.02 5.0 11.3 14.
* 86 3.42 7.6 3.89 6.6 16.5 20.
* 87 12.0 17. 4.99 8.9 26.2 31.
* 88 52.0 56. 6.06 12. 47.7 56.
* 88.1 61.8 65. 6.13 12. 51.2 60.
* 88.2 73.9 76. 6.20 13. 55.2 64.
* 88.3 89.2 90. 6.25 13. 59.7 70.
* 88.4 ---- ---- 6.28 13. 64.8 76.
* 88.5 ---- ---- 6.29 14. 70.6 82.
* 88.6 ---- ---- 6.28 14. 77.3 90.
* 88.7 ---- ---- 6.25 14. 85.1 99.
* 88.8 ---- ---- 6.19 14. 94.3 110
* 88.9 ---- ---- 6.09 15. 105. 120
* 89 ---- ---- 5.96 15. 118. 140
* 89.5 ---- ---- 4.58 16. 266. 300
* 89.9 ---- ---- 2.30 16. 1.46E3 1.7E3
* 89.99 ---- ---- 1.61 16. 1.50E4 1.7E4
* 89.999 ---- ---- 1.53 16. 1.50E5 1.7E5
*
* As the above table shows, atmospheric variable effects become important at
* large zenith angles (on the order of z = 80 and higher, but perhaps z ~ 70
* or even lower angles if you need high accuracy or the deviation of the
* atmospheric variables from standard values is extreme).
* Incidentally, sec(z) *ALWAYS* overestimates the true airmass; that is,
* you could consider sec(z) an upper bound on the value of the true airmass
* if you needed a rough estimate. (And, of course, the sec(z) estimate
* progressively gets worse with larger z.)
*
*
* IMPORTANCE OF ATMOSPHERIC VARIABLES
* ===================================
*
* For low to moderate zenith angles, the dependence of airmass on atmospheric
* variables (pressure, temperature, etc.) is generally small unless you need
* high accuracy or your observing conditions are extreme. For these zenith
* angles, you can probably safely assume the values built into the program by
* default: temperature = 15 C, pressure = 1013.25 millibars, altitude = mean
* sea level, latitude = 45 degrees, day of year = 80, wavelength = 5500
* Angstroms. (Most of these aren't technically atmospheric conditions but
* are considered here as site-dependent variables like the local temperature
* and pressure.) These values correspond in a rough way to the definition of
* the U.S. Standard Atmosphere. (The correspondence is not exact because of
* the way we calculate the characteristics of the bottom two temperature
* layers.)
* The importance of these variables in affecting the airmass calculation
* is roughly: temperature most important; pressure of medium importance;
* seasonal effects (latitude combined with day of year) not very important;
* altitude and wavelength unimportant. A note about the above sentence:
* it only refers to the effects of actually changing the PARAMETERS to the
* program; each parameter is associated with an effect, e.g., "-a" for
* altitude. However, the "real-life" effect of changing certain parameters
* is strongly coupled to changing others. That is, if the observer literally
* changes his *altitude* (for example, he moves further up a mountain), the
* locally observed *pressure* is going to change tremendously (after all,
* altimeters are in essence just barometers), as well as the observed
* temperature (it tends to get colder as you move up), so that the true
* effect of changing the observer's altitude is actually stronger than
* if the local temperature or pressure alone would change, despite what
* the first sentence says. Similarly, if the observer moves a great distance
* in latitude, the local temperature will obviously change. To make sure
* there isn't any confusion about this point: just feed in the actual values
* of the observing conditions (temperature, pressure, altitude, latitude,
* day of year, wavelength) and everything will be fine.
* The wavelength (the wavelength of light we're observing at) has such
* a weak effect because it only affects the airmass calculation indirectly,
* in the calculation of the path the light takes through the atmosphere; the
* majority of the wavelength-dependent effect enters by way of refraction,
* whose effect is to change the APPARENT zenith angle of the object, which
* of course by far is the main parameter (z) affecting airmass; non-refraction
* wavelength dependence is thus rather minor. Altitude is minor because of
* the way the bottom two temperature layers are calculated. (It only enters
* the algorithm at all through the calculation of the lower limit of the
* integration.) Temperature is the most important observing condition
* because it determines the temperature profile within the two lowest layers,
* which is of tremendous importance in calculating the density (c.f.
* Equation (5).)
* Humidity is NOT treated in this code for determining the airmass
* of dry air. This is because, while certainly the presence of water vapor
* (especially condensed water!) affects the extinction coefficient of air,
* it does not affect the airMASS very much. That is, it does not greatly
* affect the molecular weight or total density, since even saturated air does
* not carry much water vapor by mass. Furthermore, to calculate the dry
* airmass precisely, taking the presence of water vapor correctly into
* account, requires an accurate knowledge of humidity as a function of height
* in the atmosphere, and as discussed under EXTINCTION above, we don't
* usually have such information. In any case humidity's affect on airmass
* (particularly since airmass is normalized to 1 at the zenith) has got to be
* very small, for any conceivable atmospheric condition besides perhaps,
* perhaps, the innermost winds of a hurricane, and at that point other
* assumptions like spherical symmetry break down anyway. Humidity is
* implicitly included in the sense that refraction changes the apparent zenith
* angle of the object, and any good refraction algorithm includes humidity.
*
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*******************************************************************************
* Mathematical and physical constants, and constants derived from them.
* It's probably best to keep these as they are, since the temperature levels
* defined in the U.S. Standard Atmosphere were calculated using the values
* of the physical constants below.
******************************************************************************/
#define DEG2RAD (M_PI/180.) /* factor for converting degrees to radians */
#define EQUATOR 6378178. /* Earth's equatorial radius in meters */
#define POLAR (EQUATOR*(1.-1./298.32)) /* Earth's polar radius in meters */
#define GEOCENTGRAVCONST 3.9862216E14 /* G * mass of Earth (m^3 / s^2) */
#define STANDARDG0 9.80665 /* standard value of g0, in m/s^2 */
#define GASCONSTANT 8.314510 /* universal gas constant (R), in J/mol/K */
#define AIRKGPERMOLE 0.0289644 /* dry air mol.wt. 28.9644 * 0.001 g -> kg */
#define VAPORKGPERMOLE 0.0180153 /* water mol.wt. 18.0153 * 0.001 g -> kg */
#define AIRCONSTANT (GASCONSTANT/AIRKGPERMOLE) /* R for dry air, J/kg/K */
#define VAPORCONSTANT (GASCONSTANT/VAPORKGPERMOLE) /* R for vapor, J/kg/K */
#define CONSTC (AIRCONSTANT*273.15/101325.) /* constant for calculating c */
#define CONSTRHO (-(STANDARDG0)/AIRCONSTANT) /* constant for calculating rho */
#define STANDARDR02 (GEOCENTGRAVCONST/STANDARDG0) /* GM/g0, in m^2 */
#define CPVAPOR 1.81E3 /* typical heat capacity of water vapor, J/kg/K */
#define CWATER 4.20E3 /* typical specific heat of water, J/kg/K */
#define T0VAPOR 298.15 /* saturation vapor press. calibrated to this T (K) */
#define LATENTVAPOR 2443240. /* latent heat vap. at 25 C and 760 mm Hg (J/kg) */
#define P0VAPOR 3167.6 /* saturation vapor press. at 25 C (Pascals) */
#define VAPPRESSFACT1 ((CWATER-CPVAPOR)/VAPORCONSTANT) /* used in vaporrho() */
#define VAPPRESSFACT2 (((LATENTVAPOR/T0VAPOR)+(CWATER-CPVAPOR))/VAPORCONSTANT)
#define NTLEVELS 8 /* Number of different temperature regions in atmosphere */
#define K 5 /* constant used in polint() and qromb() */
#define IMAX 50 /* constant used in qromb() */
/*******************************************************************************
* Global variables (actually, constants, once they are computed in main().)
* The first two default values for betavec and tvec are replaced by values
* computed in main() (the default values come from the U.S. Standard
* Atmosphere, but we attempt a more "sophisticated" estimation of them based
* on the observed local temperature, pressure, etc.)
******************************************************************************/
static double cee, const6, const7, rhovec[NTLEVELS], bigrvec[(NTLEVELS+1)],
betavec[NTLEVELS]={-0.0065, 0., 0.0010, 0.0028, 0., -0.0020, -0.0040, 0.},
tvec[NTLEVELS]={288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 252.65,
180.65}, relhumid=-1.;
/*******************************************************************************
* vaporrho - called only if "-H" was passed on the command line.
* Calculates the wator vapor density, given the temperature (t) and
* relative humidity (rh).
******************************************************************************/
inline double vaporrho(double t, double rh)
{
return rh*(P0VAPOR/VAPORCONSTANT)*pow(T0VAPOR/t, VAPPRESSFACT1)*
exp(VAPPRESSFACT2-(VAPPRESSFACT2*T0VAPOR)/t)/t;
}
inline double frho(double rdiff, double rhobottom, double tbottom, double beta)
{
return (fabs(beta)<1.E-10) ? rhobottom*exp(CONSTRHO*rdiff/tbottom) :
rhobottom*pow(tbottom/(tbottom+beta*rdiff), 1.-CONSTRHO/beta);
}
inline double musquare(double rho)
{
return (3.+4.*cee*rho)/(3.-2.*cee*rho);
}
double columndensityint(double r)
{
unsigned long ju=NTLEVELS, jm, index=0;
double rho;
/* figure out index, the temperature layer corresponding to the given r. */
/* The algorithm is based on "locate" in Numerical Recipes, and handles */
/* r values outside the "proper" range (we simply assume the bottom (top) */
/* temp. layer continues on to negative (positive) infinity, in such cases)*/
while (ju-index>1) {
jm=(ju+index)>>1;
if (r>bigrvec[jm]) index=jm;
else ju=jm;
}
rho=frho(r-bigrvec[index], rhovec[index], tvec[index], betavec[index]);
r*=r;
return STANDARDR02*rho/(r*sqrt(1.-const6*r/musquare(rho)));
}
double vaporcolumndensityint(double r)
{
double rho, rdiff, vaprho;
if (r<bigrvec[1]) { /* if we are in the troposphere */
rdiff=r-*bigrvec;
rho=frho(rdiff, *rhovec, *tvec, *betavec);
vaprho=vaporrho(*tvec+*betavec*rdiff, relhumid);
} else { /* assume we're in the lower stratosphere */
rdiff=r-bigrvec[1];
rho=frho(rdiff, rhovec[1], tvec[1], betavec[1]);
vaprho=vaporrho(tvec[1]+betavec[1]*rdiff, const7*(bigrvec[2]-r));
}
r*=r;
return STANDARDR02*vaprho/(r*sqrt(1.-const6*r/musquare(rho)));
}
/*******************************************************************************
* trapzd - auxiliary to qromb integration routine
* Adapted from J. Percival, who in turn adapted it from the original
* Numerical Recipes routine by Press et al.
******************************************************************************/
inline double trapzd(double (*func)(double), double a, double b, int n)
{
double del, sum, tnm, x;
int j;
static double s;
static int it;
if (n<=0) {
s=0.5*(b-a)*((*func)(a)+(*func)(b));
it=1;
} else {
tnm=it;
del=(b-a)/tnm;
x=a+0.5*del;
sum=0;
for(j=0;j<it;j++) {
sum+=(*func)(x);
x+=del;
}
s=0.5*(s+(b-a)*sum/tnm);
it*=2;
}
return s;
}
/*******************************************************************************
* polint - auxiliary to qromb integration routine
* Adapted from J. Percival, who in turn adapted it from the original
* Numerical Recipes routine by Press et al.
******************************************************************************/
inline double polint(double *xa, double *ya)
{
double c[K], d[K], den, dif, dift, ho, hp, w, y;
int i, m, ns;
/* find the index of the closest table entry */
ns=0;
dif=fabs(xa[ns]);
for(i=0;i<K;i++) {
dift=fabs(xa[i]);
if (dift<dif) {
ns=i;
dif=dift;
}
c[i]=ya[i];
d[i]=ya[i];
}
y=ya[ns--]; /* first guess */
for (m=0;m<(K-1);m++) {
for(i=0;i<(K-1)-m;i++) {
ho=xa[i];
hp=xa[i+m+1];
w=c[i+1]-d[i];
den=ho-hp;
den=w/den;
c[i]=ho*den;
d[i]=hp*den;
}
y+=(2*ns+1<(K-2)-m) ? c[ns+1] : d[ns--];
}
return y;
}
/*******************************************************************************
* qromb - Romberg rule integration
* Adapted from J. Percival, who in turn adapted it from the original
* Numerical Recipes routine by Press et al.
* Changed the convergence criterion to fabs(ss-oldss)<eps.
******************************************************************************/
double qromb(double (*func)(double), double a, double b, double eps)
{
double ss=0., oldss=0., h[(IMAX+1)], s[IMAX];
int i;
h[0]=1;
for (i=0;i<IMAX;i++) {
s[i]=trapzd(func, a, b, i);
if (i>=(K-1)) {
ss=polint(h+i-(K-1), s+i-(K-1));
if (fabs((ss-oldss)/ss)<eps) return ss;
oldss=ss;
} else oldss=s[i];
h[i+1]=0.25*h[i];
}
printf("EXCEEDED MAX ITERATIONS (%d) IN qromb\n", IMAX);
return ss;
}
void usage(void)
{
printf("Usage:\n\nairmass [-a <alt>] [-d <day>] [-H <h>] [-l <lat>] [-L "
"<wavelength>] [-p <P>]\n [-t <T>] [-v] <zenith>\n\nwhere "
"<zenith> is the APPARENT zenith angle, in degrees, whose airmass you\n"
"want to calculate. Program by default uses standard atmosphere "
"values;\nto override these use any of: \n\t"
"-a <alt> : observer's altitude above sea level, in meters (default 0)"
"\n\t-d <day> : day number from beginning of year (0 = Jan 1; default "
"80)\n\t-H <h> : relative humidity at observer, in percent (default "
"n/a)\n\t-l <lat> : observer's latitude in degrees (default +45)"
"\n\t-L <wavelength> : the observing wavelength in Angstroms (default "
"5500)\n\t-p <P> : local pressure, in millibars (default 1013.25)"
"\n\t-t <T> : local temperature, in Celsius (default 15).\n"
"Additionally, you can specify the -v option to calculate the airmass"
"\npredicted by standard approximate formulas (sec(z)-1 polynomial, "
"IRAF, sec(z))\nand their error relative to the airmass computed by "
"this program.\nIf -H is used, program will also compute a *CRUDE* "
"estimate of the water vapor\ncolumn density and airmass, using the "
"specified <h>.\n");
exit(1);
}
int main(int argc, char **argv)
{
char *s, verbos=0;
int i;
double colz, col0, cosyearfrac, sinz, rmsl, bigrmsl, zfact, approximation,
rjunk, coslatsq, airmass, t0=288.15, p0=101325., lambda=5500., alt=0.,
z=-1., lat=45., daynum=80., betapoly[5]={-0.0065107, -4.5403e-06,
3.6599e-07, -2.2174e-09, 7.9392e-12}, r1poly[11]={17.204, 8.9155e-03,
-3.6420e-03, 2.5617e-05, 2.4796e-07, -1.2774e-08, 1.3017e-10,
2.0151e-12, -2.6985e-14, -1.0397e-16, 1.4849e-18}, hvec[(NTLEVELS+1)]=
{0., 11000., 20000., 32000., 47000., 52000., 61000., 79000., 88743.};
/* get command line arguments */
argc--;
while(argc--) {
if( **++argv == '-') for(s = *argv+1;*s != '\0';s++) switch(*s) {
case 'a': /* observer's altitude in meters above mean sea level */
sscanf(*++argv,"%lg",&alt);
argc--;
break;
case 'd': /* Day number from beginning of year, 0 = midnight Jan 1*/
sscanf(*++argv,"%lg",&daynum);
argc--;
break;
case 'H': /* Do water vapor calc. Relative humidity in percent */
sscanf(*++argv,"%lg",&relhumid);
relhumid*=0.01; /* convert to fraction (range 0 - 1) */
argc--;
break;
case 'l': /* observer's latitude in degrees north of Equator */
sscanf(*++argv,"%lg",&lat);
argc--;
break;
case 'L': /* observing wavelength in Angstroms */
sscanf(*++argv,"%lg",&lambda);
argc--;
break;
case 'p': /* local pressure in millibars */
sscanf(*++argv,"%lg",&p0);
p0*=100.; /* convert to Pascals */
argc--;
break;
case 't': /* local temperature in Celsius */
sscanf(*++argv,"%lg",&t0);
t0+=273.15; /* convert to Kelvins */
argc--;
break;
case 'v': /* increase verbosity: compute error in approximation */
verbos++;
break;
default:
printf("illegal option: %s\n",s);
usage();
} else sscanf(*argv,"%lg",&z);
}
if ((z<0.) || (z>= 90.)) {
printf("zenith angle not found or unacceptable value\n");
usage();
}
sinz=sin(z*DEG2RAD);
cee=CONSTC*(2.87566E-4+134.12/(lambda*lambda)+3.777E8*pow(lambda, -4.));
cosyearfrac=cos((daynum-202.)*(360.*DEG2RAD/365.));
coslatsq=cos(lat*DEG2RAD);
coslatsq*=coslatsq;
rmsl=sqrt(((POLAR*POLAR*POLAR*POLAR)+(EQUATOR*EQUATOR*EQUATOR*EQUATOR-
POLAR*POLAR*POLAR*POLAR)*coslatsq)/((POLAR*POLAR)+(EQUATOR*EQUATOR-
POLAR*POLAR)*coslatsq));
bigrmsl=(-STANDARDR02)/rmsl;
/* Calculate bigrvec, the bigr at the bottom of each temperature layer */
*bigrvec=(-STANDARDR02)/(rmsl+alt);
rjunk=r1poly[10];
for (i=9;i>=0;i--) rjunk=rjunk*lat+((i%2) ? r1poly[i]*cosyearfrac :
r1poly[i]);
bigrvec[1]=(-STANDARDR02)/(rmsl+rjunk*1000.);
for (i=2;i<(NTLEVELS+1);i++) bigrvec[i]=hvec[i]+bigrmsl;
/* Set up our temperature profile for the troposphere/lower stratosphere */
*betavec=betapoly[4];
for (i=3;i>=0;i--) *betavec=*betavec*lat+((i%2) ? betapoly[i]*cosyearfrac :
betapoly[i]);
*betavec*=(-rmsl/bigrvec[1]);
*tvec=t0;
tvec[1]=t0+*betavec*(bigrvec[1]-*bigrvec);
betavec[1]=(tvec[2]-tvec[1])/(bigrvec[2]-bigrvec[1]);
/* Compute the density at the bottom of each temperature layer */
*rhovec=p0/(AIRCONSTANT*t0);
for (i=0;i<(NTLEVELS-1);i++) rhovec[i+1]=frho(bigrvec[i+1]-bigrvec[i],
rhovec[i], tvec[i], betavec[i]);
const6=musquare(*rhovec)*sinz*sinz/(*bigrvec*(*bigrvec)); /* for z */
colz=qromb(columndensityint, *bigrvec, bigrvec[NTLEVELS], 1.E-8);
const6=0.; /* equivalent to setting z = 0. */
col0=qromb(columndensityint, *bigrvec, bigrvec[NTLEVELS], 1.E-8);
airmass=colz/col0;
printf("for z=%.15g: column density=%.8g g/cm^2; AIRMASS=%.8g\n", z,
0.1*colz, airmass);
/* if desired, compute error in using various airmass approximations */
if (verbos) {
printf("Comparisons to various approximate formulae for the airmass:\n");
/* Do the sec(z)-1 polynomial approximation */
zfact=1./cos(z*DEG2RAD)-1.;
approximation=1.+zfact*(0.9981833-zfact*(0.002875+(zfact*0.0008083)));
printf("sec(z) - 1 polynomial gives %.8g, an error of %.3g%%.\n",
approximation, 100.*fabs(approximation-airmass)/airmass);
/* Do the IRAF approximation (based on J.A. Ball, from Allen) */
zfact=750.*cos(z*DEG2RAD);
approximation=sqrt(zfact*zfact+1501.)-zfact;
printf("IRAF formula gives %.8g, an error of %.3g%%.\n",
approximation, 100.*fabs(approximation-airmass)/airmass);
/* Do the simple airmass ~ sec(z) approximation */
approximation=1./cos(z*DEG2RAD);
printf("sec(z) gives %.8g, an error of %.3g%%.\n",
approximation, 100.*fabs(approximation-airmass)/airmass);
}
/* if desired, compute the (very crude) estimate of vapor column density */
if (relhumid>0.) {
const7=relhumid/(bigrvec[2]-bigrvec[1]);
const6=musquare(*rhovec)*sinz*sinz/(*bigrvec*(*bigrvec)); /* for z */
colz=qromb(vaporcolumndensityint, *bigrvec, bigrvec[2], 1.E-8);
const6=0.; /* equivalent to setting z = 0. */
col0=qromb(vaporcolumndensityint, *bigrvec, bigrvec[2], 1.E-8);
printf("water vapor column density=%.8g g/cm^2; water vapor "
"airmass=%.8g\n", 0.1*colz, colz/col0);
}
exit(0);
}
Reference in New Issue View Git Blame Copy Permalink