Skip to main content
Changed the title to make it more interesting.
Link
pacmaninbw
  • 26.1k
  • 13
  • 47
  • 114

c, cmake, and google test Astronomical calculations in C for python bindings

added 263 characters in body
Source Link
#include "core.h"
#include <stdlib.h>
#include "array.h"

const int MONTHS_IN_YEAR = 12;
const double AVG_DAYS_PER_MONTH = 30.6001;
const double AVG_DAYS_PER_YEAR = 365.25;
const double GREGORIAN_CALENDAR =
    1582 + 10 / MONTHS_IN_YEAR + 15 / AVG_DAYS_PER_YEAR;

/*calculate julian day*/
unsigned int julian_day(int year, int month, int day) {
  double jd, gc_cor, margin;

  if (month < 3) {
    year -= 1;
    month += MONTHS_IN_YEAR;
  }

  gc_cor = 0;
  // date is greater than the GERGORIAN_CALENDAR value apply a correction
  if ((year + month / 12 + day / AVG_DAYS_PER_YEAR) >= GREGORIAN_CALENDAR) {
    margin = (double)((int)(year / 100));
    gc_cor = 2 - margin + (double)((int)(margin / 4));
  }
  //
  if (year >= 0) {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  } else {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year - 0.75)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  }
  return jd;
}

int** array_from_date_range(unsigned int* start_day,
                            unsigned int* end_day,
                            unsigned int interval) {
  unsigned int jd0, jd1, jdx;
  jd0 = julian_day(start_day[0], start_day[1], start_day[2]);
  jd1 = julian_day(end_day[0], end_day[1], end_day[2]);
  if(!interval)interval=1;
  int n = (int)floor((jd1 - jd0, i;)/interval);
  int** arr = allocate_array(n, 5);
  int count=0, i ;
  for (i = 0; i < n; i++) {
    jdx = jd0 + i + interval;+count;
    arr[i][0] = jdx;  // int function(jdx, *pyargs)
    arr[i][1] = jdx;  // int function(jdx, *pyargs)
    arr[i][2] = jdx;  // int function(jdx, *pyargs)
    arr[i][3] = jdx;  // int function(jdx, *pyargs)
    arr[i][4] = jdx;  // int function(jdx, *pyargs)
    count = i + interval;
  }
  return arr;
  free_array(arr, n);
}

#include "../core.h"
#include <Python.h>
#include <limits.h>
#include "gtest/gtest.h" 


namespace CORE {

TEST(Array2D, ArrTest) {
  unsigned int start_day[3] = {2022, 1, 1};
  unsigned int end_day[3] = {2022, 1, 15};
  int** arr = array_from_date_range(start_day, end_day, 01);
  ASSERT_EQ(arr[0][0], 2459580);
  ASSERT_EQ(arr[1][0], 2459581);

  int** arr2 = array_from_date_range(start_day, end_day, 2);
  ASSERT_EQ(arr2[0][0], 2459580);
  ASSERT_EQ(arr2[1][0], 2459582);
}

TEST(JulianDay, JdTest) {
  unsigned int jd = julian_day(2022, 1, 1);
  ASSERT_EQ(jd, 2459580);
}

}  // namespace CORE

edit:

I caught an error with my step

#include "core.h"
#include <stdlib.h>
#include "array.h"

const int MONTHS_IN_YEAR = 12;
const double AVG_DAYS_PER_MONTH = 30.6001;
const double AVG_DAYS_PER_YEAR = 365.25;
const double GREGORIAN_CALENDAR =
    1582 + 10 / MONTHS_IN_YEAR + 15 / AVG_DAYS_PER_YEAR;

/*calculate julian day*/
unsigned int julian_day(int year, int month, int day) {
  double jd, gc_cor, margin;

  if (month < 3) {
    year -= 1;
    month += MONTHS_IN_YEAR;
  }

  gc_cor = 0;
  // date is greater than the GERGORIAN_CALENDAR value apply a correction
  if ((year + month / 12 + day / AVG_DAYS_PER_YEAR) >= GREGORIAN_CALENDAR) {
    margin = (double)((int)(year / 100));
    gc_cor = 2 - margin + (double)((int)(margin / 4));
  }
  //
  if (year >= 0) {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  } else {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year - 0.75)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  }
  return jd;
}

int** array_from_date_range(unsigned int* start_day,
                            unsigned int* end_day,
                            unsigned int interval) {
  unsigned int jd0, jd1, jdx;
  jd0 = julian_day(start_day[0], start_day[1], start_day[2]);
  jd1 = julian_day(end_day[0], end_day[1], end_day[2]);
  int n = jd1 - jd0, i;
  int** arr = allocate_array(n, 5);
  for (i = 0; i < n; i++) {
    jdx = jd0 + i + interval;
    arr[i][0] = jdx;  // int function(jdx, *pyargs)
    arr[i][1] = jdx;  // int function(jdx, *pyargs)
    arr[i][2] = jdx;  // int function(jdx, *pyargs)
    arr[i][3] = jdx;  // int function(jdx, *pyargs)
    arr[i][4] = jdx;  // int function(jdx, *pyargs)

  }
  return arr;
  free_array(arr, n);
}

#include "../core.h"
#include <Python.h>
#include <limits.h>
#include "gtest/gtest.h"

namespace CORE {

TEST(Array2D, ArrTest) {
  unsigned int start_day[3] = {2022, 1, 1};
  unsigned int end_day[3] = {2022, 1, 15};
  int** arr = array_from_date_range(start_day, end_day, 0);
  ASSERT_EQ(arr[0][0], 2459580);
  ASSERT_EQ(arr[1][0], 2459581);
}

TEST(JulianDay, JdTest) {
  unsigned int jd = julian_day(2022, 1, 1);
  ASSERT_EQ(jd, 2459580);
}

}  // namespace CORE
#include "core.h"
#include <stdlib.h>
#include "array.h"

const int MONTHS_IN_YEAR = 12;
const double AVG_DAYS_PER_MONTH = 30.6001;
const double AVG_DAYS_PER_YEAR = 365.25;
const double GREGORIAN_CALENDAR =
    1582 + 10 / MONTHS_IN_YEAR + 15 / AVG_DAYS_PER_YEAR;

/*calculate julian day*/
unsigned int julian_day(int year, int month, int day) {
  double jd, gc_cor, margin;

  if (month < 3) {
    year -= 1;
    month += MONTHS_IN_YEAR;
  }

  gc_cor = 0;
  // date is greater than the GERGORIAN_CALENDAR value apply a correction
  if ((year + month / 12 + day / AVG_DAYS_PER_YEAR) >= GREGORIAN_CALENDAR) {
    margin = (double)((int)(year / 100));
    gc_cor = 2 - margin + (double)((int)(margin / 4));
  }
  //
  if (year >= 0) {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  } else {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year - 0.75)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  }
  return jd;
}

int** array_from_date_range(unsigned int* start_day,
                            unsigned int* end_day,
                            unsigned int interval) {
  unsigned int jd0, jd1, jdx;
  jd0 = julian_day(start_day[0], start_day[1], start_day[2]);
  jd1 = julian_day(end_day[0], end_day[1], end_day[2]);
  if(!interval)interval=1;
  int n = (int)floor((jd1 - jd0)/interval);
  int** arr = allocate_array(n, 5);
  int count=0, i ;
  for (i = 0; i < n; i++) {
    jdx = jd0 +count;
    arr[i][0] = jdx;  // int function(jdx, *pyargs)
    arr[i][1] = jdx;  // int function(jdx, *pyargs)
    arr[i][2] = jdx;  // int function(jdx, *pyargs)
    arr[i][3] = jdx;  // int function(jdx, *pyargs)
    arr[i][4] = jdx;  // int function(jdx, *pyargs)
    count = i + interval;
  }
  return arr;
  free_array(arr, n);
}

#include "../core.h"
#include <Python.h>
#include <limits.h>
#include "gtest/gtest.h" 


namespace CORE {

TEST(Array2D, ArrTest) {
  unsigned int start_day[3] = {2022, 1, 1};
  unsigned int end_day[3] = {2022, 1, 15};
  int** arr = array_from_date_range(start_day, end_day, 1);
  ASSERT_EQ(arr[0][0], 2459580);
  ASSERT_EQ(arr[1][0], 2459581);

  int** arr2 = array_from_date_range(start_day, end_day, 2);
  ASSERT_EQ(arr2[0][0], 2459580);
  ASSERT_EQ(arr2[1][0], 2459582);
}

TEST(JulianDay, JdTest) {
  unsigned int jd = julian_day(2022, 1, 1);
  ASSERT_EQ(jd, 2459580);
}

}  // namespace CORE

edit:

I caught an error with my step

Source Link

c, cmake, and google test for python bindings

I've started teaching myself c/c++ with the intent of being able to writing python bindings to c code. I have a c library used for astronomical calculations given a date lat lon and some other parameters.

In my implementation, I use a date range to run a start stop step for loop over the astronomical calcs (I don't actually call them).

src/core.h

#ifndef CORE_H_  // include guard
#define CORE_H_

#ifdef __cplusplus

extern "C" {
#endif

int** array_from_date_range(unsigned int* start_day,
                            unsigned int* end_day,
                            unsigned int interval);

unsigned int julian_day(int year, int month, int day);

#ifdef __cplusplus
}
#endif

#endif  // CORE_H_

src/core.c

#include "core.h"
#include <stdlib.h>
#include "array.h"

const int MONTHS_IN_YEAR = 12;
const double AVG_DAYS_PER_MONTH = 30.6001;
const double AVG_DAYS_PER_YEAR = 365.25;
const double GREGORIAN_CALENDAR =
    1582 + 10 / MONTHS_IN_YEAR + 15 / AVG_DAYS_PER_YEAR;

/*calculate julian day*/
unsigned int julian_day(int year, int month, int day) {
  double jd, gc_cor, margin;

  if (month < 3) {
    year -= 1;
    month += MONTHS_IN_YEAR;
  }

  gc_cor = 0;
  // date is greater than the GERGORIAN_CALENDAR value apply a correction
  if ((year + month / 12 + day / AVG_DAYS_PER_YEAR) >= GREGORIAN_CALENDAR) {
    margin = (double)((int)(year / 100));
    gc_cor = 2 - margin + (double)((int)(margin / 4));
  }
  //
  if (year >= 0) {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  } else {
    jd = (double)((int)(AVG_DAYS_PER_YEAR * year - 0.75)) +
         (double)((int)(AVG_DAYS_PER_MONTH * (month + 1.0))) + day + 1720994.5 +
         gc_cor;
  }
  return jd;
}

int** array_from_date_range(unsigned int* start_day,
                            unsigned int* end_day,
                            unsigned int interval) {
  unsigned int jd0, jd1, jdx;
  jd0 = julian_day(start_day[0], start_day[1], start_day[2]);
  jd1 = julian_day(end_day[0], end_day[1], end_day[2]);
  int n = jd1 - jd0, i;
  int** arr = allocate_array(n, 5);
  for (i = 0; i < n; i++) {
    jdx = jd0 + i + interval;
    arr[i][0] = jdx;  // int function(jdx, *pyargs)
    arr[i][1] = jdx;  // int function(jdx, *pyargs)
    arr[i][2] = jdx;  // int function(jdx, *pyargs)
    arr[i][3] = jdx;  // int function(jdx, *pyargs)
    arr[i][4] = jdx;  // int function(jdx, *pyargs)

  }
  return arr;
  free_array(arr, n);
}

The // int function(jdx, *pyargs) is just a placeholder my intent is to return a 2d array into pandas dataframe. Something along these lines...

pd.DataFrame(array_from_date_range((2022,1,1),(2022,1,15),0,**kwargs),columns=['jdate',...]).set_index('jdate')`

src/array.h

#ifndef ARRAY_H_  // include guard
#define ARRAY_H_

#ifdef __cplusplus
/**/
extern "C" {
#endif

int** allocate_array(int Rows, int Cols);
void free_array(int** board, int Rows);

#ifdef __cplusplus
}
#endif

#endif  // ARRAY_H_

src/array.c

#include <stdlib.h>

int** allocate_array(int nRows, int nCols) {
  // allocate Rows rows, each row is a pointer to int
  int** arr = (int**)malloc(nRows * sizeof(int*));
  int row;

  // for each row allocate Cols ints
  for (row = 0; row < nRows; row++) {
    arr[row] = (int*)malloc(nCols * sizeof(int));
  }

  return arr;
}

void free_array(int** arr, int nRows) {
  int row;
  // first free each row
  for (row = 0; row < nRows; row++) {
    free(arr[row]);
  }
  // Eventually free the memory of the pointers to the rows
  free(arr);
}

I am using a cmake and google test

CMakeLists.txt

cmake_minimum_required(VERSION 3.5)
#
project(cproject)
set(CMAKE_CXX_STANDARD 14)
# ##############################################################################
# Python
# ##############################################################################
find_package(Python 3.10 REQUIRED Interpreter Development NumPy)
find_package(PythonLibs REQUIRED)
# ##############################################################################
# Source lib
# ##############################################################################
add_library(core SHARED src/core.c src/array.c)
target_link_libraries(core PUBLIC Python::NumPy)
# ##############################################################################
# GTest
# ##############################################################################
enable_testing()
find_package(GTest CONFIG REQUIRED)
add_executable(
  core_test
  src/tests/core_test.cc
  src/tests/array_test.cc
)
target_link_libraries(core_test PUBLIC core ${PYTHON_LIBRARIES} GTest::gtest GTest::gtest_main )

add_test(UnitTests core_test)
# ##############################################################################
# Installation
# ##############################################################################
set(CMAKE_INSTALL_PREFIX /usr/local/src)
install(TARGETS core DESTINATION ${CMAKE_INSTALL_PREFIX})

src/tests/core_tests.cc

#include "../core.h"
#include <Python.h>
#include <limits.h>
#include "gtest/gtest.h"

namespace CORE {

TEST(Array2D, ArrTest) {
  unsigned int start_day[3] = {2022, 1, 1};
  unsigned int end_day[3] = {2022, 1, 15};
  int** arr = array_from_date_range(start_day, end_day, 0);
  ASSERT_EQ(arr[0][0], 2459580);
  ASSERT_EQ(arr[1][0], 2459581);
}

TEST(JulianDay, JdTest) {
  unsigned int jd = julian_day(2022, 1, 1);
  ASSERT_EQ(jd, 2459580);
}

}  // namespace CORE