# -*- coding: utf-8 -*-
"""
@author: zadavi
"""

import numpy as np
import pandas as pd   
def gmf(dmjd = None,dlat = None,dlon = None,dhgt = None,zd = None): 
    # This subroutine determines the Global Mapping Functions GMF
# Reference: Boehm, J., A.E. Niell, P. Tregoning, H. Schuh (2006),
# Global Mapping Functions (GMF): A new empirical mapping function based on numerical weather model data,
# Geoph. Res. Letters, Vol. 33, L07304, doi:10.1029/2005GL025545.
    
    # input data
# ----------
# dmjd: modified julian date
# dlat: ellipsoidal latitude in radians
# dlon: longitude in radians
# dhgt: height in m
# zd:   zenith distance in radians
    
# output data
# -----------
# gmfh: hydrostatic mapping function
# gmfw: wet mapping function
    
# Johannes Boehm, 2005 August 30
    
# ref 2006 Aug. 14: recursions for Legendre polynomials (O. Montenbruck)
# ref 2011 Jul. 21: latitude -> ellipsoidal latitude (J. Boehm)

#    File Converted to python script by Zohreh Adavi zohreh.adavi@tuwien.ac.at, February 2023
# =============================================================================  
    
    ah_mean = np.array([+ 125.17,+ 0.8503,+ 0.06936,- 6.76,+ 0.1771,+ 0.0113,+ 0.5963,+ 0.01808,+ 0.002801,- 0.001414,- 1.212,+ 0.093,+ 0.003683,+ 0.001095,+ 4.671e-05,+ 0.3959,- 0.03867,+ 0.005413,- 0.0005289,+ 0.0003229,+ 2.067e-05,+ 0.3,+ 0.02031,+ 0.0059,+ 0.0004573,- 7.619e-05,+ 2.327e-06,+ 3.845e-06,+ 0.1182,+ 0.01158,+ 0.005445,+ 6.219e-05,+ 4.204e-06,- 2.093e-06,+ 1.54e-07,- 4.28e-08,- 0.4751,- 0.0349,+ 0.001758,+ 0.0004019,- 2.799e-06,- 1.287e-06,+ 5.468e-07,+ 7.58e-08,- 6.3e-09,- 0.116,+ 0.008301,+ 0.0008771,+ 9.955e-05,- 1.718e-06,- 2.012e-06,+ 1.17e-08,+ 1.79e-08,- 1.3e-09,+ 1e-10])
    bh_mean = np.array([+ 0.0,+ 0.0,+ 0.03249,+ 0.0,+ 0.03324,+ 0.0185,+ 0.0,- 0.1115,+ 0.02519,+ 0.004923,+ 0.0,+ 0.02737,+ 0.01595,- 0.0007332,+ 0.0001933,+ 0.0,- 0.04796,+ 0.006381,- 0.0001599,- 0.0003685,+ 1.815e-05,+ 0.0,+ 0.07033,+ 0.002426,- 0.001111,- 0.0001357,- 7.828e-06,+ 2.547e-06,+ 0.0,+ 0.005779,+ 0.003133,- 0.0005312,- 2.028e-05,+ 2.323e-07,- 9.1e-08,- 1.65e-08,+ 0.0,+ 0.03688,- 0.0008638,- 8.514e-05,- 2.828e-05,+ 5.403e-07,+ 4.39e-07,+ 1.35e-08,+ 1.8e-09,+ 0.0,- 0.02736,- 0.0002977,+ 8.113e-05,+ 2.329e-07,+ 8.451e-07,+ 4.49e-08,- 8.1e-09,- 1.5e-09,+ 2e-10])
    ah_amp = np.array([- 0.2738,- 2.837,+ 0.01298,- 0.3588,+ 0.02413,+ 0.03427,- 0.7624,+ 0.07272,+ 0.0216,- 0.003385,+ 0.4424,+ 0.03722,+ 0.02195,- 0.001503,+ 0.0002426,+ 0.3013,+ 0.05762,+ 0.01019,- 0.0004476,+ 6.79e-05,+ 3.227e-05,+ 0.3123,- 0.03535,+ 0.00484,+ 3.025e-06,- 4.363e-05,+ 2.854e-07,- 1.286e-06,- 0.6725,- 0.0373,+ 0.0008964,+ 0.0001399,- 3.99e-06,+ 7.431e-06,- 2.796e-07,- 1.601e-07,+ 0.04068,- 0.01352,+ 0.0007282,+ 9.594e-05,+ 2.07e-06,- 9.62e-08,- 2.742e-07,- 6.37e-08,- 6.3e-09,+ 0.08625,- 0.005971,+ 0.0004705,+ 2.335e-05,+ 4.226e-06,+ 2.475e-07,- 8.85e-08,- 3.6e-08,- 2.9e-09,+ 0.0])
    bh_amp = np.array([+ 0.0,+ 0.0,- 0.1136,+ 0.0,- 0.1868,- 0.01399,+ 0.0,- 0.1043,+ 0.01175,- 0.00224,+ 0.0,- 0.03222,+ 0.01333,- 0.002647,- 2.316e-05,+ 0.0,+ 0.05339,+ 0.01107,- 0.003116,- 0.0001079,- 1.299e-05,+ 0.0,+ 0.004861,+ 0.008891,- 0.0006448,- 1.279e-05,+ 6.358e-06,- 1.417e-07,+ 0.0,+ 0.03041,+ 0.00115,- 0.0008743,- 2.781e-05,+ 6.367e-07,- 1.14e-08,- 4.2e-08,+ 0.0,- 0.02982,- 0.003,+ 1.394e-05,- 3.29e-05,- 1.705e-07,+ 7.44e-08,+ 2.72e-08,- 6.6e-09,+ 0.0,+ 0.01236,- 0.0009981,- 3.792e-05,- 1.355e-05,+ 1.162e-06,- 1.789e-07,+ 1.47e-08,- 2.4e-09,- 4e-10])
    aw_mean = np.array([+ 56.4,+ 1.555,- 1.011,- 3.975,+ 0.03171,+ 0.1065,+ 0.6175,+ 0.1376,+ 0.04229,+ 0.003028,+ 1.688,- 0.1692,+ 0.05478,+ 0.02473,+ 0.0006059,+ 2.278,+ 0.006614,- 0.0003505,- 0.006697,+ 0.0008402,+ 0.0007033,- 3.236,+ 0.2184,- 0.04611,- 0.01613,- 0.001604,+ 5.42e-05,+ 7.922e-05,- 0.2711,- 0.4406,- 0.03376,- 0.002801,- 0.000409,- 2.056e-05,+ 6.894e-06,+ 2.317e-06,+ 1.941,- 0.2562,+ 0.01598,+ 0.005449,+ 0.0003544,+ 1.148e-05,+ 7.503e-06,- 5.667e-07,- 3.66e-08,+ 0.8683,- 0.05931,- 0.001864,- 0.0001277,+ 0.0002029,+ 1.269e-05,+ 1.629e-06,+ 9.66e-08,- 1.015e-07,- 5e-10])
    bw_mean = np.array([+ 0.0,+ 0.0,+ 0.2592,+ 0.0,+ 0.02974,- 0.5471,+ 0.0,- 0.5926,- 0.103,- 0.01567,+ 0.0,+ 0.171,+ 0.09025,+ 0.02689,+ 0.002243,+ 0.0,+ 0.3439,+ 0.02402,+ 0.00541,+ 0.001601,+ 9.669e-05,+ 0.0,+ 0.09502,- 0.03063,- 0.001055,- 0.0001067,- 0.000113,+ 2.124e-05,+ 0.0,- 0.3129,+ 0.008463,+ 0.0002253,+ 7.413e-05,- 9.376e-05,- 1.606e-06,+ 2.06e-06,+ 0.0,+ 0.2739,+ 0.001167,- 2.246e-05,- 0.0001287,- 2.438e-05,- 7.561e-07,+ 1.158e-06,+ 4.95e-08,+ 0.0,- 0.1344,+ 0.005342,+ 0.0003775,- 6.756e-05,- 1.686e-06,- 1.184e-06,+ 2.768e-07,+ 2.73e-08,+ 5.7e-09])
    aw_amp = np.array([+ 0.1023,- 2.695,+ 0.3417,- 0.1405,+ 0.3175,+ 0.2116,+ 3.536,- 0.1505,- 0.0166,+ 0.02967,+ 0.3819,- 0.1695,- 0.07444,+ 0.007409,- 0.006262,- 1.836,- 0.01759,- 0.06256,- 0.002371,+ 0.0007947,+ 0.0001501,- 0.8603,- 0.136,- 0.03629,- 0.003706,- 0.0002976,+ 1.857e-05,+ 3.021e-05,+ 2.248,- 0.1178,+ 0.01255,+ 0.001134,- 0.0002161,- 5.817e-06,+ 8.836e-07,- 1.769e-07,+ 0.7313,- 0.1188,+ 0.01145,+ 0.001011,+ 0.0001083,+ 2.57e-06,- 2.14e-06,- 5.71e-08,+ 2e-08,- 1.632,- 0.006948,- 0.003893,+ 0.0008592,+ 7.577e-05,+ 4.539e-06,- 3.852e-07,- 2.213e-07,- 1.37e-08,+ 5.8e-09])
    bw_amp = np.array([+ 0.0,+ 0.0,- 0.08865,+ 0.0,- 0.4309,+ 0.0634,+ 0.0,+ 0.1162,+ 0.06176,- 0.004234,+ 0.0,+ 0.253,+ 0.04017,- 0.006204,+ 0.004977,+ 0.0,- 0.1737,- 0.005638,+ 0.0001488,+ 0.0004857,- 0.0001809,+ 0.0,- 0.1514,- 0.01685,+ 0.005333,- 7.611e-05,+ 2.394e-05,+ 8.195e-06,+ 0.0,+ 0.09326,- 0.01275,- 0.0003071,+ 5.374e-05,- 3.391e-05,- 7.436e-06,+ 6.747e-07,+ 0.0,- 0.08637,- 0.003807,- 0.0006833,- 3.861e-05,- 2.268e-05,+ 1.454e-06,+ 3.86e-07,- 1.068e-07,+ 0.0,- 0.02658,- 0.001947,+ 0.0007131,- 3.506e-05,+ 1.885e-07,+ 5.792e-07,+ 3.99e-08,+ 2e-08,- 5.7e-09])
    np.pi = 3.14159265359
    # reference day is 28 January
# this is taken from Niell (1996) to be consistent
    doy = dmjd - 44239 + 1 - 28
    
    # degree n and order m
    nmax = 9
    mmax = 9
    # unit vector
    x = np.cos(dlat) * np.cos(dlon)
    y = np.cos(dlat) * np.sin(dlon)
    z = np.sin(dlat)
    # Legendre polynomials
    V=pd.DataFrame(np.zeros((nmax+1,nmax+1)))
    W=pd.DataFrame(np.zeros((nmax+1,nmax+1)))
    
    
    V.loc[0,0] = 1
    W.loc[0,0] = 0
    V.loc[1,0] = z * V.loc[0,0]
    W.loc[1,0] = 0
    for n in range(2,nmax+1):
        V.loc[n ,0] = ((2 * n - 1) * z * V.loc[n-1,0] - (n - 1) * V.loc[n-2,0]) / n
        W.loc[n ,0] = 0
    
    for m in range(1,nmax+1):
        V.loc[m ,m ] = (2 * m - 1) * (x * V.loc[m-1,m-1]  - y * W.loc[m-1,m-1] )
        W.loc[m ,m ] = (2 * m - 1) * (x * W.loc[m-1,m-1] +  y * V.loc[m-1,m-1])
        if (m < nmax):
            V.loc[m + 1,m ] = (2 * m + 1) * z * V.loc[m,m]
            W.loc[m + 1,m ] = (2 * m + 1) * z * W.loc[m,m]
        for n in range(m+2,nmax+1):
            V.loc[n ,m ] = ((2 * n - 1) * z * V.loc[n-1,m] - (n + m - 1) * V.loc[n-2,m]) / (n - m)
            W.loc[n ,m ] = ((2 * n - 1) * z * W.loc[n-1,m] - (n + m - 1) * W.loc[n-2,m]) / (n - m)
    
    # (1) hydrostatic mf
    
    bh = 0.0029
    c0h = 0.062
    if (dlat < 0):
        phh = np.pi
        c11h = 0.007
        c10h = 0.002
    else:
        phh = 0
        c11h = 0.005
        c10h = 0.001
    
    ch = c0h + ((np.cos(doy / 365.25 * 2 * np.pi + phh) + 1) * c11h / 2 + c10h) * (1 - np.cos(dlat))
    ahm = 0
    aha = 0
    i = 0
    for n in range(0,nmax+1):
        for m in range(0,n+1):
            
            ahm = ahm + (ah_mean[i] * V.loc[n,m] + bh_mean[i]  * W.loc[n,m])
            aha = aha + (ah_amp[i]  * V.loc[n,m] + bh_amp[i]  * W.loc[n,m])
            i = i + 1
    
    ah = (ahm + aha * np.cos(doy / 365.25 * 2 * np.pi)) * 1e-05
    sine = np.sin(np.pi / 2 - zd)
    cose = np.cos(np.pi / 2 - zd)
    beta = bh / (sine + ch)
    gamma = ah / (sine + beta)
    topcon = (1 + ah / (1 + bh / (1 + ch)))
    gmfh = topcon / (sine + gamma)
    # height correction for hydrostatic mapping function from Niell (1996) in order to reduce the coefficients to sea level
    a_ht = 2.53e-05
    b_ht = 0.00549
    c_ht = 0.00114
    hs_km = dhgt / 1000
    beta = b_ht / (sine + c_ht)
    gamma = a_ht / (sine + beta)
    topcon = (1 + a_ht / (1 + b_ht / (1 + c_ht)))
    ht_corr_coef = 1 / sine - topcon / (sine + gamma)
    ht_corr = ht_corr_coef * hs_km
    gmfh = gmfh + ht_corr
    # (2) wet mf
    
    bw = 0.00146
    cw = 0.04391
    awm = 0
    awa = 0
    i = 0
    for n in range(0,nmax+1):
        for m in range(0,n+1):
            
            awm = awm + (aw_mean[i] * V.loc[n,m] + bw_mean[i] * W.loc[n,m])
            awa = awa + (aw_amp [i] * V.loc[n,m] + bw_amp[i] * W.loc[n,m])
            i = i + 1
    
    aw = (awm + awa * np.cos(doy / 365.25 * 2 * np.pi)) * 1e-05
    beta = bw / (sine + cw)
    gamma = aw / (sine + beta)
    topcon = (1 + aw / (1 + bw / (1 + cw)))
    gmfw = topcon / (sine + gamma)
    return gmfh,gmfw


###############################################################################
def gradient_MF(elev=None):
    
    c=0.0032
    elev=np.radians(elev)
    mfg=1/(np.sin(elev)* np.tan(elev)+c)
    mfg=mfg/1.4
    return mfg