File: wrapper.cpp

package info (click to toggle)
sgp4 2.22-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 676 kB
  • sloc: python: 2,463; cpp: 2,416; makefile: 6; xml: 5
file content (629 lines) | stat: -rw-r--r-- 22,229 bytes parent folder | download
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
#define PY_SSIZE_T_CLEAN
#include <Python.h>
#include "SGP4.h"
#include "structmember.h"

/* Whether scanf() prefers commas as the decimal point: true means that
   we cannot count on the current locale to interpret numbers correctly
   when the Vallado C++ code calls `sscanf()`. */

static bool switch_locale;

/* Satrec object that wraps a single raw SGP4 struct, and a Satrec array
   that can broadcast into NumPy arrays. */

typedef struct {
    PyObject_HEAD
    elsetrec satrec;
} SatrecObject;

typedef struct {
    PyObject_VAR_HEAD
    elsetrec satrec[0];
} SatrecArrayObject;

/* Support routine that is used to support NumPy array broadcasting for
   both individual satellite objects and also arrays. */

static PyObject *
_vectorized_sgp4(PyObject *args, elsetrec *raw_satrec_array, int imax)
{
    PyObject *jd_arg, *fr_arg, *e_arg, *r_arg, *v_arg;
    Py_buffer jd_buf, fr_buf, e_buf, r_buf, v_buf;
    PyObject *rv = NULL;

    // To prepare for "cleanup:" below.
    jd_buf.buf = fr_buf.buf = e_buf.buf = r_buf.buf = v_buf.buf = NULL;

    if (!PyArg_ParseTuple(args, "OOOOO:sgp4",
                          &jd_arg, &fr_arg, &e_arg, &r_arg, &v_arg))
        return NULL;

    if (PyObject_GetBuffer(jd_arg, &jd_buf, PyBUF_SIMPLE)) goto cleanup;
    if (PyObject_GetBuffer(fr_arg, &fr_buf, PyBUF_SIMPLE)) goto cleanup;
    if (PyObject_GetBuffer(e_arg, &e_buf, PyBUF_WRITABLE)) goto cleanup;
    if (PyObject_GetBuffer(r_arg, &r_buf, PyBUF_WRITABLE)) goto cleanup;
    if (PyObject_GetBuffer(v_arg, &v_buf, PyBUF_WRITABLE)) goto cleanup;

    if (jd_buf.len != fr_buf.len) {
        PyErr_SetString(PyExc_ValueError, "jd and fr must have the same shape");
        goto cleanup;
    }

    // This extra block allows the "goto" statements above to jump
    // across these further variable declarations.
    {
        Py_ssize_t jmax = jd_buf.len / sizeof(double);

        if ((r_buf.len != (Py_ssize_t) sizeof(double) * imax * jmax * 3) ||
            (v_buf.len != (Py_ssize_t) sizeof(double) * imax * jmax * 3) ||
            (e_buf.len != (Py_ssize_t) sizeof(uint8_t) * imax * jmax)) {
            PyErr_SetString(PyExc_ValueError, "bad output array dimension");
            goto cleanup;
        }

        double *jd = (double*) jd_buf.buf;
        double *fr = (double*) fr_buf.buf;
        double *r = (double*) r_buf.buf;
        double *v = (double*) v_buf.buf;
        uint8_t *e = (uint8_t*) e_buf.buf;

#pragma omp parallel for
        for (Py_ssize_t i=0; i < imax; i++) {
            elsetrec &satrec = raw_satrec_array[i];
            for (Py_ssize_t j=0; j < jmax; j++) {
                double t = (jd[j] - satrec.jdsatepoch) * 1440.0
                         + (fr[j] - satrec.jdsatepochF) * 1440.0;
                Py_ssize_t k1 = i * jmax + j;
                Py_ssize_t k3 = 3 * k1;
                SGP4Funcs::sgp4(satrec, t, r + k3, v + k3);
                e[k1] = (uint8_t) satrec.error;
                if (satrec.error && satrec.error < 6) {
                    r[k3] = r[k3+1] = r[k3+2] = v[k3] = v[k3+1] = v[k3+2] = NAN;
                }
            }
        }
    }

    Py_INCREF(Py_None);
    rv = Py_None;
 cleanup:
    if (jd_buf.buf) PyBuffer_Release(&jd_buf);
    if (fr_buf.buf) PyBuffer_Release(&fr_buf);
    if (r_buf.buf) PyBuffer_Release(&r_buf);
    if (v_buf.buf) PyBuffer_Release(&v_buf);
    if (e_buf.buf) PyBuffer_Release(&e_buf);
    return rv;
}

/* Details of the "Satrec" satellite object. */

static PyObject *
Satrec_twoline2rv(PyTypeObject *cls, PyObject *args)
{
    char *string1, *string2, line1[130], line2[130];
    gravconsttype whichconst = wgs72;
    double dummy;

    if (!PyArg_ParseTuple(args,
            "ss|i:twoline2rv",
            &string1, &string2, &whichconst))
        return NULL;

    // Copy both lines, since twoline2rv() might write to both buffers.
    strncpy(line1, string1, 130);
    strncpy(line2, string2, 130);

    // Truncate the line before any checksum characters, if present;
    // otherwise the C++ scanf() interprets each checksum character as
    // part of the preceding value.
    line1[68] = '\0';
    line2[68] = '\0';

    /* Allocate the new object. */
    SatrecObject *self = (SatrecObject*) cls->tp_alloc(cls, 0);
    if (!self)
        return NULL;

    /* Correct for locales that use a comma as the decimal point, since
       users report that the scanf() function on macOS is sensitive to
       locale when parsing floats.  This operation is not thread-safe,
       but we have not released the GIL. */
    float f;
    sscanf("1,5", "%f", &f);
    switch_locale = (f == 1.5);

    char *old_locale = NULL;
    if (switch_locale)
        old_locale = setlocale(LC_NUMERIC, "C");

    /* Leading spaces in a catalog number make scanf() in the official
       code consume the Classification letter as part of the catalog
       number.  (The first character of the International Designator
       then gets consumed as the Classification instead.)  But no
       parsing error is reported, which is bad for users, so let's avoid
       the situation by adding leading zeros ourselves. */
    for (int i=2; i<7; i++) {
        if (line1[i] == ' ') line1[i] = '0';
        if (line2[i] == ' ') line2[i] = '0';
    }

    /* Call the official routine. */
    SGP4Funcs::twoline2rv(line1, line2, ' ', ' ', 'i', whichconst,
                          dummy, dummy, dummy, self->satrec);

    /* Usability bonus: round the fractional day to exactly the eight
       digits specified in the TLE. */
    self->satrec.jdsatepochF = round(self->satrec.jdsatepochF * 1e8) / 1e8;

    /* To avoid having scanf() interpret the "intldesg" as zero or as
       several fields, the C++ library changes spaces to periods and
       underscores.  Let's convert them back to avoid surprising users
       and to match our Python implementation. */
    if (self->satrec.intldesg[0] == '.')
         self->satrec.intldesg[0] = ' ';
    for (int i=1; i<11; i++)
         if (self->satrec.intldesg[i] == '_')
              self->satrec.intldesg[i] = ' ';

    /* Restore previous locale. */
    if (switch_locale)
        setlocale(LC_NUMERIC, old_locale);

    return (PyObject*) self;
}

static PyObject *
Satrec_sgp4init(PyObject *self, PyObject *args)
{
    char satnum_str[6];
    int whichconst;  /* "int" rather than "gravconsttype" so we know size */
    int opsmode;     /* "int" rather than "char" because "C" needs an int */
    long int satnum;
    double epoch, bstar, ndot, nddot;
    double ecco, argpo, inclo, mo, no_kozai, nodeo;

    if (!PyArg_ParseTuple(args, "iCldddddddddd:sgp4init", &whichconst,
                          &opsmode, &satnum, &epoch, &bstar, &ndot, &nddot,
                          &ecco, &argpo, &inclo, &mo, &no_kozai, &nodeo))
        return NULL;

    // See https://www.space-track.org/documentation#tle-alpha5
    if (satnum < 100000) {
        snprintf(satnum_str, 6, "%05ld", satnum);
    } else if (satnum > 339999) {
        PyErr_SetString(PyExc_ValueError, "satellite number cannot exceed "
                        "339999, whose Alpha 5 encoding is 'Z9999'");
        return 0;
    } else {
        char c = 'A' + satnum / 10000 - 10;
        if (c > 'I') c++;
        if (c > 'O') c++;
        satnum_str[0] = c;
        snprintf(satnum_str + 1, 5, "%04ld", satnum % 10000);
    }

    elsetrec &satrec = ((SatrecObject*) self)->satrec;

    SGP4Funcs::sgp4init((gravconsttype) whichconst, opsmode, satnum_str, epoch,
                        bstar, ndot, nddot, ecco, argpo, inclo, mo, no_kozai,
                        nodeo, satrec);

    /* Populate date fields that SGP4Funcs::twoline2rv would set. */
    double whole;
    double fraction = modf(epoch, &whole);
    double whole_jd = whole + 2433281.5;

    /* Go out on a limb: if `epoch` has no decimal digits past the 8
       decimal places stored in a TLE, then assume the user is trying
       to specify an exact decimal fraction. */
    double epoch8 = epoch * 1e8;
    if (round(epoch8) == epoch8) {
        fraction = round(fraction * 1e8) / 1e8;
    }

    satrec.jdsatepoch = whole_jd;
    satrec.jdsatepochF = fraction;

    int y, m, d, H, M;
    double S, jan0jd, jan0fr /* always comes out 0.0 */;
    SGP4Funcs::invjday_SGP4(2433281.5, whole, y, m, d, H, M, S);
    SGP4Funcs::jday_SGP4(y, 1, 0, 0, 0, 0.0, jan0jd, jan0fr);
    satrec.epochyr = y % 100;
    satrec.epochdays = whole_jd - jan0jd + fraction;

    satrec.classification = 'U';  /* default to Unclassified, not '\0' */

    /* Return true, like sgp4init does; check satrec.error for an error code */
    Py_INCREF(Py_None);
    return Py_None;
}

static PyObject *
Satrec_sgp4_tsince(PyObject *self, PyObject *args)
{
    double tsince, r[3], v[3];

    if (!PyArg_ParseTuple(args, "d:sgp4_tsince", &tsince))
        return NULL;
    elsetrec &satrec = ((SatrecObject*) self)->satrec;
    SGP4Funcs::sgp4(satrec, tsince, r, v);
    if (satrec.error && satrec.error < 6)
        r[0] = r[1] = r[2] = v[0] = v[1] = v[2] = NAN;
    return Py_BuildValue("i(fff)(fff)", satrec.error,
                         r[0], r[1], r[2], v[0], v[1], v[2]);
}

static PyObject *
Satrec_sgp4(PyObject *self, PyObject *args)
{
    double jd, fr, r[3], v[3];
    if (!PyArg_ParseTuple(args, "dd:sgp4", &jd, &fr))
        return NULL;
    elsetrec &satrec = ((SatrecObject*) self)->satrec;
    double tsince = (jd - satrec.jdsatepoch) * 1440.0
                  + (fr - satrec.jdsatepochF) * 1440.0;
    SGP4Funcs::sgp4(satrec, tsince, r, v);
    if (satrec.error && satrec.error < 6)
        r[0] = r[1] = r[2] = v[0] = v[1] = v[2] = NAN;
    return Py_BuildValue("i(fff)(fff)", satrec.error,
                         r[0], r[1], r[2], v[0], v[1], v[2]);
}

static PyObject *
Satrec__sgp4(PyObject *self, PyObject *args)
{
    SatrecObject *obj = (SatrecObject*) self;
    elsetrec *raw_satrec_array = &(obj->satrec);
    Py_ssize_t imax = 1;
    return _vectorized_sgp4(args, raw_satrec_array, imax);
}

static PyMethodDef Satrec_methods[] = {
    {"twoline2rv", (PyCFunction)Satrec_twoline2rv, METH_VARARGS | METH_CLASS,
     PyDoc_STR("Initialize the record from two lines of TLE text and an optional gravity constant.")},
    {"sgp4init", (PyCFunction)Satrec_sgp4init, METH_VARARGS,
     PyDoc_STR("Initialize the record from orbital elements.")},
    {"sgp4", (PyCFunction)Satrec_sgp4, METH_VARARGS,
     PyDoc_STR("Given a modified julian date, return position and velocity.")},
    {"_sgp4", (PyCFunction)Satrec__sgp4, METH_VARARGS,
     PyDoc_STR("Given an array of modified julian dates, return position and velocity arrays.")},
    {"sgp4_tsince", (PyCFunction)Satrec_sgp4_tsince, METH_VARARGS,
     PyDoc_STR("Given minutes since epoch, return position and velocity.")},
    {NULL, NULL}
};

#define O(member) offsetof(SatrecObject, satrec.member)

static PyMemberDef Satrec_members[] = {
    /* Listed in the order they appear in a TLE record. */

    {"operationmode", T_CHAR, O(operationmode), READONLY,
     PyDoc_STR("Operation mode: 'a' legacy AFSPC, or 'i' improved.")},
    {"jdsatepoch", T_DOUBLE, O(jdsatepoch), 0,
     PyDoc_STR("Julian date of epoch, day number (see jdsatepochF).")},
    {"jdsatepochF", T_DOUBLE, O(jdsatepochF), 0,
     PyDoc_STR("Julian date of epoch, fraction of day (see jdsatepoch).")},
    {"classification", T_CHAR, O(classification), 0,
     "Usually U=Unclassified, C=Classified, or S=Secret."},
    /* intldesg: inline character array; see Satrec_getset. */
    {"epochyr", T_INT, O(epochyr), 0,
     PyDoc_STR("Year of this element set's epoch (see epochdays). Not set by sgp4init().")},
    {"epochdays", T_DOUBLE, O(epochdays), 0,
     PyDoc_STR("Day of the year of this element set's epoch (see epochyr). Not set by sgp4init().")},
    {"ndot", T_DOUBLE, O(ndot), READONLY,
     PyDoc_STR("Ballistic Coefficient in radians/minute^2.")},
    {"nddot", T_DOUBLE, O(nddot), READONLY,
     PyDoc_STR("Second Derivative of Mean Motion in radians/minute^3.")},
    {"bstar", T_DOUBLE, O(bstar), READONLY,
     PyDoc_STR("Drag Term in inverse Earth radii.")},
    {"ephtype", T_INT, O(ephtype), 0,
     PyDoc_STR("Ephemeris type (should be 0 in published TLEs).")},
    {"elnum", T_LONG, O(elnum), 0,
     PyDoc_STR("Element set number.")},
    {"inclo", T_DOUBLE, O(inclo), READONLY,
     PyDoc_STR("Inclination in radians.")},
    {"nodeo", T_DOUBLE, O(nodeo), READONLY,
     PyDoc_STR("Right ascension of ascending node in radians.")},
    {"ecco", T_DOUBLE, O(ecco), READONLY,
     PyDoc_STR("Eccentricity.")},
    {"argpo", T_DOUBLE, O(argpo), READONLY,
     PyDoc_STR("Argument of perigee in radians.")},
    {"mo", T_DOUBLE, O(mo), READONLY,
     PyDoc_STR("Mean anomaly in radians.")},
    {"no_kozai", T_DOUBLE, O(no_kozai), READONLY,
     PyDoc_STR("Mean motion in radians per minute.")},
    {"revnum", T_LONG, O(revnum), 0,
     PyDoc_STR("Integer revolution number at the epoch.")},

    /* For compatibility with the old struct members, also accept the
       plain name "no". */

    {"no", T_DOUBLE, O(no_kozai), READONLY,
     PyDoc_STR("Alias for the more carefully named ``no_kozai``.")},

    /* Derived values that do not appear explicitly in the TLE. */

    {"method", T_CHAR, O(method), READONLY,
     PyDoc_STR("Method, either 'n' near earth or 'd' deep space.")},
    {"error", T_INT, O(error), READONLY,
     PyDoc_STR("Error code (1-6) documented in sgp4()")},
    {"a", T_DOUBLE, O(a), READONLY,
     PyDoc_STR("semi-major axis")},
    {"altp", T_DOUBLE, O(altp), READONLY,
     PyDoc_STR("altitude of perigee")},
    {"alta", T_DOUBLE, O(alta), READONLY,
     PyDoc_STR("altitude of perigee")},

    /* Single averaged mean elements */

    {"am", T_DOUBLE, O(am), READONLY,
     PyDoc_STR("am: Average semi-major axis")},
    {"em", T_DOUBLE, O(em), READONLY,
     PyDoc_STR("em: Average eccentricity")},
    {"im", T_DOUBLE, O(im), READONLY,
     PyDoc_STR("im: Average inclination")},
    {"Om", T_DOUBLE, O(Om), READONLY,
     PyDoc_STR("Om: Average right ascension of ascending node")},
    {"om", T_DOUBLE, O(om), READONLY,
     PyDoc_STR("om: Average argument of perigee")},
    {"mm", T_DOUBLE, O(mm), READONLY,
     PyDoc_STR("mm: Average mean anomaly")},
    {"nm", T_DOUBLE, O(nm), READONLY,
     PyDoc_STR("nm: Average mean motion")},

    /* Gravity-constant dependent values (initialized by sgp4init() */

    {"tumin", T_DOUBLE, O(tumin), READONLY,
     PyDoc_STR("minutes in one time unit")},
    {"mu", T_DOUBLE, O(mus), READONLY,
     PyDoc_STR("Earth gravitational parameter")},
    {"radiusearthkm", T_DOUBLE, O(radiusearthkm), READONLY,
     PyDoc_STR("radius of the earth in km")},
    {"xke", T_DOUBLE, O(xke), READONLY,
     PyDoc_STR("reciprocal of tumin")},
    {"j2", T_DOUBLE, O(j2), READONLY,
     PyDoc_STR("un-normalized zonal harmonic j2 value")},
    {"j3", T_DOUBLE, O(j3), READONLY,
     PyDoc_STR("un-normalized zonal harmonic j3 value")},
    {"j4", T_DOUBLE, O(j4), READONLY,
     PyDoc_STR("un-normalized zonal harmonic j4 value")},
    {"j3oj2", T_DOUBLE, O(j3oj2), READONLY,
     PyDoc_STR("j3 divided by j2")},

    /* Other convenience variables (some required by propagation.py) */

    {"t", T_DOUBLE, O(t), READONLY,
     PyDoc_STR("Last tsince input to sgp4()")},
    {"mdot", T_DOUBLE, O(mdot), READONLY,
     PyDoc_STR("mean anomaly dot (rate)")},
    {"argpdot", T_DOUBLE, O(argpdot), READONLY,
     PyDoc_STR("argument of perigee dot (rate)")},
    {"nodedot", T_DOUBLE, O(nodedot), READONLY,
     PyDoc_STR("right ascension of ascending node dot (rate)")},
    {"gsto", T_DOUBLE, O(gsto), READONLY,
     PyDoc_STR("gsto: greenwich sidereal time")},

    {NULL}
};

#undef O

static PyObject *
get_intldesg(SatrecObject *self, void *closure)
{
    int length = 0;
    char *s = self->satrec.intldesg;
    while (length <= 7 && s[length])
        length++;
    while (length && s[length-1] == ' ')
        length--;
    return PyUnicode_FromStringAndSize(s, length);
}

static int
set_intldesg(SatrecObject *self, PyObject *value, void *closure)
{
    const char *s = PyUnicode_AsUTF8(value);
    if (!s)
        return -1;
    strncpy(self->satrec.intldesg, s, 10);
    self->satrec.intldesg[10] = '\0';
    return 0;
}

/* The official TLE code has switched `satnum` to a string, but to avoid
   a breaking change we make the string available as `satnum_str` and
   still return an integer for `satnum`. */
static PyObject *
get_satnum(SatrecObject *self, void *closure)
{
    char *s = self->satrec.satnum;
    long n;
    if (strlen(s) < 5 || s[0] <= '9')
        n = atol(s);
    else if (s[0] <= 'I')
        n = (s[0] - 'A' + 10) * 10000 + atol(s + 1);
    else if (s[0] <= 'O')
        n = (s[0] - 'A' + 9) * 10000 + atol(s + 1);
    else
        n = (s[0] - 'A' + 8) * 10000 + atol(s + 1);
    return PyLong_FromLong(n);
}

static PyObject *
get_satnum_str(SatrecObject *self, void *closure)
{
    return PyUnicode_FromString(self->satrec.satnum);
}

static int
set_satnum_str(SatrecObject *self, PyObject *value, void *closure)
{
    const char *s = PyUnicode_AsUTF8(value);
    if (!s)
        return -1;
    strncpy(self->satrec.satnum, s, 5);
    self->satrec.satnum[5] = '\0';
    return 0;
}

static PyGetSetDef Satrec_getset[] = {
    {"intldesg", (getter)get_intldesg, (setter)set_intldesg,
     PyDoc_STR("International Designator: a string of up to 7 characters"
               " from the first line of the TLE that typically provides"
               " two digits for the launch year, a 3-digit launch number,"
               " and one or two letters for which piece of the launch.")},
    {"satnum", (getter)get_satnum, NULL,
     PyDoc_STR("Satellite number from characters 3-7 of each TLE line,"
               " as an integer.")},
    {"satnum_str", (getter)get_satnum_str, (setter)set_satnum_str,
     PyDoc_STR("Satellite number from characters 3-7 of each TLE line,"
               " as a string.")},
    {NULL},
};

static PyTypeObject SatrecType = {
    PyVarObject_HEAD_INIT(NULL, 0)
    /* See the module initialization function at the bottom of this file. */
};

/* Details of the SatrecArray. */

static Py_ssize_t
Satrec_len(PyObject *self) {
    return ((SatrecArrayObject*)self)->ob_base.ob_size;
}

static PySequenceMethods SatrecArray_as_sequence = {
    Satrec_len
};

static PyObject *
SatrecArray_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
{
    PyObject *sequence;
    if (!PyArg_ParseTuple(args, "O:SatrecArray", &sequence))
        return NULL;

    Py_ssize_t length = PySequence_Length(sequence);
    if (length == -1)
        return NULL;

    return type->tp_alloc(type, length);
}


static int
SatrecArray_init(SatrecArrayObject *self, PyObject *args, PyObject *kwds)
{
    PyObject *sequence;
    if (!PyArg_ParseTuple(args, "O:SatrecArray", &sequence))
        return -1;

    Py_ssize_t length = PySequence_Length(sequence);
    if (length == -1)
        return -1;

    for (Py_ssize_t i=0; i<length; i++) {
        PyObject *item = PySequence_GetItem(sequence, i);
        if (!item)
            return -1;
        if (!PyObject_IsInstance(item, (PyObject*) &SatrecType)) {
            PyErr_Format(PyExc_ValueError, "every item must be a Satrec,"
                         " but element %d is: %R", i, item);
            Py_DECREF(item);
            return -1;
        }
        self->satrec[i] = ((SatrecObject *)item)->satrec;
        Py_DECREF(item);
    }
    return 0;
}

static PyObject *
SatrecArray_sgp4(PyObject *self, PyObject *args)
{
    SatrecArrayObject *satrec_array = (SatrecArrayObject*) self;
    elsetrec *raw_satrec_array = &(satrec_array->satrec[0]);
    Py_ssize_t imax = satrec_array->ob_base.ob_size;
    return _vectorized_sgp4(args, raw_satrec_array, imax);
}

static PyMethodDef SatrecArray_methods[] = {
    {"_sgp4", (PyCFunction)SatrecArray_sgp4, METH_VARARGS,
     PyDoc_STR("Given array of minutes since epoch, write"
               " positions, velocities, and errors to arrays.")},
    {NULL, NULL}
};

static PyTypeObject SatrecArrayType = {
    PyVarObject_HEAD_INIT(NULL, sizeof(elsetrec))
    /* See the module initialization function at the bottom of this file. */
};

/* The module that ties it all together. */
static PyModuleDef module = {
    PyModuleDef_HEAD_INIT,
    "sgp4.vallado_cpp",
    "Official C++ SGP4 implementation.",
    -1
};

#include <stdio.h>

PyMODINIT_FUNC
PyInit_vallado_cpp(void)
{
    SatrecType.tp_name = "sgp4.vallado_cpp.Satrec";
    SatrecType.tp_basicsize = sizeof(SatrecObject);
    SatrecArrayType.tp_itemsize = 0;
    SatrecType.tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE;
    SatrecType.tp_doc = "SGP4 satellite record.";
    SatrecType.tp_methods = Satrec_methods;
    SatrecType.tp_members = Satrec_members;
    SatrecType.tp_getset = Satrec_getset;
    SatrecType.tp_new = PyType_GenericNew;

    if (PyType_Ready(&SatrecType) < 0)
        return NULL;

    SatrecArrayType.tp_name = "sgp4.vallado_cpp.SatrecArray";
    SatrecArrayType.tp_basicsize = sizeof(SatrecArrayObject);
    SatrecArrayType.tp_itemsize = sizeof(elsetrec);
    SatrecArrayType.tp_as_sequence = &SatrecArray_as_sequence;
    SatrecArrayType.tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE;
    SatrecArrayType.tp_doc = "SGP4 array of satellites.";
    SatrecArrayType.tp_methods = SatrecArray_methods;
    SatrecArrayType.tp_init = (initproc) SatrecArray_init;
    SatrecArrayType.tp_new = SatrecArray_new;

    if (PyType_Ready(&SatrecArrayType) < 0)
        return NULL;

    PyObject *m = PyModule_Create(&module);
    if (m == NULL)
        return NULL;

    Py_INCREF(&SatrecType);
    if (PyModule_AddObject(m, "Satrec", (PyObject *) &SatrecType) < 0) {
        Py_DECREF(&SatrecType);
        Py_DECREF(m);
        return NULL;
    }

    Py_INCREF(&SatrecArrayType);
    if (PyModule_AddObject(m, "SatrecArray", (PyObject *) &SatrecArrayType) < 0) {
        Py_DECREF(&SatrecArrayType);
        Py_DECREF(&SatrecType);
        Py_DECREF(m);
        return NULL;
    }

    if (PyModule_AddIntConstant(m, "WGS72", (int)wgs72) ||
        PyModule_AddIntConstant(m, "WGS72OLD", (int)wgs72old) ||
        PyModule_AddIntConstant(m, "WGS84", (int)wgs84))
        return NULL;

    return m;
}