Improved kdtree and some refactoring

.gitignore

@@ -16,6 +16,7 @@ dist/
 .DS_Store
 .vscode/
 .idea/
+.copilot-project-context.md
 
 # Generated outputs
 *.ppm

Makefile

@@ -3,13 +3,12 @@
 # Unstructured data visualization tool
 
 CC = gcc
-CFLAGS = -Wall -Wextra -O2 -g
+CFLAGS = -Wall -Wextra -O2 -g -fopenmp
 # --enable-new-dtags is Linux-only, skip on macOS (Darwin)
 UNAME_S := $(shell uname -s)
-ifneq ($(UNAME_S),Darwin)
-LDFLAGS = -Wl,--enable-new-dtags
-else
-LDFLAGS =
+LDFLAGS = -fopenmp
+ifeq ($(UNAME_S),Darwin)
+LDFLAGS += -Wl,--enable-new-dtags
 endif
 
 # NetCDF - try DKRZ system installation first, fall back to system nc-config
@@ -71,7 +70,7 @@ endif
 BASE_CFLAGS = $(CFLAGS) $(NETCDF_CFLAGS)
 X11_FULL_CFLAGS = $(X11_CFLAGS) $(BASE_CFLAGS)
 
-COMMON_LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm
+COMMON_LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm -lpthread
 USHOW_LIBS = $(COMMON_LIBS) $(X11_LIBS) $(X11_RPATH)
 UTERM_LIBS = $(COMMON_LIBS)
 
@@ -179,7 +178,8 @@ SRCDIR = src
 OBJDIR = obj
 BINDIR = .
 
-COMMON_SRCS = $(SRCDIR)/kdtree.c \
+COMMON_SRCS = $(SRCDIR)/common.c \
+              $(SRCDIR)/kdtree.c \
               $(SRCDIR)/mesh.c \
               $(SRCDIR)/regrid.c \
               $(SRCDIR)/projection.c \
@@ -263,29 +263,29 @@ clean:
 	rm -f uterm
 
 # Dependencies
-$(OBJDIR)/ushow.o: $(SRCDIR)/ushow.c $(SRCDIR)/ushow.defines.h $(SRCDIR)/mesh.h \
+$(OBJDIR)/ushow.o: $(SRCDIR)/ushow.c $(SRCDIR)/us_types.h $(SRCDIR)/mesh.h \
                    $(SRCDIR)/regrid.h $(SRCDIR)/file_netcdf.h $(SRCDIR)/colormaps.h \
                    $(SRCDIR)/view.h $(SRCDIR)/interface/x_interface.h
-$(OBJDIR)/uterm.o: $(SRCDIR)/uterm.c $(SRCDIR)/ushow.defines.h $(SRCDIR)/mesh.h \
+$(OBJDIR)/uterm.o: $(SRCDIR)/uterm.c $(SRCDIR)/us_types.h $(SRCDIR)/mesh.h \
                    $(SRCDIR)/regrid.h $(SRCDIR)/file_netcdf.h $(SRCDIR)/colormaps.h \
                    $(SRCDIR)/term_render_mode.h \
                    $(SRCDIR)/view.h
 $(OBJDIR)/term_render_mode.o: $(SRCDIR)/term_render_mode.c $(SRCDIR)/term_render_mode.h
 $(OBJDIR)/kdtree.o: $(SRCDIR)/kdtree.c $(SRCDIR)/kdtree.h
-$(OBJDIR)/mesh.o: $(SRCDIR)/mesh.c $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h
+$(OBJDIR)/mesh.o: $(SRCDIR)/mesh.c $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h
 $(OBJDIR)/regrid.o: $(SRCDIR)/regrid.c $(SRCDIR)/regrid.h $(SRCDIR)/mesh.h \
-                    $(SRCDIR)/kdtree.h $(SRCDIR)/projection.h $(SRCDIR)/ushow.defines.h
-$(OBJDIR)/projection.o: $(SRCDIR)/projection.c $(SRCDIR)/projection.h $(SRCDIR)/ushow.defines.h
-$(OBJDIR)/file_netcdf.o: $(SRCDIR)/file_netcdf.c $(SRCDIR)/file_netcdf.h $(SRCDIR)/ushow.defines.h
-$(OBJDIR)/file_mitgcm.o: $(SRCDIR)/file_mitgcm.c $(SRCDIR)/file_mitgcm.h $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h
-$(OBJDIR)/colormaps.o: $(SRCDIR)/colormaps.c $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h
+                    $(SRCDIR)/kdtree.h $(SRCDIR)/projection.h $(SRCDIR)/us_types.h
+$(OBJDIR)/projection.o: $(SRCDIR)/projection.c $(SRCDIR)/projection.h $(SRCDIR)/us_types.h
+$(OBJDIR)/file_netcdf.o: $(SRCDIR)/file_netcdf.c $(SRCDIR)/file_netcdf.h $(SRCDIR)/us_types.h
+$(OBJDIR)/file_mitgcm.o: $(SRCDIR)/file_mitgcm.c $(SRCDIR)/file_mitgcm.h $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h
+$(OBJDIR)/colormaps.o: $(SRCDIR)/colormaps.c $(SRCDIR)/colormaps.h $(SRCDIR)/us_types.h
 $(OBJDIR)/view.o: $(SRCDIR)/view.c $(SRCDIR)/view.h $(SRCDIR)/file_netcdf.h \
-                  $(SRCDIR)/regrid.h $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h
+                  $(SRCDIR)/regrid.h $(SRCDIR)/colormaps.h $(SRCDIR)/us_types.h
 $(OBJDIR)/interface/x_interface.o: $(SRCDIR)/interface/x_interface.c \
                                     $(SRCDIR)/interface/x_interface.h \
                                     $(SRCDIR)/interface/colorbar.h \
                                     $(SRCDIR)/interface/timeseries_popup.h \
-                                    $(SRCDIR)/ushow.defines.h
+                                    $(SRCDIR)/us_types.h
 $(OBJDIR)/interface/colorbar.o: $(SRCDIR)/interface/colorbar.c \
                                  $(SRCDIR)/interface/colorbar.h $(SRCDIR)/colormaps.h
 $(OBJDIR)/interface/range_popup.o: $(SRCDIR)/interface/range_popup.c \
@@ -295,25 +295,25 @@ $(OBJDIR)/interface/range_utils.o: $(SRCDIR)/interface/range_utils.c \
                                     $(SRCDIR)/interface/range_utils.h
 $(OBJDIR)/interface/timeseries_popup.o: $(SRCDIR)/interface/timeseries_popup.c \
                                          $(SRCDIR)/interface/timeseries_popup.h \
-                                         $(SRCDIR)/ushow.defines.h
+                                         $(SRCDIR)/us_types.h
 
 # Zarr dependencies (when WITH_ZARR is set)
 ifdef WITH_ZARR
-$(OBJDIR)/file_zarr.o: $(SRCDIR)/file_zarr.c $(SRCDIR)/file_zarr.h $(SRCDIR)/ushow.defines.h \
+$(OBJDIR)/file_zarr.o: $(SRCDIR)/file_zarr.c $(SRCDIR)/file_zarr.h $(SRCDIR)/us_types.h \
                         $(SRCDIR)/cJSON/cJSON.h
 $(OBJDIR)/cJSON/cJSON.o: $(SRCDIR)/cJSON/cJSON.c $(SRCDIR)/cJSON/cJSON.h
 endif
 
 # GRIB dependencies (when WITH_GRIB is set)
 ifdef WITH_GRIB
-$(OBJDIR)/file_grib.o: $(SRCDIR)/file_grib.c $(SRCDIR)/file_grib.h $(SRCDIR)/ushow.defines.h
+$(OBJDIR)/file_grib.o: $(SRCDIR)/file_grib.c $(SRCDIR)/file_grib.h $(SRCDIR)/us_types.h
 endif
 
 # YAC dependencies (when WITH_YAC is set)
 ifdef WITH_YAC
 $(OBJDIR)/regrid_yac.o: $(SRCDIR)/regrid_yac.c $(SRCDIR)/regrid_yac.h \
-                         $(SRCDIR)/healpix.h $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h
-$(OBJDIR)/healpix.o: $(SRCDIR)/healpix.c $(SRCDIR)/healpix.h $(SRCDIR)/ushow.defines.h
+                         $(SRCDIR)/healpix.h $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h
+$(OBJDIR)/healpix.o: $(SRCDIR)/healpix.c $(SRCDIR)/healpix.h $(SRCDIR)/us_types.h
 endif
 
 # Print configuration

src/colormaps.h

@@ -5,7 +5,7 @@
 #ifndef COLORMAPS_H
 #define COLORMAPS_H
 
-#include "ushow.defines.h"
+#include "us_types.h"
 
 /*
  * Initialize colormaps (load built-in colormaps).

src/common.c

@@ -0,0 +1,134 @@
+/*
+ * common.c - Shared utilities for ushow and uterm
+ */
+
+#include "common.h"
+#include "kdtree.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define DEFAULT_THREADS 4
+
+void apply_thread_settings(int user_threads) {
+    if (user_threads > 0) {
+        kdtree_set_max_threads(user_threads);
+#ifdef _OPENMP
+        omp_set_num_threads(user_threads);
+#endif
+    } else if (!getenv("OMP_NUM_THREADS")) {
+#ifdef _OPENMP
+        omp_set_num_threads(DEFAULT_THREADS);
+#endif
+    }
+}
+
+int parse_time_units(const char *units, double *unit_seconds,
+                     int *y, int *mo, int *d, int *h, int *mi, double *sec) {
+    if (!units || !unit_seconds || !y || !mo || !d || !h || !mi || !sec) return 0;
+
+    const char *since = strstr(units, "since");
+    if (!since) return 0;
+
+    char unit_buf[32] = {0};
+    if (sscanf(units, "%31s", unit_buf) != 1) return 0;
+
+    /* Normalize unit token to lower case */
+    for (char *p = unit_buf; *p; ++p) {
+        if (*p >= 'A' && *p <= 'Z') *p = (char)(*p - 'A' + 'a');
+    }
+
+    if (strcmp(unit_buf, "seconds") == 0 || strcmp(unit_buf, "second") == 0 ||
+        strcmp(unit_buf, "secs") == 0 || strcmp(unit_buf, "sec") == 0 || strcmp(unit_buf, "s") == 0) {
+        *unit_seconds = 1.0;
+    } else if (strcmp(unit_buf, "minutes") == 0 || strcmp(unit_buf, "minute") == 0 ||
+               strcmp(unit_buf, "mins") == 0 || strcmp(unit_buf, "min") == 0) {
+        *unit_seconds = 60.0;
+    } else if (strcmp(unit_buf, "hours") == 0 || strcmp(unit_buf, "hour") == 0 ||
+               strcmp(unit_buf, "hrs") == 0 || strcmp(unit_buf, "hr") == 0) {
+        *unit_seconds = 3600.0;
+    } else if (strcmp(unit_buf, "days") == 0 || strcmp(unit_buf, "day") == 0) {
+        *unit_seconds = 86400.0;
+    } else {
+        return 0;
+    }
+
+    /* Parse origin date/time after "since" */
+    const char *p = since + 5;
+    while (*p == ' ') p++;
+    int n = sscanf(p, "%d-%d-%d %d:%d:%lf", y, mo, d, h, mi, sec);
+    if (n < 3) return 0;
+    if (n == 3) { *h = 0; *mi = 0; *sec = 0.0; }
+    return 1;
+}
+
+int64_t days_from_civil(int y, unsigned m, unsigned d) {
+    y -= m <= 2;
+    const int era = (y >= 0 ? y : y - 399) / 400;
+    const unsigned yoe = (unsigned)(y - era * 400);
+    const unsigned doy = (153 * (m + (m > 2 ? -3 : 9)) + 2) / 5 + d - 1;
+    const unsigned doe = yoe * 365 + yoe / 4 - yoe / 100 + doy;
+    return (int64_t)(era * 146097 + (int)doe - 719468);
+}
+
+void civil_from_days(int64_t z, int *y, unsigned *m, unsigned *d) {
+    z += 719468;
+    const int era = (z >= 0 ? z : z - 146096) / 146097;
+    const unsigned doe = (unsigned)(z - era * 146097);
+    const unsigned yoe = (doe - doe / 1460 + doe / 36524 - doe / 146096) / 365;
+    int y_tmp = (int)(yoe) + era * 400;
+    const unsigned doy = doe - (365 * yoe + yoe / 4 - yoe / 100);
+    const unsigned mp = (5 * doy + 2) / 153;
+    const unsigned d_tmp = doy - (153 * mp + 2) / 5 + 1;
+    const unsigned m_tmp = mp + (mp < 10 ? 3 : -9);
+    y_tmp += (m_tmp <= 2);
+
+    *y = y_tmp;
+    *m = m_tmp;
+    *d = d_tmp;
+}
+
+int format_time_from_units(char *out, size_t outlen, double value, const char *units) {
+    if (!out || outlen == 0) return 0;
+
+    double unit_seconds = 0.0;
+    int y, mo, d, h, mi;
+    double sec;
+    if (!parse_time_units(units, &unit_seconds, &y, &mo, &d, &h, &mi, &sec)) {
+        return 0;
+    }
+
+    int64_t days = days_from_civil(y, (unsigned)mo, (unsigned)d);
+    double total_sec = (double)days * 86400.0 + (double)h * 3600.0 + (double)mi * 60.0 + sec;
+    total_sec += value * unit_seconds;
+
+    int64_t out_days = (int64_t)(total_sec / 86400.0);
+    double rem = total_sec - (double)out_days * 86400.0;
+    if (rem < 0) {
+        rem += 86400.0;
+        out_days -= 1;
+    }
+
+    int out_y;
+    unsigned out_m, out_d;
+    civil_from_days(out_days, &out_y, &out_m, &out_d);
+
+    int out_h = (int)(rem / 3600.0);
+    rem -= out_h * 3600.0;
+    int out_mi = (int)(rem / 60.0);
+    rem -= out_mi * 60.0;
+    int out_s = (int)(rem + 0.5);
+
+    if (out_h == 0 && out_mi == 0 && out_s == 0) {
+        snprintf(out, outlen, "%04d-%02u-%02u", out_y, out_m, out_d);
+    } else {
+        snprintf(out, outlen, "%04d-%02u-%02u %02d:%02d:%02d",
+                 out_y, out_m, out_d, out_h, out_mi, out_s);
+    }
+    return 1;
+}

src/common.h

@@ -0,0 +1,41 @@
+/*
+ * common.h - Shared utilities for ushow and uterm
+ */
+
+#ifndef COMMON_H
+#define COMMON_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+/*
+ * Apply thread settings for kdtree and OpenMP.
+ *   user_threads > 0: use that value for both kdtree and omp_set_num_threads.
+ *   user_threads <= 0: leave kdtree default; if OMP_NUM_THREADS is unset,
+ *                      set omp_set_num_threads to 4.
+ */
+void apply_thread_settings(int user_threads);
+
+/*
+ * Parse a CF-convention "units since origin" string.
+ * Returns 1 on success, 0 on failure.
+ *   unit_seconds: conversion factor (e.g. 86400 for days)
+ *   y,mo,d,h,mi,sec: parsed origin date/time components
+ */
+int parse_time_units(const char *units, double *unit_seconds,
+                     int *y, int *mo, int *d, int *h, int *mi, double *sec);
+
+/*
+ * Calendar arithmetic helpers (proleptic Gregorian).
+ */
+int64_t days_from_civil(int y, unsigned m, unsigned d);
+void civil_from_days(int64_t z, int *y, unsigned *m, unsigned *d);
+
+/*
+ * Format a CF time value as "YYYY-MM-DD HH:MM:SS" (or "YYYY-MM-DD" when
+ * the time-of-day is exactly midnight).
+ * Returns 1 on success, 0 on failure.
+ */
+int format_time_from_units(char *out, size_t outlen, double value, const char *units);
+
+#endif /* COMMON_H */

src/file_grib.h

@@ -8,7 +8,7 @@
 #ifndef FILE_GRIB_H
 #define FILE_GRIB_H
 
-#include "ushow.defines.h"
+#include "us_types.h"
 
 #ifdef HAVE_GRIB
 

src/file_mitgcm.c

@@ -407,7 +407,7 @@ USFile *mitgcm_open(const char *path) {
     /* Allocate store */
     MitgcmStore *store = calloc(1, sizeof(MitgcmStore));
     if (!store) return NULL;
-    strncpy(store->directory, directory, MAX_NAME_LEN - 1);
+    snprintf(store->directory, sizeof(store->directory), "%s", directory);
     store->nx = nx;
     store->ny = ny;
     store->nz = nz;
@@ -473,7 +473,7 @@ USFile *mitgcm_open(const char *path) {
               free(store); return NULL; }
 
     f->file_type = FILE_TYPE_MITGCM;
-    strncpy(f->filename, directory, MAX_NAME_LEN - 1);
+    snprintf(f->filename, sizeof(f->filename), "%s", directory);
     f->mitgcm_data = store;
 
     return f;
@@ -487,7 +487,7 @@ USMesh *mitgcm_create_mesh(USFile *file) {
     size_t n_points = (size_t)ny * nx;
 
     /* Read XC.data and YC.data */
-    char xc_path[MAX_NAME_LEN], yc_path[MAX_NAME_LEN];
+    char xc_path[MAX_NAME_LEN + 16], yc_path[MAX_NAME_LEN + 16];
     snprintf(xc_path, sizeof(xc_path), "%s/XC.data", store->directory);
     snprintf(yc_path, sizeof(yc_path), "%s/YC.data", store->directory);
 
@@ -637,7 +637,8 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) {
             if (strcmp(prefixes[i], prefix) == 0) { found = 1; break; }
         }
         if (!found && n_prefixes < 64) {
-            strncpy(prefixes[n_prefixes++], prefix, MAX_NAME_LEN - 1);
+            snprintf(prefixes[n_prefixes], MAX_NAME_LEN, "%s", prefix);
+            n_prefixes++;
         }
     }
     closedir(dir);
@@ -659,7 +660,7 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) {
         if (!iterations || n_iters == 0) { free(iterations); continue; }
 
         /* Parse meta from first iteration to get field info */
-        char meta_path[MAX_NAME_LEN];
+        char meta_path[MAX_NAME_LEN * 2 + 32];
         snprintf(meta_path, sizeof(meta_path), "%s/%s.%010d.meta",
                  store->directory, prefixes[pi], iterations[0]);
 
@@ -683,9 +684,9 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) {
 
             /* Name */
             if (minfo.n_fields > 0 && minfo.fields[fi][0]) {
-                strncpy(var->name, minfo.fields[fi], MAX_NAME_LEN - 1);
+                snprintf(var->name, sizeof(var->name), "%s", minfo.fields[fi]);
             } else {
-                strncpy(var->name, prefixes[pi], MAX_NAME_LEN - 1);
+                snprintf(var->name, sizeof(var->name), "%s", prefixes[pi]);
             }
 
             var->mesh = m;
@@ -726,7 +727,7 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) {
             /* Private data */
             MitgcmVarData *vd = calloc(1, sizeof(MitgcmVarData));
             if (!vd) { free(var); continue; }
-            strncpy(vd->prefix, prefixes[pi], MAX_NAME_LEN - 1);
+            snprintf(vd->prefix, sizeof(vd->prefix), "%s", prefixes[pi]);
             vd->field_index = fi;
             vd->iterations = malloc(n_iters * sizeof(int));
             if (vd->iterations) {
@@ -801,7 +802,7 @@ int mitgcm_read_slice(USVar *var, size_t time_idx, size_t depth_idx, float *data
     if ((int)time_idx >= vd->n_iterations) return -1;
 
     /* Build data file path */
-    char data_path[MAX_NAME_LEN];
+    char data_path[MAX_NAME_LEN * 2 + 32];
     snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data",
              store->directory, vd->prefix, vd->iterations[time_idx]);
 
@@ -1054,7 +1055,7 @@ int mitgcm_read_timeseries(USVar *var, size_t node_idx, size_t depth_idx,
         times[t] = (double)vd->iterations[t];
 
         /* Build path and read single value */
-        char data_path[MAX_NAME_LEN];
+        char data_path[MAX_NAME_LEN * 2 + 32];
         snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data",
                  store->directory, vd->prefix, vd->iterations[t]);