# Refractive Bending and Delay Calculations: Background Information and ALMA Implementation

Last Update: JeffMangum - 16 August 2009

## Summary and Recommendations

• See "Atmospheric Refractive Signal Bending and Propagation Delay" (last update 2008/08/19) for background and specific recommendations for calculating signal bending and delay due to the atmosphere.
• CALC calculates the total delay between two "stations". This delay is composed of two components:
• Geometric Delay
• Instrumental Delay
• Atmospheric Delay
• The uncertainty in the geometric delay calculation is dominated by uncertainties in the antenna location (baseline) determination. See ALMA Memo 503 (Conway 2004) for a discussion.
• The calculation that CALC makes does not have to be perfect, only good enough such that residual phase rate in frequency and time will produce less than 1% coherence loss in the data. One only needs to make sure that the residual delay is not so large as to cause some stations to be significantly out-of-phase with respect to others in the array for sources which are not too far apart.
• Subsequent to bandpass calibration, residual delays need to be removed. To estimate "how good is good enough", assume that we don't want the delay to change by more than P degrees over the 8 GHz of bandwidth. In other words, the correction must be good enough such that phase gradients in frequency and in time over the bandpass are not introduced to the residual visibility data. Since one turn = 1 ns over 1 GHz, one turn = 125 ps of delay over 8 GHz. If we want P < 6 degrees (1/60th of a turn), then the residual delay should be less than 125/60 = 2 ps over 8 GHz at the zenith. Note that this corresponds to about 0.01 arcsec of refraction (delay*c/H0 =~ refraction, where H0 = height of troposphere).
• The delay correction made by CALC must be placed in the archival data base with sufficient accuracy and frequency so that the delay corrections applied can be interpolated for any visibility sample to an accuracy of less than 0.1 psec. For astrometric work, the total delay = calc model delay + residual delay + residual (phase/frequency) is the main observational parameter.
• Since CALC contains a rather general model for the atmospheric delay, it is not appropriate to use CALC to calculate the atmospheric delay. This can be done elsewhere (with the part of the code that calculates refraction, see "Atmospheric Refractive Signal Bending and Propagation Delay").
• Refraction calculations (as opposed to refractivity) are done elsewhere, not in CALC.
• The Lanyi (1984) model for the atmospheric delay is the most rigorous of the "mapping function" techniques for modeling delay.
• The input zenith atmospheric dry and wet delays should ultimately be calculated using ATM.
• Since the physical mechanisms which produce atmospheric delay are the same ones that produce atmospheric refraction, the same general formalism can be used for both. Note that the atmospheric refraction model proposed by Mangum (2001; ALMA Memo 366) is a variation on the continued fraction mapping function used in many of the geodesy calculations of atmospheric path delay.

## The ALMA Delay Server

The ALMA delay server is a software module which is responsible for distributing the delay correction for each antenna to the various bits of hardware in the antenna and to the correlator. The delay correction is defined as:

where is some arbitrary offset to make sure all delays are positive, is a per antenna constant representing the electronic delay, is the geometric delay derived by making a call to CALC, and is the atmospheric delay correction (which can be gotten from CALC but is (as of 2008/03/18) not used and set to zero. In the long term the atmospheric delay will be derived using ATM to calculate the zenith delay coupled to an appropriate mapping function.

## CALC

### Summary of VLA/VLBA Implementation by John Benson

The following description of CALC was derived from a description provided by John Benson:

#### History

The CALC program was written and is maintained by the VLBI group at GSFC. It is part of their off-line data reduction package for astrometric and geodetic VLBI observations. CALC is mainly paired up with a big least squares fit program (SOLVE) that solves for source positions, VLBI antenna positions, and other things. GSFC uses CALC to remodel their VLBI observations prior to running the SOLVE program. I think they run CALC-SOLVE iteratively to get the best solutions.

In 1985 we in the VLBA correlator group decided that the correlator on-line model should be a high accuracy model that would be good enough to allow phase referencing in the VLBA directly. That is, no delay model corrections would have to be made in post-processing. We chose CALC as the model to use. We also considered the JPL MODEST program. In hindsight this was a good choice, as CALC is a long-lived standard amongst astrometric model programs.

In 1985 we also recognized that besides a high accuracy on-line delay model, we needed accountability for that model. That is, the archived VLBA data needed tables that exactly described the model values that the correlator applied. If the user was going to correct his data in post-processing, the correction process needed to know what model values were used when by the correlator. For this reason the ASDM contains a set of tables that describe the delay model at various levels. They're based on the VLBA archive data model.

#### The CALC Model Philosophy

The CALC model in the VLBA correlator generates station delays with respect to the Earth's center, and CALC is run in a mode where a CALC baseline is stn A = earth center, stn B = antenna. GSFC added this mode for us. It basically turns off model components for stn A, like atmos, etc. In the VLBA correlator, the reference time (UTC) is the wave front arrival time at the earth center. That is, for a given antenna, the wave front delay from the antenna to the geocenter is the delay for that wave front. I like to think of a time tagged wave front as a pulse. The pulse hits the antenna and later hits the geocenter at time t. CALC calculates the delay. The correlator model generates antenna delays all at the same reference time, which means the delays are for the same wave front which is exactly what you want.

For the VLBA correlator in the early years we actually ran the CALC program in the on-line computer. Later on we put CALC in an off-line rpc server and that's what we do today. And in fact I often bring up a CALC server on other machines for testing, etc. So a client process sends an rpc request to a CALC server, which then receives a response containing delays, and other things. The CALC server request is for one antenna at one reference time. And the request is assembled in a simple C struct.

The CALC model itself contains individual model components where the GFSC philosophy was to calculate each model component to maximum accuracy. There is no error budget that governs the design of the model components. So you might say some items are 'over-modeled' since other items might not be as well determined and dominate the error level of the CALC delays. But CALC was part of a fitting solution, so I guess the intention was to solve for corrections to the poorly modeled components, like the atmosphere.

The CALC atmospheric model was built by Art Neill years ago, and it's based on radiosonde balloon flight data mainly from airports. It's certainly not adequate for ALMA, and can be turned off in CALC. (More information on the atmospheric model below).

The CALC subroutines are called from a wrapper that I wrote, calcmodl2.f. This includes calculations of delay corrections for curved wavefronts, and gravitational bending for objects within the solar system.

#### VLA/EVLA Refraction Calculation

The current refraction correction in the VLA is basically the same algorithm that has been in use for years (written by Barry Clark). It uses local pressure, temp, and dew point. It is not included in CALC or its wrapper, and is in fact part of the online control system.

#### CALC Server Include File

The following include file contains descriptions of the request and response values for CALC:

/*  @(#)CALCServer.h  version 1.00  created 99/07/26 10:00:00
%% include for the CALC Server RPC interface
LANGUAGE: C (include file)
ENVIRONMENT: vxWorks, Solaris, Linux
*/

#ifndef _CALC_H_RPCGEN
#define   _CALC_H_RPCGEN

#include <rpc/rpc.h>

#define MAX_STR 32

/* Client makes a CALC request through getCALC_arg arguments */
/* CALC Server Version 5 */

struct getCALC_arg {
long request_id;        /* RPC request id number, user's choice */
long date;              /* CALC model date (MJD) */
long ref_frame;         /* CALC reference frame: 0 = geocentric */
long dummy;
short int kflags[64];   /* CALC model component control flags */

double time;            /* CALC model time UTC (fraction of day) */

double a_x;             /* geocentric right-hand x coord (meters) */
double a_y;             /* geocentric right-hand y coord (meters) */
double a_z;             /* geocentric right-hand z coord (meters) */
double axis_off_a;      /* non-intersecting axis offset (meters) */
double b_x;             /* geocentric right-hand x coord (meters) */
double b_y;             /* geocentric right-hand y coord (meters) */
double b_z;             /* geocentric right-hand z coord (meters) */
double axis_off_b;      /* non-intersecting axis offset (meters) */
double ra;              /* J2000.0 coordinates (radians) */
double dec;             /* J2000.0 coordinates (radians) */
double dra;             /* J2000.0 arcsecs/year */
double ddec;            /* J2000.0 arcsecs/year */
double depoch;          /* reference date for which proper motion
* corrections are zero, mjd.fract_day */
double parallax;        /* source distance in arcsecs of annual
* parallax, = 206265.0 / distance (au) */

/* Earth Orientation Parameters */
double EOP_time[5];     /* EOP epoch date.time (MJD) */
double tai_utc[5];      /* TAI - UTC (secs) */
double ut1_utc[5];      /* UT1 - UTC (secs) */
double xpole[5];        /* earth pole offset, x (arcsec) */
double ypole[5];        /* earth pole offset, y (arcsecs) */

double pressure_a;      /* surface pressure stna (millibars)
* enter 0.0 for none availiable */
double pressure_b;      /* surface pressure stnb (millibars)
* enter 0.0 for none availiable */

char  *station_a;       /* station A name */
char  *axis_type_a;     /* station A mount type, 'altz', 'equa',
,xyns', 'xyew' */
char  *station_b;       /* station B name */
char  *axis_type_b;     /* station B mount type, 'altz', 'equa',
'xyns', 'xyew' */
char  *source;          /* source name */

};
typedef struct getCALC_arg getCALC_arg;

struct CALCRecord {
long request_id;        /* RPC request id number, returned by Server */
long date;              /* CALC model date (MJD) */
double time;            /* CALC model time UTC (fraction of day) */
double delay[2];        /* total group delay, rate (secs, sec/sec) */
double dry_atmos[4];    /* dry atmosphere delay, rate (secs, sec/sec) */
double wet_atmos[4];    /* wet atmosphere delay, rate (secs, sec/sec) */
double UV[3];           /* u, v, w coordinates in J2000.0 frame (meters) */
double riseset[2];      /* estimated rise, set times for stnb VLBA Correlator */
};
typedef struct CALCRecord CALCRecord;

struct getCALC_res {
int error;
union {
struct CALCRecord record;
char *errmsg;
} getCALC_res_u;
};
typedef struct getCALC_res getCALC_res;

#define   CALCPROG ((unsigned long)(0x20000340))
#define   CALCVERS ((unsigned long)(1))
#define   GETCALC ((unsigned long)(1))
extern  getCALC_res * getcalc_1();
extern int calcprog_1_freeresult();

/* the xdr functions */
extern bool_t xdr_getCALC_arg();
extern bool_t xdr_CALCRecord();
extern bool_t xdr_getCALC_res();

#endif /* !_CALC_H_RPCGEN */


## Atmospheric Delay Models

### References

• Niell, A. E. (1996), Journal of Geophysical Research, Vol. 101, No. B2, Pages 3227-3246: "Global Mapping Functions for the Atmospheric Delay at Radio Wavelengths". The standard reference for the derivation of a global mapping function for atmospheric delay. This derivation of the mapping function is really somewhat unique in that it attempts to analytically represent the global weather variations as a function of location (latitude) and time of year, and contains to adjustable parameters (i.e. does not require input pressure and temperature for each station). Note that Equation 4 in this paper has a typo whereby the terms which are printed as "1/term" in both the numerator and denominator should really be just "term" in both the numerator and denominator.
• Davis, J. L., Herring, T. A., Shapiro, I. I., Rogers, A. E. E., Elgered, G. (1985), Radio Science, Vol. 20, No. 6, Pages 1593-1607: "Geodesy by Radio Interferometry: Effects of Atmospheric Modeling Errors on Estimates of Baseline Length". An application of a modified Smith-Weintraub refractivity and the Niell mapping functions.
• Sovers, O. J., Fanselow, J. L., and Jacobs, C. S. (1998), Reviews of Modern Physics, Vol. 70, No. 4, Pages 1393-1454: "Astrometry and Geodesy with Radio Interferometry: Experiments, Models, Results". An excellent overview paper describing the details involved in calculating geometric and atmospheric delay. Uses the Lanyi (1984) model for the mapping function, which is a significant departure from the standard (i.e. Niell) mapping functions which derive from the Marini (1972) reduced fraction functional form.
• Lanyi, G. (1984), TDA Progress Report 42-78, Page 152-159: "Tropospheric Delay Effects in Radio Interferometry". Derivation of a new "tropospheric" (really atmospheric) mapping function which, unlike previous mapping functions, takes account of second and third order effects in the refractivity which are due to refractive bending. This derivation of the mapping function is really quite unique in that it does not fully separate the dry and wet contributions to the delay, making it a physically more exact representation. It is claimed to be more accurate than previous (i.e. Niell) mapping functions for E > 4 degrees, and the error due to the derived analytic form for the mapping function is estimated to be less than 0.02% for E > 6 degrees.
• Conway, J. (2004), ALMA Memo 503: "Antenna Position Determination: Observational Methods and Atmospheric Limits". Connection between baseline determination and delay calculations/models.
• Yan, H. and Ping, J. (1995), AJ, 110, 934: "The Generator Function Method of the Tropospheric Refraction Corrections". Another derivation of a new "tropospheric" (really atmospheric) mapping function. A cousin to existing reduced-fraction expansions of the mapping function.
• Yan, H. (1996), AJ, 112, 1312: "A New Expression for Astronomical Refraction". Related to the Yan and Ping (1995) reference above, but applied to the refraction calculation problem. Using the Yan and Ping (1995) and Yan (1996) references one can apply a unified formalism to both the atmosphere-induced refractive delay and bending problems.

### Overview

The delay experienced by an incoming signal due to its propagation through the Earth's atmosphere is given by:

where is the path through and n is the refractive index of the atmosphere. Since is very nearly unity for the Earth's atmosphere one normally uses the "refractivity" () instead of the index of refraction. Refractivity and refractive index are related as follows:

Furthermore, the atmospheric delay can be separated into contributions due to the dry and wet atmosphere:

where is the contribution due to dry air while is the contribution due to wet air. In general and are parameterized in terms of a zenith contribution to the delay which is dependent upon local atmospheric conditions () and a "mapping function" () which relate delays at an arbitrary elevation angle to that at the zenith:

Since the elevation angle is the unrefracted source elevation, refraction effects are included in the mapping functions . In the following I describe calculations of and .

### Zenith Delay

The contribution to the atmospheric delay at the zenith () is a measure of the integrated refractivity of the atmosphere at the zenith. In general, the refractivity of moist air at microwave frequencies depends upon the permanent and induced dipole moments of the molecular species that make up the atmosphere. The primary species that make up the dry atmosphere, nitrogen and oxygen, do not have permanent dipole moments, so contribute to the refractivity via their induced dipole moments. Water vapour does have a permanent dipole moment. Permanent dipole moments contribute to the refractivity as , while induced dipole moments contribute as , where is the pressure and is the temperature of the species.

A simple parameterization of the refractivity at the zenith is given by the Smith-Weintraub equation (Smith & Weintraub 1953):

where and are the partial pressures due to dry and wet air, T is the temperature of the atmosphere, and , , and are constants. The pressures and temperature are usually taken to be equivalent to the local ambient values at the antenna station under study. The dry and wet air refractivities are then given by:

The dry air contribution to this refractivity () is primarily due to oxygen and nitrogen, and is nearly in hydrostatic equilibrium. Therefore, does not depend upon the detailed behaviour of dry air pressure and temperature along the path through the atmosphere, and can be derived based on local atmospheric temperature and pressure measurements. The wet air refractivity () can be inferred from local water vapour radiometry measurements. Also, one can use an atmospheric model to derive the total atmospheric refractivity using an atmospheric model (such as ATM) with local atmospheric conditions as input.

### Mapping Function

The simplest form for the mapping function (), which relates the delay at an arbitrary elevation angle to that at the zenith, is given by the plane-parallel approximation for the Earth's atmosphere:

This simple form is in fact inadequate, which led Marini (1972) to consider corrections to this simple functional form which accounted for the Earth's curvature. Assuming an exponential atmospheric profile where the atmospheric refractivity varies exponentially with height above the antenna, Marini (1972) developed a continued fraction form for the mapping function:

where I include only the first three terms in the continued fraction. A slight modification to the Marini (1972) continued fraction functional form which forces at the zenith replaces the even numbered terms (i.e. the second, fourth, sixth, etc.) with

Chao (1974) introduced this modification by truncating the Marini (1972) form to include only two terms.

A more generalized continued-fractional form for the mapping function was developed by Yan and Ping (1995):

...where...

...is the "normalized effective zenith argument" of function which includes the "normalized effective height" of the atmosphere () defined as:

For an exponentially-distributed atmosphere:

Normally, (i.e. start the integration from the ground).

The constants , etc. in the continued fraction forms above are generally derived from analytic fits to ray-tracing results either for standard atmospheres or for observed atmospheric profiles based on radiosonde measurements. The mapping functions derived in Niell (1996) and Davis (1985) are derived in this way.

A physically more correct mapping function has been derived by Lanyi (1984). Unlike previous mapping functions, Lanyi does not fully separate the dry and wet contributions to the delay, which is a more physically correct approximation. It is based on an ideal model atmosphere whose temperature is constant from the surface to the inversion layer , then decreases linearly with height at rate (the lapse rate) from to the tropopause height , then is assumed to be constant above . This mapping function is designed then to be a semi-analytic approximation to the atmospheric delay integral that retains an explicit temperature profile that can be determined using meteorological measurements. The mapping function is expanded as a second-order polynomial in and , plus the largest third-order term). It is nonlinear in and . It also contains terms which couple and , thus including terms which arise from the bending of the signal path through the atmosphere. The functional form for the atmospheric delay in this Lanyi (1984) model is given by:

where

where

With standard values of , , (appropriate for mid-latitudes), and , .

The dry, wet, and bending contributions are expressed in terms of moments of the refractivity. The bending terms are evaluated for the ideal model atmosphere and thus give the dependence of the delay on the four parameters , , , and . Therefore, the Lanyi (1984) model relies upon accurate surface meteorological measurements at the time of the observations to which the delay model is applied.

#### Mapping Function Summary

• Differences between the various mapping functions only show up at low elevation (< 10 degrees). Since geodesists do observe at very low elevations, these differences can be significant. This is not the case for astronomers.
• The Yan and Ping (1995) form has a "cousin" used for refractive bending ("refraction correction") calculations, making it a convenient choice for both refractive delay and bending calculations.

### Antenna Height Correction to Total Atmospheric Delay

In the calculation of the zenith atmospheric delay at an antenna it is assumed that the atmospheric properties (P, T, RH) are the values measured at the focal plane of the antenna. For example, in VLBI each station has a set of associated weather measurements which are used to calculate . For a clustered array like the VLA or ALMA, the affect of the differences in antenna focal plane height above some reference point need to be accounted for.

For the VLA (not EVLA), CALC was not used to calculate the atmospheric delay. The antenna height correction was incorporated with a simple atmospheric delay correction by correcting for the path difference between each VLA antenna and a reference point at the center of the array. For the VLA case, the extra atmospheric path due to a difference in antenna height above the center-of-the-array reference point (, in ns) is given by:

where is the atmospheric refractivity, is the atmospheric zenith delay calculated using the VLA weather station (which is located near the center of the array), is the geometric w of the antenna (in ns). The first term is the antenna height correction to the zenith delay, while the second term is a simple atmospheric delay correction. For EVLA, CALC will be used to calculate both geometric and atmospheric delay. We believe (though have not confirmed) that CALC also calculates the antenna height correction (first term in the equation above) given antenna heights relative to the reference point at the center of the array. ALMA will need to include this antenna height correction term.

A simple estimate of the magnitude of the antenna height difference correction at the zenith can be gotten by assuming that the pressure P changes linearly with height. Then 52 cm of additional antenna height (the current difference in height between the two ATF antennas) out of a total atmospheric height of 8 km would correspond to:

where I have assumed P = 1053 mb. The dry term zenith atmospheric delay changes approximately like 2.3 mm/mb of pressure change. A pressure change of 0.068 mb corresponds to approximately 156 micron of path difference. This is consistent with alternate back-of-the-envelope calculations of this quantity.

### Differential Excess Atmospheric Delay Between Two Antennas

NOTE: The following is just an aside. Since CALC or any other analysis of the atmospheric delay at an antenna calculates the total integrated delay along the path of observation, the differential delay between two antennas is accounted for in any differencing calculations done during baseline determination.

The differential delay induced in an interferometer by a horizontally stratified troposphere results from the difference in zenith angle of the source at the antennas. Thompson, Moran, and Swensen (2001), pp. 516-518 discuss the atmospheric delay induced along an interferometer baseline. The excess path length is given by:

where is the refractivity at the Earth's surface, is the height above the Earth's surface, is the atmospheric scale height, is the length coordinate along the direction to the source, and E is the antenna elevation while observing the source. Note that refraction is neglected. One can relate , , , and as follows (see Figure 13.4 in TMS, page 517; reproduced below) using the cosine rule on the triangle formed by , , and :

Solving for and using elevation rather than zenith angle yields:

For the triangle which is formed by sides , , and the side which is equal to , we can write:

Since and (the height of the troposphere), . Since , . The equation for in terms of , , and then becomes:

(Thanks to Dick Thompson for filling-in some of the details of this calculation).

We can now write the expression for as follows:

Since , the second term in the equation above can be expanded with a Taylor series so that:

Integration yields:

Writing this equation in terms involving , the excess path length becomes:

Taking the derivative of with respect to E and multiplying this derivative by the baseline length divided by yields the atmospheric differential delay between two antennas separated by baseline :

...where is in m, is in km, and is in km. Note that once must calculate using a suitable atmospheric model which uses measurements of the local atmospheric pressure, temperature and relative humidity to derive the resultant differential residual delay.

I have produced a plot of this relation as a function of for baseline lengths running from 10 to 100,000 m and elevation 1 to 90 degrees. The contour levels on this plot are 10, 15, 20, 25, 50, 75, 100, 200, 400, 600, 800, 1000, 3000, and 5000 micron; right to left.

Zooming in a bit on the vertical (baseline length) axis of this plot, I have produced two zoomed versions of the above figure. One for baselines up to 1 km (contour levels: 0.25, 0.50, 0.75, 1.0, 2.5, 5.0, 7.5, 10, 25, 50, 75, 250, 500, 750, 1000 micron; right to left),

and a second plot for baselines up to 100 m (contour levels: 0.05, 0.075, 0.1, 0.25, 0.5, 0.75, 1.0, 2.5, 5.0, 7.5, 10, 25, 50 micron; right to left).
Topic revision: r28 - 2009-08-16, JeffMangum
Copyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding NRAO Public Wiki? Send feedback