365 lines
17 KiB
Bash
Executable File
365 lines
17 KiB
Bash
Executable File
#!/bin/sh
|
|
# Generates tests/fixtures/helper-cases.json: reference values for the small functions a
|
|
# port of NRLMSISE-00 is built from, so a failure can be localized to one of them rather
|
|
# than only showing up in the nine densities at the end.
|
|
#
|
|
# Fetches Dominik Brodowski's C translation of the NRL Fortran and compiles a driver that
|
|
# #includes the implementation rather than linking against it. That is the point: the shared
|
|
# state the model keeps in file-scope statics — gravity, the Earth radius, the mesosphere
|
|
# temperature nodes — is then visible to the driver, which can set it explicitly and call
|
|
# any helper at a known state. Linking cannot do that.
|
|
#
|
|
# What is committed is the numbers, which are facts about the model rather than anyone's
|
|
# code. The model itself is the work of Picone, Hedin and Drob at the US Naval Research
|
|
# Laboratory. Needs a C compiler and network access; rerun it only to regenerate.
|
|
#
|
|
# Usage: scripts/generate-helper-cases.sh
|
|
|
|
set -e
|
|
|
|
root=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd)
|
|
out="$root/tests/fixtures/helper-cases.json"
|
|
work=$(mktemp -d)
|
|
trap 'rm -rf "$work"' EXIT
|
|
|
|
base="https://raw.githubusercontent.com/magnific0/nrlmsise-00/master"
|
|
echo "fetching the reference from $base"
|
|
for f in nrlmsise-00.c nrlmsise-00.h nrlmsise-00_data.c; do
|
|
curl -sS --fail --max-time 120 -L "$base/$f" -o "$work/$f"
|
|
done
|
|
digest=$(cat "$work/nrlmsise-00.c" "$work/nrlmsise-00_data.c" | sha256sum | cut -d' ' -f1)
|
|
echo "source digest $digest"
|
|
|
|
cat > "$work/driver.c" <<'DRIVER'
|
|
/* Prints reference values for the helpers, at seventeen significant digits.
|
|
*
|
|
* The implementation is #included so that its statics are in this translation unit and can
|
|
* be set before a helper is called. printf is redefined first: dnet prints a message on a
|
|
* non-positive density, and one of the cases below is there to exercise exactly that. */
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#define printf(...) ((void)0)
|
|
#include "nrlmsise-00.c"
|
|
#undef printf
|
|
|
|
#define N(x) printf("%.17e", (double)(x))
|
|
|
|
static int first_case = 1;
|
|
|
|
static void start(const char *what) {
|
|
if (!first_case) printf(",\n");
|
|
first_case = 0;
|
|
printf(" {\"helper\": \"%s\"", what);
|
|
}
|
|
|
|
static void arg(const char *name, double v) {
|
|
printf(", \"%s\": ", name);
|
|
N(v);
|
|
}
|
|
|
|
static void array(const char *name, const double *v, int n) {
|
|
int i;
|
|
printf(", \"%s\": [", name);
|
|
for (i = 0; i < n; i++) { if (i) printf(", "); N(v[i]); }
|
|
printf("]");
|
|
}
|
|
|
|
static void end(void) { printf("}"); }
|
|
|
|
/* Latitudes to take gravity and the Earth radius at. */
|
|
static const double LATS[] = {-90.0, -60.0, -45.0, -23.5, 0.0, 23.5, 45.0, 60.0, 90.0, 30.7};
|
|
|
|
/* The node sets the model actually splines over, and a non-uniform one. */
|
|
static const double ZN1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
|
|
static const double ZN2[4] = {72.5, 55.0, 45.0, 32.5};
|
|
static const double ZN3[5] = {32.5, 20.0, 15.0, 10.0, 0.0};
|
|
|
|
int main(void) {
|
|
int i, k;
|
|
double gv, reff, tz, value;
|
|
double xs[5], ys[5], y2[5];
|
|
|
|
printf("[\n");
|
|
|
|
/* glatf: pure. */
|
|
for (i = 0; i < (int)(sizeof LATS / sizeof *LATS); i++) {
|
|
glatf(LATS[i], &gv, &reff);
|
|
start("glatf"); arg("lat", LATS[i]); arg("gravity", gv); arg("radius", reff); end();
|
|
}
|
|
|
|
/* zeta and scalh: read gsurf and re, which are set here rather than by a model run. */
|
|
glatf(45.0, &gsurf, &re);
|
|
{
|
|
const double zz[] = {0.0, 10.0, 72.5, 120.0, 400.0, 1000.0};
|
|
const double zl[] = {0.0, 32.5, 120.0};
|
|
for (i = 0; i < 6; i++) for (k = 0; k < 3; k++) {
|
|
start("zeta"); arg("gsurf", gsurf); arg("re", re);
|
|
arg("zz", zz[i]); arg("zl", zl[k]); arg("value", zeta(zz[i], zl[k])); end();
|
|
}
|
|
const double masses[] = {4.0, 16.0, 28.0, 32.0};
|
|
const double temps[] = {180.0, 300.0, 1000.0, 1500.0};
|
|
for (i = 0; i < 4; i++) for (k = 0; k < 4; k++) {
|
|
start("scalh"); arg("gsurf", gsurf); arg("re", re);
|
|
arg("alt", 400.0 + 100.0 * i); arg("xm", masses[i]); arg("temp", temps[k]);
|
|
arg("value", scalh(400.0 + 100.0 * i, masses[i], temps[k])); end();
|
|
}
|
|
}
|
|
|
|
/* ccor and ccor2: pure, and both have two clamped branches to reach. */
|
|
{
|
|
const double alts[] = {0.0, 72.5, 100.0, 200.0, 500.0, 1.0e5, -1.0e5};
|
|
const double rs[] = {-0.5, 0.0, 0.25, 2.0};
|
|
for (i = 0; i < 7; i++) for (k = 0; k < 4; k++) {
|
|
start("ccor"); arg("alt", alts[i]); arg("r", rs[k]); arg("h1", 20.0);
|
|
arg("zh", 100.0); arg("value", ccor(alts[i], rs[k], 20.0, 100.0)); end();
|
|
start("ccor2"); arg("alt", alts[i]); arg("r", rs[k]); arg("h1", 20.0);
|
|
arg("zh", 100.0); arg("h2", 5.0);
|
|
arg("value", ccor2(alts[i], rs[k], 20.0, 100.0, 5.0)); end();
|
|
}
|
|
}
|
|
|
|
/* dnet: the blend, both exponent clamps, and the non-positive guard. */
|
|
{
|
|
const double dd[] = {1.0e7, 1.0e7, 1.0e7, 1.0e-30, 1.0e30, 0.0, 1.0e7, 0.0};
|
|
const double dm[] = {1.0e7, 1.0e5, 1.0e9, 1.0e7, 1.0e7, 1.0e7, 0.0, 0.0};
|
|
for (i = 0; i < 8; i++) {
|
|
start("dnet"); arg("dd", dd[i]); arg("dm", dm[i]); arg("zhm", 4.0e4);
|
|
arg("xmm", 28.95); arg("xm", 16.0);
|
|
arg("value", dnet(dd[i], dm[i], 4.0e4, 28.95, 16.0)); end();
|
|
}
|
|
}
|
|
|
|
/* spline, splint and splini, over the model's own node sets. */
|
|
{
|
|
const double *nodes[3]; int counts[3];
|
|
nodes[0] = ZN1; counts[0] = 5;
|
|
nodes[1] = ZN2; counts[1] = 4;
|
|
nodes[2] = ZN3; counts[2] = 5;
|
|
for (i = 0; i < 3; i++) {
|
|
int n = counts[i];
|
|
/* The abscissae the model splines over are normalized geopotential heights,
|
|
* ascending, and the ordinates are inverse temperatures. */
|
|
for (k = 0; k < n; k++) {
|
|
xs[k] = (double)k / (double)(n - 1);
|
|
ys[k] = 1.0 / (200.0 + 15.0 * k - 0.4 * k * k);
|
|
}
|
|
double yd1 = -1.5e-5, yd2 = 3.0e-6;
|
|
spline(xs, ys, n, yd1, yd2, y2);
|
|
start("spline"); array("xs", xs, n); array("ys", ys, n);
|
|
arg("yd1", yd1); arg("yd2", yd2); array("y2", y2, n); end();
|
|
for (k = 0; k <= 8; k++) {
|
|
double x = -0.1 + 0.15 * k;
|
|
splint(xs, ys, y2, n, x, &value);
|
|
start("splint"); array("xs", xs, n); array("ys", ys, n);
|
|
array("y2", y2, n); arg("x", x); arg("value", value); end();
|
|
splini(xs, ys, y2, n, x, &value);
|
|
start("splini"); array("xs", xs, n); array("ys", ys, n);
|
|
array("y2", y2, n); arg("x", x); arg("value", value); end();
|
|
}
|
|
}
|
|
}
|
|
|
|
/* densu: needs gsurf, re, the tn1 nodes and the tgn1 gradients. */
|
|
glatf(60.0, &gsurf, &re);
|
|
meso_tn1[1] = 365.0; meso_tn1[2] = 240.0; meso_tn1[3] = 195.0; meso_tn1[4] = 208.0;
|
|
meso_tgn1[0] = 0.0; meso_tgn1[1] = 8.0;
|
|
{
|
|
const double alts[] = {0.0, 50.0, 72.5, 90.0, 110.0, 120.0, 200.0, 400.0, 1000.0};
|
|
const double masses[] = {0.0, 4.0, 16.0, 28.0};
|
|
const double alphas[] = {0.0, -0.38, 0.0, 0.0};
|
|
double zn1[5];
|
|
for (k = 0; k < 5; k++) zn1[k] = ZN1[k];
|
|
for (i = 0; i < 9; i++) for (k = 0; k < 4; k++) {
|
|
double tn1[5], tgn1[2];
|
|
int j;
|
|
for (j = 0; j < 5; j++) tn1[j] = meso_tn1[j];
|
|
tgn1[0] = meso_tgn1[0]; tgn1[1] = meso_tgn1[1];
|
|
tz = 0.0;
|
|
value = densu(alts[i], 3.0e11, 1250.0, 386.0, masses[k], alphas[k], &tz,
|
|
120.0, 0.02, 5, zn1, tn1, tgn1);
|
|
start("densu"); arg("gsurf", gsurf); arg("re", re);
|
|
arg("alt", alts[i]); arg("dlb", 3.0e11); arg("tinf", 1250.0); arg("tlb", 386.0);
|
|
arg("xm", masses[k]); arg("alpha", alphas[k]); arg("zlb", 120.0);
|
|
arg("s2", 0.02); array("zn1", zn1, 5);
|
|
array("tn1_in", tn1, 5); array("tgn1_in", tgn1, 2);
|
|
arg("value", value); arg("tz", tz);
|
|
/* densu writes the Bates temperature and gradient back into node zero. */
|
|
array("tn1_out", tn1, 5); array("tgn1_out", tgn1, 2); end();
|
|
}
|
|
}
|
|
|
|
/* densm: needs the tn2 and tn3 nodes and their gradients. */
|
|
{
|
|
double tn2[4] = {208.0, 260.0, 268.0, 228.0};
|
|
double tgn2[2] = {8.0, -3.0};
|
|
double tn3[5] = {228.0, 218.0, 214.0, 222.0, 288.0};
|
|
double tgn3[2] = {-3.0, -6.0};
|
|
double zn2[4], zn3[5];
|
|
const double alts[] = {0.0, 5.0, 10.0, 20.0, 32.5, 45.0, 60.0, 72.5, 80.0};
|
|
const double masses[] = {0.0, 28.95};
|
|
for (k = 0; k < 4; k++) zn2[k] = ZN2[k];
|
|
for (k = 0; k < 5; k++) zn3[k] = ZN3[k];
|
|
for (i = 0; i < 9; i++) for (k = 0; k < 2; k++) {
|
|
double a2[4], g2[2], a3[5], g3[2];
|
|
int j;
|
|
for (j = 0; j < 4; j++) a2[j] = tn2[j];
|
|
for (j = 0; j < 5; j++) a3[j] = tn3[j];
|
|
g2[0] = tgn2[0]; g2[1] = tgn2[1]; g3[0] = tgn3[0]; g3[1] = tgn3[1];
|
|
tz = 0.0;
|
|
value = densm(alts[i], 1.0e19, masses[k], &tz, 5, zn3, a3, g3, 4, zn2, a2, g2);
|
|
start("densm"); arg("gsurf", gsurf); arg("re", re);
|
|
arg("alt", alts[i]); arg("d0", 1.0e19); arg("xm", masses[k]);
|
|
array("zn3", zn3, 5); array("zn2", zn2, 4);
|
|
array("tn3", a3, 5); array("tgn3", g3, 2);
|
|
array("tn2", a2, 4); array("tgn2", g2, 2);
|
|
arg("value", value); arg("tz", tz); end();
|
|
}
|
|
}
|
|
|
|
/* globe7 and glob7s: the expansions, over every coefficient table each is called
|
|
* with. They share the LPOLY state: globe7 writes the Legendre polynomials, the local
|
|
* time harmonics, the flux departure and the magnetic activity functions, and glob7s
|
|
* reads them. So each glob7s case records the state it was called at, which is
|
|
* whatever the preceding globe7 left. */
|
|
{
|
|
struct nrlmsise_input in[6];
|
|
struct nrlmsise_flags flags;
|
|
struct ap_array aph;
|
|
int c, s;
|
|
/* The tables globe7 is called with: pt, ps, then the nine rows of pd. */
|
|
double *g7[11];
|
|
const char *g7name[11] = {"pt", "ps", "pd0", "pd1", "pd2", "pd3", "pd4",
|
|
"pd5", "pd6", "pd7", "pd8"};
|
|
/* The tables glob7s is called with: the four rows of ptl, then the ten of pma. */
|
|
double *g7s[14];
|
|
const char *g7sname[14] = {"ptl0", "ptl1", "ptl2", "ptl3", "pma0", "pma1", "pma2",
|
|
"pma3", "pma4", "pma5", "pma6", "pma7", "pma8", "pma9"};
|
|
g7[0] = pt; g7[1] = ps;
|
|
for (i = 0; i < 9; i++) g7[2 + i] = pd[i];
|
|
for (i = 0; i < 4; i++) g7s[i] = ptl[i];
|
|
for (i = 0; i < 10; i++) g7s[4 + i] = pma[i];
|
|
|
|
for (i = 0; i < 7; i++) aph.a[i] = 40.0 + 20.0 * i;
|
|
|
|
/* A spread of input states: the standard case, then the equator, local midnight,
|
|
* a storm, low flux and the southern hemisphere in the other season. */
|
|
for (i = 0; i < 6; i++) {
|
|
in[i].year = 0; in[i].doy = 172; in[i].sec = 29000.0; in[i].alt = 400.0;
|
|
in[i].g_lat = 60.0; in[i].g_long = -70.0; in[i].lst = 16.0;
|
|
in[i].f107A = 150.0; in[i].f107 = 150.0; in[i].ap = 4.0; in[i].ap_a = &aph;
|
|
}
|
|
in[1].g_lat = 0.0;
|
|
in[2].lst = 0.0; in[2].sec = 0.0; in[2].g_long = 0.0;
|
|
in[3].ap = 207.0; in[3].f107 = 230.0;
|
|
in[4].f107A = 70.0; in[4].f107 = 68.0;
|
|
in[5].doy = 15; in[5].g_lat = -45.0; in[5].g_long = 150.0; in[5].lst = 9.5;
|
|
|
|
/* Switch sets: everything on, the daily-ap mode, cross terms only on the
|
|
* seasonal terms, and the diurnal family off. */
|
|
for (s = 0; s < 4; s++) {
|
|
for (k = 0; k < 24; k++) flags.switches[k] = 1;
|
|
flags.switches[0] = 0;
|
|
if (s == 0) flags.switches[9] = -1;
|
|
if (s == 2) { flags.switches[3] = 2; flags.switches[5] = 2; flags.switches[9] = -1; }
|
|
if (s == 3) { flags.switches[7] = 0; flags.switches[8] = 0;
|
|
flags.switches[14] = 0; flags.switches[9] = 0; }
|
|
tselec(&flags);
|
|
|
|
for (c = 0; c < 6; c++) {
|
|
for (i = 0; i < 11; i++) {
|
|
double value = globe7(g7[i], &in[c], &flags);
|
|
start("globe7");
|
|
printf(", \"table\": \"%s\", \"switch_set\": %d, \"case\": %d",
|
|
g7name[i], s, c);
|
|
array("switches_in", (double[24]){
|
|
flags.switches[0], flags.switches[1], flags.switches[2],
|
|
flags.switches[3], flags.switches[4], flags.switches[5],
|
|
flags.switches[6], flags.switches[7], flags.switches[8],
|
|
flags.switches[9], flags.switches[10], flags.switches[11],
|
|
flags.switches[12], flags.switches[13], flags.switches[14],
|
|
flags.switches[15], flags.switches[16], flags.switches[17],
|
|
flags.switches[18], flags.switches[19], flags.switches[20],
|
|
flags.switches[21], flags.switches[22], flags.switches[23]}, 24);
|
|
arg("day_of_year", in[c].doy);
|
|
arg("seconds_of_day", in[c].sec);
|
|
arg("latitude_deg", in[c].g_lat);
|
|
arg("longitude_deg", in[c].g_long);
|
|
arg("local_solar_time_h", in[c].lst);
|
|
arg("f107_average", in[c].f107A);
|
|
arg("f107", in[c].f107);
|
|
arg("ap_daily", in[c].ap);
|
|
array("ap_history", aph.a, 7);
|
|
arg("value", value);
|
|
/* The state it leaves behind, which glob7s then reads. */
|
|
array("plg0", plg[0], 9); array("plg1", plg[1], 9);
|
|
array("plg2", plg[2], 9); array("plg3", plg[3], 9);
|
|
arg("stloc", stloc); arg("ctloc", ctloc);
|
|
arg("s2tloc", s2tloc); arg("c2tloc", c2tloc);
|
|
arg("s3tloc", s3tloc); arg("c3tloc", c3tloc);
|
|
arg("dfa", dfa); arg("apdf", apdf); array("apt", apt, 4);
|
|
end();
|
|
}
|
|
/* globe7 with pt again, so the state the glob7s cases record is the one a
|
|
* model run would have left. */
|
|
(void) globe7(pt, &in[c], &flags);
|
|
for (i = 0; i < 14; i++) {
|
|
double value = glob7s(g7s[i], &in[c], &flags);
|
|
start("glob7s");
|
|
printf(", \"table\": \"%s\", \"switch_set\": %d, \"case\": %d",
|
|
g7sname[i], s, c);
|
|
array("switches_in", (double[24]){
|
|
flags.switches[0], flags.switches[1], flags.switches[2],
|
|
flags.switches[3], flags.switches[4], flags.switches[5],
|
|
flags.switches[6], flags.switches[7], flags.switches[8],
|
|
flags.switches[9], flags.switches[10], flags.switches[11],
|
|
flags.switches[12], flags.switches[13], flags.switches[14],
|
|
flags.switches[15], flags.switches[16], flags.switches[17],
|
|
flags.switches[18], flags.switches[19], flags.switches[20],
|
|
flags.switches[21], flags.switches[22], flags.switches[23]}, 24);
|
|
arg("day_of_year", in[c].doy);
|
|
arg("longitude_deg", in[c].g_long);
|
|
array("plg0", plg[0], 9); array("plg1", plg[1], 9);
|
|
array("plg2", plg[2], 9); array("plg3", plg[3], 9);
|
|
arg("stloc", stloc); arg("ctloc", ctloc);
|
|
arg("s2tloc", s2tloc); arg("c2tloc", c2tloc);
|
|
arg("s3tloc", s3tloc); arg("c3tloc", c3tloc);
|
|
arg("dfa", dfa); arg("apdf", apdf); array("apt", apt, 4);
|
|
arg("value", value);
|
|
end();
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
printf("\n]\n");
|
|
return 0;
|
|
}
|
|
DRIVER
|
|
|
|
cc -O2 -o "$work/driver" "$work/driver.c" "$work/nrlmsise-00_data.c" -lm -I"$work"
|
|
mkdir -p "$(dirname "$out")"
|
|
{
|
|
echo "{"
|
|
echo " \"what\": \"Reference values for the helper functions NRLMSISE-00 is built from.\","
|
|
echo " \"model\": \"NRLMSISE-00, by Picone, Hedin and Drob at the US Naval Research Laboratory\","
|
|
echo " \"reference\": \"$base, Dominik Brodowski's C translation of the NRL Fortran\","
|
|
echo " \"reference_sha256\": \"$digest\","
|
|
echo " \"generated\": \"$(date -u +%Y-%m-%d)\","
|
|
echo " \"generated_by\": \"scripts/generate-helper-cases.sh\","
|
|
echo " \"units\": \"The model's own: km, K, cm^-3, g/cm^3, and cm/s^2 for gravity.\","
|
|
echo " \"note\": \"gsurf and re are the shared state the helper was called at, not arguments to it.\","
|
|
echo " \"cases\":"
|
|
"$work/driver" | sed 's/^/ /'
|
|
echo "}"
|
|
} > "$out"
|
|
|
|
python3 -c "
|
|
import json
|
|
from collections import Counter
|
|
cases = json.load(open('$out'))['cases']
|
|
counts = Counter(c['helper'] for c in cases)
|
|
print(f'wrote tests/fixtures/helper-cases.json: {len(cases)} cases')
|
|
for name, n in sorted(counts.items()):
|
|
print(f' {name:>8}: {n}')
|
|
"
|