import numpy as np
import struct
import warnings
import string
from astropy import log
from astropy.io import registry as io_registry
from ..spectral_cube import BaseSpectralCube
from .fits import load_fits_cube

"""
.. TODO::
    When any section length is zero, that means the following values are to be
    ignored.  No warning is needed.
"""

# Constant:
r2deg = 180/np.pi

# see sicfits.f90
_ctype_dict={'LII':'GLON',
             'BII':'GLAT',
             'VELOCITY':'VELO',
             'RA':'RA',
             'DEC':'DEC',
             'FREQUENCY': 'FREQ',
            }
_cunit_dict = {'LII':'deg',
               'BII':'deg',
               'VELOCITY':'km s-1',
               'RA':'deg',
               'DEC':'deg',
               'FREQUENCY': 'MHz',
              }
cel_types = ('RA','DEC','GLON','GLAT')

# CLASS apparently defaults to an ARC (zenithal equidistant) projection; this
# is what is output in case the projection # is zero when exporting from CLASS
_proj_dict = {0:'ARC', 1:'TAN', 2:'SIN', 3:'AZP', 4:'STG', 5:'ZEA', 6:'AIT',
              7:'GLS', 8:'SFL', }
_bunit_dict = {'k (tmb)': 'K'}

def is_lmv(origin, filepath, fileobj, *args, **kwargs):
    """
    Determine whether input is in GILDAS CLASS lmv format
    """
    return filepath is not None and filepath.lower().endswith('.lmv')

def read_lmv(lf):
    """
    Read an LMV cube file

    Specification is primarily in GILDAS image_def.f90
    """
    log.warning("CLASS LMV cube reading is tentatively supported.  "
             "Please post bug reports at the first sign of danger!")

    # lf for "LMV File"
    filetype = _read_string(lf, 12)
    #!---------------------------------------------------------------------
    #! @ private
    #!       SYCODE system code
    #!       '-'    IEEE
    #!       '.'    EEEI (IBM like)
    #!       '_'    VAX
    #!       IMCODE file code
    #!       '<'    IEEE  64 bits    (Little Endian, 99.9 % of recent computers)
    #!       '>'    EEEI  64 bits    (Big Endian, HPUX, IBM-RISC, and SPARC ...)
    #!---------------------------------------------------------------------
    imcode = filetype[6]
    if filetype[:6] != 'GILDAS' or filetype[7:] != 'IMAGE':
        raise TypeError("File is not a GILDAS Image file")

    if imcode in ('<','>'):
        if imcode =='>':
            log.warning("Swap the endianness first...")
        return read_lmv_type2(lf)
    else:
        return read_lmv_type1(lf)

def read_lmv_type1(lf):
    header = {}
    # fmt probably matters!  Default is "r4", i.e. float32 data, but could be float64
    fmt = np.fromfile(lf, dtype='int32', count=1) # 4
    # number of data blocks
    ndb = np.fromfile(lf, dtype='int32', count=1) # 5
    gdf_type = np.fromfile(lf, dtype='int32', count=1) # 6
    # Reserved Space
    reserved_fill = np.fromfile(lf, dtype='int32', count=4) # 7
    general_section_length = np.fromfile(lf, dtype='int32', count=1) # 11
    #print "Format: ",fmt," ndb: ",ndb, " fill: ",fill," other: ",unknown

    # pos 12
    naxis,naxis1,naxis2,naxis3,naxis4 = np.fromfile(lf,count=5,dtype='int32')
    header['NAXIS'] = naxis
    header['NAXIS1'] = naxis1
    header['NAXIS2'] = naxis2
    header['NAXIS3'] = naxis3
    header['NAXIS4'] = naxis4

    # We are indexing bytes from here; CLASS indices are higher by 12
    # pos 17
    header['CRPIX1'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CRVAL1'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CDELT1'] = np.fromfile(lf,count=1,dtype='float64')[0] * r2deg
    header['CRPIX2'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CRVAL2'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CDELT2'] = np.fromfile(lf,count=1,dtype='float64')[0] * r2deg
    header['CRPIX3'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CRVAL3'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CDELT3'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CRPIX4'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CRVAL4'] = np.fromfile(lf,count=1,dtype='float64')[0]
    header['CDELT4'] = np.fromfile(lf,count=1,dtype='float64')[0]
    # pos 41
    #print "Post-crval",lf.tell()
    blank_section_length = np.fromfile(lf,count=1,dtype='int32')
    if blank_section_length != 8:
        warnings.warn("Invalid section length found for blanking section")
    bval = np.fromfile(lf,count=1,dtype='float32')[0] # 42
    header['TOLERANC'] = np.fromfile(lf,count=1,dtype='int32')[0] # 43 eval = tolerance
    extrema_section_length = np.fromfile(lf,count=1,dtype='int32')[0] # 44
    if extrema_section_length != 40:
        warnings.warn("Invalid section length found for extrema section")
    vmin,vmax = np.fromfile(lf,count=2,dtype='float32') # 45
    xmin,xmax,ymin,ymax,zmin,zmax = np.fromfile(lf,count=6,dtype='int32') # 47
    wmin,wmax = np.fromfile(lf,count=2,dtype='int32') # 53
    description_section_length = np.fromfile(lf,count=1,dtype='int32')[0] # 55
    if description_section_length != 72:
        warnings.warn("Invalid section length found for description section")
    #strings = lf.read(description_section_length) # 56
    header['BUNIT']  = _read_string(lf, 12) # 56
    header['CTYPE1'] = _read_string(lf, 12) # 59
    header['CTYPE2'] = _read_string(lf, 12) # 62
    header['CTYPE3'] = _read_string(lf, 12) # 65
    header['CTYPE4'] = _read_string(lf, 12) # 68
    header['CUNIT1'] = _cunit_dict[header['CTYPE1'].strip()]
    header['CUNIT2'] = _cunit_dict[header['CTYPE2'].strip()]
    header['CUNIT3'] = _cunit_dict[header['CTYPE3'].strip()]
    header['COOSYS'] = _read_string(lf, 12) # 71
    position_section_length = np.fromfile(lf,count=1,dtype='int32') # 74
    if position_section_length != 48:
        warnings.warn("Invalid section length found for position section")
    header['OBJNAME'] = _read_string(lf, 4*3) # 75
    header['RA'] = np.fromfile(lf, count=1, dtype='float64')[0] * r2deg # 78
    header['DEC'] = np.fromfile(lf, count=1, dtype='float64')[0] * r2deg # 80
    header['GLON'] = np.fromfile(lf, count=1, dtype='float64')[0] * r2deg # 82
    header['GLAT'] = np.fromfile(lf, count=1, dtype='float64')[0] * r2deg # 84
    header['EQUINOX'] = np.fromfile(lf,count=1,dtype='float32')[0] # 86
    header['PROJWORD'] = _read_string(lf, 4) # 87
    header['PTYP'] = np.fromfile(lf,count=1,dtype='int32')[0] # 88
    header['A0'] = np.fromfile(lf,count=1,dtype='float64')[0] # 89
    header['D0'] = np.fromfile(lf,count=1,dtype='float64')[0] # 91
    header['PANG'] = np.fromfile(lf,count=1,dtype='float64')[0] # 93
    header['XAXI'] = np.fromfile(lf,count=1,dtype='float32')[0] # 95
    header['YAXI'] = np.fromfile(lf,count=1,dtype='float32')[0] # 96
    spectroscopy_section_length = np.fromfile(lf,count=1,dtype='int32') # 97
    if spectroscopy_section_length != 48:
        warnings.warn("Invalid section length found for spectroscopy section")
    header['RECVR'] = _read_string(lf, 12) # 98
    header['FRES'] = np.fromfile(lf,count=1,dtype='float64')[0] # 101
    header['IMAGFREQ'] = np.fromfile(lf,count=1,dtype='float64')[0] # 103 "FIMA"
    header['REFFREQ'] = np.fromfile(lf,count=1,dtype='float64')[0] # 105
    header['VRES'] = np.fromfile(lf,count=1,dtype='float32')[0] # 107
    header['VOFF'] = np.fromfile(lf,count=1,dtype='float32')[0] # 108
    header['FAXI'] = np.fromfile(lf,count=1,dtype='int32')[0] # 109
    resolution_section_length = np.fromfile(lf,count=1,dtype='int32')[0] # 110
    if resolution_section_length != 12:
        warnings.warn("Invalid section length found for resolution section")
    #header['DOPP'] = np.fromfile(lf,count=1,dtype='float16')[0] # 110a ???
    #header['VTYP'] = np.fromfile(lf,count=1,dtype='int16')[0] # 110b
    # integer, parameter :: vel_unk = 0      ! Unsupported referential :: planetary...)
    # integer, parameter :: vel_lsr = 1      ! LSR referential
    # integer, parameter :: vel_hel = 2      ! Heliocentric referential
    # integer, parameter :: vel_obs = 3      ! Observatory referential
    # integer, parameter :: vel_ear = 4      ! Earth-Moon barycenter referential
    # integer, parameter :: vel_aut = -1     ! Take referential from data
    header['BMAJ'] = np.fromfile(lf,count=1,dtype='float32')[0] # 111
    header['BMIN'] = np.fromfile(lf,count=1,dtype='float32')[0] # 112
    header['BPA'] = np.fromfile(lf,count=1,dtype='float32')[0] # 113
    noise_section_length = np.fromfile(lf,count=1,dtype='int32')
    if noise_section_length != 0:
        warnings.warn("Invalid section length found for noise section")
    header['NOISE'] = np.fromfile(lf,count=1,dtype='float32')[0] # 115
    header['RMS'] = np.fromfile(lf,count=1,dtype='float32')[0] # 116
    astrometry_section_length = np.fromfile(lf,count=1,dtype='int32')
    if astrometry_section_length != 0:
        warnings.warn("Invalid section length found for astrometry section")
    header['MURA'] = np.fromfile(lf,count=1,dtype='float32')[0] # 118
    header['MUDEC'] = np.fromfile(lf,count=1,dtype='float32')[0] # 119
    header['PARALLAX'] = np.fromfile(lf,count=1,dtype='float32')[0] # 120

    # Apparently CLASS headers aren't required to fill the 'value at
    # reference pixel' column
    if (header['CTYPE1'].strip() == 'RA' and header['CRVAL1'] == 0 and
        header['RA'] != 0):
        header['CRVAL1'] = header['RA']
        header['CRVAL2'] = header['DEC']


    # Copied from the type 2 reader:
    # Use the appropriate projection type
    ptyp = header['PTYP']
    for kw in header:
        if 'CTYPE' in kw:
            if header[kw].strip() in cel_types:
                n_dashes = 5-len(header[kw].strip())
                header[kw] = header[kw].strip()+ '-'*n_dashes + _proj_dict[ptyp]

    other_info = np.fromfile(lf, count=7, dtype='float32') # 121-end
    if not np.all(other_info == 0):
        warnings.warn("Found additional information in the last 7 bytes")

    endpoint = 508
    if lf.tell() != endpoint:
        raise ValueError("Header was not parsed correctly")

    data = np.fromfile(lf, count=naxis1*naxis2*naxis3, dtype='float32')

    data[data == bval] = np.nan

    # for no apparent reason, y and z are 1-indexed and x is zero-indexed
    if (wmin-1,zmin-1,ymin-1,xmin) != np.unravel_index(np.nanargmin(data),
                                                       [naxis4,naxis3,naxis2,naxis1]):
        warnings.warn("Data min location does not match that on file.  "
                      "Possible error reading data.")
    if (wmax-1,zmax-1,ymax-1,xmax) != np.unravel_index(np.nanargmax(data),
                                                       [naxis4,naxis3,naxis2,naxis1]):
        warnings.warn("Data max location does not match that on file.  "
                      "Possible error reading data.")
    if np.nanmax(data) != vmax:
        warnings.warn("Data max does not match that on file.  "
                      "Possible error reading data.")
    if np.nanmin(data) != vmin:
        warnings.warn("Data min does not match that on file.  "
                      "Possible error reading data.")

    return data.reshape([naxis4,naxis3,naxis2,naxis1]),header
    # debug
    #return data.reshape([naxis3,naxis2,naxis1]), header, hdr_f, hdr_s, hdr_i, hdr_d, hdr_d_2

def read_lmv_tofits(fileobj):
    from astropy.io import fits
    data,header = read_lmv(fileobj)
    # LMV may contain extra dimensions that are improperly labeled
    data = data.squeeze()
    bad_kws = ['NAXIS4','CRVAL4','CRPIX4','CDELT4','CROTA4','CUNIT4','CTYPE4']

    cards = [fits.header.Card(keyword=k, value=v[0], comment=v[1])
             if isinstance(v, tuple) else
             fits.header.Card(''.join(s for s in k if s in string.printable),
                              ''.join(s for s in v if s in string.printable)
                              if isinstance(v, str) else v)
             for k,v in header.items()
             if k not in bad_kws]
    Header = fits.Header(cards)
    hdu = fits.PrimaryHDU(data=data, header=Header)
    return hdu

def load_lmv_cube(fileobj, target_cls=None, use_dask=None):
    hdu = read_lmv_tofits(fileobj)
    meta = {'filename':fileobj.name}

    return load_fits_cube(hdu, meta=meta, use_dask=use_dask)

def _read_byte(f):
    '''Read a single byte (from idlsave)'''
    return np.uint8(struct.unpack('=B', f.read(4)[:1])[0])

def _read_int16(f):
    '''Read a signed 16-bit integer (from idlsave)'''
    return np.int16(struct.unpack('=h', f.read(4)[2:4])[0])

def _read_int32(f):
    '''Read a signed 32-bit integer (from idlsave)'''
    return np.int32(struct.unpack('=i', f.read(4))[0])

def _read_int64(f):
    '''Read a signed 64-bit integer '''
    return np.int64(struct.unpack('=q', f.read(8))[0])

def _read_float32(f):
    '''Read a 32-bit float (from idlsave)'''
    return np.float32(struct.unpack('=f', f.read(4))[0])

def _read_string(f, size):
    '''Read a string of known maximum length'''
    return f.read(size).decode('utf-8').strip()

def _read_float64(f):
    '''Read a 64-bit float (from idlsave)'''
    return np.float64(struct.unpack('=d', f.read(8))[0])

def _check_val(name, got,expected):
    if got != expected:
        log.warning("{2} = {0} instead of {1}".format(got, expected, name))

def read_lmv_type2(lf):
    """ See image_def.f90 """

    header = {}
    lf.seek(12)
    # DONE before integer(kind=4) :: ijtyp(3) = 0       !  1  Image Type

    # fmt probably matters!  Default is "r4", i.e. float32 data, but could be float64
    fmt = _read_int32(lf) # 4
    # number of data blocks
    ndb = _read_int64(lf) # 5
    nhb = _read_int32(lf) # 7
    ntb = _read_int32(lf) # 8
    version_gdf = _read_int32(lf) # 9
    if version_gdf != 20:
        raise TypeError("Trying to read a version-2 file, but the version"
                        " number is {0} (should be 20)".format(version_gdf))
    type_gdf = _read_int32(lf) # 10
    dim_start = _read_int32(lf) # 11
    pad_trail = _read_int32(lf) # 12

    if dim_start % 2 == 0:
        log.warning("Got even dim_start in lmv cube: this is not expected.")
    if dim_start > 17:
        log.warning("dim_start > 17 in lmv cube: this is not expected.")

    lf.seek(16*4)
    gdf_maxdims=7
    dim_words = _read_int32(lf) # 17
    if dim_words != 2*gdf_maxdims+2:
        log.warning("dim_words = {0} instead of {1}".format(dim_words,
                                                         gdf_maxdims*2+2))
    blan_start = _read_int32(lf) # 18
    if blan_start != dim_start+dim_words+2:
        log.warning("blan_star = {0} instead of {1}".format(blan_start,
                                                         dim_start+dim_words+2))

    mdim = _read_int32(lf) # 19
    ndim = _read_int32(lf) # 20
    dims = np.fromfile(lf, count=gdf_maxdims, dtype='int64')
    if np.count_nonzero(dims) != ndim:
        raise ValueError("Disagreement between ndims and number of nonzero dims.")

    header['NAXIS'] = ndim
    valid_dims = []
    for ii,dim in enumerate(dims):
        if dim != 0:
            header['NAXIS{0}'.format(ii+1)] = dim
            valid_dims.append(ii)


    blan_words = _read_int32(lf)
    if blan_words != 2:
        log.warning("blan_words = {0} instead of 2".format(blan_words))
    extr_start = _read_int32(lf)
    bval = _read_float32(lf) # blanking value
    bval_tol = _read_float32(lf) # eval = tolerance

    # FITS requires integer BLANKs
    #header['BLANK'] = bval

    extr_words = _read_int32(lf)
    if extr_words != 6:
        log.warning("extr_words = {0} instead of 6".format(extr_words))
    coor_start = _read_int32(lf)
    if coor_start != extr_start+extr_words+2:
        log.warning("coor_start = {0} instead of {1}".format(coor_start,
                                                          extr_start+extr_words+2))
    rmin = _read_float32(lf)
    rmax = _read_float32(lf)

    # position 168
    minloc = _read_int64(lf)
    maxloc = _read_int64(lf)

    # lf.seek(184)
    coor_words = _read_int32(lf)
    if coor_words != gdf_maxdims*6:
        log.warning("coor_words = {0} instead of {1}".format(coor_words,
                                                          gdf_maxdims*6))
    desc_start = _read_int32(lf)
    if desc_start != coor_start+coor_words+2:
        log.warning("desc_start = {0} instead of {1}".format(desc_start,
                                                          coor_start+coor_words+2))

    convert = np.fromfile(lf, count=3*gdf_maxdims, dtype='float64').reshape([gdf_maxdims,3])

    # conversion of "convert" to CRPIX/CRVAL/CDELT below

    desc_words = _read_int32(lf)
    if desc_words != 3*(gdf_maxdims+1):
        log.warning("desc_words = {0} instead of {1}".format(desc_words,
                                                          3*(gdf_maxdims+1)))
    null_start = _read_int32(lf)
    if null_start != desc_start+desc_words+2:
        log.warning("null_start = {0} instead of {1}".format(null_start,
                                                          desc_start+desc_words+2))
    ijuni = _read_string(lf, 12) # data unit
    ijcode = [_read_string(lf, 12) for ii in range(gdf_maxdims)]
    pad_desc = _read_int32(lf)

    if ijuni.lower() in _bunit_dict:
        header['BUNIT'] = (_bunit_dict[ijuni.lower()],
                           ijuni)
    else:
        header['BUNIT'] = ijuni

     #! The first block length is thus
     #!	s_dim-1 + (2*mdim+4) + (4) + (8) +  (6*mdim+2) + (3*mdim+5)
     #! = s_dim-1 + mdim*(2+6+3) + (4+4+2+5+8)
     #! = s_dim-1 + 11*mdim + 23
     #! With mdim = 7, s_dim=11, this is 110 spaces
     #! With mdim = 8, s_dim=11, this is 121 spaces
     #! MDIM > 8 would NOT fit in one block...
     #!
     #! Block 2: Ancillary information
     #!
     #! The same logic of Length + Pointer is used there too, although the
     #! length are fixed. Note rounding to even number for the pointer offsets
     #! in order to preserve alignement...
     #!

    lf.seek(512)
    posi_words = _read_int32(lf)
    _check_val('posi_words', posi_words, 15)

    proj_start = _read_int32(lf)
    source_name = _read_string(lf, 12)
    header['OBJECT'] = source_name
    coordinate_system = _read_string(lf, 12)

    header['RA'] = _read_float64(lf)
    header['DEC'] = _read_float64(lf)
    header['LII'] = _read_float64(lf)
    header['BII'] = _read_float64(lf)
    header['EPOCH'] = _read_float32(lf)
    #pad_posi = _read_float32(lf)
    #print pad_posi
    #raise ValueError("pad_posi should probably be 0?")

    #! PROJECTION
    #integer(kind=4) :: proj_words = 9     ! Projection length: 9 used + 1 padding
    #integer(kind=4) :: spec_start !! = proj_start + 12
    #real(kind=8) :: a0      = 0.d0        ! 89  X of projection center
    #real(kind=8) :: d0      = 0.d0        ! 91  Y of projection center
    #real(kind=8) :: pang    = 0.d0        ! 93  Projection angle
    #integer(kind=4) :: ptyp = p_none      ! 88  Projection type (see p_... codes)
    #integer(kind=4) :: xaxi = 0           ! 95  X axis
    #integer(kind=4) :: yaxi = 0           ! 96  Y axis
    #integer(kind=4) :: pad_proj
    #!

    proj_words = _read_int32(lf)
    spec_start = _read_int32(lf)
    _check_val('spec_start', spec_start, proj_start+proj_words+2)
    if proj_words == 9:
        header['PROJ_A0'] = _read_float64(lf)
        header['PROJ_D0'] = _read_float64(lf)
        header['PROJPANG'] = _read_float64(lf)
        ptyp = _read_int32(lf)
        header['PROJXAXI'] = _read_int32(lf)
        header['PROJYAXI'] = _read_int32(lf)
    elif proj_words != 0:
        raise ValueError("Invalid # of projection keywords")

    for kw in header:
        if 'CTYPE' in kw:
            if header[kw].strip() in cel_types:
                n_dashes = 5-len(header[kw].strip())
                header[kw] = header[kw].strip()+ '-'*n_dashes + _proj_dict[ptyp]

    for ii,((ref,val,inc),code) in enumerate(zip(convert,ijcode)):
        if ii in valid_dims:
            # jul14a gio/to_imfits.f90 line 284-313
            if ptyp != 0 and (ii+1) in (header['PROJXAXI'],
                                        header['PROJYAXI']):
                #! Compute reference pixel so that VAL(REF) = 0
                ref = ref - val/inc
                if (ii+1) == header['PROJXAXI']:
                    val = header['PROJ_A0']
                elif (ii+1) == header['PROJYAXI']:
                    val = header['PROJ_D0']
                else:
                    raise ValueError("Impossible state - code bug.")
                val = val*r2deg
                inc = inc*r2deg
                rota = r2deg*header['PROJPANG']
            elif code in ('RA', 'L', 'B', 'DEC', 'LII', 'BII', 'GLAT',
                          'GLON', 'LAT', 'LON'):
                val = val*r2deg
                inc = inc*r2deg
                rota = 0.0
            # These are not implemented: prefer to maintain original units (we're
            # reading in to spectral_cube after all, no need to change units until the
            # output step)
            #elseif (code.eq.'FREQUENCY') then
            #val = val*1.0d6            ! MHz to Hz
            #inc = inc*1.0d6
            #elseif (code.eq.'VELOCITY') then
            #code = 'VRAD'              ! force VRAD instead of VELOCITY for CASA
            #val = val*1.0d3            ! km/s to m/s
            #inc = inc*1.0d3

            header['CRPIX{0}'.format(ii+1)] = ref
            header['CRVAL{0}'.format(ii+1)] = val
            header['CDELT{0}'.format(ii+1)] = inc


    for ii,ctype in enumerate(ijcode):
        if ii in valid_dims:
            header['CTYPE{0}'.format(ii+1)] = _ctype_dict[ctype]
            header['CUNIT{0}'.format(ii+1)] = _cunit_dict[ctype]

    spec_words = _read_int32(lf)
    reso_start = _read_int32(lf)
    _check_val('reso_start', reso_start, proj_start+proj_words+2+spec_words+2)
    if spec_words == 14:
        header['FRES'] = _read_float64(lf)
        header['FIMA'] = _read_float64(lf)
        header['FREQ'] = _read_float64(lf)
        header['VRES'] = _read_float32(lf)
        header['VOFF'] = _read_float32(lf)
        header['DOPP'] = _read_float32(lf)
        header['FAXI'] = _read_int32(lf)
        header['LINENAME'] = _read_string(lf, 12)
        header['VTYPE'] = _read_int32(lf)
    elif spec_words != 0:
        raise ValueError("Invalid # of spectroscopic keywords")

     #! SPECTROSCOPY
     #integer(kind=4) :: spec_words  = 14   ! Spectroscopy length: 14 used
     #integer(kind=4) :: reso_start  !! = spec_words + 16
     #real(kind=8) :: fres       = 0.d0     !101  Frequency resolution
     #real(kind=8) :: fima       = 0.d0     !103  Image frequency
     #real(kind=8) :: freq       = 0.d0     !105  Rest Frequency
     #real(kind=4) :: vres       = 0.0      !107  Velocity resolution
     #real(kind=4) :: voff       = 0.0      !108  Velocity offset
     #real(kind=4) :: dopp       = 0.0      !     Doppler factor
     #integer(kind=4) :: faxi    = 0        !109  Frequency axis
     #integer(kind=4) :: ijlin(3) = 0       ! 98  Line name
     #integer(kind=4) :: vtyp    = vel_unk  ! Velocity type (see vel_... codes)

    reso_words = _read_int32(lf)
    nois_start = _read_int32(lf)
    _check_val('nois_start', nois_start, proj_start+proj_words+2+spec_words+2+reso_words+2)
    if reso_words == 3:
        header['BMAJ'] = _read_float32(lf)
        header['BMIN'] = _read_float32(lf)
        header['BPA'] = _read_float32(lf)
        #pad_reso = _read_float32(lf)
    elif reso_words != 0:
        raise ValueError("Invalid # of resolution keywords")

     #! RESOLUTION
     #integer(kind=4) :: reso_words = 3     ! Resolution length: 3 used + 1 padding
     #integer(kind=4) :: nois_start !! = reso_words + 6
     #real(kind=4) :: majo    = 0.0         !111  Major axis
     #real(kind=4) :: mino    = 0.0         !112  Minor axis
     #real(kind=4) :: posa    = 0.0         !113  Position angle
     #real(kind=4) :: pad_reso


    nois_words = _read_int32(lf)
    astr_start = _read_int32(lf)
    _check_val('astr_start', astr_start, proj_start+proj_words+2+spec_words+2+reso_words+2+nois_words+2)
    if nois_words == 2:
        header['NOISE_T'] = (_read_float32(lf), "Theoretical Noise")
        header['NOISERMS'] = (_read_float32(lf), "Measured (RMS) noise")
    elif nois_words != 0:
        raise ValueError("Invalid # of noise keywords")

     #! NOISE
     #integer(kind=4) :: nois_words = 2     ! Noise section length: 2 used
     #integer(kind=4) :: astr_start !! = s_nois + 4
     #real(kind=4) :: noise   = 0.0         ! 115 Theoretical noise
     #real(kind=4) :: rms     = 0.0         ! 116 Actual noise

    astr_words = _read_int32(lf)
    uvda_start = _read_int32(lf)
    _check_val('uvda_start', uvda_start, proj_start+proj_words+2+spec_words+2+reso_words+2+nois_words+2+astr_words+2)
    if astr_words == 3:
        header['MURA'] = _read_float32(lf)
        header['MUDEC'] = _read_float32(lf)
        header['PARALLAX'] = _read_float32(lf)
    elif astr_words != 0:
        raise ValueError("Invalid # of astrometry keywords")

     #! ASTROMETRY
     #integer(kind=4) :: astr_words = 3    ! Proper motion section length: 3 used + 1 padding
     #integer(kind=4) :: uvda_start !! = s_astr + 4
     #real(kind=4) :: mura     = 0.0        ! 118 along RA, in mas/yr
     #real(kind=4) :: mudec    = 0.0        ! 119 along Dec, in mas/yr
     #real(kind=4) :: parallax = 0.0        ! 120 in mas
     #real(kind=4) :: pad_astr
     #! real(kind=4) :: pepoch   = 2000.0     ! 121 in yrs ?

    code_uvt_last=25
    uvda_words = _read_int32(lf)
    void_start = _read_int32(lf)
    _check_val('void_start', void_start, proj_start + proj_words + 2 +
               spec_words + 2 + reso_words + 2 + nois_words + 2 + astr_words +
               2 + uvda_words + 2)
    if uvda_words == 18+2*code_uvt_last:
        version_uv = _read_int32(lf)
        nchan = _read_int32(lf)
        nvisi = _read_int64(lf)
        nstokes = _read_int32(lf)
        natom = _read_int32(lf)
        basemin = _read_float32(lf)
        basemax = _read_float32(lf)
        fcol = _read_int32(lf)
        lcol = _read_int32(lf)
        nlead = _read_int32(lf)
        ntrail = _read_int32(lf)
        column_pointer = np.fromfile(lf, count=code_uvt_last, dtype='int32')
        column_size = np.fromfile(lf, count=code_uvt_last, dtype='int32')
        column_codes = np.fromfile(lf, count=nlead+ntrail, dtype='int32')
        column_types = np.fromfile(lf, count=nlead+ntrail, dtype='int32')
        order = _read_int32(lf)
        nfreq = _read_int32(lf)
        atoms = np.fromfile(lf, count=4, dtype='int32')

    elif uvda_words != 0:
        raise ValueError("Invalid # of UV data keywords")

     #! UV_DATA information
     #integer(kind=4) :: uvda_words  = 18+2*code_uvt_last ! Length of section: 14 used
     #integer(kind=4) :: void_start  !! = s_uvda + l_uvda + 2
     #integer(kind=4) :: version_uv = code_version_uvt_current  ! 1 version number. Will allow us to change the data format
     #integer(kind=4) :: nchan = 0         ! 2 Number of channels
     #integer(kind=8) :: nvisi = 0         ! 3-4 Independent of the transposition status
     #integer(kind=4) :: nstokes = 0       ! 5 Number of polarizations
     #integer(kind=4) :: natom = 0         ! 6. 3 for real, imaginary, weight. 1 for real.
     #real(kind=4)    :: basemin = 0.      ! 7 Minimum Baseline
     #real(kind=4)    :: basemax = 0.      ! 8 Maximum Baseline
     #integer(kind=4) :: fcol              ! 9 Column of first channel
     #integer(kind=4) :: lcol              ! 10 Column of last  channel
     #! The number of information per channel can be obtained by
     #!       (lcol-fcol+1)/(nchan*natom)
     #! so this could allow to derive the number of Stokes parameters
     #! Leading data at start of each visibility contains specific information
     #integer(kind=4) :: nlead = 7           ! 11 Number of leading informations (at lest 7)
     #! Trailing data at end of each visibility may hold additional information
     #integer(kind=4) :: ntrail = 0          ! 12 Number of trailing informations
     #!
     #! Leading / Trailing information codes have been specified before
     #integer(kind=4) :: column_pointer(code_uvt_last) = code_null ! Back pointer to the columns...
     #integer(kind=4) :: column_size(code_uvt_last) = 0  ! Number of columns for each
     #! In the data, we instead have the codes for each column
     #! integer(kind=4) :: column_codes(nlead+ntrail)         ! Start column for each ...
     #! integer(kind=4) :: column_types(nlead+ntrail) /0,1,2/ ! Number of columns for each: 1 real*4, 2 real*8
     #! Leading / Trailing information codes
     #!
     #integer(kind=4) :: order = 0          ! 13  Stoke/Channel ordering
     #integer(kind=4) :: nfreq = 0          ! 14  ! 0 or = nchan*nstokes
     #integer(kind=4) :: atoms(4)           ! 15-18 Atom description
     #!
     #real(kind=8), pointer :: freqs(:) => null()     ! (nchan*nstokes) = 0d0
     #integer(kind=4), pointer :: stokes(:) => null() ! (nchan*nstokes) or (nstokes) = code_stoke
     #!
     #real(kind=8), pointer :: ref(:) => null()
     #real(kind=8), pointer :: val(:) => null()
     #real(kind=8), pointer :: inc(:) => null()

    lf.seek(1024)
    real_dims = dims[:ndim]
    data = np.fromfile(lf, count=np.prod(real_dims),
                       dtype='float32').reshape(real_dims[::-1])
    data[data==bval] = np.nan

    return data,header


io_registry.register_reader('lmv', BaseSpectralCube, load_lmv_cube)
io_registry.register_reader('class_lmv', BaseSpectralCube, load_lmv_cube)
io_registry.register_identifier('lmv', BaseSpectralCube, is_lmv)