We lay out a general framework for calculating the variation of a set of cosmological observables, down the past null cone of an arbitrarily placed observer, in a given arbitrary inhomogeneous metric. The observables include redshift, proper motions, area distance and redshift-space density. Of particular interest are observables that are zero in the spherically symmetric case, such as proper motions. The algorithm is based on the null geodesic equation and the geodesic deviation equation, and it is tailored to creating a practical numerical implementation. The algorithm provides a method for tracking which light rays connect moving objects to the observer at successive times. Our algorithm is applied to the particular case of the Szekeres metric. A numerical implementation has been created and some results will be presented in a subsequent paper. Future work will explore the range of possibilities.

The Szekeres family of inhomogeneous solutions, which are defined by six arbitrary metric functions, offers a wide range of possibilities for modelling cosmic structure. Here we present a model construction procedure for the quasispherical case using given data at initial and final times. Of the six arbitrary metric functions, the three which are common to both Szekeres and Lema\^itre-Tolman models are determined by the model construction procedure of Krasinski & Hellaby. For the remaining three functions, which are unique to Szekeres models, we derive exact analytic expressions in terms of more physically intuitive quantities - density profiles and dipole orientation angles. Using MATLAB, we implement the model construction procedure and simulate the time evolution.

This paper presents the algorithm for determining the Lemaitre-Tolman (LT) model that best fits given datasets for maximum stellar ages, and SNIa luminosities, both as functions of redshift. It then applies it to current cosmological data. Special attention must be given to the handling of the origin, and the region of the maximum diameter distances. As with a previous combination of datasets (galaxy number counts and luminosity distances versus redshift), there are relationships that must hold at the region of the maximum diameter distance, which are unlikely to be obeyed exactly by real data. We show how to make corrections that enable a self-consistent solution to be found. We address the questions of the best way to approximate discrete data with smooth functions, and how to estimate the uncertainties of the output - the 3 free functions that determine a specific LT metric. While current data does not permit any confidence in our results, we show that the method works well, and reasonable LT models do fit with or without a cosmological constant.

A number of misunderstandings about modeling the apparent accelerated expansion of the Universe, and about the `weak singularity' are clarified: 1. Of the five definitions of the deceleration parameter given by Hirata and Seljak (HS), only $q_1$ is a correct invariant measure of acceleration/deceleration of expansion. The $q_3$ and $q_4$ are unrelated to acceleration in an inhomogeneous model. 2. The averaging over directions involved in the definition of $q_4$ does not correspond to what is done in observational astronomy. 3. HS's equation (38) connecting $q_4$ to the flow invariants gives self-contradictory results when applied at the centre of symmetry of the Lema\^{\i}tre-Tolman (L-T) model. The intermediate equation (31) that determines $q_{3'}$ is correct, but approximate, so it cannot be used for determining the sign of the deceleration parameter. Even so, at the centre of symmetry of the L-T model, it puts no limitation on the sign of $q_{3'}(0)$. 4. The `weak singularity' of Vanderveld {\it et al.} is a conical profile of mass density at the centre - a perfectly acceptable configuration. 5. The so-called `critical point' in the equations of the `inverse problem' for a central observer in an L-T model is a manifestation of the apparent horizon - a common property of the past light cones in zero-lambda L-T models, perfectly manageable if the equations are correctly integrated.

We consider the spherically symmetric metric with a comoving perfect fluid and non-zero pressure -- the Lemaitre metric -- and present it in the form of a calculational algorithm. We use it to review the definition of mass, and to look at the apparent horizon relations on the observer's past null cone. We show that the introduction of pressure makes it difficult to separate the mass from other physical parameters in an invariant way. Under the usual mass definition, the apparent horizon relation, that relates the diameter distance to the cosmic mass, remains the same as in the Lemaitre-Tolman case.

An overview of some recent developments in inhomogeneous models is presented. As the volume and precision of cosmological data improves, it will become more and more essential to understand the non-linear behaviour of the Einstein field equations. This requires the study of exact inhomogeneous solutions, including their density distributions, their evolution, their geometry, and their causal structure. Observations are strongly affected by the detailed geometry and evolution of a model, and therefore interpretation of observations depends on understanding them. It is generally assumed the universe is homogeneous if averaged over large enough scales, but to actually prove this is so, will require the assumption to be relaxed, and a rigorous inhomogeneous approach to be applied. Though the \LT metric has long been used for models of spherical inhomogeneities, there have been a number of new results, including a variety of methods for creating models with specific properties, and their application to cosmic structures on several different scales. Interest in the Szekeres metrics is on the increase, and the quasi-spherical metric was recently used to model specific cosmic structures for the first time. The quasi-planar and quasi-hyperspherical metrics have been hardly studied until recent work invesigated their physical and geometric properties. There is enormous scope for work with these metrics.

The ultimate application of Einstein's field equations is to empirically determine the geometry of the Universe from its matter content, rather than simply assuming the Universe can be represented by a homogeneous model on all scales. Choosing an LTB model as the most convenient inhomogeneous model for the early stages of development, a data reduction procedure was recently validated using perfect test data. Here we simulate observational uncertainties and improve the previous numerical scheme to ensure that it will be usable with real data as soon as observational surveys are sufficiently deep and complete. Two regions require special treatment--the origin and the maximum in the areal radius. To minimize numerical errors near the origin, we use an LTB series expansion to provide the initial values for integrating the differential equations. We also use an improved method to match the numerical integration to the series expansion that bridges the region near the maximum in the areal radius. Because the mass enclosed within the maximum obeys a specific relationship, we show that it is possible to correct for a fixed systematic error in either the distance scale or the redshift-space mass density, such that the integrated values are consistent with the data at the maximum.

Recent galaxy redshift surveys have brought in a large amount of accurate cosmological data out to redshift 0.3, and future surveys are expected to achieve a high degree of completeness out to a redshift exceeding 1. Consequently, a numerical programme for determining the metric of the universe from observational data will soon become practical; and thereby realise the ultimate application of Einstein's equations. Apart from detailing the cosmic geometry, this would allow us to verify and quantify homogeneity, rather than assuming it, as has been necessary up to now, and to do that on a metric level, and not merely at the mass distribution level. This paper is the beginning of a project aimed at such a numerical implementation. The primary observational data from our past light cone consists of galaxy redshifts, apparent luminosities, angular diameters and number densities, together with source evolution functions, absolute luminosities, true diameters and masses of sources. Here we start with the simplest case, that of spherical symmetry and a dust equation of state, and execute an algorithm that determines the unknown metric functions from this data. We discuss the challenges of turning the theoretical algorithm into a workable numerical procedure, particularly addressing the origin and the maximum in the area distance. Our numerical method is tested with several artificial data sets for homogeneous and inhomogeneous models, successfully reproducing the original models. This demonstrates the basic viability of such a scheme. Although current surveys don't have sufficient completeness or accuracy, we expect this situation to change in the near future, and in the meantime there are many refinements and generalisations to be added.

We point out that the mass of the cosmos on gigaparsec scales can be measured, owing to the unique geometric role of the maximum in the areal radius. Unlike all other points on the past null cone, this maximum has an associated mass, which can be calculated with very few assumptions about the cosmological model, providing a measurable characteristic of our cosmos. In combination with luminosities and source counts, it gives the bulk mass to light ratio. The maximum is particularly sensitive to the values of the bulk cosmological parameters. In addition, it provides a key reference point in attempts to connect cosmic geometry with observations. We recommend the determination of the distance and redshift of this maximum be explicitly included in the scientific goals of the next generation of reshift surveys. The maximum in the redshift space density provides a secondary large scale characteristic of the cosmos.

We develop a non-linear relativistic model of the Shapley Concentration (SC) and its environs, including the Great Attractor (GA) and the Local Group (LG). We take the Shapley concentration as a major attractive centre, and we use the Lemaitre-Tolman model. We constrain our model with present day observations, plus the requirement that it have a physically reasonable evolution from small perturbations at last scattering. We investigate possible mass and velocity distributions, and we find that the peculiar velocity maximum near the SC is \~800 km/s inwards, the density between GA and SC must be about ~0.9 times background, the mass of the GA is probably 4 - 6 x 10^15 M_sun, the SC's contribution to the LG motion is negligible, and the value of the cosmological constant is not significant on this scale.

We develop models of void formation starting from a small initial fluctuation at recombination and growing to a realistic present day density profile in agreement with observations of voids. The model construction is an extension of previously developed algorithms for finding a Lemaitre-Tolman metric that evolves between two profiles of either density or velocity specified at two times. Of the 4 profiles of concern -- those of density and velocity at recombination and at the present day -- two can be specified and the other two follow from the derived model. We find that, in order to reproduce the present-day void density profiles, the initial velocity profile is more important than the initial density profile. Extrapolation of current CMB observations to the scales relevant to proto-voids is very uncertain. Even so, we find that it is very difficult to make both the initial density and velocity fluctuation amplitudes small enough, and still obtain a realistic void by today.

Structure formation within the Lemaitre-Tolman model is investigated in a general manner. We seek models such that the initial density perturbation within a homogeneous background has a smaller mass than the structure into which it will develop, and the perturbation then accretes more mass during evolution. This is a generalisation of the approach taken by Bonnor in 1956. It is proved that any two spherically symmetric density profiles specified on any two constant time slices can be joined by a Lemaitre-Tolman evolution, and exact implicit formulae for the arbitrary functions that determine the resulting L-T model are obtained. Examples of the process are investigated numerically.

We propose a method of testing source evolution theories that is independent of the effects of inhomogeneity, and thus complementary to other studies of evolution. It is suitable for large scale sky surveys, and the new generation of large telescopes. In an earlier paper it was shown that basic cosmological observations - luminosity versus redshift, area distance versus redshift and number counts versus redshift - cannot separate the effects of cosmic inhomogeneity, cosmic evolution and source evolution. We here investigate multicolour observations, and show that by comparing luminosity versus redshift in two or more colours, contraints can be placed on source evolution even if unknown source evolution is present, providing an important test of evolution theories that is complementary to present methods. However, number counts in different colours versus redshift are not useful in separating the effects of source evolution and inhomogeneity.

We consider a spherically symmetric distribution of dust and show that it is possible, under general physically reasonable conditions, for an overdensity to evolve to an underdensity (and vice versa). We find the conditions under which this occurs and illustrate it on a class of regular Lemaitre-Tolman-Bondi solutions. The existence of this phenomenon, if verified, would have the result that the topology of density contours, assumed fixed in standard structure formation theories, would have to change and that luminous matter would not trace the dark matter distribution so well.

This paper and the others in the series challenge the standard model of the effects of gravitational lensing on observations at large distances. We show that due to the cumulative effect of lensing, areas corresponding to an observed solid angle can be quite different than would be estimated from the corresponding Friedmann-Lema\^{\i}tre model, even when averaged over large angular scales. This paper concentrates on the specific example of spherically symmetric but spatially inhomogeneous dust universes, the Lema\^{\i}tre-Tolman-Bondi models, and shows that radial lensing significantly distorts the area distance-redshift and density-redshift relations in these exact solutions compared with the standard ones for Friedmann-Lema\^{\i}tre models. Thus inhomogeneity may introduce significant errors into distance estimates based on the standard FL relations, even after all-sky averaging. In addition a useful new gauge choice is presented for these models, solving the problem of locating the past null cone exactly.