;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_create_field_structure, fieldnames, fieldtypes, ndims
error=0
nfields=N_ELEMENTS(fieldnames)
ntypes =N_ELEMENTS(fieldtypes)

;FIELD_TYPES:
;  0 UNDEFINED Undefined
;  1 BYTE Byte
;  2 INT Integer
;  3 LONG Longword integer
;  4 FLOAT Floating point
;  5 DOUBLE Double-precision floating
;  6 COMPLEX Complex floating
;  7 STRING String
;  8 STRUCT Structure [NOT ALLOWED]
;  9 DCOMPLEX Double-precision complex
; 10 POINTER Pointer
; 11 OBJREF Object reference
; 12 UINT Unsigned Integer
; 13 ULONG Unsigned Longword Integer
; 14 LONG64 64-bit Integer
; 15 ULONG64 Unsigned 64-bit Integer 

type_names= $
   ['UNDEFINED', $
    'BYTE', $
    'INT', $
    'LONG', $
    'FLOAT', $
    'DOUBLE', $
    'COMPLEX', $
    'STRING', $
    'STRUCT', $
    'DCOMPLEX', $
    'POINTER', $
    'OBJREF', $
    'UINT', $
    'ULONG', $
    'LONG64', $
    'ULONG64']

IF nfields EQ 0 THEN BEGIN
    PRINT,' [STARDB_CREATE_FIELD_STRUCTURE] ERROR need field names and types'
    error=1
    RETURN,-1
ENDIF
IF ntypes EQ 0 THEN BEGIN
    PRINT,' [STARDB_CREATE_FIELD_STRUCTURE] setting field type to DOUBLE'
    fieldtypes='DOUBLE'
    ntypes=1
ENDIF
IF (ntypes EQ 1) AND (nfields GT 1) THEN BEGIN
    PRINT,' [STARDB_CREATE_FIELD_STRUCTURE] assuming same field type for all fields'
    fieldtypes=REPLICATE(fieldtypes,nfields)
    ntypes=nfields
ENDIF
IF (ntypes NE nfields) THEN BEGIN
    PRINT,' [STARDB_CREATE_FIELD_STRUCTURE] ERROR: #names and #types differ'
    error=1
    RETURN, -1
ENDIF
IF SIZE(fieldtypes,/TNAME) EQ 'STRING' THEN BEGIN
   xfieldtypes=BYTARR(ntypes)
   FOR itype=0,ntypes-1 DO BEGIN
      x=WHERE(type_names EQ STRUPCASE(fieldtypes[itype]))
      IF x[0] EQ -1 THEN BEGIN
         PRINT,' [STARDB_CREATE_FIELD_STRUCTURE] ERROR: field type unknown.'
         error=1
         RETURN,-1         
      ENDIF
      xfieldtypes[itype]=x
   ENDFOR
   fieldtypes=xfieldtypes
ENDIF
field=CREATE_STRUCT(fieldnames[0],(MAKE_ARRAY(1,TYPE=fieldtypes[0]))[0])
FOR ifield=1,nfields-1 DO BEGIN
    field=CREATE_STRUCT(field,fieldnames[ifield],(MAKE_ARRAY(1,TYPE=fieldtypes[ifield]))[0])
ENDFOR
IF N_ELEMENTS(ndims) EQ 0 THEN ndims=1
IF (N_ELEMENTS(ndims) NE 1) OR (ndims[0] NE 1) THEN BEGIN
    field=REPLICATE(field,ndims)
ENDIF 
RETURN,field
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_open, outfile, ERROR=error, UNIT=fileunit, COMPRESS=compress
COMMON stardb_write, unit, file
error=0
file=outfile

IF N_ELEMENTS(compress) EQ 0 THEN compress=0
IF (compress EQ 0) AND (STRLOWCASE(STRMID(STRTRIM(file,2),2,3,/REVERSE_OFFSET)) EQ '.gz') THEN compress=1
IF (compress EQ 1) AND (STRLOWCASE(STRMID(STRTRIM(file,2),2,3,/REVERSE_OFFSET)) NE '.gz') THEN file+='.gz'

OPENW,unit,file,/GET_LUN,COMPRESS=compress
fileunit=unit
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_close, ERROR=error
COMMON stardb_write, unit, file
error=0
CLOSE, unit
FREE_LUN, unit
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_head, version, ERROR=error
COMMON stardb_write, unit, file
error=0
IF N_ELEMENTS(version) EQ 0 THEN BEGIN 
    version=10000
    PRINT,'  [STARDB_WRITE_HEAD] setting version to '+STRING(version) 
ENDIF
IF (WHERE(SIZE(version,/TYPE) EQ [2,3,14]))[0] EQ -1 THEN BEGIN
    PRINT,' [STARDB_WRITE_INT] argument needs to be signed integer data type' 
    error=1
    RETURN
ENDIF
signature='ADBADBADBABDADBA'XULL
WRITEU, unit, ULONG64(signature)
WRITEU, unit, ULONG64(version)
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_tail, ERROR=error
COMMON stardb_write, unit, file
POINT_LUN, -unit, position
filesize=position+8ULL
WRITEU, unit, ULONG64(filesize)
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_str,array_in, ERROR=error
COMMON stardb_write, unit, file
error=0
IF SIZE(array_in,/TYPE) NE 7 THEN BEGIN
    PRINT,' [STARDB_WRITE_STR] argument needs to be string array' 
    error=1
    RETURN
ENDIF
len=STRLEN(array_in)
array=array_in
narray=N_ELEMENTS(array)
; align on 8 byte boundaries
FOR iarray=0,narray-1 DO BEGIN
    xlen=len[iarray] mod 8
    IF xlen NE 0 THEN BEGIN
        array[iarray]+=STRMID('       ',0,8-xlen)
    ENDIF
ENDFOR
WRITEU, unit, ULONG64(len), array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_uin,array_in, ERROR=error
COMMON stardb_write, unit, file
error=0
; IF (WHERE(SIZE(array_in,/TYPE) EQ [1,12,13,15]))[0] EQ -1 THEN BEGIN
;     PRINT,' [STARDB_WRITE_UIN] argument needs to be unsigned integer data type' 
;     error=1
;     RETURN
; ENDIF
WRITEU,unit, ULONG64(array_in)
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_int, array_in, ERROR=error
COMMON stardb_write, unit, file
error=0
; IF (WHERE(SIZE(array_in,/TYPE) EQ [2,3,14]))[0] EQ -1 THEN BEGIN
;     PRINT,' [STARDB_WRITE_INT] argument needs to be signed integer data type' 
;     error=1
;     RETURN
; ENDIF
WRITEU,unit, LONG64(array_in)
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_dbl, array_in, ERROR=error
COMMON stardb_write, unit, file
error=0
; IF (WHERE(SIZE(array_in,/TYPE) EQ [4,5]))[0] EQ -1 THEN BEGIN
;     PRINT,' [STARDB_WRITE_DBL] argument needs to be real float data type' 
;     error=1
;     RETURN
; ENDIF
ftag=' [STARDB_WRITE_DBL] '
narray=N_ELEMENTS(array_in)
status=FSTAT(unit)
nmax=2LL^26
IF (status.compress NE 0) AND (narray GT nmax) THEN BEGIN 
   PRINT,ftag+'using special algorithm to fix IDL bug for compressed files.'
   nsub=narray/nmax
   nrest=narray-nsub*nmax
   FOR i=0LL,nsub-1 DO BEGIN
      WRITEU,unit, DOUBLE(array_in[i*nmax:(i+1)*nmax-1])
   ENDFOR
   IF nrest GT 0 THEN BEGIN
      WRITEU,unit, DOUBLE(array_in[nsub*nmax:narray-1])
   ENDIF
ENDIF ELSE BEGIN
   WRITEU,unit, DOUBLE(array_in)
ENDELSE

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_stu, array_in, ERROR=error
COMMON stardb_write, unit, file
error=0
IF SIZE(array_in,/TYPE) NE 8 THEN BEGIN
    PRINT,' [STARDB_WRITE_STRUCTURE] argument needs to be of type structure' 
    error=1
    RETURN
ENDIF
ntags=N_TAGS(array_in)
; enhance to deal with strings?
FOR itag=0,ntags-1 DO BEGIN
    IF (WHERE(SIZE(array_in.(itag),/TYPE) EQ [5,13,14]))[0] EQ -1 THEN BEGIN
        PRINT,' [STARDB_WRITE_STRUCTURE] tag "'+(TAG_NAMES(array_in))[itag]+'" does not have 8 byte type'
        error=1
    ENDIF 
ENDFOR
IF error EQ 1 THEN BEGIN
    PRINT,' [STARDB_WRITE_STRUCTURE] ERROR(S): aboarding'
    RETURN
ENDIF
WRITEU,unit, array_in
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO stardb_write_file, $
                       outfile, $
                       NAME=name, $
                       COMMENTS=comments, $
                       FIELDNAMES=fieldnames, $
                       FIELDTYPES=fieldtypes, $
                       FIELDUNITS=fieldunits, $
                       FIELDFORMATS=fieldformats, $
                       FIELDFLAGS=fieldflags, $
                       ABUNDANCEE=abundancee, $
                       ABUNDANCEZ=abundancez, $
                       ABUNDANCEA=abundancea, $
                       FIELDDATA=fielddata, $
                       ABUNDANCEDATA=abundancedata, $
                       ABUNDANCE_TYPE=abundance_type, $
                       ABUNDANCE_CLASS=abundance_class, $
                       ABUNDANCE_UNIT=abundance_unit, $
                       ABUNDANCE_TOTAL=abundance_total, $
                       ABUNDANCE_NORM=abundance_norm, $
                       ABUNDANCE_DATA=abundance_data, $
                       ABUNDANCE_SUM=abundance_sum, $
                       COMPRESS=compress, $  
                       ERROR=error

starttime=SYSTIME(/SECONDS)

IF N_ELEMENTS(outfile) EQ 0 THEN BEGIN
   PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide output file name.'
   error=1
   RETURN
ENDIF

version=10100

; NEW CLASSIFICATION SCHEME
; A) abundance_type
; 0 isomere
; 1 isotope 
; 2 element
; 3 mass number (isobar)
IF N_ELEMENTS(abundance_type) EQ 0 THEN abundance_type=2

; B) abundance_class
; 0 all (raw)
; 1 rad (raw + decays)
; 2 dec (stable subset of radiso)
IF N_ELEMENTS(abundance_class) EQ 0 THEN abundance_class=2

; C) abundance_unit (should use 4)
; 0 mass (g)
; 1 mass (solar)
; 2 mol 
; 3 mass fraction (X)
; 4 mol fraction (YPS)
; 5 log epsion
; 6 [ ]
; 7 production factor (10^[])
; 8 log [Y/Si] + 6
IF N_ELEMENTS(abundance_unit) EQ 0 THEN abundance_unit=4

; D) abundance_total - for mol and mass fraction, [] and pf (should use 1)
; 0 initial (total not normalized, fallback is missing mass)
; 1 ejecta
IF N_ELEMENTS(abundance_total) EQ 0 THEN abundance_total=1

; E) abundance_norm(alization) - solar set - for [] and pf
; '' none
; 'AG89','Lod03' - which solar
IF N_ELEMENTS(abundance_norm) EQ 0 THEN abundance_norm='Lod03'

; F) abundance_data
; 0 : all ejecta (SN ejecta + wind) 
; 1 : SN ejecta, no wind 
; 2 : wind only 
; 3 : ejecta including fallback and wind
IF N_ELEMENTS(abundance_data) EQ 0 THEN abundance_data=0

; G) abundance_sum (for elements)
; 0 : mass fraction
; 1 : number fraction
IF N_ELEMENTS(abundance_sum) EQ 0 THEN abundance_sum=0

; QUESTIONS
; What would [Mg] be?
; A-2, B-2, C-6, D-1, E-'AG89', F-0
; What would log eps element be?
; A-2, B-2, C-5, D-1, E-'', F-0

ncomments=N_ELEMENTS(comments)
nstar=N_ELEMENTS(fielddata)
nfield=N_ELEMENTS(fieldnames)

CASE abundance_type OF
   0: BEGIN
      IF N_ELEMENTS(abundancee) EQ 0 THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide isomere energy levels'
         error=1
         RETURN
      ENDIF
      IF (N_ELEMENTS(abundancee) NE N_ELEMENTS(abundancez)) OR $
         (N_ELEMENTS(abundancee) NE N_ELEMENTS(abundancea)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of isomere A and Z data'
         error=1
         RETURN
      ENDIF  
   END
   1: BEGIN
      IF N_ELEMENTS(abundancez) EQ 0 THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide isotopte charge number'
         error=1
         RETURN
      ENDIF
      IF (N_ELEMENTS(abundancez) NE N_ELEMENTS(abundancea)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of isotope A and Z data'
         error=1
         RETURN
      ENDIF  
      IF N_ELEMENTS(abundancee) EQ 0 THEN abundancee=REPLICATE(0U,N_ELEMENTS(abundancez))
      IF N_ELEMENTS(abundancee) EQ 1 THEN abundancee=REPLICATE(abundancee[0],N_ELEMENTS(abundancez))
      IF (N_ELEMENTS(abundancez) NE N_ELEMENTS(abundancee)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of isotope E and Z data or skip E data'
         error=1
         RETURN
      ENDIF  
   END
   2: BEGIN
      IF N_ELEMENTS(abundancez) EQ 0 THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide element charge number'
         error=1
         RETURN
      ENDIF
      IF N_ELEMENTS(abundancea) EQ 0 THEN abundancea=REPLICATE(0U,N_ELEMENTS(abundancez))
      IF N_ELEMENTS(abundancea) EQ 1 THEN abundancea=REPLICATE(abundancea[0],N_ELEMENTS(abundancez))
      IF (N_ELEMENTS(abundancez) NE N_ELEMENTS(abundancea)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of A and Z data or skip A data'
         error=1
         RETURN
      ENDIF  
      IF N_ELEMENTS(abundancee) EQ 0 THEN abundancee=REPLICATE(0U,N_ELEMENTS(abundancez))
      IF N_ELEMENTS(abundancee) EQ 1 THEN abundancee=REPLICATE(abundancee[0],N_ELEMENTS(abundancez))
      IF (N_ELEMENTS(abundancez) NE N_ELEMENTS(abundancee)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of E and Z data or skip E data'
         error=1
         RETURN
      ENDIF  
   END
   3: BEGIN
      IF N_ELEMENTS(abundancea) EQ 0 THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide isobar neutron number'
         error=1
         RETURN
      ENDIF
      IF N_ELEMENTS(abundancez) EQ 0 THEN abundancez=REPLICATE(0U,N_ELEMENTS(abundancea))
      IF N_ELEMENTS(abundancez) EQ 1 THEN abundancez=REPLICATE(abundancez[0],N_ELEMENTS(abundancea))
      IF (N_ELEMENTS(abundancea) NE N_ELEMENTS(abundancez)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of A and Z data or skip Z data'
         error=1
         RETURN
      ENDIF  
      IF N_ELEMENTS(abundancee) EQ 0 THEN abundancee=REPLICATE(0U,N_ELEMENTS(abundancea))
      IF N_ELEMENTS(abundancee) EQ 1 THEN abundancee=REPLICATE(abundancee[0],N_ELEMENTS(abundancea))
      IF (N_ELEMENTS(abundancea) NE N_ELEMENTS(abundancee)) THEN BEGIN
         PRINT,' [STARDB_WRITE_FILE] ERROR: need to provide compatible # of E and A data or skip E data'
         error=1
         RETURN
      ENDIF  
   END
ENDCASE

nabu=N_ELEMENTS(abundancez)

IF N_ELEMENTS(name) NE 1 THEN BEGIN 
    PRINT,' [STARDB_WRITE_FILE] ERROR: name has to be scalar value'
    error=1
    RETURN
ENDIF

IF SIZE(name,/TYPE) NE 7 THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: name has to be of type "STRING"'
    error=1
    RETURN
ENDIF

IF nfield NE N_ELEMENTS(fieldtypes) THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: incompatible dimensions - types'
    error=1
    RETURN
ENDIF

IF nfield NE N_ELEMENTS(fieldunits) THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: incompatible dimensions - units'
    error=1
    RETURN
ENDIF

IF nfield NE N_ELEMENTS(fieldformats) THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: incompatible dimensions - formats'
    error=1
    RETURN
ENDIF

IF nfield NE N_ELEMENTS(fieldflags) THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: incompatible dimensions - flags'
    error=1
    RETURN
ENDIF

IF nstar*nabu NE N_ELEMENTS(abundancedata) THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: incompatible dimensions - data'
    error=1
    RETURN
ENDIF

fieldformtregexp='^[A-Z][0-9]+(\.[0-9]+)?$'
IF (WHERE(stregex(fieldformats,fieldformtregexp,/FOLD_CASE) NE 0))[0] NE -1 THEN BEGIN
    PRINT,' [STARDB_WRITE_FILE] ERROR: invalid field format.'
    PRINT,' [STARDB_WRITE_FILE] use "'+fieldformtregexp+'"'
    error=1
    RETURN
ENDIF

stardb_write_open, outfile, COMPRESS=compress
stardb_write_head, version

stardb_write_str, name
stardb_write_uin, ncomments
stardb_write_str, comments

stardb_write_uin, nstar
stardb_write_uin, nfield
stardb_write_uin, nabu

stardb_write_uin, abundance_type
stardb_write_uin, abundance_class
stardb_write_uin, abundance_unit
stardb_write_uin, abundance_total
stardb_write_str, abundance_norm 
stardb_write_uin, abundance_data
stardb_write_uin, abundance_sum

stardb_write_str, fieldnames
stardb_write_str, fieldunits
stardb_write_uin, fieldtypes
stardb_write_str, fieldformats
stardb_write_uin, fieldflags

stardb_write_uin, abundancez
stardb_write_uin, abundancea
stardb_write_uin, abundancee

stardb_write_stu, fielddata
stardb_write_dbl, abundancedata

stardb_write_tail
stardb_write_close

finishtime=SYSTIME(/SECONDS)
PRINT,' [STARDB_WRITE_FILE] data saved in '+time2human(finishtime-starttime)

END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_str,ndims, ERROR=error
COMMON stardb_read, unit
error=0
IF N_ELEMENTS(ndims) EQ 0 THEN ndims=1
narray=PRODUCT(ULONG64(ndims),/INTEGER)
len=ULON64ARR(narray)
READU, unit, len
array=STRARR(narray)
FOR iarray=0,narray-1 DO BEGIN
    ilen=len[iarray]
    xlen=ilen mod 8
    IF xlen NE 0 THEN ilen+=8-xlen
    array[iarray]=STRING('',FORMAT="(A"+STRTRIM(STRING(ilen),2)+")")
ENDFOR
READU,unit, array
FOR iarray=0,narray-1 DO BEGIN
    array[iarray]=STRMID(array[iarray],0,len[iarray])
ENDFOR
array=REFORM(array,ndims,/overwrite)
IF N_ELEMENTS(array) EQ 1 THEN array=array[0]
RETURN, array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_uin,ndims, ERROR=error
COMMON stardb_read, unit
error=0
IF N_ELEMENTS(ndims) EQ 0 THEN ndims=1
narray=PRODUCT(ULONG64(ndims),/INTEGER)
array=ULON64ARR(narray)
READU, unit, array
array=REFORM(array,ndims,/overwrite)
IF N_ELEMENTS(array) EQ 1 THEN array=array[0]
RETURN, array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_int,ndims, ERROR=error
COMMON stardb_read, unit
error=0
IF N_ELEMENTS(ndims) EQ 0 THEN ndims=1
narray=PRODUCT(ULONG64(ndims),/INTEGER)
array=LON64ARR(narray)
READU, unit, array
array=REFORM(array,ndims,/overwrite)
IF N_ELEMENTS(array) EQ 1 THEN array=array[0]
RETURN, array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_dbl,ndims, ERROR=error
COMMON stardb_read, unit
error=0
ftag=' [STARDB_READ_DBL] '
narray=PRODUCT(ULONG64(ndims),/INTEGER)
array=DBLARR(narray)
status=FSTAT(unit)
nmax=2LL^26
IF (status.compress NE 0) AND (narray GT nmax) THEN BEGIN 
   PRINT,ftag+'using special algorithm to fix IDL bug for compressed files.'
   nsub=narray/nmax
   nrest=narray-nsub*nmax
   FOR i=0LL,nsub-1 DO BEGIN
      subarray=DBLARR(nmax)
      READU, unit, subarray
      array[i*nmax:(i+1)*nmax-1]=subarray
   ENDFOR
   IF nrest GT 0 THEN BEGIN
      subarray=DBLARR(nrest)
      READU, unit, subarray
      array[nsub*nmax:narray-1]=subarray   
   ENDIF
   clear,subarray
ENDIF ELSE BEGIN
   READU, unit, array
ENDELSE
array=REFORM(array,ndims,/overwrite)
IF N_ELEMENTS(array) EQ 1 THEN array=array[0]
RETURN, array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_stu,ndims,fieldnames,fieldtypes, ERROR=error
COMMON stardb_read, unit
error=0
IF N_ELEMENTS(ndims) EQ 0 THEN ndims=1
narray=PRODUCT(ULONG64(ndims),/INTEGER)
nfields=N_ELEMENTS(fieldnames)
IF N_ELEMENTS(fieldtypes) NE nfields THEN BEGIN
    PRINT,' [STARDB_READ_STRING] need same number of elements in names and types'
    error=1
    RETURN,-1
ENDIF
FOR ifield=0,nfields-1 DO BEGIN
    IF (WHERE(fieldtypes[ifield] EQ [5,13,14]))[0] EQ -1 THEN BEGIN
        PRINT,' [STARDB_READ_STRUCTURE] ERROR: data type not yet supported'
        PRINT,' [STARDB_READ_STRUCTURE] ERROR: only supporting 8-byte scalar data types'
        error=1
        RETURN,-1
    ENDIF
ENDFOR
field=stardb_create_field_structure(fieldnames, fieldtypes)
array=REPLICATE(field,narray)
READU, unit, array
array=REFORM(array,ndims,/overwrite)
IF N_ELEMENTS(array) EQ 1 THEN array=array[0]
RETURN, array
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

FUNCTION stardb_read_tail, ERROR=error
COMMON stardb_read, unit
error=0
POINT_LUN, -unit, position
filesize=position+8ULL
savedsize=0ULL
READU, unit, savedsize
IF filesize EQ savedsize THEN BEGIN
    PRINT,' [STARDB_READ_TAIL] file integrity seems OK'
ENDIF ELSE BEGIN
    PRINT,' [STARDB_READ_TAIL] ERROR: file integrity seems broken'
    error=1
ENDELSE
RETURN,error
END

;-----------------------------------------------------------------------
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;-----------------------------------------------------------------------

PRO loadstardb, file, ERROR=error, VERBOSE=verbose, RELOAD=reload0
COMMON stardb_read, unit

@stardb_common

ftag=' [LOADSTARDB] '

starttime=SYSTIME(/SECONDS)

IF N_ELEMENTS(verbose) EQ 0 THEN verbose=2

error=0

; idea about the data format

IF N_ELEMENTS(stardb_dirtydata) EQ 0 THEN stardb_dirtydata=1

type=SIZE(file,/TYPE)
IF (type EQ 0) AND $
   (stardb_dirtydata EQ 0) AND $
   (N_ELEMENTS(stardb_filename) NE 0) THEN BEGIN
   file=stardb_filename
ENDIF

IF N_ELEMENTS(STARDB_DATADIR) NE 0 THEN BEGIN
   data_dir=STARDB_DATADIR
ENDIF ELSE BEGIN
   data_dir=GETENV('STARDB_DATADIR')
ENDELSE

IF WHERE(type EQ [1,2,3,4,5,6,12,13,14,15]) NE -1 THEN BEGIN
    file=LONG(file)
    IF file LE 0 THEN BEGIN
        PRINT,ftag+'ERROR: file number needs to be integer and bigger than zero.'
        error=1
        RETURN
    ENDIF
    IF N_ELEMENTS(STARDB_FILELIST) NE 0 THEN BEGIN
       listfile=STARDB_FILELIST
    ENDIF ELSE BEGIN
       listfile=GETENV('STARDB_FILELIST')
    ENDELSE    
    IF listfile EQ '' THEN filelist='stardb_filelist.txt'
    search_dir=''
    IF data_dir NE '' THEN search_dir=[data_dir,search_dir]
    files=data_dir+listfile
    x=WHERE(files NE '')
    IF x[0] NE -1 THEN BEGIN
        files=files[x]
        listfile=FILE_SEARCH(files,/EXPAND_ENVIRONMENT, COUNT=count)
    ENDIF ELSE BEGIN 
        count=0
    ENDELSE
    IF count EQ 0 THEN BEGIN
        PRINT,ftag+'ERROR: no file list file found.'
        error=1
        RETURN
    ENDIF
    listfile=listfile[0]
    PRINT,ftag+'using list file '+listfile
    OPENR,listfileunit,listfile,/GET_LUN
    filelist=['']
    WHILE NOT EOF(listfileunit) DO BEGIN
        s=''
        READF,listfileunit,s,FORMAT="(A)"
        filelist=[filelist,s]
    ENDWHILE
    CLOSE,listfileunit
    FREE_LUN,listfileunit
    nfiles=N_ELEMENTS(filelist)
    IF nfiles LT file THEN BEGIN
        PRINT,ftag+'ERROR: file number bigger than number of entries in list.'
        error=1
        RETURN
    ENDIF
    loadfile=filelist[file]
ENDIF
IF type EQ 7 THEN BEGIN
    loadfile=file
ENDIF
IF N_ELEMENTS(loadfile) EQ 0 THEN BEGIN
    filename='(UNDEFINED)'
    IF N_ELEMENTS(file) NE 0 THEN filemame=file
    PRINT,ftag+'ERROR: no such file ',filename
    error=1
    RETURN
ENDIF

search_dir=''
IF data_dir NE '' THEN search_dir=[data_dir,search_dir]
files=search_dir+loadfile
x=WHERE(files NE '')
IF x[0] NE -1 THEN BEGIN
    files=files[x]
    loadfiles=FILE_SEARCH(files,/EXPAND_ENVIRONMENT, COUNT=count)
ENDIF ELSE BEGIN 
    count=0
ENDELSE

IF (count EQ 0) AND (STRLOWCASE(STRMID(loadfile,2,3,/reverse)) NE '.gz') THEN BEGIN
    loadfiles=FILE_SEARCH(files+'.gz',/EXPAND_ENVIRONMENT, COUNT=count)
ENDIF
IF (count EQ 0) OR (loadfiles[0] EQ '') THEN BEGIN
    PRINT,ftag+'COULD NOT FIND FILE "'+loadfile+'"'
    loadfiles=DIALOG_PICKFILE(TITLE='Please select a stardb file', $
                              /NOCONFIRM, $
                              /READ, $
                              /MUST_EXIST, $
                              PATH=data_dir, $
                              FILTER='*')
    count=N_ELEMENTS(loadfiles)
    IF loadfiles[0] EQ '' THEN count=0
ENDIF
IF count EQ 0 THEN BEGIN
    PRINT,ftag+'ERROR: file "'+loadfile+'" not found'
    error=1
    RETURN
ENDIF
loadfile=loadfiles[0]

PRINT,ftag+'using file   '+loadfile

IF N_ELEMENTS(reload0) EQ 0 THEN reload0=0
reload=reload0
IF N_ELEMENTS(stardb_filename) EQ 0 THEN stardb_filename=''
IF stardb_filename NE loadfile THEN BEGIN 
    reload=1
ENDIF
IF N_ELEMENTS(stardb_dirtydata) EQ 0 THEN stardb_dirtydata=1
IF stardb_dirtydata EQ 1 THEN reload=1
IF reload EQ 0 THEN BEGIN 
   PRINT,ftag+'using previously loaded data'    
   RETURN
ENDIF
stardb_dirtydata=1

compress=STRLOWCASE(STRMID(STRTRIM(loadfile,2),2,3,/REVERSE_OFFSET)) EQ '.gz'

OPENR,unit,loadfile,/GET_LUN, $
  COMPRESS=compress

signature='ADBADBADBABDADBA'XULL
filesignature=0ULL
READU,unit,filesignature
IF signature NE filesignature THEN BEGIN
    CLOSE,unit
    IF verbose GE 2 THEN BEGIN
        PRINT,ftag+'swapping endian'
    ENDIF
    OPENR,unit,loadfile,/SWAP_ENDIAN,COMPRESS=compress
    READU,unit,filesignature
ENDIF

IF signature NE filesignature THEN BEGIN
    PRINT,ftag+'ERROR: not a stardb file: "'+loadfile+'"'
    error=1
    GOTO, closefiles
ENDIF

filestat=fstat(unit)

IF verbose GE 2 THEN BEGIN
    PRINT,ftag+'loading file '+loadfile
    IF compress THEN sc=' (compressed)' ELSE sc=''
    PRINT,ftag+'filesize '+byte2human(filestat.size)+sc
ENDIF

version=stardb_read_uin()

PRINT,ftag+'data version:  '+STRTRIM(STRING(version,FORMAT="(I8)"),2)

name=stardb_read_str()
PRINT,ftag+'data set name: '+name


ncomment=stardb_read_uin()
comments=stardb_read_str(ncomments)

PRINT,ftag+'=========================================================='
FOR icomment=0,ncomment-1 DO BEGIN
    PRINT,ftag+'COMMENT: '+STRTRIM(comments[icomment],1)
ENDFOR
PRINT,ftag+'=========================================================='

nstar=stardb_read_uin()
nfield=stardb_read_uin()
nabu=stardb_read_uin()

IF version LT 10100 THEN BEGIN
   iabutype=stardb_read_uin()
   IF iabutype NE 1 THEN BEGIN
      PRINT,ftag+'ERROR: currently only supporting element data (type 1)'
      error=1
      GOTO, closefiles
   ENDIF
   abundance_type=2
   abundance_class=2
   abundance_unit=7
   abundance_total=0
   abundance_norm='Lod03'
   abundance_data =0
   abundance_sum =1
ENDIF ELSE BEGIN
   abundance_type =stardb_read_uin()
   abundance_class=stardb_read_uin()
   abundance_unit =stardb_read_uin()
   abundance_total=stardb_read_uin()
   abundance_norm =stardb_read_str()
   abundance_data =stardb_read_uin()
   abundance_sum  =stardb_read_uin()
ENDELSE

PRINT,ftag+'data sets:      '+STRING(nstar, FORMAT='(I6)')
PRINT,ftag+'abundance sets:   '+STRING(nabu, FORMAT='(I4)')

PRINT,ftag+'----------------------------------------------------------'

abundance_type_names=['isomere','isotope','element','mass number (isobar)']
PRINT,ftag+'abundance type:  '+STRING(abundance_type, FORMAT='(I1)')+' - '+abundance_type_names[abundance_type]

abundance_class_names=['all (raw)','rad (raw + decays)','dec (stable subset of radiso)']
PRINT,ftag+'abundance class: '+STRING(abundance_class, FORMAT='(I1)')+' - '+abundance_class_names[abundance_class]

abundance_unit_names=['mass (g)','mass (solar)','mol','mass fraction (X)','mol fraction (YPS)','log epsion','[ ]','production factor (10^[])','log [Y/Si] + 6']
PRINT,ftag+'abundance unit:  '+STRING(abundance_unit, FORMAT='(I1)')+' - '+abundance_unit_names[abundance_unit]

abundance_total_names=['initial (total, not normalized, fallback is missing mass)','ejecta']
PRINT,ftag+'abundance total: '+STRING(abundance_total, FORMAT='(I1)')+' - '+abundance_total_names[abundance_total]

s=abundance_norm
IF s EQ '' THEN s='(NONE)'
PRINT,ftag+'abundance norm:      '+s

abundance_data_names=['all ejecta (SN ejecta + wind)','SN ejecta, no wind','wind only','ejecta including fallback and wind']
PRINT,ftag+'abundance data:  '+STRING(abundance_data, FORMAT='(I1)')+' - '+abundance_data_names[abundance_data]

abundance_sum_names=['mass fraction','number fraction']
PRINT,ftag+'abundance sum:   '+STRING(abundance_sum, FORMAT='(I1)')+' - '+abundance_sum_names[abundance_sum]


fieldnames=stardb_read_str(nfield)
fieldunits=stardb_read_str(nfield)
fieldtypes=stardb_read_uin(nfield)
fieldformats=stardb_read_str(nfield)
fieldflags=stardb_read_uin(nfield)

fieldnames=STRTRIM(fieldnames,2)
fieldunits=STRTRIM(fieldunits,2)
fieldformats=STRTRIM(fieldformats,2)

typenames=[ $
            'UNDEFINED', $
            'BYTE', $
            'INT', $
            'LONG', $
            'FLOAT', $
            'DOUBLE', $
            'COMPLEX', $
            'STRING', $
            'STRUCT', $
            'DCOMPLEX', $
            'POINTER', $
            'OBJREF', $
            'UINT', $
            'ULONG', $
            'LONG64', $
            'ULONG64']

flagnames=['parameter','property']

PRINT,ftag+'----------------------------------------------------------'
PRINT,ftag+''+STRTRIM(STRING(nfield),2)+' data fields: '

sx=STRING('',FORMAT="(A128)")
l1=MAX(STRLEN(fieldnames))
l2=MAX(STRLEN(fieldunits))
l3=MAX(STRLEN(typenames))
format="('"+ftag+"',A"+STRTRIM(STRING(l1),2)+",' ',A"+STRTRIM(STRING(l2+2),2)+",' ',A"+STRTRIM(STRING(l3+2),2)+",' ',A)"
flagnames=['parameter','property']
FOR ifield=0,nfield-1 DO BEGIN
    PRINT,fieldnames[ifield]+STRMID(sx,0,l1-STRLEN(fieldnames[ifield])), $
      '['+fieldunits[ifield]+'] '+STRMID(sx,0,l2-STRLEN(fieldunits[ifield])), $
      '('+typenames[fieldtypes[ifield]]+')'+STRMID(sx,0,l3-STRLEN(fieldtypes[ifield])), $
      '<'+flagnames[fieldflags[ifield]]+'>', $
      FORMAT=format
    IF (WHERE(typenames[fieldtypes[ifield]] EQ  $
              ['DOUBLE','LONG64','ULONG64']))[0] EQ -1 THEN BEGIN
        PRINT,ftag+'ERROR: data type not yet supported'
        PRINT,ftag+'ERROR: only supporting 8-byte scalar data types'
        error=1
        GOTO, closefiles        
    ENDIF
ENDFOR

abu_Z=stardb_read_uin(nabu)
IF version LT 10100 THEN BEGIN
   abu_A=REPLICATE(0U,nabu)
   abu_E=REPLICATE(0U,nabu)
ENDIF ELSE BEGIN
   abu_A=stardb_read_uin(nabu)
   abu_E=stardb_read_uin(nabu)
ENDELSE
abu_index=zae2idx(abu_Z,abu_A,abu_E)

field_data=stardb_read_stu(nstar,fieldnames, fieldtypes)

l1=MAX(STRLEN(fieldnames))
l2=0
l3=0
nvalues=LON64ARR(nfield)
values=DBLARR(nfield,nstar)
ivalues=LON64ARR(nfield,nstar)
FOR ifield=0,nfield-1 DO BEGIN
; values
    v=field_data.(ifield)
    vs=SORT(v)
    vu=UNIQQ(v,vs,INDICES=vx)
    vv=v[vu]
    nv=N_ELEMENTS(vv)
    values[ifield,0:nv-1]=vv
    nvalues[ifield]=nv
    vx=[-1,vx]
    FOR iv=0,nv-1 DO BEGIN
       ivalues[ifield,vs[vx[iv]+1:vx[iv+1]]]=iv
    ENDFOR
; output formatting
    ls=(STREGEX(fieldformats[ifield],'^[A-Z]([0-9]+)',/FOLD_CASE,/EXTRACT,/SUBEXPR))[1]
    flen=LONG(ls)
    l2>=flen
    l3>=STRLEN(STRTRIM(STRING(nv),2))
ENDFOR
nvalues_max=MAX(nvalues)
values= values[*,0:nvalues_max-1]

xpar=WHERE(fieldflags EQ 0)
npar=N_ELEMENTS(xpar)
IF xpar[0] EQ -1 THEN npar=0
IF npar GT 0 THEN BEGIN 
   PRINT,ftag+'----------------------------------------------------------'
   PRINT,ftag+'PARAMETER RANGES:'
   FOR ifield=0,nfield-1 DO BEGIN
      IF fieldflags[ifield] EQ 0 THEN BEGIN
         fmax=MAX(values[ifield,0:nvalues[ifield]-1],MIN=fmin)
         line=[fieldnames[ifield]+': '+STRMID(sx,0,l1-STRLEN(fieldnames[ifield])), $
               STRING(fmin,FORMAT='('+fieldformats[ifield]+')'), $
               STRING(fmax,FORMAT='('+fieldformats[ifield]+')'), $
               STRTRIM(STRING(nvalues[ifield]),2)]
         format="('"+ftag+"',"+ $
                "A"+STRTRIM(STRING(l1+2),2)+","+ $
                "A"+STRTRIM(STRING(l2),2)+","+ $
                "' ... ',"+ $
                "A"+STRTRIM(STRING(l2),2)+","+ $
                "' (',A"+STRTRIM(STRING(l3),2)+",' values)')"
         PRINT,line,FORMAT=format
      ENDIF
   ENDFOR
ENDIF
xprop=WHERE(fieldflags NE 0)
nprop=N_ELEMENTS(xprop)
IF xprop[0] EQ -1 THEN nprop=0
IF nprop GT 0 THEN BEGIN 
   PRINT,ftag+'----------------------------------------------------------'
   PRINT,ftag+'PROPERTY RANGES:'
   FOR ifield=0,nfield-1 DO BEGIN
      IF fieldflags[ifield] EQ 1 THEN BEGIN
         fmax=MAX(values[ifield,0:nvalues[ifield]-1],MIN=fmin)
         line=[fieldnames[ifield]+': '+STRMID(sx,0,l1-STRLEN(fieldnames[ifield])), $
               STRING(fmin,FORMAT='('+fieldformats[ifield]+')'), $
               STRING(fmax,FORMAT='('+fieldformats[ifield]+')'), $
               STRTRIM(STRING(nvalues[ifield]),2)]
         format="('"+ftag+"',"+ $
                "A"+STRTRIM(STRING(l1+2),2)+","+ $
                "A"+STRTRIM(STRING(l2),2)+","+ $
                "' ... ',"+ $
                "A"+STRTRIM(STRING(l2),2)+","+ $
                "' (',A"+STRTRIM(STRING(l3),2)+",' values)')"
         PRINT,line,FORMAT=format
      ENDIF
   ENDFOR
ENDIF
IF npar GT 0 THEN BEGIN 
   PRINT,ftag+'----------------------------------------------------------'
   PRINT,ftag+'PARAMETER VALUES:'
   FOR ifield=0,nfield-1 DO BEGIN
      IF fieldflags[ifield] EQ 0 THEN BEGIN
         PRINT, $
            ftag+''+fieldnames[ifield]+': '
         ls=(STREGEX(fieldformats[ifield],'^[A-Z]([0-9]+)',/FOLD_CASE,/EXTRACT,/SUBEXPR))[1]
         flen=LONG(ls)
         format=''
         FOR il=0,(50/flen)-1 DO BEGIN
            format+='," ",'+fieldformats[ifield]
         ENDFOR
         PRINT, $
            values[ifield,0:nvalues[ifield]-1], $
            FORMAT='(" [LOADSTARDB] "'+ $
            format+ $
            ')'
      ENDIF
   ENDFOR
ENDIF
maxpropvalues=100
IF nprop GT 0 THEN BEGIN 
   PRINT,ftag+'----------------------------------------------------------'
   PRINT,ftag+'PROPERTY VALUES:'
   FOR ifield=0,nfield-1 DO BEGIN
      IF fieldflags[ifield] EQ 1 THEN BEGIN
         PRINT, $
            ftag+fieldnames[ifield]+': '
         IF nvalues[ifield] GT maxpropvalues THEN BEGIN
            PRINT, ftag+'(more than '+STRTRIM(STRING(maxpropvalues),2)+' values)'
         ENDIF ELSE BEGIN
            ls=(STREGEX(fieldformats[ifield],'^[A-Z]([0-9]+)',/FOLD_CASE,/EXTRACT,/SUBEXPR))[1]
            flen=LONG(ls)
            format=''
            FOR il=0,(50/flen)-1 DO BEGIN
               format+='," ",'+fieldformats[ifield]
            ENDFOR
            PRINT, $
               values[ifield,0:nvalues[ifield]-1], $
               FORMAT='(" [LOADSTARDB] "'+ $
               format+ $
               ')'
         ENDELSE
      ENDIF
   ENDFOR
ENDIF
PRINT,ftag+'----------------------------------------------------------'
   
abu_data=stardb_read_dbl([nabu,nstar])

stardb_dirtydata=stardb_read_tail()

stardb_nvalues=nvalues
stardb_values = values
stardb_indices=ivalues 

closefiles:

CLOSE,unit
FREE_LUN,unit

abu_idx=ionidx(abu_a,abu_z)

x=WHERE(TOTAL(abu_data,1) EQ 0)
IF x[0] NE -1 THEN BEGIN
   PRINT,ftag+'ERROR: found zero data sets'
   ERROR=1
ENDIF

IF stardb_dirtydata EQ 0 THEN BEGIN 
   stardb_filename=loadfile
   PRINT,ftag+'data loaded in '+STRING(SYSTIME(/SECONDS)-starttime)+' seconds.'
ENDIF ELSE BEGIN
   clear,loadfile
   PRINT,ftag+'ERROR: Data seems bad.'
   ERORR=1
ENDELSE

END