/**
 *
 *  脱离Cesium球体进行各种坐标转换
 *
 *     有个人理解修改，注释说明个人理解
 *
 *             金磊
 *
 */

/**
 *  myScene格式及赋值
 *   myScene = {
 *          width:scene.canvas.clientWidth,
 *          height:scene.canvas.clientHeight,
 *          viewMatrix:scene.camera.viewMatrix,//([16])
 *          projectionMatrix:scene.camera.frustum.projectionMatrix,//([16])
 *          frustum:{
 *                   fovy:viewer.scene.camera.frustum.fovy,
 *                   aspectRatio:viewer.scene.camera.frustum.aspectRatio,
 *                   near:viewer.scene.camera.frustum.near,
 *                  },
 *          positionWC:viewer.scene.camera.positionWC,
 *          directionWC:viewer.scene.camera.directionWC,
 *          rightWC:viewer.scene.camera.rightWC,
 *          upWC:viewer.scene.camera.upWC
 *         }
 */

//球体参数(4490坐标系)

var Ellipsoid = {
    radii: {
        //椭球体半径
        x: 6378137.0,
        y: 6378137.0,
        z: 6356752.31414035585,
    },
    radiiSquared: {
        //椭球体半径平方
        x: 40680631590769.0,
        y: 40680631590769.0,
        z: 40408299983328.77,
    },
    oneOverRadii: {
        x: 1.0 / 6378137.0,
        y: 1.0 / 6378137.0,
        z: 1.0 / 6356752.31414035585,
    },
    oneOverRadiiSquared: {
        x: 1.0 / 40680631590769.0,
        y: 1.0 / 40680631590769.0,
        z: 1.0 / 40408299983328.77,
    },
    centerToleranceSquared: 0.1,
};
var EPSILON12 = 0.000000000001;
//屏幕坐标转经纬度
function rayEllipsoid(ray) {
    var inverseRadii = Ellipsoid.oneOverRadii;
    var q = multiplyComponents(inverseRadii, ray.origin);
    var w = multiplyComponents(inverseRadii, ray.direction);

    var q2 = q.x * q.x + q.y * q.y + q.z * q.z;
    var qw = cartesianDot(q, w);

    var difference, w2, product, discriminant, temp;

    if (q2 > 1.0) {
        // Outside ellipsoid.
        if (qw >= 0.0) {
            // Looking outward or tangent (0 intersections).
            return undefined;
        }

        // qw < 0.0.
        var qw2 = qw * qw;
        difference = q2 - 1.0; // Positively valued.
        w2 = w.x * w.x + w.y * w.y + w.z * w.z;
        product = w2 * difference;

        if (qw2 < product) {
            // Imaginary roots (0 intersections).
            return undefined;
        } else if (qw2 > product) {
            // Distinct roots (2 intersections).
            discriminant = qw * qw - product;
            temp = -qw + Math.sqrt(discriminant); // Avoid cancellation.
            var root0 = temp / w2;
            var root1 = difference / temp;
            if (root0 < root1) {
                return { start: root0, stop: root1 };
            }

            return {
                start: root1,
                stop: root0,
            };
        }
        // qw2 == product.  Repeated roots (2 intersections).
        var root = Math.sqrt(difference / w2);
        return { start: root, stop: root };
    } else if (q2 < 1.0) {
        // Inside ellipsoid (2 intersections).
        difference = q2 - 1.0; // Negatively valued.
        w2 = w.x * w.x + w.y * w.y + w.z * w.z;
        product = w2 * difference; // Negatively valued.

        discriminant = qw * qw - product;
        temp = -qw + Math.sqrt(discriminant); // Positively valued.
        return { start: 0.0, stop: temp / w2 };
    }
    // q2 == 1.0. On ellipsoid.
    if (qw < 0.0) {
        // Looking inward.
        w2 = w.x * w.x + w.y * w.y + w.z * w.z;
        return { start: 0.0, stop: -qw / w2 };
    }

    // qw >= 0.0.  Looking outward or tangent.
    return undefined;
}
//笛卡尔坐标相加
function cartesianAdd(left, right) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };
    result.x = left.x + right.x;
    result.y = left.y + right.y;
    result.z = left.z + right.z;
    return result;
}
//笛卡尔坐标相减
function subtract(left, right) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };
    result.x = left.x - right.x;
    result.y = left.y - right.y;
    result.z = left.z - right.z;
    return result;
}
function magnitude(cartesian) {
    return Math.sqrt(
        cartesian.x * cartesian.x +
            cartesian.y * cartesian.y +
            cartesian.z * cartesian.z
    );
}
//获得webgl透视图
function getPickRayPerspective(scene, windowPosition) {
    var result = {
        origin: { x: 0.0, y: 0.0, z: 0.0 },
        direction: { x: 0.0, y: 0.0, z: 0.0 },
    };

    var width = scene.width;
    var height = scene.height;

    var tanPhi = Math.tan(scene.frustum.fovy * 0.5);

    var tanTheta = scene.frustum.aspectRatio * tanPhi;
    var near = scene.frustum.near;

    //将坐标缩放到-1~1内
    var x = (2.0 / width) * windowPosition.x - 1.0;
    var y = (2.0 / height) * (height - windowPosition.y) - 1.0;

    var position = scene.positionWC;
    result.origin.x = position.x;
    result.origin.y = position.y;
    result.origin.z = position.z;
    //获取near中心点相对相机的位置
    var nearCenter = multiplyByScalar(scene.directionWC, near);
    //获取near中心点在webgl场景绝对位置
    nearCenter = cartesianAdd(position, nearCenter);

    var xDir = multiplyByScalar(scene.rightWC, x * near * tanTheta);
    var yDir = multiplyByScalar(scene.upWC, y * near * tanPhi);
    //计算传入射线与near平面交点在webgl场景内的绝对位置
    result.direction = cartesianAdd(nearCenter, xDir);
    result.direction = cartesianAdd(result.direction, yDir);
    //计算传入射线与near平面交点与相机的相对位置
    result.direction = subtract(result.direction, position);
    //除以与相机连线长度,计算出三个角度的tan值
    result.direction = normalize(result.direction);

    return result;
}
function getPoint(ray, t) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };
    result = multiplyByScalar(ray.direction, t);
    result = cartesianAdd(ray.origin, result);
    return result;
}
/**
 *
 *  scene = {
 *          width:scene.canvas.clientWidth,
 *          height:scene.canvas.clientHeight,
 *          frustum:{
 *                   fovy:viewer.scene.camera.frustum.fovy,
 *                   aspectRatio:viewer.scene.camera.frustum.aspectRatio,
 *                   near:viewer.scene.camera.frustum.near,
 *                  },
 *          positionWC:viewer.scene.camera.positionWC,
 *          directionWC:viewer.scene.camera.directionWC,
 *          rightWC:viewer.scene.camera.rightWC,
 *          upWC:viewer.scene.camera.upWC
 *          }
 *
 *
 */
function pickEllipsoid(scene, windowPosition) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    //获得获取屏幕像素在webgl内射线的透视图。
    var ray = getPickRayPerspective(scene, windowPosition);
    //计算射线与球的交点
    var intersection = rayEllipsoid(ray);
    if (!intersection) {
        return undefined;
    }
    var t = intersection.start > 0.0 ? intersection.start : intersection.stop;
    result = getPoint(ray, t);
    return result;
}
function scaleToGeodeticSurface(cartesian) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };

    var oneOverRadii = Ellipsoid.oneOverRadii;
    var oneOverRadiiSquared = Ellipsoid.oneOverRadiiSquared;
    var centerToleranceSquared = Ellipsoid._centerToleranceSquared;

    var positionX = cartesian.x;
    var positionY = cartesian.y;
    var positionZ = cartesian.z;

    var oneOverRadiiX = oneOverRadii.x;
    var oneOverRadiiY = oneOverRadii.y;
    var oneOverRadiiZ = oneOverRadii.z;

    var x2 = positionX * positionX * oneOverRadiiX * oneOverRadiiX;
    var y2 = positionY * positionY * oneOverRadiiY * oneOverRadiiY;
    var z2 = positionZ * positionZ * oneOverRadiiZ * oneOverRadiiZ;

    // Compute the squared ellipsoid norm.
    var squaredNorm = x2 + y2 + z2;
    var ratio = Math.sqrt(1.0 / squaredNorm);

    // As an initial approximation, assume that the radial intersection is the projection point.
    var intersection = multiplyByScalar(cartesian, ratio);

    // If the position is near the center, the iteration will not converge.
    if (squaredNorm < centerToleranceSquared) {
        if (!isFinite(ratio)) {
            return undefined;
        } else {
            result.x = intersection.x;
            result.y = intersection.y;
            result.z = intersection.z;
            return result;
        }
    }

    var oneOverRadiiSquaredX = oneOverRadiiSquared.x;
    var oneOverRadiiSquaredY = oneOverRadiiSquared.y;
    var oneOverRadiiSquaredZ = oneOverRadiiSquared.z;

    // Use the gradient at the intersection point in place of the true unit normal.
    // The difference in magnitude will be absorbed in the multiplier.
    var gradient = { x: 0.0, y: 0.0, z: 0.0 };
    gradient.x = intersection.x * oneOverRadiiSquaredX * 2.0;
    gradient.y = intersection.y * oneOverRadiiSquaredY * 2.0;
    gradient.z = intersection.z * oneOverRadiiSquaredZ * 2.0;

    // Compute the initial guess at the normal vector multiplier, lambda.
    var lambda =
        ((1.0 - ratio) * magnitude(cartesian)) / (0.5 * magnitude(gradient));
    var correction = 0.0;

    var func;
    var denominator;
    var xMultiplier;
    var yMultiplier;
    var zMultiplier;
    var xMultiplier2;
    var yMultiplier2;
    var zMultiplier2;
    var xMultiplier3;
    var yMultiplier3;
    var zMultiplier3;

    do {
        lambda -= correction;

        xMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredX);
        yMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredY);
        zMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredZ);

        xMultiplier2 = xMultiplier * xMultiplier;
        yMultiplier2 = yMultiplier * yMultiplier;
        zMultiplier2 = zMultiplier * zMultiplier;

        xMultiplier3 = xMultiplier2 * xMultiplier;
        yMultiplier3 = yMultiplier2 * yMultiplier;
        zMultiplier3 = zMultiplier2 * zMultiplier;

        func = x2 * xMultiplier2 + y2 * yMultiplier2 + z2 * zMultiplier2 - 1.0;

        // "denominator" here refers to the use of this expression in the velocity and acceleration
        // computations in the sections to follow.
        denominator =
            x2 * xMultiplier3 * oneOverRadiiSquaredX +
            y2 * yMultiplier3 * oneOverRadiiSquaredY +
            z2 * zMultiplier3 * oneOverRadiiSquaredZ;

        var derivative = -2.0 * denominator;

        correction = func / derivative;
    } while (Math.abs(func) > EPSILON12);

    result.x = positionX * xMultiplier;
    result.y = positionY * yMultiplier;
    result.z = positionZ * zMultiplier;
    return result;
}
function cartesianToDegrees(cartesian) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };
    var oneOverRadiiSquared = Ellipsoid.oneOverRadiiSquared;
    var p = scaleToGeodeticSurface(cartesian);

    if (p == undefined) {
        return undefined;
    }

    var n = multiplyComponents(p,oneOverRadiiSquared);
    n = normalize(n);

    var longitude = Math.atan2(n.y, n.x);
    var latitude = Math.asin(n.z);

    result.x = longitude;
    result.y = latitude;
    return result;
}
/**
 * 根据屏幕像素取经纬度
 * @param {*} x
 * @param {*} y
 */
function tempgetCoord(x, y) {
    //像素位置转迪卡尔坐标
    var car3_lc = pickEllipsoid(myScene, { x: x, y: y });
    if (car3_lc) {
        //迪卡尔坐标转弧度
        var cartographicLC = cartesianToDegrees(car3_lc);
        //转经纬度
        var point = toDegrees([
            cartographicLC.x,
            cartographicLC.y,
        ]);
        return [point.x, point.y];
    } else {
        return null;
    }
}

//经纬度转屏幕坐标

//弧度/度
var RADIANS_PER_DEGREE = Math.PI / 180.0;
//度/弧度
var DEGREES_PER_RADIAN = 180.0 / Math.PI;
//弧度转笛卡尔坐标与半径的比例
/**
 *      先计算cos(lat),结果没意义为了化简
 *      P/R = cos(lat);
 *      cos(lat) * cos(lon) = P/R * x/P = x/R
 *      cos(lat) * sin(lon) = P/R * y/P = y/R
 *      sin(lat) = z/R
 */
function geodeticSurfaceNormalCartographic(cartographic) {
    var longitude = cartographic.x;
    var latitude = cartographic.y;

    var cosLatitude = Math.cos(latitude);
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = cosLatitude * Math.cos(longitude);
    result.y = cosLatitude * Math.sin(longitude);
    result.z = Math.sin(latitude);
    return result;
}
//球体笛卡尔与坐标笛卡尔乘积
function multiplyComponents(left, right) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = left.x * right.x;
    result.y = left.y * right.y;
    result.z = left.z * right.z;

    return result;
}
function toDegrees(point) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = point[0] * DEGREES_PER_RADIAN;
    result.y = point[1] * DEGREES_PER_RADIAN;
    return result;
}
//经纬度转弧度
function toRadians(point) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = point[0] * RADIANS_PER_DEGREE;
    result.y = point[1] * RADIANS_PER_DEGREE;
    return result;
}
//计算两个笛卡尔坐标的（标量）乘积
function cartesianDot(left, right) {
    return left.x * right.x + left.y * right.y + left.z * right.z;
}
function divideByScalar(cartesian, scalar) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = cartesian.x / scalar;
    result.y = cartesian.y / scalar;
    result.z = cartesian.z / scalar;
    return result;
}
//矢量成标量
function multiplyByScalar(cartesian, scalar) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    result.x = cartesian.x * scalar;
    result.y = cartesian.y * scalar;
    result.z = cartesian.z * scalar;
    return result;
}
//经纬度转笛卡尔坐标
/**
 *      有误差，和Cesium计算方法一致，误差都一样相对位置应该没什么问题
 *      先使用了“Haversine公式”方法，最大误差应该为0.4%，
 *      然后使用有误差的值，去纠正另一个有误差的值，感觉这个思路的逻辑不太好，但应该会缩小误差。
 * */
function degreesToCartesian(point) {
    //转换坐标为弧度坐标
    var cartographic = toRadians(point);
    //弧度坐标转换为笛卡尔坐标(0~1)
    var posCartesian = geodeticSurfaceNormalCartographic(cartographic);
    /**
     *  我的理解2000坐标系是椭球，实际是有误差，最大误差为0.4%，不知道是不是理解错误
     *  实际意义不大，为了后面化简
     *  xR = (x / R) * R * R
     *  yR = (y / R) * R * R
     *  zR = (z / R) * R * R
     *  有误差，实际值  rr < R < r
     *  xR ~= (x / R) * r * r = x*r*r/R
     *  yR ~= (y / R) * r * r = y*r*r/R
     *  zR ~= (z / R) * rr * rr = z*rr*rr/R
     */
    var eldCartesian = multiplyComponents(Ellipsoid.radiiSquared, posCartesian);
    /**
     *      计算点距离球心距离，依然是上面误差的问题
     *       _______________________________
     *  d = √ xR * x/R + yR * y/R + zR * z/R
     *       _________________
     *    = √ x*x + y*y + z*z
     *
     *      实际计算公式  rr < R < r
     *       ________________________________________________
     *  d = √ x*r*r/R * x/R + y*r*r/R * y/R + z*rr*rr/R * z/R
     *       ________________________________________________
     *    = √ x*x * r*r/R*R + y*y * r*r/R*R + z*z * rr*rr/R*R
     *       ____________________________________
     *      √ x*x * r*r + y*y * r*r + z*z * rr*rr
     *    = -------------------------------------
     *                     R
     *
     */
    var gamma = Math.sqrt(cartesianDot(posCartesian, eldCartesian));
    /**
     *  猜测值,算是缩小误差  rr < R < r
     *  R ~= d
     *  xR / d = x
     *  yR / d = y
     *  zR / d = z
     *
     *
     *  实际值
     *                      ____________________________________
     *  x = (x*r*r/R) * (R/√ x*x * r*r + y*y * r*r + z*z * rr*rr)
     *             ____________________________________
     *    = x*r*r/√ x*x * r*r + y*y * r*r + z*z * rr*rr
     *                      ____________________________________
     *  y = (y*r*r/R) * (R/√ x*x * r*r + y*y * r*r + z*z * rr*rr)
     *             ____________________________________
     *    = y*r*r/√ x*x * r*r + y*y * r*r + z*z * rr*rr
     *                        ____________________________________
     *  z = (z*rr*rr/R) * (R/√ x*x * r*r + y*y * r*r + z*z * rr*rr)
     *               ____________________________________
     *    = z*rr*rr/√ x*x * r*r + y*y * r*r + z*z * rr*rr
     *
     */
    var divide = divideByScalar(eldCartesian, gamma);
    /**
     *  x/R * h
     *  y/R * h
     *  z/R * h
     */
    var multiply = multiplyByScalar(posCartesian, 0);

    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    //Cartesian3.add(divide, multiply, result);因为高度是0，所以可以不加
    return divide;
}
function multiplyByVector(matrix, cartesian) {
    var result = { x: 0.0, y: 0.0, z: 0.0, w: 0.0 };
    var vX = cartesian.x;
    var vY = cartesian.y;
    var vZ = cartesian.z;
    var vW = cartesian.w;

    var x = matrix[0] * vX + matrix[4] * vY + matrix[8] * vZ + matrix[12] * vW;
    var y = matrix[1] * vX + matrix[5] * vY + matrix[9] * vZ + matrix[13] * vW;
    var z = matrix[2] * vX + matrix[6] * vY + matrix[10] * vZ + matrix[14] * vW;
    var w = matrix[3] * vX + matrix[7] * vY + matrix[11] * vZ + matrix[15] * vW;

    result.x = x;
    result.y = y;
    result.z = z;
    result.w = w;
    return result;
}
function fromElements(x, y, z, w) {
    var result = { x: 0.0, y: 0.0, z: 0.0, w: 0.0 };
    result.x = x;
    result.y = y;
    result.z = z;
    result.w = w;
    return result;
}
function normalize(cartesian) {
    var result = {
        x: 0.0,
        y: 0.0,
        z: 0.0,
    };
    var magnitude = Math.sqrt(
        cartesian.x * cartesian.x +
            cartesian.y * cartesian.y +
            cartesian.z * cartesian.z
    );

    result.x = cartesian.x / magnitude;
    result.y = cartesian.y / magnitude;
    result.z = cartesian.z / magnitude;

    return result;
}
function worldToClip(position, viewMatrix, projectionMatrix) {
    var position4 = fromElements(position.x, position.y, position.z, 1);
    var positionEC = multiplyByVector(viewMatrix, position4);
    return multiplyByVector(projectionMatrix, positionEC);
}
function computeViewportTransformation(
    viewport,
    nearDepthRange,
    farDepthRange
) {
    var result = [];
    var x = viewport.x;
    var y = viewport.y;
    var width = viewport.width;
    var height = viewport.height;

    var halfWidth = width * 0.5;
    var halfHeight = height * 0.5;
    var halfDepth = (farDepthRange - nearDepthRange) * 0.5;

    var column0Row0 = halfWidth;
    var column1Row1 = halfHeight;
    var column2Row2 = halfDepth;
    var column3Row0 = x + halfWidth;
    var column3Row1 = y + halfHeight;
    var column3Row2 = nearDepthRange + halfDepth;
    var column3Row3 = 1.0;

    result[0] = column0Row0;
    result[1] = 0.0;
    result[2] = 0.0;
    result[3] = 0.0;
    result[4] = 0.0;
    result[5] = column1Row1;
    result[6] = 0.0;
    result[7] = 0.0;
    result[8] = 0.0;
    result[9] = 0.0;
    result[10] = column2Row2;
    result[11] = 0.0;
    result[12] = column3Row0;
    result[13] = column3Row1;
    result[14] = column3Row2;
    result[15] = column3Row3;
    return result;
}
function multiplyByPoint(matrix, cartesian) {
    var result = { x: 0.0, y: 0.0, z: 0.0 };
    var vX = cartesian.x;
    var vY = cartesian.y;
    var vZ = cartesian.z;

    var x = matrix[0] * vX + matrix[4] * vY + matrix[8] * vZ + matrix[12];
    var y = matrix[1] * vX + matrix[5] * vY + matrix[9] * vZ + matrix[13];
    var z = matrix[2] * vX + matrix[6] * vY + matrix[10] * vZ + matrix[14];

    result.x = x;
    result.y = y;
    result.z = z;
    return result;
}
function clipToGLWindowCoordinates(viewport, position) {
    var positionNDC = divideByScalar(position, position.w);

    // Viewport transform to transform from clip coordinates to window coordinates
    var viewportTransform = computeViewportTransformation(viewport, 0.0, 1.0);
    var positionWC = multiplyByPoint(viewportTransform, positionNDC);

    return { x: positionWC.x, y: positionWC.y };
}
var positionCC = { x: 0.0, y: 0.0, z: 0.0, w: 0.0 };
//笛卡尔坐标转屏幕像素
/**
 * scene = {
 *          width:scene.canvas.clientWidth,
 *          height:scene.canvas.clientHeight,
 *          viewMatrix:scene.camera.viewMatrix,([16])
 *          projectionMatrix:scene.camera.frustum.projectionMatrix([16])
 *         }
 *
 *
 * position = {
 *              x:0.0,
 *              y:0.0,
 *              z:0.0
 *             }
 *
 *
 */
function cartesianToWindow(scene, position) {
    var actualPosition = position;
    //获取画布信息
    var viewport = { x: 0.0, y: 0.0, width: 0.0, height: 0.0 };
    viewport.x = 0;
    viewport.y = 0;
    viewport.width = scene.width;
    viewport.height = scene.height;

    var cameraCentered = true;
    var result = { x: 0.0, y: 0.0 };
    //暂时只考虑三维
    if (cameraCentered) {
        // 从世界坐标转换到裁切坐标的投影矩阵
        positionCC = worldToClip(
            actualPosition,
            scene.viewMatrix,
            scene.projectionMatrix
        );
        if (positionCC.z < 0) {
            return undefined;
        }

        result = clipToGLWindowCoordinates(viewport, positionCC);
    }

    result.y = scene.height - result.y;
    return result;
}
function projection(point) {
    //经纬度转笛卡尔坐标
    var cartesian = degreesToCartesian(point);
    //迪卡尔转像素
    var pixel = cartesianToWindow(myScene, cartesian);
    return pixel;
}