/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "refellipsoid.h"

gctl::refellipsoid::refellipsoid(/* args */){}

gctl::refellipsoid::refellipsoid(refellipsoid_type_e refellipsoid)
{
    set(refellipsoid);
}

gctl::refellipsoid::refellipsoid(double R, double r, bool flat)
{
    set(R, r, flat);
}

gctl::refellipsoid::~refellipsoid(){}

void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid)
{
    if (refellipsoid == Earth) {r_ = R_ = GCTL_Earth_Radius;}
    else if (refellipsoid == MagEarth) {r_ = R_ = GCTL_Earth_RefRadius;}
    else if (refellipsoid == Moon) {r_ = R_ = GCTL_Moon_Radius;}
    else if (refellipsoid == Mars) {r_ = R_ = GCTL_Mars_Radius;}
    else if (refellipsoid == WGS84) {r_ = GCTL_WGS84_PoleRadius; R_ = GCTL_WGS84_EquatorRadius;}
    else if (refellipsoid == Ardalan2010) {r_ = GCTL_Mars_Ardalan2010_PoleRadius; R_ = GCTL_Mars_Ardalan2010_EquatorRadius;}
    else throw std::invalid_argument("Invalid reference system type for gctl::refellipsoid::set(...)");

    f_ = (R_ - r_)/R_;
    e_ = sqrt(1.0 - (r_*r_)/(R_*R_));
    return;
}

void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat)
{
    R_ = R;
    if (is_flat) r_ = R_*(1.0 - 1.0/r_or_flat);
    else r_ = r_or_flat;

    f_ = (R_ - r_)/R_;
    e_ = sqrt(1.0 - (r_*r_)/(R_*R_));
    return;
}

double gctl::refellipsoid::geodetic_radius(double geodetic_lati)
{
    return R_/sqrt(1.0 - e_*e_*sind(geodetic_lati)*sind(geodetic_lati));
}

void gctl::refellipsoid::geodetic2spherical(double geodetic_lati, 
    double geodetic_hei, double& sph_lati, double& sph_rad)
{
    double CosLat, SinLat, rc, xp, zp;
    /*
    ** Convert geodetic coordinates, (defined by the WGS-84
    ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
    ** coordinates, and then to spherical coordinates.
    */
    CosLat = cosd(geodetic_lati);
    SinLat = sind(geodetic_lati);

    // compute the local rRdius of curvature on the WGS-84 reference ellipsoid
    rc = R_/sqrt(1.0 - e_*e_*SinLat*SinLat);

    // compute ECEF Cartesian coordinates of specified point (for longitude=0)
    xp = (rc + geodetic_hei)*CosLat;
    zp = (rc*(1.0 - e_*e_) + geodetic_hei)*SinLat;

    // compute spherical rRdius and angle lambda and phi of specified point
    sph_rad = sqrt(xp*xp + zp*zp);
    sph_lati = deg(asin(zp/sph_rad));
    return;
}

void gctl::refellipsoid::spherical2geodetic(const point3ds& ps, double& geodetic_lon, double& geodetic_lati, 
    double& geodetic_hei, double eps, int cnt)
{
    point3dc pc = ps.s2c();
    xyz2geodetic(pc, geodetic_lon, geodetic_lati, geodetic_hei, eps, cnt);
    return;
}

void gctl::refellipsoid::xyz2geodetic(const point3dc& pc, double& geodetic_lon, 
    double& geodetic_lati, double& geodetic_hei, double eps, int cnt)
{
	double curB, N;
    double H = sqrt(pc.x*pc.x + pc.y*pc.y);
	double calB = atan2(pc.z, H); 
	
	int c = 0;
	do
	{
		curB = calB;
		N = R_/sqrt(1 - e_*e_*sin(curB)*sin(curB));
		calB = atan2(pc.z + N*e_*e_*sin(curB), H);
		c++;
	}
    while (abs(curB - calB)*180.0/GCTL_Pi > eps  && c < cnt);
	
	geodetic_lon = atan2(pc.y, pc.x)*180.0/GCTL_Pi;
	geodetic_lati = curB*180.0/GCTL_Pi;
	geodetic_hei = pc.z/sin(curB) - N*(1 - e_*e_);
    return;
}

void gctl::refellipsoid::geodetic2xyz(double geodetic_lon, double geodetic_lati, 
    double geodetic_hei, point3dc& pc)
{
    double L = arc(geodetic_lon);
	double B = arc(geodetic_lati);
	double H = geodetic_hei;

	double N = R_/sqrt(1 - e_*e_*sin(B)*sin(B));
	pc.x = (N + H)*cos(B)*cos(L);
	pc.y = (N + H)*cos(B)*sin(L);
	pc.z = (N*(1 - e_*e_) + H)*sin(B);
    return;
}