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

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! 🙂