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

Следствия:

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 км.

Новый полюс показан светло-розовым пятном .

В начале расчета есть искажение перетекания воды, которое исчезает в конце.

Это два оранжево-красные поднятиями воды из севера на увеличенной картинке: .

Причина - в резком сдвиге полюса. Искажение исправляется при плавном смещении.

другие углы сдвига

Сдвиг полюса на 45 градусов (из текущего положение в 40 зап.долг. 45 сев.шир.) вызывает волну высотой около 14 км:

Сдвиг на 90 градусов вызывает волну высотой около 21 км:

одинаковая сила притяжения

Если сдвинуть полюс на 17 градусов относительно текущего геоида (приплюснутой сферы) и предположить, что сила притяжения на экваторе такая же, как на полюсах, то получается волна высотой около 13 км:

Расстояния к центру Земли (минус [math]R[/math] - Эллипсоид, радиус равновеликого шара) будут такими: .