Trying to debug and verify port of c/c++ code

I’m working on some basic celestial navigation concepts, using Processing and Android. This particular problem is in Processing and is an attempt to port c/c++ code to Processing/Java.

I started with reading a paper from this link: VectorSolution

His example says the Altitude should be 22.041 and the Azimuth should be 138.2, for the assumed position and geographic position I used as input in the code I post below.

Something just isn’t right. I have a vague idea that C uses pointers and Java doesn’t worry about it until it causes an exception.

Here’s my code:

import controlP5.*;

//From NavigationalAlgorithms Vectorsolution 

/*
File: vector.cpp
 Cálculo vectorial
 Resultado test: OK
 This file contains proprietary information of Andrés Ruiz Gonzalez
 Andrés Ruiz. San Sebastian - Donostia. Gipuzkoa
 Copyright (c) 1999 - 2007
 */

double Hc;
double Z;

double GHA = 347.78;
double dec = -16.72;
double Be = 40;
double Le = 28;

//ControlP5 cp5;
//ControlP5 controlP5;
//Textarea myTextarea;
void setup() {
  //  fullscreen();

  size(127, 75, OPENGL);

  /*
  myTextarea = cp5.addTextarea("txt")
    .setPosition(450,10)
    .setSize(800,700)
    .setFont(createFont("arial",14))
    .setLineHeight(14)
    .setColor(color(255))
    .setColorBackground(color(255,0))
    .setColorForeground(color(255,100));
  */
}
void draw() {
  AltitudeAzimuth( GHA, dec, Be, Le);
  println("Altitude:  ", Hc);
  println(" Azimuth:  ", Z);
}

double Mod(double x[]) {
  return(Math.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]));
}

double[] Add(double x[], double y[]) {
  double z[] = new double[3];
  for (int i = 0; i < 3; i++) {
    z[i] = x[i] + y[i];
  }
  return   z ;
}

double[] aVector(double a, double x[]) {
  double z[] = new double[3];
  z[0] = a * x[0];
  z[1] = a * x[1];
  z[2] = a * x[2];
  return (z);
}

double[] Unit(double x[]) {
  return (aVector(1.0/Mod(x), x));
}

double Dot(double x[], double y[]) {
  return (x[0]*y[0]+x[1]*y[1]+x[2]*y[2]);
}

double[] Cross(double x[], double y[]) {
  double z[] = new double[3];
  z[0] = x[1]*y[2] - z[2]*y[1];
  z[1] = x[2]*y[0] - x[0]*y[2];
  z[2] = x[0]*y[1] - x[1]*y[0];
  return (z);
}

double[] VectorSpherical2Cartesian(double B, double L) {
  double v[] = new double[3];
  v[0] = Math.cos(B) * Math.cos(L);
  v[1] = Math.cos(B) * Math.sin(L);
  v[2] = Math.sin(B);
  println(B);
  println(L);
  return(v);
}

void AltitudeAzimuth(double GHA, double dec, double Be, double Le) {
  double GP[];
  double AP[];
  double Np[] = {0, 0, 1};
  double U1[], U2[], U3[], Az[];

  GP = VectorSpherical2Cartesian(dec, GHA);

  AP = VectorSpherical2Cartesian(Be, Le);
  println("AP  : ", AP[0], " ", AP[1], " ", AP[2]);

  Hc = 90.0 - Math.acos(Dot(AP, GP));
  println((Hc));

  U1 = Unit(Cross(AP, Np));
  println("U1  : ", U1[0], " ", U1[1], " ", U1[2]);

  U2 = Cross(U1, AP);
  println("U2  : ", U2[0], " ", U2[1], " ", U2[2]);

  U3 = Unit(Cross(AP, GP));
  println("U3  : ", U3[0], " ", U3[1], " ", U3[2]);

  Az = Cross(U3, AP);
  println("Az  : ", Az[0], " ", Az[1], " ", Az[2]);

  Z = Math.acos(Dot(U2, Az));
  //   println(Z);
  if (Dot(U1, Az) < 0.0) {
    Z = 360.0 - Z;
  }
  //   println(Z);
}
1 Like

The problem is I should get
Altitude: 22.041
Azimuth: 138.2

based on the inputs I used which came from the paper I linked.

What I get after my attempt to port the C/C++ code to Java/Processing
Altitude: 89.573
Azimuth: 0.702119

/**
 File: vector.cpp
 Cálculo vectorial
 Resultado test: OK
 This file contains proprietary information of Andrés Ruiz Gonzalez
 Andrés Ruiz. San Sebastian - Donostia. Gipuzkoa
 Copyright (c) 1999 - 2007
 */

#include <math.h>;

double Mod( double *x )
{
  return( sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]) );
}

double* Add( double *x, double *y )
{
  double *z = new double[3];

  z[0] = x[0]+y[0];
  z[1] = x[1]+y[1];
  z[2] = x[2]+y[2];

  return( z );
}

double* aVector( double a, double *x )
{
  double *z = new double[3];
  z[0] = a*x[0];
  z[1] = a*x[1];
  z[2] = a*x[2];

  return( z );
}

double* Unit( double *x )
{
  return( aVector( 1.0/Mod(x), x ) );
}

double Dot( double *x, double *y )
{
  return( x[0]*y[0]+x[1]*y[1]+x[2]*y[2] );
}

double* Cross( double *x, double *y )
{
  double *z = new double[3];

  z[0] = x[1]*y[2]-x[2]*y[1];
  z[1] = x[2]*y[0]-x[0]*y[2];
  z[2] = x[0]*y[1]-x[1]*y[0];

  return( z );
}

double* VectorSpherical2Cartesian( double B, double L )
{
  double *v = new double[3];
  // Es un vector unitario
  v[0] = COS( B )*COS( L );
  v[1] = COS( B )*SIN( L );
  v[2] = SIN( B );

  return( v );
}

void AltitudeAzimuth( double GHA, double dec, double Be, double Le, double* Hc, double* Z )
{
  double *GP; // Geographic Position of the celestial body
  double *AP; // Assumed Position
  double Np[3] = {0, 0, 1}; // North pole

  double *U1, *U2, *U3, *Az;

  GP = VectorSpherical2Cartesian( dec, GHA );
  AP = VectorSpherical2Cartesian( Be, Le );

  *Hc = 90.0 - ACOS( Dot( AP, GP ) );

  U1 = Unit( Cross( AP, Np ) );
  U2 = Cross( U1, AP );
  U3 = Unit( Cross( AP, GP ) );
  Az = Cross( U3, AP );

  *Z = ACOS( Dot( U2, Az ) );

  if ( Dot( U1, Az ) < 0.0 ) *Z = 360.0 - *Z;

  delete[] GP;
  delete[] AP;
  delete[] U1;
  delete[] U2;
  delete[] U3;
  delete[] Az;
}

double LunarDistance( double decM, double ghaM, double decB, double ghaB )
{
  double *Vm, *Vb;
  double LD;
  Vm = VectorSpherical2Cartesian( decM, ghaM );
  Vb = VectorSpherical2Cartesian( decB, ghaB );
  LD = ACOS( Dot( Vm, Vb ) );
  delete[] Vm;
  delete[] Vb;
  return( LD );
}

inline double Star2StarDistance( double dec1, double gha1, double dec2, double gha2 )
{
  return( LunarDistance( dec1, gha1, dec2, gha2 ) );
}

inline double GC_Distance( double B1, double L1, double B2, double L2 )
{
  return( 60.0* LunarDistance( dec1, gha1, dec2, gha2 ) ); // nautical miles
}

// GC Waypoints (Bx, Lx)

double B_GC( double B1, double L1, double B2, double L2, double Lx )
{
  double Bx;

  double *V1, *V2;
  double *V1xV2;

  V1 = VectorSpherical2Cartesian( B1, L1 );
  V2 = VectorSpherical2Cartesian( B2, L2 );

  V1xV2 = Cross( V1, V2 );

  Bx = -ATAN( V1xV2[0]/V1xV2[2]*COS( Lx ) + V1xV2[1]/V1xV2[2]*SIN( Lx ) );

  delete[] V1;
  delete[] V2;
  delete[] V1xV2;

  return( Bx );
}

Yeah, that’s the complete cpp code I tried to port and the cpp code from the link to the paper. There are a whole lot of *'s and I don’t know if it translates the way I tried. Thanks for posting the relevant code.

In Java, any datatype apart the 8 primitive 1s: :coffee:

is a pointer/reference type already. :wink:

3 Likes

Thanks, pointers can be eliminated as the problem. I think it has something to do with radians vs. degrees, so I added this to my initial variables:

  double GHA = Math.toRadians(347.78);
  double dec = Math.toRadians(-16.72);
  double Be = Math.toRadians(40);
  double Le = Math.toRadians(28);

It still doesn’t give me good result.

Notice Processing has a class called vector, and you can use it to store your triplets.

What is GHA and dec? Can you explain about the conversion? What does it do? You need to focus in this function: VectorSpherical2Cartesian() as it seems this is doing all the heavy work. Ok, there is also AltitudeAzimuth() which is calling the Dot() function at the end. The implementation of the Dot function there doesn’t seen to be the standard dot product used in vector algebra.

Kf

OK, think of GHA as the longitude of a star/celestial body and dec as the latitude. GHA stands for Greenwich Hourly Angle, as the earth rotates on its axis, the stars appear to rotate. If the star is projected directly onto the earth, the GHA is the longitude of where it will land. The declination(dec) is the latitude.

The conversion from Spherical to Cartesian, I think, is necessary because linear algebra works in planes.

What I’ve done is start with GHA,dec,Be,Le all in radians. VectorSpherical2Cartesian uses them and converts a x,y point into a x,y,z point using spherical trig formulas. It needs radians to use Math.cos() and Math.sin().

In the AltitudeAzimuth(), I changed Hc and Z to:

     Hc = (90.0 - Math.toDegrees(Math.acos(Dot(AP,GP))));

    Z = Math.toDegrees(Math.acos(Dot(U2,Az)));

When I check my results now, the altitude is exactly as in the paper, but azimuth is all over the place.

As for the Dot function not being the standard dot product, the wiki definition starts off with “algebraic definition” and that is how it is implemented in this code. It is a summation of the products, a1b1 + a2b2 + a3*b3. I have only seen this in my limited experience.

I still cant figure out why my version of code doesn’t give the same azimuth. One set of inputs, the azimuth is off by 4 degrees, another is off by 30 degrees, and another is only off by 2/10 degrees. The altitude agrees exactly.

I found the bug. The Dot product is fine, but the Cross product has a typo. Gots to check the simple stuff, stupid.

2 Likes

I’ve moved on to the next paper from this link: Papers on navigation

The Matrix paper is tougher to implement without an explicit code example but I’ve done it. It’s very good. I just used the formulas he provided in a function and return the CartesiantoSpherical coordinates. This time I got exactly what he wrote without hunting down typos/bugs.

One more comment on pointers/reference. I think I see it now how Java does it, after reading error messages over and over about uninitialized variables or no such variable. If I want some function to access a variable outside itself, then I can initialize the variable before setup(Processing) and that will “point” the function to that variable if needed.

I wont worry about all those *'s ever again.