Переход на матричные преобразования

Следствия:

Gradient and height crosses - Исправление ошибки карты сдвига полюса

^{комментарии в LiveJournal} Ещё недавно косинусы углов считал функцией System.Math.Cos().

Углов стало больше, использую скалярное произведение единичных векторов. Крайние условия для тех, что коллинеарны, пришлось тщательно описывать, куда же денешься.

И тут понял, что матричные преобразования компактнее, учитывают крайние условия, косинусы связанных углов считают, и проще оптимизируются.

Вот код для матриц

var aOnSurface = axisOrtohonal.Direction * b.Matrix;
aTraverse = -a * aOnSurface[1];
var aMerid = (b.Vartheta < 0 ? 1 : -1) * a * aOnSurface[2];
return aMerid;

и для сравнения, как считалось скалярным произведением

var surfaceCalm = new Plane(b.NormalCalm, b.r);
var QonAxisPlane = new Plane(axisEnd, b.Q3, Basin.O3);

// aSphere direction
var aSphereLine = surfaceCalm.IntersectionWith(QonAxisPlane);

var b3unit = new Line3D(Basin.O3, b.Q3).Direction;
// lays in surfaceCalm plane, directed to equator of AxisOfRotation if Math.Abs used
var aSphere = //Math.Abs
    (a * AxisOfRotation.DotProduct(b3unit));// axisOrtohonal.Direction.DotProduct(aSphereLine.Direction)); 

var aMeridianLine = surfaceCalm.IntersectionWith(b.MeridianCalm); //new Plane(OzEnd, Q3, O3);
var aTraverseLine = surfaceCalm.IntersectionWith(b.TraverseCalm);
Assert.AreEqual(0, aMeridianLine.Direction.DotProduct(aTraverseLine.Direction), .000000001);

aTraverse = Math.Abs(aSphere * aSphereLine.Direction.DotProduct(aTraverseLine.Direction));

var planeOZ = new Plane(Basin.Oz);
var planeAxis = new Plane(AxisOfRotation);

//directed to equator of Oz if Math.Abs used
double aMeridian;
/*if (aSphereLine.IsCollinear(aMeridianLine, .1))
{
    aMeridian = aSphere;
}
else*/
{
    var dotProduct = aSphereLine.Direction.DotProduct(aMeridianLine.Direction);
    aMeridian = Math.Abs(aSphere * dotProduct);
}


// if (AxisOfRotation != Basin.Oz)
var spin = new Plane(Basin.OzEnd, b3unit.ToPoint3D(), axisEnd).Normal.DotProduct(b3unit);

// "north" hemisphere of AxisOfRotation
if (b3unit.DotProduct(AxisOfRotation) > 0) 
{
    if (spin > 0)
    {
        aTraverse = -aTraverse;
    }
}
else
{
    if (spin < 0)
    {
        aTraverse = -aTraverse;
    }
}

//aMeridian<0 if Q3 between planes or inside of cones (bug when angle is near 90) of OZ and AxisOfRotation
if (planeOZ.SignedDistanceTo(b.Q3) * planeAxis.SignedDistanceTo(b.Q3) < 0)
{
    aMeridian = -aMeridian;
}
else
{
    var coneAxis = new UnitVector3D((Basin.Oz + AxisOfRotation).ToVector());
    if (new UnitVector3D(b.Q3.ToVector()).DotProduct(coneAxis) > coneAxis.DotProduct(Basin.Oz))
    {
        //inside cone
        aMeridian = -aMeridian;
    }
}
return aMeridian;

диаграмма векторов

результаты

^{комментарии в LiveJournal} одолел-таки скалярные произведения векторов

в модели океана без дна сдвиг полюса на 17 градусов (из текущего положение в 40 зап.долг. 73 сев.шир.) вызывает волну относительно текущего геоида высотой около 5 км

полюс показан розовыми точками

этот пост будет изменяться...