worldmagmodel.cpp

                    coeff += (decimal_date - MagneticModel.epoch) * get_secular_var_coeff_h(sum_index);
            }
        }

        return coeff;
    }

    double WorldMagModel::get_secular_var_coeff_g(int index)
    {
        if (index >= WMM_NUMTERMS)
            return 0;

        return CoeffFile[index][4];
    }

    double WorldMagModel::get_secular_var_coeff_h(int index)
    {
        if (index >= WMM_NUMTERMS)
            return 0;

        return CoeffFile[index][5];
    }

    int WorldMagModel::DateToYear(int month, int day, int year)
    {
        // Converts a given calendar date into a decimal year

        int temp = 0;	// Total number of days
        int MonthDays[13] = { 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
        int ExtraDay = 0;

        if ((year % 4 == 0 && year % 100 != 0) || (year % 400 == 0))
            ExtraDay = 1;
        MonthDays[2] += ExtraDay;

        /******************Validation********************************/

        if (month <= 0 || month > 12)
            return -1;  // error

        if (day <= 0 || day > MonthDays[month])
            return -2;  // error

        /****************Calculation of t***************************/
        for (int i = 1; i <= month; i++)
            temp += MonthDays[i - 1];
        temp += day;

        decimal_date = year + (temp - 1) / (365.0 + ExtraDay);

        return 0;   // OK
    }

    void WorldMagModel::GeodeticToSpherical(WMMtype_CoordGeodetic *CoordGeodetic, WMMtype_CoordSpherical *CoordSpherical)
    {
        // Converts Geodetic coordinates to Spherical coordinates
        // Convert geodetic coordinates, (defined by the WGS-84
        // reference ellipsoid), to Earth Centered Earth Fixed Cartesian
        // coordinates, and then to spherical coordinates.

        double CosLat = cos(DEG2RAD(CoordGeodetic->phi));
        double SinLat = sin(DEG2RAD(CoordGeodetic->phi));

        // compute the local radius of curvature on the WGS-84 reference ellipsoid
        double rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);

        // compute ECEF Cartesian coordinates of specified point (for longitude=0)
        double xp = (rc + CoordGeodetic->HeightAboveEllipsoid) * CosLat;
        double zp = (rc * (1.0 - Ellip.epssq) + CoordGeodetic->HeightAboveEllipsoid) * SinLat;

        // compute spherical radius and angle lambda and phi of specified point
        CoordSpherical->r = sqrt(xp * xp + zp * zp);
        CoordSpherical->phig = RAD2DEG(asin(zp / CoordSpherical->r));	// geocentric latitude
        CoordSpherical->lambda = CoordGeodetic->lambda;	// longitude
    }

}