Description and use of the POLynomial ANalysis subroutine


for the calculation of real-height profiles from sweep-frequency ionograms.

This outline is for the Version of September 1986. It supplements the information contained in the report Ionogram analysis with the generalised program POLAN, obtainable as report UAG-93

World Data Center A, NOAA E/GC2, 325 Broadway, Boulder, Colorado 80303

If problems arise, run one data set with LIST = 3, and mail all output to:

J.E. Titheridge, Physics Dept., University of Auckland, New Zealand.

NEW TO THIS VERSION (September 1986)


POLAN is called with frequency, height data in the arrays FV, HT.

In the data arrays, intermediate layers are terminated by a scaled (or zero) value for the critical frequency, with:-

The final layer is terminated by at least 2 null points, with h=f=0.
Data can be terminated without a peak by using a final frequency of -1.0.

Data for the extraordinary ray, if any, precedes the o-ray data for each layer. This is because x-ray data is used only to calculate the (start or valley) corrections to be made at the beginning of the calculation for that layer.

The format for input data is best seen by study of the examples in the test file POLRUN.DAT.


FB gives the gyrofrequency at the ground in MHz, for an inverse cube variation. If you have only the gyrofrequency FH at a height h km, the ground value is obtained from FB = FH * (1. + h/6371.2)**3.
To use a gyrofrequency (FH, say) which is independent of height, set FB = -FH.

DIP is the magnetic dip angle IN DEGREES. Use of a negative value for DIP suppresses the physical checks which are normally applied to the calculated profile, so that the result obtained is the best mathematical (but possibly non-physical) fit to the virtual-height data. [Some physically based equations are still included in start and valley calculations, unless AMODE is negative.]

START normally gives a model height at 0.5 MHz. Typical values are:

noon   sunset-2/rise+2hr   set/rise    set+1hr   set+2    set+4 to rise-1
85km 88km(E layer) 90(E)/80(F) 100 km 130 km 150 km.

A preferred procedure is to calculate model values of START from the equations given in J. Atmosph. Terr. Phys. 48, 435-446, 1986.

Use of START = 0.0 makes some allowance for underlying ionisation based on a limited extrapolation of the first few virtual heights.

With initial x-ray data, START is taken to give the gyrofrequency height for underlying ionisation calculations; the values listed above are still suitable for this purpose. The x-ray data is used to calculate a slab start correction from 0.3*fmin (adding points at 0.3, 0.6 and 0.8 *fmin).

[Alternative procedures can be obtained using non-standard values of START:-
START between 0 and 44 defines the plasma frequency for a model start.
Start = -1.0 uses a direct start, from the first scaled point.
Start < -1.0 for x-starts to use a polynomial from (-Start -1.0) MHz. ]

THE final three parameters - AMODE, VALLEY and LIST, are zero for most work.

AMODE sets the type of analysis, as listed below. Zero uses mode 6. Use AMODE+10. for 12-point integrals, for high accuracy at large dip angles (this is done automatically, at DIP > 60, when AMODE=0).

Values of AMODE greater than 29.0 are used to specify the number of polynomial constants to be used to describe each ionospheric layer; thus

Setting AMODE negative causes physical relations to be omitted from the start and valley calculations.

VALLEY= 0.0 or 1.0 uses a valley width equal to the initial default value of twice the local scale height. The initial default depth is 0.05 MHz. The calculated depth is scaled according to (calculated width)**2.

Alternative solutions may be obtained as follows:


LIST negative suppresses most trace output below the first peak. LIST= -10 suppresses all output, even the normal layer summaries.


The arrays FV, HT contain the calculated frequencies and real heights.

N gives the number of calculated real-height data points.

The peak of the last layer is at FC = fv(N-3), Hmax = ht(N-3).
Points at N-2, N-1 and N in the output arrays are extrapolated heights at 0.5, 1.0 and 1.5 scale heights above the peak (calculated using the Chapman expression with a scale height gradient of 0.2).

QQ returns the real-height coefficients, for single-polynomial calculations, as described under Section 4.2 below.



MODE is obtained from the input parameter AMODE, modified to the range 1 to 10, and is used to select the type of analysis as summarised below. All Modes include an estimated start correction, a Chapman-layer peak, and a model valley between layers.


The basic parameters which define the type of analysis depend on the parameter MODE, and are obtained from the arrays given below.

       |-------- First step ------|    |------- Following steps --| 
MODE= 1, 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
NT = 1, 2, 3, 4, 4, 5, 6, 6, 7,73, 1, 2, 3, 4, 5, 6, 6, 6, 7,73
NV = 1, 2, 3, 4, 5, 7, 8,10,12,40, 1, 1, 2, 3, 4, 5, 7, 8,13,40
NR = 0, 0, 0, 1, 1, 2, 2, 3, 5, 2, 0,-1,-1, 1,-2,-3,-3,-4,-6,-3
NH = 1, 1, 2, 3, 3, 4, 5, 6, 8,28, 1, 1, 1, 1, 1, 1, 2, 2, 3,28


These use a defined number of real-height coefficients for each layer, and return all profile parameters in the array QQ.The order of the analysis is set by the parameter AMODE, as follows.

QQ returns the real-height parameters which describe the profile, for single-polynomial modes of analysis (unless QQ(1) was set equal to -1.0 by the calling program).

The returned value of QQ(1) gives the total number of stored values (mq). Starting at QQ(2), the parameters returned for each layer are:
FA, HA, nq, q1, q2, .. qn, devn, FM, FC, HM, and SH.

nq is the number of polynomial coefficients (q1 to qn) used for this layer. This is normally equal to the number of coefficients requested in AMODE.

HA is the true height at FA, after any start or valley adjustments, so the real-height profile is
h = HA + q1.(f-FA) + q2.(f-FA)^2 + ... qn.(f-FA)^nq.

devn is the rms deviation (in km) of the fit to the virtual height data.

FC, HM andSH are the constants which define the Chapman-layer peak; this joins the polynomial section at the frequency FM (the highest scaled frequency for the layer).

For a 2nd (or 3rd) layer, FA, HA give the new real-height origin at the top of the valley region. Thus FA is equal to the previous FC, and the valley width is W = HA - HM in km. The valley depth D, (in MHz) can be obtained from the width using equations (14) of the report UAG-93, which give

D = 0.008 W**2/(20 + W) MHz, followed by D = D.FC/(D + FC).

The end point of the data in QQ is verified by a value QQ(mq+1) = -99 for a normal exit, and -98. for an error exit.



Analysis can proceed with a minimum of 2 scaled virtual heights (or 1 height and a critical frequency) for each layer. If the number of data points NV is less than the number of polynomial terms NT (as specified by AMODE), NT is automatically decreased.

Calculate one polynomial, with NT terms, from the point FA = fv(K) ,HA = ht(K), to fit the next NV virtual and NR real heights. (The fitted real heights include one point below HA, if NR is negative.)

The real-height origin (FA,HA) is at K = KR, in the data arrays FV, HT. The corresponding virtual height is at K = KV.

With x-ray data (-ve frequencies), at the start or after a peak, recalculate HA to include the correction for underlying or valley ionisation. Calculate a further NH real heights, and set KR = KR + NH; KV = KV + NH.

Repeat this loop, calculating successive overlapping real-height sections, until a critical frequency (or end-of-layer) is found in the range KV+1 to KV+NV+1. Then calculate real heights at the remaining scaled frequencies and determine a least-squares Chapman-layer peak.


(numbered according to the corresponding section in the program POLAN.) SECTION 2.2 Count initial x-rays. Check frequency sequencing. Check for cusp, peak, or end of data.


SECTION 2.3 Subtract the group retardation due to the last calculated real-height section. This modifies all the virtual heights at f > FA (where FA = fv(KR)), and increases the index LK (which gives the point up to which the group retardation has been removed) to KR.

SECTION 3. Set up equations for the next profile step. Check for the occurrence of a valley; if this is required, set the valley flag HVAL and set initial values for the width and depth. Set up equations in the matrix B. For start calculations using x-ray data, or for any valley calculations, add suitably weighted equations specifying desired physical properties of the solution.

SECTION 4. Solve the set of simultaneous equations in the array B.

Check that the solution satisfies basic physical constraints. If it does not, obtain a new least-squares solution with the limiting constraints imposed (in the subroutine ADJUST).

For an x-start or valley calculation, iterate the solution as required to ensure the use of a correct gyrofrequency height, and the correct relation between depth and width of the valley. For an o-ray valley, loop once to adjust the valley depth.

SECTION 5. Calculate and store the real heights. Set KRM as the index for the highest calculated real height.

SECTION 6. Least-squares fitting of a Chapman layer peak.

Calculate the critical frequency and the scale height of a layer peak, by an iterative fit to the real-height gradients at the last few calculated points (as in Radio Science 20, 247, 1985). Determine the height of the peak by fitting the peak shape to a weighted mean of the last few calculated real heights. Adjust the last real height to agree exactly with the Chapman peak (sept'86).

SECTION 7. Go to section 2, to restart for a new layer.

If there are no further data, extrapolate 3 points for the topside ionosphere (assuming a Chapman layer with a scale height gradient of 0.2 km/km), store constants relating to the last layer peak, and return.

J.E. Titheridge, November 1987.

True height analysis | Ionosondes | WDC

4-DEC-1997 Chris Davis