text:
# define a primitive type for IBM floating number
primitive type IBMFloat32 <: AbstractFloat 32 end

ieeeOfPieces(fr::UInt32, exp::Int32, sgn::UInt32) = reinterpret(Float32, convert(UInt32,fr >>> 9) | convert(UInt32,exp << 23) | sgn) :: Float32

"""
   convert IBM Float32 to IEEE Float32
"""
function convert(::Type{Float32}, ibm::IBMFloat32)

    local fr::UInt32 = ntoh(reinterpret(UInt32, ibm))
    local sgn::UInt32 = fr & 0x80000000 # save sign
    fr <<= 1 # shift sign out
    local exp::Int32 = convert(Int32,fr >>> 25) # save exponent
    fr <<= 7 # shift exponent out

    if (fr == convert(UInt32,0))
       zero(Float32)
    else
       # normalize the signficand
       local norm::UInt32 = leading_zeros(fr)
       fr <<= norm
       exp = (exp << 2) - 130 - norm

       # exp <= 0 --> ieee(0,0,sgn)
       # exp >= 255 --> ieee(0,255,sgn)
       # else -> ieee(fr<<1, exp, sgn)
       local clexp::Int32 = exp & convert(Int32,0xFF)
       ieeeOfPieces(clexp == exp ? fr << 1 : convert(UInt32,0), clexp, sgn)
    end
end

"""
   the mapping between edbcb and ascii code, as the text header of segy file is
encode with edbcb format, so we need first translate it into ascii format,
then convert to character.
"""
const edbcb2ascii_table = [0, 1, 2, 3, 156, 9, 134,127,151,141,142, 11, 12, 13, 14, 15,
       16, 17, 18, 19,157,133,  8,135, 24, 25,146,143, 28, 29, 30, 31,
       128,129,130,131,132, 10, 23, 27,136,137,138,139,140,  5,  6,  7,
       144,145, 22,147,148,149,150,  4,152,153,154,155, 20, 21,158, 26,
       32,160,161,162,163,164,165,166,167,168, 91, 46, 60, 40, 43, 33,
       38,169,170,171,172,173,174,175,176,177, 93, 36, 42, 41, 59, 94,
       45, 47,178,179,180,181,182,183,184,185,124, 44, 37, 95, 62, 63,
       186,187,188,189,190,191,192,193,194, 96, 58, 35, 64, 39, 61, 34,
       195, 97, 98, 99,100,101,102,103,104,105,196,197,198,199,200,201,
       202,106,107,108,109,110,111,112,113,114,203,204,205,206,207,208,
       209,126,115,116,117,118,119,120,121,122,210,211,212,213,214,215,
       216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,
       123, 65, 66, 67, 68, 69, 70, 71, 72, 73,232,233,234,235,236,237,
       125, 74, 75, 76, 77, 78, 79, 80, 81, 82,238,239,240,241,242,243,
       92,159, 83, 84, 85, 86, 87, 88, 89, 90,244,245,246,247,248,249,
       48, 49, 50, 51, 52, 53, 54, 55, 56, 57,250,251,252,253,254,255]

"""
   create the mapping table which maps ascii code to edbcb code.
"""
function create_ascii2edbcb_table(e2a::Vector{Int64})

    a2e = zeros(Int64, 256)
    for i = 1 : length(e2a)
       	a2e[e2a[i]+1] = i-1
    end
    return a2e
end
const ascii2edbcb_table = create_ascii2edbcb_table(edbcb2ascii_table);

"""
   conversions between ascii and edbcb table
"""
function conversion_between_edbcb_ascii(a::Vector{UInt8}; flag="edbcb2ascii")

    r = zeros(UInt, length(a))
    if flag == "edbcb2ascii"

       idx = 1
       for c in a
           c1 = edbcb2ascii_table[c+1]
           r[idx] = c1
           idx = idx + 1
       end

    elseif flag == "ascii2edbcb"
       idx = 1
       for c in a
           c1 = ascii2edbcb_table[c+1]
           r[idx] = c1
           idx = idx + 1
       end

    end

    return convert(Vector{UInt8}, r)
end

"""
   read text header
"""
function read_text_header(path::String)

    # read the text header as UInt8
	  fid = open(path, "r")
		thdr = read(fid, 3200)
		close(fid)

    # convert edbcb coding to ascii coding
    thdr = conversion_between_edbcb_ascii(thdr; flag="edbcb2ascii")

    # convert each element to character
		thdr = convert(Vector{Char}, thdr)

		return thdr
end

"""
   print text header
"""
function print_text_header(thdr)

    for i = 1 : 40
        il = (i-1)*80+1
        iu = il+80-1
        println(join(thdr[il:iu]))
    end
    return nothing
end

"""
   put the text header to segy file
"""
function put_text_header(path::String, thdr::Vector{Char})

    #convert string of Char to a vector of UInt
    a = convert(Vector{UInt8}, thdr)

    # convert ascii coding to edbcb coding
    a = conversion_between_edbcb_ascii(a; flag="ascii2edbcb")

    # write to file
    fid = open(path, "w")
    write(fid, a)
    close(fid)

		return nothing
end

"""
   composite type of file header
*   ntrpe : number of traces per ensemble
*   dt    : sample interval in microseconds
*   dtfr  : sample interval in microseconds in original field recording
*   ns    : number of samples per data trace
*   nsfr  : number of samples per data tracf in original field recording
*   fmtc  : data sample format code(1=>IBM; 5=>IEEE)
*   expf  : expected number of trace per trace ensemble (such as CMP fold)
*   trsc  : trace sorting code(1=>as recorded, 2=>CDP, 3=>single fold profile, 4=>post stack
           5=>CSP, 6=>CRP, 7=>common offset, 8=>CMP, 9=>common conversion point)
*   sfs   : sweep frequency at start (Hz)
*   sfe   : sweep frequency at end (Hz)
*   slen  : sweep length (ms)
*   styp  : sweep type code (1=>linear, 2=>parabolic, 3=>exponential)
*   tnumsc: trace number sweep channel
*   stalens: sweep trace taper length in milliseconds at start
*   stalene: sweep trace taper length in milliseconds at end
*   taty   : taper type(1=>linear, 2=>cos^2, 3=>other)
*   corr   : correlated trace(1=>no, 2=>yes)
*   rgc    : gain recovered(1=>yes, 2=>no)
*   arm    : amplitude recover method(1=>none, 2=>spherical divergence, 3=>AGC, 4=>other)
*   unit   : measurement system(1=>meters, 2=>Feet)
*   pol    : impulse polarity(1=>increase pressure and upward movement is negative number,
                              2=>increase pressure and upward movement is positive number)
*   vpol   : Vibratory polarity code
*   fvn    : file reversion number
*   fltf   : fixed trace length flag
*   netfh  : number of extended text file header(0=> no extende, -1=>variable, positive=>exactly that many)
"""
mutable struct FileHeader
# corresponding to 3200-3259
	  jobid  :: Int32
	  linnum :: Int32
	  renum  :: Int32
	  ntrpe  :: Int16
	  natrpe :: Int16
	  dt     :: Int16
	  dtfr   :: Int16
	  ns     :: Int16
	  nsfr   :: Int16
	  fmtc   :: Int16
	  expf   :: Int16
	  trsc   :: Int16
	  vsumc  :: Int16
	  sfs    :: Int16
	  sfe    :: Int16
	  slen   :: Int16
	  styp   :: Int16
	  tnumsc :: Int16
	  stalens:: Int16
	  stalene:: Int16
	  taty   :: Int16
	  corr   :: Int16
	  rgc    :: Int16
	  arm    :: Int16
	  unit   :: Int16
	  pol    :: Int16
	  vpol   :: Int16
# 3260-3501 is unsigned, need to skip 242 byte
	  fltf   :: Int16 # 3502-3503
	  netfh  :: Int16 # 3504-3505
# then skip 94 byte to the start point of trace header
end

# record the start position of file header item
const file_header_location = Dict(
      :jobid   => 3200,
      :linnum  => 3204,
      :renum   => 3208,
      :ntrpe   => 3212,
      :natrpe  => 3214,
      :dt      => 3216,
      :dtfr    => 3218,
      :ns      => 3220,
      :nsfr    => 3222,
      :fmtc    => 3224,
      :expf    => 3226,
      :trsc    => 3228,
      :vsumc   => 3230,
      :sfs     => 3232,
      :sfe     => 3234,
      :slen    => 3236,
      :styp    => 3238,
      :tnumsc  => 3240,
      :stalens => 3242,
      :stalene => 3244,
      :taty    => 3246,
      :corr    => 3248,
      :rgc     => 3250,
      :arm     => 3252,
      :unit    => 3254,
      :pol     => 3256,
      :vpol    => 3258,
# 3260-3501 is unsigned, need to skip 242 byte
      :fltf    => 3502,  # 3502-3503
      :netfh   => 3504); # 3504-3505
  # then skip 94 byte to the start point of trace header

function init_FileHeader()
	  fhdr = FileHeader(0,0,0,0,0,0,0,0,0,0,
	                    0,0,0,0,0,0,0,0,0,0,
	                    0,0,0,0,0,0,0,0,0,)
	  return fhdr
end

function show(io::IO, fhdr::FileHeader)
    @printf("number of traces per ensemble                               : %6d\n", fhdr.ntrpe)
    @printf("sample interval in microseconds                             : %6d\n", fhdr.dt   )
    @printf("sample interval in microseconds in original field recording : %6d\n", fhdr.dtfr )
    @printf("number of samples per data trace                            : %6d\n", fhdr.ns   )
    @printf("number of samples per data tracf in original field recording: %6d\n", fhdr.nsfr )

    if fhdr.fmtc == 1
       @printf("data sample format                                          : IBM\n"            )
    elseif fhdr.fmtc == 5
       @printf("data sample format                                          : IEEE\n"            )
    end

    @printf("expected number of trace per trace ensemble                 : %6d\n", fhdr.expf )

    if fhdr.trsc == 1
       @printf("trace sorting code                                          : as recorded\n"      )
    elseif fhdr.trsc == 2
       @printf("trace sorting code                                          : CDP\n"      )
    elseif fhdr.trsc == 3
       @printf("trace sorting code                                          : single fold profile\n")
    elseif fhdr.trsc == 4
       @printf("trace sorting code                                          : post stack\n"      )
    elseif fhdr.trsc == 5
       @printf("trace sorting code                                          : common shot point\n")
    elseif fhdr.trsc == 6
       @printf("trace sorting code                                          : common receiver point\n")
    elseif fhdr.trsc == 7
       @printf("trace sorting code                                          : common offset profile\n")
    elseif fhdr.trsc == 8
       @printf("trace sorting code                                          : CMP\n"                  )
    elseif fhdr.trsc == 9
       @printf("trace sorting code                                          : common conversion point\n")
    end

    @printf("sweep frequency at start (Hz)                               : %6d\n", fhdr.sfs  )
    @printf("sweep frequency at end (Hz)                                 : %6d\n", fhdr.sfe  )
    @printf("sweep length (ms)                                           : %6d\n", fhdr.slen )

    if fhdr.styp == 1
       @printf("sweep type                                                  : linear\n"           )
    elseif fhdr.styp == 2
       @printf("sweep type                                                  : parabolic\n"        )
    elseif fhdr.styp == 3
       @printf("sweep type                                                  : exponential\n"      )
    end

    if fhdr.corr == 1
       @printf("correlated trace                                            : NO\n"      )
    elseif fhdr.corr == 2
       @printf("correlated trace                                            : YES\n"      )
    end

    if fhdr.rgc == 1
       @printf("gain recovered                                              : YES\n"      )
    elseif fhdr.rgc == 2
       @printf("gain recovered                                              : NO\n"      )
    end

    if fhdr.arm == 1
       @printf("amplitue recover method                                     : None\n"      )
    elseif fhdr.arm == 2
       @printf("amplitue recover method                                     : spherical divergence\n")
    elseif fhdr.arm == 3
       @printf("amplitue recover method                                     : AGC\n"              )
    elseif fhdr.arm == 4
       @printf("amplitue recover method                                     : other\n"            )
    end

    if fhdr.unit == 1
       @printf("measurement system                                          : Meters\n"          )
    elseif fhdr.unit == 2
       @printf("measurement system                                          : Feet\n"          )
    end

    if fhdr.pol == 1
       @printf("polarity : increase pressure and upward movement use negative number\n"        )
    elseif fhdr.pol == 2
       @printf("polarity : increase pressure and upward movement use positive number\n"        )
    end

    if fhdr.netfh == 0
       @printf("number of extended text file header                         : none\n"             )
    elseif fhdr.netfh == -1
       @printf("number of extended text file header                         : variable\n"         )
    elseif fhdr.netfh > 0
       @printf("number of extended text file header                         : %6d\n", fhdr.netfh)
    end

    return nothing
end

"""
   read 400 byte file header of SEGY file.
"""
function read_file_header(path::String; swap_bytes=true)

    stream = open(path, "r")

    # seek to the start point of file header
		seek(stream, file_header_location[:jobid])

    # initialize a empty file header
		fhdr = init_FileHeader()

    # read file header items
    # There are unsigned space in segy file, for loop can't be used
    fhdr.jobid   = read(stream, typeof(fhdr.jobid))
    fhdr.linnum  = read(stream, typeof(fhdr.linnum))
    fhdr.renum   = read(stream, typeof(fhdr.renum))
    fhdr.ntrpe   = read(stream, typeof(fhdr.ntrpe))
    fhdr.natrpe  = read(stream, typeof(fhdr.natrpe))
    fhdr.dt      = read(stream, typeof(fhdr.dt))
    fhdr.dtfr    = read(stream, typeof(fhdr.dtfr))
    fhdr.ns      = read(stream, typeof(fhdr.ns))
    fhdr.nsfr    = read(stream, typeof(fhdr.nsfr))
    fhdr.fmtc    = read(stream, typeof(fhdr.fmtc))
    fhdr.expf    = read(stream, typeof(fhdr.expf))
    fhdr.trsc    = read(stream, typeof(fhdr.trsc))
    fhdr.vsumc   = read(stream, typeof(fhdr.vsumc))
    fhdr.sfs     = read(stream, typeof(fhdr.sfs))
    fhdr.sfe     = read(stream, typeof(fhdr.sfe))
    fhdr.slen    = read(stream, typeof(fhdr.slen))
    fhdr.styp    = read(stream, typeof(fhdr.styp))
    fhdr.tnumsc  = read(stream, typeof(fhdr.tnumsc))
    fhdr.stalens = read(stream, typeof(fhdr.stalens))
    fhdr.stalene = read(stream, typeof(fhdr.stalene))
    fhdr.taty    = read(stream, typeof(fhdr.taty))
    fhdr.corr    = read(stream, typeof(fhdr.corr))
    fhdr.rgc     = read(stream, typeof(fhdr.rgc))
    fhdr.arm     = read(stream, typeof(fhdr.arm))
    fhdr.unit    = read(stream, typeof(fhdr.unit))
    fhdr.pol     = read(stream, typeof(fhdr.pol))
    fhdr.vpol    = read(stream, typeof(fhdr.vpol))

    skip(stream, 242) # skip the unsigned space
    fhdr.fltf    = read(stream, typeof(fhdr.fltf))
    fhdr.netfh   = read(stream, typeof(fhdr.netfh))

    # convert endian
    if swap_bytes
    	 for field in fieldnames(FileHeader)
    			 setfield!(fhdr, field, bswap(getfield(fhdr, field)))
    	 end
    end

    # skip to the start point of trace header
    # skip(stream, 94)
		close(stream)

	  return fhdr
end

"""
   write 400 byte file header of SEGY file.
"""
function put_file_header(path::String, fhdr::FileHeader; swap_bytes=true)

    # the start point for file header
    stream = open(path, "a+")
    seek(stream, file_header_location[:jobid])

    # initialize a temporary file header
    # keep the original file header untouch
    fh = init_FileHeader()
    if swap_bytes
       for field in fieldnames(FileHeader)
           setfield!(fh, field, bswap(getfield(fhdr, field)))
       end
    else
       for field in fieldnames(FileHeader)
           setfield!(fh, field, getfield(fhdr, field))
       end
    end

    # write file header
    write(stream, fh.jobid)
    write(stream, fh.linnum)
    write(stream, fh.renum)
    write(stream, fh.ntrpe)
    write(stream, fh.natrpe)
    write(stream, fh.dt)
    write(stream, fh.dtfr)
    write(stream, fh.ns)
    write(stream, fh.nsfr)
    write(stream, fh.fmtc)
    write(stream, fh.expf)
    write(stream, fh.trsc)
    write(stream, fh.vsumc)
    write(stream, fh.sfs)
    write(stream, fh.sfe)
    write(stream, fh.slen)
    write(stream, fh.styp)
    write(stream, fh.tnumsc)
    write(stream, fh.stalens)
    write(stream, fh.stalene)
    write(stream, fh.taty)
    write(stream, fh.corr)
    write(stream, fh.rgc)
    write(stream, fh.arm)
    write(stream, fh.unit)
    write(stream, fh.pol)
    write(stream, fh.vpol)

    # assign 0 to 242 unsigned bytes
    write(stream, zeros(Int8, 242))

    # write the last two field
    write(stream, fh.fltf)
    write(stream, fh.netfh)

    # assign 0 to 94 unsigned bytes
    write(stream, zeros(Int8, 94))
    close(stream)

    return nothing
end

"""
* tracl: trace number within line
* tracr: trace number within file
* fldr : original field record number
* tracf: trace number within the original field record
* ep   : energy source point number
* cdp  : Ensemble number (i.e. CDP, CMP, CRP,)
* cdpt : trace number within the ensemble
* trid : trace identification code (-1=other, 0=Unknown, 1=seismic, 2=dead, 3=Dummy, 4=Time break, 5=Uphole, 6=sweep
         7=timing, 8=waterbreak, 9=near field gun signature, 10=far field gun signature, 11=pressure sensor, 12=vertical component,
         13=cross-line component, 14=in-line component, 15=rotated vertical component, 16=rotated Transversal, 17=rotated radial component,
         18=vibrator reation mass, 19=vibrator baseplate, 19=vibrator estimated ground force, 21=vibrator reference
         22=time-velocity pairs, 23=optional use)
* nva    : number of vertically summed traces yield this trace
* nhs    : number of horizontally summed traces yield this trace
* duse   : data use (1=production, 2=Test)
* offset : distance between source and receiver
* gelev  : receiver elevation(above datum is positive, under is negative)
* selev  : surface elevation at source
* sdepth : source depth below surface (positive number)
* gdel   : datum elevation at receiver group
* sdel   : datum elevation at source
* swdep  : water depth at source
* gwdep  : water depth at receiver
* scalel : scalar applied to all evevations and depths (1,10,100,1000, positive used as multiplier, negative used as divisor)
* scalco : scalar applied to all coordinate (same as above case)
* sx     : source x
* sy     : source y
* gx     : receiver x
* gy     : receiver y
* counit : coordinate unit (1=length, 2=seconds of arc, 3=decimal degrees, 4=DMS)
* wevel  : weathering velocity (ft/s or meter/s)
* swevel : subweathing velocity
* sut    : uphole time at source in milliseconds
* gut    : uphole time at receiver in milliseconds
* sstat  : source static correction in milliseconds
* gstat  : group  static correction in milliseconds
* tstat  : total static applied in milliseconds
* laga   : lag time A
* lagb   : lag time B
* delrt  : delay recording time in millliseconds
* muts   : mute time - start time in milliseconds
* mute   : mute time - end time in milliseconds
* ns     : number of samples in this trace
* dt     : sample interval in microseconds
* gain   : gain type of field instruments (1=fixed, 2=binary, 3=floating point)
* igc    : instrument gain constant (dB)
* igi    : instrument early or initial gain (dB)
* corr   : correlated (1=no 2=yes)
* sfs    : sweep frequency at start
* sfe    : sweep frequency at end
* slen   : sweep length
* styp   : sweep type (1=linear 2=parabolic 3=exponential 4=other)
* stas   : sweep trace taper length at start in milliseconds
* stae   : sweep trace taper length at end in milliseconds
* tatyp  : taper type (1=linear, 2=cos^2, 3=other)
* afilf  : alias filter frequency (Hz), if applied
* afils  : alias filter slope
* nofilf : Notch filter frequency
* nofils : Notch filter slope
* lcf    : low-cut frequency
* hcf    : high-cut frequency
* lcs    : low-cut slope
* hcs    : high-cut slope
* year
* day
* hour
* minute
* sec
* timbas : timme basis code
* trwf   : trace weighting factor
* grnors : geophone group number of roll switch position one
* grnofr : geophone group number of roll switch position one within original field
* grnlof : geophone group number of last trace within original field
* gaps   : gap size
* otrav  : over travel
* CDPx   : CDP x coordinate, scala for coordinate need to apply
* CDPy   : CDP y
* inlineNum : for 3D post-stack data, inline number
* crosslineNum : for 3D post-stack data, crossline number
* shotNum: shot point number for poststack case
* scalarOnShotNum : applied to shot point number (=0, not scaled)
* didf   : device identifier
* scalart: scalar applied to acquisition time
* stype  : source type(1=vv, 2=vc, 3=vi, 4=iv, 5=ic, 6=ii, 7=div, 8=dic, 9=dii{v=vibrator v=vertical c=crossline i=inline i=impulse di=distributed impulse})
          source energy direction and source measurement, source measurement unit need to be handled.
"""
mutable struct TraceHeader
		tracl   :: Int32
		tracr   :: Int32
		fldr    :: Int32
		tracf   :: Int32
		ep      :: Int32
		cdp     :: Int32
		cdpt    :: Int32
    trid    :: Int16
		nva     :: Int16
		nhs     :: Int16
		duse    :: Int16
		offset  :: Int32
		gelev   :: Int32
		selev   :: Int32
		sdepth  :: Int32
		gdel    :: Int32
		sdel    :: Int32
		swdep   :: Int32
		gwdep   :: Int32
		scalel  :: Int16
		scalco  :: Int16
		sx      :: Int32
		sy      :: Int32
		gx      :: Int32
		gy      :: Int32
		counit  :: Int16
		wevel   :: Int16
		swevel  :: Int16
		sut     :: Int16
		gut     :: Int16
		sstat   :: Int16
		gstat   :: Int16
		tstat   :: Int16
		laga    :: Int16
		lagb    :: Int16
		delrt   :: Int16
		muts    :: Int16
		mute    :: Int16
		ns      :: Int16
		dt      :: Int16
		gain    :: Int16
		igc     :: Int16
		igi     :: Int16
		corr    :: Int16
		sfs     :: Int16
		sfe     :: Int16
		slen    :: Int16
		styp    :: Int16
		stas    :: Int16
		stae    :: Int16
		tatyp   :: Int16
		afilf   :: Int16
		afils   :: Int16
		nofilf  :: Int16
		nofils  :: Int16
		lcf     :: Int16
		hcf     :: Int16
		lcs     :: Int16
		hcs     :: Int16
		year    :: Int16
		day     :: Int16
		hour    :: Int16
		minute  :: Int16
		sec     :: Int16
		timbas  :: Int16
		trwf    :: Int16
		grnors  :: Int16
		grnofr  :: Int16
		grnlof  :: Int16
		gaps    :: Int16
		otrav   :: Int16
		CDPx    :: Int32
		CDPy    :: Int32
		inlineNum       :: Int32
		crosslineNum    :: Int32
		shotNum         :: Int32
		scalarOnShotNum :: Int16
		tvmu            :: Int16
    didf            :: Int16
		scalart         :: Int16
		stype           :: Int16
end

"""
   write trace's header into binary file
"""
function write_traces_header(path_hdr::String, thdr::Vector{TraceHeader}; print_interval=10000)

     # open file for writting
     fp = open(path_hdr, "w")

     # total number of traces
     num = length(thdr)
     write(fp, num)

     # loop over traces
     for i = 1 : num
         for field in fieldnames(TraceHeader)
             write(fp, getfield(thdr[i], field))
         end

         if (i % print_interval) == 0
            println("$i")
         end
     end
     close(fp)
end

"""
   read traces header from binary file
"""
function read_traces_header(path_hdr::String; print_interval=10000)

    fp = open(path_hdr, "r")

    # total number of traces
    num = read(fp, Int64)
    println("total number of traces: $num")

    # allocate memory
    thdr = Vector{TraceHeader}(undef, num)
    for i = 1 : num
        h = init_TraceHeader()
        for field in fieldnames(TraceHeader)
            setfield!(h, field, read(fp, typeof(getfield(h, field))))
        end
        thdr[i] = h

        if (i % print_interval) == 0
           println("$i")
        end
    end

    return thdr
end

"""
   extract the value of selected fields of traces header, the field is specified
by fields and the output is an 2D array of Int32.
"""
function extract_field_trace_header(thdr::Vector{TraceHeader}, fields::Vector{Symbol};
         istart=1, iend=10000, idx=Vector{Int}(undef, 0), print_interval=10000)

    # number of fields will be extracted
    num_field = length(fields)

    # number of traces
    if length(idx) == 0
       num_traces = iend - istart + 1
       idx = collect(istart : iend)
    else
       num_traces = length(idx)
    end

    # allocate memory to save the output
    dout = zeros(Int32, num_traces, num_field)

    # loop over the selected traces
    icount = 1
    for i in idx

        # loop over the fields
        for j = 1 : num_field
            dout[icount, j] = getfield(thdr[i], fields[j])
        end
        icount = icount + 1

        #print info
        if (icount % print_interval) == 0
           println("finished $icount traces")
        end
    end

    # return the value
    return dout
end

"""
   print the value of selected fields of trace header
"""
function print_field_trace_header(thdr::Vector{TraceHeader}, fields::Vector{Symbol};
         istart=1, iend=10000, idx=Vector{Int}(undef,0))

    # number of fields will be printed
    num_field = length(fields)

    # number of traces
    if length(idx) == 0
       idx = collect(istart : iend)
    end

    # loop over the selected traces
    for i in idx

        @printf("trace: %8d ", i)
        # loop over the field
        for j = 1 : num_field
            print(fields[j])
            val = getfield(thdr[i], fields[j])
            @printf(": %8d ", val)
        end
        @printf("\n")
    end
    return nothing
end

"""
   save the byte location of trace header
"""
const trace_header_location = Dict(
    :tracl  => 0,
    :tracr  => 4,
    :fldr   => 8,
    :tracf  => 12,
    :ep     => 16,
    :cdp    => 20,
    :cdpt   => 24,
    :trid   => 28,
    :nva    => 30,
    :nhs    => 32,
    :duse   => 34,
    :offset => 36,
    :gelev  => 40,
    :selev  => 44,
    :sdepth => 48,
    :gdel   => 52,
    :sdel   => 56,
    :swdep  => 60,
    :gwdep  => 64,
    :scalel => 68,
    :scalco => 70,
    :sx     => 72,
    :sy     => 76,
    :gx     => 80,
    :gy     => 84,
    :counit => 88,
    :wevel  => 90,
    :swevel => 92,
    :sut    => 94,
    :gut    => 96,
    :sstat  => 98,
    :gstat  => 100,
    :tstat  => 102,
    :laga   => 104,
    :lagb   => 106,
    :delrt  => 108,
    :muts   => 110,
    :mute   => 112,
    :ns     => 114,
    :dt     => 116,
    :gain   => 118,
    :igc    => 120,
    :igi    => 122,
    :corr   => 124,
    :sfs    => 126,
    :sfe    => 128,
    :slen   => 130,
    :styp   => 132,
    :stas   => 134,
    :stae   => 136,
    :tatyp  => 138,
    :afilf  => 140,
    :afils  => 142,
    :nofilf => 144,
    :nofils => 146,
    :lcf    => 148,
    :hcf    => 150,
    :lcs    => 152,
    :hcs    => 154,
    :year   => 156,
    :day    => 158,
    :hour   => 160,
    :minute => 162,
    :sec    => 164,
    :timbas => 166,
    :trwf   => 168,
    :grnors => 170,
    :grnofr => 172,
    :grnlof => 174,
    :gaps   => 176,
    :otrav  => 178,
    :CDPx   => 180,
    :CDPy   => 184,
    :inlineNum       => 188,
    :crosslineNum    => 192,
    :shotNum         => 196,
    :scalarOnShotNum => 200,
    :tvmu            => 202,
# 10 bytes gap
    :didf            => 214,
    :scalart         => 216,
    :stype           => 218)
# 20 bytes gap

# there are 81 fields
function init_TraceHeader()
	  h = TraceHeader(
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0,0,0,0,0,0,0,0,0,0,
	  0)

	  return h
end

# show the content of trace header
function show(io::IO, h::TraceHeader)
    @printf("offset(distance between source and receiver)                : %6d\n", h.offset)
    @printf("receiver elevation (above datum is positive)                : %6d\n", h.gelev )
    @printf("source elevation (above datum is positive)                  : %6d\n", h.selev )
    @printf("source depth below surface(positve number)                  : %6d\n", h.sdepth)
    @printf("datum elevation at receiver group                           : %6d\n", h.gdel  )
    @printf("datum elevation at source location                          : %6d\n", h.sdel  )
    @printf("scalar applied to all evevations and depths(1,10,100...)    : %6d\n", h.scalel)
    @printf("scalar applied to all coordinate                            : %6d\n", h.scalco)
    @printf("source x                                                    : %6d\n", h.sx    )
    @printf("source y                                                    : %6d\n", h.sy    )
    @printf("receiver x                                                  : %6d\n", h.gx    )
    @printf("receiver y                                                  : %6d\n", h.gy    )

    if h.counit == 1
       @printf("coordinate unit                                             : length\n")
    elseif h.counit == 2
       @printf("coordinate unit                                             : seconds of arc\n")
    elseif h.counit == 3
       @printf("coordinate unit                                             : decimal degrees\n")
    elseif h.counit == 4
       @printf("coordinate unit                                             : DMS\n")
    end

    @printf("source static correction in milliseconds                    : %6d\n", h.sstat  )
    @printf("group  static correction in milliseconds                    : %6d\n", h.gstat  )
    @printf("total static applied in milliseconds                        : %6d\n", h.tstat  )
    @printf("number of samples in this trace                             : %6d\n", h.ns     )
    @printf("sample interval in microseconds                             : %6d\n", h.dt     )
    @printf("CDP x coordinate, scala for coordinate need to apply        : %6d\n", h.CDPx   )
    @printf("CDP y coordinate, scala for coordinate need to apply        : %6d\n", h.CDPy   )
    @printf("for 3D post-stack data, inline number                       : %6d\n", h.inlineNum   )
    @printf("for 3D post-stack data, crossline number                    : %6d\n", h.crosslineNum)

    return nothing
end

"""
   read the header of whole traces, return vector of TraceHeader
"""
function extract_traces_header(path::String; swap_bytes=true, print_interval=10000)

		file_header = read_file_header(path; swap_bytes=swap_bytes)
		println("===================================================================\n")
		println("=======================Print file header===========================\n")
    print(file_header)

    if file_header.netfh == -1
			 error("Does not support variable text file header")
    elseif file_header.netfh == 0
       file_hsize = 3600
    elseif file_header.netfh > 0
       file_hsize = 400 + 3200 * file_header.netfh
		end

    # compute the number of traces in this file
    file_size = filesize(path)
    num_samples = Int64(file_header.ns)
    trace_data_size = num_samples * 4
    num_traces  = round(Int64, (file_size - file_hsize) / (240 + trace_data_size))

    @printf("number of traces in this file                               : %6d\n", num_traces)

    fid = open(path, "r")
    seek(fid, file_hsize)

    vector_trace_header = Vector{TraceHeader}(undef, num_traces)
		for i = 1 : num_traces
        h = init_TraceHeader()
        h.tracl  = read(fid, Int32)
        h.tracr  = read(fid, Int32)
        h.fldr   = read(fid, Int32)
        h.tracf  = read(fid, Int32)
        h.ep     = read(fid, Int32)
        h.cdp    = read(fid, Int32)
        h.cdpt   = read(fid, Int32)
        h.trid   = read(fid, Int16)
        h.nva    = read(fid, Int16)
        h.nhs    = read(fid, Int16)
        h.duse   = read(fid, Int16)
        h.offset = read(fid, Int32)
        h.gelev  = read(fid, Int32)
        h.selev  = read(fid, Int32)
        h.sdepth = read(fid, Int32)
        h.gdel   = read(fid, Int32)
        h.sdel   = read(fid, Int32)
        h.swdep  = read(fid, Int32)
        h.gwdep  = read(fid, Int32)
        h.scalel = read(fid, Int16)
        h.scalco = read(fid, Int16)
        h.sx     = read(fid, Int32)
        h.sy     = read(fid, Int32)
        h.gx     = read(fid, Int32)
        h.gy     = read(fid, Int32)
        h.counit = read(fid, Int16)
        h.wevel  = read(fid, Int16)
        h.swevel = read(fid, Int16)
        h.sut    = read(fid, Int16)
        h.gut    = read(fid, Int16)
        h.sstat  = read(fid, Int16)
        h.gstat  = read(fid, Int16)
        h.tstat  = read(fid, Int16)
        h.laga   = read(fid, Int16)
        h.lagb   = read(fid, Int16)
        h.delrt  = read(fid, Int16)
        h.muts   = read(fid, Int16)
        h.mute   = read(fid, Int16)
        h.ns     = read(fid, Int16)
        h.dt     = read(fid, Int16)
        h.gain   = read(fid, Int16)
        h.igc    = read(fid, Int16)
        h.igi    = read(fid, Int16)
        h.corr   = read(fid, Int16)
        h.sfs    = read(fid, Int16)
        h.sfe    = read(fid, Int16)
        h.slen   = read(fid, Int16)
        h.styp   = read(fid, Int16)
        h.stas   = read(fid, Int16)
        h.stae   = read(fid, Int16)
        h.tatyp  = read(fid, Int16)
        h.afilf  = read(fid, Int16)
        h.afils  = read(fid, Int16)
        h.nofilf = read(fid, Int16)
        h.nofils = read(fid, Int16)
        h.lcf    = read(fid, Int16)
        h.hcf    = read(fid, Int16)
        h.lcs    = read(fid, Int16)
        h.hcs    = read(fid, Int16)
        h.year   = read(fid, Int16)
        h.day    = read(fid, Int16)
        h.hour   = read(fid, Int16)
        h.minute = read(fid, Int16)
        h.sec    = read(fid, Int16)
        h.timbas = read(fid, Int16)
        h.trwf   = read(fid, Int16)
        h.grnors = read(fid, Int16)
        h.grnofr = read(fid, Int16)
        h.grnlof = read(fid, Int16)
        h.gaps   = read(fid, Int16)
        h.otrav  = read(fid, Int16)
        h.CDPx            = read(fid, Int32)
        h.CDPy            = read(fid, Int32)
        h.inlineNum       = read(fid, Int32)
        h.crosslineNum    = read(fid, Int32)
        h.shotNum         = read(fid, Int32)
        h.scalarOnShotNum = read(fid, Int16)
        h.tvmu            = read(fid, Int16)
        skip(fid, 10) # skip 10 bytes
        h.didf            = read(fid, Int16)
        h.scalart         = read(fid, Int16)
        h.stype           = read(fid, Int16)
        skip(fid, 20 + trace_data_size) # skip 20 bytes and trace data to the begining of next trace

        if swap_bytes
           for field in fieldnames(TraceHeader)
               setfield!(h, field, bswap(getfield(h, field)))
           end
        end
        vector_trace_header[i] = h

        if (i % print_interval) == 0
           println("trace: $i")
        end
	  end

    close(fid)
    return vector_trace_header

end

"""
   read the segy file, return FileHeader, vector of TraceHeader and matrix of data
"""
function read_segy_file(path::String; swap_bytes=true, data_format="NULL")

    text_header = read_text_header(path)
    println("===================================================================\n")
    println("=======================Print text header===========================\n")
    print_text_header(text_header)

    file_header = read_file_header(path; swap_bytes=swap_bytes)
    println("===================================================================\n")
    println("=======================Print file header===========================\n")
    print(file_header)

    if file_header.netfh == -1
       error("Does not support variable text file header")
    elseif file_header.netfh == 0
       file_hsize = 3600
    elseif file_header.netfh > 0
       file_hsize = 400 + 3200 * file_header.netfh
    end

    # compute the number of traces in this file
    file_size = filesize(path)
    num_samples = Int64(file_header.ns)
    trace_data_size = num_samples * 4
    num_traces  = round(Int64, (file_size - file_hsize) / (240 + trace_data_size))
    @printf("number of traces in this file                               : %6d\n", num_traces)

    # determine the data format
    if data_format == "NULL"
       if file_header.fmtc == 1
          data_format = "IBM"
          tmp_trace   = Vector{IBMFloat32}(undef, num_samples)
          @printf("data format determined from file header                  : IBM\n")
       elseif file_header.fmtc == 5
          data_format = "IEEE"
          @printf("data format determined from file header                  : IEEE\n")
          tmp_trace   = Vector{Float32}(undef, num_samples)
       end
    else
       @printf("data format determined from keyword input                : %s\n", data_format)
    end

    # the start location of the first trace
    fid = open(path, "r")
    seek(fid, file_hsize)

    # allocate space for trace header and data
    trace_header = Vector{TraceHeader}(undef, num_traces)
    data         = zeros(Float32, num_samples, num_traces)
    idx          = 1

    for i = 1 : num_traces

        # read each item of trace header
        h = init_TraceHeader()
        h.tracl  = read(fid, Int32)
        h.tracr  = read(fid, Int32)
        h.fldr   = read(fid, Int32)
        h.tracf  = read(fid, Int32)
        h.ep     = read(fid, Int32)
        h.cdp    = read(fid, Int32)
        h.cdpt   = read(fid, Int32)
        h.trid   = read(fid, Int16)
        h.nva    = read(fid, Int16)
        h.nhs    = read(fid, Int16)
        h.duse   = read(fid, Int16)
        h.offset = read(fid, Int32)
        h.gelev  = read(fid, Int32)
        h.selev  = read(fid, Int32)
        h.sdepth = read(fid, Int32)
        h.gdel   = read(fid, Int32)
        h.sdel   = read(fid, Int32)
        h.swdep  = read(fid, Int32)
        h.gwdep  = read(fid, Int32)
        h.scalel = read(fid, Int16)
        h.scalco = read(fid, Int16)
        h.sx     = read(fid, Int32)
        h.sy     = read(fid, Int32)
        h.gx     = read(fid, Int32)
        h.gy     = read(fid, Int32)
        h.counit = read(fid, Int16)
        h.wevel  = read(fid, Int16)
        h.swevel = read(fid, Int16)
        h.sut    = read(fid, Int16)
        h.gut    = read(fid, Int16)
        h.sstat  = read(fid, Int16)
        h.gstat  = read(fid, Int16)
        h.tstat  = read(fid, Int16)
        h.laga   = read(fid, Int16)
        h.lagb   = read(fid, Int16)
        h.delrt  = read(fid, Int16)
        h.muts   = read(fid, Int16)
        h.mute   = read(fid, Int16)
        h.ns     = read(fid, Int16)
        h.dt     = read(fid, Int16)
        h.gain   = read(fid, Int16)
        h.igc    = read(fid, Int16)
        h.igi    = read(fid, Int16)
        h.corr   = read(fid, Int16)
        h.sfs    = read(fid, Int16)
        h.sfe    = read(fid, Int16)
        h.slen   = read(fid, Int16)
        h.styp   = read(fid, Int16)
        h.stas   = read(fid, Int16)
        h.stae   = read(fid, Int16)
        h.tatyp  = read(fid, Int16)
        h.afilf  = read(fid, Int16)
        h.afils  = read(fid, Int16)
        h.nofilf = read(fid, Int16)
        h.nofils = read(fid, Int16)
        h.lcf    = read(fid, Int16)
        h.hcf    = read(fid, Int16)
        h.lcs    = read(fid, Int16)
        h.hcs    = read(fid, Int16)
        h.year   = read(fid, Int16)
        h.day    = read(fid, Int16)
        h.hour   = read(fid, Int16)
        h.minute = read(fid, Int16)
        h.sec    = read(fid, Int16)
        h.timbas = read(fid, Int16)
        h.trwf   = read(fid, Int16)
        h.grnors = read(fid, Int16)
        h.grnofr = read(fid, Int16)
        h.grnlof = read(fid, Int16)
        h.gaps   = read(fid, Int16)
        h.otrav  = read(fid, Int16)
        h.CDPx            = read(fid, Int32)
        h.CDPy            = read(fid, Int32)
        h.inlineNum       = read(fid, Int32)
        h.crosslineNum    = read(fid, Int32)
        h.shotNum         = read(fid, Int32)
        h.scalarOnShotNum = read(fid, Int16)
        h.tvmu            = read(fid, Int16)
        skip(fid, 10) # skip 10 bytes
        h.didf            = read(fid, Int16)
        h.scalart         = read(fid, Int16)
        h.stype           = read(fid, Int16)
        skip(fid, 20) # skip 20 bytes to the begining of trace data

        if swap_bytes
           for field in fieldnames(TraceHeader)
               setfield!(h, field, bswap(getfield(h, field)))
           end
        end
        trace_header[i] = h

        # read trace data
        read!(fid, tmp_trace)
        if data_format == "IEEE"
           copyto!(data, idx, tmp_trace, 1, num_samples)
		    elseif data_format == "IBM"
           copyto!(data, idx, convert(Vector{Float32}, tmp_trace), 1, num_samples)
		    end
        idx = idx + num_samples
    end

    close(fid)
    return text_header, file_header, trace_header, data

end

"""
   write segy file with inputed text header, file header, vector of trace header and data,
only support save the data in IEEE format (different to the default IBM format).
"""
function write_segy_file(path::String, text_header::Vector{Char}, fhdr::FileHeader,
         thdr::Union{Vector{TraceHeader}, TraceHeader}, data; swap_bytes=true, data_format="IEEE")

    # write text header
    put_text_header(path, text_header)

    # check number of samples per trace
    num_samples = size(data, 1)
    if num_samples != fhdr.ns
       fhdr.ns = num_samples
    end

    # write file header
    fhdr.fmtc = 5  # make sure the data format is IEEE
    put_file_header(path, fhdr)

    # convert data to Float32
    if eltype(data) != Float32
       data = convert(Array{Float32}, data)
    end

    # the start location trace header
    file_hsize = 3600
    fp = open(path, "a+")
    seek(fp, file_hsize)

    # write only 1 trace
    if typeof(thdr) == TraceHeader
       if prod(size(data)[2:end]) != 1
          error("input more than one trace")
       end
       write_one_trace(fp, thdr, data)

    # write multiple traces
   elseif typeof(thdr) == Vector{TraceHeader}
       num_traces = length(thdr)
       if num_traces != prod(size(data)[2:end])
          error("number of trace header not equal to the number of data trace")
       end
       data  = reshape(data, num_samples, num_traces)
       for i = 1 : num_traces
           write_one_trace(fp, thdr[i], data[:,i])
       end

    # wronge input
    else
       error("input trace header is wrong")
    end

    close(fp)

    return nothing
end

"""
   write one trace header and data to segy file
"""
function write_one_trace(fid::IOStream, h::TraceHeader, d::Vector{Float32})

    #write the field of one trace header with swaped byte
    write(fid, bswap(h.tracl          ))
    write(fid, bswap(h.tracr          ))
    write(fid, bswap(h.fldr           ))
    write(fid, bswap(h.tracf          ))
    write(fid, bswap(h.ep             ))
    write(fid, bswap(h.cdp            ))
    write(fid, bswap(h.cdpt           ))
    write(fid, bswap(h.trid           ))
    write(fid, bswap(h.nva            ))
    write(fid, bswap(h.nhs            ))
    write(fid, bswap(h.duse           ))
    write(fid, bswap(h.offset         ))
    write(fid, bswap(h.gelev          ))
    write(fid, bswap(h.selev          ))
    write(fid, bswap(h.sdepth         ))
    write(fid, bswap(h.gdel           ))
    write(fid, bswap(h.sdel           ))
    write(fid, bswap(h.swdep          ))
    write(fid, bswap(h.gwdep          ))
    write(fid, bswap(h.scalel         ))
    write(fid, bswap(h.scalco         ))
    write(fid, bswap(h.sx             ))
    write(fid, bswap(h.sy             ))
    write(fid, bswap(h.gx             ))
    write(fid, bswap(h.gy             ))
    write(fid, bswap(h.counit         ))
    write(fid, bswap(h.wevel          ))
    write(fid, bswap(h.swevel         ))
    write(fid, bswap(h.sut            ))
    write(fid, bswap(h.gut            ))
    write(fid, bswap(h.sstat          ))
    write(fid, bswap(h.gstat          ))
    write(fid, bswap(h.tstat          ))
    write(fid, bswap(h.laga           ))
    write(fid, bswap(h.lagb           ))
    write(fid, bswap(h.delrt          ))
    write(fid, bswap(h.muts           ))
    write(fid, bswap(h.mute           ))
    write(fid, bswap(h.ns             ))
    write(fid, bswap(h.dt             ))
    write(fid, bswap(h.gain           ))
    write(fid, bswap(h.igc            ))
    write(fid, bswap(h.igi            ))
    write(fid, bswap(h.corr           ))
    write(fid, bswap(h.sfs            ))
    write(fid, bswap(h.sfe            ))
    write(fid, bswap(h.slen           ))
    write(fid, bswap(h.styp           ))
    write(fid, bswap(h.stas           ))
    write(fid, bswap(h.stae           ))
    write(fid, bswap(h.tatyp          ))
    write(fid, bswap(h.afilf          ))
    write(fid, bswap(h.afils          ))
    write(fid, bswap(h.nofilf         ))
    write(fid, bswap(h.nofils         ))
    write(fid, bswap(h.lcf            ))
    write(fid, bswap(h.hcf            ))
    write(fid, bswap(h.lcs            ))
    write(fid, bswap(h.hcs            ))
    write(fid, bswap(h.year           ))
    write(fid, bswap(h.day            ))
    write(fid, bswap(h.hour           ))
    write(fid, bswap(h.minute         ))
    write(fid, bswap(h.sec            ))
    write(fid, bswap(h.timbas         ))
    write(fid, bswap(h.trwf           ))
    write(fid, bswap(h.grnors         ))
    write(fid, bswap(h.grnofr         ))
    write(fid, bswap(h.grnlof         ))
    write(fid, bswap(h.gaps           ))
    write(fid, bswap(h.otrav          ))
    write(fid, bswap(h.CDPx           ))
    write(fid, bswap(h.CDPy           ))
    write(fid, bswap(h.inlineNum      ))
    write(fid, bswap(h.crosslineNum   ))
    write(fid, bswap(h.shotNum        ))
    write(fid, bswap(h.scalarOnShotNum))
    write(fid, bswap(h.tvmu           ))

    # write 5 Int16(0) to fill the 10 byte gap
    write(fid, zeros(Int8, 10))

    # write another 3 fields
    write(fid, bswap(h.didf           ))
    write(fid, bswap(h.scalart        ))
    write(fid, bswap(h.stype          ))

    # write 5 Int32(0) to fill the 20 byte gap
    write(fid, zeros(Int8, 20))

    # write the seismic data
    write(fid, d)

    # flush the IOBuffer
    flush(fid)

    return nothing
end

"""
   extract unique locations of sources or receivers
"""
function extract_sensors_location(thdr::Vector{TraceHeader}; tol = 1.0, print_interval=10000)

    num_traces = length(thdr)

    # the vector to save unique source location
    sx = [thdr[1].sx]
    sy = [thdr[1].sy]

    # the vector to save unique receiver location
    gx = [thdr[1].gx]
    gy = [thdr[1].gy]

    # the vector to save number of traces per shot
    num_receiver = Vector{Int64}(undef, 0)

    # temporary variable for sensors' location difference
    res = zeros(Float32, 2)

    # the distance between sensors
    dis = Inf32

    for i = 2 : num_traces

        # source location
        dis = Inf32
        for j = 1 : length(sx)
            res[1] = sx[j] - thdr[i].sx
            res[2] = sy[j] - thdr[i].sy
            tmp = norm(res)
            if tmp < dis
               dis = tmp
            end
        end
        if dis > tol
           push!(sx, thdr[i].sx)
           push!(sy, thdr[i].sy)
        end

        # receiver location
        dis = Inf
        for j = 1 : length(gx)
            res[1] = gx[j] - thdr[i].gx
            res[2] = gy[j] - thdr[i].gy
            tmp = norm(res)
            if tmp < dis
               dis = tmp
            end
        end
        if dis > tol
           push!(gx, thdr[i].gx)
           push!(gy, thdr[i].gy)
        end

        if (i % print_interval) == 0
           println("trace number: $i")
        end
    end

    return sx, sy, gx, gy
end

"""
   create acquisition geometry from trace header, currently only work for data sorted in shot gather
num_component used for multi-component seismic data.
"""
function create_acquisition_geometry(traces_header::Vector{TraceHeader}; num_component=1)

    num_traces = length(traces_header)

    # determine the start and end trace index for each shot
    trace_start = [1]
    trace_end   = Vector{Int64}(undef, 0)
    for i  = 2 : num_traces
        if traces_header[i].sx != traces_header[i-1].sx || traces_header[i].sy != traces_header[i-1].sy
           push!(trace_end, i-1)
           push!(trace_start, i)
        end
    end
    push!(trace_end, num_traces)

    # determine the number of traces for each shot
    num_shots = length(trace_start)
    num_traces_per_shot = zeros(Int64, num_shots)
    for i = 1 : num_shots
        tmp = trace_end[i] - trace_start[i] + 1
        num_traces_per_shot[i] = convert(Int64, tmp / num_component)
    end

    # get the source coordinate
    data_format = typeof(traces_header[1].sx)
    source_x = zeros(data_format, num_shots)
    source_y = zeros(data_format, num_shots)

    # get the receiver coordinate for each shot
    data_format = typeof(traces_header[1].gx)
    receiver_x = Vector{Vector{data_format}}(undef, num_shots)
    receiver_y = Vector{Vector{data_format}}(undef, num_shots)

    for i = 1 : num_shots
        idx = trace_start[i]
        source_x[i] = traces_header[idx].sx
        source_y[i] = traces_header[idx].sy
        tmp_x = zeros(data_format, num_traces_per_shot[i])
        tmp_y = zeros(data_format, num_traces_per_shot[i])
        count = 0
        for j = trace_start[i] : trace_end[i]
            if (j % num_component) == 0
               count = count + 1
               tmp_x[count] = traces_header[j].gx
               tmp_y[count] = traces_header[j].gy
            end
        end
        receiver_x[i] = tmp_x
        receiver_y[i] = tmp_y
    end

    return source_x, source_y, receiver_x, receiver_y
end

"""
   create a lookup table to facilitate the reading of gathers(can be shot gather or CMP gathers)
"""
function create_lookup_table(path_txt::String, path::String; swap_bytes=true, data_format="NULL", print_interval=100000)

    # print file header
    file_header = read_file_header(path; swap_bytes=swap_bytes)
    println("===================================================================\n")
    println("=======================Print file header===========================\n")
    print(file_header)

    # compute the size of file header in terms of number of byte
    if file_header.netfh == -1
       error("Does not support variable text file header")
    elseif file_header.netfh == 0
       file_hsize = 3600
    elseif file_header.netfh > 0
       file_hsize = 400 + 3200 * file_header.netfh
    end
    @printf("the size of file header is                                  : %6d\n", file_hsize)

    # compute the number of traces in this file
    file_size = filesize(path)
    num_samples = Int64(file_header.ns)
    trace_data_size = num_samples * 4
    num_traces  = round(Int64, (file_size - file_hsize) / (240 + trace_data_size))
    @printf("number of traces in this file                               : %6d\n", num_traces)

    # determine the data format
    if data_format == "NULL"
       if file_header.fmtc == 1
          data_format = "IBM"
          @printf("data format determined from file header                     : IBM\n")
       elseif file_header.fmtc == 5
          data_format = "IEEE"
          @printf("data format determined from file header                     : IEEE\n")
       end
    else
       @printf("data format determined from keyword input                   : %s\n", data_format)
    end

    # skip text header and file header
    fid = open(path, "r")
    seek(fid, file_hsize)

    # read the header of all traces
    trace_header = Vector{TraceHeader}(undef, num_traces)
    for i = 1 : num_traces
        h = init_TraceHeader()
        h.tracl  = read(fid, Int32)
        h.tracr  = read(fid, Int32)
        h.fldr   = read(fid, Int32)
        h.tracf  = read(fid, Int32)
        h.ep     = read(fid, Int32)
        h.cdp    = read(fid, Int32)
        h.cdpt   = read(fid, Int32)
        h.trid   = read(fid, Int16)
        h.nva    = read(fid, Int16)
        h.nhs    = read(fid, Int16)
        h.duse   = read(fid, Int16)
        h.offset = read(fid, Int32)
        h.gelev  = read(fid, Int32)
        h.selev  = read(fid, Int32)
        h.sdepth = read(fid, Int32)
        h.gdel   = read(fid, Int32)
        h.sdel   = read(fid, Int32)
        h.swdep  = read(fid, Int32)
        h.gwdep  = read(fid, Int32)
        h.scalel = read(fid, Int16)
        h.scalco = read(fid, Int16)
        h.sx     = read(fid, Int32)
        h.sy     = read(fid, Int32)
        h.gx     = read(fid, Int32)
        h.gy     = read(fid, Int32)
        h.counit = read(fid, Int16)
        h.wevel  = read(fid, Int16)
        h.swevel = read(fid, Int16)
        h.sut    = read(fid, Int16)
        h.gut    = read(fid, Int16)
        h.sstat  = read(fid, Int16)
        h.gstat  = read(fid, Int16)
        h.tstat  = read(fid, Int16)
        h.laga   = read(fid, Int16)
        h.lagb   = read(fid, Int16)
        h.delrt  = read(fid, Int16)
        h.muts   = read(fid, Int16)
        h.mute   = read(fid, Int16)
        h.ns     = read(fid, Int16)
        h.dt     = read(fid, Int16)
        h.gain   = read(fid, Int16)
        h.igc    = read(fid, Int16)
        h.igi    = read(fid, Int16)
        h.corr   = read(fid, Int16)
        h.sfs    = read(fid, Int16)
        h.sfe    = read(fid, Int16)
        h.slen   = read(fid, Int16)
        h.styp   = read(fid, Int16)
        h.stas   = read(fid, Int16)
        h.stae   = read(fid, Int16)
        h.tatyp  = read(fid, Int16)
        h.afilf  = read(fid, Int16)
        h.afils  = read(fid, Int16)
        h.nofilf = read(fid, Int16)
        h.nofils = read(fid, Int16)
        h.lcf    = read(fid, Int16)
        h.hcf    = read(fid, Int16)
        h.lcs    = read(fid, Int16)
        h.hcs    = read(fid, Int16)
        h.year   = read(fid, Int16)
        h.day    = read(fid, Int16)
        h.hour   = read(fid, Int16)
        h.minute = read(fid, Int16)
        h.sec    = read(fid, Int16)
        h.timbas = read(fid, Int16)
        h.trwf   = read(fid, Int16)
        h.grnors = read(fid, Int16)
        h.grnofr = read(fid, Int16)
        h.grnlof = read(fid, Int16)
        h.gaps   = read(fid, Int16)
        h.otrav  = read(fid, Int16)
        h.CDPx            = read(fid, Int32)
        h.CDPy            = read(fid, Int32)
        h.inlineNum       = read(fid, Int32)
        h.crosslineNum    = read(fid, Int32)
        h.shotNum         = read(fid, Int32)
        h.scalarOnShotNum = read(fid, Int16)
        h.tvmu            = read(fid, Int16)
        skip(fid, 10) # skip 10 bytes
        h.didf            = read(fid, Int16)
        h.scalart         = read(fid, Int16)
        h.stype           = read(fid, Int16)
        # skip to the begining of the header of next trace bytes
        skip(fid, 20 + trace_data_size)

        # swap byte
        if swap_bytes
           for field in fieldnames(TraceHeader)
               setfield!(h, field, bswap(getfield(h, field)))
           end
        end

        # include in the vector
        trace_header[i] = h

        if (i % print_interval) == 0
           println("$i")
        end
    end
    close(fid)

    # count the number of shot gathers
    trace_start = [1]
    trace_end   = [ ]
    for i  = 2 : num_traces
        if trace_header[i].sx != trace_header[i-1].sx || trace_header[i].sy != trace_header[i-1].sy
           push!(trace_end, i-1)
           push!(trace_start, i)
        end
    end
    push!(trace_end, num_traces)
    num_shots = length(trace_start)

    # form a shot gather loopup table
    lookup = hcat(collect(1:num_shots), trace_start, trace_end)

    # create the lookup table into a text file
    fid = open(path_txt, "w")
    write(fid, join(["file_header_size = " "$file_hsize" "\n"]))
    write(fid, join(["num_samples      = " "$num_samples" "\n"]))
    write(fid, join(["num_traces       = " "$num_traces" "\n"]))
    write(fid, join(["num_shots        = " "$num_shots" "\n"]))
    write(fid, join(["swap_bytes       = " "$swap_bytes" "\n"]))
    write(fid, join(["data_format      = " data_format "\n"]))

    # write the trace index table
    writedlm(fid,  lookup)
    close(fid)
end

"""
   read the lookup table for shot gathers
   example: (file_header_size, num_samples, num_traces, swap_bytes, data_format, shot_idx) = read_lookup_table(path);
"""
function read_lookup_table(path::String)

    fid = open(path, "r")

    # read the file_header_size
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    file_header_size = parse(Int32, tmp[loc+2:end])

    # read the number of samples per trace
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    num_samples = parse(Int32, tmp[loc+2:end])

    # read the number of traces
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    num_traces = parse(Int32, tmp[loc+2:end])

    # read the number of shots
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    num_shots = parse(Int32, tmp[loc+2:end])

    # read the flag of byte swap
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    swap_bytes = parse(Bool, tmp[loc+2:end])

    # read the data format of the data sampale
    tmp = readline(fid); loc = findfirst(isequal('='), tmp);
    data_format = tmp[loc+2:end]

    # close the file
    close(fid)

    # read shot start and end index
    shot_idx = readdlm(path, skipstart=6)
    shot_idx = convert(Matrix{Int32}, shot_idx)

    if num_shots != size(shot_idx, 1)
       error("the information on number of shots doesn't match")
    end

    return file_header_size, num_samples, num_traces, swap_bytes, data_format, shot_idx

end

"""
   read one shot gather based on look up table
"""
function read_one_shot(path_sgy::String, path_txt::String, idx::Integer)

    # read the look up table
    (file_header_size, num_samples, num_traces, swap_bytes, data_format, shot_idx) = read_lookup_table(path_txt)

    # the size of a trace(header plus data)
    trace_size = 240 + num_samples * 4

    # start and end trace of this shot
    trace_start   = shot_idx[idx, 2]
    trace_end     = shot_idx[idx, 3]
    num_receivers = trace_end - trace_start + 1

    # all memory to save data of one trace
    if data_format == "IBM"
       tmp_trace = Vector{IBMFloat32}(undef, num_samples)
    elseif data_format == "IEEE"
       tmp_trace = Vector{Float32}(undef, num_samples)
    end

    # the start position
    location = file_header_size + (trace_start-1) * trace_size
    fid = open(path_sgy, "r")
    seek(fid, location)

    # read the shot gather
    trace_header = Vector{TraceHeader}(undef, num_receivers)
    data         = zeros(Float32, num_samples, num_receivers)
    idx          = 1
    for i = 1 : num_receivers
        h = init_TraceHeader()
        h.tracl  = read(fid, Int32)
        h.tracr  = read(fid, Int32)
        h.fldr   = read(fid, Int32)
        h.tracf  = read(fid, Int32)
        h.ep     = read(fid, Int32)
        h.cdp    = read(fid, Int32)
        h.cdpt   = read(fid, Int32)
        h.trid   = read(fid, Int16)
        h.nva    = read(fid, Int16)
        h.nhs    = read(fid, Int16)
        h.duse   = read(fid, Int16)
        h.offset = read(fid, Int32)
        h.gelev  = read(fid, Int32)
        h.selev  = read(fid, Int32)
        h.sdepth = read(fid, Int32)
        h.gdel   = read(fid, Int32)
        h.sdel   = read(fid, Int32)
        h.swdep  = read(fid, Int32)
        h.gwdep  = read(fid, Int32)
        h.scalel = read(fid, Int16)
        h.scalco = read(fid, Int16)
        h.sx     = read(fid, Int32)
        h.sy     = read(fid, Int32)
        h.gx     = read(fid, Int32)
        h.gy     = read(fid, Int32)
        h.counit = read(fid, Int16)
        h.wevel  = read(fid, Int16)
        h.swevel = read(fid, Int16)
        h.sut    = read(fid, Int16)
        h.gut    = read(fid, Int16)
        h.sstat  = read(fid, Int16)
        h.gstat  = read(fid, Int16)
        h.tstat  = read(fid, Int16)
        h.laga   = read(fid, Int16)
        h.lagb   = read(fid, Int16)
        h.delrt  = read(fid, Int16)
        h.muts   = read(fid, Int16)
        h.mute   = read(fid, Int16)
        h.ns     = read(fid, Int16)
        h.dt     = read(fid, Int16)
        h.gain   = read(fid, Int16)
        h.igc    = read(fid, Int16)
        h.igi    = read(fid, Int16)
        h.corr   = read(fid, Int16)
        h.sfs    = read(fid, Int16)
        h.sfe    = read(fid, Int16)
        h.slen   = read(fid, Int16)
        h.styp   = read(fid, Int16)
        h.stas   = read(fid, Int16)
        h.stae   = read(fid, Int16)
        h.tatyp  = read(fid, Int16)
        h.afilf  = read(fid, Int16)
        h.afils  = read(fid, Int16)
        h.nofilf = read(fid, Int16)
        h.nofils = read(fid, Int16)
        h.lcf    = read(fid, Int16)
        h.hcf    = read(fid, Int16)
        h.lcs    = read(fid, Int16)
        h.hcs    = read(fid, Int16)
        h.year   = read(fid, Int16)
        h.day    = read(fid, Int16)
        h.hour   = read(fid, Int16)
        h.minute = read(fid, Int16)
        h.sec    = read(fid, Int16)
        h.timbas = read(fid, Int16)
        h.trwf   = read(fid, Int16)
        h.grnors = read(fid, Int16)
        h.grnofr = read(fid, Int16)
        h.grnlof = read(fid, Int16)
        h.gaps   = read(fid, Int16)
        h.otrav  = read(fid, Int16)
        h.CDPx            = read(fid, Int32)
        h.CDPy            = read(fid, Int32)
        h.inlineNum       = read(fid, Int32)
        h.crosslineNum    = read(fid, Int32)
        h.shotNum         = read(fid, Int32)
        h.scalarOnShotNum = read(fid, Int16)
        h.tvmu            = read(fid, Int16)
        skip(fid, 10) # skip 10 bytes
        h.didf            = read(fid, Int16)
        h.scalart         = read(fid, Int16)
        h.stype           = read(fid, Int16)
        skip(fid, 20) # skip 20 bytes to the begining of trace data

        if swap_bytes
           for field in fieldnames(TraceHeader)
               setfield!(h, field, bswap(getfield(h, field)))
           end
        end
        trace_header[i] = h

        # read seismic data
        read!(fid, tmp_trace)
        if data_format == "IEEE"
           copyto!(data, idx, tmp_trace, 1, num_samples)
  		  elseif data_format == "IBM"
           copyto!(data, idx, convert(Vector{Float32}, tmp_trace), 1, num_samples)
  		  end
        idx = idx + num_samples

    end

    return trace_header, data
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays

sz=[3, 3, 3]
sz=[2, 2, 2]
# sz=[5, 6, 7]

# V,(VV,EV,FV,CV) = Lar.cuboidGrid([3,3,3], true)
V,(VV,EV,FV,CV) = Lar.cuboidGrid(sz, true)

# Plasm.view( Plasm.numbering(.6)((V,[VV, EV, FV])) )
Plasm.view( Plasm.numbering(.6)((V,[VV, EV, FV])) )

ext: jl


text:
basic_opts = [
("add",".+","s+xi"),
("sub",".-","s-xi"),
("mul",".*","s*xi"),
("div","./","s/xi"),
("pow",".^","pow(s,xi)"),
("max","max","(xi>s?xi:s)"),
("min","min","(xi<s?xi:s)"),
("eq",".==","xi==s"),
("gt",".>","xi>s"),
("ge",".>=","xi>=s"),
("lt",".<","xi<s"),
("le",".<=","xi<=s"),
]

function basic_cuda_src(f, j=f, ex="$f(s,xi)"; BLK=256, THR=256)
"""
__global__ void _$(f)_32_01(int n, float s, float *x, float *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    float xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_32_01(int n, float s, float *x, float *y) {
    _$(f)_32_01<<<$BLK,$THR>>>(n,s,x,y);
  }
}
__global__ void _$(f)_64_01(int n, double s, double *x, double *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    double xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_64_01(int n, double s, double *x, double *y) {
    _$(f)_64_01<<<$BLK,$THR>>>(n,s,x,y);
  }
}
"""
end

for a in basic_opts
    isa(a,Tuple) || (a=(a,))
    print(basic_cuda_src(a...))
end


text:
function format_header(header::String)
    tmp = ""
    for i in 1:80:length(header)
        tmp *= i+79 < length(header) ? header[i:i+79] : header[i:end]
        tmp *= "\n"
    end
    return tmp
end
function format_label(label::Tuple{Tuple{Real,Real},String})
    return ". (" * string.(label[1][1]) * "," * string.(label[1][2]) * ")\t" * label[2]
end

function form_image(img::AstroImage, indx::Integer)
    rimg = render(img,indx)
    img24 = Matrix(UInt32[reinterpret(UInt32, convert(RGB24, rimg[i,j])) for i = 1:size(rimg,1), j = size(rimg,2):-1:1]')
    fv = XY(0.0..size(img24,2), 0.0..size(img24,1))
    zr = Signal(ZoomRegion(fv, fv))
    surf = Cairo.CairoRGBSurface(img24)

    return fv,zr,surf
end

function get_button(name::String)
    btn =  ToolButton(name)
    set_gtk_property!(btn, :label, name)
    set_gtk_property!(btn, :is_important, true)
    return btn
end

get_label(name::String) = Label(name)

function header_window(header::WCSTransform)
    header = WCS.to_header(header)
    new_win = Window("Header",600,500)
    txt = get_label("Header")
    push!(new_win,txt)
    set_gtk_property!(txt, :label, format_header(header))
    return new_win
end

function label_window(img_labels::Array{<:Tuple{Tuple{Real,Real},String},1})
    new_win = Window("Labels",600,500)
    text = ""
    for (idx, i) in enumerate(img_labels)
        text *= string.(idx) * format_label(i) * "\n"
    end
    txt = get_label("Label")
    push!(new_win,txt)
    set_gtk_property!(txt, :label, text)
    return new_win
end

function close_window(closebtn::GtkToolButtonLeaf, win::GtkWindow)
    signal_connect(closebtn, :clicked) do widget
       destroy(win)
    end
end

function draw_canvas(c::GtkReactive.Canvas{UserUnit}, zr::Signal,surf::Cairo.CairoSurfaceImage)
   draw(c, zr) do widget, r
        ctx = getgc(widget)
        set_coordinates(ctx, r)
        rectangle(ctx, BoundingBox(r.currentview))
        set_source(ctx, surf)
        fill(ctx)
        end
end

function map_mouse_motion(img::AstroImage, indx::Int, c::GtkReactive.Canvas{UserUnit}, g::GtkGrid, header::WCSTransform)
    map(c.mouse.motion) do btn
        xu, yu = btn.position.x.val, btn.position.y.val
        world_coord = pix_to_world(header, [xu, yu])
        GAccessor.text(g[1,3] ,"Pixel Coordinate X = $(round(xu, digits = 2))")
        GAccessor.text(g[2,3] ,"Pixel Coordinate Y = $(round(yu, digits = 2))")
        GAccessor.text(g[1,4] ,"World Coordinate X = $(round(world_coord[1], digits = 2))"*"°")
        GAccessor.text(g[2,4] ,"World Coordinate Y = $(round(world_coord[2], digits = 2))"*"°")
        if xu >= 1 && yu >= 1
            GAccessor.text(g[1,5] ,"Value = $(img.data[indx][Int(round(xu)),Int(round(yu))])")
        else
            GAccessor.text(g[1,5] , "Value = ")
        end
    end
end


text: # https://arxiv.org/abs/1901.01780


text:
abstract type OrthogonalPolyBasis end

struct LegendrePolyBasis{T,A} <: OrthogonalPolyBasis
    n::T
    weights::A

    function LegendrePolyBasis(n; wType = Float64)
        weights_helper = zeros(wType, n, n)
        weights_helper[1, 1] = one(wType)
        weights_helper[2, 2] = one(wType)
        for i = 3:n
            real_index = i - 1
            two_below = ((real_index - 1) / real_index) .* weights_helper[i-2, :]
            one_below =
                ((2 * real_index - 1) / real_index) .*
                [[zero(wType)]; weights_helper[i-1, :][Base.OneTo(n - 1)]]
            weights_helper[i, :] .= one_below .- two_below
        end
        weights =
            [Tuple(weights_helper[i, j] for j in Base.OneTo(i)) for i in Base.OneTo(n)]
        return new{typeof(n),typeof(weights)}(n, weights)
    end
end

function (b::LegendrePolyBasis)(x, i::Int64)
    @unpack weights = b
    @evalpoly(x, weights[i+1]...)
end

struct HermitePolyBasis{T,A} <: OrthogonalPolyBasis
    n::T
    weights::A

    function HermitePolyBasis(n; wType = Int64)
        weights_helper = zeros(wType, n, n)
        weights_helper[1, 1] = one(wType)
        weights_helper[2, 2] = 2
        for i = 3:n
            real_index = i - 1
            weights_helper[i, 1] = -weights_helper[i-1, 2]
            for j = 2:i-1
                weights_helper[i, j] =
                    2 * weights_helper[i-1, j-1] - j * weights_helper[i-1, j+1]
            end
            weights_helper[i, i] = 2 * weights_helper[i-1, i-1]
        end
        weights =
            [Tuple(weights_helper[i, j] for j in Base.OneTo(i)) for i in Base.OneTo(n)]
        return new{typeof(n),typeof(weights)}(n, weights)
    end
end

function (b::HermitePolyBasis)(x, i::Int64)
    @unpack weights = b
    @evalpoly(x, weights[i+1]...)
end

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
# This file was generated, do not modify it. # hide
is_active(col) = !occursin(r"source|generation", string(col))
active_cols = [col for col in names(data) if is_active(col)]
select!(data, active_cols);


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP

"""
    generate_initial_conditions!(instance, optimizer)

Generates feasible initial conditions for the given instance, by constructing
and solving a single-period mixed-integer optimization problem, using the given
optimizer. The instance is modified in-place.
"""
function generate_initial_conditions!(
    instance::UnitCommitmentInstance,
    optimizer,
)::Nothing
    G = instance.units
    B = instance.buses
    t = 1
    mip = JuMP.Model(optimizer)

    # Decision variables
    @variable(mip, x[G], Bin)
    @variable(mip, p[G] >= 0)

    # Constraint: Minimum power
    @constraint(mip, min_power[g in G], p[g] >= g.min_power[t] * x[g])

    # Constraint: Maximum power
    @constraint(mip, max_power[g in G], p[g] <= g.max_power[t] * x[g])

    # Constraint: Production equals demand
    @constraint(
        mip,
        power_balance,
        sum(b.load[t] for b in B) == sum(p[g] for g in G)
    )

    # Constraint: Must run
    for g in G
        if g.must_run[t]
            @constraint(mip, x[g] == 1)
        end
    end

    # Objective function
    function cost_slope(g)
        mw = g.min_power[t]
        c = g.min_power_cost[t]
        for k in g.cost_segments
            mw += k.mw[t]
            c += k.mw[t] * k.cost[t]
        end
        if mw < 1e-3
            return 0.0
        else
            return c / mw
        end
    end
    @objective(mip, Min, sum(p[g] * cost_slope(g) for g in G))

    JuMP.optimize!(mip)

    for g in G
        if JuMP.value(x[g]) > 0
            g.initial_power = JuMP.value(p[g])
            g.initial_status = 24
        else
            g.initial_power = 0
            g.initial_status = -24
        end
    end
    return
end


text:
########################
# helper functions
########################

argerr(s::AbstractString) = throw(ArgumentError(s))
daterr(s::AbstractString) = throw(DataError(s))
parserr(s::AbstractString) = throw(Meta.ParseError(s))

"""
    listgroups()

Return a list of available groups (array of `Symbol`).
"""
function listgroups()
    groups = Symbol[]
    append!(groups, sort!(collect(keys(SUBSETS))))
    append!(groups, sort!(collect(keys(MATRIXCLASS))))
    append!(groups, sort!(collect(keys(usermatrixclass))))
    groups
end

#######################
# matrix group
#######################

# write one property association
function propline(io::IO, propname, matnames)
    write(io, repr(propname))
    write(io, " => [")
    for str in matnames
        write(io, repr(str))
        write(io, ", ")
    end
    write(io, "],\n")
end

# add, remove, or replace complete user group
function modgroup(prop::Symbol, mats::Union{Nothing,Vector{<:AbstractString}})
    prop in keys(MATRIXCLASS) && daterr("$prop can not be modified.")

    user = abspath(MY_DEPOT_DIR, "group.jl")
    s = read(user, String)          # read complete file into s
    rg = Regex(repr(prop) * r"\W*=>\W*(\[.*\]\W*,\W*\n)".pattern)
    ppos = findfirst(rg, s)         # locate the prop in user.jl to remove.
    if ppos !== nothing
        start_char = first(ppos) - 1    # the start of the line
        end_char = last(ppos)           # the end of the line
    else
        ppos = findnext(r"\);", s, 1)
        start_char = ppos !== nothing ? first(ppos) - 1 : length(s)
        end_char = start_char
    end
    if mats !== nothing
        mats = sort(mats)
    end
    open(user, "w") do io
        write(io, s[1:start_char])
        if mats !== nothing
            propline(io, prop, mats)
        end
        write(io, s[end_char+1:end])
    end
    if mats !== nothing
        usermatrixclass[prop] = mats
    else
        delete!(usermatrixclass, prop)
    end
    return nothing
end

"add a group to Matrix Depot"
macro addgroup(ex)
    name = Symbol(ex.args[1])
    esc(modgroup(name, eval(ex.args[2])))
end

"add or replace group in Matrix Depot"
macro modifygroup(ex)
    name = Symbol(ex.args[1])
    esc(modgroup(name, eval(ex.args[2])))
end

"remove an added group from Matrix Depot"
macro rmgroup(ex)
    name = Symbol(ex)
    esc(modgroup(name, nothing))
end

################################
# user defined matrix generators
################################

abstract type MatrixGenerator end
abstract type FunctionName <: MatrixGenerator end
abstract type Group <: MatrixGenerator end

function include_generator(::Type{FunctionName}, fn::AbstractString, f::Function)
    (haskey(MATRIXDICT, fn) ? MATRIXDICT : USERMATRIXDICT)[fn] = f
end

function addtogroup(dir::Dict, groupname::Symbol, f::Function)
    if groupname in keys(dir)
        fn = fname(f)
        gr = dir[groupname]
        fn in gr || push!(gr, fn)
        true
    else
        false
    end
end
function include_generator(::Type{Group}, groupname::Symbol, f::Function)
    addtogroup(MATRIXCLASS, groupname, f) ||
    addtogroup(usermatrixclass, groupname, f) ||
    argerr("$(groupname) is not a group in MatrixDepot, use
              @addgroup to add this group")
end

"return the name of the function `f` as a string."
function fname(f::Function)
    for (key, value) in MATRIXDICT
        value == f && return key
    end
    for (key, value) in USERMATRIXDICT
        value == f && return key
    end
    "unknown-function"
end

"""
listdir([db,] p::AbstractString)

list directories and the number of matrices contained in them.
get an overview of the count of names down directories.
return a list with summary information for directories in matrix name space.
The input argument is split into 2 patterns by the first double slash `"//"`.
The whole string (with multiple slashes reduced to single slashes) determines
a subset of all matrix names. They are then grouped by the first pattern and
for each different group value the number of names in the subset is counted.
A final `/` is replaced by `"//**"`.

E.g.
+ `listdir("/*")`     - count names without a `/`.
+ `listdir("/k*")`    - count names without `/` starting with `k*`.
+ `listdir("*//*")`   - count names with one directory part (sp-collection)
+ `listdir("*/*//*")` - count names with two directory parts (mm-collection)
+ `listdir("*//*/*")` - count names with two directory parts (mm-collection)
+ `listdir("Har*//*/*")` - restrict to directories starting with "Har"
+ `listdir("Har*/*//*")` - all subdirectoreis of the previous
"""
listdir(p::AbstractString, xp::Pattern=()) = listdir(MATRIX_DB, p, xp)
function listdir(db::MatrixDatabase, p::AbstractString, xp::Pattern)
    check_symbols(xp)
    r = findfirst(r"/+", p)
    if r !== nothing && first(r) == 1
        p = p[last(r)+1:end]
        depth = 0
    else
        m = match(r"^(([^/]+/)+)(/|$)", p)
        depth = m !== nothing ? count(x == '/' for x in m.captures[1]) : -1
    end
    p = replace(p, r"//+" => '/')
    endswith(p, '/') && ( p = string(p, "**") )
    r = shell_to_regex(p, false)
    if depth >= 0
        length(p) == 1 && ( p = ".*/" )
        listdir(db, r & xp, depth)
    else
        argerr("pattern '$p' needs '//' in the middle or '/' at start or end")
    end
end
function listdir(db::MatrixDatabase, r::Pattern, depth::Int)
    result = Dict{AbstractString, Int}()
    f(x, n) = string(join(x[1:n], '/'), "/*" ^ max(length(x) - n, 0))
    for name in mdlist(r)
        li = split(name, '/')
        if length(li) >= depth
            key = f(li, depth)
            result[key] = get!(result, key, 0) + 1
        end
    end
    sort!([string(k, " - (", v, ")") for (k,v) in result])
end

"""
    listdata([db,] p::Pattern)
Return an array of `MatrixData` objects according to matched patterns.
"""
listdata(p::Pattern) = listdata(MATRIX_DB, p)
listdata(db::MatrixDatabase, p::Pattern) = mdata.(mdlist(db, p))

"""
    mdlist([db,] p::Pattern)
return a vector of full matrix names where name or alias match given pattern.
`p` can be one of the following:
+ a plain string (without characters `*` and `?`) which must match exactly
+ a string containing `*` and `?` acting like a shell path pattern
+ a regular expression
+ an integer matching equivalent to the alias string `"#\$p"`
+ a range of integers
+ a group name expressed as a symbol e.g. `:local`, `:all`, `:illcond`, `posdef`
+ the name of a predicate function `f(::MatrixData)::Bool`, e.g. `symmetric`, ...
+ a vector of patterns meaning the union
+ a tuple of patterns meaning the intersection
"""
mdlist(p::Pattern) = mdlist(MATRIX_DB, p)
is_all(res::Vector) = length(res) == 1 && res[1] == ""

function mdlist(db::MatrixDatabase, p::Pattern)
    res = list!(db, [""], p)
    is_all(res) ? list_all(db) : res
end
function list!(db::MatrixDatabase, res::Vector{String}, r::Regex)
    isempty(res) && return res
    if is_all(res)
        empty!(res)
        for name in keys(db.data)
            if match(r, name) !== nothing
                push!(res, name)
            end
        end
        sort!(res)
    else
        for i = 1:length(res)
            name = res[i]
            if match(r, name) === nothing
                res[i] = ""
            end
        end
        for i = length(res):-1:1
            if res[i] == ""
                deleteat!(res, i)
            end
        end
    end
    res
end

function list!(db::MatrixDatabase, res::Vector{String}, p::Symbol)
    isempty(res) && return res
    x = if haskey(SUBSETS, p)
        SUBSETS[p](db)
    elseif haskey(MATRIXCLASS, p)
        MATRIXCLASS[p]
    elseif haskey(usermatrixclass, p)
        usermatrixclass[p]
    else
        argerr("unknown group name '$p'")
        # EMPTY_PATTERN
    end
    if is_all(res)
        empty!(res)
        append!(res, sort(x))
    else
        intersect!(res, x)
    end
end

"""
    shell_to_regex

return a regular expression if shell pattern characters `"*?]"` are contained in
string, otherwise return string.
If no `'/'` is contained in p and p is not "*", insert `"(.*/)?"` in regular expression.
"""
function shell_to_regex(p::AbstractString, retain_pure::Bool)
    regex(p) = Regex(string('^', p, '$'))
    p = replace(p, r"//+" => '/')
    # p = p == "*" || '/' in p ? p : string("(**/)\x03", p)
    if occursin(r"[*?.]", p)
        p = replace(p, "**" => "\x01\x02")
        p = replace(p, '*' => "[^/]*")
        p = replace(p, '?' => "[^/]")
        p = replace(p, '.' => "[.]")
        p = replace(p, '\x01' => '.')
        p = replace(p, '\x02' => '*')
        p = replace(p, '\x03' => '?')
        regex(p)
    else
        retain_pure ? p : regex(p)
    end
end

function singlist!(db::MatrixDatabase, res::Vector{String}, p::AbstractString)
    if is_all(res)
        if haskey(db.data, p)
            res[1] = p
            res
        else
            empty!(res)
        end
    else
        if p in res
            empty!(res)
            push!(res, p)
        else
            empty!(res)
        end
    end
end

function list!(db::MatrixDatabase, res::Vector{String}, p::AbstractString)
    isempty(res) && return res
    r = shell_to_regex(p, true)
    r isa Regex ? list!(db, res, r) : singlist!(db, res, r)
end

list!(db::MatrixDatabase, res::Vector{String}, p::Alias) = list!(db, res, aliasresolve(db, p))

# If res is symbolically [""], populate it with the set of all available names
function resall!(db::MatrixDatabase, res::Vector{String})
    if is_all(res)
        x = list_all(db)
        resize!(res, length(x))
        copyto!(res, x)
    end
    res
end

function list!(db::MatrixDatabase, res::Vector{String}, p::Not)
    isempty(res) && return res
    cres = copy(res)
    resall!(db, res)
    setdiff!(res, list!(db, cres, p.pattern))
end

# logical OR
function list!(db::MatrixDatabase, res::Vector{String}, r::AbstractVector)
    isempty(r) && return empty!(res)
    check_symbols(r)
    isempty(res) && return res
    length(r) == 1 && return list!(db, res, r[1])
    cres = copy(res)
    list!(db, res, r[1])
    for y in r[2:end]
        union!(res, list!(db, copy(cres), y))
    end
    sort!(res)
end

list!(db::MatrixDatabase, res::Vector{String}, ::Tuple{}) = res

# logical AND
function list!(db::MatrixDatabase, res::Vector{String}, r::Tuple)
    isempty(res) && return res
    check_symbols(r)
    for y in sort_by_type(r)
        list!(db, res, y)
        isempty(res) && break
    end
    res
end

function list!(db::MatrixDatabase, res::Vector{String}, pred::Function)
    isempty(res) && return res
    resall!(db, res)
    filter!(k -> pred(db.data[k]), res)
end

# return a vector with re-arranged the contents of itr
# the elements of its must be of type Pattern
function sort_by_type(itr)
    f(p::AbstractString) = true
    f(p::Symbol) = true
    f(p::AbstractVector) = all(f.(p))
    f(p::Tuple) = all(f.(p))
    f(p::Pattern) = false
    vall = collect(itr)
    vend = filter(!f, vall)
    filter!(f, vall)
    append!(vall, vend)
    vall
end

## internal list special cases
list_all(db::MatrixDatabase) = sort!(String.(collect(keys(db.data))))
list_local(db::MatrixDatabase) = union(collect(keys(MATRIXDICT)), keys(USERMATRIXDICT))
list_builtin(db::MatrixDatabase) = collect(keys(MATRIXDICT))
list_user(db::MatrixDatabase) = collect(keys(USERMATRIXDICT))

const SUBSETS = Dict(
                     :local => list_local,
                     :builtin => list_builtin,
                     :user => list_user,
                     :all => list_all,
)

function verify_loaded(db::MatrixDatabase, data::RemoteMatrixData)
    if isempty(data.metadata) || !isfile(matrixfile(data))
        loadmatrix(data)
    end
    data
end
verify_loaded(db::MatrixDatabase, data::MatrixData) = data

function verify_loadinfo(data::RemoteMatrixData)
    file = matrixfile(data)
    if !isfile(file)
        file = matrixinfofile(data)
        if !isfile(file)
            loadinfo(data)
        end
    end
    file
end

mdatav(db::MatrixDatabase, p::Pattern) = verify_loaded(db, mdata(db, p))

"""
    load([db,] pattern)

Load data from remote repository for all problems matching pattern.

Return the number of successfully loaded matrices.
"""
load(p::Pattern) = load(MATRIX_DB, p)
load(db::MatrixDatabase, p::Pattern) = _load(db, loadmatrix, p)

"""
    loadinfo([db,] pattern)
"""
loadinfo(p::Pattern) = loadinfo(MATRIX_DB, p)
loadinfo(db::MatrixDatabase, p::Pattern) = _load(db, loadinfo, p)

function _load(db::MatrixDatabase, loadfunc::Function, p::Pattern)
    check_symbols(p)
    n = 0
    for name in mdlist(p)
        try
            n += loadfunc(db.data[name])
        catch ex
            ex isa InterruptException && rethrow()
            @warn "could not load $name: $ex"
        end
    end
    n
end

"""
    mdopen([db,] pattern)
    mdopen(f, [db,] pattern)

Return `MatrixDescriptor` object, which can be used with data access functions.

Make sure that data files are loaded.
Keeps a cache of already delivered matrices and metadata.
If the pattern has not a unique resolution, an error is thrown.
"""
mdopen(p::Pattern, args...) = mdopen(MATRIX_DB, p, args...)
function mdopen(db::MatrixDatabase, p::Pattern, args...)
    _mdopen(mdatav(db, p), args...)
end

mdopen(f::Function, p::Pattern, args...) = mdopen(f, MATRIX_DB, p, args...)
function mdopen(f::Function, db::MatrixDatabase, p::Pattern, args...)
    data = _mdopen(mdatav(db, p), args...)
    f(data)
end

"""
    mdata(db, pattern)

Return unique `MatrixData` object according to pattern.
"""
mdata(p::Pattern) = mdata(MATRIX_DB, p)
function mdata(db::MatrixDatabase, p::Pattern)
    check_symbols(p)
    li = mdlist(db, p)
    length(li) == 0 && daterr("no matrix according to $p found")
    length(li) > 1  && daterr("pattern not unique: $p -> $li")
    db.data[li[1]]
end

"""
    metadata([db, ], Union{MatrixDescriptor,MatrixData})

Return copy of list of metadata names.
"""
metadata(mdesc::MatrixDescriptor) = metadata(mdesc.data)
metadata(data::RemoteMatrixData) = copy(data.metadata)
metadata(data::MatrixData) = String[]

_mdopen(data::RemoteMatrixData)= MatrixDescriptor(data)
function _mdopen(data::GeneratedMatrixData, args...)
    md = MatrixDescriptor(data, args...)
    md.A # trigger generation of data and fill cache
    md
end

###
# convenience API
###

"""
    matrixdepot(p::Pattern, args...)

Return matrix according to pattern or local matrix according to name and arguments.

If not loaded, load remote matrix first.
`p` must be a unique pattern (match only one name). The presence of arguments makes
sense only if the pattern matches the name of a generated (=local) matrix.

Only the matrix part is delivered, also in the local cases, where the underlying
function returns a structure containing matrix and vectors.
Use `md = mdopen; md.A, md.b ...`
to access those objects.
"""
matrixdepot(p::Pattern, args...) = matrixdepot(MATRIX_DB, p, args...)
function matrixdepot(db::MatrixDatabase, p::Pattern, args...)
    mdopen(db, p, args...) do md
        md.A
    end
end


text:
# Floating point data format invariants

const SIGN_BIT = 0x80000000
const SIGN_BIT_64 = UInt64(SIGN_BIT)

const bmask16 = 0xFFFF
const bmask32 = 0xFFFFFFFF

const UNO = one(UInt32)
const UNO64 = one(UInt64)

# VAX floating point data formats (see VAX Architecture Reference Manual)

const VAX_F_SIGN_BIT      = SIGN_BIT
const VAX_F_EXPONENT_MASK = Int32(0x7F800000)
const VAX_F_EXPONENT_SIZE = UInt32(8)
const VAX_F_EXPONENT_BIAS = Int32(128)
const VAX_F_MANTISSA_MASK = UInt32(0x007FFFFF)
const VAX_F_MANTISSA_SIZE = UInt32(23)
const VAX_F_HIDDEN_BIT    = UInt32( UNO << VAX_F_MANTISSA_SIZE )

const VAX_D_EXPONENT_MASK = Int64(VAX_F_EXPONENT_MASK)
const VAX_D_EXPONENT_SIZE = UInt64(VAX_F_EXPONENT_SIZE)
const VAX_D_EXPONENT_BIAS = Int64(VAX_F_EXPONENT_BIAS)
const VAX_D_MANTISSA_MASK = UInt64(VAX_F_MANTISSA_MASK)
const VAX_D_MANTISSA_SIZE = UInt64(VAX_F_MANTISSA_SIZE)
const VAX_D_HIDDEN_BIT    = UInt64(VAX_F_HIDDEN_BIT)

const VAX_G_EXPONENT_MASK = Int64(0x7FF00000)
const VAX_G_EXPONENT_SIZE = UInt64(11)
const VAX_G_EXPONENT_BIAS = Int64(1024)
const VAX_G_MANTISSA_MASK = UInt64(0x000FFFFF)
const VAX_G_MANTISSA_SIZE = UInt64(20)
const VAX_G_HIDDEN_BIT    = UInt64( UNO << VAX_G_MANTISSA_SIZE )

# IEEE floating point data formats (see Alpha Architecture Reference Manual)

const IEEE_S_SIGN_BIT      = SIGN_BIT
const IEEE_S_EXPONENT_MASK = Int32(0x7F800000)
const IEEE_S_EXPONENT_SIZE = UInt32(8)
const IEEE_S_EXPONENT_BIAS = Int32(127)
const IEEE_S_MANTISSA_MASK = UInt32(0x007FFFFF)
const IEEE_S_MANTISSA_SIZE = UInt32(23)
const IEEE_S_HIDDEN_BIT    = UInt32( UNO << IEEE_S_MANTISSA_SIZE )

const IEEE_T_EXPONENT_MASK = Int64(0x7FF00000)
const IEEE_T_EXPONENT_SIZE = UInt64(11)
const IEEE_T_EXPONENT_BIAS = Int64(1023)
const IEEE_T_MANTISSA_MASK = UInt64(0x000FFFFF)
const IEEE_T_MANTISSA_SIZE = UInt64(20)
const IEEE_T_HIDDEN_BIT    = UInt64( 1 << IEEE_T_MANTISSA_SIZE )


text:
#Adding Julia Scripts to Jupyter on Windows

Pkg.add("IJulia")


text:
#Reference:
#Batoz J L, Tahar M B. Evaluation of a new quadrilateral thin plate bending element [J]. International Journal for Numerical Methods in Engineering, 1982, 18(11): 1655―1677.
function K_DKQ(elm::Quad)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc.-center)
    x₂,y₂,z₂=T*(elm.node2.loc.-center)
    x₃,y₃,z₃=T*(elm.node3.loc.-center)
    x₄,y₄,z₄=T*(elm.node4.loc.-center)
    K=Matrix{Float64}(undef,12,12)
    J=Matrix{Float64}(undef,2,2)
    j=Matrix{Float64}(undef,2,2)

    function BtDB(r)
        ξ,η=r[1],r[2]
        x=[x₁,x₂,x₃,x₄]
        y=[y₁,y₂,y₃,y₄]

        xᵢⱼ=[x[i]-x[j] for (i,j) in zip([1,2,3,4],[2,3,4,1])]
        yᵢⱼ=[y[i]-y[j] for (i,j) in zip([1,2,3,4],[2,3,4,1])]
        lᵢⱼ²=xᵢⱼ.^2 .+ yᵢⱼ.^2

        a₅,a₆,a₇,a₈=-xᵢⱼ ./ lᵢⱼ²
        b₅,b₆,b₇,b₈=0.75*xᵢⱼ.*yᵢⱼ ./ lᵢⱼ²
        c₅,c₆,c₇,c₈=(0.25*xᵢⱼ.^2 .- 0.5*yᵢⱼ.^2) ./ lᵢⱼ²
        d₅,d₆,d₇,d₈=-yᵢⱼ ./ lᵢⱼ²
        e₅,e₆,e₇,e₈=(-0.25*xᵢⱼ.^2 .+ 0.5*yᵢⱼ.^2) ./ lᵢⱼ²

        dHˣdξ = [1.5*(0.5*a₈*(1 - η^2) - 1.0*ξ*(1 - η)*a₅), -1.0*ξ*(1 - η)*b₅, 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 - η)*c₅,1.5*(0.5*a₆*(1 - η^2) + 1.0*ξ*(1 - η)*a₅), 0.5*b₆*(1 - η^2), -0.5*c₆*(1 - η^2) + 1.0*ξ*(1 - η)*c₅ - 0.25*(1 + η), 1.5*(-0.5*a₆*(1 - η^2) - 1.0*ξ*(1 + η)*a₇), -1.0*ξ*(1 + η)*b₇, 0.25 + 0.25*η - 0.5*c₆*(1 - η^2) + 1.0*ξ*(1 + η)*c₇, 1.5*(-0.5*a₈*(1 - η^2) + 1.0*ξ*(1 + η)*a₇), -0.5*b₈*(1 - η^2), 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 + η)*c₇]
        dHˣdη = [1.5*(0.5*a₈*(1 - η^2) - 1.0*ξ*(1 - η)*a₅), -1.0*ξ*(1 - η)*b₅, 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 - η)*c₅,1.5*(0.5*a₆*(1 - η^2) + 1.0*ξ*(1 - η)*a₅), 0.5*b₆*(1 - η^2), -0.5*c₆*(1 - η^2) + 1.0*ξ*(1 - η)*c₅ - 0.25*(1 + η), 1.5*(-0.5*a₆*(1 - η^2) - 1.0*ξ*(1 + η)*a₇), -1.0*ξ*(1 + η)*b₇, 0.25 + 0.25*η - 0.5*c₆*(1 - η^2) + 1.0*ξ*(1 + η)*c₇, 1.5*(-0.5*a₈*(1 - η^2) + 1.0*ξ*(1 + η)*a₇), -0.5*b₈*(1 - η^2), 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 + η)*c₇]
        dHʸdξ = [1.5*(0.5*a₈*(1 - η^2) - 1.0*ξ*(1 - η)*a₅), -1.0*ξ*(1 - η)*b₅, 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 - η)*c₅,1.5*(0.5*a₆*(1 - η^2) + 1.0*ξ*(1 - η)*a₅), 0.5*b₆*(1 - η^2), -0.5*c₆*(1 - η^2) + 1.0*ξ*(1 - η)*c₅ - 0.25*(1 + η), 1.5*(-0.5*a₆*(1 - η^2) - 1.0*ξ*(1 + η)*a₇), -1.0*ξ*(1 + η)*b₇, 0.25 + 0.25*η - 0.5*c₆*(1 - η^2) + 1.0*ξ*(1 + η)*c₇, 1.5*(-0.5*a₈*(1 - η^2) + 1.0*ξ*(1 + η)*a₇), -0.5*b₈*(1 - η^2), 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 + η)*c₇]
        dHʸdη = [1.5*(0.5*a₈*(1 - η^2) - 1.0*ξ*(1 - η)*a₅), -1.0*ξ*(1 - η)*b₅, 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 - η)*c₅,1.5*(0.5*a₆*(1 - η^2) + 1.0*ξ*(1 - η)*a₅), 0.5*b₆*(1 - η^2), -0.5*c₆*(1 - η^2) + 1.0*ξ*(1 - η)*c₅ - 0.25*(1 + η), 1.5*(-0.5*a₆*(1 - η^2) - 1.0*ξ*(1 + η)*a₇), -1.0*ξ*(1 + η)*b₇, 0.25 + 0.25*η - 0.5*c₆*(1 - η^2) + 1.0*ξ*(1 + η)*c₇, 1.5*(-0.5*a₈*(1 - η^2) + 1.0*ξ*(1 + η)*a₇), -0.5*b₈*(1 - η^2), 0.25 + 0.25*η + 0.5*c₈*(1 - η^2) + 1.0*ξ*(1 + η)*c₇]

        x₁₂,x₂₃,x₃₄,x₄₁=xᵢⱼ
        y₁₂,y₂₃,y₃₄,y₄₁=xᵢⱼ
        x₂₁=-x₁₂
        x₃₂=-x₂₃
        y₂₁=-y₁₂
        y₃₂=-y₂₃
        # J[1,1]=0.25*(x₂₁+x₃₄+η*(x₁₂+x₃₄))
        # J[2,1]=0.25*(x₃₂+x₄₁+ξ*(x₁₂+x₃₄))
        # J[1,2]=0.25*(y₂₁+y₃₄+η*(y₁₂+y₃₄))
        # J[2,2]=0.25*(y₃₂+y₄₁+ξ*(y₁₂+y₃₄))
        J[1,1]=-0.25*(1 - η)*x₁ + 0.25*(1 - η)*x₂ + 0.25*(1 + η)*x₃ - 0.25*(1 + η)*x₄
        J[1,2]=-0.25*(1 - η)*y₁ + 0.25*(1 - η)*y₂ + 0.25*(1 + η)*y₃ - 0.25*(1 + η)*y₄
        J[2,1]=-0.25*(1 - ξ)*x₁ + 0.25*(1 - ξ)*x₄ - 0.25*(1 + ξ)*x₂ + 0.25*(1 + ξ)*x₃
        J[2,2]=-0.25*(1 - ξ)*y₁ + 0.25*(1 - ξ)*y₄ - 0.25*(1 + ξ)*y₂ + 0.25*(1 + ξ)*y₃

        detJ=det(J)
        j[1,1]=J[1,1]/detJ
        j[2,1]=-J[2,1]/detJ
        j[1,2]=-J[1,2]/detJ
        j[2,2]=J[2,2]/detJ
        @show detJ
        B=[j[1,1]*dHˣdξ+j[1,2]*dHˣdη j[2,1]*dHʸdξ+j[2,2]*dHʸdη j[1,1]*dHʸdξ+j[1,2]*dHʸdη+j[2,1]*dHˣdξ+j[2,2]*dHˣdξ]'
        D₀=E₀*t^3/12/(1-ν₀^2)
        D=D₀*[1 ν₀ 0;
              ν₀ 1 0;
              0  0 (1-ν₀)/2]
        return transpose(B)*D*B*detJ
    end
    K=hcubature(BtDB,[-1,-1],[1,1],initdiv=2)[1]
    #12x12 to 24x24
    I=1:12
    J=[3,4,5,9,10,11,15,16,17,21,22,23]
    L=sparse(I,J,1.,12,24)
    return L'*K*L
end


text:
using FermiCG
using Printf
using NPZ
using JSON
using Pkg
using LinearAlgebra

testdir = joinpath(dirname(pathof(FermiCG)), "..", "test")
println(testdir)

ket_a = ConfigString(no=4, ne=4)

filepath = joinpath(testdir, "data_h4/ints_0b.npy")
ints_0b = npzread(filepath);

filepath = joinpath(testdir, "data_h4/ints_1b.npy")
ints_1b = npzread(filepath);

filepath = joinpath(testdir, "data_h4/ints_2b.npy")
ints_2b = npzread(filepath);

filepath = joinpath(testdir, "data_h4/problem.json")
data = JSON.parsefile(filepath)

println(data)
spin = data["spin"]
n_elec = data["n_elec"]
n_elec_a = round(Int,(spin + n_elec)/2)
n_elec_b = n_elec - n_elec_a

ham 	= InCoreInts(ints_0b, ints_1b, ints_2b)
#display(rand(Float64, (4,4)) )

#using Plots
#N = size(ham.h2)[1]
#println(size(ham.h2))
#heatmap(ham.h1, color = :greys)
#gui()
#print("what?")

# problem = ElectronicProblem(no=size(ints_1b,1), na=n_elec_a, nb=n_elec_b)
# display(problem)

atoms = []
push!(atoms,Atom(1,"H",[0,0,0]))
push!(atoms,Atom(2,"H",[0,0,1]))
push!(atoms,Atom(3,"H",[0,0,2]))
push!(atoms,Atom(4,"H",[0,0,3]))

mol     = Molecule(0,1,atoms)
display(mol)

mf = FermiCG.pyscf_do_scf(mol,"6-31g")

FermiCG.pyscf_write_molden(mol,"sto-3g",mf.mo_coeff)
FermiCG.pyscf_write_molden(mf,filename="2.molden")
FermiCG.pyscf_fci(ham,2,2)

n_orbs = size(mf.mo_coeff)[2]
orb_indices = [i for i = 1:n_orbs]
emb_indices = [i for i = 1:2]
display(emb_indices)
emb_orbs = mf.mo_coeff[:,emb_indices];
emb_density = 2*emb_orbs * emb_orbs';
display(emb_density)
PS = emb_density*mf.get_ovlp();
n = tr(PS)

#ints = FermiCG.pyscf_build_ints(mf.mol,mf.mo_coeff[:,[3,4]], emb_density);
#FermiCG.pyscf_fci(ints,0,0)

ints = FermiCG.pyscf_build_ints(mf.mol,mf.mo_coeff);
FermiCG.pyscf_fci(ints,2,2)

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
@testset "ldg" begin
    ir = sprint(io->CUDA.code_llvm(io, CUDA.pointerref_ldg, Tuple{Core.LLVMPtr{Int,AS.Global},Int,Val{1}}))
    @test occursin("@llvm.nvvm.ldg", ir)
end

capability(device()) >= v"3.2" && @testset "unsafe_cached_load" begin

@testset for T in (Int8, UInt16, Int32, UInt32, Int64, UInt64, Int128, Float32, Float64)
    d_a = CuArray(ones(T))
    d_b = CuArray(zeros(T))
    @test Array(d_a) != Array(d_b)

    ptr_a = reinterpret(Core.LLVMPtr{T,AS.Global}, pointer(d_a))
    ptr_b = reinterpret(Core.LLVMPtr{T,AS.Global}, pointer(d_b))

    let ptr_a=ptr_a, ptr_b=ptr_b #JuliaLang/julia#15276
        @on_device unsafe_store!(ptr_b, unsafe_cached_load(ptr_a))
    end

    @test Array(d_a) == Array(d_b)
end

@testset "Const" begin
    function kernel(a, b, i)
        @inbounds b[i] = Base.Experimental.Const(a)[i]
        return
    end

    buf = IOBuffer()

    a = CuArray([0])
    b = CuArray([0])
    @device_code_ptx io=buf @cuda kernel(a, b, 1)
    @test Array(a) == Array(b)

    asm = String(take!(copy(buf)))
    @test occursin("ld.global.nc", asm)
end

end


text:
@testset "Luxor trading logic"
  @testset "Market state"
    fm = TradingLogic.luxormktstate
    @test fm(120.0, 50.0) == :trendup
    @test fm(20.0, 50.0) == :trenddown
    @test fm(20.0, 20.0) == :undefined
    @test fm(NaN, 20.0) == :undefined
    @test fm(20.0, NaN) == :undefined
    @test fm(NaN, NaN) == :undefined
  end
  @testset "Target position"
    tq = 100
    ft(mkt, posnow) = TradingLogic.luxorposlogic(
      mkt, 155.0, 150.0, 10.0, tq, [posnow])
    stlim = Array(Float64, 0)

    # enter long
    @test ft(:trendup, 0) == (tq, [165.0, 155.0])
    # enter long partial fill
    p = round(Int64, tq/2)
    @test ft(:trendup, p) == (tq - p, [165.0, 155.0])

    # enter short
    @test ft(:trenddown, 0) == (-tq, [140.0, 150.0])
    # enter short partial fill
    p = -round(Int64, tq/2)
    @test ft(:trenddown, p) == (-abs(tq-p), [140.0, 150.0])

    # hold position in line with the market state
    @test ft(:trendup, tq) == (0, stlim)
    @test ft(:trenddown, -tq) == (0, stlim)

    # exit long position
    @test ft(:trenddown, tq) == (-tq, stlim)
    p = round(Int64, tq/2)
    @test ft(:trenddown, p) == (-p, stlim)

    # exit short position
    @test ft(:trendup, -tq) == (tq, stlim)
    p = -round(Int64, tq/2)
    @test ft(:trendup, p) == (abs(p), stlim)

    # undefined: wait
    @test ft(:undefined, 0) == (0, stlim)
    @test ft(:undefined, -tq) == (0, stlim)
    @test ft(:undefined, tq) == (0, stlim)
  end
end

@testset "Luxor strategy backtesting"
  @testset "GBPUSD vs. quantstrat"
    mafast = 10
    maslow = 40
    targetqty = 100
    pthresh = 10.0

    ### TODO change data input
    s_ohlc = Reactive.Signal((Dates.DateTime(ohlc.timestamp[1]),
                              vec(ohlc.values[1,:])))
    ohlc_inds = Dict{Symbol,Int64}()
    ohlc_inds[:open] = 1
    ohlc_inds[:high] = 2
    ohlc_inds[:low] = 3
    ohlc_inds[:close] = 4

    # backtest at next-open price
    # quantstrat fills tracsactions at next open on enter-signal
    s_pnow = Reactive.map(s -> s[2][ohlc_inds[:open]], s_ohlc, typ=Float64)
    blotter = TradingLogic.emptyblotter()

    s_status = TradingLogic.runtrading!(
      blotter, true, s_ohlc, ohlc_inds, s_pnow, 0,
      TradingLogic.luxortarget, mafast, maslow, pthresh, targetqty)

    s_perf = TradingLogic.tradeperfcurr(s_status)

    for i = 2:length(ohlc)
      push!(s_ohlc, (Dates.DateTime(ohlc.timestamp[i]),
                     vec(ohlc.values[i,:])))
    end

    TradingLogic.printblotter(STDOUT, blotter)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    println(vt)
    println(perfm)

    ### TODO
    @pending 0 == 1
  end
end