Arc points calculation in spherical coordinates for DMX Moving Head

From time to time in the process of developing software engineers come across various mathematical issues. Usually those issues are easy fixes and can be resolved applying elementary mathematics, but sometimes we are facing complex and intriguing task. Please, see below a case study, I hope you find it interesting.

While helping to develop Gunnebo Business Solution product, our customer Bjorn Nostdahl requested to create custom control system to be used for the set of different DMX Lights. You know, those kinds of lights that are being used for huge concerts, TV shows or night clubs.

One of those lights is DMX Moving Head, it rotates almost to any angle, changes focus, colors, intensity, etc. The customer requested the sun movement simulation on a ceiling, starting from sunrise till sunset to impress conference participants.

Imagine, we would need to calculate light path from point A to point B along the arc. Whereas light can be rotated by 360 degrees vertically and horizontally, anyway we should use spherical coordinate system for points forming the arc.

To start with, let’s try to recall spherical coordinates formulas

Where r > 0, 0 ≤ θ ≤ 2π , 0 ≤ φ ≤ π

Pic.1. Spherical coordinate system

In our case there are three points in spherical coordinates: starting point A, target point B and center point C. All these points belong to the sphere. We don’t care about sphere radius and can use unit sphere for calculation purposes because ray of light does not have physical distance.

The first step is converting all points into Cartesian coordinate system. We’ll use formula 1 for this purpose. JavaScript implementation for Cartesian conversion is the following:

// spherical to cartesian coordinates
// phi is in [0,2*pi]
// theta is in [0,pi]
function toCartesian(point) {
   var _point = {};
   var r = 10;
   _point.x = r * Math.sin(point.theta) * Math.cos(point.phi);
   _point.y = r * Math.sin(point.theta) * Math.sin(point.phi);
   _point.z = r * Math.cos(point.theta);
   return _point;
}

Ac = toCartesian(A);
Bc = toCartesian(B);

Next we would need to find normal vector for a plane, where points A and B belong. On this stage normal vector will be the equal unit vector vn, which will cross circle center point C. We can calculate normal vector the following way:

function moduleVector(v){
   return Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

function normalizeVector(v){
   var _vectorModule = moduleVector(v);
   return {x:v.x/_vectorModule,y:v.y/_vectorModule,z:v.z/_vectorModule};
}

vn = normalizeVector({x:Cc.x, y:Cc.y, z:Cc.z});

Also we need to find the real center point Cc’, which will be a cross point of normal vector vn and the plane where points Ac and Bc belong.

function vectorScalar(v1, v2){
   return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
// Find center point C0 what is on vector OC and on ABC plane
// A,B - points on the sphere
// C - projection of the center point on the sphere
function calcCenterPoint(centerOnSphere, p) {
  var moduleSqr=vectorScalar(centerOnSphere,centerOnSphere);
  var сoef=1-(moduleSqr-vectorScalar(centerOnSphere,p))/moduleSqr;
  return {
       x: centerOnSphere.x * coef,
       y: centerOnSphere.y * coef,
       z: centerOnSphere.z * coef};
}
Cc' = calcCenterPoint(Cc, Ac);

Finally, there are three points describing the plane. We’ll save Cc’ to movement_point  variable and will subtract point Cc’ from every point we have to simplify calculation. We need movement_point point for reverse transformation. After such transformations the point  Cc’ will be positioned in coordinate system center O(0,0,0). And the most interesting point – rotating coordinate system to make AcBcCc’ plane equal Oxy plane. We will use the following formulas (2) for that purpose:

Please, find JavaScript implementation below:

const OZ = {x:0, y:0, z:-1};
const epsilon = Math.exp(-16);
function getVectorMultiplication(p1p2Vect, p1p3Vect){
   var i_v=p1p2Vect.y*p1p3Vect.z-p1p2Vect.z*p1p3Vect.y;
   var j_v=-(p1p2Vect.x*p1p3Vect.z-p1p2Vect.z*p1p3Vect.x);
   var k_v = p1p2Vect.x*p1p3Vect.y-p1p2Vect.y*p1p3Vect.x;
   return {x:i_v, y:j_v, z:k_v};
}

function vectorsAngle(v1, v2){
   return Math.acos(vectorScalar(v1,v2)/(moduleVector(v1)*moduleVector(v2)));
}

function rotatePoint_xy_quarterion(p, normalVector){
   var angle = vectorsAngle(OZ,normalVector);
   // find case when angle near 0 or PI
   var an = (angle + 2*Math.PI)%(Math.PI);
   if(an < epsilon || an > Math.PI - epsilon){
       return p;
   }
   var b = getVectorMultiplication(OZ,normalVector);
   b = normalizeVector(b);

   var q0 = Math.cos(angle/2);
   var q1 = Math.sin(angle/2)*b.x;
   var q2 = Math.sin(angle/2)*b.y;
   var q3 = Math.sin(angle/2)*b.z;

   var rm = [
       [q0*q0+q1*q1-q2*q2-q3*q3,2*(q1*q2-q0*q3),2*(q1*q3+q0*q2)],
       [2*(q1*q2+q0*q3),q0*q0-q1*q1+q2*q2-q3*q3,2*(q2*q3-q0*q1)],
       [2*(q1*q3-q0*q2),2*(q2*q3+q0*q1),q0*q0-q1*q1-q2*q2+q3*q3]];
   var res = {};
   res.x = p.x*rm[0][0] + p.y*rm[1][0] + p.z*rm[2][0];
   res.y = p.x*rm[0][1] + p.y*rm[1][1] + p.z*rm[2][1];
   res.z = p.x*rm[0][2] + p.y*rm[1][2] + p.z*rm[2][2];

   return res;
}
Acr= rotatePoint_xy_quarterion(Ac,vn);
Bcr = rotatePoint_xy_quarterion(Bc ,vn);
Ccr = rotatePoint_xy_quarterion(Cc’,vn);
backVector = rotatePoint_xy_quarterion(OZ,vn);

At this stage we can ignore z coordinate and use 2D formulas because point Сс’ is in the coordinate system center and points Ac, Bc belong to Oxy plane. Z coordinate equals zero. Also we calculated backVector. It’s necessary for coordinate system reverse rotation.

As a next step we would need to calculate t-parameter value for points Ac and Bc:

function pointsDistanceInCartesian(point1, point2){
   var xDelta = point1.x - point2.x;
   var yDelta = point1.y - point2.y;
   var zDelta = point1.z - point2.z;
   return Math.sqrt(xDelta*xDelta+yDelta*yDelta+zDelta*zDelta);
}
var radius = pointsDistanceInCartesian(Acr,Ccr);

var Acr_t = (Math.acos(Acr.x/radius)+2*Math.PI)%(2*Math.PI);
var Bcr_t = (Math.acos(Bcr.x/radius)+2*Math.PI)%(2*Math.PI);

Now, as both sector angle created by points Acr, Bcr, Ccr  and angle step are defined, we could easily calculate circle points coordinates. We’ll find sector angle and iteration step.

var iterations = 100;
var sectAngle = Bcr_t - Acr_t;
var angleDelta=Math.abs(sectAngle)<=Math.PI?sectAngle:-(sectAngle-Math.PI);
var angleStep = angleDelta/iterations;

The final step is to find point coordinates for each iteration and to make reverse transformation.

// spherical coordinates => polar coordinates
// circle plane is parallel with Oxy plane => z=0; theta = pi/2
// sin(theta) = sin(pi/2) = 1;
// cos(theta) = cos(pi/2) = 0;
// x = r*sin(theta)*cos(phi) = r*cos(phi)
// y = r*sin(theta)*sin(phi) = r*sin(phi)
// z = r*cos(theta) = 0;
function getCirclePoint(t, radius){
   return {x:radius * Math.cos(t), y:radius*Math.sin(t), z:0};
}
// convert cartesian to spherical coordinates
function toSpherical(p) {
   var _p = {};
   _p.r = Math.sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
   if(p.x == 0 && p.y == 0){
       _p.theta = p.z >= 0 ? 0 : Math.PI;
       _p.phi = 0;
   } else {
       _p.theta = Math.acos(p.z / _p.r);
       _p.phi = (Math.atan2(p.y, p.x) + 2*Math.PI)%(2*Math.PI);
   }
   return _p;
}
var timeoutID;
for (var i = 0; i <= iterations; i++) {
   // calculate next t value
   var t = Acr_t + i*angleStep;
   // calculate point coordinates
   var iterPoint = getCirclePoint(t,radius);
   // rotate coordinaеe system back 
   iterPoint = rotatePoint_xy_quarterion(iterPoint, backVector);
   // move points back
   iterPoint = movePointInCartesian(iterPoint, movement_point,1);
  
   var p = toSpherical({x:iterPoint.x,y:iterPoint.y,z:iterPoint.z});
}

That’s it! In this article I outlined how to find Arc point coordinates for circle on the sphere. I hope you will find it helpful one day! 🙂