def coord_list_mapping_pbc(subset, superset, atol=1e-8):
    """
    Gives the index mapping from a subset to a superset.
    Subset and superset cannot contain duplicate rows

    Args:
        subset, superset: List of frac_coords

    Returns:
        list of indices such that superset[indices] = subset
    """
    c1 = np.array(subset)
    c2 = np.array(superset)

    diff = c1[:, None, :] - c2[None, :, :]
    diff -= np.round(diff)
    inds = np.where(np.all(np.abs(diff) < atol, axis=2))[1]

    # verify result (its easier to check validity of the result than
    # the validity of inputs)
    test = c2[inds] - c1
    test -= np.round(test)
    if not np.allclose(test, 0):
        if not is_coord_subset_pbc(subset, superset):
            raise ValueError("subset is not a subset of superset")
    if not test.shape == c1.shape:
        raise ValueError("Something wrong with the inputs, likely duplicates " "in superset")
    return inds

def ReadBPLASMA(file_name,BNORM,Ns):
    #Read the BPLASMA output file from MARS-F
    #Return BM1, BM2, BM3
    BPLASMA = num.loadtxt(open(file_name))
 
    Nm1 = BPLASMA[0,0]
    n = num.round(BPLASMA[0,2])
    Mm = num.round(BPLASMA[1:Nm1+1,0])
    Mm.resize([len(Mm),1])


    BM1 = BPLASMA[Nm1+1:,0] + BPLASMA[Nm1+1:,1]*1j
    BM2 = BPLASMA[Nm1+1:,2] + BPLASMA[Nm1+1:,3]*1j
    BM3 = BPLASMA[Nm1+1:,4] + BPLASMA[Nm1+1:,5]*1j

    BM1 = num.reshape(BM1,[Ns,Nm1],order='F')
    BM2 = num.reshape(BM2,[Ns,Nm1],order='F')
    BM3 = num.reshape(BM3,[Ns,Nm1],order='F')

    BM1 = BM1[0:Ns,:]*BNORM
    BM2 = BM2[0:Ns,:]*BNORM
    BM3 = BM3[0:Ns,:]*BNORM

    #NEED TO KNOW WHY THIS SECTION IS INCLUDED - to do with half grid???!!
    #BM2[1:,:] = BM2[0:-1,:] Needed to comment out to compare with RZPlot3
    #BM3[1:,:] = BM3[0:-1,:]

    return BM1, BM2, BM3,Mm

def bandpass_filter(files, lowpass_freq, highpass_freq, fs):
    """Bandpass filter the input files

    Parameters
    ----------
    files: list of 4d nifti files
    lowpass_freq: cutoff frequency for the low pass filter (in Hz)
    highpass_freq: cutoff frequency for the high pass filter (in Hz)
    fs: sampling rate (in Hz)
    """
    out_files = []
    for filename in filename_to_list(files):
        path, name, ext = split_filename(filename)
        out_file = os.path.join(os.getcwd(), name + '_bp' + ext)
        img = nb.load(filename)
        timepoints = img.shape[-1]
        F = np.zeros((timepoints))
        lowidx = int(timepoints / 2) + 1
        if lowpass_freq > 0:
            lowidx = np.round(float(lowpass_freq) / fs * timepoints)
        highidx = 0
        if highpass_freq > 0:
            highidx = np.round(float(highpass_freq) / fs * timepoints)
        F[highidx:lowidx] = 1
        F = ((F + F[::-1]) > 0).astype(int)
        data = img.get_data()
        if np.all(F == 1):
            filtered_data = data
        else:
            filtered_data = np.real(np.fft.ifftn(np.fft.fftn(data) * F))
        img_out = nb.Nifti1Image(filtered_data, img.affine, img.header)
        img_out.to_filename(out_file)
        out_files.append(out_file)
    return list_to_filename(out_files)

def _af_majo(age, smic55, af_nbenf, _P, _option={'age': ENFS, 'smic55': ENFS}):
    '''
    Allocations familiales - majoration pour âge
    'fam'
    '''
    # TODO: Date d'entrée en vigueur de la nouvelle majoration
    # enfants nés après le "1997-04-30"
    bmaf = _P.fam.af.bmaf
    P_af = _P.fam.af
    P = _P.fam.af.maj_age
    af_maj1 = round(bmaf * P.taux1, 2)
    af_maj2 = round(bmaf * P.taux2, 2)

    ageaine = age_aine(age, smic55, P_af.age1, P_af.age2)

    def age_sf_aine(age, ag1, ag2, ageaine):
        dum = (ag1 <= ageaine) & (ageaine <= ag2)
        return nb_enf(age, smic55, ag1, ag2) - dum * 1


    nbenf_maj1 = ( (af_nbenf == 2)*age_sf_aine(age, P.age1, P.age2 - 1, ageaine)
                   + nb_enf(age, smic55, P.age1, P.age2 - 1)*(af_nbenf >= 3)  )
    nbenf_maj2 = ( (af_nbenf == 2)*age_sf_aine(age, P.age2, P_af.age2, ageaine)
                   + nb_enf(age, smic55, P.age2, P_af.age2)*(af_nbenf >= 3)  )

    af_majo = nbenf_maj1 * af_maj1 + nbenf_maj2 * af_maj2

    return 12*af_majo # annualisé

def _apje(br_pf, age, smic55, isol, biact, _P, _option={'age': ENFS, 'smic55': ENFS}):
    '''
    Allocation pour jeune enfant
    '''
    # TODO: APJE courte voir doc ERF 2006
    P = _P.fam
    nbenf = nb_enf(age, smic55, 0, P.apje.age - 1)
    bmaf = P.af.bmaf
    bmaf_n_2 = P.af.bmaf_n_2
    base = round(P.apje.taux * bmaf, 2)
    base2 = round(P.apje.taux * bmaf_n_2, 2)

    plaf_tx = (nbenf > 0) + P.apje.plaf_tx1 * min_(nbenf, 2) + P.apje.plaf_tx2 * max_(nbenf - 2, 0)
    majo = isol | biact
    plaf = P.apje.plaf * plaf_tx + P.apje.plaf_maj * majo
    plaf2 = plaf + 12 * base2

    apje = (nbenf >= 1) * ((br_pf <= plaf) * base
                            + (br_pf > plaf) * max_(plaf2 - br_pf, 0) / 12.0)

    # Pour bénéficier de cette allocation, il faut que tous les enfants du foyer soient nés, adoptés, ou recueillis en vue d’une adoption avant le 1er janvier 2004, et qu’au moins l’un d’entre eux ait moins de 3 ans.
    # Cette allocation est verséE du 5ème mois de grossesse jusqu’au mois précédant le 3ème anniversaire de l’enfant.

    # Non cumul APE APJE CF
    #  - L’allocation parentale d’éducation (APE), sauf pour les femmes enceintes.
    #    L’APJE est alors versée du 5ème mois de grossesse jusqu’à la naissance de l’enfant.
    #  - Le CF
    return 12*apje  # annualisé

	def on_draw1(self):
		"""
		Устанавливается режим тридэ. Очищается буфер. Загружается единичная матрица.
		Поворачивается всё по оси x (направление вверх/вниз), затем поворачивается
		по оси y (влево/вправо). Затем вызывается функция опреления высоты игрока
		над уровнем моря. И после этого мир двигается куда нужно.
		Порядок матриц ВАЖЕН. ОЧЕНЬ.
		Затем выбирается материал ящиков, рисуются ящики, выбирается материал
		стен комнаты, рисуется комната.
		Перенастраивается в 2д, пишется фпс.
		"""
		self.setup3d()
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
		glLoadIdentity()
		glRotatef(self.player.xrot,1,0,0)
		glRotatef(self.player.yrot,0,1,0)
		glLightfv(GL_LIGHT1, GL_POSITION, vec(0,1,-1))
		glTranslatef(-self.player.xpos,-self.player.ypos-self.player.height,self.player.zpos)
# 		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, vec(1,0.3,0,0))
# 		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, vec(0,0,0,0))
		glPointSize(20)
		glColor3f(0.5,0,0.2)
		glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, vec(0.75,0.52,0.15,0))
		self.batch.draw()
		glColor3f(0,0,0)
		glEnable(GL_TEXTURE_2D)
		glBindTexture(self.boxtexture.target,self.boxtexture.id)
		self.batch_box.draw()
		glDisable(GL_TEXTURE_2D)
		glColor3f(0.5,0,0.2)
		model.draw()
# 		glTranslatef(0,-20,0)
# 		robot[self.robot_fr].draw()
# 		glTranslatef(-5,0,0)
# 		robot[self.robot_fr].draw()
# 		glTranslatef(10,0,0)
# 		robot[self.robot_fr].draw()
# 		glTranslatef(-15,0,0)
# 		robot[self.robot_fr].draw()
# 		glTranslatef(20,0,0)
# 		robot[self.robot_fr].draw()
# 		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, vec(0,0.3,0,0))
		self.setup2d()
			# лейблы с координатами. Округлено до трёх знаков после запятой, всё в одну строку чтобы показывалось.
		pltxt = "X: "+str(np.round(self.player.xpos,3))+" Y: "+str(np.round(self.player.ypos,3))+" Z: "+str(np.round(self.player.zpos,3))
		self.coords.text = pltxt
		self.coords.y = self.height-30
		plrottxt = "Xrot: "+str(np.round(self.player.xrot,3))+" Yrot: "+str(np.round(self.player.yrot,3))+" Zrot: "+str(np.round(self.player.zrot,3))
		self.rot.text = plrottxt
		self.rot.y = self.height-50
		self.pljfw.text = self.player.jumping
		self.pljfw.y = self.height-70
		self.times.text = str(np.round(self.time,3))
		self.times.y = self.height-90
		self.times.draw()
		self.pljfw.draw()
		self.coords.draw()
		self.rot.draw()
		self.fps.draw()
		self.picspr.draw()

def test_maskandscale():
    t = np.linspace(20, 30, 15)
    t[3] = 100
    tm = np.ma.masked_greater(t, 99)
    fname = pjoin(TEST_DATA_PATH, 'example_2.nc')
    with netcdf_file(fname, maskandscale=True) as f:
        Temp = f.variables['Temperature']
        assert_equal(Temp.missing_value, 9999)
        assert_equal(Temp.add_offset, 20)
        assert_equal(Temp.scale_factor, np.float32(0.01))
        found = Temp[:].compressed()
        del Temp  # Remove ref to mmap, so file can be closed.
        expected = np.round(tm.compressed(), 2)
        assert_allclose(found, expected)

    with in_tempdir():
        newfname = 'ms.nc'
        f = netcdf_file(newfname, 'w', maskandscale=True)
        f.createDimension('Temperature', len(tm))
        temp = f.createVariable('Temperature', 'i', ('Temperature',))
        temp.missing_value = 9999
        temp.scale_factor = 0.01
        temp.add_offset = 20
        temp[:] = tm
        f.close()

        with netcdf_file(newfname, maskandscale=True) as f:
            Temp = f.variables['Temperature']
            assert_equal(Temp.missing_value, 9999)
            assert_equal(Temp.add_offset, 20)
            assert_equal(Temp.scale_factor, np.float32(0.01))
            expected = np.round(tm.compressed(), 2)
            found = Temp[:].compressed()
            del Temp
            assert_allclose(found, expected)

    def _to_dense(self):
        """ Convert the sparse [onset, duration, amplitude] representation
        typical of event files to a dense matrix where each row represents
        a fixed unit of time. """
        end = int((self.events['onset'] + self.events['duration']).max())

        targ_hz, orig_hz = self.target_hz, self.orig_hz
        len_ts = end * targ_hz
        conditions = self.events['condition'].unique().tolist()
        n_conditions = len(conditions)
        ts = np.zeros((len_ts, n_conditions))

        _events = self.events.copy().reset_index()
        _events[['onset', 'duration']] = \
            _events[['onset', 'duration']] * targ_hz / orig_hz

        cond_index = [conditions.index(c) for c in _events['condition']]
        ev_end = np.round(_events['onset'] + _events['duration']).astype(int)
        onsets = np.round(_events['onset']).astype(int)

        for i, row in _events.iterrows():
            ts[onsets[i]:ev_end[i], cond_index[i]] = row['amplitude']

        self.data = pd.DataFrame(ts, columns=conditions)
        onsets = np.arange(len(ts)) / self.target_hz
        self.data.insert(0, 'onset', onsets)

    def imageCoCenter(self, inst, algo):

        x1, y1, tmp = getCenterAndR_ef(self.image)
        if algo.debugLevel >= 3:
            print('imageCoCenter: (x1,y1)=(%8.2f,%8.2f)\n' % (x1, y1))

        stampCenterx1 = inst.sensorSamples / 2. + 0.5
        stampCentery1 = inst.sensorSamples / 2. + 0.5
        radialShift = 3.5 * algo.upReso * \
            (inst.offset / 1e-3) * (10e-6 / inst.pixelSize)

        radialShift = radialShift * self.fldr / 1.75
        if (self.fldr > 1.75):
            radialShift = 0

        if self.fldr != 0:
            I1c = self.fieldX / self.fldr
            I1s = self.fieldY / self.fldr
        else:
            I1c = 0
            I1s = 0

        stampCenterx1 = stampCenterx1 + radialShift * I1c
        stampCentery1 = stampCentery1 + radialShift * I1s

        self.image = np.roll(self.image, int(
            np.round(stampCentery1 - y1)), axis=0)
        self.image = np.roll(self.image, int(
            np.round(stampCenterx1 - x1)), axis=1)

def test_D_infinity_flat_closed_upper():
    mg = RasterModelGrid((5, 4), xy_spacing=(1, 1))
    z = mg.add_zeros("node", "topographic__elevation")
    z[mg.core_nodes] -= 1
    mg.set_closed_boundaries_at_grid_edges(
        bottom_is_closed=True,
        left_is_closed=True,
        right_is_closed=True,
        top_is_closed=True,
    )

    fd = FlowDirectorDINF(mg)
    fd.run_one_step()

    node_ids = np.arange(mg.number_of_nodes)
    true_recievers = -1 * np.ones(fd.receivers.shape)
    true_recievers[:, 0] = node_ids

    true_proportions = np.zeros(fd.proportions.shape)
    true_proportions[:, 0] = 1

    assert_array_equal(fd.receivers, true_recievers)
    assert_array_equal(
        np.round(fd.proportions, decimals=6), np.round(true_proportions, decimals=6)
    )

def query_form(filename="merged_table.ipac"):

    table = Table.read(os.path.join(app.config['DATABASE_FOLDER'], filename), format='ascii.ipac')
    
    tolerance=1.1

    min_values=[np.round(min(table['SurfaceDensity'])/tolerance,4),np.round(min(table['VelocityDispersion'])/tolerance,4),np.round(min(table['Radius'])/tolerance,4)]
    max_values=[np.round(max(table['SurfaceDensity'])*tolerance,1),np.round(max(table['VelocityDispersion'])*tolerance,1),np.round(max(table['Radius'])*tolerance,1)]

    usetable = table[use_column_names]

    best_matches = {difflib.get_close_matches(vcn, usetable.colnames,  n=1,
                                              cutoff=0.4)[0]: vcn
                    for vcn in use_column_names
                    if any(difflib.get_close_matches(vcn, usetable.colnames, n=1, cutoff=0.4))
                   }

    best_column_names = [best_matches[colname] if colname in best_matches else 'Ignore'
                         for colname in usetable.colnames]

    return render_template("query_form.html", table=table, usetable=usetable,
                           use_units=use_units, filename=filename,
                           use_column_names=use_column_names,
                           best_column_names=best_column_names,
                           min_values=min_values,
                           max_values=max_values
                          )

 def round_values(self):
     """
     PURPOSE:
     To round the input parameters for lookup in the Claret table.
     """
     self.rteff = np.round(self.teff/250) * 250
     self.rlogg = np.round(self.logg/0.5) * 0.5

    def function(self, simulation, period):
        period = period.start.offset('first-of', 'month').period('month')
        age_holder = simulation.compute('age', period)
        smic55_holder = simulation.compute('smic55', period)
        af_nbenf = simulation.calculate('af_nbenf', period)
        P = simulation.legislation_at(period.start).fam.af

        age = self.split_by_roles(age_holder, roles = ENFS)
        smic55 = self.split_by_roles(smic55_holder, roles = ENFS)
        # TODO: Date d'entrée en vigueur de la nouvelle majoration
        # enfants nés après le "1997-04-30"
        bmaf = P.bmaf
        P_maj = P.maj_age
        af_maj1 = round(bmaf * P_maj.taux1, 2)
        af_maj2 = round(bmaf * P_maj.taux2, 2)
        ageaine = age_aine(age, smic55, P.age1, P.age2)

        def age_sf_aine(age, ag1, ag2, ageaine):
            dum = (ag1 <= ageaine) & (ageaine <= ag2)
            return nb_enf(age, smic55, ag1, ag2) - dum * 1

        nbenf_maj1 = (
            (af_nbenf == 2) * age_sf_aine(age, P_maj.age1, P_maj.age2 - 1, ageaine) +
            nb_enf(age, smic55, P_maj.age1, P_maj.age2 - 1) * (af_nbenf >= 3)
            )
        nbenf_maj2 = (
            (af_nbenf == 2) * age_sf_aine(age, P_maj.age2, P.age2, ageaine) +
            nb_enf(age, smic55, P_maj.age2, P.age2) * (af_nbenf >= 3)
            )
        return period, nbenf_maj1 * af_maj1 + nbenf_maj2 * af_maj2

def _as_timedelta64_scalar(time, unit=None):
    unit_args = [unit] if unit else []
    flt_unit = unit if unit else 's'
    # turn 'H:M:S.ms', 'M:S.ms', 'S.ms' into floating point seconds
    if isinstance(time, string_types):# and ':' in time:
        time = [float(t) for t in time.lstrip('T').split(':')][::-1]
        if len(time) > 1 and unit is not None:
            raise ValueError("When giving time as a string, units are automatic")
        if len(time) > 3:
            raise ValueError("Timedelta as string only goes up to hours")
        t_flt = 0.0
        for factor, t in zip([1, 60, 60 * 60], time):
            t_flt += factor * t
        time = t_flt
        flt_unit = 's'
    # turn floating point time into integer with the correct unit
    if is_datetime_like(time):
        time = as_datetime64(time) - as_datetime64(np.timedelta64(0, 's'))
    elif isinstance(time, (np.timedelta64, timedelta)):
        time = np.timedelta64(time).astype(_format_unit(unit, base=DELTA_BASE))
    elif isinstance(time, (int, float, np.integer, np.floating)):
        orig_time, orig_flt_unit = time, flt_unit
        unit_idx = TIME_UNITS.index(flt_unit)
        while not np.isclose(time, int(np.round(time)), rtol=1e-4, atol=1e-18):
            if unit_idx <= 0:
                raise ValueError("Floating point time {0} [{1}] is too precise "
                                 "for any time unit?".format(orig_time, orig_flt_unit))
            unit_idx -= 1
            time *= TIME_SCALE[unit_idx]
            flt_unit = TIME_UNITS[unit_idx]
        time = np.timedelta64(int(np.round(time)), flt_unit)
        unit, unit_args = flt_unit, [flt_unit]
    return np.timedelta64(time, *unit_args)

def get_histogram(series, bins, bins_decimals=0, bins_is_percent=False,
                  block_count=100):
    """Creates a text-based histogram.

    Args:
        series: pandas.Series of numeric values.
        bins: List of boundaries between bins in ascending order.
        bins_decimals: Number of decimals to use for bins in format string.
        bins_is_percent: Whether to print a '%' character for bins.
        block_count: Total number of block characters in histogram.
    """
    histogram = ''
    buckets = series.groupby(pd.cut(series, bins)).count()
    scaled_bins = 100 * bins if bins_is_percent else bins

    # Find the max string length for an individual bin value so that right
    # alignment works properly.
    max_bin_value_len = len(str(int(np.round(max(abs(scaled_bins)))))) + (
        (bins_decimals + 1) if bins_decimals > 0 else 0) + (
            1 if min(scaled_bins) < 0 else 0)

    format_str = '  '.join(['{:' + str(max_bin_value_len) + '.' + str(
        bins_decimals) + ('f}%' if bins_is_percent else 'f}')] * 2) + (
            '  {:<' + str(len(str(buckets.max()))) + '}  {}\n')
    for i in range(buckets.size):
        # Due to rounding exact number of blocks may vary.
        histogram += format_str.format(
            scaled_bins[i],
            scaled_bins[i + 1],
            buckets[i],
            ''.join(['*'] * np.round(
                block_count * buckets[i] / series.size)))
    return histogram

 def _scaling_fun(arr, scale):
     """
     scales and rounds -- does it all in place.
     """
     arr *= scale
     np.round(arr, out=arr)
     return arr

def fitting(d0, d1):
    idx_list=[]
    pos_list=[]
    for tp in['beta', 'sw']:
        e=0
        for net, sl in zip(['Net_0', 'Net_1'], [slice(2,4), slice(0,4)]):
            z=d0[tp][net]['mean_rates'][:,sl]
            target=d1[tp][net]['mean_rates'][sl]
            target=numpy.array([target]*z.shape[0])
            if e is 0:
                e=z-target
            else:
                e=numpy.concatenate((z-target, e),axis=1)
        
        e**=2
        e=numpy.sqrt(numpy.mean(e, axis=1))
        idx=numpy.argsort(e)
#         idx_list.append(idx)
#         l=[]
#         for i, _id in enumerate(idx_list[-2]):
#             j=list(idx_list[-1]).index(_id)
#             l.append([i,j])
#         l=numpy.array(l)
#         pos_list.append(l)
#         e=numpy.mean(l,axis=1)
        idx=numpy.argsort(e)
        
#         pp(list(l[idx,:]))
            
#         print idx
#         print e[idx]
        print tp
        for _id in idx[:100]:
            print d0[tp]['Net_0']['ylabels'][_id], d1[tp]['Net_0']['mean_rates'][:], numpy.round(d0[tp]['Net_0']['mean_rates'][_id,:],1),e[_id]
            print d0[tp]['Net_1']['ylabels'][_id], d1[tp]['Net_1']['mean_rates'][:], numpy.round(d0[tp]['Net_1']['mean_rates'][_id,:],1)

        def on_press(event):

            try:
                self.lasttime = self.curtime
            except AttributeError:
                self.lasttime = 0

            self.curtime = time.time()
            if np.round(self.curtime, 2) != np.round(self.lasttime, 2):
                tdiff = self.curtime - self.lasttime
                print "time diff", tdiff
                print "you pressed", event.button, event.xdata, event.ydata
                if tdiff < 0.25:
                    tdiff = 0
                    x1, x2 = self.axes.get_xlim()
                    y1, y2 = self.axes.get_ylim()
                    print x1, x2, y1, y2
                    # ~ self.axes.set_xlim(x1/2,x2/2)
                    # ~ self.axes.set_ylim(y1/2,y2/2)
                    # self.axes.set_xlim(event.xdata-(x2-x1)/2, event.xdata+(x2-x1)/2)
                    if event.button == 1:  # zoom in
                        self.axes.set_xlim(event.xdata - (x2 - x1) / 4, event.xdata + (x2 - x1) / 4)
                        self.axes.set_ylim(event.ydata - (y2 - y1) / 4, event.ydata + (y2 - y1) / 4)
                    if event.button == 3:  # zoom out
                        self.axes.set_xlim(event.xdata - (x2 - x1) * 1, event.xdata + (x2 - x1) * 1)
                        self.axes.set_ylim(event.ydata - (y2 - y1) * 1, event.ydata + (y2 - y1) * 1)

                    # self.axes.set_xlim(np.mean(event.xdata,x1), np.mean(event.xdata,x2))
                    # self.axes.set_ylim(np.mean(event.ydata,y1), np.mean(event.ydata,y2))

                    self.canvas.draw()

def get_facet_values(facet, ra, dec, root="facet", default=0):
    """
    Extract the value from a fits facet file
    """
    import numpy as np
    from astropy.io import fits
    from astropy.wcs import WCS

    # TODO: Check astropy version
    # TODO: Check facet is a fits file

    with fits.open(facet) as f:
        shape = f[0].data.shape

        w = WCS(f[0].header)
        freq = w.wcs.crval[2]
        stokes = w.wcs.crval[3]

        xe, ye, _1, _2 = w.all_world2pix(ra, dec, freq, stokes, 1)
        x, y = np.round(xe).astype(int), np.round(ye).astype(int)

        # Dummy value for points out of the fits area
        x[(x < 0) | (x >= shape[-1])] = -1
        y[(y < 0) | (y >= shape[-2])] = -1

        data = f[0].data[0,0,:,:]

        values = data[y, x]

        # Assign the default value to NaNs and points out of the fits area
        values[(x == -1) | (y == -1)] = default
        values[np.isnan(values)] = default

        #TODO: Flexible format for other data types ?
        return np.array(["{}_{:.0f}".format(root, val) for val in values])

def Mie_ab(m,x):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#Mie_ab
  mx = m*x
  nmax = np.round(2+x+4*(x**(1/3)))
  nmx = np.round(max(nmax,np.abs(mx))+16)
  n = np.arange(1,nmax+1)
  nu = n + 0.5

  sx = np.sqrt(0.5*np.pi*x)
  px = sx*jv(nu,x)

  p1x = np.append(np.sin(x), px[0:int(nmax)-1])
  chx = -sx*yv(nu,x)

  ch1x = np.append(np.cos(x), chx[0:int(nmax)-1])
  gsx = px-(0+1j)*chx
  gs1x = p1x-(0+1j)*ch1x

  # B&H Equation 4.89
  Dn = np.zeros(int(nmx),dtype=complex)
  for i in range(int(nmx)-1,1,-1):
    Dn[i-1] = (i/mx)-(1/(Dn[i]+i/mx))

  D = Dn[1:int(nmax)+1] # Dn(mx), drop terms beyond nMax
  da = D/m+n/x
  db = m*D+n/x

  an = (da*px-p1x)/(da*gsx-gs1x)
  bn = (db*px-p1x)/(db*gsx-gs1x)

  return an, bn