From 77191da3d05810079200000a6cc60113bc898fc3 Mon Sep 17 00:00:00 2001 From: Chris Mackey Date: Mon, 4 Mar 2024 12:01:12 -0800 Subject: [PATCH] fix(north): Add a class to convert magnetic to true north --- MANIFEST.in | 1 + ladybug/WMM.COF | 93 ++++++++++++ ladybug/north.py | 295 +++++++++++++++++++++++++++++++++++++++ tests/assets/txt/WMM.COF | 93 ++++++++++++ tests/north_test.py | 47 +++++++ 5 files changed, 529 insertions(+) create mode 100644 ladybug/WMM.COF create mode 100644 ladybug/north.py create mode 100644 tests/assets/txt/WMM.COF create mode 100644 tests/north_test.py diff --git a/MANIFEST.in b/MANIFEST.in index 942a8a3e..cfd4ab6e 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,4 +1,5 @@ include ladybug/config.json +include ladybug/WMM.COF recursive-exclude tests * recursive-exclude docs * recursive-exclude .github * diff --git a/ladybug/WMM.COF b/ladybug/WMM.COF new file mode 100644 index 00000000..7940efb4 --- /dev/null +++ b/ladybug/WMM.COF @@ -0,0 +1,93 @@ + 2020.0 WMM-2020 12/10/2019 + 1 0 -29404.5 0.0 6.7 0.0 + 1 1 -1450.7 4652.9 7.7 -25.1 + 2 0 -2500.0 0.0 -11.5 0.0 + 2 1 2982.0 -2991.6 -7.1 -30.2 + 2 2 1676.8 -734.8 -2.2 -23.9 + 3 0 1363.9 0.0 2.8 0.0 + 3 1 -2381.0 -82.2 -6.2 5.7 + 3 2 1236.2 241.8 3.4 -1.0 + 3 3 525.7 -542.9 -12.2 1.1 + 4 0 903.1 0.0 -1.1 0.0 + 4 1 809.4 282.0 -1.6 0.2 + 4 2 86.2 -158.4 -6.0 6.9 + 4 3 -309.4 199.8 5.4 3.7 + 4 4 47.9 -350.1 -5.5 -5.6 + 5 0 -234.4 0.0 -0.3 0.0 + 5 1 363.1 47.7 0.6 0.1 + 5 2 187.8 208.4 -0.7 2.5 + 5 3 -140.7 -121.3 0.1 -0.9 + 5 4 -151.2 32.2 1.2 3.0 + 5 5 13.7 99.1 1.0 0.5 + 6 0 65.9 0.0 -0.6 0.0 + 6 1 65.6 -19.1 -0.4 0.1 + 6 2 73.0 25.0 0.5 -1.8 + 6 3 -121.5 52.7 1.4 -1.4 + 6 4 -36.2 -64.4 -1.4 0.9 + 6 5 13.5 9.0 -0.0 0.1 + 6 6 -64.7 68.1 0.8 1.0 + 7 0 80.6 0.0 -0.1 0.0 + 7 1 -76.8 -51.4 -0.3 0.5 + 7 2 -8.3 -16.8 -0.1 0.6 + 7 3 56.5 2.3 0.7 -0.7 + 7 4 15.8 23.5 0.2 -0.2 + 7 5 6.4 -2.2 -0.5 -1.2 + 7 6 -7.2 -27.2 -0.8 0.2 + 7 7 9.8 -1.9 1.0 0.3 + 8 0 23.6 0.0 -0.1 0.0 + 8 1 9.8 8.4 0.1 -0.3 + 8 2 -17.5 -15.3 -0.1 0.7 + 8 3 -0.4 12.8 0.5 -0.2 + 8 4 -21.1 -11.8 -0.1 0.5 + 8 5 15.3 14.9 0.4 -0.3 + 8 6 13.7 3.6 0.5 -0.5 + 8 7 -16.5 -6.9 0.0 0.4 + 8 8 -0.3 2.8 0.4 0.1 + 9 0 5.0 0.0 -0.1 0.0 + 9 1 8.2 -23.3 -0.2 -0.3 + 9 2 2.9 11.1 -0.0 0.2 + 9 3 -1.4 9.8 0.4 -0.4 + 9 4 -1.1 -5.1 -0.3 0.4 + 9 5 -13.3 -6.2 -0.0 0.1 + 9 6 1.1 7.8 0.3 -0.0 + 9 7 8.9 0.4 -0.0 -0.2 + 9 8 -9.3 -1.5 -0.0 0.5 + 9 9 -11.9 9.7 -0.4 0.2 + 10 0 -1.9 0.0 0.0 0.0 + 10 1 -6.2 3.4 -0.0 -0.0 + 10 2 -0.1 -0.2 -0.0 0.1 + 10 3 1.7 3.5 0.2 -0.3 + 10 4 -0.9 4.8 -0.1 0.1 + 10 5 0.6 -8.6 -0.2 -0.2 + 10 6 -0.9 -0.1 -0.0 0.1 + 10 7 1.9 -4.2 -0.1 -0.0 + 10 8 1.4 -3.4 -0.2 -0.1 + 10 9 -2.4 -0.1 -0.1 0.2 + 10 10 -3.9 -8.8 -0.0 -0.0 + 11 0 3.0 0.0 -0.0 0.0 + 11 1 -1.4 -0.0 -0.1 -0.0 + 11 2 -2.5 2.6 -0.0 0.1 + 11 3 2.4 -0.5 0.0 0.0 + 11 4 -0.9 -0.4 -0.0 0.2 + 11 5 0.3 0.6 -0.1 -0.0 + 11 6 -0.7 -0.2 0.0 0.0 + 11 7 -0.1 -1.7 -0.0 0.1 + 11 8 1.4 -1.6 -0.1 -0.0 + 11 9 -0.6 -3.0 -0.1 -0.1 + 11 10 0.2 -2.0 -0.1 0.0 + 11 11 3.1 -2.6 -0.1 -0.0 + 12 0 -2.0 0.0 0.0 0.0 + 12 1 -0.1 -1.2 -0.0 -0.0 + 12 2 0.5 0.5 -0.0 0.0 + 12 3 1.3 1.3 0.0 -0.1 + 12 4 -1.2 -1.8 -0.0 0.1 + 12 5 0.7 0.1 -0.0 -0.0 + 12 6 0.3 0.7 0.0 0.0 + 12 7 0.5 -0.1 -0.0 -0.0 + 12 8 -0.2 0.6 0.0 0.1 + 12 9 -0.5 0.2 -0.0 -0.0 + 12 10 0.1 -0.9 -0.0 -0.0 + 12 11 -1.1 -0.0 -0.0 0.0 + 12 12 -0.3 0.5 -0.1 -0.1 +999999999999999999999999999999999999999999999999 +999999999999999999999999999999999999999999999999 diff --git a/ladybug/north.py b/ladybug/north.py new file mode 100644 index 00000000..debcba3d --- /dev/null +++ b/ladybug/north.py @@ -0,0 +1,295 @@ +# coding=utf-8 +"""Module converting between magnetic and true North.""" +from __future__ import division + +import math +import os + + +class WorldMagneticModel(object): + """World Magnetic Model (WMM) that can convert from magnetic to true North. + + Args: + cof_file: The full path to a .COF file containing the coefficients that + form the inputs for the World Magnetic Model (WMM). A new set of coefficients + is published roughly every 5 years as the magnetic poles continue to + move. If None, coefficients will be derived from the WMM.COF file contained + within this Python package, which should be for the most recent model. + If not, the most recent coefficients are available at + https://www.ncei.noaa.gov/products/world-magnetic-model/wmm-coefficients + """ + + def __init__(self, cof_file=None): + """Initialize WorldMagneticModel.""" + # use the .COF file in this package if not specified + if cof_file is None: + cof_file = os.path.join(os.path.dirname(__file__), 'WMM.COF') + + # parse the coefficients from the file contents + wmm = [] + with open(cof_file) as wmm_file: + for line in wmm_file: + vals = line.strip().split() + if len(vals) == 3: + self.epoch = float(vals[0]) + self.model = vals[1] + self.modeldate = vals[2] + elif len(vals) == 6: + line_dict = { + 'n': int(float(vals[0])), + 'm': int(float(vals[1])), + 'gnm': float(vals[2]), + 'hnm': float(vals[3]), + 'dgnm': float(vals[4]), + 'dhnm': float(vals[5]) + } + wmm.append(line_dict) + + # set the global constants used by the WMM + z = [0.0] * 15 + self.maxord = 12 + self.tc = [] + for _ in range(14): + self.tc.append(z[0:13]) + self.sp = z[0:14] + self.cp = z[0:14] + self.cp[0] = 1.0 + self.pp = z[0:13] + self.pp[0] = 1.0 + self.p = [] + for _ in range(14): + self.p.append(z[0:14]) + self.p[0][0] = 1.0 + self.dp = [] + for _ in range(14): + self.dp.append(z[0:13]) + self.a = 6378.137 + self.b = 6356.7523142 + self.re = 6371.2 + self.a2 = self.a*self.a + self.b2 = self.b*self.b + self.c2 = self.a2-self.b2 + self.a4 = self.a2*self.a2 + self.b4 = self.b2*self.b2 + self.c4 = self.a4 - self.b4 + + self.c = [] + self.cd = [] + for _ in range(14): + self.c.append(z[0:14]) + self.cd.append(z[0:14]) + + # adjust C and CD from the file contents + for wmmnm in wmm: + m = wmmnm['m'] + n = wmmnm['n'] + gnm = wmmnm['gnm'] + hnm = wmmnm['hnm'] + dgnm = wmmnm['dgnm'] + dhnm = wmmnm['dhnm'] + if (m <= n): + self.c[m][n] = gnm + self.cd[m][n] = dgnm + if (m != 0): + self.c[n][m - 1] = hnm + self.cd[n][m - 1] = dhnm + + # convert schmidt normalized gauss coefficients to un-normalized + self.snorm = [] + for _ in range(13): + self.snorm.append(z[0:13]) + self.snorm[0][0] = 1.0 + self.k = [] + for _ in range(13): + self.k.append(z[0:13]) + self.k[1][1] = 0.0 + self.fn = [float(i) for i in range(14)] + self.fm = [float(i) for i in range(13)] + for n in range(1, self.maxord + 1): + self.snorm[0][n] = self.snorm[0][n - 1] * (2.0 * n - 1) / n + j = 2.0 + m = 0 + d_1 = 1 + d_2 = (n - m + d_1) / d_1 + while (d_2 > 0): + self.k[m][n] = (((n - 1) * (n - 1)) - (m * m)) / \ + ((2.0 * n - 1) * (2.0 * n - 3.0)) + if (m > 0): + flnmj = ((n - m + 1.0) * j) / (n + m) + self.snorm[m][n] = self.snorm[m - 1][n] * math.sqrt(flnmj) + j = 1.0 + self.c[n][m - 1] = self.snorm[m][n] * self.c[n][m - 1] + self.cd[n][m - 1] = self.snorm[m][n] * self.cd[n][m - 1] + self.c[m][n] = self.snorm[m][n] * self.c[m][n] + self.cd[m][n] = self.snorm[m][n] * self.cd[m][n] + d_2 = d_2 - 1 + m = m + d_1 + + # set default lat, lon and altitude, which will be overwritten as used + self.otime = -1000.0 + self.oalt = -1000.0 + self.olat = -1000.0 + self.olon = -1000.0 + + def magnetic_declination(self, latitude=0, longitude=0, elevation=0, year=2025): + """Compute the magnetic declination using the World Magnetic Model (WMM). + + Magnetic declination is the difference between magnetic North and true North at + a given location on the globe (expressed in terms of degrees). The function + here uses the same method that underlies the NOAA Magnetic Declination + calculator. + + Args: + latitude: A number between -90 and 90 for the latitude of the location + in degrees. (Default: 0 for the equator). + longitude: A number between -180 and 180 for the longitude of the location + in degrees (Default: 0 for the prime meridian). + elevation: A number for elevation of the location in meters. (Default: 0). + year: A number for the year in which the magnetic declination is + being evaluated. Decimal values are accepted. (Default: 2025). + + Returns: + A number for the magnetic declination in degrees. + """ + # set up the properties given the location and year information + alt = elevation / 1000 # to kilometers + dt = year - self.epoch + glat = latitude + glon = longitude + rlat = math.radians(glat) + rlon = math.radians(glon) + srlon = math.sin(rlon) + srlat = math.sin(rlat) + crlon = math.cos(rlon) + crlat = math.cos(rlat) + srlat2 = srlat * srlat + crlat2 = crlat * crlat + self.sp[1] = srlon + self.cp[1] = crlon + + # convert from geodetic to spherical coordinates + if (alt != self.oalt or glat != self.olat): + q = math.sqrt(self.a2 - self.c2 * srlat2) + q1 = alt * q + q2 = ((q1 + self.a2) / (q1 + self.b2)) * ((q1 + self.a2) / (q1 + self.b2)) + ct = srlat / math.sqrt(q2 * crlat2 + srlat2) + st = math.sqrt(1.0 - (ct * ct)) + r2 = (alt * alt) + 2.0 * q1 + (self.a4 - self.c4 * srlat2) / (q * q) + r = math.sqrt(r2) + d = math.sqrt(self.a2 * crlat2 + self.b2 * srlat2) + ca = (alt + d) / r + sa = self.c2 * crlat * srlat / (r * d) + + if (glon != self.olon): + for m in range(2, self.maxord + 1): + self.sp[m] = self.sp[1] * self.cp[m - 1] + self.cp[1] * self.sp[m - 1] + self.cp[m] = self.cp[1] * self.cp[m - 1] - self.sp[1] * self.sp[m - 1] + + aor = self.re / r + ar = aor * aor + br = bt = bp = bpp = 0.0 + for n in range(1, self.maxord + 1): + ar = ar * aor + m = 0 + D4 = (n + m + 1) + while D4 > 0: + # compute un-normalized associated polynomials and derivatives + if (alt != self.oalt or glat != self.olat): + if (n == m): + self.p[m][n] = st * self.p[m - 1][n - 1] + self.dp[m][n] = \ + st * self.dp[m - 1][n - 1] + ct * self.p[m - 1][n - 1] + + elif (n == 1 and m == 0): + self.p[m][n] = ct * self.p[m][n-1] + self.dp[m][n] = ct * self.dp[m][n - 1] - st * self.p[m][n - 1] + + elif (n > 1 and n != m): + if (m > n - 2): + self.p[m][n - 2] = 0 + if (m > n - 2): + self.dp[m][n - 2] = 0.0 + self.p[m][n] = \ + ct * self.p[m][n - 1] - self.k[m][n] * self.p[m][n - 2] + self.dp[m][n] = ct * self.dp[m][n - 1] - \ + st * self.p[m][n - 1]-self.k[m][n] * self.dp[m][n - 2] + + # time adjust the gauss coefficients + if (year != self.otime): + self.tc[m][n] = self.c[m][n] + dt * self.cd[m][n] + if (m != 0): + self.tc[n][m - 1] = self.c[n][m - 1] + dt * self.cd[n][m - 1] + + # accumulate terms of the spherical harmonic expansions + par = ar * self.p[m][n] + if (m == 0): + temp1 = self.tc[m][n] * self.cp[m] + temp2 = self.tc[m][n] * self.sp[m] + else: + temp1 = self.tc[m][n] * self.cp[m] + self.tc[n][m - 1] * self.sp[m] + temp2 = self.tc[m][n] * self.sp[m] - self.tc[n][m - 1] * self.cp[m] + + bt = bt - ar * temp1 * self.dp[m][n] + bp = bp + (self.fm[m] * temp2 * par) + br = br + (self.fn[n] * temp1 * par) + + # special case: north/south geographic poles + if (st == 0.0 and m == 1): + if (n == 1): + self.pp[n] = self.pp[n - 1] + else: + self.pp[n] = ct * self.pp[n - 1] - self.k[m][n] * self.pp[n - 2] + parp = ar * self.pp[n] + bpp = bpp + (self.fm[m] * temp2 * parp) + D4 = D4 - 1 + m = m + 1 + + if (st == 0.0): + bp = bpp + else: + bp = bp/st + + # rotate magnetic vector components from spherical to geodetic coordinates + bx = -bt * ca - br * sa + by = bp + declination = math.degrees(math.atan2(by, bx)) + + # set the location attributes that the model has been aligned to + self.otime = year + self.oalt = alt + self.olat = glat + self.olon = glon + + return declination + + def magnetic_to_true_north(self, location, magnetic_north=0, year=2025): + """Compute true North from a magnetic North vector. + + Args: + location: A Ladybug Location object that will be used to determine the + magnetic declination. + magnetic_north: A number between -360 and 360 for the counterclockwise + difference between the North and the positive Y-axis in degrees. + 90 is West and 270 is East (Default: 0). + year: A number for the year in which the magnetic declination is + being evaluated. Decimal values are accepted. (Default: 2025). + + Returns: + A number between -360 and 360 for the true North angle in degrees. + """ + declination = self.magnetic_declination( + location.latitude, location.longitude, location.elevation, year) + true_north = magnetic_north + declination + if true_north > 360: + true_north = true_north - 360 + elif true_north < -360: + true_north = 360 + true_north + return true_north + + def ToString(self): + """Overwrite .NET ToString.""" + return self.__repr__() + + def __repr__(self): + """Return WorldMagneticModel as a string.""" + return 'WorldMagneticModel: {}'.format(self.epoch) diff --git a/tests/assets/txt/WMM.COF b/tests/assets/txt/WMM.COF new file mode 100644 index 00000000..cd86b31a --- /dev/null +++ b/tests/assets/txt/WMM.COF @@ -0,0 +1,93 @@ + 2015.0 WMM-2015 12/15/2014 + 1 0 -29438.5 0.0 10.7 0.0 + 1 1 -1501.1 4796.2 17.9 -26.8 + 2 0 -2445.3 0.0 -8.6 0.0 + 2 1 3012.5 -2845.6 -3.3 -27.1 + 2 2 1676.6 -642.0 2.4 -13.3 + 3 0 1351.1 0.0 3.1 0.0 + 3 1 -2352.3 -115.3 -6.2 8.4 + 3 2 1225.6 245.0 -0.4 -0.4 + 3 3 581.9 -538.3 -10.4 2.3 + 4 0 907.2 0.0 -0.4 0.0 + 4 1 813.7 283.4 0.8 -0.6 + 4 2 120.3 -188.6 -9.2 5.3 + 4 3 -335.0 180.9 4.0 3.0 + 4 4 70.3 -329.5 -4.2 -5.3 + 5 0 -232.6 0.0 -0.2 0.0 + 5 1 360.1 47.4 0.1 0.4 + 5 2 192.4 196.9 -1.4 1.6 + 5 3 -141.0 -119.4 0.0 -1.1 + 5 4 -157.4 16.1 1.3 3.3 + 5 5 4.3 100.1 3.8 0.1 + 6 0 69.5 0.0 -0.5 0.0 + 6 1 67.4 -20.7 -0.2 0.0 + 6 2 72.8 33.2 -0.6 -2.2 + 6 3 -129.8 58.8 2.4 -0.7 + 6 4 -29.0 -66.5 -1.1 0.1 + 6 5 13.2 7.3 0.3 1.0 + 6 6 -70.9 62.5 1.5 1.3 + 7 0 81.6 0.0 0.2 0.0 + 7 1 -76.1 -54.1 -0.2 0.7 + 7 2 -6.8 -19.4 -0.4 0.5 + 7 3 51.9 5.6 1.3 -0.2 + 7 4 15.0 24.4 0.2 -0.1 + 7 5 9.3 3.3 -0.4 -0.7 + 7 6 -2.8 -27.5 -0.9 0.1 + 7 7 6.7 -2.3 0.3 0.1 + 8 0 24.0 0.0 0.0 0.0 + 8 1 8.6 10.2 0.1 -0.3 + 8 2 -16.9 -18.1 -0.5 0.3 + 8 3 -3.2 13.2 0.5 0.3 + 8 4 -20.6 -14.6 -0.2 0.6 + 8 5 13.3 16.2 0.4 -0.1 + 8 6 11.7 5.7 0.2 -0.2 + 8 7 -16.0 -9.1 -0.4 0.3 + 8 8 -2.0 2.2 0.3 0.0 + 9 0 5.4 0.0 0.0 0.0 + 9 1 8.8 -21.6 -0.1 -0.2 + 9 2 3.1 10.8 -0.1 -0.1 + 9 3 -3.1 11.7 0.4 -0.2 + 9 4 0.6 -6.8 -0.5 0.1 + 9 5 -13.3 -6.9 -0.2 0.1 + 9 6 -0.1 7.8 0.1 0.0 + 9 7 8.7 1.0 0.0 -0.2 + 9 8 -9.1 -3.9 -0.2 0.4 + 9 9 -10.5 8.5 -0.1 0.3 + 10 0 -1.9 0.0 0.0 0.0 + 10 1 -6.5 3.3 0.0 0.1 + 10 2 0.2 -0.3 -0.1 -0.1 + 10 3 0.6 4.6 0.3 0.0 + 10 4 -0.6 4.4 -0.1 0.0 + 10 5 1.7 -7.9 -0.1 -0.2 + 10 6 -0.7 -0.6 -0.1 0.1 + 10 7 2.1 -4.1 0.0 -0.1 + 10 8 2.3 -2.8 -0.2 -0.2 + 10 9 -1.8 -1.1 -0.1 0.1 + 10 10 -3.6 -8.7 -0.2 -0.1 + 11 0 3.1 0.0 0.0 0.0 + 11 1 -1.5 -0.1 0.0 0.0 + 11 2 -2.3 2.1 -0.1 0.1 + 11 3 2.1 -0.7 0.1 0.0 + 11 4 -0.9 -1.1 0.0 0.1 + 11 5 0.6 0.7 0.0 0.0 + 11 6 -0.7 -0.2 0.0 0.0 + 11 7 0.2 -2.1 0.0 0.1 + 11 8 1.7 -1.5 0.0 0.0 + 11 9 -0.2 -2.5 0.0 -0.1 + 11 10 0.4 -2.0 -0.1 0.0 + 11 11 3.5 -2.3 -0.1 -0.1 + 12 0 -2.0 0.0 0.1 0.0 + 12 1 -0.3 -1.0 0.0 0.0 + 12 2 0.4 0.5 0.0 0.0 + 12 3 1.3 1.8 0.1 -0.1 + 12 4 -0.9 -2.2 -0.1 0.0 + 12 5 0.9 0.3 0.0 0.0 + 12 6 0.1 0.7 0.1 0.0 + 12 7 0.5 -0.1 0.0 0.0 + 12 8 -0.4 0.3 0.0 0.0 + 12 9 -0.4 0.2 0.0 0.0 + 12 10 0.2 -0.9 0.0 0.0 + 12 11 -0.9 -0.2 0.0 0.0 + 12 12 0.0 0.7 0.0 0.0 +999999999999999999999999999999999999999999999999 +999999999999999999999999999999999999999999999999 diff --git a/tests/north_test.py b/tests/north_test.py new file mode 100644 index 00000000..b93c42e1 --- /dev/null +++ b/tests/north_test.py @@ -0,0 +1,47 @@ +# coding=utf-8 +from __future__ import division + +import pytest +import datetime + +from ladybug.north import WorldMagneticModel + + +def test_world_magnetic_model_defaults(): + """Test the defaults of the WorldMagneticModel.""" + wmm_obj = WorldMagneticModel() + dec_val = wmm_obj.magnetic_declination(80, 0, 0, 2015) + + assert dec_val == pytest.approx(-3.7947, rel=1e-2) + + +def test_world_magnetic_model_results(): + """Test that the WorldMagneticModel returns correct results.""" + d1 = datetime.date(2015, 1, 1) + d2 = datetime.date(2017, 7, 2) + y1 = d1.year + ((d1 - datetime.date(d1.year, 1, 1)).days / 365.0) + y2 = d2.year + ((d2 - datetime.date(d2.year, 1, 1)).days / 365.0) + + test_values = ( + # date, alt, lat, lon, var + (y1, 0, 80, 0, -3.85), + (y1, 0, 0, 120, 0.57), + (y1, 0, -80, 240, 69.81), + (y1, 100000.0, 80, 0, -4.27), + (y1, 100000.0, 0, 120, 0.56), + (y1, 100000.0, -80, 240, 69.22), + (y2, 0, 80, 0, -2.75), + (y2, 0, 0, 120, 0.32), + (y2, 0, -80, 240, 69.58), + (y2, 100000.0, 80, 0, -3.17), + (y2, 100000.0, 0, 120, 0.32), + (y2, 100000.0, -80, 240, 69.00), + ) + + wmm_2015 = './tests/assets/txt/WMM.COF' + wmm_obj = WorldMagneticModel(wmm_2015) + + for values in test_values: + dec_val = wmm_obj.magnetic_declination( + values[2], values[3], values[1], values[0]) + assert dec_val == pytest.approx(values[4], rel=1e-2)