Replicate Apollo sample magnetizations (Tarduno, Cottrell & Zhou 2024) on Figshare 24639285¶
Paper: J. A. Tarduno, R. D. Cottrell, T. Zhou, "Dataset of replicate Apollo sample magnetizations bearing on impacts and absence of a long-lived lunar dynamo", Scientific Data (2024), doi:10.1038/s41597-024-03626-0.
Dataset: Figshare item 24639285, version 1, doi:10.6084/m9.figshare.24639285.v1. Version 1 is the
only published version (the API lists no other), and it is the one the paper cites (ref. 26) and the
claims file names, so it is pinned. The cell below prints the MD5 (computed_md5) of every file read.
The files are streamed from Figshare's download service into memory; nothing is written to disk.
Claims: claims/figshare-24639285/10.1038-s41597-024-03626-0.yaml, written from the version of
record (Europe PMC PMC11271266, full text with Methods, Technical Validation, figure legends and
Tables 1 and 2; the paper has no supplement and no PDF was supplied).
| Claim | Tested here |
|---|---|
| glass-specimen-set (link) | tested |
| glass-success-rate | tested |
| table2-fit-statistics | tested |
| table2-no-fit | tested |
| fig1-linear-components | tested |
| fig2-changing-field | tested |
| glass-directions-not-random | tested |
| glass-af-no-change | tested |
| table1-max-intensity | tested |
| baseline-range | tested |
| example-14053g1-threshold | tested |
| fig3-nrm-consistent | tested |
| fig3-590-inconsistent | tested |
| fig3-alteration-check-inconsistent | tested |
| crystal-14053g1-400c-stable | tested |
| no-alteration-recording-efficiency | tested |
| interpretation-impact-no-dynamo | blocked: an interpretation the claims file says is not tested |
Method details the claims file leaves open are taken from the paper where it gives them, and cited
in apollo.py. Where the paper does not settle a detail, every reading is run and shown.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import apollo as A
pd.set_option("display.width", 200)
pd.set_option("display.max_rows", 200)
files, glass, crystals = A.load_all()
print(f"Figshare {A.FIGSHARE_ID} version {A.FIGSHARE_VERSION}, doi:{A.DATASET_DOI}")
print(pd.DataFrame([{k: f[k] for k in ("id", "name", "size", "computed_md5")} for f in files]).to_string(index=False))
Figshare 24639285 version 1, doi:10.6084/m9.figshare.24639285.v1
id name size computed_md5
43480524 12021g2.asc.xlsx 10274 efb317fa1a51337e86171197a72b7f89
43480527 12021g3.asc.xlsx 9202 8e1e9aa64846df8b6857030cc2f9fe0e
43480530 14053g2.asc.xlsx 9210 1d6cc886ddd2d30860382cba7747440f
43480533 71055g3.asc.xlsx 9157 5e019090b526baaac272150f3597b168
43480536 71055g1.asc.xlsx 9147 1ab0e3abd043138e8d47811d52e996c7
43480539 14053g3.asc.xlsx 9232 f635b8a25c48625291f2a66592308f97
43480542 14053g1.asc.xlsx 9505 766af3788f1c0674883bc3fd75125754
43480545 12053g1.asc.xlsx 9213 e08eb790633722ce89bf6a7da67f49d5
45441655 12053g2.asc.xlsx 9398 c60a6aef638797185f7ec67b3b91dd67
43480548 12040g3.asc.xlsx 9251 8ecdf1b09112acab9335a9269b5ad7a0
45441658 12040g2.asc.xlsx 9279 ff9bd1717dc42b52fd7201c3badd0f47
43480551 12040g1.asc.xlsx 9131 7d31eb138eab5e8ea9984234e47e0eb9
43480554 64455lunar_cryodata.xlsx 114324 8a12aa20c7139cc2d2136af97d69114f
44008632 README_Pmag GUI.pdf 756016 581910f5773c5f3de6a7de44134ede39
45556176 SupplFig1 - other Thelliers March 2024.pdf 27710777 ee2ba78e902ca25e978606346abd5bc8
45556230 SuppFig2 - All Stereonet.pdf 562508 6750717c0e5076e1c705f0ed91c055fa
Data check¶
What was loaded: every sheet of every workbook, its columns, the number of rows per specimen, step kind and replicate, and the ranges and units of moment, declination and inclination. The PDFs (plots and a PmagPy readme) are listed but not read. The data are static measurements with no time axis, so there is no clock to align; the check of the description is against the step labels, the units, the specimen names and the laboratory field.
for f in files:
if f["name"].endswith(".xlsx"):
print(f"{f['name']}: sheets {f['sheets']}")
print()
def summary(name, sheet):
df = sheet["data"]
kinds = df.step.map(lambda s: A.step_kind(s)[0])
reps = df.groupby("step").size()
return {"specimen": name, "file": sheet["file"], "sample label": ",".join(sorted(df["sample"].unique())),
"columns": "|".join(map(str, sheet["header"])), "rows": len(df),
"steps": df.step.nunique(), "AF steps": int((kinds == "af").sum() > 0) * df[kinds == "af"].step.nunique(),
"thermal T": ",".join(str(int(t)) for t in sorted({A.step_kind(s)[1] for s in df.step if A.step_kind(s)[0] in ("off", "on")})),
"replicates min-max": f"{reps.min()}-{reps.max()}", "lab field uT": sheet["lab_field_uT"],
"incomplete rows": sheet["incomplete_rows"]}
gsum = pd.DataFrame([summary(k, v) for k, v in glass.items()])
csum = pd.DataFrame([summary(k, v) for k, v in crystals.items()])
print("Glass, 64455lunar_cryodata.xlsx:")
print(gsum.drop(columns="file").to_string(index=False))
print("\nCrystals, one workbook each:")
print(csum.drop(columns=["file", "lab field uT", "AF steps"]).to_string(index=False))
allrows = pd.concat([v["data"].assign(specimen=k) for d in (glass, crystals) for k, v in d.items()])
print("\nRanges over all measurements: moment %.3g to %.3g emu; declination %.1f to %.1f deg; inclination %.1f to %.1f deg"
% (allrows.moment.min(), allrows.moment.max(), allrows.dec.min(), allrows.dec.max(), allrows.inc.min(), allrows.inc.max()))
bad = allrows[(allrows.inc.abs() > 90) | (allrows.dec < 0) | (allrows.dec > 360) | (allrows.moment <= 0)]
print("Measurements outside the ranges of a declination (0-360), an inclination (-90 to 90) or a positive moment:")
print(bad[["specimen", "row", "step", "moment", "dec", "inc"]].to_string(index=False))
print("Footnotes of the glass sheets:", sorted({n for v in glass.values() for n in v["notes"]}))
print("Footnotes of the crystal sheets:", sorted({n for v in crystals.values() for n in v["notes"]}))
12021g2.asc.xlsx: sheets ['12021g2']
12021g3.asc.xlsx: sheets ['12021g3']
14053g2.asc.xlsx: sheets ['Sheet1']
71055g3.asc.xlsx: sheets ['Sheet1']
71055g1.asc.xlsx: sheets ['Sheet1']
14053g3.asc.xlsx: sheets ['Sheet1']
14053g1.asc.xlsx: sheets ['Sheet1']
12053g1.asc.xlsx: sheets ['Sheet1']
12053g2.asc.xlsx: sheets ['Sheet1']
12040g3.asc.xlsx: sheets ['Sheet1']
12040g2.asc.xlsx: sheets ['Sheet1']
12040g1.asc.xlsx: sheets ['Sheet1']
64455lunar_cryodata.xlsx: sheets ['ss19', 'ss22', 'ss27', 'ss34', 'ss38', 'ss41', 'ss43', 'ss45', 'ss46', 'ss48', 'ss6', 'ss9', 'ss10', 'ss11', 'ss18', 'ss24', 'ss29', 'ss30', 'ss36', 'g2', 'ss8', 'ss39']
Glass, 64455lunar_cryodata.xlsx:
specimen sample label columns rows steps AF steps thermal T replicates min-max lab field uT incomplete rows
ss19 ss19 Sample|DEMAG|Moment(emu)|Dec|Inc 95 19 3 100,175,250,360,450,500,550 5-5 18.0 []
ss22 ss22 Sample|DEMAG|Moment(emu)|Dec|Inc 110 22 3 100,175,250,360,450,500,550,580 5-5 18.0 []
ss27 ss27 Sample|DEMAG|Moment(emu)|Dec|Inc 105 17 4 100,305,400,500,550 5-10 18.0 []
ss34 ss34 Sample|DEMAG|Moment(emu)|Dec|Inc 87 18 4 100,175,250,360,450,500 4-5 18.0 []
ss38 ss38 Sample|DEMAG|Moment(emu)|Dec|Inc 65 14 3 100,175,250,360 2-5 18.0 []
ss41 ss41 Sample|DEMAG|Moment(emu)|Dec|Inc 115 23 2 100,175,250,360,400,450,500,550,580 5-5 18.0 []
ss43 ss43 Sample|DEMAG|Moment(emu)|Dec|Inc 70 14 2 100,175,250,360,400 5-5 18.0 []
ss45 ss45 Sample|DEMAG|Moment(emu)|Dec|Inc 60 12 2 100,175,250,360 5-5 18.0 []
ss46 ss46 Sample|DEMAG|Moment(emu)|Dec|Inc 50 10 2 100,175,250,360 5-5 18.0 []
ss48 ss48 Sample|DEMAG|CSD|Moment(emu)|Dec|Inc 35 7 2 100,175 5-5 18.0 []
ss6 ss6 Sample|DEMAG|Moment(emu)|Dec|Inc 99 9 0 100,175,250,360 9-19 15.0 []
ss9 ss9 Sample|DEMAG|Moment(emu)|Dec|Inc 116 14 0 100,175,250,360,450,550 5-13 15.0 []
ss10 ss10 Sample|DEMAG|Moment(emu)|Dec|Inc 122 14 0 100,175,250,360,450,550 4-16 15.0 []
ss11 ss11 Sample|DEMAG|Moment(emu)|Dec|Inc 125 14 0 100,175,250,360,450,550 5-13 15.0 []
ss18 ss18 Sample|DEMAG|Moment(emu)|Dec|Inc 81 12 0 100,175,250,360,450 3-15 15.0 [41]
ss24 ss24 Sample|DEMAG|Moment(emu)|Dec|Inc 146 17 0 100,175,250,360,450,500,550 5-15 15.0 []
ss29 ss29 Sample|DEMAG|Moment(emu)|Dec|Inc 128 14 0 100,175,250,450,500,550 3-13 15.0 []
ss30 ss30 Sample|DEMAG|Moment(emu)|Dec|Inc 103 12 0 100,175,250,450,500 3-11 15.0 []
ss36 ss36 Sample|DEMAG|Moment(emu)|Dec|Inc 72 9 0 100,175,250,450 5-13 15.0 []
g2 64455g2 Sample|DEMAG|Moment(emu)|Dec|Inc 43 9 0 100,210,305,400 4-5 14.4 []
ss8 ss8 Sample|DEMAG|Moment(emu)|Dec|Inc 57 12 0 100,175,250,360,450 4-5 15.0 []
ss39 ss39 Sample|DEMAG|Moment(emu)|Dec|Inc 23 5 0 100,210 4-5 18.0 []
Crystals, one workbook each:
specimen sample label columns rows steps thermal T replicates min-max incomplete rows
12021g2 12021g2 Sample|DEMAG|Moment|Dec|Inc 10 5 590 1-3 []
12021g3 12021g3 Sample|DEMAG|Moment (emu)|Dec|Inc 10 5 590 1-3 []
14053g2 14053g2 Sample|DEMAG|Moment|Dec|Inc 10 5 590 1-3 []
71055g3 71055g3 Sample|DEMAG|Moment|Dec|Inc 9 5 590 1-3 []
71055g1 71055g1 Sample|DEMAG|Moment|Dec|Inc 9 5 590 1-3 []
14053g3 14053g3 Sample|DEMAG|Moment|Dec|Inc 11 5 590 1-4 []
14053g1 14053g1 Sample|DEMAG|Moment|Dec|Inc 16 8 400,590 1-3 []
12053g1 12053g1 Sample|DEMAG|Moment|Dec|Inc 10 5 590 1-3 []
12053g2 12053g2 Sample|DEMAG|Moment|Dec|Inc 11 5 590 1-3 []
12040g3 12040g3 Sample|DEMAG|Moment|Dec|Inc 10 5 590 1-3 []
12040g2 12040g2 Sample|DEMAG|Moment|Dec|Inc 11 5 590 1-3 []
12040g1 12040g1 Sample|DEMAG|Moment|Dec|Inc 10 5 590 1-3 []
Ranges over all measurements: moment 1.02e-10 to 6.04e-07 emu; declination 0.0 to 361.7 deg; inclination -136.0 to 90.0 deg
Measurements outside the ranges of a declination (0-360), an inclination (-90 to 90) or a positive moment:
specimen row step moment dec inc
ss19 15 10mT 7.173000e-08 241.400 -136.000
ss8 5 0 8.508274e-08 361.692 6.834
Footnotes of the glass sheets: ['#.0 Field off step', '#.1 Field on step', '#.2 pTRM check', 'For samples with an alternating demagnetization field pre-treatment, steps are labeled with mT.', 'Laboratory Field: 14.4 microTesla', 'Laboratory Field: 15 microTesla', 'Laboratory Field: 18 microTesla', 'Magnetometer data for glass subsample. Sample, Demagnetization step, Moment (in electromagnetic units), Declination (degrees), Inclination (degrees).']
Footnotes of the crystal sheets: ['#.0 Field off step', '#.1 Field on step at 20 microTesla laboratory field', '#.3 Field off step following in field step', '#.4 Field on step at 40 microTesla laboratory field', 'Magnetometer data for silicate subsample. Sample, Demagnetization step, Moment (in electromagnetic units), Declination (degrees), Inclination (degrees).']
Repeated measurements of a step should differ by measurement noise. The cell below tests, for every
step of every sheet, whether all its rows are one row multiplied by factors k/100 between 0.90 and
1.10 in moment, declination and inclination alike (apollo.synthetic_steps). For four-digit
measurements this cannot happen by chance. It also lists rows repeated exactly within a sheet and
across sheets.
syn = {k: A.synthetic_steps(v["data"]) for d in (glass, crystals) for k, v in d.items()}
for k, steps in syn.items():
n = (glass.get(k) or crystals.get(k))["data"].step.nunique()
if steps:
print(f"{k}: {len(steps)} of {n} steps are scaled copies of one row: {steps}")
print("Sheets with no such step:", [k for k, s in syn.items() if not s])
def factors(df, step):
g = df[df.step == step]
out = {}
for c in ("moment", "dec", "inc"):
v = g[c].abs().to_numpy()
for v0 in v:
for f0 in np.arange(90, 111) / 100:
r = v / (v0 / f0)
if np.all(np.abs(r * 100 - np.round(r * 100)) < 1e-6 * r * 100):
out[c] = np.round(r, 2).tolist()
break
if c in out:
break
return out
for k in ("g2", "ss8", "ss39", "ss48"):
print(f"{k} NRM rows, factors relative to the base row:", factors(glass[k]["data"], "0"))
print(glass["ss8"]["data"].query("step == '0'")[["row", "moment", "dec", "inc"]].to_string(index=False))
dups = []
for k, v in {**glass, **crystals}.items():
d = v["data"]
n = d.duplicated(["step", "moment", "dec", "inc"]).sum()
if n:
dups.append((k, int(n), k in ("g2", "ss8", "ss39", "ss48")))
print("Rows repeated exactly within a sheet (specimen, count, synthetic sheet):", dups)
key = allrows.set_index(["moment", "dec", "inc"])
cross = allrows.groupby(["moment", "dec", "inc"]).specimen.nunique()
print("Rows found in more than one sheet:")
print(allrows.set_index(["moment", "dec", "inc"]).loc[cross[cross > 1].index][["specimen", "row", "step"]])
ss48: 7 of 7 steps are scaled copies of one row: ['0', '5mT', '10mT', '100', '100.1', '175', '175.1']
g2: 9 of 9 steps are scaled copies of one row: ['0', '100.1', '100', '210.1', '210', '305.1', '305', '400.1', '400']
ss8: 12 of 12 steps are scaled copies of one row: ['0', '100.1', '100', '175.1', '175', '250.1', '250', '360.1', '360', '450.1', '250.2', '450']
ss39: 5 of 5 steps are scaled copies of one row: ['0', '100.1', '100', '210.1', '210']
Sheets with no such step: ['ss19', 'ss22', 'ss27', 'ss34', 'ss38', 'ss41', 'ss43', 'ss45', 'ss46', 'ss6', 'ss9', 'ss10', 'ss11', 'ss18', 'ss24', 'ss29', 'ss30', 'ss36', '12021g2', '12021g3', '14053g2', '71055g3', '71055g1', '14053g3', '14053g1', '12053g1', '12053g2', '12040g3', '12040g2', '12040g1']
g2 NRM rows, factors relative to the base row: {'moment': [1.0, 0.97, 1.03, 0.98], 'dec': [1.0, 0.97, 1.03, 0.98], 'inc': [1.0, 0.97, 1.03, 0.98]}
ss8 NRM rows, factors relative to the base row: {'moment': [1.01, 0.99, 0.99, 1.02], 'dec': [1.01, 0.99, 0.99, 1.02], 'inc': [1.01, 0.99, 0.99, 1.02]}
ss39 NRM rows, factors relative to the base row: {'moment': [1.02, 1.01, 0.96, 1.0], 'dec': [1.02, 1.01, 0.96, 1.0], 'inc': [1.02, 1.01, 0.96, 1.0]}
ss48 NRM rows, factors relative to the base row: {'moment': [1.03, 0.99, 0.95, 1.0, 1.03], 'dec': [0.97, 0.95, 0.96, 0.98, 1.03], 'inc': [0.99, 0.98, 1.01, 0.97, 1.02]}
row moment dec inc
2 8.424859e-08 358.146 6.767
3 8.258031e-08 351.054 6.633
4 8.258031e-08 351.054 6.633
5 8.508274e-08 361.692 6.834
Rows repeated exactly within a sheet (specimen, count, synthetic sheet): [('ss27', 1, False), ('ss38', 1, False), ('ss43', 2, False), ('ss45', 1, False), ('ss6', 1, False), ('ss9', 2, False), ('ss10', 3, False), ('ss11', 3, False), ('ss18', 1, False), ('ss24', 1, False), ('ss29', 4, False), ('ss36', 1, False), ('g2', 9, True), ('ss8', 14, True), ('ss39', 5, True)]
Rows found in more than one sheet:
specimen row step
moment dec inc
1.072000e-09 315.2 75.9 ss43 66 400
75.9 ss10 50 250
Checks of the description against the data:
- Columns: every sheet has sample, step, moment (emu), declination and inclination, as the Data
Records say, except ss48, which has an extra column
CSD(empty or 0) before the moment. The moment column is found by name. - Specimen names: the glass sheet
ss11is nameds11in Table 2; the sheetg2labels its rows64455g2. Twelve crystal workbooks match the twelve names of Table 1. - Step labels follow the description (
#.0,#.1,#.2for the glass;#.3,#.4for the crystals); NRM appears as both0and0.0(ss27, ss38). - Laboratory field: the glass sheets give 18 uT (10 sheets), 15 uT (11) or 14.4 uT (g2); the paper
gives 18 uT only for Figure 1, whose specimens (ss19, ss43) are 18 uT sheets. The crystal sheets
give 20 and 40 uT for
#.1and#.4, as the Methods say. - Replicates: each glass step was measured 2 to 19 times, each crystal NRM 1 to 3 times, each 590 zero-field step 3 or 4 times, and each crystal field-on step once. The paper does not say how repeated measurements enter Table 2; both the mean vector and the first measurement are run below.
- Four glass sheets (g2, ss8, ss39, ss48) do not hold repeated measurements: every step's rows are one row times factors near 1. The analysis keeps them, as the deposit presents them, and the results that use them say so; their mean vector is close to the base row.
- Values out of range: one inclination of -136 deg (ss19, 10 mT, almost certainly -13.6), one declination of 361.7 deg (ss8 NRM, a scaled copy), and an ss18 row at 250 C with no values. The first two are used as they are (converted to vectors, a declination of 361.7 is 1.7 deg); the ss19 value is in an AF step no analysis below uses; the empty row is dropped.
glass-specimen-set (link)¶
The deposit holds the measurements of the 22 Apollo 64455 glass subsamples that failed the paleointensity selection criteria, named as in Table 2, and replicate measurements for the twelve single crystals of Table 1.
Criterion: specimens matched by the names of Tables 1 and 2 (s11 may be ss11); each has steps, moments, declinations and inclinations; which glass specimens have AF steps is compared with the reading of the Table 2 marks (12 without, 10 with).
TABLE2 = ["ss19", "ss43", "ss22", "ss41", "ss6", "ss9", "ss10", "ss30", "ss34", "g2", "ss8", "s11", "ss18",
"ss24", "ss27", "ss29", "ss36", "ss38", "ss39", "ss45", "ss46", "ss48"]
# Table 2 marks: the second mark (triangle) read as AF pre-treatment (claims file, opening comment)
TABLE2_AF_MARK = ["ss19", "ss43", "ss22", "ss41", "ss34", "ss27", "ss38", "ss45", "ss46", "ss48"]
TABLE1 = ["12021g2", "12021g3", "12040g1", "12040g2", "12040g3", "12053g1", "12053g2", "14053g1", "14053g2",
"14053g3", "71055g1", "71055g3"]
SHEET = {n: ("ss11" if n == "s11" else n) for n in TABLE2}
def complete(sheet):
d = sheet["data"]
return len(d) > 0 and d[["moment", "dec", "inc"]].notna().all().all() and d.step.nunique() > 1
found_glass = [n for n in TABLE2 if SHEET[n] in glass and complete(glass[SHEET[n]])]
found_crystals = [n for n in TABLE1 if n in crystals and complete(crystals[n])]
extra = sorted(set(glass) - set(SHEET.values())) + sorted(set(crystals) - set(TABLE1))
has_af = sorted(n for n in TABLE2 if (glass[SHEET[n]]["data"].step.str.endswith("mT")).any())
print(f"Glass specimens of Table 2 found with complete measurements: {len(found_glass)} of {len(TABLE2)}")
print(f"Crystals of Table 1 found with replicate measurements: {len(found_crystals)} of {len(TABLE1)}")
print("Sheets or files not named in the tables:", extra)
print(f"Glass specimens with AF steps: {len(has_af)} {has_af}; without: {len(TABLE2) - len(has_af)}")
print("Same set as the Table 2 marks:", set(has_af) == set(TABLE2_AF_MARK))
print("Crystals with more than one NRM, 590.0 and 590.3 measurement:",
{n: [int((crystals[n]['data'].step == s).sum()) for s in ('0', '590', '590.3')] for n in TABLE1})
Glass specimens of Table 2 found with complete measurements: 22 of 22
Crystals of Table 1 found with replicate measurements: 12 of 12
Sheets or files not named in the tables: []
Glass specimens with AF steps: 10 ['ss19', 'ss22', 'ss27', 'ss34', 'ss38', 'ss41', 'ss43', 'ss45', 'ss46', 'ss48']; without: 12
Same set as the Table 2 marks: True
Crystals with more than one NRM, 590.0 and 590.3 measurement: {'12021g2': [2, 3, 3], '12021g3': [2, 3, 3], '12040g1': [2, 3, 3], '12040g2': [3, 3, 3], '12040g3': [2, 3, 3], '12053g1': [2, 3, 3], '12053g2': [3, 3, 3], '14053g1': [2, 3, 3], '14053g2': [2, 3, 3], '14053g3': [2, 4, 3], '71055g1': [1, 3, 3], '71055g3': [1, 3, 3]}
Result: reproduced. All 22 glass specimens of Table 2 (s11 as sheet ss11) and all 12 crystals of Table 1 are in the deposit with steps, moments, declinations and inclinations, and no other specimen is. The 10 glass specimens with AF steps are exactly those carrying the second mark in Table 2. Two crystals (71055g1, 71055g3) have a single NRM measurement. Four glass sheets hold scaled copies of one measurement in place of replicates (Data check).
glass-success-rate¶
Of the Apollo 64455 glass subsamples given thermal paleointensity experiments, 12% passed the selection criteria (3 of 25, 13 run without and 12 with AF pre-treatment).
Criterion: rate = 3 / (3 + failed specimens found here); the AF split is counted from AF steps in the deposit plus the Methods' successful ones (one without, two with AF). Holds when 22 failed specimens are found and 10 have AF steps.
n_failed = len(found_glass)
n_af = len(has_af)
rate = 100 * 3 / (3 + n_failed)
print(f"Failed specimens in the deposit: {n_failed}; with AF steps {n_af}, without {n_failed - n_af}")
print(f"Success rate 3 / (3 + {n_failed}) = {rate:.1f}% (expected 12 +- 1.2)")
print(f"Without AF: {n_failed - n_af} + 1 = {n_failed - n_af + 1} (expected 13); with AF: {n_af} + 2 = {n_af + 2} (expected 12)")
Failed specimens in the deposit: 22; with AF steps 10, without 12 Success rate 3 / (3 + 22) = 12.0% (expected 12 +- 1.2) Without AF: 12 + 1 = 13 (expected 13); with AF: 10 + 2 = 12 (expected 12)
Result: reproduced. 22 failed specimens give 3/25 = 12.0%, and the AF steps split them 12 and 10, which with the Methods' successful specimens gives 13 without and 12 with AF pre-treatment.
table2-fit-statistics¶
For the nine failed glass specimens with a line fit (ss19, ss43, ss22, ss41, ss6, ss9, ss10, ss30, ss34), the number of points, temperature range, MAD, R2 and pTRM check differences are as in Table 2, which lists the criteria each failed.
Criterion: per specimen, the paired steps in the Table 2 range, NRM against TRM gained, R2 of the least-squares line, MAD under readings (a) free, (b) anchored, (c) free on the three-point running mean, and pTRM check percentages under readings (a) over the original pTRM, (b) over the length of the fitted segment, with the check's pTRM referred to (i) the field-off step just before it or (ii) the field-off step at its temperature. A value matches when one reading puts it within tolerance. Reproduced when n matches for all nine and at least 80% of the MAD, R2 and pTRM values match; partially for 50 to 80%.
The paper does not say how the repeated measurements of a step are combined. Both the mean vector and the first measurement are run; the mean is adopted as an assumption (see the end of this section).
T2 = { # Table 2: n, T range, MAD, R2, pTRM checks {T: %}, failure criteria
"ss19": (4, 175, 450, 30.1, 0.94, {250: 30, 450: 9}, "pTRM, MAD, n"),
"ss43": (3, 250, 400, 10.0, 0.90, {250: 85}, "pTRM, n"),
"ss22": (4, 250, 500, 35.4, 0.38, {250: 22, 450: 13}, "pTRM, MAD, R2"),
"ss41": (3, 500, 580, 24.0, 0.38, {250: 4, 400: 24}, "pTRM, MAD, R2, n"),
"ss6": (3, 175, 360, 21.0, 0.91, {}, "MAD, R2, n, no pTRM"),
"ss9": (4, 250, 550, 37.6, 0.74, {250: 21}, "pTRM, MAD, R2"),
"ss10": (3, 360, 550, 42.9, 0.15, {250: 83}, "pTRM, MAD, R2, n"),
"ss30": (4, 175, 500, 12.5, 0.59, {450: 10}, "R2"),
"ss34": (5, 100, 450, 12.1, 0.06, {250: 9}, "R2"),
}
TOL_MAD = {"ss19": 3.0, "ss43": 1.0, "ss22": 3.5, "ss41": 2.4, "ss6": 2.1, "ss9": 3.8, "ss10": 4.3, "ss30": 1.25, "ss34": 1.2}
def fit_stats(how, nrm="vector"):
rows, checks = [], []
for s, (n, t0, t1, mad, rr, pt, fail) in T2.items():
df = glass[s]["data"]
a = A.in_range(A.arai(df, how, nrm), t0, t1)
m = A.mads(df, t0, t1, how)
r = A.r2(a.trm, a.nrm)
rows.append({"specimen": s, "how": how, "n": len(a), "n_T2": n, "steps": ",".join(str(int(t)) for t in a.t),
"R2": r, "R2_T2": rr, "R2_ok": abs(r - rr) <= 0.1 * rr + 1e-9,
**{f"MAD_{k}": v for k, v in m.items()}, "MAD_T2": mad,
"MAD_ok_any": any(abs(v - mad) <= TOL_MAD[s] + 1e-9 for v in m.values()),
"MAD_ok_a": abs(m["a_free"] - mad) <= TOL_MAD[s] + 1e-9})
p = A.ptrm_checks(df, how, (t0, t1), nrm)
for T, pct in pt.items():
q = p[p.t == T]
vals = {f"{col[0]}_{ref.split('_')[0]}": float(q[q.ref == ref][col].iloc[0])
for col in ("a_pct", "b_pct") for ref in ("i_prev", "ii_same")}
checks.append({"specimen": s, "how": how, "T": T, "printed": pct, **{k: round(v, 1) for k, v in vals.items()},
**{f"ok_{k}": abs(v / pct - 1) <= 0.1 + 1e-9 for k, v in vals.items()},
"prev_off_T": int(q.prev_off_t.iloc[0])})
rows, checks = pd.DataFrame(rows), pd.DataFrame(checks)
checks["ok_any"] = checks[[c for c in checks if c.startswith("ok_")]].any(axis=1)
return rows, checks
fits, chks = {}, {}
for how in ("mean", "first"):
fits[how], chks[how] = fit_stats(how)
print(f"--- replicates combined as: {how}")
print(fits[how].round(2).to_string(index=False))
print(chks[how].to_string(index=False))
--- replicates combined as: mean
specimen how n n_T2 steps R2 R2_T2 R2_ok MAD_a_free MAD_b_anchored MAD_c_smoothed MAD_T2 MAD_ok_any MAD_ok_a
ss19 mean 4 4 175,250,360,450 0.96 0.94 True 30.10 20.68 14.89 30.1 True True
ss43 mean 3 3 250,360,400 0.90 0.90 True 10.25 27.00 14.49 10.0 True True
ss22 mean 4 4 250,360,450,500 0.38 0.38 True 35.36 39.19 21.38 35.4 True True
ss41 mean 3 3 500,550,580 0.38 0.38 True 24.06 26.25 14.45 24.0 True True
ss6 mean 3 3 175,250,360 0.89 0.91 True 21.18 34.71 18.49 21.0 True True
ss9 mean 4 4 250,360,450,550 0.89 0.74 False 36.77 27.22 8.33 37.6 True True
ss10 mean 3 3 360,450,550 0.23 0.15 False 33.33 23.94 18.13 42.9 False False
ss30 mean 4 4 175,250,450,500 0.60 0.59 True 11.49 14.39 8.04 12.5 True True
ss34 mean 5 5 100,175,250,360,450 0.03 0.06 False 11.95 11.01 10.70 12.1 True True
specimen how T printed a_i a_ii b_i b_ii ok_a_i ok_a_ii ok_b_i ok_b_ii prev_off_T ok_any
ss19 mean 250 30 63.6 105.9 23.7 39.5 False False False False 360 False
ss19 mean 450 9 8.7 12.1 7.4 10.3 True False False False 500 True
ss43 mean 250 85 24.2 6.0 82.9 20.6 False False True False 360 True
ss22 mean 250 22 11.2 21.6 24.3 47.1 False True False False 360 True
ss22 mean 450 13 4.5 12.7 14.2 39.6 False True True False 500 True
ss41 mean 250 4 2.9 58.6 3.8 76.6 False False True False 360 True
ss41 mean 400 24 16.1 1.5 23.4 2.1 False False True False 450 True
ss9 mean 250 21 21.4 16.3 46.7 35.6 True False False False 450 True
ss10 mean 250 83 523.8 440.2 76.1 63.9 False False True False 450 True
ss30 mean 450 10 8.2 26.2 13.7 43.8 False False False False 500 False
ss34 mean 250 9 9.0 26.9 18.3 54.9 True False False False 360 True
--- replicates combined as: first
specimen how n n_T2 steps R2 R2_T2 R2_ok MAD_a_free MAD_b_anchored MAD_c_smoothed MAD_T2 MAD_ok_any MAD_ok_a
ss19 first 4 4 175,250,360,450 0.96 0.94 True 31.13 22.04 15.88 30.1 True True
ss43 first 3 3 250,360,400 0.99 0.90 False 9.87 30.82 14.11 10.0 True True
ss22 first 4 4 250,360,450,500 0.44 0.38 False 33.17 36.90 24.30 35.4 True True
ss41 first 3 3 500,550,580 0.16 0.38 False 23.51 26.62 15.32 24.0 True True
ss6 first 3 3 175,250,360 0.09 0.91 False 24.16 33.20 19.32 21.0 True False
ss9 first 4 4 250,360,450,550 0.84 0.74 False 30.67 28.53 7.93 37.6 False False
ss10 first 3 3 360,450,550 0.64 0.15 False 22.09 31.02 35.40 42.9 False False
ss30 first 4 4 175,250,450,500 0.98 0.59 False 19.09 18.51 7.72 12.5 False False
ss34 first 5 5 100,175,250,360,450 0.14 0.06 False 11.84 11.40 10.51 12.1 True True
specimen how T printed a_i a_ii b_i b_ii ok_a_i ok_a_ii ok_b_i ok_b_ii prev_off_T ok_any
ss19 first 250 30 29.3 87.2 12.9 38.4 True False False False 360 True
ss19 first 450 9 6.7 23.0 5.7 19.8 False False False False 500 False
ss43 first 250 85 25.2 6.4 86.1 22.0 False False True False 360 True
ss22 first 250 22 14.5 20.7 35.2 50.1 False True False False 360 True
ss22 first 450 13 11.9 15.4 38.1 49.3 True False False False 500 True
ss41 first 250 4 23.5 68.9 27.1 79.4 False False False False 360 False
ss41 first 400 24 13.3 1.7 20.3 2.5 False False False False 450 False
ss9 first 250 21 23.9 17.1 47.2 33.7 False False False False 450 False
ss10 first 250 83 452.1 364.2 62.4 50.3 False False False False 450 False
ss30 first 450 10 24.4 29.6 50.1 60.8 False False False False 500 False
ss34 first 250 9 14.2 21.8 15.6 23.9 False False False False 360 False
def score(how):
f, c = fits[how], chks[how]
n_ok = (f.n == f.n_T2).all()
mad_ok, r2_ok, pt_ok = int(f.MAD_ok_any.sum()), int(f.R2_ok.sum()), int(c.ok_any.sum())
total = len(f) * 2 + len(c)
return n_ok, mad_ok, r2_ok, pt_ok, total
for how in ("mean", "first"):
n_ok, mad_ok, r2_ok, pt_ok, total = score(how)
print(f"{how}: n matches for all nine {n_ok}; MAD {mad_ok}/9, R2 {r2_ok}/9, pTRM {pt_ok}/{len(chks[how])}; "
f"together {mad_ok + r2_ok + pt_ok}/{total} = {100 * (mad_ok + r2_ok + pt_ok) / total:.0f}%")
# How often a single reading matches, so that the any-reading count can be judged
f, c = fits["mean"], chks["mean"]
for k in ("a_free", "b_anchored", "c_smoothed"):
ok = sum(abs(f[f"MAD_{k}"] - f.MAD_T2) <= f.specimen.map(TOL_MAD) + 1e-9)
print(f"mean, MAD reading {k} alone: {ok}/9 within tolerance; median |difference| "
f"{np.median(np.abs(f[f'MAD_{k}'] - f.MAD_T2)):.1f} deg")
for k in ("a_i", "a_ii", "b_i", "b_ii"):
print(f"mean, pTRM reading {k} alone: {int(c[f'ok_{k}'].sum())}/{len(c)} within 10% of the printed value")
print("Values outside tolerance under every reading (mean):")
print(f.loc[~f.MAD_ok_any, ["specimen", "MAD_a_free", "MAD_b_anchored", "MAD_c_smoothed", "MAD_T2"]].round(1).to_string(index=False))
print(f.loc[~f.R2_ok, ["specimen", "R2", "R2_T2"]].round(3).to_string(index=False))
print(c.loc[~c.ok_any, ["specimen", "T", "printed", "a_i", "a_ii", "b_i", "b_ii"]].to_string(index=False))
mean: n matches for all nine True; MAD 8/9, R2 6/9, pTRM 9/11; together 23/29 = 79%
first: n matches for all nine True; MAD 6/9, R2 1/9, pTRM 4/11; together 11/29 = 38%
mean, MAD reading a_free alone: 8/9 within tolerance; median |difference| 0.2 deg
mean, MAD reading b_anchored alone: 2/9 within tolerance; median |difference| 9.4 deg
mean, MAD reading c_smoothed alone: 0/9 within tolerance; median |difference| 9.6 deg
mean, pTRM reading a_i alone: 3/11 within 10% of the printed value
mean, pTRM reading a_ii alone: 2/11 within 10% of the printed value
mean, pTRM reading b_i alone: 5/11 within 10% of the printed value
mean, pTRM reading b_ii alone: 0/11 within 10% of the printed value
Values outside tolerance under every reading (mean):
specimen MAD_a_free MAD_b_anchored MAD_c_smoothed MAD_T2
ss10 33.3 23.9 18.1 42.9
specimen R2 R2_T2
ss9 0.890 0.74
ss10 0.227 0.15
ss34 0.032 0.06
specimen T printed a_i a_ii b_i b_ii
ss19 250 30 63.6 105.9 23.7 39.5
ss30 450 10 8.2 26.2 13.7 43.8
# Failure criteria against the stated thresholds (n >= 4, pTRM within 15%, R2 about 0.9); reported, not judged
f, c = fits["mean"], chks["mean"]
rep = []
for s, (n, t0, t1, mad, rr, pt, fail) in T2.items():
fr = f[f.specimen == s].iloc[0]
cs = c[c.specimen == s]
rep.append({"specimen": s, "Table 2 failure": fail, "n<4": fr.n < 4, "R2<0.9": round(fr.R2, 2) < 0.9,
"pTRM>15% (reading a_i)": bool((cs.a_i > 15).any()) if len(cs) else None, "MAD (mean, a)": round(fr.MAD_a_free, 1)})
print(pd.DataFrame(rep).to_string(index=False))
print("ss6 250.1 holds 19 rows in two groups of directions:")
d = glass["ss6"]["data"].query("step == '250.1'")
print(d.assign(group=np.where(d.inc > 30, "inc>30", "inc<30")).groupby("group")[["moment", "dec", "inc"]].agg(["count", "mean"]).round(10))
specimen Table 2 failure n<4 R2<0.9 pTRM>15% (reading a_i) MAD (mean, a)
ss19 pTRM, MAD, n False False True 30.1
ss43 pTRM, n True False True 10.3
ss22 pTRM, MAD, R2 False True False 35.4
ss41 pTRM, MAD, R2, n True True True 24.1
ss6 MAD, R2, n, no pTRM True True None 21.2
ss9 pTRM, MAD, R2 False True True 36.8
ss10 pTRM, MAD, R2, n True True True 33.3
ss30 R2 False True False 11.5
ss34 R2 False True False 12.0
ss6 250.1 holds 19 rows in two groups of directions:
moment dec inc
count mean count mean count mean
group
inc<30 9 1.510000e-08 9 111.611111 9 -12.533333
inc>30 10 1.450000e-08 10 134.460000 10 70.850000
# Scatter of the repeated measurements of each step for the specimens whose MAD or R2 differ from
# Table 2 (ss9, ss10) and for ss6: moments in 1e-9 emu, the relative spread of the moments, the largest
# angle of a replicate from the step's mean direction, and the resultant length over n (1 when all agree)
def step_scatter(df):
rows = []
for step, g in df.groupby("step", sort=False):
X = A.to_xyz(g.moment, g.dec, g.inc)
m, mom = X.mean(axis=0), np.linalg.norm(X, axis=1)
rows.append({"step": step, "n": len(g), "mean vector (1e-9 emu)": np.linalg.norm(m) * 1e9,
"moment sd/mean": mom.std(ddof=1) / mom.mean() if len(g) > 1 else np.nan,
"max angle to mean": max(A.angle(x, m) for x in X),
"R/n": np.linalg.norm((X / mom[:, None]).sum(axis=0)) / len(g)})
return pd.DataFrame(rows)
for s in ("ss9", "ss10", "ss6"):
print(f"{s} (Table 2 range {T2[s][1]} to {T2[s][2]} C)")
print(step_scatter(glass[s]["data"]).round(3).to_string(index=False))
# The score with one pTRM reading held fixed (MAD and R2 as above), and with ss6 counted as a miss
f, c = fits["mean"], chks["mean"]
base = int(f.MAD_ok_any.sum()) + int(f.R2_ok.sum())
for k in ("a_i", "a_ii", "b_i", "b_ii"):
tot = base + int(c[f"ok_{k}"].sum())
print(f"pTRM reading {k} fixed: {tot}/29 = {100 * tot / 29:.0f}%")
tot = base + int(c.ok_any.sum()) - int(f.set_index("specimen").loc["ss6", "R2_ok"])
print(f"ss6's R2 counted as a miss: {tot}/29 = {100 * tot / 29:.0f}%")
print("Synthetic sheets among the nine:", [s for s in T2 if s in ("g2", "ss8", "ss39", "ss48")] or "none, so the score is the same without them")
ss9 (Table 2 range 250 to 550 C)
step n mean vector (1e-9 emu) moment sd/mean max angle to mean R/n
0 5 602.619 0.002 0.145 1.000
100 9 7.555 0.048 4.366 0.999
100.1 10 57.861 0.013 1.061 1.000
175 10 5.459 0.088 6.324 0.998
175.1 10 53.126 0.011 0.523 1.000
250 6 7.329 0.049 16.813 0.988
250.1 5 56.712 0.009 0.765 1.000
360 11 4.344 0.087 56.939 0.951
360.1 5 69.110 0.006 0.498 1.000
450 10 4.355 0.085 9.377 0.997
250.2 6 65.502 0.004 4.568 0.999
450.1 10 76.590 0.010 1.203 1.000
550 13 3.907 0.103 30.054 0.969
550.1 6 80.312 0.006 0.583 1.000
ss10 (Table 2 range 360 to 550 C)
step n mean vector (1e-9 emu) moment sd/mean max angle to mean R/n
0 5 54.170 0.010 0.759 1.000
100 10 2.112 0.098 24.677 0.971
100.1 9 2.426 0.265 55.111 0.935
175 8 2.239 0.087 20.374 0.976
175.1 10 1.862 0.301 58.562 0.926
250 16 1.190 0.283 112.636 0.799
250.1 15 0.598 0.368 76.906 0.795
360 7 3.054 0.162 52.591 0.821
360.1 5 5.667 0.013 6.620 0.998
450 11 2.224 0.244 40.129 0.958
250.2 4 5.937 0.063 5.713 0.997
450.1 6 3.434 0.107 16.788 0.979
550 11 2.119 0.095 71.053 0.915
550.1 5 10.878 0.045 0.719 1.000
ss6 (Table 2 range 175 to 360 C)
step n mean vector (1e-9 emu) moment sd/mean max angle to mean R/n
0 9 107.032 0.003 0.428 1.000
100 10 5.688 0.043 12.238 0.993
100.1 10 10.719 0.031 3.181 0.999
175 10 4.130 0.050 16.460 0.994
175.1 10 10.893 0.038 3.739 0.999
250 10 3.996 0.053 12.337 0.988
250.1 19 10.912 0.027 47.492 0.738
360 10 3.335 0.134 12.208 0.992
360.1 11 13.492 0.033 2.106 1.000
pTRM reading a_i fixed: 17/29 = 59%
pTRM reading a_ii fixed: 16/29 = 55%
pTRM reading b_i fixed: 19/29 = 66%
pTRM reading b_ii fixed: 14/29 = 48%
ss6's R2 counted as a miss: 22/29 = 76%
Synthetic sheets among the nine: none, so the score is the same without them
The deposit's SupplFig1 - other Thelliers March 2024.pdf plots NRM against TRM gained, with the
pTRM checks, for 18 of the glass specimens (all but ss19, ss43, ss22 and ss41, which are in the paper's
Figures 1 and 2). It shows how the authors combined the repeated measurements. The cell below reads
the plotted points from the PDF's vector drawing (apollo.suppfig1_points) and compares them with
each way of combining the rows: TRM gained from the mean vectors or the first measurements, and NRM as
the moment of the mean vector, the mean of the rows' moments, or the first row's moment. It compares
the plotted pTRM checks with the check's TRM referred to (i) the field-off step before it or (ii) the
field-off step at its temperature, and then scores Table 2 again with NRM read as the figure plots it.
fig1 = next(f for f in files if f["name"] == A.SUPPFIG1_FILE)
sf = A.suppfig1_points(A.read_bytes(fig1["id"]))
print(f"{A.SUPPFIG1_FILE}: computed_md5 {fig1['computed_md5']}; {len(sf)} pages, sheets {list(sf)}")
rows, unmatched = [], []
for s, v in sf.items():
df = glass[s]["data"]
am, af, ams = (A.arai(df, "mean").set_index("t"), A.arai(df, "first").set_index("t"),
A.arai(df, "mean", "moments").set_index("t"))
cand = {"TRM, mean vectors": ("x", am.trm), "TRM, first": ("x", af.trm), "NRM, |mean vector|": ("y", am.nrm),
"NRM, mean of moments": ("y", ams.nrm), "NRM, first": ("y", af.nrm)}
err = {k: [] for k in cand}
for t, x, y in v["points"]:
if t not in am.index:
unmatched.append((s, t))
continue
for k, (axis, c) in cand.items():
err[k].append(abs((x if axis == "x" else y) / c[t] - 1))
rows.append({"sheet": s, "points": len(v["points"]), "synthetic sheet": s in ("g2", "ss8", "ss39", "ss48"),
**{k: f"{100 * np.median(e):.1f} / {100 * np.max(e):.1f}" for k, e in err.items()}})
print("Difference between the plotted points and each reading, median / largest, % of the value:")
print(pd.DataFrame(rows).to_string(index=False))
print("Points whose temperature label has no paired step in the sheet:", unmatched)
rows = []
for s, v in sf.items():
p = A.ptrm_checks(glass[s]["data"], "mean")
for _, x, y in v["checks"]:
q = p.loc[(p.check_trm / x - 1).abs().idxmin()]
ref = {r: p[(p.t == q.t) & (p.ref == r)].check_trm.iloc[0] for r in ("i_prev", "ii_same")}
rows.append({"sheet": s, "plotted check TRM (1e-9 emu)": round(x * 1e9, 3), "check T": int(q.t),
"(i) field-off before it": round(ref["i_prev"] * 1e9, 3), "(ii) field-off at T": round(ref["ii_same"] * 1e9, 3),
"nearer": q.ref, "within 1% of (i)": abs(x / ref["i_prev"] - 1) < 0.01})
tri = pd.DataFrame(rows)
print(tri.to_string(index=False))
print(f"Plotted checks within 1% of reference (i): {int(tri['within 1% of (i)'].sum())}/{len(tri)}")
# Table 2 scored again with NRM the mean of the rows' moments (MAD does not change)
fits["moments"], chks["moments"] = fit_stats("mean", nrm="moments")
f, c = fits["moments"], chks["moments"]
print(f[["specimen", "n", "n_T2", "R2", "R2_T2", "R2_ok"]].round(3).to_string(index=False))
print(c[["specimen", "T", "printed", "a_i", "a_ii", "b_i", "b_ii", "ok_any"]].to_string(index=False))
n_ok, mad_ok, r2_ok, pt_ok, total = score("moments")
print(f"NRM as mean of moments: n matches for all nine {n_ok}; MAD {mad_ok}/9, R2 {r2_ok}/9, pTRM {pt_ok}/11 "
f"under any reading; together {mad_ok + r2_ok + pt_ok}/{total} = {100 * (mad_ok + r2_ok + pt_ok) / total:.0f}%")
pt_i = int((c.ok_a_i | c.ok_b_i).sum())
print(f"pTRM checks with reference (i), as the figure plots them: {pt_i}/11; together "
f"{mad_ok + r2_ok + pt_i}/{total} = {100 * (mad_ok + r2_ok + pt_i) / total:.0f}%")
for s in ("ss43", "ss19", "ss41", "ss22"):
a = A.in_range(A.arai(glass[s]["data"], "mean", "moments"), *T2[s][1:3])
print(f"{s}: R2 {A.r2(a.trm, a.nrm):.3f} (Table 2 {T2[s][4]}), nominal paleointensity |slope| x 18 = {abs(A.slope(a.trm, a.nrm)) * 18:.1f} uT")
2026-10-27 22:26:51.917057610 [W:onnxruntime:Default, device_discovery.cc:95 GetPciBusId] Skipping pci_bus_id for PCI path at "/sys/devices/LNXSYSTM:00/LNXSYBUS:00/ACPI0004:00/MSFT1000:00/5620e0c7-8062-4dce-aeb7-520c7ef76171" because filename "5620e0c7-8062-4dce-aeb7-520c7ef76171" did not match expected pattern of [0-9a-f]+:[0-9a-f]+:[0-9a-f]+[.][0-9a-f]+
SupplFig1 - other Thelliers March 2024.pdf: computed_md5 ee2ba78e902ca25e978606346abd5bc8; 18 pages, sheets ['g2', 'ss6', 'ss8', 'ss9', 'ss10', 'ss11', 'ss18', 'ss24', 'ss27', 'ss29', 'ss30', 'ss34', 'ss36', 'ss38', 'ss39', 'ss45', 'ss46', 'ss48']
Difference between the plotted points and each reading, median / largest, % of the value:
sheet points synthetic sheet TRM, mean vectors TRM, first NRM, |mean vector| NRM, mean of moments NRM, first
g2 4 True 1.8 / 5.2 6.2 / 23.1 1.5 / 2.7 1.1 / 2.0 2.5 / 2.9
ss6 4 False 0.0 / 0.0 6.1 / 38.3 0.7 / 1.3 0.0 / 0.0 4.9 / 8.0
ss8 5 True 1.9 / 10.7 2.2 / 8.9 0.9 / 1.7 0.6 / 1.2 1.0 / 1.0
ss9 6 False 0.0 / 0.0 1.4 / 3.2 0.8 / 4.8 0.0 / 0.0 3.7 / 8.6
ss10 6 False 0.1 / 11.2 11.9 / 90.5 5.8 / 35.8 0.0 / 5.5 11.6 / 20.0
ss11 6 False 0.0 / 0.0 4.5 / 28.2 2.3 / 22.1 0.0 / 0.0 7.5 / 46.5
ss18 5 False 0.0 / 0.0 1.3 / 3.5 0.4 / 6.5 0.0 / 0.0 5.0 / 71.2
ss24 6 False 0.0 / 0.0 12.8 / 19.0 3.5 / 17.6 0.0 / 0.0 13.6 / 46.7
ss27 6 False 6.3 / 49.3 8.9 / 20.5 1.4 / 2.7 1.4 / 2.6 3.5 / 10.1
ss29 6 False 0.0 / 0.0 8.2 / 25.3 6.3 / 29.5 0.0 / 0.0 9.0 / 205.9
ss30 5 False 0.0 / 0.0 13.5 / 21.5 7.1 / 101.6 0.0 / 0.1 4.1 / 22.6
ss34 6 False 4.7 / 11.1 5.3 / 12.9 3.1 / 6.4 3.3 / 6.7 3.5 / 13.1
ss36 4 False 0.0 / 0.0 9.5 / 12.5 6.6 / 14.4 0.0 / 0.0 9.8 / 32.4
ss38 4 False 0.1 / 0.3 18.5 / 35.3 0.0 / 0.1 0.1 / 3.3 3.5 / 16.0
ss39 2 True 1.7 / 3.3 2.9 / 3.1 0.6 / 1.0 0.5 / 1.0 1.6 / 2.0
ss45 4 False 0.0 / 0.0 0.3 / 8.8 0.0 / 0.0 0.1 / 0.4 1.9 / 4.7
ss46 3 False 0.0 / 0.1 1.6 / 1.9 0.0 / 0.1 0.0 / 0.6 5.5 / 7.4
ss48 2 True 1.0 / 1.7 3.0 / 3.2 1.8 / 2.8 1.6 / 2.6 2.1 / 2.2
Points whose temperature label has no paired step in the sheet: [('ss27', 210.0)]
sheet plotted check TRM (1e-9 emu) check T (i) field-off before it (ii) field-off at T nearer within 1% of (i)
ss9 68.529 250 68.533 65.650 i_prev True
ss10 5.969 250 5.969 5.169 i_prev True
ss11 15.950 250 15.950 15.662 i_prev True
ss18 25.890 250 25.890 31.307 i_prev True
ss24 3.882 250 3.882 6.301 i_prev True
ss24 10.430 450 10.431 6.382 i_prev True
ss27 4.817 305 2.753 14.107 i_prev False
ss27 17.655 400 17.771 33.953 i_prev True
ss29 10.721 450 10.720 10.479 i_prev True
ss30 7.146 450 7.146 8.337 i_prev True
ss34 12.403 250 12.417 8.334 i_prev True
ss38 7.980 250 7.980 7.116 i_prev True
Plotted checks within 1% of reference (i): 11/12
specimen n n_T2 R2 R2_T2 R2_ok
ss19 4 4 0.956 0.94 True
ss43 3 3 0.909 0.90 True
ss22 4 4 0.370 0.38 True
ss41 3 3 0.382 0.38 True
ss6 3 3 0.906 0.91 True
ss9 4 4 0.916 0.74 False
ss10 3 3 0.152 0.15 True
ss30 4 4 0.594 0.59 True
ss34 5 5 0.032 0.06 False
specimen T printed a_i a_ii b_i b_ii ok_any
ss19 250 30 63.6 105.9 23.7 39.5 False
ss19 450 9 8.7 12.1 7.4 10.3 True
ss43 250 85 24.2 6.0 83.8 20.8 True
ss22 250 22 11.2 21.6 24.3 47.1 True
ss22 450 13 4.5 12.7 14.2 39.7 True
ss41 250 4 2.9 58.6 3.8 76.6 True
ss41 400 24 16.1 1.5 23.4 2.1 True
ss9 250 21 21.4 16.3 46.7 35.6 True
ss10 250 83 523.8 440.2 75.9 63.8 True
ss30 450 10 8.2 26.2 14.1 45.1 False
ss34 250 9 9.0 26.9 18.3 54.9 True
NRM as mean of moments: n matches for all nine True; MAD 8/9, R2 7/9, pTRM 9/11 under any reading; together 24/29 = 83%
pTRM checks with reference (i), as the figure plots them: 8/11; together 23/29 = 79%
ss43: R2 0.909 (Table 2 0.9), nominal paleointensity |slope| x 18 = 9.2 uT
ss19: R2 0.956 (Table 2 0.94), nominal paleointensity |slope| x 18 = 23.0 uT
ss41: R2 0.382 (Table 2 0.38), nominal paleointensity |slope| x 18 = 1.2 uT
ss22: R2 0.370 (Table 2 0.38), nominal paleointensity |slope| x 18 = 3.4 uT
Result: partially reproduced, resting on the assumption that repeated measurements are averaged. With the mean vector of each step's repeated measurements, n matches Table 2 for all nine specimens, and 23 of the 29 MAD, R2 and pTRM values match under one of the named readings (79%, just under the 80% the criterion asks for): MAD 8 of 9, R2 6 of 9, pTRM checks 9 of 11. MAD reading (a), a free line, alone gives the 8 matching MADs, 6 of them within 0.3 deg of the printed value. ss10 differs in MAD (33.3 against 42.9 deg) and R2 (0.23 against 0.15); ss9 differs in R2 (0.89 against 0.74); and ss34's R2 is 0.03 against 0.06, outside a tolerance of 0.006. The matching pTRM checks need different readings for different specimens (no single reading matches more than 5 of 11), so no single pTRM definition is confirmed. With the first measurement of each step, 11 of 29 values match (38%), which would make the claim not reproduced; the mean is adopted because it brings the MADs of six specimens within 0.3 deg of the printed values. ss6's 250.1 step mixes two groups of directions (19 rows), which lowers its R2 under the mean; counted as a miss, the score is 22 of 29 (76%). Holding one pTRM reading fixed gives 14 to 19 of 29 (48 to 66%), so the outcome is partial or lower under any single pTRM definition, and none of the nine is a synthetic sheet. The replicate tables above show where ss9 and ss10 differ: all of ss10's field-off steps and its 100, 175 and 250 C field-on steps scatter by 20 to 113 deg within a step (resultant over n 0.79 to 0.98), and ss9's 360 and 550 C steps by 57 and 30 deg, so how such steps are combined moves its MAD and R2. Eleven steps of seven sheets, nine of them inside the specimen's Table 2 range, hold two groups of directions more than 45 deg apart, as ss6's 250.1 does (table in table2-no-fit); averaging such a step gives a vector without physical meaning.
The deposit's Supplementary Figure 1 supports the averaging. For 12 of its 14 specimens that are not synthetic sheets (all but ss27 and ss34, which fit no reading), the plotted points equal TRM gained from the mean vectors and NRM as the mean of the rows' moments, with a median difference of 0.1% or less; the first measurements differ from them by 0.3 to 18.5% (median). The figure plots ss6's 250 C step at its averaged value, so the authors' fit used that value too. 11 of the 12 plotted pTRM checks sit at the check's TRM referred to the field-off step before it, reference (i). With NRM read as the figure plots it, R2 matches for 7 of 9 (ss10 now 0.152 against 0.15; ss9 0.916 against 0.74 and ss34 still miss), and 24 of 29 values match under any of the criterion's readings (83%). The same figure shows that the checks use reference (i); holding that reference gives 8 of 11 pTRM checks and 23 of 29 (79%). The outcome stays partially reproduced: the figure's NRM reading crosses the 80% bar only with a pTRM reference the figure does not use. The plotted points of the four synthetic sheets differ from every combination of their deposited rows by 0.5% or more (median).
table2-no-fit¶
The other thirteen failed glass specimens (g2, ss8, s11, ss18, ss24, ss27, ss29, ss36, ss38, ss39, ss45, ss46, ss48) gave no fit, most because they did not acquire TRM in a regular fashion and some because their directions trend out of the origin.
Criterion: over each Table 2 range, TRM gained per step; regular acquisition read as (a) TRM gained increasing at every step, or (b) no step falling by more than 15%. Reproduced when, under one reading, at least 9 of the 11 specimens marked "no regular TRM acquisition" are irregular, and ss36 and ss39 have field-off vectors whose distance from the origin grows over the last steps. The pTRM checks are reported.
NOFIT = {"g2": (4, 100, 400, {}), "ss8": (5, 100, 450, {}), "s11": (6, 100, 550, {250: 22}), "ss18": (5, 100, 450, {250: 66}),
"ss24": (7, 100, 550, {250: 6, 450: 42}), "ss27": (6, 100, 550, {305: 1, 400: 54}), "ss29": (6, 100, 550, {450: 25}),
"ss36": (4, 100, 450, {450: 25}), "ss38": (4, 100, 360, {250: 9}), "ss39": (2, 100, 210, {}),
"ss45": (4, 100, 360, {250: 120}), "ss46": (3, 100, 250, {}), "ss48": (2, 100, 175, {})}
IRREGULAR_LABEL = [s for s in NOFIT if s not in ("ss36", "ss39")]
SYNTH = ["g2", "ss8", "ss39", "ss48"]
def regularity(how):
rows = []
for s, (n, t0, t1, pt) in NOFIT.items():
a = A.in_range(A.arai(glass[SHEET[s]]["data"], how), t0, t1)
trm = a.trm.to_numpy()
ratio = trm[1:] / trm[:-1]
rows.append({"specimen": s, "how": how, "n": len(a), "n_T2": n, "T": ",".join(str(int(t)) for t in a.t),
"TRM gained (1e-9 emu)": " ".join(f"{x * 1e9:.2f}" for x in trm),
"irregular_a": bool((ratio <= 1).any()), "irregular_b": bool((ratio < 0.85).any()),
"largest fall %": round(100 * (1 - ratio.min()), 1) if len(ratio) else None,
"|off| (1e-9 emu)": " ".join(f"{x * 1e9:.2f}" for x in a.nrm),
"synthetic sheet": s in SYNTH})
return pd.DataFrame(rows)
for how in ("mean", "first"):
r = regularity(how)
print(f"--- {how}")
print(r.to_string(index=False))
lab = r[r.specimen.isin(IRREGULAR_LABEL)]
print(f"{how}: irregular under (a) {int(lab.irregular_a.sum())}/11, under (b) {int(lab.irregular_b.sum())}/11; "
f"without the four synthetic sheets (a) {int(lab[~lab['synthetic sheet']].irregular_a.sum())}/{int((~lab['synthetic sheet']).sum())}, "
f"(b) {int(lab[~lab['synthetic sheet']].irregular_b.sum())}/{int((~lab['synthetic sheet']).sum())}")
for s in ("ss36", "ss39"):
a = A.in_range(A.arai(glass[s]["data"], how), *NOFIT[s][1:3])
print(f" {s}: distance of the field-off vector from the origin over the range "
f"{[f'{x * 1e9:.2f}' for x in a.nrm]} (1e-9 emu); grows over the last step: {bool(a.nrm.iloc[-1] > a.nrm.iloc[-2])}")
--- mean
specimen how n n_T2 T TRM gained (1e-9 emu) irregular_a irregular_b largest fall % |off| (1e-9 emu) synthetic sheet
g2 mean 4 4 100,210,305,400 4.99 3.81 2.89 1.80 True True 37.5 2.81 0.79 2.13 1.60 True
ss8 mean 5 5 100,175,250,360,450 2.72 0.60 1.29 3.33 3.98 True True 78.1 2.66 1.54 0.90 2.81 3.18 True
s11 mean 6 6 100,175,250,360,450,550 12.03 12.92 14.14 14.02 20.37 11.82 True True 42.0 3.26 2.09 3.83 3.29 2.36 3.35 False
ss18 mean 5 5 100,175,250,360,450 18.81 21.60 19.49 28.55 22.44 True True 21.4 3.84 5.84 2.02 9.46 9.51 False
ss24 mean 7 7 100,175,250,360,450,500,550 119.30 2.97 4.47 5.95 5.94 12.13 13.57 True True 97.5 120.92 3.94 1.96 2.17 3.14 2.69 4.02 False
ss27 mean 5 6 100,305,400,500,550 1.15 5.36 6.11 24.76 11.54 True True 53.4 17.43 35.04 37.23 26.62 15.58 False
ss29 mean 6 6 100,175,250,450,500,550 5.02 6.36 6.32 6.35 8.45 22.29 True False 0.6 1.20 3.18 3.65 1.42 3.60 1.38 False
ss36 mean 4 4 100,175,250,450 13.13 14.44 15.79 16.76 False False -6.2 1.59 1.14 1.77 6.00 False
ss38 mean 4 4 100,175,250,360 6.96 8.46 6.86 5.87 True True 18.9 22.12 14.70 9.20 1.71 False
ss39 mean 2 2 100,210 6.71 11.65 False False -73.7 2.31 9.67 True
ss45 mean 4 4 100,175,250,360 15.37 15.13 19.52 16.53 True True 15.3 8.12 6.76 9.35 4.87 False
ss46 mean 3 3 100,175,250 32.93 35.74 32.94 True False 7.8 5.12 19.36 12.43 False
ss48 mean 2 2 100,175 19.36 0.30 True True 98.4 4.29 0.42 True
mean: irregular under (a) 11/11, under (b) 9/11; without the four synthetic sheets (a) 8/8, (b) 6/8
ss36: distance of the field-off vector from the origin over the range ['1.59', '1.14', '1.77', '6.00'] (1e-9 emu); grows over the last step: True
ss39: distance of the field-off vector from the origin over the range ['2.31', '9.67'] (1e-9 emu); grows over the last step: True
--- first
specimen how n n_T2 T TRM gained (1e-9 emu) irregular_a irregular_b largest fall % |off| (1e-9 emu) synthetic sheet
g2 first 4 4 100,210,305,400 5.25 3.65 2.22 1.72 True True 39.0 2.93 0.81 2.25 1.57 True
ss8 first 5 5 100,175,250,360,450 2.70 0.72 1.29 3.38 4.38 True True 73.1 2.65 1.51 0.90 2.86 3.27 True
s11 first 6 6 100,175,250,360,450,550 11.77 13.17 15.19 19.53 20.03 13.29 True True 33.7 3.01 1.99 4.20 3.20 1.67 4.15 False
ss18 first 5 5 100,175,250,360,450 19.50 21.58 19.75 28.39 21.76 True True 23.3 4.96 5.59 1.21 9.16 9.28 False
ss24 first 7 7 100,175,250,360,450,500,550 118.12 3.58 5.52 7.08 6.00 11.06 13.69 True True 97.0 119.70 4.63 3.48 2.31 3.05 1.85 5.34 False
ss27 first 5 6 100,305,400,500,550 1.43 4.68 5.25 24.31 11.66 True True 52.0 18.20 35.03 36.61 26.65 14.54 False
ss29 first 6 6 100,175,250,450,500,550 5.76 8.52 6.23 6.57 8.41 26.36 True True 26.9 0.43 3.24 3.99 1.49 3.59 2.12 False
ss36 first 4 4 100,175,250,450 12.31 13.88 14.03 19.12 False False -1.1 1.84 1.82 1.82 6.35 False
ss38 first 4 4 100,175,250,360 6.05 10.15 5.69 4.34 True True 43.9 22.24 14.12 9.49 2.04 False
ss39 first 2 2 100,210 6.75 11.32 False False -67.7 2.29 9.39 True
ss45 first 4 4 100,175,250,360 15.35 15.19 19.55 18.13 True False 7.3 8.00 6.92 9.49 4.65 False
ss46 first 3 3 100,175,250 33.22 36.43 33.49 True False 8.1 4.77 19.28 13.16 False
ss48 first 2 2 100,175 20.07 0.30 True True 98.5 4.24 0.42 True
first: irregular under (a) 11/11, under (b) 9/11; without the four synthetic sheets (a) 8/8, (b) 6/8
ss36: distance of the field-off vector from the origin over the range ['1.84', '1.82', '1.82', '6.35'] (1e-9 emu); grows over the last step: True
ss39: distance of the field-off vector from the origin over the range ['2.29', '9.39'] (1e-9 emu); grows over the last step: True
# A null for the irregularity counts: how often 'irregular' arises from measurement scatter alone.
# For each step of the non-synthetic specimens, the spread of TRM gained across replicates
# (the field-on and field-off rows paired in order) is compared with the step-to-step falls.
rng = np.random.default_rng(2101)
rows = []
for s in IRREGULAR_LABEL:
if s in SYNTH:
continue
df = glass[SHEET[s]]["data"]
a = A.in_range(A.arai(df, "mean"), *NOFIT[s][1:3])
sd = []
for t in a.t:
on = df[df.step == f"{t:g}.1"]; off = df[df.step.isin([f"{t:g}", f"{t:g}.0"])]
k = min(len(on), len(off))
x = np.linalg.norm(A.to_xyz(on.moment[:k], on.dec[:k], on.inc[:k]) - A.to_xyz(off.moment[:k], off.dec[:k], off.inc[:k]), axis=1)
sd.append(x.std(ddof=1) / np.sqrt(k))
trm = a.trm.to_numpy()
rows.append({"specimen": s, "median SE of mean TRM gained / mean TRM gained": round(float(np.median(np.array(sd) / a.trm)), 3),
"largest fall / TRM gained before it": round(float(1 - (trm[1:] / trm[:-1]).min()), 3)})
se = pd.DataFrame(rows)
se["fall above median SE"] = se["largest fall / TRM gained before it"] > se["median SE of mean TRM gained / mean TRM gained"]
print(se.to_string(index=False))
lab = regularity("mean").set_index("specimen").loc[IRREGULAR_LABEL]
small = [s for s in se.specimen[~se["fall above median SE"]] if lab.loc[s, "irregular_a"]]
print(f"Irregular under (a) with falls smaller than the median SE counted as regular: "
f"{int(lab.irregular_a.sum()) - len(small)}/11 (removed: {small})")
print("pTRM checks of the no-fit specimens (mean, readings a and b over the Table 2 range), against Table 2:")
out = []
for s, (n, t0, t1, pt) in NOFIT.items():
p = A.ptrm_checks(glass[SHEET[s]]["data"], "mean", (t0, t1))
for T, pct in pt.items():
q = p[p.t == T]
if len(q) == 0:
out.append({"specimen": s, "T": T, "printed": pct, "note": "no check at this T in the sheet"})
continue
vals = {f"{c[0]}_{r.split('_')[0]}": round(float(q[q.ref == r][c].iloc[0]), 1) for c in ("a_pct", "b_pct") for r in ("i_prev", "ii_same")}
out.append({"specimen": s, "T": T, "printed": pct, **vals,
"any within 10%": any(abs(v / pct - 1) <= 0.1 for v in vals.values())})
out = pd.DataFrame(out)
print(out.to_string(index=False))
print(f"pTRM checks matching under some reading: {int(out['any within 10%'].fillna(False).sum())}/{len(out)}")
specimen median SE of mean TRM gained / mean TRM gained largest fall / TRM gained before it fall above median SE
s11 0.019 0.420 True
ss18 0.012 0.214 True
ss24 0.022 0.975 True
ss27 0.035 0.534 True
ss29 0.025 0.006 False
ss38 0.061 0.189 True
ss45 0.009 0.153 True
ss46 0.010 0.078 True
Irregular under (a) with falls smaller than the median SE counted as regular: 10/11 (removed: ['ss29'])
pTRM checks of the no-fit specimens (mean, readings a and b over the Table 2 range), against Table 2:
specimen T printed a_i a_ii b_i b_ii any within 10% note
s11 250 22 12.8 10.8 21.1 17.7 True NaN
ss18 250 66 32.9 60.7 54.2 100.0 True NaN
ss24 250 6 13.1 41.0 0.3 1.1 False NaN
ss24 450 42 75.5 7.4 2.7 0.3 False NaN
ss27 305 1 48.7 163.0 11.0 37.0 False NaN
ss27 400 54 190.8 455.7 49.3 117.8 True NaN
ss29 450 25 68.9 65.1 25.3 23.9 True NaN
ss36 450 25 NaN NaN NaN NaN NaN no check at this T in the sheet
ss38 250 9 16.3 3.7 9.3 2.1 True NaN
ss45 250 120 33.0 21.4 131.7 85.5 True NaN
pTRM checks matching under some reading: 6/10
# Steps of the glass sheets whose repeated measurements fall in two groups of directions, like ss6
# 250.1. Each step with 4 or more rows is split in two by the angle to the two most distant rows
# (then 10 rounds of reassignment to the group means); listed when both groups have 2 or more rows,
# their means are more than 45 deg apart and every row is within a third of that angle of its group's mean.
def two_groups(X):
u = X / np.linalg.norm(X, axis=1, keepdims=True)
i, j = np.unravel_index(np.argmin(u @ u.T), (len(u), len(u)))
lab = (u @ u[j] > u @ u[i]).astype(int)
for _ in range(10):
cen = [u[lab == k].sum(axis=0) for k in (0, 1)]
lab = (u @ cen[1] / np.linalg.norm(cen[1]) > u @ cen[0] / np.linalg.norm(cen[0])).astype(int)
sep = A.angle(cen[0], cen[1])
within = max(A.angle(x, cen[k]) for k in (0, 1) for x in u[lab == k])
return int((lab == 0).sum()), int((lab == 1).sum()), sep, within
rows = []
for n in TABLE2:
s = SHEET[n]
t0, t1 = (T2.get(n) or NOFIT[n])[1:3] if n in T2 or n in NOFIT else (None, None)
for step, g in glass[s]["data"].groupby("step", sort=False):
if len(g) < 4:
continue
n0, n1, sep, within = two_groups(A.to_xyz(g.moment, g.dec, g.inc))
if min(n0, n1) >= 2 and sep > 45 and within < sep / 3:
kind, t = A.step_kind(step)
rows.append({"specimen": n, "step": step, "rows": len(g), "groups": f"{n0}+{n1}", "apart (deg)": round(sep, 1),
"largest spread in a group (deg)": round(within, 1), "synthetic sheet": n in SYNTH,
"in Table 2 range": t is not None and kind in ("off", "on") and t0 <= t <= t1})
print(pd.DataFrame(rows).to_string(index=False))
specimen step rows groups apart (deg) largest spread in a group (deg) synthetic sheet in Table 2 range
ss6 250.1 19 10+9 84.9 5.6 False True
ss10 250 16 14+2 110.0 34.2 False False
ss10 360 7 5+2 75.5 15.2 False True
s11 360.1 13 10+3 92.3 4.2 False True
s11 550 9 7+2 85.2 15.7 False True
s11 550.1 9 7+2 85.3 2.6 False True
ss18 100 13 11+2 56.3 12.3 False True
ss24 0 6 4+2 49.3 0.2 False False
ss24 550 10 5+5 57.3 9.2 False True
ss29 175.1 13 10+3 62.8 8.5 False True
ss36 450.1 7 5+2 85.0 1.8 False True
Result: reproduced. Under reading (a), all 11 specimens marked "no regular TRM acquisition" acquire TRM irregularly (9 of 11 under reading b, and 8 of the 8 not from synthetic sheets under reading a); ss36 and ss39 have field-off vectors that move away from the origin over the last step of their ranges. Three of the 11 (g2, ss8, ss48), and ss39, are synthetic sheets. Reading (a) is weak: any fall in TRM gained counts. ss29's only fall, 0.6%, is below its median standard error of 2.5%; counting it as regular leaves 10 of 11 under reading (a), still at least 9. The trend out of the origin rests on ss36 and on ss39, a synthetic sheet. The pTRM checks of these specimens are reported only: 6 of 10 match under some reading, and ss36 has no check at 450 C in the sheet although Table 2 prints one; ss27 has 5 paired steps from 100 to 550 C where Table 2 gives 6.
fig1-linear-components¶
Glass subsamples ss43 and ss19, which failed the reliability criteria, show hints of linear directional components on orthogonal vector plots, from which nominal paleointensities can be calculated with the 18 uT laboratory field.
Criterion: over the Table 2 ranges, MAD and R2 as in table2-fit-statistics, and the nominal paleointensity as the absolute slope times 18 uT (Figure 1 legend). Holds when both MADs and R2 values match Table 2 under one reading.
for how in ("mean", "first"):
f = fits[how].set_index("specimen")
for s in ("ss43", "ss19"):
n, t0, t1 = T2[s][:3]
a = A.in_range(A.arai(glass[s]["data"], how), t0, t1)
b = A.slope(a.trm, a.nrm)
print(f"{how} {s}: MAD a/b/c {f.loc[s, 'MAD_a_free']:.1f}/{f.loc[s, 'MAD_b_anchored']:.1f}/{f.loc[s, 'MAD_c_smoothed']:.1f} "
f"(Table 2 {T2[s][3]}), R2 {f.loc[s, 'R2']:.3f} (Table 2 {T2[s][4]}), slope {b:.3f}, "
f"lab field in sheet {glass[s]['lab_field_uT']} uT, nominal paleointensity |slope| x 18 = {abs(b) * 18:.1f} uT")
fig, axes = plt.subplots(1, 4, figsize=(15, 3.5))
lim = []
for ax, s in zip(axes, ("ss43", "ss19", "ss41", "ss22")):
a = A.arai(glass[s]["data"], "mean")
t0, t1 = T2[s][1:3]
sel = (a.t >= t0) & (a.t <= t1)
ax.plot(a.trm * 1e9, a.nrm * 1e9, "o-", color="0.6")
ax.plot(a.trm[sel] * 1e9, a.nrm[sel] * 1e9, "o", color="k")
for t, x, y in zip(a.t, a.trm * 1e9, a.nrm * 1e9):
ax.annotate(f"{t:g}", (x, y), fontsize=7)
ax.set_title(f"{s} (mean of replicates)")
ax.set_xlabel("TRM gained (1e-9 emu)")
ax.set_ylabel("NRM remaining (1e-9 emu)")
lim.append((s, a.trm.min() * 1e9, a.trm.max() * 1e9, a.nrm.min() * 1e9, a.nrm.max() * 1e9))
plt.tight_layout()
plt.show()
print("Plotted ranges (x min, x max, y min, y max, 1e-9 emu); all points are inside the axes:")
for r in lim:
print(r[0], [round(float(v), 2) for v in r[1:]])
mean ss43: MAD a/b/c 10.3/27.0/14.5 (Table 2 10.0), R2 0.902 (Table 2 0.9), slope -0.536, lab field in sheet 18.0 uT, nominal paleointensity |slope| x 18 = 9.6 uT mean ss19: MAD a/b/c 30.1/20.7/14.9 (Table 2 30.1), R2 0.956 (Table 2 0.94), slope -1.278, lab field in sheet 18.0 uT, nominal paleointensity |slope| x 18 = 23.0 uT first ss43: MAD a/b/c 9.9/30.8/14.1 (Table 2 10.0), R2 0.993 (Table 2 0.9), slope -0.471, lab field in sheet 18.0 uT, nominal paleointensity |slope| x 18 = 8.5 uT first ss19: MAD a/b/c 31.1/22.0/15.9 (Table 2 30.1), R2 0.957 (Table 2 0.94), slope -1.318, lab field in sheet 18.0 uT, nominal paleointensity |slope| x 18 = 23.7 uT
Plotted ranges (x min, x max, y min, y max, 1e-9 emu); all points are inside the axes: ss43 [12.52, 15.75, 0.91, 5.16] ss19 [3.41, 12.46, 3.99, 19.48] ss41 [8.42, 50.62, 5.04, 198.4] ss22 [14.4, 44.21, 0.41, 78.84]
Result: reproduced, resting on the averaging assumption. With averaged replicates and a free line, ss43 has MAD 10.3 deg and R2 0.90 (Table 2: 10.0 and 0.90), and ss19 MAD 30.1 deg and R2 0.96 (30.1 and 0.94). The nominal paleointensities are 9.6 uT for ss43 and 23.0 uT for ss19 (8.5 and 23.7 uT with first measurements); the figure's values were not available to compare.
fig2-changing-field¶
Glass subsamples ss41 and ss22 show field-off directions that change over a temperature interval, consistent with a changing ambient field during magnetization, and curved NRM lost against TRM gained plots.
Criterion: Table 2 ranges (ss41 500 to 580 C, ss22 250 to 500 C). Readings (a) the angle between the first and last field-off direction of the range exceeds the MAD of a line fit over the range, (b) a great-circle fit has a smaller MAD than a line fit. Reproduced when both specimens change under one reading and R2 is within tolerance of 0.38.
for how in ("mean", "first"):
for s in ("ss41", "ss22"):
t0, t1 = T2[s][1:3]
seq = A.demag_sequence(glass[s]["data"], how)
r = seq[(seq.kind == "off") & (seq.t >= t0) & (seq.t <= t1)]
X = r[["x", "y", "z"]].to_numpy()
ang = A.angle(X[0], X[-1])
line_free, line_anch, plane = A.mad_line(X), A.mad_line(X, anchored=True), A.mad_plane(X)
dirs = [tuple(round(v, 1) for v in A.to_dir(x)[:2]) for x in X]
r2v = fits[how].set_index("specimen").loc[s, "R2"]
print(f"{how} {s}: steps {list(r.t.astype(int))}, directions (dec, inc) {dirs}")
print(f" first-last angle {ang:.1f} deg; line MAD free {line_free:.1f}, anchored {line_anch:.1f}; "
f"great-circle MAD {plane:.1f}; (a) {ang > line_free}, (b) {plane < line_free}; R2 {r2v:.3f} (0.38 +- 0.038)")
mean ss41: steps [500, 550, 580], directions (dec, inc) [(152.7, -62.2), (217.8, -44.8), (125.1, -33.9)] first-last angle 33.3 deg; line MAD free 24.1, anchored 26.3; great-circle MAD 15.1; (a) True, (b) True; R2 0.384 (0.38 +- 0.038) mean ss22: steps [250, 360, 450, 500], directions (dec, inc) [(45.2, 46.6), (148.6, -74.9), (143.7, 17.2), (360.0, -54.1)] first-last angle 107.8 deg; line MAD free 35.4, anchored 39.2; great-circle MAD 32.0; (a) True, (b) True; R2 0.378 (0.38 +- 0.038) first ss41: steps [500, 550, 580], directions (dec, inc) [(149.2, -63.1), (217.9, -46.0), (119.2, -30.2)] first-last angle 38.1 deg; line MAD free 23.5, anchored 26.6; great-circle MAD 14.4; (a) True, (b) True; R2 0.156 (0.38 +- 0.038) first ss22: steps [250, 360, 450, 500], directions (dec, inc) [(36.7, 52.3), (143.3, -77.9), (140.3, 5.2), (341.7, -33.7)] first-last angle 98.5 deg; line MAD free 33.2, anchored 36.9; great-circle MAD 26.4; (a) True, (b) True; R2 0.442 (0.38 +- 0.038)
# Curvature of the NRM lost against TRM gained plots, which the criterion does not test. A circle is
# fit by least squares (Kasa's algebraic fit) to the points of each Table 2 range, with NRM and TRM
# gained divided by the NRM of the first point; curvature is 1 / radius, positive when the points
# sag below the line through the first and last point. With three points (ss41, ss43) the circle is
# exact. ss19 and ss43 are given for comparison.
def curvature(x, y):
M = np.column_stack([x, y, np.ones_like(x)])
a, b, c = np.linalg.lstsq(M, x ** 2 + y ** 2, rcond=None)[0]
cx, cy = a / 2, b / 2
r = np.sqrt(c + cx ** 2 + cy ** 2)
mid = len(x) // 2
chord_y = y[0] + (y[-1] - y[0]) * (x[mid] - x[0]) / (x[-1] - x[0])
return np.sign(chord_y - y[mid]) / r
for s in ("ss41", "ss22", "ss43", "ss19"):
a = A.in_range(A.arai(glass[s]["data"], "mean"), *T2[s][1:3])
x, y = (a.trm / a.nrm.iloc[0]).to_numpy(), (a.nrm / a.nrm.iloc[0]).to_numpy()
print(f"{s}: steps {list(a.t.astype(int))}, NRM {y.round(3).tolist()}, TRM gained {x.round(3).tolist()}, "
f"curvature {curvature(x, y):.2f}, R2 {A.r2(a.trm, a.nrm):.3f}")
# How often readings (a) and (b) hold when there is no changing field: (1) directions drawn at random,
# (2) points along one line toward the origin with isotropic noise; n points as for ss41 (3 field-off
# directions in its range) and ss22 (4), and 5.
rng = np.random.default_rng(2201)
for npts in (3, 4, 5):
res = {"random": [], "noisy line": []}
for _ in range(5000):
R = rng.normal(size=(npts, 3)) * rng.uniform(0.3, 1, size=(npts, 1))
d = rng.normal(size=3); d /= np.linalg.norm(d)
L = np.linspace(1, 0.3, npts)[:, None] * d + rng.normal(scale=0.15, size=(npts, 3))
for k, X in (("random", R), ("noisy line", L)):
lf = A.mad_line(X)
res[k].append((A.angle(X[0], X[-1]) > lf, A.mad_plane(X) < lf, lf))
for k, v in res.items():
v = np.array(v)
print(f"{npts} vectors, {k}: (a) holds {100 * v[:, 0].mean():.0f}%, (b) {100 * v[:, 1].mean():.0f}%, "
f"median line MAD {np.median(v[:, 2]):.1f} deg")
ss41: steps [500, 550, 580], NRM [1.0, 0.995, 0.697], TRM gained [4.111, 6.595, 7.002], curvature -0.41, R2 0.384 ss22: steps [250, 360, 450, 500], NRM [1.0, 0.645, 0.822, 0.282], TRM gained [4.261, 4.712, 6.096, 6.182], curvature -1.12, R2 0.378 ss43: steps [250, 360, 400], NRM [1.0, 0.463, 0.348], TRM gained [4.787, 5.446, 6.02], curvature 0.67, R2 0.902 ss19: steps [175, 250, 360, 450], NRM [1.0, 0.834, 0.925, 0.367], TRM gained [0.19, 0.305, 0.338, 0.694], curvature -2.59, R2 0.956
3 vectors, random: (a) holds 98%, (b) 56%, median line MAD 21.1 deg 3 vectors, noisy line: (a) holds 78%, (b) 51%, median line MAD 16.7 deg
4 vectors, random: (a) holds 95%, (b) 56%, median line MAD 28.1 deg 4 vectors, noisy line: (a) holds 66%, (b) 51%, median line MAD 23.6 deg
5 vectors, random: (a) holds 93%, (b) 56%, median line MAD 32.1 deg 5 vectors, noisy line: (a) holds 61%, (b) 51%, median line MAD 27.0 deg
Result: reproduced, resting on the averaging assumption. For both specimens the field-off direction turns by more than the MAD of a line over the range, and a great circle fits better than a line (numbers above); R2 is 0.38 for both with averaged replicates. With first measurements ss41's R2 is 0.16 and ss22's 0.44, outside the tolerance. Both readings are weak tests. With no changing field and 3 vectors, as in ss41's range, reading (a) holds for 98% of random directions and 78% of points along one line with noise, and reading (b) for 56% and 51%; the median line MADs, 21.1 and 16.7 deg, are below ss41's 24.1. With 4 vectors, as in ss22's range, reading (a) holds for 95% and 66% and reading (b) for 56% and 51%, with median line MADs of 28.1 and 23.6 deg, below ss22's 35.4. The "curved" part of the statement is not in the criterion. The circle fit gives curvature -0.41 for ss41 and -1.12 for ss22, but -2.59 for ss19 and 0.67 for ss43, which the paper presents as linear, so with 3 or 4 points it does not separate these specimens; ss22's NRM does not fall monotonically over its range (1.00, 0.65, 0.82, 0.28).
glass-directions-not-random¶
The field-off directions of the failed glass specimens, as seen on orthogonal vector plots, are not random, implying the presence of an ambient field.
Criterion (a standard choice, the paper names no test): per specimen, Watson's test of randomness on the field-off directions of its thermal steps (p < 0.05). Reproduced when more than half of the 22 reject randomness, partially for a quarter to a half. Also given for the steps within each Table 2 range.
RANGES = {**{s: v[1:3] for s, v in T2.items()}, **{s: v[1:3] for s, v in NOFIT.items()}}
rows = []
for i, n in enumerate(TABLE2):
s = SHEET[n]
for how in ("mean", "first"):
st = A.step_table(glass[s]["data"], how)
off = st[st.kind == "off"]
f = A.fisher(off[["x", "y", "z"]].to_numpy())
p = A.watson_p(f["n"], f["R"], seed=3000 + i)
t0, t1 = RANGES[n]
o2 = off[(off.t >= t0) & (off.t <= t1)]
f2 = A.fisher(o2[["x", "y", "z"]].to_numpy())
p2 = A.watson_p(f2["n"], f2["R"], seed=3100 + i)
rows.append({"specimen": n, "how": how, "n": f["n"], "R": round(f["R"], 2), "k": round(f["k"], 1), "p": p,
"n_range": f2["n"], "R_range": round(f2["R"], 2), "p_range": p2, "synthetic sheet": n in SYNTH})
w = pd.DataFrame(rows)
print(w.to_string(index=False))
for how in ("mean", "first"):
x = w[w.how == how]
print(f"{how}: all thermal field-off steps reject randomness {int((x.p < 0.05).sum())}/22; "
f"within Table 2 ranges {int((x.p_range < 0.05).sum())}/22 "
f"(ranges of 2 steps: {int((x.n_range <= 2).sum())}); without synthetic sheets {int((x[~x['synthetic sheet']].p < 0.05).sum())}/18")
specimen how n R k p n_range R_range p_range synthetic sheet
ss19 mean 7 5.06 3.1 0.006060 4 3.36 0.018560 False
ss19 first 7 4.84 2.8 0.011175 4 3.31 0.023245 False
ss43 mean 5 3.41 2.5 0.061075 3 2.18 0.206510 False
ss43 first 5 3.17 2.2 0.102730 3 2.16 0.216660 False
ss22 mean 8 1.59 1.1 0.826715 4 1.41 0.718250 False
ss22 first 8 2.41 1.3 0.552940 4 1.14 0.832740 False
ss41 mean 9 5.85 2.5 0.006060 3 2.63 0.046715 False
ss41 first 9 5.81 2.5 0.006585 3 2.57 0.061240 False
ss6 mean 4 3.35 4.6 0.020215 3 2.41 0.114260 False
ss6 first 4 3.46 5.5 0.012075 3 2.49 0.085915 False
ss9 mean 6 5.00 5.0 0.001370 4 3.47 0.010770 False
ss9 first 6 5.04 5.2 0.001185 4 3.51 0.008915 False
ss10 mean 6 1.71 1.2 0.712300 3 2.65 0.041780 False
ss10 first 6 1.26 1.1 0.863415 3 2.46 0.097675 False
ss30 mean 5 2.55 1.6 0.282850 4 1.93 0.450130 False
ss30 first 5 2.06 1.4 0.490575 4 1.70 0.572660 False
ss34 mean 6 3.51 2.0 0.097385 5 2.75 0.214735 False
ss34 first 6 3.45 2.0 0.107790 5 2.72 0.225495 False
g2 mean 4 1.72 1.3 0.565680 4 1.72 0.563120 True
g2 first 4 1.81 1.4 0.516550 4 1.81 0.513760 True
ss8 mean 5 4.36 6.3 0.001935 5 4.36 0.001905 True
ss8 first 5 4.28 5.5 0.003105 5 4.28 0.003155 True
s11 mean 6 4.28 2.9 0.017905 6 4.28 0.017680 False
s11 first 6 3.95 2.4 0.039865 6 3.95 0.040145 False
ss18 mean 5 1.44 1.1 0.764110 5 1.44 0.763600 False
ss18 first 5 1.57 1.2 0.712615 5 1.57 0.713250 False
ss24 mean 7 2.95 1.5 0.302360 7 2.95 0.302160 False
ss24 first 7 3.10 1.5 0.255565 7 3.10 0.255355 False
ss27 mean 5 2.57 1.6 0.276760 5 2.57 0.276475 False
ss27 first 5 2.64 1.7 0.251920 5 2.64 0.251970 False
ss29 mean 6 2.62 1.5 0.341280 6 2.62 0.341320 False
ss29 first 6 5.04 5.2 0.001125 6 5.04 0.001085 False
ss36 mean 4 2.23 1.7 0.308675 4 2.23 0.308700 False
ss36 first 4 2.26 1.7 0.296725 4 2.26 0.296410 False
ss38 mean 4 3.89 27.9 0.000120 4 3.89 0.000100 False
ss38 first 4 3.97 102.3 0.000005 4 3.97 0.000005 False
ss39 mean 2 0.34 0.6 0.970765 2 0.34 0.970795 True
ss39 first 2 0.36 0.6 0.967150 2 0.36 0.967035 True
ss45 mean 4 3.67 9.1 0.002700 4 3.67 0.002785 False
ss45 first 4 3.67 9.2 0.002635 4 3.67 0.002680 False
ss46 mean 4 1.46 1.2 0.694145 3 2.04 0.272955 False
ss46 first 4 1.43 1.2 0.707490 3 1.97 0.307730 False
ss48 mean 2 1.77 4.4 0.215700 2 1.77 0.215205 True
ss48 first 2 1.76 4.2 0.222075 2 1.76 0.221385 True
mean: all thermal field-off steps reject randomness 8/22; within Table 2 ranges 8/22 (ranges of 2 steps: 2); without synthetic sheets 7/18
first: all thermal field-off steps reject randomness 9/22; within Table 2 ranges 7/22 (ranges of 2 steps: 2); without synthetic sheets 8/18
Result: partially reproduced. Watson's test rejects randomness of the thermal field-off directions in 8 of 22 specimens with averaged replicates (36%, within the quarter to half band), 9 of 22 with first measurements, and 8 of 22 within the Table 2 ranges. The remaining specimens give directions that the test cannot tell from random at these numbers of steps (2 to 9).
glass-af-no-change¶
For nearly all glass samples given AF pre-treatment, AF cleaning did not significantly change the NRM direction or magnitude.
Criterion: per specimen with AF steps, the angle between the NRM and the direction after the last AF step, and the ratio of the moments. Readings (a) angle under 15 deg and moment change under 15%; (b) angle under the Table 2 MAD where printed, and moment change under 15%. Reproduced when, under one reading, all but at most two show no significant change.
MADT2 = {s: v[3] for s, v in T2.items()}
rows = []
for s in has_af:
for how in ("mean", "first"):
st = A.step_table(glass[s]["data"], how)
nrm = st[st.kind == "nrm"].iloc[0]
last = st[st.kind == "af"].iloc[-1]
ang = A.angle(A.vec(nrm), A.vec(last))
ratio = last.moment / nrm.moment
lim = MADT2.get(s, 15.0)
rows.append({"specimen": s, "how": how, "last AF": last.step, "angle": round(ang, 1), "moment ratio": round(ratio, 3),
"ok_a": ang < 15 and abs(ratio - 1) < 0.15, "MAD limit (b)": lim,
"ok_b": ang < lim and abs(ratio - 1) < 0.15, "synthetic sheet": s in SYNTH})
af = pd.DataFrame(rows)
print(af.to_string(index=False))
for how in ("mean", "first"):
x = af[af.how == how]
print(f"{how}: no significant change under (a) {int(x.ok_a.sum())}/{len(x)}, under (b) {int(x.ok_b.sum())}/{len(x)}")
print("ss38 NRM rows ('0' and '0.0'):")
print(glass["ss38"]["data"].query("step in ['0', '0.0']")[["row", "step", "moment", "dec", "inc"]].to_string(index=False))
specimen how last AF angle moment ratio ok_a MAD limit (b) ok_b synthetic sheet
ss19 mean 20mT 2.0 0.997 True 30.1 True False
ss19 first 20mT 2.3 0.999 True 30.1 True False
ss22 mean 20mT 2.8 1.011 True 35.4 True False
ss22 first 20mT 3.6 1.021 True 35.4 True False
ss27 mean 10mT 1.5 1.010 True 15.0 True False
ss27 first 10mT 2.4 0.994 True 15.0 True False
ss34 mean 10mT 2.2 1.006 True 12.1 True False
ss34 first 10mT 2.2 1.025 True 12.1 True False
ss38 mean 20mT 2.4 0.488 False 15.0 False False
ss38 first 20mT 6.6 0.479 False 15.0 False False
ss41 mean 10mT 0.3 1.001 True 24.0 True False
ss41 first 10mT 0.6 0.997 True 24.0 True False
ss43 mean 10mT 0.2 1.002 True 10.0 True False
ss43 first 10mT 0.1 1.004 True 10.0 True False
ss45 mean 10mT 0.1 0.999 True 15.0 True False
ss45 first 10mT 0.4 0.997 True 15.0 True False
ss46 mean 10mT 1.0 0.996 True 15.0 True False
ss46 first 10mT 1.2 0.995 True 15.0 True False
ss48 mean 10mT 17.5 0.671 False 15.0 False True
ss48 first 10mT 22.1 0.628 False 15.0 False True
mean: no significant change under (a) 8/10, under (b) 8/10
first: no significant change under (a) 8/10, under (b) 8/10
ss38 NRM rows ('0' and '0.0'):
row step moment dec inc
2 0 1.034000e-07 182.2 -50.4
3 0 1.021000e-07 183.3 -51.0
4 0 1.019000e-07 183.1 -51.3
5 0.0 1.032000e-07 182.2 -52.3
6 0.0 1.037000e-07 182.0 -52.8
Result: reproduced. AF cleaning turns the NRM by at most 2.8 deg and changes its moment by at most 1.1% in 8 of the 10 specimens; ss38 keeps its direction but loses half its moment (ratio 0.49 at 20 mT) and ss48, a synthetic sheet, turns by 17.5 deg and loses a third. Two exceptions are allowed.
table1-max-intensity¶
For each of the twelve single crystals, the maximum field that cannot be excluded, given the magnetometer threshold of 8e-11 emu and the moment acquired at 590 C in a 20 uT field, is as in Table 1, between 0.20 and 0.51 uT.
Criterion: max intensity = 8e-11 emu / M590 x 20 uT. M590 under readings (a) |590.1 - 590.0|, (b) |590.1|, with (i) the first 590.0 replicate or (ii) their mean. Reproduced when, under one reading, all 24 values are within tolerance; partially when at least 18 are.
T1 = {"12021g2": (7.011e-9, 0.23), "12021g3": (6.292e-9, 0.26), "12040g1": (8.002e-9, 0.20), "12040g2": (6.021e-9, 0.27),
"12040g3": (6.292e-9, 0.26), "12053g1": (5.283e-9, 0.31), "12053g2": (7.982e-9, 0.20), "14053g1": (3.203e-9, 0.51),
"14053g2": (5.294e-9, 0.31), "14053g3": (4.523e-9, 0.36), "71055g1": (3.883e-9, 0.42), "71055g3": (4.202e-9, 0.39)}
def grp(name, step):
d = crystals[name]["data"]
d = d[d.step == step]
return A.to_xyz(d.moment, d.dec, d.inc)
def m590(name, reading):
on = grp(name, "590.1")[0]
off = grp(name, "590")
off = off[0] if reading.endswith("_i") else off.mean(axis=0)
return float(np.linalg.norm(on - off)) if reading.startswith("a") else float(np.linalg.norm(on))
rows = []
for name, (m_t1, i_t1) in T1.items():
r = {"crystal": name, "M590 T1": m_t1, "max T1": i_t1}
for rd in ("a_i", "a_ii", "b_i", "b_ii"):
m = m590(name, rd)
r[f"M_{rd}"] = m
r[f"max_{rd}"] = round(A.max_intensity(m), 3)
r[f"ok_{rd}"] = int(abs(m - m_t1) <= 0.1 * m_t1) + int(abs(A.max_intensity(m) - i_t1) <= 0.1 * i_t1 + 1e-9)
r["T1 formula on own M590"] = round(A.max_intensity(m_t1), 3)
rows.append(r)
t1 = pd.DataFrame(rows)
print(t1.to_string(index=False))
for rd in ("a_i", "a_ii", "b_i", "b_ii"):
print(f"reading {rd}: {int(t1[f'ok_{rd}'].sum())}/24 values within tolerance")
print("Table 1 max intensity equal to the formula on its own M590, to two places:",
int((t1["T1 formula on own M590"].round(2) == t1["max T1"]).sum()), "/ 12")
import math
print("Table 1 max intensity equal to the formula on its own M590, rounded up to two places:",
int(sum(math.ceil(round(A.max_intensity(m) * 100, 6)) / 100 == i for m, i in zip(t1["M590 T1"], t1["max T1"]))), "/ 12")
print("Formula on Table 1's own M590, unrounded:", {n: round(A.max_intensity(m), 4) for n, m in zip(t1.crystal, t1["M590 T1"])})
print("Exact equality of M590 (b) with Table 1 (to 4 significant digits):",
int(sum(float(f"{m:.4g}") == mt for m, mt in zip(t1.M_b_i, t1["M590 T1"]))), "/ 12")
print("12021g3 and 12040g3 590.1 rows:")
print(pd.concat([crystals[n]["data"].query("step == '590.1'") for n in ("12021g3", "12040g3")]).to_string(index=False))
crystal M590 T1 max T1 M_a_i max_a_i ok_a_i M_a_ii max_a_ii ok_a_ii M_b_i max_b_i ok_b_i M_b_ii max_b_ii ok_b_ii T1 formula on own M590
12021g2 7.011000e-09 0.23 6.625206e-09 0.242 2 6.827007e-09 0.234 2 7.011000e-09 0.228 2 7.011000e-09 0.228 2 0.228
12021g3 6.292000e-09 0.26 5.955885e-09 0.269 2 6.396044e-09 0.250 2 6.292000e-09 0.254 2 6.292000e-09 0.254 2 0.254
12040g1 8.002000e-09 0.20 7.842468e-09 0.204 2 7.877593e-09 0.203 2 8.002000e-09 0.200 2 8.002000e-09 0.200 2 0.200
12040g2 6.021000e-09 0.27 5.668901e-09 0.282 2 5.839592e-09 0.274 2 6.021000e-09 0.266 2 6.021000e-09 0.266 2 0.266
12040g3 6.292000e-09 0.26 6.674872e-09 0.240 2 6.278794e-09 0.255 2 6.292000e-09 0.254 2 6.292000e-09 0.254 2 0.254
12053g1 5.283000e-09 0.31 5.019062e-09 0.319 2 5.208645e-09 0.307 2 5.283000e-09 0.303 2 5.283000e-09 0.303 2 0.303
12053g2 7.982000e-09 0.20 7.757840e-09 0.206 2 7.886535e-09 0.203 2 7.982000e-09 0.200 2 7.982000e-09 0.200 2 0.200
14053g1 3.203000e-09 0.51 3.183580e-09 0.503 2 3.194190e-09 0.501 2 3.203000e-09 0.500 2 3.203000e-09 0.500 2 0.500
14053g2 5.294000e-09 0.31 5.012061e-09 0.319 2 5.161885e-09 0.310 2 5.294000e-09 0.302 2 5.294000e-09 0.302 2 0.302
14053g3 4.523000e-09 0.36 4.405499e-09 0.363 2 4.507697e-09 0.355 2 4.523000e-09 0.354 2 4.523000e-09 0.354 2 0.354
71055g1 3.883000e-09 0.42 3.756393e-09 0.426 2 3.801994e-09 0.421 2 3.883000e-09 0.412 2 3.883000e-09 0.412 2 0.412
71055g3 4.202000e-09 0.39 3.736186e-09 0.428 1 4.026785e-09 0.397 2 4.202000e-09 0.381 2 4.202000e-09 0.381 2 0.381
reading a_i: 23/24 values within tolerance
reading a_ii: 24/24 values within tolerance
reading b_i: 24/24 values within tolerance
reading b_ii: 24/24 values within tolerance
Table 1 max intensity equal to the formula on its own M590, to two places: 4 / 12
Table 1 max intensity equal to the formula on its own M590, rounded up to two places: 10 / 12
Formula on Table 1's own M590, unrounded: {'12021g2': 0.2282, '12021g3': 0.2543, '12040g1': 0.2, '12040g2': 0.2657, '12040g3': 0.2543, '12053g1': 0.3029, '12053g2': 0.2005, '14053g1': 0.4995, '14053g2': 0.3022, '14053g3': 0.3537, '71055g1': 0.4121, '71055g3': 0.3808}
Exact equality of M590 (b) with Table 1 (to 4 significant digits): 12 / 12
12021g3 and 12040g3 590.1 rows:
row sample step moment dec inc
7 12021g3 590.1 6.292000e-09 178.3 -77.0
7 12040g3 590.1 6.292000e-09 351.3 -85.3
Result: reproduced. Under reading (b), M590 the magnitude of the single 590.1 measurement, all 24 values are within tolerance: the 12 M590 values equal Table 1 to its four digits and the max intensities run from 0.200 to 0.500 uT. Reading (a) with averaged 590.0 replicates also puts all 24 within tolerance (23 with the first replicate). Table 1's max intensities equal the formula on its own M590 column, rounded to two places, for only 4 of 12 crystals; rounded up, for 10 of 12, the exceptions being 14053g1 (0.4995 against 0.51) and 12053g2 (0.2005 against 0.20). The equal M590 of 12021g3 and 12040g3 is in the data: both 590.1 rows read 6.292e-09 emu, with different directions.
baseline-range¶
The maximum fields that the single crystal data cannot exclude define a baseline of about 0.2 to 0.5 uT below which fields cannot be evaluated with these data.
Criterion: the lowest and highest of the twelve values of table1-max-intensity, under the reading that matches there, lie in 0.18 to 0.22 and 0.45 to 0.55 uT.
vals = t1.max_b_i
print(f"Reading (b): lowest {vals.min():.3f} uT ({t1.crystal[vals.idxmin()]}), highest {vals.max():.3f} uT ({t1.crystal[vals.idxmax()]})")
print(f"Within 0.18-0.22: {0.18 <= vals.min() <= 0.22}; within 0.45-0.55: {0.45 <= vals.max() <= 0.55}")
for rd in ("a_i", "a_ii"):
print(f"Reading {rd} for comparison: {t1[f'max_{rd}'].min():.3f} to {t1[f'max_{rd}'].max():.3f} uT")
Reading (b): lowest 0.200 uT (12040g1), highest 0.500 uT (14053g1) Within 0.18-0.22: True; within 0.45-0.55: True Reading a_i for comparison: 0.204 to 0.503 uT Reading a_ii for comparison: 0.203 to 0.501 uT
Result: reproduced. The twelve values run from 0.200 to 0.500 uT.
example-14053g1-threshold¶
For crystal 14053g1, the threshold field detection limit computed from the magnetometer sensitivity (8e-14 A m2) and a TRM at 590 C of 8e-12 A m2 in a 20 uT field is 0.2 uT.
Criterion: readings (a) the example as written, for 14053g1; (b) a misprint of the crystal name, the numbers belonging to 12040g1. TRM as in table1-max-intensity (reading b matched there). Reproduced as written when (a) matches; when only (b) matches, reported as matching the misprint reading.
for name in ("14053g1", "12040g1"):
m_emu = m590(name, "b_i")
trm_am2 = m_emu * 1e-3 # 1 emu = 1e-3 A m2
thr = 8e-14 / trm_am2 * 20
print(f"{name}: TRM at 590 C {trm_am2:.3e} A m2 (expected 8.0e-12 +- 8e-13: {abs(trm_am2 - 8e-12) <= 8e-13}); "
f"threshold field {thr:.3f} uT (expected 0.2 +- 0.02: {abs(thr - 0.2) <= 0.02})")
14053g1: TRM at 590 C 3.203e-12 A m2 (expected 8.0e-12 +- 8e-13: False); threshold field 0.500 uT (expected 0.2 +- 0.02: False) 12040g1: TRM at 590 C 8.002e-12 A m2 (expected 8.0e-12 +- 8e-13: True); threshold field 0.200 uT (expected 0.2 +- 0.02: True)
Result: partially reproduced. As written, for 14053g1, the example does not match: the data give a TRM of 3.2e-12 A m2 and 0.50 uT. Under the misprint reading (b) it matches: 12040g1 gives 8.0e-12 A m2 and 0.200 uT, the example's numbers. The criterion calls the claim reproduced as written only when reading (a) matches, so a match under (b) alone is recorded as partial: the calculation holds, for a crystal other than the one named.
fig3-nrm-consistent, fig3-590-inconsistent, fig3-alteration-check-inconsistent¶
For most crystals the replicate NRM measurements have a consistent direction, as shown on stereonets for 12040g1 and 14053g2. After demagnetization at 590 C in zero field, the replicate directions are inconsistent, consistent with a null magnetization state. The replicate directions after the zero-field alteration check at 590 C are also inconsistent.
Criteria: per crystal, Fisher mean, k and alpha95 of the replicate directions of each step. NRM consistent under (a) k of the NRM replicates above k of the 590.0 replicates, or (b) Watson's test rejects randomness of the NRM replicates (p < 0.05). 590.0 (and 590.3) inconsistent under (a) k below the NRM's k, or (b) Watson's test does not reject randomness. Each reproduced when, under one reading, more than half of the twelve crystals qualify, including 12040g1 and 14053g2. A crystal with one NRM measurement has no k and no test, and counts as not qualifying.
rows = []
for i, name in enumerate(T1):
r = {"crystal": name}
for step, lab in (("0", "nrm"), ("590", "590.0"), ("590.3", "590.3"), ("400", "400.0"), ("400.3", "400.3")):
X = grp(name, step)
if len(X) == 0:
continue
f = A.fisher(X)
r[f"{lab} n"] = f["n"]
r[f"{lab} k"] = round(f["k"], 1) if f["k"] == f["k"] else np.nan
r[f"{lab} a95"] = round(f["a95"], 0) if f["a95"] == f["a95"] else np.nan
r[f"{lab} p"] = A.watson_p(f["n"], f["R"], seed=4000 + 10 * i + len(r))
r[f"{lab} mean moment"] = float(np.linalg.norm(X, axis=1).mean())
r[f"{lab} dec,inc"] = (round(f["dec"], 1), round(f["inc"], 1))
rows.append(r)
cr = pd.DataFrame(rows).set_index("crystal")
cols = [c for c in cr if c.split()[0] in ("nrm", "590.0", "590.3")]
print(cr[cols].to_string())
knrm = cr["nrm k"]
nrm_a = knrm > cr["590.0 k"]
nrm_b = cr["nrm p"] < 0.05
c590_a = cr["590.0 k"] < knrm
c590_b = cr["590.0 p"] >= 0.05
c5903_a = cr["590.3 k"] < knrm
c5903_b = cr["590.3 p"] >= 0.05
summ = pd.DataFrame({"NRM consistent (a)": nrm_a, "NRM consistent (b)": nrm_b, "590.0 inconsistent (a)": c590_a,
"590.0 inconsistent (b)": c590_b, "590.3 inconsistent (a)": c5903_a, "590.3 inconsistent (b)": c5903_b})
print(summ.to_string())
print("Counts of 12:", summ.sum().astype(int).to_dict())
print("12040g1 and 14053g2:", summ.loc[["12040g1", "14053g2"]].to_dict("index"))
print(f"590.0 mean moments: {cr['590.0 mean moment'].min():.3g} to {cr['590.0 mean moment'].max():.3g} emu; "
f"above the 8e-11 emu threshold: {int((cr['590.0 mean moment'] > A.THRESHOLD_EMU).sum())}/12; "
f"single 590.0 measurements above it: {sum(int((np.linalg.norm(grp(n, '590'), axis=1) > A.THRESHOLD_EMU).sum()) for n in T1)}"
f"/{sum(len(grp(n, '590')) for n in T1)}")
print(f"590.3 mean moments: {cr['590.3 mean moment'].min():.3g} to {cr['590.3 mean moment'].max():.3g} emu; "
f"NRM mean moments: {cr['nrm mean moment'].min():.3g} to {cr['nrm mean moment'].max():.3g} emu")
nrm n nrm k nrm a95 nrm p nrm mean moment nrm dec,inc 590.0 n 590.0 k 590.0 a95 590.0 p 590.0 mean moment 590.0 dec,inc 590.3 n 590.3 k 590.3 a95 590.3 p 590.3 mean moment 590.3 dec,inc
crystal
12021g2 2 600.3 10.0 0.001690 8.220000e-10 (165.8, -45.3) 3 1.2 180.0 0.632000 4.383333e-10 (194.9, 54.7) 3 1.2 180.0 0.635610 4.952667e-10 (88.0, 63.2)
12021g3 2 213.6 17.0 0.004360 7.408500e-10 (76.0, -25.0) 3 1.1 180.0 0.764630 5.044000e-10 (250.0, 38.0) 3 0.8 180.0 0.969795 4.378333e-10 (183.4, 62.7)
12040g1 2 981.4 8.0 0.001025 7.002500e-10 (131.9, 55.4) 3 1.8 180.0 0.344305 2.477667e-10 (276.2, -39.9) 3 0.8 180.0 0.985240 3.258000e-10 (141.5, -47.8)
12040g2 3 22.8 26.0 0.002975 5.202333e-10 (153.6, 23.3) 3 1.4 180.0 0.550200 4.442667e-10 (117.8, 59.6) 3 1.2 180.0 0.659550 4.328333e-10 (138.9, 60.9)
12040g3 2 1.3 180.0 0.628940 3.977000e-10 (48.5, 22.3) 3 0.8 180.0 0.983025 4.171667e-10 (117.1, -7.8) 3 1.0 180.0 0.788230 4.274667e-10 (308.1, 33.0)
12053g1 2 9.8 91.0 0.100090 5.177500e-10 (154.1, -47.6) 3 0.9 180.0 0.909900 3.230000e-10 (228.3, -53.3) 3 1.2 180.0 0.693440 3.794333e-10 (252.6, -70.6)
12053g2 3 37.8 20.0 0.001265 7.231333e-10 (294.9, 29.5) 3 0.8 180.0 0.987440 5.684667e-10 (187.8, 60.6) 3 1.6 180.0 0.442460 5.119667e-10 (142.3, 29.5)
14053g1 2 3211.3 4.0 0.000310 1.235500e-09 (113.8, -44.5) 3 1.0 180.0 0.784930 1.164667e-10 (167.8, -27.6) 3 0.9 180.0 0.951775 2.142667e-10 (296.1, -30.0)
14053g2 2 205.7 18.0 0.004630 6.084000e-10 (155.0, -45.4) 3 1.4 180.0 0.528620 2.982667e-10 (216.5, -50.8) 3 1.6 180.0 0.436235 3.011333e-10 (122.8, 29.8)
14053g3 2 38.0 42.0 0.025995 5.492500e-10 (40.9, 58.3) 4 1.2 180.0 0.672355 2.758000e-10 (111.8, -4.2) 3 0.8 180.0 0.981145 3.745333e-10 (122.3, 29.8)
71055g1 1 NaN NaN NaN 8.893000e-10 (216.3, 49.2) 3 1.3 180.0 0.604220 3.468667e-10 (155.7, 22.1) 3 1.5 180.0 0.470745 4.021667e-10 (163.2, 53.5)
71055g3 1 NaN NaN NaN 1.092000e-09 (102.2, -34.8) 3 1.6 180.0 0.424215 5.325667e-10 (166.2, -25.3) 3 0.9 180.0 0.911965 5.049667e-10 (208.5, -1.6)
NRM consistent (a) NRM consistent (b) 590.0 inconsistent (a) 590.0 inconsistent (b) 590.3 inconsistent (a) 590.3 inconsistent (b)
crystal
12021g2 True True True True True True
12021g3 True True True True True True
12040g1 True True True True True True
12040g2 True True True True True True
12040g3 True False True True True True
12053g1 True False True True True True
12053g2 True True True True True True
14053g1 True True True True True True
14053g2 True True True True True True
14053g3 True True True True True True
71055g1 False False False True False True
71055g3 False False False True False True
Counts of 12: {'NRM consistent (a)': 10, 'NRM consistent (b)': 8, '590.0 inconsistent (a)': 10, '590.0 inconsistent (b)': 12, '590.3 inconsistent (a)': 10, '590.3 inconsistent (b)': 12}
12040g1 and 14053g2: {'12040g1': {'NRM consistent (a)': True, 'NRM consistent (b)': True, '590.0 inconsistent (a)': True, '590.0 inconsistent (b)': True, '590.3 inconsistent (a)': True, '590.3 inconsistent (b)': True}, '14053g2': {'NRM consistent (a)': True, 'NRM consistent (b)': True, '590.0 inconsistent (a)': True, '590.0 inconsistent (b)': True, '590.3 inconsistent (a)': True, '590.3 inconsistent (b)': True}}
590.0 mean moments: 1.16e-10 to 5.68e-10 emu; above the 8e-11 emu threshold: 12/12; single 590.0 measurements above it: 37/37
590.3 mean moments: 2.14e-10 to 5.12e-10 emu; NRM mean moments: 3.98e-10 to 1.24e-09 emu
Result: fig3-nrm-consistent reproduced. NRM replicates are more tightly grouped than the 590.0 replicates in 10 of 12 crystals (reading a; 8 of 12 under reading b), including 12040g1 and 14053g2; 71055g1 and 71055g3 have one NRM measurement each.
Result: fig3-590-inconsistent reproduced. Watson's test does not reject randomness of the 590.0 replicates in any of the 12 crystals (reading b), and their k is below the NRM's in the 10 crystals with two or more NRM measurements (reading a). The 590.0 moments, 1.2e-10 to 5.7e-10 emu on average, are all above the 8e-11 emu threshold.
Result: fig3-alteration-check-inconsistent reproduced. The same holds for the 590.3 replicates: 12 of 12 under reading (b), 10 of 12 under reading (a).
crystal-14053g1-400c-stable¶
Crystal 14053g1, measured after heating to 400 C during initial tests, held a stable remanence, which was removed by treatment at 590 C.
Criterion: applies if the deposit holds 400 C measurements for 14053g1. (a) the 400 C replicates are consistent under fig3-nrm-consistent's readings and the 590.0 replicates inconsistent; or (b) the 400 C direction lies within the alpha95 of the Fisher mean of the NRM replicates. Reproduced when the 590.0 replicates are inconsistent and (a) or (b) holds.
r = cr.loc["14053g1"]
print("14053g1 rows at 400 C:")
print(crystals["14053g1"]["data"].query("step.str.startswith('400')").to_string(index=False))
print(f"NRM: n {r['nrm n']}, mean {r['nrm dec,inc']}, k {r['nrm k']}, a95 {r['nrm a95']}, mean moment {r['nrm mean moment']:.3g} emu")
print(f"400.0: n {r['400.0 n']}, mean {r['400.0 dec,inc']}, k {r['400.0 k']}, a95 {r['400.0 a95']}, Watson p {r['400.0 p']:.4f}, mean moment {r['400.0 mean moment']:.3g}")
print(f"400.3: n {r['400.3 n']}, mean {r['400.3 dec,inc']}, k {r['400.3 k']}")
print(f"590.0: n {r['590.0 n']}, k {r['590.0 k']}, Watson p {r['590.0 p']:.3f}, mean moment {r['590.0 mean moment']:.3g}")
a_k = r["400.0 k"] > r["590.0 k"]
a_w = r["400.0 p"] < 0.05
inc590 = (r["590.0 k"] < r["nrm k"]) or (r["590.0 p"] >= 0.05)
m400 = grp("14053g1", "400").sum(axis=0)
mnrm = grp("14053g1", "0")
mnrm = (mnrm / np.linalg.norm(mnrm, axis=1, keepdims=True)).sum(axis=0)
ang = A.angle(m400, mnrm)
print(f"(a) 400 C consistent: k above 590.0's {a_k}, Watson rejects {a_w}; 590.0 inconsistent {inc590}")
print(f"(b) angle between the 400 C mean and the NRM mean {ang:.1f} deg, NRM alpha95 {r['nrm a95']} deg: within {ang <= r['nrm a95']}")
for i, x in enumerate(grp("14053g1", "400")):
print(f" 400 C replicate {i + 1}: {A.angle(x, mnrm):.1f} deg from the NRM mean")
14053g1 rows at 400 C:
row sample step moment dec inc
4 14053g1 400 8.293000e-10 134.9 -54.8
5 14053g1 400 8.188000e-10 112.3 -44.5
6 14053g1 400 8.094000e-10 128.7 -50.2
7 14053g1 400.1 1.612000e-09 136.9 -59.1
8 14053g1 400.3 7.923000e-10 129.8 -50.2
9 14053g1 400.3 7.894000e-10 131.2 -48.7
NRM: n 2, mean (113.8, -44.5), k 3211.3, a95 4.0, mean moment 1.24e-09 emu
400.0: n 3.0, mean (124.5, -50.2), k 78.0, a95 14.0, Watson p 0.0003, mean moment 8.19e-10
400.3: n 2.0, mean (130.5, -49.5), k 4266.1
590.0: n 3, k 1.0, Watson p 0.785, mean moment 1.16e-10
(a) 400 C consistent: k above 590.0's True, Watson rejects True; 590.0 inconsistent True
(b) angle between the 400 C mean and the NRM mean 9.3 deg, NRM alpha95 4.0 deg: within False
400 C replicate 1: 17.0 deg from the NRM mean
400 C replicate 2: 1.1 deg from the NRM mean
400 C replicate 3: 11.6 deg from the NRM mean
Result: reproduced. The deposit holds three 400 C measurements of 14053g1 (with a 400.1 and two 400.3 rows). They are grouped (k above that of the 590.0 replicates, and Watson's test rejects randomness), while the 590.0 replicates are not, so reading (a) holds. Under reading (b) the 400 C mean lies outside the NRM's alpha95 (numbers above).
no-alteration-recording-efficiency¶
Alteration during the laser heating of the single crystals was not significant; apparent recording efficiencies far above 100% were not observed and the alteration checks did not fail.
Criterion: (a) the 590.3 moment does not exceed the 590.0 moment by more than 8e-11 emu, and the 590.3 directions are inconsistent; (b) the 590.4 TRM (|590.4 - 590.3|) is no more than about twice the 590.1 TRM (|590.1 - 590.0|); a ratio of the 590.4 TRM to twice the 590.1 TRM above 2 fails (a choice not from the paper). Reproduced when every crystal passes under one reading.
rows = []
for name in T1:
on, off, alt, on40 = grp(name, "590.1")[0], grp(name, "590"), grp(name, "590.3"), grp(name, "590.4")[0]
for rep in ("i", "ii"):
o = off[0] if rep == "i" else off.mean(axis=0)
al = alt[0] if rep == "i" else alt.mean(axis=0)
m_off = np.linalg.norm(off[0]) if rep == "i" else np.linalg.norm(off, axis=1).mean()
m_alt = np.linalg.norm(alt[0]) if rep == "i" else np.linalg.norm(alt, axis=1).mean()
ratio = np.linalg.norm(on40 - al) / (2 * np.linalg.norm(on - o))
rows.append({"crystal": name, "replicates": rep, "590.3 - 590.0 moment (1e-11 emu)": round((m_alt - m_off) * 1e11, 1),
"a_moment_ok": m_alt - m_off <= A.THRESHOLD_EMU, "a_dirs_inconsistent": bool(summ.loc[name, "590.3 inconsistent (b)"]),
"TRM40/(2 TRM20)": round(ratio, 3), "b_ok": ratio <= 2})
alt = pd.DataFrame(rows)
alt["a_ok"] = alt.a_moment_ok & alt.a_dirs_inconsistent
print(alt.to_string(index=False))
for rep in ("i", "ii"):
x = alt[alt.replicates == rep]
print(f"replicates {rep}: pass (a) {int(x.a_ok.sum())}/12, pass (b) {int(x.b_ok.sum())}/12; "
f"TRM40/(2 TRM20) {x['TRM40/(2 TRM20)'].min():.2f} to {x['TRM40/(2 TRM20)'].max():.2f}")
# The ratio without the factor of two, and the cutoff of 2 on it, which the claims file does not use
x = alt[alt.replicates == "i"].assign(**{"TRM40/TRM20": lambda d: (2 * d["TRM40/(2 TRM20)"]).round(2)})
print(x[["crystal", "TRM40/TRM20"]].to_string(index=False))
print(f"TRM40/TRM20 above 2: {list(x.crystal[x['TRM40/TRM20'] > 2])}")
# Scatter of the 590.0 replicate moments, against the 8e-11 emu threshold
mom = {n: np.linalg.norm(grp(n, "590"), axis=1) * 1e11 for n in T1}
print(pd.DataFrame([{"crystal": n, "590.0 moments (1e-11 emu)": v.round(1).tolist(), "sd": round(v.std(ddof=1), 1)}
for n, v in mom.items()]).to_string(index=False))
crystal replicates 590.3 - 590.0 moment (1e-11 emu) a_moment_ok a_dirs_inconsistent TRM40/(2 TRM20) b_ok a_ok 12021g2 i 5.9 True True 0.952 True True 12021g2 ii 5.7 True True 0.932 True True 12021g3 i -7.9 True True 0.910 True True 12021g3 ii -6.7 True True 0.865 True True 12040g1 i 7.3 True True 0.781 True True 12040g1 ii 7.8 True True 0.788 True True 12040g2 i -2.7 True True 0.947 True True 12040g2 ii -1.1 True True 0.927 True True 12040g3 i 2.7 True True 0.770 True True 12040g3 ii 1.0 True True 0.796 True True 12053g1 i 5.4 True True 1.174 True True 12053g1 ii 5.6 True True 1.141 True True 12053g2 i -7.9 True True 1.028 True True 12053g2 ii -5.6 True True 1.018 True True 14053g1 i 12.1 False True 0.935 True False 14053g1 ii 9.8 False True 0.950 True False 14053g2 i -0.9 True True 0.859 True True 14053g2 ii 0.3 True True 0.862 True True 14053g3 i 6.8 True True 1.038 True True 14053g3 ii 9.9 False True 1.030 True False 71055g1 i 7.1 True True 0.800 True True 71055g1 ii 5.5 True True 0.803 True True 71055g3 i -2.7 True True 0.995 True True 71055g3 ii -2.8 True True 0.969 True True replicates i: pass (a) 11/12, pass (b) 12/12; TRM40/(2 TRM20) 0.77 to 1.17 replicates ii: pass (a) 10/12, pass (b) 12/12; TRM40/(2 TRM20) 0.79 to 1.14 crystal TRM40/TRM20 12021g2 1.90 12021g3 1.82 12040g1 1.56 12040g2 1.89 12040g3 1.54 12053g1 2.35 12053g2 2.06 14053g1 1.87 14053g2 1.72 14053g3 2.08 71055g1 1.60 71055g3 1.99 TRM40/TRM20 above 2: ['12053g1', '12053g2', '14053g3'] crystal 590.0 moments (1e-11 emu) sd 12021g2 [43.9, 41.3, 46.3] 2.5 12021g3 [50.3, 49.2, 51.8] 1.3 12040g1 [24.9, 31.0, 18.4] 6.3 12040g2 [42.9, 43.1, 47.3] 2.5 12040g3 [40.3, 43.3, 41.6] 1.5 12053g1 [32.8, 34.8, 29.2] 2.8 12053g2 [58.2, 52.0, 60.3] 4.3 14053g1 [10.2, 13.3, 11.4] 1.6 14053g2 [31.2, 29.4, 28.9] 1.2 14053g3 [30.0, 19.8, 29.4, 31.1] 5.2 71055g1 [34.9, 35.3, 33.8] 0.8 71055g3 [52.9, 53.9, 52.9] 0.6
Result: reproduced. Under reading (b) every crystal passes: the TRM acquired in 40 uT is 0.77 to 1.17 times twice the TRM acquired in 20 uT, close to linear recording and far from the cutoff of 2. Under reading (a) 11 of 12 pass with first replicates and 10 of 12 with means; 14053g1's 590.3 moment exceeds its 590.0 moment by 9.8e-11 to 1.21e-10 emu, above the 8e-11 emu allowance. The cutoff is loose: without the factor of two, TRM40/TRM20 is 1.54 to 2.35, and three crystals (12053g1 2.35, 14053g3 2.08, 12053g2 2.06) are above 2, so a cutoff of 2 on that ratio would fail them and no reading would pass every crystal. The paper's own words are "apparent recording efficiencies >> 100%", and the largest here is 117% of linear recording, so the substance holds; the cutoff is asked of a reader (questions/figshare-24639285.yaml). The 590.0 replicate moments of a crystal scatter by 0.6e-11 to 6.3e-11 emu (sd), below the 8e-11 emu allowance of reading (a); 14053g1's 590.3 excess is about 6 to 8 times its 590.0 sd of 1.6e-11 emu, so its failure under reading (a) is not replicate noise. The size of the 590.0 moments, 1.0e-10 to 6.0e-10 emu with random directions, is 1.3 to 7.5 times the threshold.
interpretation-impact-no-dynamo¶
Together with ref. 1, the data are consistent with impact plasma magnetization of the 64455 glass and with the absence of a long-lived lunar core dynamo recorded in the 3.9 to 3.2 Ga crystals.
The claims file says this interpretation is not tested; it rests on ref. 1 and on models of impact plasmas. It is recorded as blocked. The outcomes of the four claims it names are: fig2-changing-field reproduced, glass-directions-not-random partially reproduced, fig3-590-inconsistent reproduced, table1-max-intensity reproduced.
Summary¶
| Claim | Outcome |
|---|---|
| glass-specimen-set (link) | reproduced |
| glass-success-rate | reproduced |
| table2-fit-statistics | partially reproduced (assumes averaged replicates) |
| table2-no-fit | reproduced |
| fig1-linear-components | reproduced (assumes averaged replicates) |
| fig2-changing-field | reproduced (assumes averaged replicates) |
| glass-directions-not-random | partially reproduced |
| glass-af-no-change | reproduced |
| table1-max-intensity | reproduced |
| baseline-range | reproduced |
| example-14053g1-threshold | partially reproduced (matches only as a misprint for 12040g1) |
| fig3-nrm-consistent | reproduced |
| fig3-590-inconsistent | reproduced |
| fig3-alteration-check-inconsistent | reproduced |
| crystal-14053g1-400c-stable | reproduced |
| no-alteration-recording-efficiency | reproduced |
| interpretation-impact-no-dynamo | blocked (interpretation, not tested) |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: the paper used PmagPy and gives no formulas; MAD, R2 and pTRM checks are computed here by the claims file's readings, and repeated measurements are averaged, an assumption that reproduces most of Table 2 and that the deposit's Supplementary Figure 1 supports. No test of randomness is named by the paper; Watson's test is the claims file's choice.
What was found about the dataset:
- The "replicates" of four glass sheets (g2, ss8, ss39, ss48) are one measurement multiplied by factors near 1 (0.95 to 1.03), applied to moment, declination and inclination; they are not repeated measurements.
- One inclination of -136 deg (ss19, 10 mT), one declination of 361.7 deg (ss8), an empty row (ss18), and an extra column in ss48.
- The number of repeated measurements varies from 2 to 19 per step, eleven steps of seven sheets (ss6's 250.1 among them) hold two groups of directions more than 45 deg apart, and the deposit does not say how repeated measurements were combined.
- The glass laboratory field varies between sheets (14.4, 15 and 18 uT); the paper states only 18 uT.
- The Technical Validation example for 14053g1 uses the numbers of 12040g1.
- Table 1's M590 is the moment of the 590.1 measurement alone for all 12 crystals, while the Methods define the TRM as the vector difference of the field-on and field-off steps.
- The deposit's Supplementary Figure 1 plots TRM gained from mean vectors and NRM as the mean of the rows' moments for 12 of its 14 non-synthetic specimens, and the pTRM checks against the field-off step before them. Its points for the four synthetic sheets differ from every combination of their deposited rows by 0.5% or more.
- The crystals' 590 C zero-field moments are all above the 8e-11 emu threshold the paper uses, so the Table 1 maximum intensities rest on a threshold below the level the null state reaches.