diff --git a/bta_print_header/Makefile b/bta_print_header/Makefile index 041c303..5408d01 100644 --- a/bta_print_header/Makefile +++ b/bta_print_header/Makefile @@ -1,7 +1,7 @@ # run `make DEF=...` to add extra defines PROGRAM := bta_print_header LDFLAGS := -fdata-sections -ffunction-sections -Wl,--gc-sections -Wl,--discard-all -LDFLAGS += -lusefull_macros +LDFLAGS += -lusefull_macros -lm -lerfa SRCS := $(wildcard *.c) DEFINES := $(DEF) -D_GNU_SOURCE -D_XOPEN_SOURCE=1111 OBJDIR := mk diff --git a/bta_print_header/airmass.inc b/bta_print_header/airmass.inc new file mode 100644 index 0000000..ebca45d --- /dev/null +++ b/bta_print_header/airmass.inc @@ -0,0 +1,1157 @@ +/******************************************************************************* + * + * 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 ] [-d ] [-H ] [-l ] [-L ] + * [-p

] [-t ] [-v] + * + * where: + * is the APPARENT zenith angle (not the true angle in the + * absence of atmosphere) at which you want the airmass; + * -a sets the altitude of the observer to meters above mean + * sea level (the default is 0 meters); + * -d 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 , 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, . 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 sets the observer's latitude to degrees (negative = south; + * default is 45 degrees); + * -L 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

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 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 +#include +#include + +/******************************************************************************* + * 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=(K-1)) { + ss=polint(h+i-(K-1), s+i-(K-1)); + if (fabs((ss-oldss)/ss)] [-d ] [-H ] [-l ] [-L " + "] [-p

]\n [-t ] [-v] \n\nwhere " + " 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 : observer's altitude above sea level, in meters (default 0)" + "\n\t-d : day number from beginning of year (0 = Jan 1; default " + "80)\n\t-H : relative humidity at observer, in percent (default " + "n/a)\n\t-l : observer's latitude in degrees (default +45)" + "\n\t-L : the observing wavelength in Angstroms (default " + "5500)\n\t-p

: local pressure, in millibars (default 1013.25)" + "\n\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 .\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); +} diff --git a/bta_print_header/am.c b/bta_print_header/am.c new file mode 100644 index 0000000..86cf60d --- /dev/null +++ b/bta_print_header/am.c @@ -0,0 +1,75 @@ +static int main(int argc, char *argv[]) __attribute((unused)); +static void usage(void) __attribute((unused)); +#define inline +#include "airmass.inc" +#undef inline + +#include "bta_site.h" + +void calc_airmass( + // in + double daynum, // Day number from beginning of year, 0 = midnight Jan 1 + double h0,// Relative humidity in percent + double p0, // local pressure in mmHg + double t0, // temperature in kelvins + double z, // zenith distance in degr + // out + double *am, // AIRMASS + double *acd, // column density + double *wcd, // water vapor column density + double *wam // water vapor airmass + ){ + int i; + double colz, col0, cosyearfrac, sinz, rmsl, bigrmsl, + rjunk, coslatsq, lambda=5500., alt=TELALT, + lat=TELLAT, 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.}; + p0 *= 101325. / 760.; + relhumid = h0 / 100.; + 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); + *am = colz/col0; *acd = 0.1*colz; + if(h0 < 0.1){ *wcd = 0.; *wam = 0.; return;} + 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); + *wcd = 0.1*colz; *wam = colz/col0; +} diff --git a/bta_print_header/bta_print.c b/bta_print_header/bta_print.c index 0825164..1b74b13 100644 --- a/bta_print_header/bta_print.c +++ b/bta_print_header/bta_print.c @@ -1,6 +1,6 @@ /* * This file is part of the btaprinthdr project. - * Copyright 2022 Edward V. Emelianov . + * Copyright 2022 Edward V. Emelianov , Vladimir S. Shergin . * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -16,45 +16,98 @@ * along with this program. If not, see . */ +#include #include #include #include #include #include +#include +#include #include "bta_print.h" #include "bta_shdata.h" +#include "bta_site.h" +// rad to time sec +#ifndef ERFA_DR2S +#define ERFA_DR2S 1.3750987083139757010431557155385240879777313391975e4 +#endif +// seconds of time in full circle +#define S24 (24.*3600.) + +const double longitude = TELLONG * 3600.; +const double Fi = TELLAT * 3600.; + +const double cos_fi=0.7235272793; +const double sin_fi=0.6902957888; + +/* +const double cos_fi = 0.723527277857; // Cos of SAO latitude +const double sin_fi = 0.690295790366; // Sin --- "" ----- +*/ + +/** + * convert apparent coordinates (nowadays) to mean (JD2000) + * appRA, appDecl in seconds + * r, d in seconds + */ +void calc_mean(double appRA, double appDecl, double *r, double *dc){ + double ra=0., dec=0., utc1, utc2, tai1, tai2, tt1, tt2, fd, eo, ri; + int y, m, d, H, M; + DBG("appRa: %g'', appDecl'': %g", appRA, appDecl); + appRA *= ERFA_DS2R; + appDecl *= ERFA_DAS2R; +#define ERFA(f, ...) do{if(f(__VA_ARGS__)){WARNX("Error in " #f); goto rtn;}}while(0) + // 1. convert system JDate to UTC + ERFA(eraJd2cal, JDate, 0., &y, &m, &d, &fd); + fd *= 24.; + H = (int)fd; + fd = (fd - H)*60.; + M = (int)fd; + fd = (fd - M)*60.; + ERFA(eraDtf2d, "UTC", y, m, d, H, M, fd, &utc1, &utc2); + ERFA(eraUtctai, utc1, utc2, &tai1, &tai2); + ERFA(eraTaitt, tai1, tai2, &tt1, &tt2); + eraAtic13(appRA, appDecl, tt1, tt2, &ri, &dec, &eo); + ra = eraAnp(ri + eo); + ra *= ERFA_DR2S; + dec *= ERFA_DR2AS; +#undef ERFA +rtn: + if(r) *r = ra; + if(dc) *dc = dec; +} static char buf[1024]; char *time_asc(double t){ - int h, m; - double s; - h = (int)(t/3600.); - m = (int)((t - (double)h*3600.)/60.); - s = t - (double)h*3600. - (double)m*60.; - h %= 24; - if(s>59) s=59; - sprintf(buf, "%d:%02d:%04.1f", h,m,s); - return buf; + int h, m; + double s; + h = (int)(t/3600.); + m = (int)((t - (double)h*3600.)/60.); + s = t - (double)h*3600. - (double)m*60.; + h %= 24; + if(s>59) s=59; + sprintf(buf, "%d:%02d:%04.1f", h,m,s); + return buf; } char *angle_asc(double a){ - char s; - int d, min; - double sec; - if (a >= 0.) - s = '+'; - else{ - s = '-'; a = -a; - } - d = (int)(a/3600.); - min = (int)((a - (double)d*3600.)/60.); - sec = a - (double)d*3600. - (double)min*60.; - d %= 360; - if(sec>59.9) sec=59.9; - sprintf(buf, "%c%d:%02d:%04.1f", s,d,min,sec); - return buf; + char s; + int d, min; + double sec; + if (a >= 0.) + s = '+'; + else{ + s = '-'; a = -a; + } + d = (int)(a/3600.); + min = (int)((a - (double)d*3600.)/60.); + sec = a - (double)d*3600. - (double)min*60.; + d %= 360; + if(sec>59.9) sec=59.9; + sprintf(buf, "%c%d:%02d:%04.1f", s,d,min,sec); + return buf; } /** @@ -73,9 +126,9 @@ static int printhdr(int fd, const char *key, const char *val, const char *cmnt){ key = tk; } if(cmnt){ - snprintf(tmp, 81, "%-8s= %-21s / %s", key, val, cmnt); + snprintf(tmp, 80, "%-8s= %-21s / %s", key, val, cmnt); }else{ - snprintf(tmp, 81, "%-8s= %s", key, val); + snprintf(tmp, 80, "%-8s= %s", key, val); } size_t l = strlen(tmp); tmp[l] = '\n'; @@ -87,10 +140,89 @@ static int printhdr(int fd, const char *key, const char *val, const char *cmnt){ return 0; } +static void calc_AZ(double alpha, double delta, double stime, double *az, double *zd){ + double ha = (stime - alpha) * 15.; + if(ha < 0.) ha += ERFA_TURNAS; + ha *= ERFA_DAS2R; + double A, alt; + eraHd2ae(ha, delta*ERFA_DAS2R, TELLAT*ERFA_DD2R, &A, &alt); + if(az) *az = (A - M_PI) * ERFA_DR2AS; // BTA counts Az from south, not north + if(zd) *zd = (M_PI_2 - alt) * ERFA_DR2AS; +} + +extern void calc_airmass( + // in + double daynum, // Day number from beginning of year, 0 = midnight Jan 1 + double relhumid,// Relative humidity in percent + double p0, // local pressure in mmHg + double t0, // temperature in kelvins + double z, // zenith distance in degr + // out + double *am, // AIRMASS + double *acd, // column density + double *wcd, // water vapor column density + double *wam // water vapor airmass +); + + +#if 0 +//double coeff[8] = { -97.4, 13.0, -11.7, -3.6, -6.2, -279.9, -33.3, 8.2}; +void pos_corr (double Delta, double A, double Z, double P, + double *pcA, double *pcZ) { + double ar, dr, zr, pr; + + dr = Delta; + ar = A; + zr = Z; + pr = P; +double *coeff = (double*)sdt->pc_coeff; + *pcA = coeff[0] + coeff[1]/tan(zr) + coeff[2]/sin(zr) - + coeff[3]*sin(ar)/tan(zr) + coeff[4]*cos(dr)*cos(pr)/sin(zr); + + *pcZ = coeff[5] + coeff[6]*sin(zr) + coeff[7]*cos(zr) + + coeff[3]*cos(ar) + coeff[4]*cos_fi*sin(ar); + +} + +void calc_AZP(double alpha, double delta, double stime, + double *az, double *zd, double *pa){ + double sin_t,cos_t, sin_d,cos_d, sin_z,cos_z, sin_p,cos_p; + double t, d, z, a, p , x, y; + + t = (stime - alpha) * 15.; + red("st=%g, alpha=%g, delta=%g, t=%g\n", stime/3600., alpha/3600., delta/3600., t/3600./15.); + if (t < 0.) + t += ERFA_TURNAS; + t *= ERFA_DAS2R; + d = delta * ERFA_DAS2R; + sincos(t, &sin_t, &cos_t); + sincos(d, &sin_d, &cos_d); + + cos_z = cos_fi * cos_d * cos_t + sin_fi * sin_d; + z = acos(cos_z); + sin_z = sqrt(1. - cos_z*cos_z); + + y = cos_d * sin_t; + x = cos_d * sin_fi * cos_t - cos_fi * sin_d; + a = atan2(y, x); + + sin_p = sin_t * cos_fi / sin_z; + cos_p = (sin_fi * cos_d - sin_d * cos_fi * cos_t) / sin_z; + p = atan2(sin_p, cos_p); + if (p < 0.0) + p += 2.0*M_PI; +double pca, pcz; +pos_corr(d, a, z, p, &pca, &pcz); +red("pca=%g, pcz=%g; tel_cor_A=%g, tel_cor_Z=%g\n", pca, pcz, tel_cor_A, tel_cor_Z); + *zd = z * ERFA_DR2AS; + *az = a * ERFA_DR2AS; + *pa = p * ERFA_DR2AS; +} +#endif + #define WRHDR(k, v, c) do{if(printhdr(fd, k, v, c)){goto returning;}}while(0) int print_header(_U_ const char *path){ int ret = FALSE; - double dtmp; char val[23], comment[71]; #define COMMENT(...) do{snprintf(comment, 70, __VA_ARGS__);}while(0) #define VAL(fmt, x) do{snprintf(val, 22, fmt, x);}while(0) @@ -107,33 +239,99 @@ int print_header(_U_ const char *path){ } fchmod(fd, 0644); WRHDR("TELESCOP", "'BTA 6m telescope'", "Telescope name"); - WRHDR("ORIGIN", "'SAO RAS'", "Organization responsible for the data"); - dtmp = S_time-EE_time; - COMMENT("Sidereal time, seconds: %s", time_asc(dtmp)); - VALD(dtmp); + WRHDR("ORIGIN", "'SAO RAS, Russia'", "Organization responsible for the data"); + VALD(TELLAT); + COMMENT("Telescope lattitude (degr): %s", angle_asc(TELLAT*3600.)); + WRHDR("SITELAT", val, comment); + VALD(TELLONG); + COMMENT("Telescope longitude (degr): %s", angle_asc(TELLONG*3600.)); + WRHDR("SITELONG", val, comment); + VAL("%.1f", TELALT); + WRHDR("SITEALT", val, "Telescope altitude (m)"); + double sidtm = S_time-EE_time; + if(sidtm < 0.) sidtm += S24; + else if(sidtm > S24){ + int x = (int)(sidtm / S24); + sidtm -= S24 * (double)x; + } + COMMENT("Sidereal time, seconds: %s", time_asc(sidtm)); + VALD(sidtm); WRHDR("ST", val, comment); COMMENT("Universal time, seconds: %s", time_asc(M_time)); VALD(M_time); WRHDR("UT", val, comment); VALD(JDate); WRHDR("JD", val, "Julian date"); + {time_t t_now = time(NULL); + struct tm *tm_ut = gmtime(&t_now); + double l = ERFA_DTY; int y = 1900 + tm_ut->tm_year; + if((0 == y%400) || ((0 == y%4)&&(0 != y%100))) l += 1.; + double dtmp = y + (double)tm_ut->tm_yday / l; + VALD(dtmp); + WRHDR("EQUINOX", val, "Coordinates epoch");} switch(Tel_Focus){ - default: - case Prime: - VALS("Prime"); break; - case Nasmyth1: + default: + case Prime: + VALS("Prime"); break; + case Nasmyth1: VALS("Nasmyth1"); break; - case Nasmyth2: + case Nasmyth2: VALS("Nasmyth2"); break; - } + } WRHDR("FOCUS", val, "Observation focus"); + char *telmode = "Undefined"; + if(Tel_Hardware == Hard_Off) telmode = "Off"; + else if(Tel_Mode != Automatic) telmode = "Manual"; + else{ + switch (Sys_Mode){ + default: + case SysStop : telmode = "Stopping"; break; + case SysWait : telmode = "Waiting"; break; + case SysPointAZ : + case SysPointAD : telmode = "Pointing"; break; + case SysTrkStop : + case SysTrkStart: + case SysTrkMove : + case SysTrkSeek : telmode = "Seeking"; break; + case SysTrkOk : telmode = "Tracking"; break; + case SysTrkCorr : telmode = "Correction";break; + case SysTest : telmode = "Testing"; break; + } + } + VALS(telmode); + WRHDR("TELMODE", val, "Telescope working mode"); + VAL("%.2f", val_F); WRHDR("VAL_F", val, "Focus value of telescope (mm)"); -#define RA(ra, dec, text, pref) do{VALD(ra*15./3600.); COMMENT(text " R.A. (degr): %s", angle_asc(15.*ra)); WRHDR(pref "RA", val, comment); \ - VALD(dec/3600.); COMMENT(text " Decl (degr): %s", angle_asc(dec)); WRHDR(pref "DEC", val, comment);}while(0) - RA(CurAlpha, CurDelta, "Current object", ""); - RA(SrcAlpha, SrcDelta, "Source", "S_"); - RA(val_Alp, val_Del, "Telescope", "T_"); + double a2000, d2000; + calc_mean(InpAlpha, InpDelta, &a2000, &d2000); + VALD(a2000 * 15. / 3600.); + COMMENT("Input R.A. for J2000 (deg): %s", time_asc(a2000)); + WRHDR("RA_INP0", val, comment); + VALD(d2000 / 3600.); + COMMENT("Input Decl. for J2000 (deg): %s", angle_asc(d2000)); + WRHDR("DEC_INP0", val, comment); + calc_mean(CurAlpha, CurDelta, &a2000, &d2000); + VALD(a2000 * 15. / 3600.); + COMMENT("Telescope R.A. for J2000 (deg): %s", time_asc(a2000)); + WRHDR("RA_0", val, comment); + VALD(d2000 / 3600.); + COMMENT("Telescope Decl. for J2000 (deg): %s", angle_asc(d2000)); + WRHDR("DEC_0", val, comment); +#define RA(ra, dec, text, pref) do{VALD(ra*15./3600.); COMMENT(text " R.A. (degr): %s", time_asc(ra)); WRHDR("RA" pref, val, comment); \ + VALD(dec/3600.); COMMENT(text " Decl (degr): %s", angle_asc(dec)); WRHDR("DEC" pref, val, comment);}while(0) + RA(InpAlpha, InpDelta, "Input", "_INP"); + RA(CurAlpha, CurDelta, "Current object", "_OBJ"); + RA(SrcAlpha, SrcDelta, "Source", "_SRC"); + RA(val_Alp, val_Del, "Telescope", ""); +#undef RA +#define AZ(a, z, text, pref) do{VALD(a/3600.); COMMENT(text " Az (degr): %s", angle_asc(a)); WRHDR("A" pref, val, comment); \ + VALD(z/3600.); COMMENT(text " ZD (degr): %s", angle_asc(z)); WRHDR("Z" pref, val, comment);}while(0) + AZ(InpAzim, InpZdist, "Input", "_INP"); + AZ(tag_A, tag_Z, "Current object", "_OBJ"); + AZ(val_A, val_Z, "Telescope", ""); +#undef AZ +/* VALD(tag_A/3600.); COMMENT("Target Az (degr): %s", angle_asc(tag_A)); WRHDR("TARG_A", val, comment); @@ -146,25 +344,133 @@ int print_header(_U_ const char *path){ VALD(val_Z/3600.); COMMENT("Telescope ZD (degr): %s", angle_asc(val_Z)); WRHDR("Z", val, comment); + */ + double t = sidtm * ERFA_DS2R - val_Alp * ERFA_DS2R; + if(t < 0) t += 2*M_PI; + double P = ERFA_DR2D * eraHd2pa(t, val_Del * ERFA_DAS2R, TELLAT * ERFA_DD2R); + if(P < 0.) P += 360.; + VALD(P); + COMMENT("Parallactic angle (degr): %s", angle_asc(P*3600.)); + WRHDR("PARANGLE", val, comment); VALD(tag_P/3600.); - COMMENT("Target rotation angle (degr): %s", angle_asc(tag_P)); - WRHDR("TARG_R", val, comment); + COMMENT("Target par. angle (degr): %s", angle_asc(tag_P)); + WRHDR("TAGANGLE", val, comment); VALD(val_P/3600.); - COMMENT("Telescope rotation angle (degr): %s", angle_asc(val_P)); + COMMENT("Current P2 rot. angle (degr): %s", angle_asc(val_P)); WRHDR("ROTANGLE", val, comment); VALD(val_D/3600.); COMMENT("Dome Az (degr): %s", angle_asc(val_D)); WRHDR("DOME_A", val, comment); + if(Sys_Mode==SysTrkSeek||Sys_Mode==SysTrkOk||Sys_Mode==SysTrkCorr||Sys_Mode==SysTrkStart||Sys_Mode==SysTrkMove){ + double curA,curZ,srcA,srcZ; + double corAlp,corDel,corA,corZ; + corAlp = CurAlpha-SrcAlpha; + corDel = CurDelta-SrcDelta; + if(corAlp > 23*3600.) corAlp -= 24*3600.; + if(corAlp < -23*3600.) corAlp += 24*3600.; + calc_AZ(SrcAlpha, SrcDelta, S_time, &srcA, &srcZ); + calc_AZ(CurAlpha, CurDelta, S_time, &curA, &curZ); + corA = curA-srcA; + corZ = curZ-srcZ; + VALD(corAlp); + WRHDR("RACORR", val, "RA correction (current - source)"); + VALD(corDel); + WRHDR("DECCORR", val, "DEC correction (current - source)"); + VALD(corA); + WRHDR("ACORR", val, "A correction (current - source)"); + VALD(corZ); + WRHDR("ZCORR", val, "Z correction (current - source)"); + } +/* + double az,zd;//,pa; + calc_AZ(val_Alp, val_Del, sidtm, &az, &zd); + green("calc_AZ: AZ=%g, ZD=%g\n", az/3600., zd/3600.); + calc_AZx(val_Alp, val_Del, sidtm, &az, &zd); + green("calc_AZx: AZ=%g, ZD=%g\n", az/3600., zd/3600.); + calc_AZP(val_Alp, val_Del, sidtm, &az, &zd, &pa); + green("AZ=%g, ZD=%g, PA=%g; sidtm=%g\n", az/3600.,zd/3600.,pa/3600., sidtm/3600.); +*/ + if(UsePCorr) VALS("true"); else VALS("false"); + WRHDR("USEPCORR", val, "P.corr.sys.: K0..K7 (real = measured - PCS)"); + if(UsePCorr){ + char kx[3]; + const char *descr[8] = { + "A0 - PCS Azimuth zero", + "L - PCS horizontal axe inclination", + "k - PCS collimation error", + "F - PCS lattitude error of vert. axe", + "dS - PCS time error", + "Z0 - PCS zenith zero", + "d - PCS tube bend [sin(Z)]", + "d1 - PCS tube bend [cos(Z)]" + }; + for(int i = 0; i < 8; ++i){ + VAL("%.2f", PosCor_Coeff[i]); + snprintf(kx, 3, "K%d", i); + WRHDR(kx, val, descr[i]); + } +#if 0 + // SKN: + // dA = K0 + K1/tg(Z) + K2/sin(Z) - K3*sin(A)/tg(Z) + K4 *cos(delta)*cos(P)/sin(Z) + // dZ = K5 + K6*sin(Z) + K7*cos(Z) + K3*cos(A) + K4*cos(phi)*sin(A) + double A = val_A * ERFA_DAS2R, Z = val_Z * ERFA_DAS2R; + double sinZ, cosZ, sinA, cosA, tgZ = tan(Z); + sincos(Z, &sinZ, &cosZ); + sincos(A, &sinA, &cosA); + double dA = PosCor_Coeff[0] + PosCor_Coeff[1]/tgZ + PosCor_Coeff[2]/sinZ - + PosCor_Coeff[3]*sinA/tgZ + + PosCor_Coeff[4]*cos(CurDelta*ERFA_DAS2R)*cos(P*ERFA_DD2R)/sinZ; + double dZ = PosCor_Coeff[5] + PosCor_Coeff[6]*sinZ + PosCor_Coeff[7]*cosZ + + PosCor_Coeff[3]*cosA + PosCor_Coeff[4]*cos_fi*sinA; + red("dA=%g, dZ=%g; tel_cor_A=%g, tel_cor_Z=%g\n", dA, dZ, tel_cor_A, tel_cor_Z); +#endif + VAL("%.1f", pos_cor_A); + WRHDR("PCSDA_O", val, "Pointing correction for object A (arcsec)"); + VAL("%.1f", pos_cor_Z); + WRHDR("PCSDZ_O", val, "Pointing correction for object Z (arcsec)"); + VAL("%.1f", tel_cor_A); + WRHDR("PCSDA_T", val, "Pointing correction for telescope A (arcsec)"); + VAL("%.1f", tel_cor_Z); + WRHDR("PCSDZ_T", val, "Pointing correction for telescope Z (arcsec)"); + } + VAL("%.1f", refract_Z); + WRHDR("REFR_O", val, "Refraction for object position (arcsec)"); + VAL("%.1f", tel_ref_Z); + WRHDR("REFR_T", val, "Refraction for telescope position (arcsec)"); + {double refa, refb, tz = tan(val_Z*ERFA_DAS2R); + eraRefco(val_B*1.33322, val_T1, val_Hmd/100., 0.45, &refa, &refb); + double dz = (refa*tz + refb*tz*tz*tz)*ERFA_DR2AS; + VAL("%.2f", dz); + WRHDR("REFR_T_E", val, "RREFR_T by ERFA eraRefco()");} + #define T(hdr, t, text) do{VAL("%.1f", t); COMMENT(text " temperature (degC)"); WRHDR(hdr, val, comment);}while(0) T("OUTTEMP", val_T1, "Outern"); T("DOMETEMP", val_T2, "In-dome"); T("MIRRTEMP", val_T3, "Mirror"); +#undef T VAL("%.1f", val_B); WRHDR("PRESSURE", val, "Atm. pressure (mmHg)"); VAL("%.1f", val_Wnd); WRHDR("WIND", val, "Wind speed (m/s)"); - VAL("%.0f", val_Hmd); + VAL("%.1f", val_Hmd); WRHDR("HUMIDITY", val, "Relative humidity (%)"); + /* + * Airmass calculation + * by Reed D. Meyer + */ + {double am, acd, wcd, wam; + time_t t_now = time(NULL); + struct tm *tm_loc; + tm_loc = localtime(&t_now); + calc_airmass(tm_loc->tm_yday, val_Hmd, val_B, val_T1+273.15, val_Z/3600., &am, &acd, &wcd, &wam); + VALD(am); + WRHDR("AIRMASS", val, "Air mass by Reed D. Meyer"); + VALD(wam); + WRHDR("WVAM", val, "Water vapour air mass by Reed D. Meyer"); + VALD(acd); + WRHDR("ATMDENS", val, "Atm. column density by Reed D. Meyer (g/cm^2)"); + VALD(wcd); + WRHDR("WVDENS", val, "WV column density by Reed D. Meyer (g/cm^2)");} ret = TRUE; returning: close(fd); diff --git a/bta_print_header/bta_site.h b/bta_print_header/bta_site.h new file mode 100644 index 0000000..1afcb06 --- /dev/null +++ b/bta_print_header/bta_site.h @@ -0,0 +1,38 @@ +/* + * This file is part of the btaprinthdr project. + * Copyright 2025 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +/* + * SAO longitude 41 26 29.175 + * SAO latitude 43 39 12.7 + * SAO altitude 2070 + * BTA focal ratio 24.024 m + */ +#ifndef TELLAT + #define TELLAT (43.6535278) +#endif +#ifndef TELLONG + #define TELLONG (41.4414375) +#endif +#ifndef TELALT + #define TELALT (2070.0) +#endif +#ifndef TELFOCUS + #define TELFOCUS (24.024) +#endif diff --git a/bta_print_header/btaprinthdr.creator.user b/bta_print_header/btaprinthdr.creator.user index 8e8c7d5..0344072 100644 --- a/bta_print_header/btaprinthdr.creator.user +++ b/bta_print_header/btaprinthdr.creator.user @@ -1,6 +1,6 @@ - + EnvironmentId @@ -8,7 +8,7 @@ ProjectExplorer.Project.ActiveTarget - 0 + 0 ProjectExplorer.Project.EditorSettings @@ -28,7 +28,7 @@ QmlJSGlobal - 2 + 2 KOI8-R false 4 @@ -37,6 +37,8 @@ true true 1 + 0 + false true false 1 @@ -45,29 +47,56 @@ 0 8 true + false 2 true false true + *.md, *.MD, Makefile false + true + true ProjectExplorer.Project.PluginSettings + + true + false + true + true + true + true + + + 0 + true true + + true + true + Builtin.DefaultTidyAndClazy + 4 + true + + + + true + ProjectExplorer.Project.Target.0 + Desktop Desktop Desktop {91347f2c-5221-46a7-80b1-0a054ca02f79} - 0 - 0 - 0 + 0 + 0 + 0 /tmp/bgd/C/bta_control_net-x86_64/bta_print_header @@ -75,14 +104,10 @@ all - false - - - false true GenericProjectManager.GenericMakeStep - 1 + 1 Сборка Сборка ProjectExplorer.BuildSteps.Build @@ -92,117 +117,60 @@ clean - true - - - false true GenericProjectManager.GenericMakeStep - 1 + 1 Очистка Очистка ProjectExplorer.BuildSteps.Clean 2 false + + false По умолчанию GenericProjectManager.GenericBuildConfiguration - 1 + 1 - 0 + 0 Развёртывание Развёртывание ProjectExplorer.BuildSteps.Deploy 1 + + false ProjectExplorer.DefaultDeployConfiguration - 1 - + 1 - dwarf - - cpu-cycles - - - 250 - - -e - cpu-cycles - --call-graph - dwarf,4096 - -F - 250 - - -F true - 4096 - false - false - 1000 - true - - false - false - false - false - true - 0.01 - 10 - true - kcachegrind - 1 - 25 - - 1 true true - true - valgrind - - 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 10 - 11 - 12 - 13 - 14 - + + 2 - + false + -e cpu-cycles --call-graph dwarf,4096 -F 250 + ProjectExplorer.CustomExecutableRunConfiguration - - false - - false + false true - false - false true - - - 1 + 1 ProjectExplorer.Project.TargetCount - 1 + 1 ProjectExplorer.Project.Updater.FileVersion diff --git a/bta_print_header/btaprinthdr.files b/bta_print_header/btaprinthdr.files index 332ff40..e4c29d2 100644 --- a/bta_print_header/btaprinthdr.files +++ b/bta_print_header/btaprinthdr.files @@ -1,7 +1,10 @@ +airmass.inc +am.c bta_print.c bta_print.h bta_shdata.c bta_shdata.h +bta_site.h cmdlnopts.c cmdlnopts.h main.c