# $Id: datalib.py,v 1.8 2020/06/20 17:53:06 andrei Exp $
# -*- coding: koi8-r -*-
"""First attempt to systematize some functions for accessing and analyzing data
available either locally or on the internet.
"""
from __future__ import print_function
from __future__ import division

# some prerequisites
from builtins import chr
from builtins import range
from builtins import object
#from past.utils import old_div
import sys
import os
import glob
import time
import datetime as dt
import warnings as warn
from string import Template as tpl
import math
import bisect
import numpy as np
from scipy.interpolate import interp1d

# CDF library path
os.putenv("CDF_LIB", "/home/andrei/cdf/lib")
import spacepy.datamodel as dm
from spacepy import pycdf
import spacepy.time as spt

def dataroot(prj=None):
  """Default local data root directory for all projects"""
  if prj is None:
    prj = ''
  prjl = prj.lower()
  if 'pwing' in prjl or 'isee' in prjl:
    root = '/media/andrei/My Passport/data'
  else:
    root = '/home/andrei/data'
  return root

def known_projects():
  """Returns a list or set of known project names"""
# ÄÒÕÇÏÊ ÐÏÄÈÏÄ: ÐÏÓÍÏÔÒÅÔØ ×ÓÅ ËÁÔÁÌÏÇÉ × dataroot?
# ôÒÅÔÉÊ ÐÏÄÈÏÄ: ÐÒÏÞÉÔÁÔØ ÓÐÅÃÉÁÌØÎÙÊ ÆÁÊÌ (× dataroot)?
  proj_set = ["CLUSTER", "THEMIS", "RBSP", "MMS", 'ERG', 'SGO', 'ISEE', 'PWING', 'NOAA', 'EISCAT']
  return set(proj_set)

def is_known(prnam):
  """Tells whether the specified name denotes a known project"""
  return (prnam.upper() in known_projects())

def datadir(prj=None, drng=None, root=None, **kw):
  """Data directory structure/formats for local files in project prj for date drng; to be replaced by reading a configuration file?
  
  The following keyword arguments will be recognized:
  
  root: root local directory for all data (absolute)

  pdir: project data directory (relative to root)
  
  sat: spacecraft or project station
  
  suite: suite of instruments

  instr: the instrument
  
  prd: product

  coord: coordinate system

  ftyp: file type, e.g., 'cdf', 'h5'
  
  ldirtmpl: template for local data directory (relative to root/pdir)

  rdirtmpl: template for local data directory (relative to remdir)

  remdir: remote root directory for given type of data (absolute)

  """
  # recarg = {'prj': prj, 'drng': drng}
  # for a in recarg:
  #   if recarg[a] is None and a in kw:
  if prj is None:
    if 'prj' in kw:
      prj = kw['prj']
    else:
      warn.warn('Project must be specified')
      return None
  dr = dataroot(prj)
  ds = dict(root=dr)
  if drng is None:
    if 'drng' in kw:
      drng = kw['drng']
    else:
      warn.warn('Date/time range must be specified')
      return None
  if(is_known(prj.upper())):
    prjl = prj.lower()
    ds['prj'] = prjl
    if(prjl in ('eiscat', 'mms', 'rbsp', 'themis', 'erg', 'sgo', 'pwing', 'pwing')):
      ds['pdir'] = prj.upper()
    elif(prjl == 'cluster'):
      ds['pdir'] = 'Cluster'
    elif(prjl in ('noaa', 'isee')):
      ds['pdir'] = prjl
  else:
    warn.warn(prj + ': Unknown project')
    return ds
  # ÓÐÉÓÏË ÄÏÐÕÓÔÉÍÙÈ ËÌÀÞÅÊ
  valid_kw = ['sat', 'suite', 'instr', 'prd', 'coord', 'lev', 'rel', 'ver', 'ftyp']
  # ÓÐÉÓÏË ÎÅÏÂÈÏÄÉÍÙÈ ËÌÀÞÅÊ
  req_kw = []
  # ÓÐÉÓÏË ËÌÀÞÅÊ, ÄÌÑ ËÏÔÏÒÙÈ (ÍÏÇÕÔ ÂÙÔØ) ÎÕÖÎÙ ÚÁÇÌÁ×ÎÙÅ ×ÁÒÉÁÎÔÙ
  cap_kw = ['sat', 'suite', 'instr']
  # ÓÐÉÓÏË ËÌÀÞÅÊ, ËÏÔÏÒÙÅ ÎÅ ÎÕÖÎÏ ÐÒÅÏÂÒÁÚÏ×Ù×ÁÔØ
  asis_kw = ['prd', 'coord']
  for k in list(kw):
    if k in valid_kw:
      if kw[k] is not None:
        if str is bytes:
          ctup = (str, unicode) #py2k
        else:
          ctup = str
        if isinstance(kw[k], ctup) and (k not in asis_kw):
          ds[k] = kw[k].lower()
          if k in cap_kw:
            ds[k.upper()] = kw[k].upper()
        else:
          ds[k] = kw[k]
#          print ('kw[{}] = {}'.format(k, kw[k]))
  for k in req_kw:
    if k not in ds:
      ds[k] = ''
  ### ÐÏ ÕÍÏÌÞÁÎÉÀ ÉÝÅÍ ÎÁ CDAWeb
   site = 'http://cdaweb.gsfc.nasa.gov/istp_public/data'
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ RBSP
  if prjl == 'rbsp':
    ds['ftyp'] = 'cdf'
    ### ÐÒÅÄ×ÁÒÉÔÅÌØÎÙÅ ÐÒÏ×ÅÒËÉ É ÎÁÚÎÁÞÅÎÉÅ ÎÅÎÁÚÎÁÞÅÎÎÏÇÏ...
    if ds['instr'] == 'hope' or ds['instr'] == 'mageis':
      if 'rel' not in kw:
        ds['rel'] = 4
        warn.warn('Using default release {} for {}'.format(ds['rel'], ds['INSTR']))
      ds['rel'] = '{:02}'.format(ds['rel'])
      if 'lev' not in kw:
        ds['lev'] = 3
        warn.warn('Using default level {} for {}'.format(ds['lev'], ds['INSTR']))
      if ds['instr'] == 'hope' and 'prd' not in kw:
        ds['prd'] = 'pa'
        ds['PRD'] = 'PA'
        warn.warn('Using default product {} for {}'.format(ds['prd'], ds['INSTR']))
      if 'suite' not in kw or kw['suite'] != 'ect':
        ds['suite'] = 'ect'
        warn.warn('Using default suite {} for {}'.format(ds['suite'], ds['INSTR']))
    if ds['instr'] == 'emfisis':
      site = 'https://emfisis.physics.uiowa.edu/Flight'
      rdirtmps = '$site/$PRJ-$SAT/L$lev/$Y/$mon/$day'
      ldirtmps = '$root/$pdir/$SAT/L$lev/$Y/$mon/$day'
      # ÎÅ ×ÓÅÇÄÁ ×ÅÒÎÏ! îÏ ÎÁ ÐÒÁËÔÉËÅ ÐÏËÁ ÇÏÄÉÔÓÑ
      if 'prd' in list(ds):
        prd = ds['prd']
        if 'lev' not in ds and ('WFR-spectral-matrix' in prd or 'continuous' in prd or 'HFR' in prd):
          ds['lev'] = 2
        elif 'magnetometer' in prd:
          if 'uvw' in prd:
            ds['lev'] = 2
          else:
            ds['lev'] = 3
        elif 'svd' in prd or 'dcf' in prd or 'density' in prd or 'NURD' in prd or 'HFR-UHR' in prd:
          ds['lev'] = 4
          # ðÏËÁ ÐÏ ÕÍÏÌÞÁÎÉÀ ÇÒÕÚÉÍ ÆÁÊÌÙ, ÐÏÌÕÞÅÎÎÙÅ Ó ÐÏÍÏÝØÀ octave
          if 'svd' in prd:
            if 'ftyp' not in kw:
              ds['ftyp'] = 'h5'
              ds['ver'] = None
              if 'svd*' not in prd:
                ds['prd'] = ds['prd'].replace('svd', 'svd*')
          if 'HFR-UHR' in prd:
            ds['ftyp'] = 'txt'
        if 'waveform' in prd:
          filetmps = '$prj-${sat}_${prd}_$instr-L${lev}_$Y$mon${day}T${H}'
        elif 'NURD' in prd:
          ldirtmps = '$root/$pdir/$SAT/L$lev/$Y'
          filetmps = '$prj-${sat}_$Y$mon${day}_${prd}_$instr-L${lev}'
        else:
          filetmps = '$prj-${sat}_${prd}_$instr-L${lev}_$Y$mon${day}'
    elif ds['instr'] == 'efw':
      site = 'http://rbsp.space.umn.edu/data'
      if 'lev' not in kw:
        ds['lev'] = 3
        warn.warn('Using default level {} for {}'.format(ds['lev'], ds['INSTR']))
      if ds['lev'] == 3:
        rdirtmps = '$site/$prj/${prj}${sat}/l$lev/$Y'
        ldirtmps = '$root/$pdir/$SAT/L$lev/$Y'
        filetmps = '${prj}${sat}_$instr-[Ll]${lev}_$Y$mon${day}'
    elif ds['instr'] == 'hope':
      if ds['PRD'] == 'PA':
        ds['remprd'] = 'pitchangle'
        rdirtmps = '$site/$prj/$prj$sat/l$lev/$suite/$instr/$remprd'
      else:
        rdirtmps = '$site/$prj/$prj$sat/l$lev/$suite/$instr/$prd'
      rdirtmps = rdirtmps + 'rel$rel/$Y'
      ldirtmps = '$root/$pdir/$SAT/$suite/$instr/l$lev/$prd'
      filetmps = '${prj}${sat}_rel${rel}_${suite}-$instr-[${PRD}${prd}]*-[Ll]${lev}_$Y$mon${day}'
    elif ds['instr'] == 'mageis':
      rdirtmps = '$site/$prj/$prj$sat/l$lev/$suite/$instr/$prd/rel$rel/$Y'
      ldirtmps = '$root/$pdir/$SAT/$suite/$instr/level$lev'
      filetmps = '${prj}${sat}_rel${rel}_${suite}-${instr}-[Ll]${lev}_$Y$mon${day}'
    elif ds['instr'] == 'ephem':
      site = 'https://emfisis.physics.uiowa.edu/Flight'
      rdirtmps = '$site/$PRJ-$SAT/LANL/MagEphem/$Y'
      ldirtmps = '$root/$pdir/$SAT/MagEphem/$Y'
      filetmps = '$prj${sat}_def_MagEphem_${prd}_$Y$mon$day'
      ds['ftyp'] = 'h5'
      if 'prd' not in list(ds):
        ds['prd'] = 'TS04D'
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ ERG
  elif prjl == 'erg':
    if 'ftyp' not in ds:
      ds['ftyp'] = 'cdf'
    site = 'https://ergsc.isee.nagoya-u.ac.jp/data/ergsc/satellite'
    filetmps = '${prj}_${instr}_l${lev}_${prd}_$Y$mon${day}'
    rdirtmps = '$site/$prj/$instr/l$lev/$prd/$Y/$mon'
    ldirtmps = '$root/$pdir/$instr/l$lev/$prd/$Y'
    if 'suite' in list(ds.keys()) and ds['suite'] == 'pwe':
     ldirtmps = '$root/$pdir/$suite/$instr/l$lev/$Y/$mon'
      rdirtmps = '$site/$prj/$suite/$instr/l$lev/$prd/$Y/$mon'
      filetmps = '${prj}_${suite}_${instr}_l${lev}_${prd}_$Y$mon${day}'
    if ds['prd'] == '8sec':
      ds['lev'] = 2
    elif '64hz' in ds['prd'] and 'polariz' not in ds['prd']:
      if 'coord' not in ds:
        raise KeyError('Coord system must be given for hires {} data'.format(kw['instr']))
      ds['lev'] = 2
      ds['baseprd'] = '64hz'
      if 'suite' in ds and len(ds['suite']):
        ldirtmps = '$root/$pdir/$suite/$instr/l$lev/$baseprd/$Y/$mon'
        filetmps = '${prj}_${suite}_${instr}_l${lev}_${prd}_${coord}_$Y$mon${day}'
      else:
        ldirtmps = '$root/$pdir/$instr/l$lev/$baseprd/$Y/$mon'
        filetmps = '${prj}_${instr}_l${lev}_${prd}_${coord}_$Y$mon${day}${H}'
    if 'svd' in ds['prd'] or 'polariz' in ds['prd']:
      ds['lev'] = 4
      if 'instr' in ds:
        filetmps = '${prj}_${instr}_l${lev}_${prd}_$Y$mon${day}T${H}'
      else:
        ldirtmps = '$root/$pdir/l$lev/$prd/$Y/$mon'
        filetmps = '${prj}_l${lev}_${prd}_$Y$mon${day}T${H}'
      if 'ftyp' not in kw:
        ds['ftyp'] = 'h5'
        ds['ver'] = None
        if 'svd*' not in ds['prd']:
          ds['prd'] = ds['prd'].replace('svd', 'svd*')
    elif ds['prd'] == 'orb':
      rdirtmps = '$site/$prj/$prd/l$lev/$Y'
      ldirtmps = '$root/$pdir/$prd/l$lev/$Y'
      filetmps = '${prj}_${prd}_l${lev}_$Y$mon${day}'
    elif 'ep' in ds['instr']:
      if ds['prd'] == 'omniflux':
        ds['lev'] = 2
      elif 'pad' in ds['prd'] or 'pa' in ds['prd']:
        ds['lev'] = 3
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ NOAA
  elif prjl == 'noaa':
    if 'sat' not in kw:
      raise TypeError('sat num must be given')
    # NB! Should check for valid dates and satellite numbers
    if 'ftyp' not in kw:
      ds['ftyp'] = 'cdf'
    if 'instr' not in kw:
      ds['instr'] = 'poes'
    if 'prd' not in kw:
      ds['prd'] = 'sem2_fluxes-2sec'
    filetmps = '${prj}${sat}_${instr}-${prd}_$Y$mon${day}'
    ldirtmps = '$root/$pdir/$pdir$sat/$Y'
    rdirtmps = '$root/$pdir/$pdir$sat/$prd/$Y'
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ Cluster
  elif prjl == 'cluster':
    if 'ftyp' not in kw:
      ds['ftyp'] = 'cdf'
    filetmps = '${SAT}_${INSTR}-${prd}_$Y$mon${day}'
    ldirtmps = '$root/$pdir/$INSTR/$Y'
    rdirtmps = '$root/$pdir/$pdir$sat/$prd/$Y'
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ SGO(VLF)
  elif prjl == 'sgo':
    if 'ftyp' not in kw:
      ds['ftyp'] = 'mat'
    if 'instr' not in kw:
      ds['instr'] = 'vlf'
    ldirtmps = '$root/$pdir/$instr/$ftyp/$Y$mon${day}'
    rdirtmps = ldirtmps
    if 'ffmt' not in kw:
      filetmps = 'comp${comp}Matrix$Y$mon${day}$H$M$S_D$Y$mon${day}_ql1_${dur}${tu}'
    else:
      filetmps = kw['ffmt']
  ### óÔÒÕËÔÕÒÁ ÄÁÎÎÙÈ ÄÌÑ PWING(VLF)
  ### prj=isee/pwing, suite=VLF, instr=stn, prd=40k/spec
  ### $droot/$prj/$suite/$instr/$stn/$prd/$Y
  ### VLFNo_yyyymmddhhnnss_FRQk_CH.STN
  elif prjl == 'pwing' or prjl == 'isee':
    prj = 'isee'
    if 'ftyp' not in kw:
      ds['ftyp'] = 'cdf'
    if 'suite' not in kw:
      ds['suite'] = 'VLF'
    filetmps = '${suite}03_$Y$mon${day}${H}${M}${S}_${prd}_03'
    ldirtmps = '$root/$pdir/$pdir$sat/$Y'
    rdirtmps = '$root/$pdir/$pdir$sat/$prd/$Y'
  ### ÐÒÏÅËÔ ÎÅ ××ÅÄÅÎ; ÐÏÐÒÏÂÏ×ÁÔØ CDAWEB???
  else:
    site = None
    ldirtmpl = tpl('')
    filetmpl = tpl('')
    rdirtmpl = tpl('')
    warn.warn('Unknown project? ds={}'.format(ds))
    pass
  ############################
  ds['site'] = site
  if 'debug' in kw and kw['debug']:
    print('ds =', ds)
    print('rdirtmps =', rdirtmps)
    print('ldirtmps =', ldirtmps)
    print('filetmps =', filetmps)
  ds['rdirtmpl'] = tpl(rdirtmps)
  ds['ldirtmpl'] = tpl(ldirtmps)
  if 'ver' in ds and ds['ver']:
    filetmps = filetmps + '_v${ver}'
  ds['filetmpl'] = tpl(filetmps)
  if drng is None:
    if 'drng' in kw:
      drng = kw['drng']
    else:
      warn.warn('Date empty, assume a limited output')
      return ds
  tmr = datrng2t(drng)
  ts = tmr['tstart']
  ds['Y'] = ts.tm_year
  ds['mon'] = '{:02}'.format(ts.tm_mon)
  ds['day'] = '{:02}'.format(ts.tm_mday)
  ds['H'] = '{:02}'.format(ts.tm_hour)
  ds['M'] = '{:02}'.format(ts.tm_min)
  ds['S'] = '{:02}'.format(ts.tm_sec)
  if 'debug' in kw and kw['debug']:
    print('ds =', ds)
  return ds

def fload(prj=None, drng=None, ds=None, **kw):
    """Load a data file by patterns in ds OR create a new one

    No checking of ds yet!!!
    """
    if ds is None:
      ds = datadir(prj, drng, **kw)
      if ds is None:
        warn.warn('Data structure undefined')
        return None
    dn = ds['ldirtmpl'].safe_substitute(ds)
    fp = ds['filetmpl'].safe_substitute(ds)
    if 'ver' in ds and ds['ver']:
      ptn = glob.glob('{}/{}.{}'.format(dn, fp, ds['ftyp']))
    else:
      ptn = glob.glob('{}/{}*.{}'.format(dn, fp, ds['ftyp']))
    if 'debug' in kw and kw['debug']:
      print('dn = ' + dn)
      print('fp = ' + fp)
      print('ptn: ', ptn)
    # ÐÒÏ×ÅÒËÁ ÆÁÊÌÁ
    if len(ptn):
      ptn.sort()
      # ÅÓÌÉ ÎÕÖÎÏ ÚÁÇÒÕÚÉÔØ ×ÓÅ ÆÁÊÌÙ ÐÏ ÛÁÂÌÏÎÕ:
      if 'loadall' in kw and kw['loadall']:
        fn = ptn
      else: # ÂÅÒÅÍ Ó ËÏÎÃÁ - ÄÕÍÁÅÍ, ÞÔÏ ÜÔÏ ÓÁÍÙÊ ÎÏ×ÙÊ
        fn = [ptn[-1]]
    elif 'new' in kw and kw['new']:
      if not os.path.exists(dn):
        try:
          os.makedirs(dn)
        except:
          warn.warn('Creating directory {} failed'.format(dn))
          return None
      if 'create' in kw and kw['create']:
        if ds['ftyp'] == 'cdf':
          fn = '{}/{}.cdf'.format(dn, fp)
          if os.path.exists(fn):
            if 'overwrite' in kw and kw['overwrite']:
              os.remove(fn)
            else:
              warn.warn('File exists {}'.format(fn))
              return None
          if 'skel' in kw:
            skel = kw['skel']
          else:
            skel = ''
          ### ÓÏÚÄÁÅÍ ÆÁÊÌ 
          pycdf.lib.set_backward(False)
          df = pycdf.CDF(fn, skel)
          df.compress(pycdf.const.GZIP_COMPRESSION, 5)
          return df
        else:
          warn.warn('Creation of .{} file not supported'.format(ds['ftyp']))
          raise TypeError
      else:
        warn.warn('Directory {} created'.format(dn))
        return None
    else:
      return None
    dfl = []
    for i in range(len(fn)):
      dfl.append(dataload(fn[i], **kw))
    if 'debug' in kw and kw['debug']:
      print('{} files in the list'.format(len(fn)))
    if len(dfl) == 1:
      df = dfl[0]
    else:
      df = dfl
    return df

def dataload(fi, **kw):
  """Load data by filename, single name is expected"""
  fn, ext = os.path.splitext(fi)
  if ext == '.cdf':
      df = pycdf.CDF(fi)
      if 'overwrite' in kw and kw['overwrite']:
        df.readonly(False)
        if 'debug' in kw and kw['debug']:
          warn.warn('Overwriting: {}'.format(fn))
  elif ext in ('.h5', '.nc'):
    df = dm.fromHDF5(fi)
  elif ext in ('.sav', '.tplot'):
    from scipy import io as sio
    df = sio.readsav(fi)
  else:
    warn.warn('Unknown file type .{}'.format(ext))
  return df

def datrng2t(drng, fmt='%Y%m%dT%H:%M', sep='-'):
  """Just convert date range to time structures
  """
  d = drng.partition(sep)
  if d[1] != sep:
    warn.warn('{} not found in {}?'.format(sep, drng))
  tms = time.strptime(d[0], fmt)
  r = dict(tstart=tms)
  if len(d[2]):
    de = '{0:.9s}{1}'.format(d[0], d[2])
    tme = time.strptime(de, fmt)
    r['tend'] = tme
  return r

def drng2dt(drng, fmt='%Y%m%dT%H:%M', sep='-', ret_t='dt', debug=0):
  """Convert date range to datetime objects
  """
  d = drng.partition(sep)
  if d[1] != sep and debug:
    warn.warn('{} not found in {}?'.format(sep, drng))
  if ret_t == 'dt':
    fret = dt.datetime.strptime
  else:
    fret = time.strptime
  tms = fret(d[0], fmt)
  r = dict(start=tms)
  if len(d[2]):
    de = '{0:.9s}{1}'.format(d[0], d[2])
    tme = fret(de, fmt)
    r['end'] = tme
  return r

def dr2itr(epoch, dtr):
    """Find the indices of epoch (datetime) array enclosing given interval
    """
    if isinstance(dtr, str):
      dtrng = list(drng2dt(dtr).values())
    elif isinstance(dtr, list):
      dtrng = dtr
    elif isinstance(dtr, dict):
      dtrng = list(dtr.values())
    if isinstance(epoch[0], dt.datetime):
      tt = epoch
    elif isinstance(epoch[0], float):
      tt = spt.Ticktock(epoch, 'UNX')
    else:
      tt = spt.Ticktock(epoch).UTC
    istart = bisect.bisect_left(tt, dtrng[0])
    if istart:
      istart -= 1
    itr = dict(start=istart)
    iend = bisect.bisect_right(tt, dtrng[1])
    if iend < len(epoch) - 1:
      iend += 1
    itr['end'] = iend
    return itr

def set_tlim(drng, ax, fmt='%Y%m%dT%H:%M', sep='-'):
  """Set datetime limits based on drng
  """
  dtr = drng2dt(drng, fmt, sep)
  ts = dtr['start']
  tend = dtr['end']
  ax.set_xlim(ts, tend)

# ÄÌÑ ÉÎÔÅÒÐÏÌÑÃÉÉ ÎÕÖÎÙ "ÏÂÙÞÎÙÅ" ÓÅËÕÎÄÙ
def ep2t(epoch, ep0=None):
  """Convert epoch times to float seconds, e.g., for interpolation

  epoch can be any type convertible to spacepy.time.Ticktock 
  """
  tt = spt.Ticktock(epoch)
  t = tt.UNX
  if ep0 is not None:
    tt0 = spt.Ticktock(ep0)
    ep0 = tt0.UTC[0].replace(hour=0, minute=0, second=0, microsecond=0)
    t -= spt.Ticktock(ep0).UNX
  return t

def ep2dt(epoch, fmt='%Y-%m-%dT%H:%M:%S'):
  """Convert epoch times to datetime values, e.g., for plotting
  """
  nt = epoch.shape[0]
  d = np.zeros(nt, dtype='O')
  for i in range(nt):
    d[i] = dt.datetime.strptime(epoch[i], fmt)
  return d

## éÎÔÅÒÐÏÌÉÒÕÅÍ ÉÚÍÅÒÅÎÎÙÅ ×ÅÌÉÞÉÎÙ (ÍÁÇÎÉÔÎÏÅ ÐÏÌÅ) ÎÁ ÎÕÖÎÕÀ ÓÅÔËÕ ×ÒÅÍÅÎ
def ep_interp(ep0, v0, ep, kind='cubic', debug=0):
    """ep0[nt0] - Epochs for which v0 is known;

    v0(nt0, ncomp) - initial value as 2D matrix

    ep - Epochs to which we want to interpolate each row of v0

    No type checking at the moment
    
    ep0 and ep can be any types convertible to spacepy.time.Ticktock but NOT
    Ticktock itself!!

    """
    if debug:
      start = time.time()
      print('Interpolating...')
    if debug:
      end = time.time()
      wt1 = end-start
      start = end
    t0 = ep2t(ep0)
    t = ep2t(ep)
    if debug:
      end = time.time()
      wt2 = end-start
      start = end
    nt = ep.shape[0]
    ff = interp1d(t0, v0, kind=kind, axis=0, fill_value='extrapolate')
    v = ff(t)
    if debug:
      end = time.time()
      wt3 = end-start
      print('t1 = {}, t2 = {}, t3 = {}'.format(wt1, wt2, wt3))
    return v

def h5ep(d, v='Epoch', asis=False):
    """
    ÷ÏÓÓÔÁÎÏ×ÉÍ ÓÔÒÏÞÎÕÀ ÐÅÒÅÍÅÎÎÕÀ ÉÚ ËÒÉ×ÏÇÏ ÉÍÐÏÒÔÁ HDF5 ÉÚ Octave

    d - ÏÂßÅËÔ SpaceData, × ËÏÔÏÒÏÍ ÎÁÄÏ ÎÁÊÔÉ v (='Epoch') É ×ÏÓÓÔÁÎÏ×ÉÔØ
    """
    if v not in list(d):
        raise TypeError
    ee = d[v]['value']
    if ee.shape[1] > 1:
      ee = ee[:,0]
    l = len(ee)
    r = ''
    for ii in list(range(l)):
      c = ee[ii] #[0]
      r = r + chr(c)
    if asis:
      return r
    else:
      t = dt.datetime.strptime(str.rstrip(str(r)), '%Y-%m-%dT%H:%M:%S.%f')
      return t

## ÁÔÔÒÉÂÕÔÙ CDF ÆÁÊÌÏ× ÐÏ ÕÍÏÌÞÁÎÉÀ
def defattr(a, **kw):
  """Default attribute dict: add it to a"""
  if type(a) is not dict:
    raise TypeError('input must be dict')
  a['Author'] = 'Andrei Demekhov'
  a['CreateDate'] = dt.datetime.now().strftime('%Y%m%d')
  a['Generated_by'] = 'Polar Geophysical Institute'
  a['Generation_date'] = dt.datetime.now().strftime('%Y%m%d')
  a['Rules_of_use'].append('For collaborative use only: contact the author')
  a['TEXT'].append('Generated by a Python script for internal use')
  # valid_kw = ['module']
  # for k in valid_kw:
  #   if k in kw:
  return a

## ×ÙÞÉÓÌÉÍ ÐÉÔÞ-ÕÇÌÙ ÄÌÑ ËÁÖÄÏÊ ÔÏÞËÉ, ÐÅÒÅÍÎÏÖÁÑ ÍÁÔÒÉÃÙ ÐÏ×ÏÒÏÔÏ×
def erg_gse2pa(df, mcdf=None, *args):
    """Calculate pitch-angles from 3dflux gse data

    Input data will be just 3dflux CDF and magnetic field CDF [and optionally
    datetime range, not implemented]

    Assume that the MF resolution is known and is 8 seconds. 

    No variable check yet.
    """
    dn, fname = os.path.split(df.pathname)
    fn, ext = os.path.splitext(fname)
    fn = str(fn)
    if mcdf is None:
      #try to find magnetic field file based on df - not implemented yet
      pass
    print('Pitch-angle calculation...')
    start = time.time()
    if 'mepe' in fn:
      anam = 'FEDU_Angle_gse'
    elif 'lepe' in fn:
      anam = 'FEDU_Angle_GSE'
    else: # 'mepi' in fn or  'lepi' in fn
      anam = 'FIDU_Angle_gse'
    if anam not in df:
      raise KeyError
    if 'lepi' not in fn:
      agse = df[anam]
    else:
      ### calculations required for Angle_gse!!!
      raise TypeError
    if 'epoch' in df:
      ep = df['epoch']
    elif 'Epoch' in df:
      ep = df['Epoch']
    else:
      raise KeyError
    bgse = mcdf['mag_8sec_gse']
    epm = mcdf['epoch_8sec']
    b = ep_interp(epm, bgse, ep)
    btm = mcdf['magt_8sec']
    bt = ep_interp(epm, btm, ep)
    nt = ep.shape[0]
    # agse[nt, 2, n_detect, nspin]
    napd = agse.shape[2]
    nspin = agse.shape[3]
    pa = np.zeros((nt, nspin, napd))
    ii = 0.0
    for i in range(nt):
        ct = b[i,2] / bt[i]
        st = np.sqrt(1 - ct*ct)
        phi = np.arctan2(b[i,1], b[i,0])
        cphi = np.cos(phi)
        sphi = np.sin(phi)
        th_mat = np.array([[ct, 0, -st], [0, 1, 0], [st, 0, ct]])
        phi_mat = np.array([[cphi, sphi, 0], [-sphi, cphi, 0], [0, 0, 1]])
        T = np.matmul(th_mat, phi_mat)
        rad = np.deg2rad(agse[i])
        ct = np.sin(rad[0])
        st = np.cos(rad[0])
        cphi = np.cos(rad[1])
        sphi = np.sin(rad[1])
        for ia in range(napd):
          for isp in range(nspin):
            v0 = np.array((st[ia,isp]*cphi[ia,isp],
                           st[ia,isp]*sphi[ia,isp],
                           ct[ia,isp]))
            v = np.matmul(T, v0)
            pa[i, isp, ia] = np.arccos(v[2])
        ii = ii + 5.0/nt
        if(int(ii)):
          ip = float(i*100) / nt
          print('{}% i={}'.format(ip, i))
          ii = 0
    end = time.time()
    print('Pitch-angle calculation time = {}'.format(end-start))
    return np.rad2deg(pa)

def erg_ep_pad(df, mcdf=None, pa=None, npa=36, **kw):
    """Calculate the MEP(e,i)/LEP(e,i) PAD from ERG 3dflux data

    Valid keywords: (drng, save:format/bool); 
      let = which particle to treat
    """
    # first we need to check which variables we need (depends on the detector)
    fn = str(df.pathname)
    if 'epe' in fn:
      flets = ['e']
    elif 'epi' in fn:
      flets = ['p', 'o', 'opp', 'o2p', 'he', 'he2']
    # ÔÉÐ ÞÁÓÔÉÃÙ ÚÁÄÁÎ kw['let'] ÉÌÉ ÐÏ ÕÍÏÌÞÁÎÉÀ (ÉÍÑ ÐÒÉÂÏÒÁ, e/p)
    if 'let' in kw:
      if kw['let'] in flets:
        plet = kw['let']
      else:
        raise KeyError
    else:
      plet = flets[0]
    f = 'F{}DU'.format(plet.upper())
    fe = 'F{}EDU'.format(plet.upper())
    e = 'F{}DU_Energy'.format(plet.upper())
    # check for the needed variables; Energy just for checking
    for k in [f, fe, e]:
        if k not in list(df):
            warn.warn('No {} in {}'.format(k, df.pathname))
            if k != fe:
              raise KeyError
    if pa is None:
        # ÎÁÄÏ ×ÙÞÉÓÌÉÔØ pa Ó ÐÏÍÏÝØÀ gse2pa
        pa = erg_gse2pa(df, mcdf)
    start = time.time()
    # ÄÁÌÅÅ pa ÕÖÅ ÅÓÔØ
    print('Making {} pitch-angle distribution...'.format(f))
    nt = pa.shape[0]
    if 'mep' in fn: # (nt, nspin(sectors), nE, napd(anodes))
      flux = df[f]
      eflux = df[fe]
    # NB! ÄÌÑ ÉÏÎÏ× ÔÏÖÅ!
    elif 'lep' in fn: # (nt, nE, napd(anodes), nspin(sectors))
      if 'lepi' in fn and 'lepinarrowchannels' in kw:
        fsub = f + '_sub'
        fesub = fe + '_sub'
        f0 = np.concatenate((df[f], df[fsub]), axis=2)
        if fe in list(df):
          fe0 = np.concatenate((df[fe], df[fesub]), axis=2)
      else:
        for k in [f, fe]:
          if k in list(df.keys()):
            f0 = df[k]
      flux = np.moveaxis(f0, [1, 2, -1], [-2, -1, 1])
      if fe in list(df):
        eflux = np.moveaxis(fe0, [1, 2, -1], [-2, -1, 1])
    nspin = flux.shape[1]
    nE = flux.shape[2]
    bins = np.linspace(0, 180, npa+1)
    pad = np.zeros((nt, nE, npa))
    if fe in list(df):
      epad = np.zeros((nt, nE, npa))
    else:
      epad = 0
    pa_bins = np.digitize(pa, bins) - 1
    ii = 0.0
    for ia in range(npa):
      ii = ii + 6.0/npa
      if(int(ii)):
        ip = float(ia*100) / npa
        print('{}% i={}'.format(ip, ia))
        ii = 0
      msk = pa_bins == ia
      for iE in range(nE):
        m1 = msk & (flux[:, :, iE, :] > 0)
        ff = np.ma.array(flux[:, :, iE, :], mask = ~m1)
        pad[:, iE, ia] = np.nanmean(ff, axis=(1,2))
        if fe in list(df):
          m1 = msk & (eflux[:, :, iE, :] > 0)
          ff = np.ma.array(eflux[:, :, iE, :], mask = ~m1)
          epad[:, iE, ia] = np.nanmean(ff, axis=(1,2))
    end = time.time()
    print('PAD calculation time = {}'.format(end-start))
    if 'save' in kw:
      if type(kw['save']) is bool:
        sext = 'cdf'
      elif type(kw['save']) is str:
        sext = kw['save']
    return pad, epad, bins

### óÏÈÒÁÎÉÍ ðõò ÉÌÉ ÐÏÔÏË Ó ÐÉÔÞ-ÕÇÌÁÍÉ (ÉÌÉ ÔÏÌØËÏ ÐÉÔÞ-ÕÇÌÙ) × CDF
def addpad(df3d, dfm, dfpad, let='e', pa=None, npa=36, mkpad=True):
    """Add the calculated MEP(e,i)/LEP(e,i) PAD from ERG 3dflux data 
    to the CDF PA file, or just save pitch-angle data (not working)
    
      let = which particle to treat
    """
    k = 'F{}DU'.format(let.upper())
    if k not in list(df3d.keys()):
      warn.warn('No {} in {}'.format(k, df3d.pathname))
      return
    if let == 'e':
        iORe = 'e'
    else:
        iORe = 'i'
    sfdu = 'F{}DU'.format(let.upper())
    sen = '{}_Energy'.format(sfdu)
    sfedu = 'F{}EDU'.format(let.upper())
    sangl = 'F{}DU_Angle'.format(iORe.upper())
    sanglab = '{}_LABL'.format(sangl)
## ÇÅÎÅÒÉÒÕÅÍ ðõò?
    if mkpad:
        pad, epad, bins = erg_ep_pad(df3d, dfm, pa, npa, **dict(let=let))
        vars4cp = [sen, '{}_LABL'.format(sen)]
        for v in vars4cp:
            if v in df3d:
                dfpad.clone(df3d[v])
    else:
        # vars4cp = ['FEDU', 'FEDU_Energy', 'FEDU_Energy_LABL',
        #            'FEDU_Sector', 'FEDU_Sector_LABL',
        #            'FEDU_Phase', 'FEDU_Phase_LABL', # LEPe
        #            'FEDU_APDno', 'FEDU_APD_LABL', 'epoch', 'Epoch',
        #            'epoch_sp', 'spin_phase', 'mode_reduction', 'Reduct_Flag']
        pass
        ### Now add the real things
    if 'lepe' in str(df3d.pathname.lower()):
      # make bin centers from bin margins
        senbins = '{}_bins'.format(sen)
        dfpad[sen].rename(senbins)
        dfpad[senbins].attrs['FIELDNAM'] = '{} bin margins'.format(sen)
        dfpad[sen] = 0.5 * (df3d[sen][:,0,:] + df3d[sen][:,1,:])
        dfpad[sen].attrs = dfpad[senbins].attrs
        dfpad[sen].attrs['FIELDNAM'] = '{} bin centers'.format(sen)
        e_un = 'eV'
    else:
        e_un = 'keV'
    un = 'cm!U-2!Ns!U-1!Nsr!U-1!N'
    #del dfpad['Epoch'].attrs['DISPLAY_TYPE']
    if mkpad:
        # Energy
        a2del = ['DISPLAY_TYPE', 'FILLVAL']
        for a in a2del:
          if a in dfpad[sen].attrs:
            del dfpad[sen].attrs[a]
        a4cp = ['AVG_TYPE', 'CATDESC', 'FIELDNAM', 'FILLVAL',
                'FORMAT', 'MONOTON', 'LABLAXIS', 'SCALEMAX','SCALEMIN',
                'SCALETYP', 'VALIDMAX', 'VALIDMIN', 'VAR_TYPE']
        dfpad[sfdu] = pad
        dfpad[sfdu].attrs['UNITS'] = '{}{}!U-1!N'.format(un, e_un)
        vlst = [sfdu]
        if sfedu in df3d:
            dfpad[sfedu] = epad
            dfpad[sfedu].attrs['UNITS'] = '{} {}{}!U-1!N'.format(e_un, un, e_un)
            vlst.append(sfedu)
        for v in vlst:
            for a in a4cp:
                if a in df3d[v].attrs:
                    dfpad[v].attrs[a] = df3d[v].attrs[a]
            dfpad[v].attrs['DISPLAY_TYPE'] = 'spectrogram'
            dfpad[v].attrs['DEPEND_0'] = 'Epoch'
            dfpad[v].attrs['DEPEND_1'] = sen
            dfpad[v].attrs['DEPEND_2'] = sangl
            if 'LABL_PTR_2' in df3d[v].attrs:
                dfpad[v].attrs['LABL_PTR_1'] = df3d[v].attrs['LABL_PTR_2']
            if sanglab in df3d[v].attrs:
                dfpad[v].attrs['LABL_PTR_2'] = sanglab
        if 'lepe' in str(df3d.pathname.lower()):
          pass
        # Angle
        dfpad[sangl] = bins[0:bins.size-1] + 0.5 * (bins[1]-bins[0])
        dfpad[sangl].attrs = dfpad[sen].attrs
        dfpad[sangl].attrs['CATDESC'] = 'Pitch angles for {}, middle of each bin'.format(sfdu)
        dfpad[sangl].attrs['FIELDNAM'] = 'Pitch angle'
        dfpad[sangl].attrs['LABLAXIS'] = 'Pitch angle'
        dfpad[sangl].attrs['MONOTON'] = 'INCREASE'
        dfpad[sangl].attrs['SCALEMAX'] = 180
        dfpad[sangl].attrs['UNITS'] = 'deg'
        dfpad[sangl].attrs['VALIDMAX'] = 180
        del dfpad[sangl].attrs['SCALETYP']
        # Angle labels
        dfpad[sanglab] = np.empty(bins.size-1, dtype='S7')
        for i in range(bins.size-1):
            dfpad[sanglab][i] = '{:>3d} deg'.format(int(bins[i]))
    else:
        dfpad[sangl] = pa
        a4cp = ['AVG_TYPE', 'FORMAT', 'MONOTON', 'SCALEMAX','SCALEMIN',
                'SCALETYP','UNITS', 'VALIDMAX', 'VALIDMIN']
        if 'FEDU_Angle_gse' in df3d:
            agse = df3d['FEDU_Angle_gse']
        if 'FEDU_Angle_GSE' in df3d:
            agse = df3d['FEDU_Angle_GSE']
        for a in a4cp:
            dfpad[sangl].attrs[a] = agse.attrs[a]
        dfpad[sangl].attrs['CATDESC'] = sfdu + ' pitch angles'
        dfpad[sangl].attrs['DEPEND_0'] = 'Epoch'
        dfpad[sangl].attrs['DEPEND_1'] = dfpad['FEDU'].attrs['DEPEND_1']
        dfpad[sangl].attrs['DEPEND_2'] = dfpad['FEDU'].attrs['DEPEND_3']
        dfpad[sangl].attrs['DISPLAY_TYPE'] = 'time_series'
        dfpad[sangl].attrs['FIELDNAM'] = sfdu + ' pitch angles'
        dfpad[sangl].attrs['FILLVAL'] = -1.0E+31
        dfpad[sangl].attrs['LABLAXIS'] = 'Pitch-angle'
        dfpad[sangl].attrs['VAR_TYPE'] = 'data'
        dfpad[sangl].attrs['VAR_NOTES'] = sfdu + ' pitch-angles calculated from GSE angles'

def erg_lepe_emerge(l3f, debug=0):
    """combine LEPe data for low energy-resolution mode"""
    ep = l3f['Epoch']
    e0 = l3f['FEDU_Energy']
    pad0 = l3f['FEDU']
    nt = e0.shape[0]
    nE = e0.shape[1]
    npa = pad0.shape[2]
    ## ÚÁÇÏÔÏ×ËÁ pad[nt, nE, npa]
    pad = np.zeros(pad0.shape)
    ## ÚÁÇÏÔÏ×ËÉ ÄÌÑ ÐÏÔÏËÏ× Ó ÌÉÎÅÊÎÏÊ É ÜËÓÐÏÎÅÎÃÉÁÌØÎÏÊ ÉÎÔÅÒÐÏÌÑÃÉÅÊ [nt, nE, npa]
    F_int = np.zeros(pad0.shape)
    F_logint = np.zeros(pad0.shape)
    ## ÚÁÇÏÔÏ×ËÁ energy[nt, nE]
    energy = np.zeros(e0.shape)
    e_int = np.zeros(e0.shape)
    igot_e = False
    if debug:
      start = time.time()
      print('Regularizing LEPe data...')
      ii = 0.0
    for it in range(nt):
        eax = e0[it]
        eu = np.unique(eax)
        nEu = len(eu)
        if not (nEu == 4 or nEu == 8):
            if not igot_e:
                igot_e = True
                ## ÎÁ ÜÔÕ ÏÓØ ÂÕÄÅÍ ÉÎÔÅÒÐÏÌÉÒÏ×ÁÔØ ÄÁÎÎÙÅ Ó ÎÉÚËÉÍ ÒÁÚÒÅÛÅÎÉÅÍ
                eax_int = eax
            pad[it] = pad0[it]
            F_int[it] = pad0[it]
#            F_logint[it] = pad0[it]
#            Fsum[it] = pad0[it]
            energy[it] = e0[it]
            e_int[it] = e0[it]
        else:
            # ÎÉÚËÏÅ ÒÁÚÒÅÛÅÎÉÅ
            pad1 = np.zeros((nEu, npa))
#            pad2 = np.zeros((nEu, npa))
            for iE in range(nEu):
                # ÚÁËÒÙ×ÁÅÍ ÎÅÎÕÖÎÙÅ ÚÎÁÞÅÎÉÑ ÜÎÅÒÇÉÉ
                msk = np.broadcast_to(eax != eu[iE], (npa, nE)).T
#                cc = pad0[it] == np.broadcast_to(np.nan, (nE, npa))
                cc = np.isnan(pad0[it])
                msk = msk | cc
                padm = np.ma.array(pad0[it, :, :], mask=msk)
                pad1[iE, :] = np.nanmean(padm, axis=0)
#                pad2[iE, :] = np.ma.sum(padm, axis=0)
            pad[it, list(range(nEu)), :] = pad1
#            Fsum[it, range(nEu), :] = pad2
            energy[it, list(range(nEu))] = eu
#            Fint = interp.PchipInterpolator(np.log(eu), np.log(pad1), axis=0)
#            F_logint[it] = np.exp(Fint(np.log(eax_int)))
            Fint = interp.PchipInterpolator(eu, pad1, axis=0)
            F_int[it] = Fint(eax_int)
            e_int[it] = eax_int
        if debug:
            ii = ii + 5.0/nt
            if(int(ii)):
                ip = float(it*100) / nt
                print('{}% i={}'.format(ip, it))
                ii = 0
    ## ÔÅÐÅÒØ ÚÁÍÅÎÉÍ ÕÓÒÅÄÅÎÎÎÙÅ ÄÁÎÎÙÅ ÎÉÚËÏÇÏ ÒÁÚÒÅÛÅÎÉÑ ÎÁ ÉÎÔÅÒÐÏÌÉÒÏ×ÁÎÎÙÅ
    if debug:
        end = time.time()
        print('Calculation time = {}'.format(end-start))
    if l3f.readonly():
        l3f.readonly(False)
    ss = ', sorted'
    if 'FEDU_srt' not in l3f:
        l3f['FEDU_srt'] = pad
        l3f['FEDU_srt'].attrs = l3f['FEDU'].attrs
        l3f['FEDU_srt'].attrs['CATDESC'] += ss
        l3f['FEDU_Energy_srt'] = energy
        l3f['FEDU_Energy_srt'].attrs = l3f['FEDU_Energy'].attrs
        l3f['FEDU_Energy_srt'].attrs['CATDESC'] += ss
        l3f['FEDU_srt'].attrs['DEPEND_1'] = 'FEDU_Energy_srt'
        l3f['FEDU_sum'] = Fsum
        l3f['FEDU_sum'].attrs = l3f['FEDU'].attrs
        l3f['FEDU_sum'].attrs['CATDESC'] += ', summed'
        l3f['FEDU_sum'].attrs['DEPEND_1'] = 'FEDU_Energy_srt'
    if 'FEDU_int' not in l3f:
        l3f['FEDU_int'] = F_int
        l3f['FEDU_int'].attrs = l3f['FEDU'].attrs
        l3f['FEDU_int'].attrs['CATDESC'] += ', interpolated'
        l3f['FEDU_logint'] = F_logint
        l3f['FEDU_logint'].attrs = l3f['FEDU'].attrs
        l3f['FEDU_logint'].attrs['CATDESC'] += ', logarithm interpolated'
        l3f['FEDU_Energy_int'] = e_int
        l3f['FEDU_Energy_int'].attrs = l3f['FEDU_Energy'].attrs
        l3f['FEDU_Energy_int'].attrs['CATDESC'] += ', interpolated'
        l3f['FEDU_int'].attrs['DEPEND_1'] = 'FEDU_Energy_int'
        l3f['FEDU_logint'].attrs['DEPEND_1'] = 'FEDU_Energy_int'

class Dataset(object):
  """
  A class specifying attributes for a given data set.
  
  Input characteristics: 

  date range, project, subproject/spacecraft, instrument suite, instrument, product

  auxiliary inputs/dependent attributes:
    local_repository, remote_repository, data_format,
  .. 
  """
  pass

# -----------------------------------------------
# End of dataset class
# -----------------------------------------------

class Project(dict):
  """
  Project class specifying dictionary of attributes specific for a given project.
  
  Important characteristics: 

  subprojects/spacecraft, instrument suites, instruments

  local_repository, remote_repository,
  data_format,

  .. 
  """

  def __init__(self, name):
    """
    """

# -----------------------------------------------
# End of Project class
# -----------------------------------------------
                                                                                                                                                       1018,1        Bot