Sunday, January 2, 2022

How to write a basic C program to calculate ecliptic longitudinal position of astrological planets

In the title I speak of "astrological" planets. As you know, Sun and Moon aren't planets in astronomy, neither is Pluto, which has been declassified as a planet in recent times. In Astrology, planet means "moving object" and this definition includes Sun, Moon, the solar system planets plus Pluto. Asteroids like Vesta, Chiron, Juno and many other are called asteroids as well in Astrology. Some object are instead considered "fixed", the Zodiac as a whole, and the most visible stars. Other objects are tipically astrological, such as the Ascendant, House cusps, the midpoints, the Arabian parts and some other less frequently used elements. All the related calculus techniques are based on astrological assumptions.

In this post we'll gather all those elements which have the purpose to let us calculate, within the available accuracy and precision, the longitudinal position of planets, i.e. the one and only spatial coordinate that astrology take into account. When you see the zodiac circle in a natal theme, don't think of it as a planar representation of universe. Instead, think of it as a monodimentional representation, like a segment, which has been coiled on itself to join the head and the tail. I'll talk about this argument later, about the inconsistency of the actual methods used to calculate house cusps.

Basic C program and dependencies

Jumping to the start, we must gather a few essential resources:

  1. One bsp file, from this site https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/

  2. de441.bsp is huge, but we can extract a subset. Don't use your phone to download it, it's more than 3GB large, use your PC.

  3. naif0012.tls file, is a generic kernel which contains the annual leapseconds as they are established by IERS (International Earth Rotation and Reference System Service).You can download it at this link. If you want to download it directly in your phone, use the following command under Termux in your Termux home directory:

    wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls

  4. The file mytrueepoch.tf, we use it to define a non inertial system of coordinates, in essence we refer longitudinal coordinates to an estimated vernal point as in traditional occidental tropical astrology. You can create a file in your home Termux directory in your phone:

    touch mytrueepoch.tf
    

    then copy the following code inside of it:

    KPL/FK 
    
    	private Kernel FK, dynamic, geocentric, time-based. It 
    integrates precession, nutation and ecliptic obliquity models 
    
    \begindata
    
       FRAME_MYTRUEEPOCH = 1987654
       FRAME_1987654_NAME = 'MYTRUEEPOCH'
       FRAME_1987654_CLASS = 5
       FRAME_1987654_CLASS_ID = 1987654
       FRAME_1987654_CENTER = 399
       FRAME_1987654_RELATIVE = 'J2000'
       FRAME_1987654_DEF_STYLE = 'PARAMETERIZED'
       FRAME_1987654_FAMILY = 'MEAN_ECLIPTIC_AND_EQUINOX_OF_DATE'
       FRAME_1987654_PREC_MODEL = 'EARTH_IAU_1976'
       FRAME_1987654_OBLIQ_MODEL = 'EARTH_IAU_1980'
       FRAME_1987654_ROTATION_STATE= 'ROTATING'
    
    \begintext
    
  5. obviously, the library file libcspice.so that you have generated from the jpl toolkit

Once collected the resources (the other important one is in the cspice dir) you are almost ready. Let's go back to the big file you downloaded in yor PC: de441.bsp

As large as it is, it's of no use. We must extract a subset based on a useful date interval (from an astrological point of view, of course) then transfer that newly created bsp file to your phone memory.

While the download is going on, you can take a little time to install a graphical scp. For Windows your winning candidate is WinSCP, simple and fast, for Linux, if you don't like scp and the command line, you can use the connection to server ssh that your file manager (nautilus, caja, nemo, etc) provides.

When the download of the monster de441.bsp ends, you have got in your download directory the mother ephemeris, from which you can extract any time interval you like. For my astrological purposes, an interval of time between 1500 AD and 2100 AD is largely adequate. Now we need some executable files from the cspice collection. You can't use what we have downloaded before. As we have seen, the executables are compiled with cgywin o aarch64 in mind. You may need a 64bit executable for windows of linux or macos, browse again to the jpl toolkit C sources page and download the tarball or the zip file that fits best your architecture. At present I'm using debian, so i'm going to download the PC, Linux, gCC, 64bit tarball cspice.tar.gz. If you are on Windows download the PC, Windows, MS Visual C, 64bit cspice.zip file. If your PC still has a 32bit architecture, choose the appropriate source.

If you unzip the cspice.tar.gz via tar or the cspice.zip via winzip or similar, you will see a replica of what you have already seen inside your smartphone. Change dir to cspice/exe, where you will find the already seen brief or brief.exe and another useful utility called spkmerge or spkmerge.exe. These are compiled for your PC architecture, so they'll surely work.

Copy brief(.exe) and spkmerge(.exe) to the directory where you have that giant 3gb file named de441.bsp

open a terminal in that directory and type:

./brief de441.bsp (in Linux)
brief.exe de441.bsp (in Windows)

You should see something like:

debian@debian11:~/Desktop/cspice$ ./brief de441.bsp
 
BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066
 
 
Summary for: de441.bsp
 
Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
        VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
        EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
        MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
        JUPITER BARYCENTER (5)  SUN (10)
        Start of Interval (ET)              End of Interval (ET)
        -----------------------------       -----------------------------
        13201 B.C. MAY 06 00:00:00.000      17191 MAR 15 00:00:00.000
 
debian11@debian11:~/Desktop/cspice$ 

the listed bodies are what serve our purpose. The time interval, as you can appreciate, is humongous, we don't need all that

Let's create an appropriate file to order the extraction of a subinterval, more adequate to our purposes.

Create a plain text file, let's call it extract.txt ad type inside of it:


   LEAPSECONDS_KERNEL  = naif0012.tls

   SPK_KERNEL          = from1501ADto2100AD.bsp
     SOURCE_SPK_KERNEL = de441.bsp
       BEGIN_TIME      = 01 JAN 1501 00:00:00.000
       END_TIME        = 31 DEC 2100 23:59:59.000

As you can see, I give the name of the leapsecond kernel (copy it to the same directory), the name of the destination bsp file (don't use dexxx, you are not allowed, according to the jpl policy) the source bsp file and the begin and end time required.

We're ready, type:

spkmerge extract.txt

In a fraction of time you will see a new file named from1501ADto2100AD.bsp, it's a subset of the big file. Using brief you can inspect it:

debian11@debian11:~/Desktop/cspice$ ./brief from1501ADto2100AD.bsp 
 
BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066
 
 
Summary for: from1501ADto2100AD.bsp
 
Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
        VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
        EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
        MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
        JUPITER BARYCENTER (5)  SUN (10)
        Start of Interval (ET)              End of Interval (ET)
        -----------------------------       -----------------------------
        1501 JAN 01 00:00:41.184            2101 JAN 01 00:01:08.183
 
debian11@debian11:~/Desktop/cspice$ 

It's rather large (64mb) but handier than the 3GB one. Using WinSCP o the the linux file manager, copy it to your phone home directory.

Creating a minimal C source file

I'm not an expert at C programming, but I'll try to write something simple and useful

Just a few preprocessor instructions, some functions and a main()

But first let's give a glance to our phone home directory. Access it using ssh from you PC, as we have seen before. It should be something like this:

~ $ ls
cspice      from1501ADto2100AD.bsp  mytrueepoch.tf  storage
cspice.tar  importCSpice.csh        naif0012.tls 	libcspice.so
~ $ 

Let's move our dependency files in a new subdirectory.

~ $ ls
clear   cspice.tar              importCSpice.csh  naif0012.tls
cspice  from1501ADto2100AD.bsp  mytrueepoch.tf    storage
~ $ mkdir build
~ $ mv from1501ADto2100AD.bsp naif0012.tls mytrueepoch.tf libcspice.so build/
~ $ ls
build  cspice  cspice.tar  importCSpice.csh  storage
~ $ cd build
~/build $ ls
from1501ADto2100AD.bsp  mytrueepoch.tf  libcspice.so	naif0012.tls
~/build $ 

Now we need a helper C file with some useful functions. This first example will only include two functions, one for converting radians to degrees (or hours) and decimal fractions, and another for converting the result in a more readeable form, with degrees, minutes and seconds.

We'll call the first r2d and the second ddd2dms. If you have never written a C script before, I'll give some useful hints

Let's begin writing two different files, one header and one source. We'll call the first tools.h and the second tools.c

tools.h

/*In most cases we want to avoid repeated callings to the header files.
  The following preprocessors instructions let you compile the header only once
  and only if not called before.
 */

#ifndef TOOLS_H
#define TOOLS_H

/* The following lines are preprocessor include instructions.
   The first one takes all the declarations contained in our cspice library and
   makes them available to the program.
   The second one calls the stdio component of the standard library. We use it
   mainly for printing the results of calculations.
   The third one supports mathematical calculations.
*/

#include "SpiceUsr.h"
#include <stdio.h>
#include <math.h>

/* Then we declare a structure, which is a complex set of data under one
   new variable type, to represent sign, degrees, minutes and seconds in every
   sexagesimal value.
*/

typedef struct _dms {
	int sign;
	int units;
	int minutes;
	int seconds;
} dms;

// Next we create a function which converts radians to degrees

double r2d (double rads);

// Then a functions which takes a decimal  and converts it to a sexagesimal.

dms ddd2dms(double dec_units);

#endif

tools.c

#include "tools.h"

double r2d (double rads){
		if (rads < 0){
			rads = rads + 2*M_PI;
		}
		return rads/M_PI*180.0;
}

dms ddd2dms(double dec_units){

	dms d2s;
	d2s.sign = '+';
    if (dec_units < 0){
        d2s.sign = '-';
        dec_units *= -1;
    }
    double total_seconds = ceil(dec_units * 3600.0);
    double seconds = fmod(total_seconds, 60.0);
    int total_minutes = (int)((total_seconds - seconds)/60.0);
    double minutes = fmod(total_minutes, 60.0);
    int units = (int)((total_minutes - minutes)/60.0);
    d2s.units = units;
    d2s.minutes = (int)minutes;
    d2s.seconds = (int)seconds;
    return d2s;
}

You can see that the header file explicitly recalls the SpiceUsr.h header of the cspice library. Let's assume that the compiler and the linker know where to find it (it's not true, but we'll get to it a few lines ahead). Apart from the SpiceUsr.h header, we call the stdio.h and math.h from the standard library. The first is used mainly for communications with the standard io, the second for calling mathematical methods. In the source file, tools.c, we simply call tools.h which automatically calls the other library components. As you can see, the header file contains only the declarations, the definitions are in the source file.

Now, just for testing, a little C source file, named test_tools.c

test_tools.c

#include "tools.h"

int main(){

	double radians = 12.3456;
	double decimals = r2d(radians);
	dms sexas;
	printf("radians: %f\n", radians);
	printf("decimals: %f\n", decimals);
	sexas = ddd2dms(decimals);
	printf("sexagesimal: %c %d %d %d\n", sexas.sign, \
		sexas.units, sexas.minutes, sexas.seconds);
	/* Let's call the spice builtin variable for pi and print it
	   just to be sure that the spice library is available
	   Give a glance to "https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/info/mostused.html#Y"
       to learn more about Spice constants
       The jpl toolkit signature for pi is <SpiceDouble pi_c ( void )>
    */
	printf("Pi value: %f\n", pi_c());
	return 0;
}

This source file contains the main routine, which defines a double variable with a sample value, converts it from radians to degrees, than converts the result to a string with sign, degrees, minutes and seconds as a sequence of integer values.

It may seem that the job is done. Actually, even if the cspice library is in the same directory, the SpiceUsr.h file is in another directory, we must tell the compiler where to find it.

Let's write a little script, containing what is needed (I could have written a Makefile, we'll do it in the future). Let's call it test_tools.sh where sh is the standard extension for bash files.

Before going on, let's download the package termux-elf-cleaner from the repository. It removes some boring warnings that come out when we execute the compiled file.

pkg install termux-elf-cleaner

Below a possible test_tools.sh file

test_tools.sh

gcc -L. -g -Wall -O2 -Wl,-rpath=. \
    -o test_tools test_tools.c tools.c \
    -lm -ldl -lcspice \
    -I/data/data/com.termux/files/home/cspice/include/

termux-elf-cleaner test_tools

Let's analyze this script, I have broken it in a few lines to make it easier to understand

  1. gcc is the calling to the compiler (and to the linker in this case)
  2. -L. tells the compiler where to find the shared library in use, in this case libcspice.so in the same directory
  3. -g tells the compiler to generate debugging infos
  4. -Wall for maximizing warnings
  5. -O2 for a very good optimization
  6. -Wl,-rpath=. to embed the location of the shared library in the executable itself, in this case the working directory. It's important for giving access to the shared library during runtime execution. Other methods are available.
  7. -lm -ldl -lcspice are links to math library, libdl is the dynamic link library, lcspice stands for load the libcpice.so library from where you know, lib and .so are conventionally omitted.
  8. The capital I is the location where lies the unfamous SpiceUsr.h header. In this case I prefer to give the full path. If you are wandering how you can obtain that long string, try typing in your command line:
    pwd
    which stays for "print the working directory", and you'll get the full path. Of course, you must modify it to obtain the correct location of the header.

Let's make this script executable using chmod test_tools.sh followed by Enter

Then we launch its execution: ./test_tools.sh (remember to begin with ./, it's not a standard shell command.

~/build $ ./test_tools.sh
termux-elf-cleaner: Replacing unsupported DF_1_* flags 134217729 with 1 in 'test_tools'
~/build $ 

Ok, now we can launch the test_tools executable (if it say Permission denied, it means you have to chmod +x it, as for the .sh file before):

./test_tools

This is what you get:

~/build $ ./test_tools
radians: 12.345600
decimals: 707.350776
sexagesimal: + 707 21 3
Pi value: 3.141593
~/build $ 

We have the demonstration that our new born functions work fine and that the SpiceUsr.h header is correctly addressed. Let me know in your comments. See you in the next post. Happy New Year to all of you.

How to create a virtual linux machine with qemu under Debian or Ubuntu with near native graphics performance

It's been a long time since my latest post. I know, I'm lazy. But every now and then I like to publish something that other people c...