A major occupation of the astrophysicist is that of making attempts to understand the physical structure of heavenly bodies for which he has quantitative observations. In most cases his procedure is indirect: he sets up a mathematical model of the relevant features of the body of interest, and compares its predictions with the observations. If they do not fit sufficiently well, he tries to make changes to improve the agreement. Unfortunately, there are many cases in which the mathematical analysis of even the simplest models is so involved that drastic simplifications are needed to obtain an answer. Such analyses often contain inconsistencies, and one has frequently little confidence in the precision of their results.
The situation we have described is typified by many analyses of the structure of the Venus cloud layer (1,2,3). It seemed to us to be worthwhile to attempt a study, in much greater depth than one could contemplate with the desk methods hitherto employed, by making use of the great power of the modern computer. The data provided by Spinrad (3) on the strengths of absorption lines in the P-branch of the carbon dioxide rotation-vibration band at A 7820 Angstroms at various Venus phases seemed the most suitable. The lines are formed purely by diffuse reflection of sunlight; our aim is to construct a model of this process and to determine the optimum parameters of the atmospheric models which fit the observations best. (It remains to be seen whether the data are good enough to support this procedure. However, our main aim is to demonstrate that our methods provide a practical solution of the problem, and perhaps to stimulate more accurate observation should this prove to be necessary.)
There are, as we have said, two parts to the problem: first, the construction of atmospheric models; second, an efficient and economical optimization procedure. We have written up our completed work on the first part (4) and are currently studying the second problem.
We make the following assumptions:
A. Atmosphere
Six parameters are required to completely specify the model, namely the temperature (To) and pressure (Po) near the cloud top, the acceleration due to gravity (g), the temperature at the ground (Tg), the temperature lapse rate (L-1 in km-1), and the mean molecular weight (μ). Estimates are able of the first four of these, guesses for the last two (5). The composition is assumed to be in the region of 20% CO2, 80% N2 by volume.
B. Optical Properties
Absorption primarily due to CO2. The line with rotational quantum number J, centred on the frequency vJ, has an absorption coefficient for radiation of frequency v at temperature T, pressure p given by
(1)
where , β is a molecular constant.
Cloud particles have number density proportional to gas density. They are assumed to scatter radiation isotropically with an albedo near unity (cf. (2)). The scattering cross-section is known to vary in a complex manner as the frequency changes, but we assume it is constant over the spectral range of interest here.
Combining assumptions A and B, we obtain an optical thickness variation with pressure τv(p), and the behaviour of the albedo for single scattering ωv(τv) (the ratio of scattering to (scattering + absorption) cross-sections for frequency v at optical height τv ). The diffusely scattered intensity emerging at the top of the atmosphere can be obtained by computing the solution of the well-known integro-differential equation for the diffuse reflection function Sv (μ, μ0;τ) ) at the top of the clouds, τ = τ (p0), where
where
(2)
prescribed.
Our solution of (1) depends on reducing it to a set of ordinary non-linear differential equations in τ by approximating the r.h.s. of (2) by a quadrature formula. A new method of quadrature, based on a quadratic spline fit to Sv (μ, μ';τ) as μ' varies, is chosen to facilitate the subsequent computation of integrals over awkward ranges of μ and μ0. With the range 0 ≤ μ,μ0 ≤ 1 split into N intervals, we then have to solve ½N(N+1) simultaneous equations (symmetry in μ and μ0). Error bounds have been obtained and verified; with N = 8, the solution for
is in error by less than ½%.
A code has been constructed to perform these calculations on Atlas in the Hartran language. For N = 8, the numerical solution of (1) at a fixed value of v requires about 4 sec for 70-100 steps in Τ. Thus a typical run with, say, 23 frequencies in each of 5 lines takes of order 10 minutes. The complete results of the calculation are filed on magnetic tape to eliminate tedious data handling. Subsequent processing, for example to compute equivalent widths at the observed Venus phase angles, depends on solutions stored in this way. Recovery of the data is expedited by using the block addressing facility on Ampex tape in conjunction with the Hartran BLOCK COMMON feature. Machine-coded subroutines are used to transmit buffered data in 512-word blocks to and from tape. An index of the current state of the tape is maintained in block 1 and is updated at the end of each machine session. Data for many models can be accommodated on a single tape, a useful feature for optimization calculations.
Provision has been made for the automatic production of graphs via IBM 1inch magnetic tape and Benson-Lehner electroplotter. As this doubles the running time, the facility has been little used.
Two problems arise in attempting to find models that will best fit the observations. One is the selection of an efficient, economical, optimization procedure; the other is the selection of convergence criteria. For the procedure, we are attempting to devise a least-square strategy which minimizes the number of models to be constructed, since it is this part of the problem that takes most of the machine time. At present, the most favoured scheme is a combination of that due to Powell (6) with that due to Marquardt (7). The former scheme is a method for minimizing sums of squares without using derivatives. The latter attempts to secure rapid convergence from an initially poor guess by slightly warping the surface on which minimization is sought.
For the convergence criterion, we remark that the observational data are subject to fluctuations arising from errors of measurement and from day to day physical changes in the Venus atmosphere. This defines a confidence region in which we accept all models giving predictions within a given "distance" of the observations. The "distance" is determined by statistical arguments. We hope that such tests will enable us to decide whether the data are sufficiently well-defined to distinguish between models of the type we have described. Later, we hope to consider models of more general physical type and to use the tests to discriminate between them.
1. J. W. Chamberlain and G. P. Kuiper, Astrophys. J. 124, 399, 1956.
2. J.W. Chamberlain, Astrophys.J. 141, 1184, 1965.
3. H. Spinrad, Proc. Astron. Soc. Pacific 74, 187, 1962.
4. 1. P. Grant and V. M. Burke, J. Computational Phys., to be published.
5. L. D. Kaplan, J. Quant. Spectrosc. and Rad. Transf. 3, 537, 1963; and Jet Propulsion Lab. Report no. 32-379, Dec. 1962.
6. M.J.D. Powell, Comp.J., 7, 303, 1965.
7. D.W. Marquardt, SIAM J., 11, 431, 1963.
