From aecf192b4a4941dd9062bc2d39f7cc26621b6ac2 Mon Sep 17 00:00:00 2001 From: Hanno Rein Date: Tue, 9 Apr 2024 19:50:01 -0400 Subject: [PATCH] first version working --- dev_tools/ephem_slicer/Makefile | 2 + dev_tools/ephem_slicer/main.c | 247 ++++++++++++++++++++++++++++++++ 2 files changed, 249 insertions(+) create mode 100644 dev_tools/ephem_slicer/Makefile create mode 100644 dev_tools/ephem_slicer/main.c diff --git a/dev_tools/ephem_slicer/Makefile b/dev_tools/ephem_slicer/Makefile new file mode 100644 index 0000000..024d4af --- /dev/null +++ b/dev_tools/ephem_slicer/Makefile @@ -0,0 +1,2 @@ +all: + gcc main.c -std=c99 -lc -I../../src/ -I../../../rebound/src/ -o ephem_slicer diff --git a/dev_tools/ephem_slicer/main.c b/dev_tools/ephem_slicer/main.c new file mode 100644 index 0000000..f8432be --- /dev/null +++ b/dev_tools/ephem_slicer/main.c @@ -0,0 +1,247 @@ +// Ephem slicer +// (c) Hanno Rein 2024 +// A tool to trim down the DE440 ephemeris file, covering a smaller baseline +// +#include +#include +#include +#include + +#include "planets.h" // for jpl_s struct + +// https://gist.github.com/dgoguerra/7194777 +static const char *humanSize(uint64_t bytes) { + char *suffix[] = {"B", "KB", "MB", "GB", "TB"}; + char length = sizeof(suffix) / sizeof(suffix[0]); + + int i = 0; + double dblBytes = bytes; + + if (bytes > 1024) { + for (i = 0; (bytes / 1024) > 0 && ibeg, sizeof(double)); // Start JD + ret += read(fd, &jpl->end, sizeof(double)); // End JD + ret += read(fd, &jpl->inc, sizeof(double)); // Days per block + ret += read(fd, &jpl->num, sizeof(int32_t)); // Number of constants + ret += read(fd, &jpl->cau, sizeof(double)); // AU to km + ret += read(fd, &jpl->cem, sizeof(double)); // Ratio between Earth/Moon + + // number of coefficients for all components + for (int p = 0; p < JPL_N; p++){ + jpl->ncm[p] = 3; + } + // exceptions: + jpl->ncm[JPL_NUT] = 2; // nutations + jpl->ncm[JPL_TDB] = 1; // TT-TDB + + for (int p = 0; p < 12; p++) { // Columns 1-12 of Group 1050 + ret += read(fd, &jpl->off[p], sizeof(int32_t)); + ret += read(fd, &jpl->ncf[p], sizeof(int32_t)); + ret += read(fd, &jpl->niv[p], sizeof(int32_t)); + } + + ret += read(fd, &jpl->ver, sizeof(int32_t)); // Version. e.g. 440 + ret += read(fd, &jpl->off[12], sizeof(int32_t)); // Columns 13 of Group 1050 + ret += read(fd, &jpl->ncf[12], sizeof(int32_t)); + ret += read(fd, &jpl->niv[12], sizeof(int32_t)); + + // Get all the constant names + jpl->str = calloc(jpl->num, sizeof(char *)); + + // retrieve the names of the first 400 constants + lseek(fd, 0x00FC, SEEK_SET); + for (int p = 0; p < 400; p++) { // Group 1040 + jpl->str[p] = calloc(8, sizeof(char)); + read(fd, jpl->str[p], 6); + } + + // read the remaining constant names + lseek(fd, 0x0B28, SEEK_SET); + for (int p = 400; p < jpl->num; p++) { + jpl->str[p] = calloc(8, sizeof(char)); + read(fd, jpl->str[p], 6); + } + + for (int p = 13; p < 15; p++) { // Columns 14 and 15 of Group 1050 + ret += read(fd, &jpl->off[p], sizeof(int32_t)); + ret += read(fd, &jpl->ncf[p], sizeof(int32_t)); + ret += read(fd, &jpl->niv[p], sizeof(int32_t)); + } + + // adjust for correct indexing (ie: zero based) + for (int p = 0; p < JPL_N; p++){ + jpl->off[p] -= 1; + } + + // save file size, and determine 'kernel size' or 'block size' (=8144 bytes for DE440/441) + jpl->rec = sizeof(double) * 2; + + for (int p = 0; p < JPL_N; p++){ + jpl->rec += sizeof(double) * jpl->ncf[p] * jpl->niv[p] * jpl->ncm[p]; + } + + unsigned int nblocks = (jpl->end-jpl->beg)/jpl->inc; + + printf("Original file (%s):\n", argv[1]); + printf("Days per block: %.4f\n", jpl->inc); + printf("Number of blocks: %u\n", nblocks); + printf("Block size (bytes): %zu\n", jpl->rec); + printf("Start (JD): %.4f ", jpl->beg); + printf_jd(jpl->beg); + printf("\nEnd (JD): %.4f ", jpl->end); + printf_jd(jpl->end); + printf("\nHeader size: %s\n", humanSize((uint64_t)(jpl->rec*2))); + printf("Data size: %s\n", humanSize((uint64_t)(jpl->rec*nblocks))); + + if (argc<4){ + printf("No new file created.\n"); + exit(0); + } + int first_block = atoi(argv[2]); + int num_blocks = atoi(argv[3]); + char* outfilename = malloc(sizeof(char)*1024); + strcpy(outfilename, "linux_"); + double beg = jpl->beg+jpl->inc*first_block; + double end = jpl->beg+jpl->inc*(first_block+num_blocks); + + + + short int beg_year, beg_month, beg_day; + double beg_hour; + short int end_year, end_month, end_day; + double end_hour; + cal_date(beg, &beg_year, &beg_month, &beg_day, &beg_hour); + cal_date(end, &end_year, &end_month, &end_day, &end_hour); + + + if (beg_year<0) { + sprintf(outfilename + strlen(outfilename), "m%04d", -beg_year); + }else{ + sprintf(outfilename + strlen(outfilename), "p%04d", beg_year); + } + if (end_year<0) { + sprintf(outfilename + strlen(outfilename), "m%04d.440", -end_year); + }else{ + sprintf(outfilename + strlen(outfilename), "p%04d.440", end_year); + } + + + FILE* fdo; + if ((fdo = fopen(outfilename, "w")) < 0){ + fprintf(stderr, "Cannot open output file.\n"); + exit(1); + } + + // copy header + lseek(fd, 0, SEEK_SET); + char* buf = malloc(2*jpl->rec); + read(fd, buf, 2*jpl->rec); + fwrite(buf, 2*jpl->rec, 1, fdo); + + printf("\n\nNew file (%s):\n", outfilename); + printf("Start (JD): %.4f ", beg); + printf_jd(beg); + printf("\nEnd (JD): %.4f ", end); + printf_jd(end); + printf("\nHeader size: %s\n", humanSize((uint64_t)(jpl->rec*2))); + printf("Data size: %s\n", humanSize((uint64_t)(jpl->rec*num_blocks))); + + // Copy data blocks + for (int blk = first_block; blkrec, SEEK_SET); + read(fd, buf, jpl->rec); + fwrite(buf, jpl->rec, 1, fdo); + } + + // Overwrite header dates + fseek(fdo, 0x0A5C, SEEK_SET); + fwrite(&beg, sizeof(double), 1, fdo); + fwrite(&end, sizeof(double), 1, fdo); + + if (close(fd) < 0) { + fprintf(stderr, "Error while closing input file.\n"); + exit(1); + } + if (fclose(fdo) < 0) { + fprintf(stderr, "Error while closing output file.\n"); + exit(1); + } + + return 0; + +} +