def magsyn(gh, sv, b, date, itype, alt, colat, elong):
    """
# Computes x, y, z, and f for a given date and position, from the
# spherical harmonic coefficients of the International Geomagnetic
# Reference Field (IGRF).
# From Malin and Barraclough (1981), Computers and Geosciences, V.7, 401-405.
#
# Input:
#       date  = Required date in years and decimals of a year (A.D.)
#       itype = 1, if geodetic coordinates are used, 2 if geocentric
#       alt   = height above mean sea level in km (if itype = 1)
#       alt   = radial distance from the center of the earth (itype = 2)
#       colat = colatitude in degrees (0 to 180)
#       elong = east longitude in degrees (0 to 360)
#               gh        = main field values for date (calc. in igrf subroutine)
#               sv        = secular variation coefficients (calc. in igrf subroutine)
#               begin = date of dgrf (or igrf) field prior to required date
#
# Output:
#       x     - north component of the magnetic force in nT
#       y     - east component of the magnetic force in nT
#       z     - downward component of the magnetic force in nT
#       f     - total magnetic force in nT
#
#       NB: the coordinate system for x,y, and z is the same as that specified
#       by itype.
#
# Modified 4/9/97 to use DGRFs from 1945 to 1990 IGRF
# Modified 10/13/06 to use  1995 DGRF, 2005 IGRF and sv coefficient
# for extrapolation beyond 2005. Coefficients from Barton et al. PEPI, 97: 23-26
# (1996), via web site for NOAA, World Data Center A. Modified to use
#degree and
# order 10 as per notes in Malin and Barraclough (1981).
# coefficients for DGRF 1995 and IGRF 2005 are from http://nssdcftp.gsfc.nasa.gov/models/geomagnetic/igrf/fortran_code/
# igrf subroutine calculates
# the proper main field and secular variation coefficients (interpolated between
# dgrf values or extrapolated from 1995 sv values as appropriate).
    """
#
#       real gh(120),sv(120),p(66),q(66),cl(10),sl(10)
#               real begin,dateq
    p = np.zeros((66), 'f')
    q = np.zeros((66), 'f')
    cl = np.zeros((10), 'f')
    sl = np.zeros((10), 'f')
    begin = b
    t = date - begin
    r = alt
    one = colat * 0.0174532925
    ct = np.cos(one)
    st = np.sin(one)
    one = elong * 0.0174532925
    cl[0] = np.cos(one)
    sl[0] = np.sin(one)
    x, y, z = 0.0, 0.0, 0.0
    cd, sd = 1.0, 0.0
    l, ll, m, n = 1, 0, 1, 0
    if itype != 2:
        #
        # if required, convert from geodectic to geocentric
        a2 = 40680925.0
        b2 = 40408585.0
        one = a2 * st * st
        two = b2 * ct * ct
        three = one + two
        rho = np.sqrt(three)
        r = np.sqrt(alt * (alt + 2.0 * rho) +
                    old_div((a2 * one + b2 * two), three))
        cd = old_div((alt + rho), r)
        sd = (a2 - b2) / rho * ct * st / r
        one = ct
        ct = ct * cd - st * sd
        st = st * cd + one * sd
    ratio = old_div(6371.2, r)
    rr = ratio * ratio
#
# compute Schmidt quasi-normal coefficients p and x(=q)
    p[0] = 1.0
    p[2] = st
    q[0] = 0.0
    q[2] = ct
    for k in range(1, 66):
        if n < m:   # else go to 2
            m = 0
            n = n + 1
            rr = rr * ratio
            fn = n
            gn = n - 1
# 2
        fm = m
        if k != 2:  # else go to 4
            if m == n:   # else go to 3
                one = np.sqrt(1.0 - old_div(0.5, fm))
                j = k - n - 1
                p[k] = one * st * p[j]
                q[k] = one * (st * q[j] + ct * p[j])
                cl[m - 1] = cl[m - 2] * cl[0] - sl[m - 2] * sl[0]
                sl[m - 1] = sl[m - 2] * cl[0] + cl[m - 2] * sl[0]
            else:
                # 3
                gm = m * m
                one = np.sqrt(fn * fn - gm)
                two = old_div(np.sqrt(gn * gn - gm), one)
                three = old_div((fn + gn), one)
                i = k - n
                j = i - n + 1
                p[k] = three * ct * p[i] - two * p[j]
                q[k] = three * (ct * q[i] - st * p[i]) - two * q[j]
#
# synthesize x, y, and z in geocentric coordinates.
# 4
        #print (l,ll,t,rr)
        one = (gh[l - 1] + sv[ll + l - 1] * t) * rr
        if m != 0:  # else go to 7
            two = (gh[l] + sv[ll + l] * t) * rr
            three = one * cl[m - 1] + two * sl[m - 1]
            x = x + three * q[k]
            z = z - (fn + 1.0) * three * p[k]
            if st != 0.0:  # else go to 5
                y = y + (one * sl[m - 1] - two * cl[m - 1]) * fm * p[k] / st
            else:
                # 5
                y = y + (one * sl[m - 1] - two * cl[m - 1]) * q[k] * ct
            l = l + 2
        else:
            # 7
            x = x + one * q[k]
            z = z - (fn + 1.0) * one * p[k]
            l = l + 1
        m = m + 1
#
# convert to coordinate system specified by itype
    one = x
    x = x * cd + z * sd
    z = z * cd - one * sd
    f = np.sqrt(x * x + y * y + z * z)
#
    return x, y, z, f