Ice Permittivity Model

In A summary of the complex dieletric permittivity of ice in the megahertz range and its applications for radar sounding of polar ice sheets Fujita et Al. present a model for epsilon double prime. Their data ( with depth added ) is available here.

They fit to the following equation:

However when fitting curves through the supplied data points it became clear that fitting for the C parameter was going to be difficult. Instead it was fit to:

The best fit found is implemented by the following python code: 

# x_in is the ice temp in deg C
def fit_matsuoka_a(x_in):
    temp = 0.0
    # coefficients
    A = -6.3293731455351277E-05
    B = 1.6797867246920113E-08
    C = -1.4744880953559558E-12
    D = 9.1142297584796472E-02
    E = -1.3061526747480247E+01
    F = 1.1626675440193469E+03
    G = -4.2940915550767430E+04
    temp = A + B * math.pow(x_in, 2.0) + \
        C * math.pow(x_in, 4.0) + \
        D / math.pow(x_in, 2.0) + \
        E / math.pow(x_in, 4.0) + \
        F / math.pow(x_in, 6.0) + \
        G / math.pow(x_in, 8.0)
    
    return temp

# x_in is the ice temp in deg C
def fit_matsuoka_b(x_in):
    temp = 0.0
    # coefficients
    A = 3.6656269531875900E-05
    B = -6.2610754337689568E-09
    C = 4.5947045594069011E-13
    D = 5.3986512806528447E-03
    E = -6.1224071031577942E-01
    F = -2.5506123420729601E+01
    G = 4.9422024952720685E+03
    temp = A + B * math.pow(x_in, 2.0) + \
        C * math.pow(x_in, 4.0) + \
        D / math.pow(x_in, 2.0) + \
        E / math.pow(x_in, 4.0) + \
        F / math.pow(x_in, 6.0) + \
        G / math.pow(x_in, 8.0)
    
    return temp

# x_in is the ice temp in deg C
def fit_matsuoka_c(x_in):
    temp = 0.0
    # coefficients
    A = 9.0628357600078143E-08
    B = 1.7535861695933664E-12
    C = -4.6909178338196244E-16
    D = -1.5145302144421151E-05
    E = 2.8959559913410331E-03
    F = -1.7486054687356639E-01
    temp = A + B * math.pow(x_in, 2.0) + \
        C * math.pow(x_in, 4.0) + \
        D / math.pow(x_in, 2.0) + \
        E / math.pow(x_in, 4.0) + \
        F / math.pow(x_in, 6.0)
    return temp