How to use the pyuvdata.utils._fits_gethduaxis function in pyuvdata

self.Nspws = 1
                self.spw_array = np.array([0])

                # the axis number for phase center depends on if the spw exists
                self.phase_center_ra_degrees = np.float(vis_hdr.pop("CRVAL5"))
                self.phase_center_dec_degrees = np.float(vis_hdr.pop("CRVAL6"))

            # get shapes
            self.Nfreqs = vis_hdr.pop("NAXIS4")
            self.Npols = vis_hdr.pop("NAXIS3")
            self.Nblts = vis_hdr.pop("GCOUNT")

            self.freq_array = uvutils._fits_gethduaxis(vis_hdu, 4)
            self.freq_array.shape = (self.Nspws,) + self.freq_array.shape
            self.channel_width = vis_hdr.pop("CDELT4")
            self.polarization_array = np.int32(uvutils._fits_gethduaxis(vis_hdu, 3))

            # other info -- not required but frequently used
            self.object_name = vis_hdr.pop("OBJECT", None)
            self.telescope_name = vis_hdr.pop("TELESCOP", None)
            self.instrument = vis_hdr.pop("INSTRUME", None)
            latitude_degrees = vis_hdr.pop("LAT", None)
            longitude_degrees = vis_hdr.pop("LON", None)
            altitude = vis_hdr.pop("ALT", None)
            self.x_orientation = vis_hdr.pop("XORIENT", None)
            blt_order_str = vis_hdr.pop("BLTORDER", None)
            if blt_order_str is not None:
                self.blt_order = tuple(blt_order_str.split(", "))
                if self.blt_order == ("bda",):
                    self._blt_order.form = (1,)
            self.history = str(vis_hdr.get("HISTORY", ""))
            if not uvutils._check_history_version(

self.data_array = data

                # Note: This axis is called STOKES by analogy with the equivalent
                # uvfits axis
                # However, this is confusing because it is NOT a true Stokes axis,
                #   it is really the polarization axis.
                if (
                    primary_header.pop("CTYPE" + str(ax_nums["feed_pol"]))
                    .lower()
                    .strip()
                    == "stokes"
                ):
                    self.Npols = primary_header.pop("NAXIS" + str(ax_nums["feed_pol"]))

                self.polarization_array = np.int32(
                    uvutils._fits_gethduaxis(primary_hdu, ax_nums["feed_pol"])
                )
                self._set_power()
            elif self.beam_type == "efield":
                self._set_efield()
                if n_dimensions < n_efield_dims:
                    raise ValueError(
                        "beam_type is efield and data dimensionality is too low"
                    )
                complex_arrs = np.split(data, 2, axis=0)
                self.data_array = np.squeeze(
                    complex_arrs[0] + 1j * complex_arrs[1], axis=0
                )
                if (
                    primary_header.pop("CTYPE" + str(ax_nums["feed_pol"]))
                    .lower()
                    .strip()

self.delay_array = data[:, :, :, :, :, 0]
                self.quality_array = data[:, :, :, :, :, 1]

                sechdu = fname[hdunames["FLAGS"]]
                flag_data = sechdu.data
                if sechdu.header["NAXIS1"] == 2:
                    self.flag_array = flag_data[:, :, :, :, :, 0].astype("bool")
                    self.input_flag_array = flag_data[:, :, :, :, :, 1].astype("bool")
                else:
                    self.flag_array = flag_data[:, :, :, :, :, 0].astype("bool")

                # generate frequency array from flag data unit
                # (no freq axis in primary).
                self.Nfreqs = sechdu.header["NAXIS4"]
                self.freq_array = uvutils._fits_gethduaxis(sechdu, 4)
                self.freq_array.shape = (self.Nspws,) + self.freq_array.shape

                spw_array = uvutils._fits_gethduaxis(sechdu, 5) - 1

                if not np.allclose(spw_array, self.spw_array):
                    raise ValueError(
                        "Spectral window values are different in FLAGS HDU than"
                        " in primary HDU"
                    )

                time_array = uvutils._fits_gethduaxis(sechdu, 3)
                if not np.allclose(
                    time_array,
                    self.time_array,
                    rtol=self._time_array.tols[0],
                    atol=self._time_array.tols[0],

"primary HDU"
                        )
                else:
                    basisvec_ax_nums = reg_basisvec_ax_nums
                    basisvec_coord_list = [
                        basisvec_header[
                            "CTYPE" + str(basisvec_ax_nums["img_ax1"])
                        ].lower(),
                        basisvec_header[
                            "CTYPE" + str(basisvec_ax_nums["img_ax2"])
                        ].lower(),
                    ]
                    basisvec_axis1_array = uvutils._fits_gethduaxis(
                        basisvec_hdu, basisvec_ax_nums["img_ax1"]
                    )
                    basisvec_axis2_array = uvutils._fits_gethduaxis(
                        basisvec_hdu, basisvec_ax_nums["img_ax2"]
                    )

                    # if units aren't defined they are degrees by FITS convention
                    # convert degrees to radians because UVBeam uses radians
                    axis1_units = basisvec_header.pop(
                        "CUNIT" + str(basisvec_ax_nums["img_ax1"]), "deg"
                    )
                    if axis1_units == "deg":
                        basisvec_axis1_array = np.deg2rad(basisvec_axis1_array)
                    elif axis1_units != "rad":
                        raise ValueError(
                            "Units of first axis array in BASISVEC HDU are not"
                            " 'deg' or 'rad'."
                        )
                    axis2_units = basisvec_header.pop(

"Spectral window values are different in FLAGS HDU than"
                        " in primary HDU"
                    )

                time_array = uvutils._fits_gethduaxis(sechdu, 3)
                if not np.allclose(
                    time_array,
                    self.time_array,
                    rtol=self._time_array.tols[0],
                    atol=self._time_array.tols[0],
                ):
                    raise ValueError(
                        "Time values are different in FLAGS HDU than in primary HDU"
                    )

                jones_array = uvutils._fits_gethduaxis(sechdu, 2)
                if not np.allclose(
                    jones_array,
                    self.jones_array,
                    rtol=self._jones_array.tols[0],
                    atol=self._jones_array.tols[0],
                ):
                    raise ValueError(
                        "Jones values are different in FLAGS HDU than in primary HDU"
                    )

            self.extra_keywords = uvutils._get_fits_extra_keywords(hdr)

            # get total quality array if present
            if "TOTQLTY" in hdunames:
                totqualhdu = fname[hdunames["TOTQLTY"]]
                self.total_quality_array = totqualhdu.data

hdunames = uvutils._fits_indexhdus(hdu_list)  # find the rest of the tables

            # First get everything we can out of the header.
            self._set_phased()
            # check if we have an spw dimension
            if vis_hdr["NAXIS"] == 7:
                if vis_hdr["NAXIS5"] > 1:
                    raise ValueError(
                        "Sorry.  Files with more than one spectral"
                        "window (spw) are not yet supported. A "
                        "great project for the interested student!"
                    )

                self.Nspws = vis_hdr.pop("NAXIS5")

                self.spw_array = np.int32(uvutils._fits_gethduaxis(vis_hdu, 5)) - 1

                # the axis number for phase center depends on if the spw exists
                self.phase_center_ra_degrees = np.float(vis_hdr.pop("CRVAL6"))
                self.phase_center_dec_degrees = np.float(vis_hdr.pop("CRVAL7"))
            else:
                self.Nspws = 1
                self.spw_array = np.array([0])

                # the axis number for phase center depends on if the spw exists
                self.phase_center_ra_degrees = np.float(vis_hdr.pop("CRVAL5"))
                self.phase_center_dec_degrees = np.float(vis_hdr.pop("CRVAL6"))

            # get shapes
            self.Nfreqs = vis_hdr.pop("NAXIS4")
            self.Npols = vis_hdr.pop("NAXIS3")
            self.Nblts = vis_hdr.pop("GCOUNT")

ax_nums = hpx_primary_ax_nums
                self.nside = primary_header.pop("NSIDE", None)
                self.ordering = primary_header.pop("ORDERING", None)
                data_Npixels = primary_header.pop("NAXIS" + str(ax_nums["pixel"]))
                if data_Npixels != self.Npixels:
                    raise ValueError(
                        "Number of pixels in HPX_IND extension does "
                        "not match number of pixels in data array"
                    )
            else:
                ax_nums = reg_primary_ax_nums
                self.Naxes1 = primary_header.pop("NAXIS" + str(ax_nums["img_ax1"]))
                self.Naxes2 = primary_header.pop("NAXIS" + str(ax_nums["img_ax2"]))

                self.axis1_array = uvutils._fits_gethduaxis(
                    primary_hdu, ax_nums["img_ax1"]
                )
                self.axis2_array = uvutils._fits_gethduaxis(
                    primary_hdu, ax_nums["img_ax2"]
                )

                # if units aren't defined they are degrees by FITS convention
                # convert degrees to radians because UVBeam uses radians
                axis1_units = primary_header.pop(
                    "CUNIT" + str(ax_nums["img_ax1"]), "deg"
                )
                if axis1_units == "deg":
                    self.axis1_array = np.deg2rad(self.axis1_array)
                elif axis1_units != "rad":
                    raise ValueError(
                        'Units of first axis array are not "deg" or "rad".'

)

            if (
                self.Nspws is None or self.Naxes_vec is None
            ) and self.beam_type == "power":
                if self.Nspws is None:
                    self.Nspws = 1
                    self.spw_array = np.array([0])
                if self.Naxes_vec is None:
                    self.Naxes_vec = 1

                # add extra empty dimensions to data_array as appropriate
                while len(self.data_array.shape) < n_efield_dims - 1:
                    self.data_array = np.expand_dims(self.data_array, axis=0)

            self.freq_array = uvutils._fits_gethduaxis(primary_hdu, ax_nums["freq"])
            self.freq_array.shape = (self.Nspws,) + self.freq_array.shape
            # default frequency axis is Hz, but check for corresonding CUNIT
            freq_units = primary_header.pop("CUNIT" + str(ax_nums["freq"]), "Hz")
            if freq_units != "Hz":
                freq_factor = {"kHz": 1e3, "MHz": 1e6, "GHz": 1e9}
                if freq_units in freq_factor.keys():
                    self.freq_array = self.freq_array * freq_factor[freq_units]
                else:
                    raise ValueError("Frequency units not recognized.")

            self.history = str(primary_header.get("HISTORY", ""))
            if not uvutils._check_history_version(
                self.history, self.pyuvdata_version_str
            ):
                self.history += self.pyuvdata_version_str

"Frequency values are different in TOTQLTY HDU than"
                            " in primary HDU"
                        )

                time_array = uvutils._fits_gethduaxis(totqualhdu, 2)
                if not np.allclose(
                    time_array,
                    self.time_array,
                    rtol=self._time_array.tols[0],
                    atol=self._time_array.tols[0],
                ):
                    raise ValueError(
                        "Time values are different in TOTQLTY HDU than in primary HDU"
                    )

                jones_array = uvutils._fits_gethduaxis(totqualhdu, 1)
                if not np.allclose(
                    jones_array,
                    self.jones_array,
                    rtol=self._jones_array.tols[0],
                    atol=self._jones_array.tols[0],
                ):
                    raise ValueError(
                        "Jones values are different in TOTQLTY HDU than in primary HDU"
                    )

            else:
                self.total_quality_array = None

        if run_check:
            self.check(
                check_extra=check_extra, run_check_acceptability=run_check_acceptability