# coding: cp1251
__author__ = 'Fliegende Hollander'
from math import sqrt, sin, cos, tan, fabs, radians, floor, pi
#========================================================================================
# Вспомогательные процедуры пересчета координат
#========================================================================================
# Расчет параметров геоида
def geoid_params():
    """
    Процедура вычисления характеристик геоида
    :return:    Набор характеристик геоида
    """
    a = 6378137.0
    b = 6356752.3142
    f    = (a - b) / a
    e_2  = f * (2 - f)
    e_21 = e_2 / (1 - e_2)
    n    = f / (2 - f)
    A = a * (1 - n + 5 * (n ** 2 - n ** 3) / 4. + 81 * (n ** 4 - n ** 5) / 64.)
    B = 3 * a * (n - n ** 2 + 7 * (n ** 3 - n ** 4) / 8. + 55 * (n ** 5) / 64.) / 2.
    C = 15 * a * (n ** 2 - n ** 3 + 3 * (n ** 4 - n ** 5) / 4.) / 16.
    D = 35 * a * (n ** 3 - n ** 4 + 11 * (n ** 5) / 16.) /48.
    E = 315 * a * (n ** 4 - n ** 5) / 512.

    return a, b, f, e_2, e_21, n, A, B, C, D, E

# Расчет степение синуса угла (для ускорения расчета
def sins(fi):
    """
    Процедура вычисления степеней синуса широты (для ускорения вычислений)
    :param fi:    Широта в радианах
    :return:      Массив степеней синуса
    """
    s = []
    for i in range(1, 9):
        s.append(sin(fi) ** i)

    return s

# Расчет степеней косинуса угла
def coss(fi):
    """
    Процедура вычисления степеней косинуса широты (для ускорения вычислений)
    :param fi:    Широта в радианах
    :return:      Массив степеней косинуса
    """
    s = []
    for i in range(1, 9):
        s.append(cos(fi) ** i)

    return s

# Расчет степеней тангенса угла
def tans(fi):
    """
    Процедура вычисления степеней тангенса широты (для ускорения вычислений)
    :param fi:    Широта в радианах
    :return:      Массив степеней тангенса
    """
    s = []
    for i in range(1, 9):
        s.append(tan(fi) ** i)

    return s

# Расчет длины меридианной дуги для заданной широты
def S(fi):
    """
    Расчет длины меридианной дуги
    :param fi:  Широта точки в радианах
    :return:    Длина меридиональной дуги
    """
    geoid = geoid_params()
    A = geoid[6]
    B = geoid[7]
    C = geoid[8]
    D = geoid[9]
    E = geoid[10]
    S = A * fi - B * sin(2 * fi) + C * sin(4 * fi) - D * sin(6 * fi) + E * sin(8 * fi)

    return S

# Расчет широты основаня перпендикуляра от точки к центральному меридиану зоны
def UTM_footpoint(N, hemisphere, eps = 0.000000001):
    """
    Вычисление широты основания перпендикуляра от точки до центрального меридиана зоны
    :param N:      Координата N точки
    :param eps:    Точность вычислений
    :return:       Footpoint-широта заданной точки
    """
    k0 = 0.9996
    fi = 0
    dN  = 1
    # Максимальное значение дуги
    Nm  = S(radians(84.)) * k0 if hemisphere == 'N' else S(radians(80.)) * k0
    # Итерационной вычисление искомой широты
    while fabs(dN) > eps:
        N0 = S(fi) * k0
        dN = N - N0
        fi += dN / Nm

    return fi

# Вычисление термов для прямого преобразования WGS84 -> UTM
def T_WGS84_UTM(fi):
    """
    Вычисляет термы для прямого преобразования WGS84 в UTM
    :param fi:      Широта в радианах
    :return:        Массив термов
    """
    k0 = 0.9996
    g = geoid_params()
    s = sins(fi)
    c = coss(fi)
    t = tans(fi)
    nu  = g[0] / sqrt(1 - g[3] * s[1])
    T1  = S(fi) * k0
    T2  = nu * s[0] * c[0] * k0 / 2.
    T3  = nu * s[0] * c[2] * k0 * (5 - t[1] + 9 * g[4] * c[1] + 4 * g[4] ** 2 * c[3]) / 24.
    T4  = nu * s[0] * c[4] * k0 * (61 - 58 * t[1] + t[3] * 270 * g[4] * c[1] - 330 * t[1] * g[4] * c[1] + \
          445 * g[4] ** 2 *c[3] + 324 * g[4] ** 3 * c[5] - 680 * t[1] * g[4] ** 2 * c[3] + \
          88 * g[4] ** 4 * c[7] - 600 * t[1] * g[4] ** 3 * c[5] - 192 *t[1] * g[4] ** 4 * c[7]) / 720.
    T5  = nu * s[0] * c[6] * k0 * (1385 - 3111 * t[1] + 543 * t[3] - t[5]) / 40320.
    T6  = nu * c[0] * k0
    T7  = nu * c[2] * k0 * (1 - t[1] + g[4] * c[1]) / 6.
    T8  = nu * c[4] * k0 * (5 - 18 * t[1] + t[3] + 14 * g[4] * c[1] - 58 *t[1] * g[4] * c[1] + \
          13 * g[4] ** 2 * c[3] + 4 * g[4] ** 3 * c[5] - 64 * t[1] * g[4] ** 2 * c[3] - 24 * t[1] * g[4] ** 3 * c[5]) / 120.
    T9  = nu * c[6] * k0 * (61 - 479 * t[1] + 179 * t[3] - t[5]) / 5040.

    return T1, T2, T3, T4, T5, T6, T7, T8, T9

# Вычисление термов для обратного преобразования UTM -> WGS84
def T_UTM_WGS84(N, hemisphere):
    """
    Вычисляет термы для обратного преобразования UTM и WGS84
    :param N:           Северное склонение точки
    :param hemisphere:  Полушарие
    :return:            Массив термов
    """
    k0 = 0.9996
    fi1 = UTM_footpoint(N, hemisphere)

    g = geoid_params()
    s = sins(fi1)
    c = coss(fi1)
    t = tans(fi1)
    nu  = g[0] / sqrt(1 - g[3] * s[1])
    ro  = g[0] * (1 - g[3]) / sqrt((1 - g[3] * s[1]) ** 3)
    T10 = t[0] / (2 * ro * nu * k0 ** 2)
    T11 = t[0] * (5 + 3 * t[1] + g[4] * c[1] - 4 * g[4] ** 2 * c[3] - 9 * t[1] * g[4] * c[1]) / (24 * ro * nu ** 3 * k0 ** 4)
    T12 = t[0] * (61 + 90 * t[1] + 46 * g[4] * c[1] + 45 * t[3] - 252 * t[1] * g[4] * c[1] - \
          3 * g[4] ** 2 * c[3] + 100 * g[4] ** 3 * c[5] - 66 * t[1] * g[4] ** 2 * c[3] - \
          90 * t[3] * g[4] * c[1] + 88 * g[4] ** 4 * c[7] + 225 * t[3] * g[4] ** 2 * c[3] + \
          84 * t[1] * g[4] ** 3 * c[5] - 192 * t[1] * g[4] ** 4 * c[7]) / (720 * ro * nu ** 5 * k0 ** 6)
    T13 = t[1] * (1385 + 3633 * t[1] + 4095 * t[3] + 1575 * t[5]) / (40320 * ro * nu ** 7 * k0 ** 8)
    T14 = 1. / (nu * c[0] * k0)
    T15 = (1 + 2 * t[1] + g[4] * c[1]) / (6 * nu ** 3 * c[0] * k0 ** 3)
    T16 = (5 + 6 * g[4] * c[1] + 28 * t[1] - 3 * g[4] ** 2 * c[3] + 8 * t[1] * g[4] * c[1] + \
           24 * t[3] - 4 * g[4] ** 3 * c[5] + 4 * t[1] * g[4] ** 2 * c[3] + 24 * t[1] * g[4] ** 3 * c[5]) / (120 * nu ** 5 * c[0] * k0 ** 5)
    T17 = (61 + 6621 * t[1] + 1320 * t[3] + 720 * t[5]) / (5040 * nu ** 7 * c[0] * k0 ** 7)

    return T10, T11, T12, T13, T14, T15, T16, T17, fi1

#========================================================================================
# Процедуры пересчета координат WGS84 <-> UTM
#========================================================================================
# Прямое преобразование
def WGS84_UTM(fi, hemisphere, lo, easting):
    """
    Процедура прямого преобразования координат WGS84 -> UTM
    :param fi:          Широта в радианах
    :param hemisphere:  Полушарие N - северное, S - южное
    :param lo:          Долгота в радианах
    :param easting:     Положение относительно нулевого меридиана E - долгота восточная, W - западная
    :return:            Координаты точки в проекции Universal Transverse Mercator, номер зоны
    """
    FN = 0 if hemisphere == 'N' else 10000000
    FE = 500000
    zone = (30 + floor(lo / 6) + 1) if easting == "E" else (floor((180 - lo) / 6) + 1)
    lo_0 = radians((zone - 1) * 6 - 180 + 3) if easting == 'E' else radians(180 - (zone - 1) * 6 + 3)
    dlo  = radians(lo) - lo_0

    T = T_WGS84_UTM(radians(fi))

    N = FN + (T[0] + dlo ** 2 * T[1] + dlo ** 4 * T[2] + dlo ** 6 * T[3] + dlo ** 8 * T[4])
    E = FE + (dlo * T[5] + dlo ** 3 * T[6] + dlo ** 5 * T[7] + dlo ** 7 * T[8])

    return N, E, zone

# Обратное преобразование
def UTM_WGS84(N, E, zone, hemisphere):
    """
    Процедура обратного преобразования UTM -> WGS84
    :param N:           Северное склонение в метрах
    :param E:           Восточное склонение в метрах
    :param zone:        Номер зоны UTM
    :param hemisphere:  Полушарие - северное или южное
    :return:            Геодезические широта и долгота в радианах
    """
    lo_0 = radians((zone - 1) * 6 - 180 + 3) if zone > 30 else radians(180 - (zone - 1) * 6 + 3)
    FE = 500000
    DE = E - FE
    T = T_UTM_WGS84(N, hemisphere)

    lat = T[8] - DE ** 2 * T[0] + DE ** 4 * T[1] - DE ** 6 * T[2] + DE ** 8 * T[3]
    lon = lo_0 + DE * T[4] - DE ** 3 * T[5] + DE ** 5 * T[6] - DE ** 7 * T[7]

    return lat, lon

def MIN_MAX_EST(zone):
    """
    Вычисление максимального и минимального восточного смещеня заданной зоны
    :param zone:    Номер зоны
    :return:        Массив максимальных и минимальных восточных смещений для зоны от широты
    """
    lon_0   = radians((zone - 1) * 6 - 180 + 3) if zone > 30 else radians(180 - (zone - 1) * 6 + 3)
    lon_min = lon_0 - radians(3.)
    lon_max = lon_0 + radians(3.)
    fi      = radians(0.)
    dfi     = radians(1.)
    M = []
    for i in range(85):
        N, Emin, zone = WGS84_UTM(fi, 'N', lon_min, 'E')
        N, Emax, zone = WGS84_UTM(fi, 'N', lon_max, 'E')
        M.append([(180 * fi) / pi, Emin, Emax])
        fi += dfi

    return M

if __name__ == "__main__":
    # Тестовый расчет по известным координатам опорной точки
    lat = 49.4345794666
    lon = 27.0021217333

    N, E, zone = WGS84_UTM(lat, 'N', lon, 'E')
    print "NORTHING IS : ", N
    print "EASTING IS  : ", E
    print "ZONE IS     : ", zone

    lat_1, lon_1 = UTM_WGS84(N, E, zone, 'N')
    print "-----------------------------------------------------------------"
    print "LATITUDE IS  : ", 180 * lat_1 / pi
    print "LONGITUDE IS : ", 180 * lon_1 / pi

    # Тест вычисления пределов зоны по восточному смещению
    M = MIN_MAX_EST(zone)
    print "*****************************************************************"
    print " Zone ", zone, " boundaries"
    print "*****************************************************************"
    for i in range(len(M)):
        print "LAT = ", M[i][0], " \t", "Emin = ", M[i][1], "\t", "Emax = ", M[i][2]


