-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathtest_geomag.c
135 lines (102 loc) · 3.01 KB
/
test_geomag.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/**
* @file
* @author Nicholas DeCicco <[email protected]>
*/
#include <stdio.h>
#include <math.h>
#include <netcdf.h>
#include "geomag.h"
#define N_LAT 40
#define N_LON 90
#define LIST_VARIABLES(VAR) \
VAR(d) \
VAR(i) \
VAR(h) \
VAR(f) \
VAR(x) \
VAR(y) \
VAR(z) \
VAR(ddot) \
VAR(idot) \
VAR(hdot) \
VAR(fdot) \
VAR(xdot) \
VAR(ydot) \
VAR(zdot)
#define NUM_VARIABLES 14
#define DEF_VAR(x) double x[N_LON][N_LAT];
#define DEFINE_VARIABLES LIST_VARIABLES(DEF_VAR)
#define DEF_NC_VAR(x) int x ## _var;
#define DEFINE_NC_VARS LIST_VARIABLES(DEF_NC_VAR)
#define LON_DIM 0
#define LAT_DIM 1
#define TRYNC(x) if ((status = (x)) != NC_NOERR) goto ncErr;
void print_model(BFieldModel const*const model);
#ifdef TEST_EMBEDDED
extern BFieldModel model;
#endif
int main()
{
#ifndef TEST_EMBEDDED
BFieldModel model;
#endif
BField bfield[N_LON][N_LAT];
// Define arrays for each variable
DEFINE_VARIABLES
// Define
DEFINE_NC_VARS
int ncid;
const double alt = 500.0;
const double date = julday(10, 10, 2016);
double epsilon = 10; /* We can't get the H field at the poles, so only go
within 'epsilon' (degrees) of them. */
int ii, ij;
int dims[2];
int status;
double lat[N_LAT], lon[N_LON];
int lat_var, lon_var;
#ifndef TEST_EMBEDDED
if (!read_model(&model, "IGRF12.COF")) {
fprintf(stderr, "Fatal: failed to open coefficient file\n");
return 0;
}
#endif
#define COPY_VAR(x) \
x[ij][ii] = bfield[ij][ii].x;
for (ii = 0; ii < N_LAT; ii++) {
for (ij = 0; ij < N_LON; ij++) {
lat[ii] = -90.0 + epsilon +
((double) ii)/((double) (N_LAT-1))*(180.0 - 2.0*epsilon);
lon[ij] = -180.0 + ((double) ij)/((double) (N_LON-1))*360.0;
get_field_components(bfield[ij]+ii, &model, alt, kUnitsKilometers,
kCoordSysGeodetic, lat[ii], lon[ij], date);
LIST_VARIABLES(COPY_VAR);
}
}
if ((status = nc_create("test.nc", NC_NETCDF4, &ncid)) != NC_NOERR)
goto ncErr_dontclose;
TRYNC(nc_def_dim(ncid, "lat", (long) N_LAT, dims+LAT_DIM));
TRYNC(nc_def_dim(ncid, "lon", (long) N_LON, dims+LON_DIM));
TRYNC(nc_def_var(ncid, "lat", NC_DOUBLE, 1, dims+LAT_DIM, &lat_var));
TRYNC(nc_def_var(ncid, "lon", NC_DOUBLE, 1, dims+LON_DIM, &lon_var));
const char degrees_north[] = "degrees_north";
const char degrees_east[] = "degrees_east";
TRYNC(nc_put_att_text(ncid, lat_var, "units", sizeof(degrees_north)/sizeof(char), degrees_north));
TRYNC(nc_put_att_text(ncid, lon_var, "units", sizeof(degrees_east)/sizeof(char), degrees_east));
TRYNC(nc_put_var_double(ncid, lat_var, lat));
TRYNC(nc_put_var_double(ncid, lon_var, lon));
#define CREATE_NC_VAR(x) \
TRYNC(nc_def_var(ncid, #x, NC_DOUBLE, 2, dims, & x ## _var));
LIST_VARIABLES(CREATE_NC_VAR);
#define FILL_VAR(x) \
TRYNC(nc_put_var_double(ncid, x ## _var, (double*) x));
LIST_VARIABLES(FILL_VAR);
nc_close(ncid);
return 0;
ncErr:
nc_close(ncid);
ncErr_dontclose:
fprintf(stderr, "NetCDF error occured: \"%s\"", nc_strerror(status));
return 1;
}