More Sanity Checking of the LANL GPS Charged-Particle Dataset

Prior posts in this series may be found at:

• Original data
• Stages 1 to 4
• Revised stages 1 to 5  
• Stages 6 to 10 
• Stages 11 to 15
• Stages 16 to 20  
• Stages 21 to 25
• Stages 26 to 30 

In this post I continue looking at the values of fields for the CXD experiment carried on ns53 through ns73.

Sanity check for L_shell

The values for L_shell include some NA values (created during earlier processing), representing locations for which the L-shell calculation is not meaningful. It is fairly debatable whether these should be removed, so as to ensure that the user cannot make the mistake of failing to check for NA values. On balance, though, it seems better to leave these values in the dataset, since the actual data from the satellite are not tainted.

Ignoring the NA values, the values of L_shell lie in the following ranges:

L_shell
Satellite	Minimum	Maximum
53	4.079050e+00	6.885292e+01
54	4.073716e+00	7.250067e+01
55	4.067056e+00	7.125280e+01
56	4.081161e+00	7.127625e+01
57	4.113699e+00	7.339611e+01
58	4.088591e+00	7.147417e+01
59	4.058853e+00	7.154358e+01
60	4.083813e+00	7.114283e+01
61	4.061549e+00	7.002652e+01
62	4.082852e+00	6.885380e+01
63	4.076490e+00	7.115145e+01
64	4.105269e+00	7.008071e+01
65	4.068625e+00	6.727069e+01
66	4.084220e+00	6.957386e+01
67	4.108083e+00	7.065966e+01
68	4.108754e+00	7.211983e+01
69	4.106989e+00	6.798032e+01
70	4.100414e+00	6.645398e+01
71	4.103767e+00	7.137340e+01
72	4.103526e+00	7.117196e+01
73	4.096220e+00	6.786675e+01

These look reasonable.

Sanity check for L_LGM_TS04IGRF

The values for L_LGM_TS04IGRF include some NA values (created during earlier processing), representing locations for which the L-shell calculation is not meaningful. It is fairly debatable whether these should be removed, so as to ensure that the user cannot make the mistake of failing to check for NA values. On balance, though, it seems better to leave these values in the dataset, since the actual data from the satellite are not tainted.

Ignoring the NA values, the values of L_LGM_TS04IGRF lie in the following ranges:

L_LGM_TS04IGRF
Satellite	Minimum	Maximum
53	3.947573e+00	6.927239e+01
54	3.901199e+00	6.919588e+01
55	3.936277e+00	6.879617e+01
56	3.792127e+00	7.015551e+01
57	4.008310e+00	7.110403e+01
58	3.900022e+00	7.005057e+01
59	4.007891e+00	6.847356e+01
60	4.003237e+00	7.221038e+01
61	3.902738e+00	6.978320e+01
62	3.912530e+00	7.011196e+01
63	3.971503e+00	6.971264e+01
64	4.024707e+00	6.769095e+01
65	4.015073e+00	7.188151e+01
66	4.007551e+00	6.856100e+01
67	4.050529e+00	6.987180e+01
68	3.995247e+00	6.942418e+01
69	4.048874e+00	6.902460e+01
70	4.046164e+00	6.815424e+01
71	4.071240e+00	6.757528e+01
72	3.999839e+00	6.450151e+01
73	4.044781e+00	6.871470e+01

These look reasonable.

Sanity check for L_LGM_OP77IGRF

The values for L_LGM_OP77IGRF include some NA values (created during earlier processing), representing locations for which the L-shell calculation is not meaningful. It is fairly debatable whether these should be removed, so as to ensure that the user cannot make the mistake of failing to check for NA values. On balance, though, it seems better to leave these values in the dataset, since the actual data from the satellite are not tainted. I note that the description of this field claims that it is "not currently filled". This is obviously not true.

Ignoring the NA values, the values of L_LGM_OP77IGRF lie in the following ranges:

L_LGM_OP77IGRF
Satellite	Minimum	Maximum
53	4.113965e+00	3.660894e+01
54	4.106457e+00	3.617539e+01
55	4.100884e+00	3.591780e+01
56	4.113720e+00	3.660868e+01
57	4.146636e+00	3.850798e+01
58	4.120502e+00	3.787992e+01
59	4.091215e+00	3.613692e+01
60	4.117398e+00	3.684020e+01
61	4.094553e+00	3.684505e+01
62	4.114965e+00	3.617936e+01
63	4.108885e+00	3.639508e+01
64	4.138934e+00	3.603204e+01
65	4.103389e+00	3.535154e+01
66	4.116871e+00	3.776598e+01
67	4.144117e+00	3.279393e+01
68	4.145199e+00	2.837767e+01
69	4.138738e+00	3.545018e+01
70	4.132139e+00	3.260918e+01
71	4.136635e+00	3.643418e+01
72	4.137440e+00	3.532261e+01
73	4.127833e+00	3.632726e+01

These look (mostly) reasonable; but note that the maximum valid values are considerably less than for the other magnetic field models, and also the low maximum for ns68, which may indicate an issue worthy of more investigation with the data for that satellite.

Sanity check for L_LGM_T89CDIP

The values for L_LGM_T89CDIP include some NA values (created during earlier processing), representing locations for which the L-shell calculation is not meaningful. It is fairly debatable whether these should be removed, so as to ensure that the user cannot make the mistake of failing to check for NA values. On balance, though, it seems better to leave these values in the dataset, since the actual data from the satellite are not tainted.

Ignoring the NA values, the values of L_LGM_T89CDIP lie in the following ranges:

L_LGM_T89CDIP
Satellite	Minimum	Maximum
53	4.156665e+00	7.046841e+01
54	4.134309e+00	7.297807e+01
55	4.141375e+00	7.265581e+01
56	4.147715e+00	7.362871e+01
57	4.171050e+00	7.021322e+01
58	4.160699e+00	6.858366e+01
59	4.138150e+00	7.196883e+01
60	4.132394e+00	7.063522e+01
61	4.135012e+00	7.077780e+01
62	4.164926e+00	6.988359e+01
63	4.155301e+00	7.192335e+01
64	4.165912e+00	6.925883e+01
65	4.153215e+00	6.977235e+01
66	4.163230e+00	6.761487e+01
67	4.174354e+00	6.961524e+01
68	4.172899e+00	6.893099e+01
69	4.178059e+00	7.078880e+01
70	4.174908e+00	7.211150e+01
71	4.179994e+00	7.007543e+01
72	4.171566e+00	7.102184e+01
73	4.177882e+00	6.814721e+01

These look reasonable.

Sanity check for bfield_ratio

The definition of bfield_ratio is:

Column	Variable name	type	Dim.	description
33	bfield_ratio	double	1	Bsatellite/Bequator

 This appears to mean:

Column	Variable name	type	Dim.	description
33	bfield_ratio	double	1	Ratio of magnetic field at the satellite to the magnetic field along the field line to the equator.

The recorded value of this field is not obvious (especially since it is naturally highly dependent on the field model one uses).

Looking at the values of the bfield_ratio, one quickly sees many NA values and, worse, values of the form -n.nnnnnne+99 that make no sense at all. It is therefore not obvious that this field was ever subject to any kind of quality control.

In private communication with the LANL team, I discovered that:

1. There is no in situ satellite instrument for measuring magnetic fields;
2. All fields related to magnetic measurement are calculated, not measured.
3. The primary reference for the calculation is: https://github.com/drsteve/LANLGeoMag.

It does not, therefore, seem sensible to retain these fields from the records (i.e., in particular, bfield_ratio, b_sattelite [which I remove with some relief, for obvious reasons] and b_equator [removed with relief for reasons that are almost as obvious if one looks at the official documentation for the field]. If a user wishes to calculate the values of equivalent fields for himself, he may do so, using the then-current versions of the routines cited above or similar routines. (Although it is obviously troubling that apparently it was possible to obtain nonsense results from the published routines, at least as of the date when the original 1.3 dataset was released. So, if a user calculates magnetic values for himself, he is urged to check their sanity before relying on them.)

Stage 31: Remove bfield_ratio

for file in ns[567]*
  do 
    awk '{$33=""; print $0}' $file | tr -s " " | sed 's/ $//' > ../gps-stage-31/$file
  done

The records for ns53 to ns73 now look like this:

Column	Variable name	type	Dim.	description
1	decimal_day	double	1	GPS time: a decimal number in the range [1, 367) in leap years or [1, 366) otherwise, representing the day of the year (1-Jan 00:00 to 31-Dec 24:00).
2	Geographic_Latitude	double	1	Latitude of satellite (°, N +ve)
3	Geographic_Longitude	double	1	Longitude of satellite (°, E +ve, measured from Greenwich meridian)
4	Rad_Re	double	1	Distance from centre of Earth, in units of Earth radii.
5-15	rate_electron_measured	double	11	Measured rate (Hz) in each of the 11 CXD electron channels
16-20	rate_proton_measured	double	5	Measured rate (Hz) in each of the 5 CXD proton channels (P1-P5)
21	LEP_thresh	double	1	LEP threshold in E1 channels, in keV
22	collection_interval	int	1	dosimeter collection period (seconds)
23	year	int	1	year (e.g. 2015)
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0) / (days in year)
25	SVN	int	1	Satellite Vehicle Number
26	b_coord_radius	double	1	Distance from dipole axis, in units of Earth radii.
27	b_coord_height	double	1	Distance from dipole equatorial plane, in units of Earth radii (N +ve).
28	magnetic_longitude	double	1	Magnetic longitude (degrees)
29	L_shell	double	1	L shell: McIlwain calculation according to model with T89 External Field, IGRF Internal Field.
30	L_LGM_TS04IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, TS04 External Field, IGRF Internal Field.
31	L_LGM_OP77IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, OP77 External Field, IGRF Internal Field (not currently filled)
32	L_LGM_T89CDIP	double	1	LanlGeoMag L-shell McIlwain calculation, T89 External Field, Centered Dipole Internal Field
33	local_time	double	1	magnetic local time (0-24 hours)
34	utc_lgm	double	1	UTC (0-24 hours)
35	b_sattelite	double	1	B field at satellite (gauss)
36	b_equator	double	1	B field at equator (on this field line I think) (gauss)
37-47	electron_background	double	11	estimated background in electron channels E1-E11 (Hz)
48-52	proton_background	double	5	estimated background in proton channels P1-P5 (Hz)
53	proton_activity	int	1	=1 if there is significant proton activity
54	proton_temperature_fit	double	1	characteristic momentum -- R0 in the expression given above (MeV/c)
55	proton_density_fit	double	1	N0 parameter in fit to proton flux ((protons/(cm2 sec sr MeV))
56	electron_temperature_fit	double	1	electron temperature from a one Maxwellian fit (MeV)
57	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
58-62	model_counts_proton_fit_pf	double	5	P1-P5 rate from proton fit (using proton_temperature_fit, proton_density_fit)
63-73	model_counts_electron_fit	double	11	E1-E11 rates from the 9-parameter electron flux model
74-79	proton_integrated_flux_fit	double	6	integral of proton flux (based on fit) above 10, 15.85, 25.11, 30, 40, 79.43 MeV (proton kinetic energy)
80-109	integral_flux_instrument	double	30	(based on 9 parameter fit) integral of electron flux above integral_flux_energy[i] particles/(cm2 sec)
110-139	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
140-154	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
155-169	electron_diff_flux	double	15	(based on 9 parameter fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))
170-178	Efitpars	double	9	fit parameters for 9 parameter electron fit

Stage 32: Remove b_sattelite

  for file in ns[567]*
  do 
    awk '{$35=""; print $0}' $file | tr -s " " | sed 's/ $//' > ../gps-stage-32/$file
  done

The records for ns53 to ns73 now look like this:

Column	Variable name	type	Dim.	description
1	decimal_day	double	1	GPS time: a decimal number in the range [1, 367) in leap years or [1, 366) otherwise, representing the day of the year (1-Jan 00:00 to 31-Dec 24:00).
2	Geographic_Latitude	double	1	Latitude of satellite (°, N +ve)
3	Geographic_Longitude	double	1	Longitude of satellite (°, E +ve, measured from Greenwich meridian)
4	Rad_Re	double	1	Distance from centre of Earth, in units of Earth radii.
5-15	rate_electron_measured	double	11	Measured rate (Hz) in each of the 11 CXD electron channels
16-20	rate_proton_measured	double	5	Measured rate (Hz) in each of the 5 CXD proton channels (P1-P5)
21	LEP_thresh	double	1	LEP threshold in E1 channels, in keV
22	collection_interval	int	1	dosimeter collection period (seconds)
23	year	int	1	year (e.g. 2015)
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0) / (days in year)
25	SVN	int	1	Satellite Vehicle Number
26	b_coord_radius	double	1	Distance from dipole axis, in units of Earth radii.
27	b_coord_height	double	1	Distance from dipole equatorial plane, in units of Earth radii (N +ve).
28	magnetic_longitude	double	1	Magnetic longitude (degrees)
29	L_shell	double	1	L shell: McIlwain calculation according to model with T89 External Field, IGRF Internal Field.
30	L_LGM_TS04IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, TS04 External Field, IGRF Internal Field.
31	L_LGM_OP77IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, OP77 External Field, IGRF Internal Field (not currently filled)
32	L_LGM_T89CDIP	double	1	LanlGeoMag L-shell McIlwain calculation, T89 External Field, Centered Dipole Internal Field
33	local_time	double	1	magnetic local time (0-24 hours)
34	utc_lgm	double	1	UTC (0-24 hours)
35	b_equator	double	1	B field at equator (on this field line I think) (gauss)
36-46	electron_background	double	11	estimated background in electron channels E1-E11 (Hz)
47-51	proton_background	double	5	estimated background in proton channels P1-P5 (Hz)
52	proton_activity	int	1	=1 if there is significant proton activity
53	proton_temperature_fit	double	1	characteristic momentum -- R0 in the expression given above (MeV/c)
54	proton_density_fit	double	1	N0 parameter in fit to proton flux ((protons/(cm2 sec sr MeV))
55	electron_temperature_fit	double	1	electron temperature from a one Maxwellian fit (MeV)
56	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
57-61	model_counts_proton_fit_pf	double	5	P1-P5 rate from proton fit (using proton_temperature_fit, proton_density_fit)
62-72	model_counts_electron_fit	double	11	E1-E11 rates from the 9-parameter electron flux model
73-78	proton_integrated_flux_fit	double	6	integral of proton flux (based on fit) above 10, 15.85, 25.11, 30, 40, 79.43 MeV (proton kinetic energy)
79-108	integral_flux_instrument	double	30	(based on 9 parameter fit) integral of electron flux above integral_flux_energy[i] particles/(cm2 sec)
109-138	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
139-153	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
154-168	electron_diff_flux	double	15	(based on 9 parameter fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))
169-177	Efitpars	double	9	fit parameters for 9 parameter electron fit

Stage 33: Remove b_equator

  for file in ns[567]*
  do 
    awk '{$35=""; print $0}' $file | tr -s " " | sed 's/ $//' > ../gps-stage-32/$file
  done

The records for ns53 to ns73 now look like this:

Column	Variable name	type	Dim.	description
1	decimal_day	double	1	GPS time: a decimal number in the range [1, 367) in leap years or [1, 366) otherwise, representing the day of the year (1-Jan 00:00 to 31-Dec 24:00).
2	Geographic_Latitude	double	1	Latitude of satellite (°, N +ve)
3	Geographic_Longitude	double	1	Longitude of satellite (°, E +ve, measured from Greenwich meridian)
4	Rad_Re	double	1	Distance from centre of Earth, in units of Earth radii.
5-15	rate_electron_measured	double	11	Measured rate (Hz) in each of the 11 CXD electron channels
16-20	rate_proton_measured	double	5	Measured rate (Hz) in each of the 5 CXD proton channels (P1-P5)
21	LEP_thresh	double	1	LEP threshold in E1 channels, in keV
22	collection_interval	int	1	dosimeter collection period (seconds)
23	year	int	1	year (e.g. 2015)
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0) / (days in year)
25	SVN	int	1	Satellite Vehicle Number
26	b_coord_radius	double	1	Distance from dipole axis, in units of Earth radii.
27	b_coord_height	double	1	Distance from dipole equatorial plane, in units of Earth radii (N +ve).
28	magnetic_longitude	double	1	Magnetic longitude (degrees)
29	L_shell	double	1	L shell: McIlwain calculation according to model with T89 External Field, IGRF Internal Field.
30	L_LGM_TS04IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, TS04 External Field, IGRF Internal Field.
31	L_LGM_OP77IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, OP77 External Field, IGRF Internal Field (not currently filled)
32	L_LGM_T89CDIP	double	1	LanlGeoMag L-shell McIlwain calculation, T89 External Field, Centered Dipole Internal Field
33	local_time	double	1	magnetic local time (0-24 hours)
34	utc_lgm	double	1	UTC (0-24 hours)
35-45	electron_background	double	11	estimated background in electron channels E1-E11 (Hz)
46-50	proton_background	double	5	estimated background in proton channels P1-P5 (Hz)
51	proton_activity	int	1	=1 if there is significant proton activity
52	proton_temperature_fit	double	1	characteristic momentum -- R0 in the expression given above (MeV/c)
53	proton_density_fit	double	1	N0 parameter in fit to proton flux ((protons/(cm2 sec sr MeV))
54	electron_temperature_fit	double	1	electron temperature from a one Maxwellian fit (MeV)
55	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
56-60	model_counts_proton_fit_pf	double	5	P1-P5 rate from proton fit (using proton_temperature_fit, proton_density_fit)
61-71	model_counts_electron_fit	double	11	E1-E11 rates from the 9-parameter electron flux model
72-77	proton_integrated_flux_fit	double	6	integral of proton flux (based on fit) above 10, 15.85, 25.11, 30, 40, 79.43 MeV (proton kinetic energy)
78-107	integral_flux_instrument	double	30	(based on 9 parameter fit) integral of electron flux above integral_flux_energy[i] particles/(cm2 sec)
108-137	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
138-152	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
153-167	electron_diff_flux	double	15	(based on 9 parameter fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))
168-176	Efitpars	double	9	fit parameters for 9 parameter electron fit

Sanity check for local_time

The definition of local_time is:

Column	Variable name	type	Dim.	description
33	local_time	double	1	magnetic local time (0-24 hours)

 There are a couple of problems with this definition:

1. The field name gives no indication that refers to a magnetic value, so should be renamed magnetic_local_time;
2. It provides no indication of exactly how magnetic local time is defined (and several mutually exclusive definitions are extant).

However, providing that one does not ascribe great accuracy to this value, it might conceivably be useful, so it is not unreasonable to retain it, with some emendations.

Firstly, we redefine it slightly (although unfortunately being forced to retain the ambiguity as to its precise meaning):

Column	Variable name	type	Dim.	description
33	magnetic_local_time	double	1	magnetic local time [0-24) hours

The actual values of this field lie in the following ranges:

magnetic_local_time
Satellite	Minimum	Maximum
53	7.845492e-06	2.399998e+01
54	7.956284e-06	2.400000e+01
55	3.351804e-05	2.399994e+01
56	1.586655e-05	2.399998e+01
57	4.068317e-05	2.399998e+01
58	2.519977e-05	2.399999e+01
59	4.977282e-06	2.400000e+01
60	7.915747e-06	2.399999e+01
61	2.649532e-06	2.400000e+01
62	1.991213e-05	2.400000e+01
63	3.556485e-06	2.400000e+01
64	1.172278e-05	2.399996e+01
65	1.809704e-05	2.399994e+01
66	4.205324e-05	2.399998e+01
67	4.227346e-05	2.399994e+01
68	2.200456e-05	2.399982e+01
69	3.367094e-05	2.399992e+01
70	1.883996e-05	2.399991e+01
71	1.952489e-04	2.399984e+01
72	1.620355e-04	2.399977e+01
73	6.868977e-05	2.399989e+01

There are obvious issues with these values:

1. They are not normalised to [0-24)
2. The numbers are reported to different precision, depending on their magnitude (and, for small numbers, they are reported to an unjustifiable precision).

These matters are easily fixed by the following script, do-stage-34.py:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

field_nr = 32    # magnetic_local_time (wrt 0)

with open(filename) as records:
  for line in records:
    fields = line.split()
    value_str = fields[field_nr]
    value = float(value_str)
  
    if (value >= 24):
      value = value - 24
    
    output_format = '{:8.5f}'
    
    new_value_str = output_format.format(value)
  
    fields[field_nr] = new_value_str
    newline = " ".join(fields)
 
    print newline

Stage 34: Reformat values of magnetic_local_time

Execute do-stage-34.py on all the relevant satellites:

for file in ns[567]*
do  
  ./do-stage-34.py $file > ../gps-stage-34/$file
done

The maxima and minima now have these (reasonable) values:

magnetic_local_time
Satellite	Minimum	Maximum
53	0.00001	23.99998
54	0.00000	23.99998
55	0.00003	23.99994
56	0.00002	23.99998
57	0.00004	23.99998
58	0.00003	23.99999
59	0.00000	23.99999
60	0.00001	23.99999
61	0.00000	23.99999
62	0.00000	23.99996
63	0.00000	23.99999
64	0.00001	23.99996
65	0.00002	23.99994
66	0.00004	23.99998
67	0.00004	23.99994
68	0.00002	23.99982
69	0.00003	23.99992
70	0.00002	23.99991
71	0.00020	23.99984
72	0.00016	23.99977
73	0.00007	23.99989

Sanity check for utc_lgm

The definition of utc_lgm [it is unclear why the field is not called simply utc] is:

Column	Variable name	type	Dim.	description
34	utc_lgm	double	1	UTC (0-24 hours)

 We offer a small correction:

Column	Variable name	type	Dim.	description
34	utc_lgm	double	1	UTC [0-24) hours

The actual values of this field lie in the following ranges:

utc_lgm
Satellite	Minimum	Maximum
53	5.567111e-03	2.396917e+01
54	5.567111e-03	2.398610e+01
55	1.916133e-02	2.399556e+01
56	6.108889e-03	2.397277e+01
57	1.248933e-02	2.399583e+01
58	2.527778e-02	2.398555e+01
59	5.563556e-03	2.399222e+01
60	1.949778e-03	2.395944e+01
61	9.177333e-03	2.398278e+01
62	9.993333e-03	2.398639e+01
63	7.219556e-03	2.399306e+01
64	1.554756e-02	2.399083e+01
65	2.001156e-02	2.399666e+01
66	3.249156e-02	2.399500e+01
67	1.055556e-02	2.398972e+01
68	1.111778e-02	2.395138e+01
69	1.638756e-02	2.399332e+01
70	4.445378e-02	2.399693e+01
71	1.221156e-02	2.399445e+01
72	2.527778e-02	2.396389e+01
73	4.111778e-02	2.397445e+01

The numbers are obviously reported to a precision that is a function of their magnitude. This is easily fixed by the following script, do-stage-35.py:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

field_nr = 33    # utc_lgm (wrt 0)

with open(filename) as records:
  for line in records:
    fields = line.split()
    value_str = fields[field_nr]
    value = float(value_str)
    
    output_format = '{:8.5f}'
    
    new_value_str = output_format.format(value)
  
    fields[field_nr] = new_value_str
    newline = " ".join(fields)
 
    print newline

Note, however, that the minima are orders of magnitude greater than the originally reported values for the field magnetic_local_time. There is  a slight inconsistency in the mechanisms for reporting these two types of time. It is hard to see how this will make any difference when actually using the dataset, but I note it in passing.

Stage 35: Reformat values of utc_lgm

Execute do-stage-34.py on all the relevant satellites:

for file in ns[567]*
do  
  ./do-stage-35.py $file > ../gps-stage-35/$file
done

The maxima and minima now have these values:

utc_lgm
Satellite	Minimum	Maximum
53	0.00557	23.96917
54	0.00557	23.98610
55	0.01916	23.99556
56	0.00611	23.97277
57	0.01249	23.99583
58	0.02528	23.98555
59	0.00556	23.99222
60	0.00195	23.95944
61	0.00918	23.98278
62	0.00999	23.98639
63	0.00722	23.99306
64	0.01555	23.99083
65	0.02001	23.99666
66	0.03249	23.99500
67	0.01056	23.98972
68	0.01112	23.95138
69	0.01639	23.99332
70	0.04445	23.99693
71	0.01221	23.99445
72	0.02528	23.96389
73	0.04112	23.97445

A checkpoint copy of the stage 35 files is available here, with MD5 checksum b85d41e615126e3912bd50aa55a4af54.

The data table for ns41 and ns48 still looks like this:

Column	Variable name	type	Dim.	Description
1	decimal_day	double	1	GPS time -- a number from 1 (1-Jan 00:00) to 366 (31-Dec 24:00) or 367 in leap years
2	Geographic_Latitude	double	1	Latitude of satellite (deg)
3	Geographic_Longitude	double	1	Longitude of satellite (deg)
4	Rad_Re	double	1	(radius of satellite)/Rearth
5-12	rate_electron_measured	double	8	Measured rate (Hz) in each of the 8 BDD electron channels (E1-E8)
13-20	rate_proton_measured	double	8	Measured rate (Hz) in each of the 8 BDD proton channels (P1-P8)
21	collection_interval	int	1	dosimeter collection period (seconds)
22	year	int	1	year (e.g. 2015)
23	decimal_year	double	1	decimal year = year + (decimal_day-1.0)/(days in year)
24	svn_number	int	1	SVN number of satellite
25	b_coord_radius	double	1	radius from earth's dipole axis (earth radii)
26	b_coord_height	double	1	height above the earth's dipole equatorial plane (earth radii)
27	magnetic_longitude	double	1	Magnetic longitude (degrees)
28	L_shell	double	1	L_shell (earth radii) -- I do not clearly understand the origin of the calculation, but it seems to be a dipole field/T-89
29	bfield_ratio	double	1	Bsatellite/Bequator
30	local_time	double	1	magnetic local time (0-24 hours)
31	b_sattelite	double	1	B field at satellite (gauss)
32	b_equator	double	1	B field at equator (on this field line I think) (gauss)
33-40	electron_background	double	8	estimated background in electron channels E1-E8 (Hz)
41-48	proton_background	double	8	estimated background in proton channels P1-P8 (Hz)
49	proton_activity	int	1	=1 if there is significant proton activity
50	electron_temperature	double	1	electron temperature from a one Maxwellian fit (MeV)
51	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
52-59	model_counts_electron_fit	double	8	E1-E8 rates from the 2-parameter Maxwellian fit to the electron data
60-67	dtc_counts_electron	double	8	Dead time corrected electron rates (from data, not fit)
68-97	integral_flux_instrument	double	30	(based on 2 parameter Maxwellian fit) integral of electron flux above integral_flux_energy[i] particles/(cm2sec)
98-127	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
128-142	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
143-157	electron_diff_flux	double	15	(based on 2 parameter Maxwellian fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))

And for the remaining satellites we now have (I have slightly edited the description of the magnetic_longitude and L_LGM_OP77IGRF fields):

Column	Variable name	type	Dim.	description
1	decimal_day	double	1	GPS time: a decimal number in the range [1, 367) in leap years or [1, 366) otherwise, representing the day of the year (1-Jan 00:00 to 31-Dec 24:00).
2	Geographic_Latitude	double	1	Latitude of satellite (°, N +ve)
3	Geographic_Longitude	double	1	Longitude of satellite (°, E +ve, measured from Greenwich meridian)
4	Rad_Re	double	1	Distance from centre of Earth, in units of Earth radii.
5-15	rate_electron_measured	double	11	Measured rate (Hz) in each of the 11 CXD electron channels
16-20	rate_proton_measured	double	5	Measured rate (Hz) in each of the 5 CXD proton channels (P1-P5)
21	LEP_thresh	double	1	LEP threshold in E1 channels, in keV
22	collection_interval	int	1	dosimeter collection period (seconds)
23	year	int	1	year (e.g. 2015)
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0) / (days in year)
25	SVN	int	1	Satellite Vehicle Number
26	b_coord_radius	double	1	Distance from dipole axis, in units of Earth radii.
27	b_coord_height	double	1	Distance from dipole equatorial plane, in units of Earth radii (N +ve).
28	magnetic_longitude	double	1	Magnetic longitude (°)
29	L_shell	double	1	L shell: McIlwain calculation according to model with T89 External Field, IGRF Internal Field.
30	L_LGM_TS04IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, TS04 External Field, IGRF Internal Field.
31	L_LGM_OP77IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, OP77 External Field, IGRF Internal Field
32	L_LGM_T89CDIP	double	1	LanlGeoMag L-shell McIlwain calculation, T89 External Field, Centered Dipole Internal Field
33	magnetic_local_time	double	1	magnetic local time [0-24) hours
34	utc_lgm	double	1	UTC [0-24) hours
35-45	electron_background	double	11	estimated background in electron channels E1-E11 (Hz)
46-50	proton_background	double	5	estimated background in proton channels P1-P5 (Hz)
51	proton_activity	int	1	=1 if there is significant proton activity
52	proton_temperature_fit	double	1	characteristic momentum -- R0 in the expression given above (MeV/c)
53	proton_density_fit	double	1	N0 parameter in fit to proton flux ((protons/(cm2 sec sr MeV))
54	electron_temperature_fit	double	1	electron temperature from a one Maxwellian fit (MeV)
55	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
56-60	model_counts_proton_fit_pf	double	5	P1-P5 rate from proton fit (using proton_temperature_fit, proton_density_fit)
61-71	model_counts_electron_fit	double	11	E1-E11 rates from the 9-parameter electron flux model
72-77	proton_integrated_flux_fit	double	6	integral of proton flux (based on fit) above 10, 15.85, 25.11, 30, 40, 79.43 MeV (proton kinetic energy)
78-107	integral_flux_instrument	double	30	(based on 9 parameter fit) integral of electron flux above integral_flux_energy[i] particles/(cm2 sec)
108-137	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
138-152	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
153-167	electron_diff_flux	double	15	(based on 9 parameter fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))
168-176	Efitpars	double	9	fit parameters for 9 parameter electron fit

More on Reporting Bust Rates in Contests

Here I began to address the problem of how to report and compare bust rates in contests in a defensible, objective manner.

At the end of that post, we decided that: rather than simply quoting some kind of "bust rate", a range should be quoted for each station, representing, say, the 99% confidence limit for the rate.

If we do that, and reorder the stations in order of decreasing upper limit (which seems like the most reasonable ordering: it means that we are 99.5% sure that the actual bust rate is less than this number), then we find: and then we follow this text with a table.

Let us apply this methodology to an example contest (CQ WW SSB for 2005, as that is the first CQ WW contest for which public logs are available), and look at the stations that appear to have the best rate of copying. What we tabulate is the 99.5% confidence limit for the probability of a bust, $p_{bust}$, as in this table:

2005 CQ WW SSB -- 99.5% confidence upper bound to $p_{bust}$
Position	Call	$p_{995}$
1	OH5BM	0.0005
2	9A2EU	0.0008
3	N5AU	0.0012
4	UA1OMX	0.0015
5	V31MQ	0.0018
6	UR6IJ	0.0019
7	KF2O	0.0020
8	DL3BRA	0.0021
9	EA1JO	0.0022
10	W3YY	0.0022

Let's expand the table so as to include the verified number of QSOs, $Q_v$, and the number of busts, $B$, for each station:

2005 CQ WW SSB -- 99.5% confidence upper bound to $p_{bust}$
Position	Call	$p_{995}$	$Q_v$	$B$
1	OH5BM	0.0005	1904	0
2	9A2EU	0.0008	1175	0
3	N5AU	0.0012	802	0
4	UA1OMX	0.0015	649	0
5	V31MQ	0.0018	563	0
6	UR6IJ	0.0019	533	0
7	KF2O	0.0020	488	0
8	DL3BRA	0.0021	470	0
9	EA1JO	0.0022	457	0
10	W3YY	0.0022	447	0

It should be comforting to observe that by ordering the stations in increasing order of $p_{995}$, we have also sorted it in decreasing order of $Q_v$.

However, note that none of the listed stations have any busts. What happens when we look at stations that have different numbers of busts? Let us look at an example further down the table:

2005 CQ WW SSB -- 99.5% confidence upper bound to $p_{bust}$
Position	Call	$p_{995}$	$Q_v$	$B$
4025	EC7ALM	0.6096	19	6
4026	GM8KSJ	0.6307	7	1

Let's now look at the actual probability distribution of $p_{bust}$ for these two stations (the vertical lines are at the two values of $p_{995}$):

It is obvious from this plot that the order of these two stations should be reversed in any meaningful table that purports to order stations in increasing order of estimates of $p_{bust}$. What is not so obvious is that the order of the two stations depends on the particular confidence limit chosen. For the 99.5% limit, the order is as stated, but (for example) for the 98% limit, the order is reversed. This arbitrariness is obviously as unacceptable as the basic notion that somehow GM8KSJ has a higher value of $p_{bust}$ than does EC7ALM.

So we need a different approach: although quoting a particular confidence limit is useful for defining $p_{bust}$ for a single station, it is inappropriate for ordering multiple stations, particularly in the case when different numbers of busts are involved (because the shapes of the probability curves differ markedly as $B$ varies).

A better way to order two stations is to determine the weighted mean value of $p_{bust}$ as distributed according to probability function of each station. This automatically results in an ordering such that if one selects a large number of values of $p_{bust}$ distributed according to the two relevant probability functions, the mean value for the first station will be less than the mean value for the second station.

[It is probably worth noting that, because of the lack of symmetry in the probability curves, this is not the same as the first station being necessarily more likely to have a lower value of  $p_{bust}$ than is the second. If this isn't obvious, just plot the difference between the probability curves for two stations with different values of $B$, such as EC7ALM and GM8KSJ.]

Note also that this procedure will leave the order of stations with equal numbers of busts unchanged, which is comforting.

Applying this procedure, our top-ten table now looks like this:

2005 CQ WW SSB -- weighted mean values of $p_{bust}$
Position	Call	weighted mean	$Q_v$	$B$
1	OH5BM	0.0005	1904	0
2	9A2EU	0.0008	1175	0
3	N5AU	0.0012	802	0
4	UA1OMX	0.0015	649	0
5	V31MQ	0.0017	563	0
6	UR6IJ	0.0018	533	0
7	ES5RY	0.0019	1068	1
8	KF2O	0.0020	488	0
9	DL3BRA	0.0021	470	0
10	EA1JO	0.0022	457	0

This looks reasonable (to me, anyway).

Further Sanity Checking of the LANL GPS Charged-Particle Dataset

Prior posts in this series may be found at:

• Original data
• Stages 1 to 4
• Revised stages 1 to 5  
• Stages 6 to 10 
• Stages 11 to 15
• Stages 16 to 20  
• Stages 21 to 25

In this post I continue looking at the values of fields for the CXD experiment carried on ns53 through ns73.

Sanity check for collection_interval

The values for collection_interval seem reasonable. (Recall that we processed these values earlier, at stage 4 and stage 10.)

Sanity check for year

The values for year seem reasonable.

Sanity check for decimal_year

The documentation for decimal_year states:

Column	Variable name	type	Dim.	description
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0)/(days in year)

But if we look at the values in the files, we see an immediate problem. The values are formatted as, for example:
  2.007879e+03
That is, the smallest increment is 0.000001e+03 years, which is 0.001 years, or somewhat less than nine hours. This is obviously useless.

The decimal_year field is redundant, as it can be derived, as stated in the documentation, from the year field and the decimal_day field. But if the decimal_year field is going to be present, it should at least have the same value as the result of the calculation in the documentation. There is obvious use for this field if it has an accurate value, so that brings us to define stage 26:

Stage 26: Recalculate decimal_year

The values of decimal_day are recorded to 10^-6 day. So decimal_year should be recorded to (roughly) 10^-9 year. We can accomplish the reformatting easily, by applying the following script, do-stage-26.py,  to each file:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

with open(filename) as records:
  for line in records:
    fields = line.split()
    decimal_day_str = fields[0]
    year_str = fields[22]
    decimal_day_value = float(decimal_day_str)
    year_value = int(year_str)
   
    days_in_year = 365
   
    if ( (year_value == 2004) or (year_value == 2008) or (year_value == 2012) or (year_value == 2016)):
      days_in_year = 366
     
    decimal_year_value = year_value + (decimal_day_value - 1.0) / days_in_year
     
    decimal_year_str = '{:14.9f}'.format(decimal_year_value)
   
    fields[23] = decimal_year_str
    newline = " ".join(fields)
  
    print newline

Applying this to each file:

 for file in ns[567]*
 do 
   ./do-stage-26.py $file | tr -s " " | sed 's/ $//' > ../gps-stage-26/$file
 done

The resulting maximum and minimum values of decimal_year look like this:

decimal_year
Satellite	Minimum	Maximum
53	2005.769864123	2016.999996964
54	2001.169868310	2016.999995825
55	2007.879454718	2017.000000000
56	2003.145206589	2016.999995066
57	2008.051914464	2017.000000000
58	2006.939732433	2016.999995825
59	2004.256832087	2016.999999620
60	2004.562844536	2016.999993169
61	2004.887981525	2016.999998104
62	2010.465755041	2016.999998325
63	2011.578088408	2016.999999653
64	2014.164385844	2016.923493000
65	2012.784155790	2016.999996205
66	2013.435620655	2016.999996964
67	2014.394527493	2016.999996396
68	2014.624659342	2016.999995003
69	2014.912332318	2016.999999716
70	2016.120224178	2016.999992601
71	2015.276714230	2016.999999874
72	2015.583565671	2016.999995825
72	2015.871238110	2016.999997628

Sanity check for SVN_number

The documentation for SVN_number states:

Column	Variable name	type	Dim.	description
25	SVN_number	int	1	SVN number of satellite

The values for SVN_number seem to be correct, except for the obvious silliness that the field name, when expanded, is Satellite_Vehicle_Number_number, and the description includes the nonsensical term SVN number. Consequently I simply redefine this field as:

Column	Variable name	type	Dim.	description
25	SVN	int	1	Satellite Vehicle Number

Sanity check for b_coord_radius

The maximum and minimum values for b_coord_radius are:

b_coord_radius
Satellite	Minimum	Maximum
53	6.828054e-01	4.200857e+00
54	8.312735e-01	4.211072e+00
55	8.557934e-01	4.198010e+00
56	7.165368e-01	4.202026e+00
57	7.242072e-01	4.171993e+00
58	7.577640e-01	4.188587e+00
59	NA	4.214008e+00
60	NA	4.205339e+00
61	NA	4.223202e+00
62	1.304388e+00	4.193189e+00
63	6.851737e-01	4.209296e+00
64	1.461084e+00	4.209824e+00
65	1.361361e+00	4.202571e+00
66	1.049070e+00	4.182467e+00
67	7.584448e-01	4.211562e+00
68	7.835127e-01	4.211309e+00
69	7.074734e-01	4.169802e+00
70	8.317451e-01	4.207506e+00
71	1.880229e+00	4.171600e+00
72	1.366415e+00	4.170103e+00
72	7.233613e-01	4.213620e+00

There are two obvious problems with these values: some satellites contain invalid (NA) values; and the precision of the field varies as a function of the value of the field.

We can remove all the records that contain the an invalid value easily enough, which gives us stage 27.

Stage 27: Remove records with invalid values of b_coord_radius

for file in ns[567]*
do
  awk '$26!="NA" { print $0 }' $file > ../gps-stage-27/$file
done

This gives us:

b_coord_radius
Satellite	Minimum	Maximum
53	6.828054e-01	4.200857e+00
54	8.312735e-01	4.211072e+00
55	8.557934e-01	4.198010e+00
56	7.165368e-01	4.202026e+00
57	7.242072e-01	4.171993e+00
58	7.577640e-01	4.188587e+00
59	6.926215e-01	4.214008e+00
60	8.069324e-01	4.205339e+00
61	7.764630e-01	4.223202e+00
62	1.304388e+00	4.193189e+00
63	6.851737e-01	4.209296e+00
64	1.461084e+00	4.209824e+00
65	1.361361e+00	4.202571e+00
66	1.049070e+00	4.182467e+00
67	7.584448e-01	4.211562e+00
68	7.835127e-01	4.211309e+00
69	7.074734e-01	4.169802e+00
70	8.317451e-01	4.207506e+00
71	1.880229e+00	4.171600e+00
72	1.366415e+00	4.170103e+00
72	7.233613e-01	4.213620e+00

Stage 28: Reformat values of b_coord_radius

As noted above, the accuracy of the numbers reported varies as a function of the value of the field. This is not a desirable situation. Since values greater than unity are reported with a precision of a millionth of a terrestrial radius, there seems no point in reporting values at greater precision (i.e., less than about 6m), especially since such a precision surely vastly exceeds the accuracy of the magnetic field model on which the numbers are based.

Consequently, we can reformat the values so that they have a consistent precision of one millionth of the terrestrial radius. We define a script called do-stage-28.py:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

with open(filename) as records:
  for line in records:
    fields = line.split()
    value_str = fields[25]      # b_coord_radius
    value = float(value_str)
     
    output_format = '{:8.6f}'
         
    new_value_str = output_format.format(value)
   
    fields[25] = new_value_str
    newline = " ".join(fields)
  
    print newline

And operate on the files:

for file in ns[567]*
do 
  ./do-stage-28.py $file | tr -s " " | sed 's/ $//' > ../gps-stage-28/$file
done

Now we have:

b_coord_radius
Satellite	Minimum	Maximum
53	0.682805	4.200857
54	0.831273	4.211072
55	0.855793	4.198010
56	0.716537	4.202026
57	0.724207	4.171993
58	0.757764	4.188587
59	0.692622	4.214008
60	0.806932	4.205339
61	0.776463	4.223202
62	1.304388	4.193189
63	0.685174	4.209296
64	1.461084	4.209824
65	1.361361	4.202571
66	1.049070	4.182467
67	0.758445	4.211562
68	0.783513	4.211309
69	0.707473	4.169802
70	0.831745	4.207506
71	1.880229	4.171600
72	1.366415	4.170103
72	0.723361	4.213620

Sanity check for b_coord_height

The maximum and minimum values for b_coord_height are:

b_coord_height
Satellite	Minimum	Maximum
53	-4.105137e+00	4.082754e+00
54	-4.011254e+00	4.119900e+00
55	-4.097133e+00	3.981823e+00
56	-4.091484e+00	4.110250e+00
57	-4.076309e+00	4.113386e+00
58	-4.051769e+00	4.112167e+00
59	-4.092374e+00	4.107748e+00
60	-3.965876e+00	4.114727e+00
61	-4.095934e+00	4.074463e+00
62	-3.918038e+00	3.950894e+00
63	-4.087862e+00	4.103668e+00
64	-3.801516e+00	3.899932e+00
65	-3.946078e+00	3.755339e+00
66	-3.898498e+00	4.027861e+00
67	-4.095983e+00	4.105099e+00
68	-4.039725e+00	4.085420e+00
69	-4.077435e+00	4.103233e+00
70	-4.077986e+00	4.069141e+00
71	-3.667803e+00	3.713735e+00
72	-3.727146e+00	3.938838e+00
72	-4.134635e+00	4.132903e+00

These require the same kind of massaging as b_coord_radius.

Stage 29: Reformat values of b_coord_height

We define do-stage-29.py:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

with open(filename) as records:
  for line in records:
    fields = line.split()
    value_str = fields[26]      # b_coord_height
    value = float(value_str)
    
    output_format = '{:8.6f}'
    
    new_value_str = output_format.format(value)
  
    fields[26] = new_value_str
    newline = " ".join(fields)
 
    print newline

and apply it to the files:

for file in ns[567]*
do 
  ./do-stage-29.py $file | tr -s " " | sed 's/ $//' > ../gps-stage-29/$file
done

with the result:

b_coord_height
Satellite	Minimum	Maximum
53	-4.105137	4.082754
54	-4.011254	4.119900
55	-4.097133	3.981823
56	-4.091484	4.110250
57	-4.076309	4.113386
58	-4.051769	4.112167
59	-4.092374	4.107748
60	-3.965876	4.114727
61	-4.095934	4.074463
62	-3.918038	3.950894
63	-4.087862	4.103668
64	-3.801516	3.899932
65	-3.946078	3.755339
66	-3.898498	4.027861
67	-4.095983	4.105099
68	-4.039725	4.085420
69	-4.077435	4.103233
70	-4.077986	4.069141
71	-3.667803	3.713735
72	-3.727146	3.938838
72	-4.134635	4.132903

Sanity check for magnetic_longitude

The maximum and minimum values for magnetic_longitude are:

magnetic_longitude
Satellite	Minimum	Maximum
53	4.185175e-03	3.599997e+02
54	1.637927e-04	3.599992e+02
55	3.244386e-05	3.600000e+02
56	3.396840e-04	3.599996e+02
57	5.629913e-04	3.599997e+02
58	1.443975e-03	3.599997e+02
59	2.364732e-05	3.600000e+02
60	3.059084e-04	3.599989e+02
61	6.888084e-04	3.599998e+02
62	9.731220e-04	3.599994e+02
63	1.367992e-03	3.599995e+02
64	3.209574e-04	3.599982e+02
65	1.819354e-03	3.599992e+02
66	8.708465e-04	3.599990e+02
67	1.153824e-05	3.599956e+02
68	6.599287e-04	3.599986e+02
69	5.731681e-03	3.599997e+02
70	3.041434e-02	3.599747e+02
71	5.421115e-04	3.599966e+02
72	1.369661e-03	3.599986e+02
72	6.566200e-04	3.599996e+02

We apply changes similar to stage 23, as follows.

Stage 30: reformat values of magnetic_longitude

We define do-stage-30.py:

#!/usr/bin/env python
# -*- coding: utf8 -*-

import re
import sys

filename = sys.argv[1]

field_nr = 27    # magnetic_longitude

with open(filename) as records:
  for line in records:
    fields = line.split()
    value_str = fields[field_nr]      # longitude
    value = float(value_str)
   
    if value < 0:
      value += 360
     
    if value >= 360:
      value -= 360
     
    new_value_str = '{:9.4f}'.format(value)
   
    fields[field_nr] = new_value_str
    newline = " ".join(fields)
  
    print newline

And apply it:

for file in ns[567]*
do 
  ./do-stage-30.py $file | tr -s " " | sed 's/ $//' > ../gps-stage-30/$file
done

With the result:

magnetic_longitude
Satellite	Minimum	Maximum
53	0.0042	359.9997
54	0.0002	359.9992
55	0.0000	359.9998
56	0.0003	359.9996
57	0.0006	359.9997
58	0.0014	359.9997
59	0.0000	359.9999
60	0.0003	359.9989
61	0.0007	359.9998
62	0.0010	359.9994
63	0.0014	359.9995
64	0.0003	359.9982
65	0.0018	359.9992
66	0.0009	359.9990
67	0.0000	359.9956
68	0.0007	359.9986
69	0.0057	359.9997
70	0.0304	359.9747
71	0.0005	359.9966
72	0.0014	359.9986
72	0.0007	359.9996

We  checkpoint the stage 30 dataset. The checkpoint file has the MD5 checksum: eba369c61f7974a636358758e1035ae9.

The data table for ns41 and ns48 still looks like this:

Column	Variable name	type	Dim.	Description
1	decimal_day	double	1	GPS time -- a number from 1 (1-Jan 00:00) to 366 (31-Dec 24:00) or 367 in leap years
2	Geographic_Latitude	double	1	Latitude of satellite (deg)
3	Geographic_Longitude	double	1	Longitude of satellite (deg)
4	Rad_Re	double	1	(radius of satellite)/Rearth
5-12	rate_electron_measured	double	8	Measured rate (Hz) in each of the 8 BDD electron channels (E1-E8)
13-20	rate_proton_measured	double	8	Measured rate (Hz) in each of the 8 BDD proton channels (P1-P8)
21	collection_interval	int	1	dosimeter collection period (seconds)
22	year	int	1	year (e.g. 2015)
23	decimal_year	double	1	decimal year = year + (decimal_day-1.0)/(days in year)
24	svn_number	int	1	SVN number of satellite
25	b_coord_radius	double	1	radius from earth's dipole axis (earth radii)
26	b_coord_height	double	1	height above the earth's dipole equatorial plane (earth radii)
27	magnetic_longitude	double	1	Magnetic longitude (degrees)
28	L_shell	double	1	L_shell (earth radii) -- I do not clearly understand the origin of the calculation, but it seems to be a dipole field/T-89
29	bfield_ratio	double	1	Bsatellite/Bequator
30	local_time	double	1	magnetic local time (0-24 hours)
31	b_sattelite	double	1	B field at satellite (gauss)
32	b_equator	double	1	B field at equator (on this field line I think) (gauss)
33-40	electron_background	double	8	estimated background in electron channels E1-E8 (Hz)
41-48	proton_background	double	8	estimated background in proton channels P1-P8 (Hz)
49	proton_activity	int	1	=1 if there is significant proton activity
50	electron_temperature	double	1	electron temperature from a one Maxwellian fit (MeV)
51	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
52-59	model_counts_electron_fit	double	8	E1-E8 rates from the 2-parameter Maxwellian fit to the electron data
60-67	dtc_counts_electron	double	8	Dead time corrected electron rates (from data, not fit)
68-97	integral_flux_instrument	double	30	(based on 2 parameter Maxwellian fit) integral of electron flux above integral_flux_energy[i] particles/(cm2sec)
98-127	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
128-142	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
143-157	electron_diff_flux	double	15	(based on 2 parameter Maxwellian fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))

And for the remaining satellites we now have (I have clarified the descriptions of the b_coord_radius and b_coord_height fields):

Column	Variable name	type	Dim.	description
1	decimal_day	double	1	GPS time: a decimal number in the range [1, 367) in leap years or [1, 366) otherwise, representing the day of the year (1-Jan 00:00 to 31-Dec 24:00).
2	Geographic_Latitude	double	1	Latitude of satellite (°, N +ve)
3	Geographic_Longitude	double	1	Longitude of satellite (°, E +ve, measured from Greenwich meridian)
4	Rad_Re	double	1	Distance from centre of Earth, in units of Earth radii.
5-15	rate_electron_measured	double	11	Measured rate (Hz) in each of the 11 CXD electron channels
16-20	rate_proton_measured	double	5	Measured rate (Hz) in each of the 5 CXD proton channels (P1-P5)
21	LEP_thresh	double	1	LEP threshold in E1 channels, in keV
22	collection_interval	int	1	dosimeter collection period (seconds)
23	year	int	1	year (e.g. 2015)
24	decimal_year	double	1	decimal year = year + (decimal_day-1.0) / (days in year)
25	SVN	int	1	Satellite Vehicle Number
26	b_coord_radius	double	1	Distance from dipole axis, in units of Earth radii.
27	b_coord_height	double	1	Distance from dipole equatorial plane, in units of Earth radii (N +ve).
28	magnetic_longitude	double	1	Magnetic longitude (degrees)
29	L_shell	double	1	L shell: McIlwain calculation according to model with T89 External Field, IGRF Internal Field.
30	L_LGM_TS04IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, TS04 External Field, IGRF Internal Field.
31	L_LGM_OP77IGRF	double	1	LanlGeoMag L-shell McIlwain calculation, OP77 External Field, IGRF Internal Field (not currently filled)
32	L_LGM_T89CDIP	double	1	LanlGeoMag L-shell McIlwain calculation, T89 External Field, Centered Dipole Internal Field
33	bfield_ratio	double	1	Bsatellite/Bequator
34	local_time	double	1	magnetic local time (0-24 hours)
35	utc_lgm	double	1	UTC (0-24 hours)
36	b_sattelite	double	1	B field at satellite (gauss)
37	b_equator	double	1	B field at equator (on this field line I think) (gauss)
38-48	electron_background	double	11	estimated background in electron channels E1-E11 (Hz)
49-53	proton_background	double	5	estimated background in proton channels P1-P5 (Hz)
54	proton_activity	int	1	=1 if there is significant proton activity
55	proton_temperature_fit	double	1	characteristic momentum -- R0 in the expression given above (MeV/c)
56	proton_density_fit	double	1	N0 parameter in fit to proton flux ((protons/(cm2 sec sr MeV))
57	electron_temperature_fit	double	1	electron temperature from a one Maxwellian fit (MeV)
58	electron_density_fit	double	1	electron number density from a one Maxwellian fit (cm-3)
59-63	model_counts_proton_fit_pf	double	5	P1-P5 rate from proton fit (using proton_temperature_fit, proton_density_fit)
64-74	model_counts_electron_fit	double	11	E1-E11 rates from the 9-parameter electron flux model
75-80	proton_integrated_flux_fit	double	6	integral of proton flux (based on fit) above 10, 15.85, 25.11, 30, 40, 79.43 MeV (proton kinetic energy)
81-110	integral_flux_instrument	double	30	(based on 9 parameter fit) integral of electron flux above integral_flux_energy[i] particles/(cm2 sec)
111-140	integral_flux_energy	double	30	energies for the integral of integral_flux_instrument (MeV)
141-155	electron_diff_flux_energy	double	15	energies for the fluxes in electron_diff_flux_energy (MeV)
156-170	electron_diff_flux	double	15	(based on 9 parameter fit) electron flux at energies electron_diff_flux[i] (particle/(cm2 sr MeV sec))
171-179	Efitpars	double	9	fit parameters for 9 parameter electron fit

Finally, here are the number of records for each satellite at the end of the processing for stage 30:

Satellite	Stage 30 Records
ns41	1,990,340
ns48	1,105,300
ns53	1,331,017
ns54	1,938,603
ns55	1,054,569
ns56	1,680,369
ns57	1,082,028
ns58	1,174,629
ns59	1,513,323
ns60	1,492,497
ns61	1,467,283
ns62	774,976
ns63	651,078
ns64	343,643
ns65	480,017
ns66	446,139
ns67	327,174
ns68	304,964
ns69	262,021
ns70	110,260
ns71	220,694
ns72	181,519
ns73	145,292

Switching From Nouveau to the Proprietary NVIDIA driver in Debian Jessie

A couple of weeks ago, I started to experience random crashes on my 64-bit jessie desktop machine. Generally, these took the form of either a frozen desktop or a sudden blank screen, along with a complete lack of responsiveness to either mouse or keyboard.

Usually there was no obvious associated entry in any system log, but at last a series of messages appeared in the syslog file, starting with this one:

Jul 17 13:55:05 homebrew kernel: [24064.296254] nouveau E[ PFIFO][0000:01:00.0] write fault at 0x000029d000 [PTE] from GR/GPC0/GPCCS on channel 0x003fbad000 [Xorg[2071]] 

This was followed by several more messages that appeared to be related, the last of which was:

Jul 17 13:58:23 homebrew kernel: [24262.075187] nouveau E[ DRM] GPU lockup - switching to software fbcon

(On this occasion, although the desktop was non-responsive after the first message, I could still ssh into the machine, and shut it down cleanly from the ssh session, which is why there is a period of several minutes between these two messages.)

This suggested that the cause lay in the nouveau video driver, so I decided to switch to the proprietary NVIDIA driver. This turned out not to be as easy as one might expect, since there didn't seem to be a single place that defines the complete procedure  in detail. Hence this post.

Here are the steps that I followed:

1. Install the nvidia-driver package.

2. Install the nvidia-xconfig package.

3. Run nvidia-xconfig.

This complained about the lack of an xorg.conf file, but generated a default one with an nvidia entry for the driver.  There were several other errors, but rebooting at this point resulted in a system that booted and ran X.

So far so good, but during the boot sequence I noticed that the text on the system console was enormous. Similarly, if I switched to the console once the system had booted, the text appeared to be about 80x24, which is quite obnoxious on a 27-inch monitor.

Following the instructions at:

https://wiki.archlinux.org/index.php/GRUB/Tips_and_tricks#Setting_the_framebuffer_resolution

I added two lines to the file /etc/default/grub:

GRUB_GFXMODE=1280x1024x16,1024x768,auto 
GRUB_GFXPAYLOAD_LINUX=keep

DO NOT DO THIS.

After executing
  grub-mkconfig -o /boot/grub/grub.cfg
and rebooting, although the text on the console looked much better, I no longer had any X-based desktop. Switching to :0 merely gave me a blank screen. So I restored the grub.cfg file to the original version.

The above-named URL provides a deprecated mechanism for changing the console font, so that's what I ended up using. In particular, I changed one line of the /etc/default/grub file to read:

GRUB_CMDLINE_LINUX_DEFAULT="quiet vga=794"

and executed: 
  grub-mkconfig -o /boot/grub/grub.cfg

According to this documentation, this gives a 1280x1024 16-bit console, which is a somewhat lower resolution than I had with the nouveau driver, but is vastly better than the resolution without this line in the grub configuration file.

Now everything is working to my satisfaction. The only quirk I see is that at boot time, there is a LOT of disk activity for about 30 seconds after the desktop starts. I'm not sure what the reason for this might be, but at the end of it I have a fully-functioning system with a KDE desktop on :0 and i3 on :1, and can switch to a reasonable-looking console at will.

The best news is that, at least so far, I have experienced no system crashes since switching to the proprietary driver.

Most-Logged Stations in CQ WW CW 2016

The public CQ WW CW logs allow us easily to tabulate the stations that appear in the largest number of entrants' logs. For 2016, the ten stations with the largest number of appearances were:

Callsign	Appearances	% logs
HK1NA	10,277	59
9A1A	9,921	63
TK0C	9,682	61
CN2R	9,675	63
PJ2T	9,569	53
CR3W	9,354	60
LZ9W	9,091	59
CN2AA	8,912	59
P33W	8,661	56
EF8R	8,409	57

The first column in the table is the callsign. The second column is the total number of times that the call appears in other stations' logs. That is, if a station worked HK1NA on six bands, that will increment the value in the second column of the HK1NA row by six. The third column is the percentage of logs that contain the callsign at least once.

Tables for prior years are available here.

We can also list the cumulative data for the ten-year span from 2007 to 2016:

Callsign	Appearances	% logs
LZ9W	91,382	67
PJ2T	84,001	58
9A1A	81,583	61
DF0HQ	80,298	63
PJ4A	72,262	55
D4C	71,571	49
W3LPL	69,447	52
LX7I	69,254	55
K3LR	68,908	53
CR3L	64,426	48

Additional Information in Augmented Logs for CQ WW, 2005 to 2016

Now available are new augmented versions of the public logs for CQ WW CW and SSB for the period 2005 to 2016.

The cleaned logs are the result of processing the QSO: lines from the entrants' submitted Cabrillo files to ensure that all fields contain valid values and all the data match the format required in the rules. Any line containing illegal data in a field (for example, a zone number greater than 40, or a date/time stamp that is outside the contest period) has simply been removed. Also, only the QSO: lines are retained, so that each line in the file can be processed easily. The MD5 checksum for the file of cleaned logs is: 1b47059d1f2431b55d89a5eb954a05cc.

The augmented logs contain the same information as the cleaned logs, with the addition of some useful information on each line. The MD5 checksum for the compressed (~800 MB) file of augmented logs is: 7a728987fb8637ab8c156df3fa27d582. The information added to each line now includes two new fields: the callsign copied by the second party in the case that the second party bust the cull of the first party; amd the correct callsign of the second party in the case that the first party bust the second party's call.

In all, the addition fields in the augmented file comprise:

1. The letter "A" or "U" indicating "assisted" or "unassisted"
2. A four-digit number representing the time if the contact in minutes measured from the start of the contest. (I realise that this can be calculated from the other information on the line, but it saves a lot of time to have the number readily available in the file without having to calculate it each time.)
3. Band
4. A set of eleven flags, each -- apart from column k -- encoded as T/F: 
   • a. QSO is confirmed by a log from the second party 
   • b. QSO is a reverse bust (i.e., the second party appears to have bust the call of the first party) 
   • c. QSO is an ordinary bust (i.e., the first party appears to have bust the call of the second party) 
   • d. the call of the second party is unique 
   • e. QSO appears to be a NIL 
   • f. QSO is with a station that did not send in a log, but who did make 20 or more QSOs in the contest 
   • g. QSO appears to be a country mult 
   • h. QSO appears to be a zone mult 
   • i. QSO is a zone bust (i.e., the received zone appears to be a bust)
   • j. QSO is a reverse zone bust (i.e. the second party appears to have bust the zone of the first party)
   • k. This entry has three possible values rather than just T/F:
     • T: QSO appears to be made during a run by the first party
     • F: QSO appears not to be made during a run by the first party
     • U: the run status is unknown because insufficient frequency information is available in the first party's log 
5. If the QSO is a reverse bust, the call logged by the second party; otherwise, the placeholder "-"
6. If the QSO is an ordinary bust, the correct call that should have been logged by the first party; otherwise, the placeholder "-"
7. If the QSO is a reverse zone bust, the zone logged by the second party; otherwise, the placeholder "-"
8.  If the QSO is an ordinary zone bust, the correct zone that should have been logged by the first party; otherwise, the placeholder "-"

Notes:

• The encoding of some of the flags requires subjective decisions to be made as to whether the flag should be true or false; consequently, and because CQ has yet to understand the importance of making their scoring code public, the value of a flag for a specific QSO line in some circumstances might not match the value that CQ would assign. (Also, CQ has more data available in the form of check logs, which are not made public.)
• I made no attempt to deduce the run status of a QSO in the second party's log (if such exists), regardless of the status in the first party's log. This allows one cleanly to perform correct statistical analyses anent the number of QSOs made by running stations merely by excluding QSOs marked with a U in column k.
• No attempt is made to detect the case in which both participants of a QSO bust the other station's call. This is a problematic situation because of the relatively high probability of a false positive unless both stations log the frequency as opposed to the band. (Also, on bands on which split-frequency QSOs are common, the absence of both transmit and receive frequency is a problem.) Because of the likelihood of false positives, it seems better, given the presumed rarity of double-bust QSOs, that no attempt be made to mark them.
• The entries for the zones in the case of zone or reverse zone busts are normalised to two-digit values.