text:
module ERFA

import Base.getindex

export
    ERFA_DPI,
    ERFA_D2PI,
    ERFA_DR2D,
    ERFA_DD2R,
    ERFA_DR2AS,
    ERFA_DAS2R,
    ERFA_DS2R,
    ERFA_TURNAS,
    ERFA_DMAS2R,
    ERFA_DTY,
    ERFA_DAYSEC,
    ERFA_DJ,
    ERFA_DJC,
    ERFA_DJM,
    ERFA_DJ00,
    ERFA_DJM0,
    ERFA_DJM00,
    ERFA_DJM77,
    ERFA_TTMTAI,
    ERFA_DAU,
    ERFA_CMPS,
    ERFA_AULT,
    ERFA_DC,
    ERFA_ELG,
    ERFA_ELB,
    ERFA_TDB0,
    ERFA_SRS,
    ERFA_WGS84,
    ERFA_GRS80,
    ERFA_WGS72

export
    eraASTROM,
    eraLDBODY,
    eraA2af,
    eraA2tf,
    eraAb,
    eraAf2a,
    eraAnp,
    eraAnpm,
    eraApcg,
    eraApcg13,
    eraApci,
    eraApci13,
    eraApco,
    eraApco13,
    eraApcs,
    eraApcs13,
    eraAper,
    eraAper13,
    eraApio,
    eraApio13,
    eraAtci13,
    eraAtciq,
    eraAtciqn,
    eraAtciqz,
    eraAtco13,
    eraAtic13,
    eraAticq,
    eraAticqn,
    eraAtio13,
    eraAtioq,
    eraAtoc13,
    eraAtoi13,
    eraAtoiq,
    eraBi00,
    eraBp00,
    eraBp06,
    eraBpn2xy,
    eraCal2jd,
    eraC2s,
    eraC2i00a,
    eraC2i00b,
    eraC2i06a,
    eraC2ibpn,
    eraC2ixy,
    eraC2ixys,
    eraC2t00a,
    eraC2t00b,
    eraC2t06a,
    eraC2tcio,
    eraC2teqx,
    eraC2tpe,
    eraC2txy,
    eraCp,
    eraCpv,
    eraCr,
    eraDat,
    eraD2dtf,
    eraD2tf,
    eraDtdb,
    eraDtf2d,
    eraEceq06,
    eraEcm06,
    eraEe00,
    eraEe00a,
    eraEe00b,
    eraEe06a,
    eraEect00,
    eraEform,
    eraEo06a,
    eraEors,
    eraEpb,
    eraEpb2jd,
    eraEpj,
    eraEpj2jd,
    eraEpv00,
    eraEqec06,
    eraEqeq94,
    eraEra00,
    eraFad03,
    eraFae03,
    eraFaf03,
    eraFaju03,
    eraFal03,
    eraFalp03,
    eraFama03,
    eraFame03,
    eraFane03,
    eraFaom03,
    eraFapa03,
    eraFasa03,
    eraFaur03,
    eraFave03,
    eraFk52h,
    eraFk5hip,
    eraFk5hz,
    eraFw2m,
    eraFw2xy,
    eraG2icrs,
    eraGc2gd,
    eraGc2gde,
    eraGd2gc,
    eraGd2gce,
    eraGmst00,
    eraGmst06,
    eraGmst82,
    eraGst00a,
    eraGst00b,
    eraGst06,
    eraGst06a,
    eraGst94,
    eraH2fk5,
    eraHfk5z,
    eraIcrs2g,
    eraIr,
    eraJd2cal,
    eraJdcalf,
    eraLd,
    eraLdn,
    eraLdsun,
    eraLteceq,
    eraLtecm,
    eraLteqec,
    eraLtp,
    eraLtpb,
    eraLtpecl,
    eraLtpequ,
    eraNum00a,
    eraNum00b,
    eraNum06a,
    eraNumat,
    eraNut00a,
    eraNut00b,
    eraNut06a,
    eraNut80,
    eraNutm80,
    eraObl06,
    eraObl80,
    eraP06e,
    eraP2pv,
    eraP2s,
    eraPap,
    eraPas,
    eraPb06,
    eraPdp,
    eraPfw06,
    eraPlan94,
    eraPm,
    eraPmp,
    eraPmat00,
    eraPmat06,
    eraPmat76,
    eraPmpx,
    eraPmsafe,
    eraPn,
    eraPn00,
    eraPn00a,
    eraPn00b,
    eraPn06,
    eraPn06a,
    eraPnm00a,
    eraPnm00b,
    eraPnm06a,
    eraPnm80,
    eraPom00,
    eraPpp,
    eraPpsp,
    eraPr00,
    eraPrec76,
    eraPv2p,
    eraPv2s,
    eraPvdpv,
    eraPvm,
    eraPvmpv,
    eraPvppv,
    eraPvstar,
    eraPvtob,
    eraPvu,
    eraPvup,
    eraPvxpv,
    eraPxp,
    eraRefco,
    eraRm2v,
    eraRv2m,
    eraRx,
    eraRy,
    eraRz,
    eraRxp,
    eraRxpv,
    eraRxr,
    eraS00,
    eraS00a,
    eraS00b,
    eraS06,
    eraS06a,
    eraS2c,
    eraS2p,
    eraS2pv,
    eraS2xpv,
    eraSp00,
    eraSepp,
    eraSeps,
    eraStarpm,
    eraStarpv,
    eraSxp,
    eraSxpv,
    eraTaitt,
    eraTaiut1,
    eraTaiutc,
    eraTcbtdb,
    eraTcgtt,
    eraTdbtcb,
    eraTdbtt,
    eraTf2a,
    eraTf2d,
    eraTr,
    eraTrxp,
    eraTrxpv,
    eraTttai,
    eraTttcg,
    eraTttdb,
    eraTtut1,
    eraUt1tai,
    eraUt1tt,
    eraUt1utc,
    eraUtctai,
    eraUtcut1,
    eraXy06,
    eraXys00a,
    eraXys00b,
    eraXys06a,
    eraZp,
    eraZpv,
    eraZr

include("../deps/deps.jl")
include("erfa_common.jl")

function getindex(A::Array_3_Cdouble, i::Integer)
    if i === 1
        return A.d1
    elseif i === 2
        return A.d2
    elseif i === 3
        return A.d3
    end
end

function getindex(A::Array_3_Array_3_Cdouble, i::Integer)
    if i === 1
        return getindex(A.d1,1)
    elseif i === 2
        return getindex(A.d1,2)
    elseif i === 3
        return getindex(A.d1,3)
    elseif i === 4
        return getindex(A.d2,1)
    elseif i === 5
        return getindex(A.d2,2)
    elseif i === 6
        return getindex(A.d2,3)
    elseif i === 7
        return getindex(A.d3,1)
    elseif i === 8
        return getindex(A.d3,2)
    elseif i === 9
        return getindex(A.d3,3)
    end
end

function eraASTROM(pmt::Float64, eb::Array{Float64}, eh::Array{Float64}, em::Float64, v::Array{Float64}, bm1::Float64, bpn::Array{Float64}, along::Float64, phi::Float64, xpl::Float64, ypl::Float64, sphi::Float64, cphi::Float64, diurab::Float64, eral::Float64, refa::Float64, refb::Float64)
    eraASTROM(pmt,
              Array_3_Cdouble(eb[1], eb[2], eb[3]),
              Array_3_Cdouble(eh[1], eh[2], eh[3]),
              em,
              Array_3_Cdouble(v[1], v[2], v[3]),
              bm1,
              Array_3_Array_3_Cdouble(Array_3_Cdouble(bpn[1],bpn[2],bpn[3]),
                                      Array_3_Cdouble(bpn[4],bpn[5],bpn[6]),
                                      Array_3_Cdouble(bpn[7],bpn[8],bpn[9])),
              along,
              phi,
              xpl,
              ypl,
              sphi,
              cphi,
              diurab,
              eral,
              refa,
              refb)
end

function eraAb(pnat::Array{Cdouble},v::Array{Cdouble},s::Cdouble,bm1::Cdouble)
    ppr = zeros(3)
    ccall((:eraAb,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,Ptr{Cdouble}),
          pnat,v,s,bm1,ppr)
    ppr
end

function eraApcg(date1::Cdouble,date2::Cdouble,ebpv::Array{Cdouble},ehp::Array{Cdouble})
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApcg,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{eraASTROM}),
          date1,date2,ebpv,ehp,&astrom)
    astrom
end

function eraApcg13(date1::Cdouble,date2::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApcg13,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM}),
          date1,date2,&astrom)
    astrom
end

function eraApci(date1::Cdouble,date2::Cdouble,ebpv::Array{Cdouble},ehp::Array{Cdouble},x::Cdouble,y::Cdouble,s::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApci,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,Cdouble,Ptr{eraASTROM}),
          date1,date2,ebpv,ehp,x,y,s,&astrom)
    astrom
end

function eraApci13(date1::Cdouble,date2::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    eo = [0.]
    ccall((:eraApci13,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble}),
          date1,date2,&astrom,eo)
    astrom, eo[1]
end

function eraApco(date1::Cdouble,date2::Cdouble,ebpv::Array{Cdouble},ehp::Array{Cdouble},x::Cdouble,y::Cdouble,s::Cdouble,theta::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,sp::Cdouble,refa::Cdouble,refb::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApco,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM}),
          date1,date2,ebpv,ehp,x,y,s,theta,elong,phi,hm,xp,yp,sp,refa,refb,&astrom)
    astrom
end

function eraApco13(utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    eo = [0.]
    i = ccall((:eraApco13,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble}),
              utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,&astrom,eo)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    astrom, eo[1]
end

function eraApcs(date1::Cdouble,date2::Cdouble,pv::Array{Cdouble},ebpv::Array{Cdouble},ehp::Array{Cdouble})
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApcs,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{eraASTROM}),
          date1,date2,pv,ebpv,ehp,&astrom)
    astrom
end

function eraApcs13(date1::Cdouble,date2::Cdouble,pv::Array{Cdouble})
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApcs13,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{eraASTROM}),
          date1,date2,pv,&astrom)
    astrom
end

function eraAper(theta::Cdouble,astrom::eraASTROM)
    ccall((:eraAper,liberfa),Void,
          (Cdouble,Ptr{eraASTROM}),
          theta,&astrom)
    astrom
end

function eraAper13(ut11::Cdouble,ut12::Cdouble,astrom::eraASTROM)
    ccall((:eraAper13,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM}),
          ut11,ut12,&astrom)
    astrom
end

function eraApio(sp::Cdouble,theta::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,refa::Cdouble,refb::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    ccall((:eraApio,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM}),
          sp,theta,elong,phi,hm,xp,yp,refa,refb,&astrom)
    astrom
end

function eraApio13(utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    astrom = eraASTROM(0.0,zeros(3),zeros(3),0.0,zeros(3),0.0,zeros((3,3)),0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)
    i = ccall((:eraApio13,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM}),
              utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,&astrom)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    astrom
end

function eraAtci13(rc::Cdouble,dc::Cdouble,pr::Cdouble,pd::Cdouble,px::Cdouble,rv::Cdouble,date1::Cdouble,date2::Cdouble)
    ri = [0.]
    di = [0.]
    eo = [0.]
    ccall((:eraAtci13,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          rc,dc,pr,pd,px,rv,date1,date2,ri,di,eo)
    ri[1], di[1], eo[1]
end

function eraAtciq(rc::Cdouble,dc::Cdouble,pr::Cdouble,pd::Cdouble,px::Cdouble,rv::Cdouble,astrom::eraASTROM)
    ri = [0.]
    di = [0.]
    ccall((:eraAtciq,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble},Ptr{Cdouble}),
          rc,dc,pr,pd,px,rv,&astrom,ri,di)
    ri[1], di[1]
end

function eraAtciqn(rc::Cdouble,dc::Cdouble,pr::Cdouble,pd::Cdouble,px::Cdouble,rv::Cdouble,astrom::eraASTROM,b::Array{eraLDBODY})
    ri = [0.]
    di = [0.]
    n = length(b)
    ccall((:eraAtciqn,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{eraASTROM},Cint,Ptr{eraLDBODY},Ptr{Cdouble},Ptr{Cdouble}),
          rc,dc,pr,pd,px,rv,&astrom,n,b,ri,di)
    ri[1], di[1]
end

function eraAtciqz(rc::Cdouble,dc::Cdouble,astrom::eraASTROM)
    ri = [0.]
    di = [0.]
    ccall((:eraAtciqz,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble},Ptr{Cdouble}),
          rc,dc,&astrom,ri,di)
    ri[1], di[1]
end

function eraAtco13(rc::Cdouble,dc::Cdouble,pr::Cdouble,pd::Cdouble,px::Cdouble,rv::Cdouble,utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    aob = [0.]
    zob = [0.]
    hob = [0.]
    dob = [0.]
    rob = [0.]
    eo = [0.]
    i = ccall((:eraAtco13,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              rc,dc,pr,pd,px,rv,utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,aob,zob,hob,dob,rob,eo)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    aob[1],zob[1],hob[1],dob[1],rob[1],eo[1]
end

function eraAtic13(ri::Cdouble,di::Cdouble,date1::Cdouble,date2::Cdouble)
    rc = [0.]
    dc = [0.]
    eo = [0.]
    ccall((:eraAtic13,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          ri,di,date1,date2,rc,dc,eo)
    rc[1],dc[1],eo[1]
end

function eraAticq(ri::Cdouble,di::Cdouble,astrom::eraASTROM)
    rc = [0.]
    dc = [0.]
    ccall((:eraAticq,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble},Ptr{Cdouble}),
          ri,di,&astrom,rc,dc)
    rc[1],dc[1]
end

function eraAticqn(ri::Cdouble,di::Cdouble,astrom::eraASTROM,b::Array{eraLDBODY})
    rc = [0.]
    dc = [0.]
    n = length(b)
    ccall((:eraAticqn,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM},Cint,Ptr{eraLDBODY},Ptr{Cdouble},Ptr{Cdouble}),
          ri,di,&astrom,n,b,rc,dc)
    rc[1],dc[1]
end

function eraAtio13(ri::Cdouble,di::Cdouble,utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    aob = [0.]
    zob = [0.]
    hob = [0.]
    dob = [0.]
    rob = [0.]
    i = ccall((:eraAtio13,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              ri,di,utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,aob,zob,hob,dob,rob)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    aob[1],zob[1],hob[1],dob[1],rob[1]
end

function eraAtioq(ri::Cdouble,di::Cdouble,astrom::eraASTROM)
    aob = [0.]
    zob = [0.]
    hob = [0.]
    dob = [0.]
    rob = [0.]
    ccall((:eraAtioq,liberfa),Void,
          (Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          ri,di,&astrom,aob,zob,hob,dob,rob)
    aob[1],zob[1],hob[1],dob[1],rob[1]
end

function eraAtoc13(typeofcoordinates::Char,ob1::Cdouble,ob2::Cdouble,utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    rc = [0.]
    dc = [0.]
    if !(typeofcoordinates in ('R', 'r', 'H', 'h', 'A', 'a'))
        typeofcoordinates = 'A'
    end
    i = ccall((:eraAtoc13,liberfa),Cint,
              (Ptr{Char},Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
              &typeofcoordinates,ob1,ob2,utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,rc,dc)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    rc[1],dc[1]
end

function eraAtoi13(typeofcoordinates::Char,ob1::Cdouble,ob2::Cdouble,utc1::Cdouble,utc2::Cdouble,dut1::Cdouble,elong::Cdouble,phi::Cdouble,hm::Cdouble,xp::Cdouble,yp::Cdouble,phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    ri = [0.]
    di = [0.]
    i = ccall((:eraAtoi13,liberfa),Cint,
              (Ptr{Char},Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
              &typeofcoordinates,ob1,ob2,utc1,utc2,dut1,elong,phi,hm,xp,yp,phpa,tk,rh,wl,ri,di)
    if i == -1
        error("unacceptable date")
    elseif i == +1
        warn("dubious year")
    end
    ri[1],di[1]
end

function eraAtoiq(typeofcoordinates::Char,ob1::Cdouble,ob2::Cdouble,astrom::eraASTROM)
    ri = [0.]
    di = [0.]
    ccall((:eraAtoiq,liberfa),
          Void,(Ptr{Char},Cdouble,Cdouble,Ptr{eraASTROM},Ptr{Cdouble},Ptr{Cdouble}),
          &typeofcoordinates,ob1,ob2,&astrom,ri,di)
    ri[1],di[1]
end

function eraBi00()
    dpsibi = [0.]
    depsbi = [0.]
    dra = [0.]
    ccall((:eraBi00,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          dpsibi,depsbi,dra)
    dpsibi[1], depsbi[1], dra[1]
end

function eraBpn2xy(rbpn::Array{Cdouble})
    x = [0.]
    y = [0.]
    ccall((:eraBpn2xy,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          rbpn,x,y)
    x[1], y[1]
end

function eraC2ibpn(date1::Cdouble,date2::Cdouble,rbpn::Array{Cdouble})
    rc2i = zeros((3,3))
    ccall((:eraC2ibpn,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          date1,date2,rbpn,rc2i)
    rc2i
end

function eraC2s(p::Array{Cdouble})
    theta = [0.]
    phi = [0.]
    ccall((:eraC2s,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          p,theta,phi)
    theta[1], phi[1]
end

function eraCal2jd(iy::Integer, imo::Integer, id::Integer)
    r1 = [0.]
    r2 = [0.]
    i = ccall((:eraCal2jd,liberfa), Cint,
              (Cint,Cint,Cint,Ptr{Float64},Ptr{Float64}),
              iy, imo, id, r1, r2)
    @assert i == 0
    r1[1], r2[1]
end

function eraCp(p::Array{Cdouble})
    c = zeros(3)
    ccall((:eraCp,liberfa),Void,(Ptr{Cdouble},Ptr{Cdouble}),
          p,c)
    c
end

function eraCpv(pv::Array{Cdouble})
    c = zeros((2,3))
    ccall((:eraCpv,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          pv,c)
    c
end

function eraCr(p::Array{Cdouble})
    r = zeros((3,3))
    ccall((:eraCr,liberfa),Void,(Ptr{Cdouble},Ptr{Cdouble}),
          p,r)
    r
end

function eraDtdb(date1::Cdouble,date2::Cdouble,ut::Cdouble,elong::Cdouble,u::Cdouble,v::Cdouble)
    ccall((:eraDtdb,liberfa),Cdouble,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble),
          date1,date2,ut,elong,u,v)
end

function eraEform(n::Integer)
    a = [0.]
    f = [0.]
    i = ccall((:eraEform,liberfa),Cint,
              (Cint,Ptr{Cdouble},Ptr{Cdouble}),
              n,a,f)
    if i == -1
        error("illegal identifier")
    end
    a[1], f[1]
end

function eraEors(rnpb::Array{Cdouble},s::Cdouble)
    ccall((:eraEors,liberfa),Cdouble,
          (Ptr{Cdouble},Cdouble),
          rnpb,s)
end

function eraJd2cal(d1::Real, d2::Real)
    iy = Int32[0]
    imo = Int32[0]
    id = Int32[0]
    fd = [0.]
    i = ccall((:eraJd2cal,liberfa), Cint,
              (Float64,Float64,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Float64}),
              d1, d2, iy, imo, id, fd)
    @assert i == 0
    iy[1], imo[1], id[1], fd[1]
end

function eraJdcalf(ndp::Integer, d1::Real, d2::Real)
    iymdf = Int32[0, 0, 0, 0]
    i = ccall((:eraJdcalf,liberfa), Cint,
              (Cint,Float64,Float64,Ptr{Cint}),
              ndp, d1, d2, iymdf)
    @assert i == 0
    iymdf[1], iymdf[2], iymdf[3], iymdf[4]
end

function eraDat(iy::Integer, im::Integer, id::Integer, fd::Real)
    d = [0.]
    i = ccall((:eraDat, liberfa), Cint,
              (Cint, Cint, Cint, Float64, Ptr{Float64}),
              iy, im, id, fd, d)
    @assert i == 0
    d[1]
end

function eraD2dtf(scale::AbstractString, ndp::Integer, d1::Real, d2::Real)
    iy = Int32[0]
    imo = Int32[0]
    id = Int32[0]
    ihmsf = Int32[0, 0, 0, 0]
    i = ccall((:eraD2dtf,liberfa), Cint,
              (Ptr{Cchar},Cint,Float64,Float64,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
              scale, ndp, d1, d2, iy, imo, id, ihmsf)
    @assert i == 0
    iy[1], imo[1], id[1], ihmsf[1], ihmsf[2], ihmsf[3], ihmsf[4]
end

function eraDtf2d(scale::AbstractString, iy::Integer, imo::Integer, id::Integer, ih::Integer, imi::Integer, sec::Real)
    r1 = [0.]
    r2 = [0.]
    i = ccall((:eraDtf2d,liberfa), Cint,
              (Ptr{Cchar},Cint,Cint,Cint,Cint,Cint,Float64,Ptr{Float64},Ptr{Float64}),
              scale, iy, imo, id, ih, imi, sec, r1, r2)
    @assert i == 0
    r1[1], r2[1]
end

function eraEpv00(date1::Float64, date2::Float64)
    pvh = zeros((2,3))
    pvb = zeros((2,3))
    i = ccall((:eraEpv00, liberfa),
              Cint,
              (Float64, Float64, Ptr{Float64}, Ptr{Float64}),
              date1, date2, pvh, pvb)
    if i == 1
        warn("date outside the range 1900-2100 AD")
    end
    pvh, pvb
end

function eraFk5hip()
    r5h = zeros((3,3))
    s5h = zeros(3)
    ccall((:eraFk5hip,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          r5h,s5h)
    r5h,s5h
end

function eraFk5hz(r5::Cdouble,d5::Cdouble,date1::Cdouble,date2::Cdouble)
    rh = [0.]
    dh = [0.]
    ccall((:eraFk5hz,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          r5,d5,date1,date2,rh,dh)
    rh[1],dh[1]
end

function eraFw2xy(gamb::Cdouble,phib::Cdouble,psi::Cdouble,eps::Cdouble)
    x =[0.]
    y = [0.]
    ccall((:eraFw2xy,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          gamb,phib,psi,eps,x,y)
    x[1], y[1]
end

function eraGc2gd(n::Integer,xyz::Array{Cdouble})
    elong = [0.]
    phi = [0.]
    height = [0.]
    i = ccall((:eraGc2gd,liberfa),Cint,
              (Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              n,xyz,elong,phi,height)
    if i == -1
        error("illegal identifier")
    elseif i == -2
        error("internal error")
    end
    elong[1],phi[1],height[1]
end

function eraGc2gde(a::Cdouble,f::Cdouble,xyz::Array{Cdouble})
    elong = [0.]
    phi = [0.]
    height = [0.]
    i = ccall((:eraGc2gde,liberfa),Cint,
              (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              a,f,xyz,elong,phi,height)
    if i == -1
        error("illegal f")
    elseif i == -2
        error("internal a")
    end
    elong[1],phi[1],height[1]
end

function eraGd2gc(n::Integer,elong::Cdouble,phi::Cdouble,height::Cdouble)
    xyz = zeros(3)
    i = ccall((:eraGd2gc,liberfa),Cint,
              (Cint,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
              n,elong,phi,height,xyz)
    if i == -1
        error("illegal identifier")
    elseif i == -2
        error("illegal case")
    end
    xyz
end

function eraGd2gce(a::Cdouble,f::Cdouble,elong::Cdouble,phi::Cdouble,height::Cdouble)
    xyz = zeros(3)
    i = ccall((:eraGd2gce,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
              a,f,elong,phi,height,xyz)
    if i == -1
        error("illegal case")
    end
    xyz
end

function eraGst06(uta::Cdouble,utb::Cdouble,tta::Cdouble,ttb::Cdouble,rnpb::Array{Cdouble})
    ccall((:eraGst06,liberfa),Cdouble,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
          uta,utb,tta,ttb,rnpb)
end

function eraHfk5z(rh::Cdouble,dh::Cdouble,date1::Cdouble,date2::Cdouble)
    r5 = [0.]
    d5 = [0.]
    dr5 = [0.]
    dd5 = [0.]
    ccall((:eraHfk5z,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          rh,dh,date1,date2,r5,d5,dr5,dd5)
    r5[1],d5[1],dr5[1],dd5[1]
end

function eraLDBODY(bm::Cdouble, dl::Cdouble, pv::Array{Float64})
    p = Array_3_Cdouble(pv[1], pv[2], pv[3])
    v = Array_3_Cdouble(pv[4], pv[5], pv[6])
    eraLDBODY(bm, dl, Array_2_Array_3_Cdouble(p, v))
end

function eraLd(bm::Cdouble,p::Array{Cdouble},q::Array{Cdouble},e::Array{Cdouble},em::Cdouble,dlim::Cdouble)
    p1 = zeros(3)
    ccall((:eraLd,liberfa),Void,
          (Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,Ptr{Cdouble}),
          bm,p,q,e,em,dlim,p1)
    p1
end

function eraLdn(l::Array{eraLDBODY}, ob::Array{Float64}, sc::Array{Float64})
    sn = zeros(3)
    n = length(l)
    ccall((:eraLdn, liberfa),Void,
          (Cint, Ptr{eraLDBODY}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          n, l, ob, sc, sn)
    sn
end

function eraLdsun(p::Array{Cdouble},e::Array{Cdouble},em::Cdouble)
    p1 = zeros(3)
    ccall((:eraLdsun,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Cdouble,Ptr{Cdouble}),
          p,e,em,p1)
    p1
end

function eraPmpx(rc::Cdouble,dc::Cdouble,pr::Cdouble,pd::Cdouble,px::Cdouble,rv::Cdouble,pmt::Cdouble,vob::Array{Cdouble})
    pco = zeros(3)
    ccall((:eraPmpx,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          rc,dc,pr,pd,px,rv,pmt,vob,pco)
    pco
end

function eraNumat(epsa::Real, dpsi::Real, deps::Real)
    rmatn = zeros((3,3))
    ccall((:eraNumat,liberfa),Void,
          (Float64, Float64, Float64, Ptr{Float64}),
          epsa, dpsi, deps, rmatn)
    rmatn
end

function eraP06e(date1::Cdouble,date2::Cdouble)
    eps0 = [0.]
    psia = [0.]
    oma = [0.]
    bpa = [0.]
    bqa = [0.]
    pia = [0.]
    bpia = [0.]
    epsa = [0.]
    chia = [0.]
    za = [0.]
    zetaa = [0.]
    thetaa = [0.]
    pa = [0.]
    gam = [0.]
    phi = [0.]
    psi = [0.]
    ccall((:eraP06e,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          date1,date2,eps0,psia,oma,bpa,bqa,pia,bpia,epsa,chia,za,zetaa,thetaa,pa,gam,phi,psi)
    eps0[1],psia[1],oma[1],bpa[1],bqa[1],pia[1],bpia[1],epsa[1],chia[1],za[1],zetaa[1],thetaa[1],pa[1],gam[1],phi[1],psi[1]
end

function eraP2s(p::Array{Cdouble})
    theta = [0.]
    phi = [0.]
    r = [0.]
    ccall((:eraP2s,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          p,theta,phi,r)
    theta[1], phi[1], r[1]
end

function eraP2pv(p::Array{Cdouble})
    pv = zeros((2,3))
    ccall((:eraP2pv,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          p,pv)
    pv
end

function eraPb06(date1::Cdouble,date2::Cdouble)
    bzeta = [0.]
    bz = [0.]
    btheta = [0.]
    ccall((:eraPb06,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          date1,date2,bzeta,bz,btheta)
    bzeta[1],bz[1],btheta[1]
end

function eraPfw06(date1::Cdouble,date2::Cdouble)
    gamb = [0.]
    phib = [0.]
    psib = [0.]
    epsa = [0.]
    ccall((:eraPfw06,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          date1,date2,gamb,phib,psib,epsa)
    gamb[1],phib[1],psib[1],epsa[1]
end

function eraPlan94(date1::Float64, date2::Float64, np::Integer)
    pv = zeros((2,3))
    i = ccall((:eraPlan94, liberfa),Cint,
              (Float64, Float64, Cint, Ptr{Float64}),
              date1, date2, np, pv)
    if i == -1
        error("illegal np,  not in range(1,8) for planet")
    elseif i == 1
        warn("year outside range(1000:3000)")
        return pv
    elseif i == 2
        error("computation failed to converge")
    elseif i == 0
        return pv
    end
end

function eraPm(p::Array{Cdouble})
    ccall((:eraPm,liberfa),Cdouble,(Ptr{Cdouble},),p)
end

function eraPmsafe(ra1::Cdouble,dec1::Cdouble,pmr1::Cdouble,pmd1::Cdouble,px1::Cdouble,rv1::Cdouble,ep1a::Cdouble,ep1b::Cdouble,ep2a::Cdouble,ep2b::Cdouble)
    ra2 = [0.]
    dec2 = [0.]
    pmr2 = [0.]
    pmd2 = [0.]
    px2 = [0.]
    rv2 = [0.]
    i = ccall((:eraPmsafe,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              ra1,dec1,pmr1,pmd1,px1,rv1,ep1a,ep1b,ep2a,ep2b,ra2,dec2,pmr2,pmd2,px2,rv2)
    if i == -1
        error("system error")
    elseif i == 1
        warn("distance overridden")
    elseif i == 2
        warn("excessive velocity")
    elseif i == 4
        error("solution didn't converge")
    end
    ra2[1],dec2[1],pmr2[1],pmd2[1],px2[1],rv2[1]
end

function eraPn(p::Array{Cdouble})
    r = [0.]
    u = zeros(3)
    ccall((:eraPn,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          p,r,u)
    r[1],u
end

function eraPpsp(a::Array{Cdouble},s::Cdouble,b::Array{Cdouble})
    apsb = zeros(3)
    ccall((:eraPpsp,liberfa),Void,
          (Ptr{Cdouble},Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          a,s,b,apsb)
    apsb
end

function eraPrec76(ep01::Cdouble,ep02::Cdouble,ep11::Cdouble,ep12::Cdouble)
    zeta = [0.]
    z = [0.]
    theta = [0.]
    ccall((:eraPrec76,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          ep01,ep02,ep11,ep12,zeta,z,theta)
    zeta[1],z[1],theta[1]
end

function eraPv2s(pv::Array{Cdouble})
    theta = [0.]
    phi = [0.]
    r = [0.]
    td = [0.]
    pd = [0.]
    rd = [0.]
   ccall((:eraPv2s,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          pv,theta,phi,r,td,pd,rd)
    theta[1], phi[1], r[1], td[1], pd[1], rd[1]
end

function eraPv2p(pv::Array{Cdouble})
    p = zeros(3)
    ccall((:eraPv2p,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          pv,p)
    p
end

function eraPvdpv(a::Array{Cdouble},b::Array{Cdouble})
    adb = zeros(2)
    ccall((:eraPvdpv,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          a,b,adb)
    adb
end

function eraPvm(pv::Array{Cdouble})
    s = [0.]
    r = [0.]
    ccall((:eraPvm,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          pv,r,s)
    r[1], s[1]
end

function eraPvstar(pv::Array{Cdouble})
    ra = [0.]
    dec = [0.]
    pmr = [0.]
    pmd = [0.]
    px = [0.]
    rv = [0.]
    i = ccall((:eraPvstar,liberfa),Cint,
              (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              pv,ra,dec,pmr,pmd,px,rv)
    if i == -1
        error("superluminal speed")
    elseif i == -2
        warn("null position vector")
        return ra[1],dec[1],pmr[1],pmd[1],px[1],rv[1]
    end
    ra[1],dec[1],pmr[1],pmd[1],px[1],rv[1]
end

function eraPvtob(elong::Cdouble,phi::Cdouble,height::Cdouble,xp::Cdouble,yp::Cdouble,sp::Cdouble,theta::Cdouble)
    pv = zeros((2,3))
    ccall((:eraPvtob,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
          elong,phi,height,xp,yp,sp,theta,pv)
    pv
end

function eraRefco(phpa::Cdouble,tk::Cdouble,rh::Cdouble,wl::Cdouble)
    refa = [0.]
    refb = [0.]
    ccall((:eraRefco,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          phpa,tk,rh,wl,refa,refb)
    refa[1],refb[1]
end

function eraPvu(dt::Cdouble,pv::Array{Cdouble})
    upv = zeros((2,3))
    ccall((:eraPvu,liberfa),Void,
          (Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          dt,pv,upv)
    upv
end

function eraPvup(dt::Cdouble,pv::Array{Cdouble})
    p = zeros(3)
    ccall((:eraPvup,liberfa),Void,
          (Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          dt,pv,p)
    p
end

function eraRm2v(r::Array{Cdouble})
    w = zeros(3)
    ccall((:eraRm2v,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          r,w)
    w
end

function eraRv2m(w::Array{Cdouble})
    r = zeros((3,3))
    ccall((:eraRv2m,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          w,r)
    r
end

function eraRxr(a::Array{Cdouble},b::Array{Cdouble})
    atb = zeros((3,3))
    ccall((:eraRxr,liberfa),
          Void,
          (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
          a,b,atb)
    atb
end

function eraS2c(theta::Cdouble,phi::Cdouble)
    c = zeros(3)
    ccall((:eraS2c,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble}),
          theta,phi,c)
    c
end

function eraS2p(theta::Cdouble,phi::Cdouble,r::Cdouble)
    p = zeros(3)
    ccall((:eraS2p,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
          theta,phi,r,p)
    p
end

function eraS2pv(theta::Cdouble,phi::Cdouble,r::Cdouble,td::Cdouble,pd::Cdouble,rd::Cdouble)
    pv = zeros((2,3))
    ccall((:eraS2pv,liberfa),Void,
          (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
          theta,phi,r,td,pd,rd,pv)
    pv
end

function eraS2xpv(s1::Cdouble,s2::Cdouble,pv::Array{Cdouble})
    spv = zeros((2,3))
    ccall((:eraS2xpv,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          s1,s2,pv,spv)
    spv
end

function eraStarpm(ra1::Cdouble,dec1::Cdouble,pmr1::Cdouble,pmd1::Cdouble,px1::Cdouble,rv1::Cdouble,ep1a::Cdouble,ep1b::Cdouble,ep2a::Cdouble,ep2b::Cdouble)
    ra2 = [0.]
    dec2 = [0.]
    pmr2 = [0.]
    pmd2 = [0.]
    px2 = [0.]
    rv2 = [0.]
    i = ccall((:eraStarpm,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
              ra1,dec1,pmr1,pmd1,px1,rv1,ep1a,ep1b,ep2a,ep2b,ra2,dec2,pmr2,pmd2,px2,rv2)
    if i == -1
        error("system error")
    elseif i == 1
        warn("distance overridden")
        return ra2[1],dec2[1],pmr2[1],pmd2[1],px2[1],rv2[1]
    elseif i == 2
        warn("excessive velocity")
        return ra2[1],dec2[1],pmr2[1],pmd2[1],px2[1],rv2[1]
    elseif i == 4
        error("solution didn't converge")
    end
    ra2[1],dec2[1],pmr2[1],pmd2[1],px2[1],rv2[1]
end

function eraStarpv(ra::Cdouble,dec::Cdouble,pmr::Cdouble,pmd::Cdouble,px::Cdouble,rv::Cdouble)
    pv = zeros((2,3))
    i = ccall((:eraStarpv,liberfa),Cint,
              (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
              ra,dec,pmr,pmd,px,rv,pv)
    if i == 1
        warn("distance overridden")
        return pv
    elseif i == 2
        warn("excessive speed ")
        return pv
    elseif i == 4
        error("solution didn't converge")
    end
    pv
end

function eraSxp(s::Cdouble,p::Array{Cdouble})
    sp = zeros(3)
    ccall((:eraSxp,liberfa),Void,
          (Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          s,p,sp)
    sp
end

function eraSxpv(s::Cdouble,pv::Array{Cdouble})
    spv = zeros((2,3))
    ccall((:eraSxpv,liberfa),Void,
          (Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          s,pv,spv)
    spv
end

function eraTr(r::Array{Cdouble})
    rt = zeros((3,3))
    ccall((:eraTr,liberfa),Void,
          (Ptr{Cdouble},Ptr{Cdouble}),
          r,rt)
    rt
end

function eraXy06(date1::Real,date2::Real)
    x = [0.]
    y = [0.]
    ccall((:eraXy06,liberfa),Void,
          (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
          date1,date2,x,y)
    x[1], y[1]
end

for f in (:eraAf2a,
          :eraTf2a,
          :eraTf2d)
    @eval begin
        function ($f)(s::Char,ideg::Integer,iamin::Integer,asec::Real)
            rad = [0.]
            i = ccall(($(Expr(:quote,f)),liberfa),Cint,
                       (Cchar,Cint,Cint,Cdouble,Ptr{Cdouble}),
                       s,ideg,iamin,asec,rad)
            @assert i == 0
            rad[1]
        end
    end
end

function eraIr()
    r = zeros((3,3))
    ccall((:eraIr,liberfa),Void,(Ptr{Cdouble},),r)
    r
end

function eraZp()
    p = zeros(3)
    #ccall((:eraZp,liberfa),Void,(Ptr{Cdouble},),p)
    p
end

function eraZpv()
    pv = zeros((2,3))
    #ccall((:eraZpv,liberfa),Void,(Ptr{Cdouble},),pv)
    pv
end

function eraZr()
    r = zeros((3,3))
    #ccall((:eraZr,liberfa),Void,(Ptr{Cdouble},),r)
    r
end

for f in (:eraA2af,
          :eraA2tf,
          :eraD2tf)
    @eval begin
        function ($f)(ndp::Integer, a::Float64)
            s = Vector{UInt8}(1)
            s[1] = '+'
            i = Int32[0, 0, 0, 0]
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cint, Float64, Ptr{UInt8}, Ptr{Cint}),
                  ndp, a, s, i)
            Char(s[1]), i[1], i[2], i[3], i[4]
        end
    end
end

for f in (:eraAnp,
          :eraAnpm)
    @eval begin
        function ($f)(a::Cdouble)
            ccall(($(Expr(:quote,f)),liberfa),Cdouble,(Cdouble,),a)
        end
    end
end

for f in (:eraBp00,
          :eraBp06)
    @eval begin
        function ($f)(a::Float64, b::Float64)
            rb = zeros((3,3))
            rp = zeros((3,3))
            rbp = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),
                  Void,
                  (Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                  a, b, rb, rp, rbp)
            rb, rp, rbp
        end
    end
end

for f in (:eraC2ixys,
          :eraPom00)
    @eval begin
        function ($f)(x::Cdouble, y::Cdouble, s::Cdouble)
            r = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
                  x, y, s, r)
            r
        end
    end
end

for f in (:eraC2tcio,
          :eraC2teqx)
    @eval begin
        function ($f)(rc2i::Array{Cdouble}, era::Cdouble, rpom::Array{Cdouble})
            rc2t = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Ptr{Cdouble},Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
                  rc2i,era,rpom,rc2t)
            rc2t
        end
    end
end

for f in (:eraC2ixy,
          :eraFw2m)
    @eval begin
        function ($f)(x::Cdouble, y::Cdouble, s::Cdouble, t::Cdouble)
            r = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
                  x, y, s, t, r)
            r
        end
    end
end

for f in (:eraC2t00a,
          :eraC2t00b,
          :eraC2t06a)
    @eval begin
        function ($f)(tta::Cdouble,ttb::Cdouble,uta::Cdouble,utb::Cdouble,xp::Cdouble,yp::Cdouble)
            rc2t =zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
                  tta,ttb,uta,utb,xp,yp,rc2t)
            rc2t
        end
    end
end

for f in (:eraC2tpe,
          :eraC2txy)
    @eval begin
        function ($f)(tta::Cdouble,ttb::Cdouble,uta::Cdouble,utb::Cdouble,x::Cdouble,y::Cdouble,xp::Cdouble,yp::Cdouble)
            rc2t =zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble}),
                  tta,ttb,uta,utb,x,y,xp,yp,rc2t)
            rc2t
        end
    end
end

for f in (:eraC2i00a,
          :eraC2i00b,
          :eraC2i06a,
          :eraEcm06,
          :eraNum00a,
          :eraNum00b,
          :eraNum06a,
          :eraNutm80,
          :eraPmat00,
          :eraPmat06,
          :eraPmat76,
          :eraPnm00a,
          :eraPnm00b,
          :eraPnm06a,
          :eraPnm80)
    @eval begin
        function ($f)(a::Float64, b::Float64)
            r = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Float64, Float64, Ptr{Float64}),
                  a, b, r)
            r
        end
    end
end

for f in (:eraPn00,
          :eraPn06)
    @eval begin
        function ($f)(date1::Float64, date2::Float64, dpsi::Float64, deps::Float64)
            epsa = [0.]
            rb = zeros((3,3))
            rp = zeros((3,3))
            rbp = zeros((3,3))
            rn = zeros((3,3))
            rbpn = zeros((3,3))
            ccall(($(Expr(:quote,f)), liberfa),Void,
                  (Float64, Float64, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                  date1, date2, dpsi, deps, epsa, rb, rp, rbp, rn, rbpn)
            epsa[1], rb, rp, rbp, rn, rbpn
        end
    end
end

for f in (:eraPmp,
          :eraPpp,
          :eraPxp)
    @eval begin
        function ($f)(a::Array{Cdouble},b::Array{Cdouble})
            ab = zeros(3)
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  a,b,ab)
            ab
        end
    end
end

for f in (:eraPvmpv,
          :eraPvppv,
          :eraPvxpv)
    @eval begin
        function ($f)(a::Array{Cdouble},b::Array{Cdouble})
            ab = zeros((2,3))
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  a,b,ab)
            ab
        end
    end
end

for f in (:eraPn00a,
          :eraPn00b,
          :eraPn06a)
    @eval begin
        function ($f)(date1::Float64, date2::Float64)
            dpsi = [0.]
            deps = [0.]
            epsa = [0.]
            rb = zeros((3,3))
            rp = zeros((3,3))
            rbp = zeros((3,3))
            rn = zeros((3,3))
            rbpn = zeros((3,3))
            ccall(($(Expr(:quote,f)), liberfa),Void,
                  (Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                  date1, date2, dpsi, deps, epsa, rb, rp, rbp, rn, rbpn)
            dpsi[1], deps[1], epsa[1], rb, rp, rbp, rn, rbpn
        end
    end
end

for f in (:eraNut00a,
          :eraNut00b,
          :eraNut06a,
          :eraNut80,
          :eraPr00,
          :eraG2icrs,
          :eraIcrs2g)
    @eval begin
        function ($f)(a::Float64, b::Float64)
            r1 = [0.]
            r2 = [0.]
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Float64, Float64, Ptr{Float64}, Ptr{Float64}),
                  a, b, r1, r2)
            r1[1], r2[1]
        end
    end
end

for f in (:eraFk52h,
          :eraH2fk5)
    @eval begin
        function ($f)(ra::Cdouble,dec::Cdouble,dra::Cdouble,ddec::Cdouble,px::Cdouble,rv::Cdouble)
            r = [0.]
            d = [0.]
            dr = [0.]
            dd = [0.]
            p = [0.]
            v = [0.]
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  ra,dec,dra,ddec,px,rv,r,d,dr,dd,p,v)
            r[1],d[1],dr[1],dd[1],p[1],v[1]
        end
    end
end

for f in (:eraFad03,
          :eraFae03,
          :eraFaf03,
          :eraFaju03,
          :eraFal03,
          :eraFalp03,
          :eraFama03,
          :eraFame03,
          :eraFane03,
          :eraFaom03,
          :eraFapa03,
          :eraFasa03,
          :eraFaur03,
          :eraFave03)
    @eval ($f)(d::Real) = ccall(($(Expr(:quote,f)),liberfa), Float64, (Float64,), d)
end

for f in (:eraEe00,
          :eraGmst00,
          :eraGmst06,
          :eraGst00a,
          :eraGst06a,
          :eraS00,
          :eraS06)
    @eval ($f)(d1::Real,d2::Real,t1::Real,t2::Real) = ccall(($(Expr(:quote,f)),liberfa), Float64, (Float64,Float64,Float64,Float64), d1,d2,t1,t2)
end

for f in (:eraEe00a,
          :eraEe00b,
          :eraEe06a,
          :eraEect00,
          :eraEo06a,
          :eraEpb,
          :eraEpj,
          :eraEqeq94,
          :eraEra00,
          :eraGmst82,
          :eraGst00b,
          :eraGst94,
          :eraObl06,
          :eraObl80,
          :eraS00a,
          :eraS00b,
          :eraS06a,
          :eraSp00)
    @eval ($f)(d1::Real, d2::Real) = ccall(($(Expr(:quote,f)),liberfa), Float64, (Float64,Float64), d1, d2)
end

for f in (:eraTaitt,
          :eraTaiutc,
          :eraTcbtdb,
          :eraTcgtt,
          :eraTdbtcb,
          :eraTttai,
          :eraTttcg,
          :eraUtctai)
    @eval begin
        function ($f)(a::Float64, b::Float64)
            r1 = [0.]
            r2 = [0.]
            i = ccall(($(Expr(:quote,f)),liberfa), Cint,
                      (Float64, Float64, Ptr{Float64}, Ptr{Float64}),
                      a, b, r1, r2)
            @assert i == 0
            r1[1], r2[1]
        end
    end
end

for f in (:eraEpb2jd,
          :eraEpj2jd)
    @eval begin
        function ($f)(d::Real)
            r1 = [0.]
            r2 = [0.]
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Float64, Ptr{Float64}, Ptr{Float64}),
                  d, r1, r2)
            r1[1], r2[1]
        end
    end
end

for f in (:eraTaiut1,
          :eraTdbtt,
          :eraTttdb,
          :eraTtut1,
          :eraUt1tai,
          :eraUt1tt,
          :eraUt1utc,
          :eraUtcut1)
    @eval begin
        function ($f)(a::Float64, b::Float64, c::Float64)
            r1 = [0.]
            r2 = [0.]
            i = ccall(($(Expr(:quote,f)),liberfa), Cint,
                      (Float64, Float64, Float64, Ptr{Float64}, Ptr{Float64}),
                      a, b, c, r1, r2)
            @assert i == 0
            r1[1], r2[1]
        end
    end
end

for f in (:eraPap,
          :eraPdp,
          :eraSepp)
    @eval begin
        function ($f)(a::Array{Cdouble}, b::Array{Cdouble})
            ccall(($(Expr(:quote,f)),liberfa),Cdouble,(Ptr{Cdouble},Ptr{Cdouble}),a,b)
        end
    end
end

for f in (:eraPas,
          :eraSeps)
    @eval begin
        function ($f)(al::Cdouble,ap::Cdouble,bl::Cdouble,bp::Cdouble)
            ccall(($(Expr(:quote,f)),liberfa),Cdouble,(Cdouble,Cdouble,Cdouble,Cdouble),al,ap,bl,bp)
        end
    end
end

for f in (:eraRxp,
          :eraTrxp)
    @eval begin
        function ($f)(r::Array{Cdouble}, p::Array{Cdouble})
            rp = zeros(3)
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  r,p,rp)
            rp
        end
    end
end

for f in (:eraRx,
          :eraRy,
          :eraRz)
    @eval begin
        function ($f)(a::Cdouble, r::Array{Cdouble})
            ccall(($(Expr(:quote,f)),liberfa),Void,
                  (Cdouble,Ptr{Cdouble}),
                  a,r)
            r
        end
    end
end

for f in (:eraRxpv,
          :eraTrxpv)
    @eval begin
        function ($f)(r::Array{Cdouble}, p::Array{Cdouble})
            rp = zeros((2,3))
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  r,p,rp)
            rp
        end
    end
end

for f in (:eraXys00a,
          :eraXys00b,
          :eraXys06a)
    @eval begin
        function ($f)(date1::Float64, date2::Float64)
            x = [0.]
            y = [0.]
            s = [0.]
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
                  date1,date2,x,y,s)
            x[1], y[1], s[1]
        end
    end
end

for f in (:eraLtecm,
          :eraLtp,
          :eraLtpb)
    @eval begin
        function ($f)(epj::Cdouble)
            rp = zeros((3,3))
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Cdouble,Ptr{Cdouble}),
                  epj,rp)
            rp
        end
    end
end

for f in (:eraLtpecl,
          :eraLtpequ)
    @eval begin
        function ($f)(epj::Cdouble)
            vec = zeros(3)
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Cdouble,Ptr{Cdouble}),
                  epj,vec)
            vec
        end
    end
end

for f in (:eraEceq06,
          :eraEqec06)
    @eval begin
        function ($f)(date1::Cdouble,date2::Cdouble,d1::Cdouble,d2::Cdouble)
            r1 = [0.0]
            r2 = [0.0]
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
                  date1, date2, d1, d2, r1,r2)
            r1[1], r2[1]
        end
    end
end

for f in (:eraLteceq,
          :eraLteqec)
    @eval begin
        function ($f)(epj::Cdouble,d1::Cdouble,d2::Cdouble)
            r1 = [0.0]
            r2 = [0.0]
            ccall(($(Expr(:quote,f)),liberfa), Void,
                  (Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}),
                  epj, d1, d2, r1,r2)
            r1[1], r2[1]
        end
    end
end

end # module


text:
function create_graph(start_node, end_node, link_length)
	@assert length(start_node)==length(end_node)

	no_node = max(maximum(start_node), maximum(end_node))
	no_arc = length(start_node)

	graph = Graph(no_node)
	distmx = Inf*ones(no_node, no_node)
	for i=1:no_arc
		add_edge!(graph, start_node[i], end_node[i])
		distmx[start_node[i], end_node[i]] = link_length[i]
	end
	return graph, distmx
end

function get_shortest_path(start_node::Array, end_node::Array, link_length::Array, origin::Int, destination::Int)
	@assert length(start_node)==length(end_node)
	@assert length(start_node)==length(link_length)

	graph, distmx = create_graph(start_node, end_node, link_length)

	state = dijkstra_shortest_paths(graph, origin, distmx)

	path = get_path(state, origin, destination)
	x = get_vector(path, start_node, end_node)

	len = state.dists[destination]

	return path, x, len
end

function get_path(state, origin, destination)
	_current = destination
	_rpath = _current
	_parent = -1

	while _parent != origin
		_parent = state.parents[_current]
		_rpath = [_rpath; _parent]
		_current = _parent
	end

	_path = zeros(Int, length(_rpath))
	for i=length(_rpath):-1:1
		_path[length(_rpath)-i+1] = _rpath[i]
	end
	return _path
end

function get_vector(path, start_node, end_node)
	x = zeros(Int, length(start_node))

	for i=1:length(path)-1
		st = path[i]
		en = path[i+1]

		for j=1:length(start_node)
			if start_node[j]==st && end_node[j]==en
				x[j] = 1
				break
			end
		end

	end
	return x
end

ext: jl


text:
# ===================================================================
#                       ERRORS
# ===================================================================

abstract type FEPrepException <: Exception end

struct FEPrepError <: FEPrepException
    msg::AbstractString
end

Base.showerror(io::IO, e::FEPrepError) = begin
    print(io, "FEPrep Error :: " * e.msg)
end

struct MeshPeriodicityError <: FEPrepException
    node_tag::Int
    group::AbstractString
end

Base.showerror(io::IO, e::MeshPeriodicityError) = begin
    print(io, "MeshPeriodicityError: Unable to find the corresponding node
    for a node, with tag ", e.node_tag, ", located in group ", e.group)
end

# ===================================================================
#                       INTEGRATION POINTS
# ===================================================================

abstract type AbstractIntegrationPoint end

struct IntegrationPoint1D <: AbstractIntegrationPoint
    wt::Float64
    r::Float64
end

struct IntegrationPoint2D <: AbstractIntegrationPoint
    wt::Float64
    r::Float64
    s::Float64
end

struct IntegrationPoint3D <: AbstractIntegrationPoint
    wt::Float64
    r::Float64
    s::Float64
    t::Float64
end

"""It returns tuple of integration points for a given finite elemenet type. One needs to provide either the order of gauss quadrature or number of integration points or both depending on the finite element type.

1. `num_ip`    : Triangular or Tetrahedral regions
2. `gq_order`  : Rectangular or Brick regions
3. `num_ip` and `gq_order` : Prism regions

"""
function get_integration_points(
    fele_type::DataType;
    gq_order::Int64 = -1,
    num_ip::Int64 = -1
)::Tuple{Vararg{AbstractIntegrationPoint}}
    @assert gq_order > 0 || num_ip > 0 "Gauss Quadrature order or Number of integration points depending on the element type."
    if fele_type == T2D2
        @assert gq_order > 0 "Line elements requires order of integration > 1 and it depends on the desired level of accuracy but $gq_order is given"
        return get_1D_integration_points(gq_order)
    elseif fele_type == CPS3
        @assert num_ip > 0 "Triangular region requires nuber of integration points. Valid numbers are 1, 3, 4"
        return get_triangular_region_integration_points(num_ip)
    elseif fele_type == CPS4
        @assert gq_order > 0 "Rectangular regions requires order of integration > 1 and it depends on the desired level of accuracy but $gq_order is given"
        return get_rectangular_region_integration_points(gq_order)
    elseif fele_type == C3D4
        @assert num_ip > 0 "Tetrahdral region requires number of integration points. Valid numbers are 1, 4, 5"
        return get_tetrahedral_region_integration_points(num_ip)
    elseif fele_type == C3D8
        @assert gq_order > 0 "Brick element regions require order of integration > 1 and it depends on the desired level of accuracy but $gq_order is given"
        return get_brick_region_integration_points(gq_order)
    elseif fele_type == C3D6
        @assert (num_ip > 0 && gq_order > 0) "Wedge element regions require order of integration > 1 and number of integration points > 1 which depends on the desired level of accuracy. Given order of integration = $gq_order, number of integration points = $num_ip"
        return get_wedge_region_integration_points(gq_order, num_ip)
    else
        @error "Encounteres unknown finite element type :: $fele_type"
    end
end

"""
Source:: https://pomax.github.io/bezierinfo/legendre-gauss.html
It contains data up to an order of 64.
"""
function get_1D_integration_points(
    order::Int64
)::Tuple{Vararg{IntegrationPoint1D}}
    wt::Vector{Float64} = []
    ξ::Vector{Float64} = []
    if order == 1
        wt, ξ = ([2.0,], [0.0,])
    elseif order == 2
        wt, ξ = ([1.0, 1.0], [-1 / sqrt(3), 1 / sqrt(3)])
    elseif order == 3
        wt, ξ = ([5 / 9, 8 / 9, 5 / 9], [-sqrt(0.6), 0.0, sqrt(0.6)])
    elseif order == 4
        wt, ξ = ([0.3478548451374538, 0.6521451548625461, 0.6521451548625461, 0.3478548451374538], [-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526])
    end
    return tuple((IntegrationPoint1D(aw, aξ) for (aw, aξ) in zip(wt, ξ))...)
end

function get_rectangular_region_integration_points(
    order::Int64,
)::Tuple{Vararg{IntegrationPoint2D}}
    ip_1d::Tuple{Vararg{IntegrationPoint1D}} = get_1D_integration_points(order)
    return tuple(
        (IntegrationPoint2D(a1Dpnti.wt * a1Dpntj.wt, a1Dpnti.r, a1Dpntj.r) for a1Dpnti in ip_1d for a1Dpntj in ip_1d)...
    )
end

function get_brick_region_integration_points(
    order::Int64,
)::Tuple{Vararg{IntegrationPoint3D}}
    ip_1d::Tuple{Vararg{IntegrationPoint1D}} = get_1D_integration_points(order)
    return tuple(
        (IntegrationPoint3D(a1Dpnti.wt * a1Dpntj.wt * a1Dpntk.wt, a1Dpnti.r, a1Dpntj.r, a1Dpntk.r) for a1Dpnti in ip_1d for a1Dpntj in ip_1d for a1Dpntk in ip_1d)...
    )
end

function get_wedge_region_integration_points(
    order::Int64,
    num_ip::Int64,
)::Tuple{Vararg{IntegrationPoint3D}}
    return tuple(
        IntegrationPoint3D(1 / 6, 1 / 6, 1 / 6, -1 / sqrt(3.0)),
        IntegrationPoint3D(1 / 6, 1 / 6, 1 / 6, 1 / sqrt(3.0)),
        IntegrationPoint3D(1 / 6, 1 / 6, 2 / 3, -1 / sqrt(3.0)),
        IntegrationPoint3D(1 / 6, 1 / 6, 2 / 3, 1 / sqrt(3.0)),
        IntegrationPoint3D(1 / 6, 2 / 3, 1 / 6, -1 / sqrt(3.0)),
        IntegrationPoint3D(1 / 6, 2 / 3, 1 / 6, 1 / sqrt(3.0)),
    )
end

"""
returns integration Points of triangular regions, given the number of gauss/integration points.
Source:: Concepts and applications of Finite Element Analysis by Robert D. Cook, Fourth Edition, Page Number: 267
"""
function get_triangular_region_integration_points(
    num_ip::Int64,  # number of integration points
)::Tuple{Vararg{IntegrationPoint2D}}
    wt::Vector{Float64} = []
    ξ::Vector{Float64} = []
    η::Vector{Float64} = []
    if num_ip == 1
        return (IntegrationPoint2D(1 / 2, 1 / 3, 1 / 3),)
    elseif num_ip == 3
        return (
            IntegrationPoint2D(1 / 6, 2 / 3, 1 / 6),
            IntegrationPoint2D(1 / 6, 1 / 6, 1 / 6),
            IntegrationPoint2D(1 / 6, 1 / 6, 2 / 3),
        )
    elseif num_ip == 4
        return (
            IntegrationPoint2D(-27 / 96, 1 / 3, 1 / 3),
            IntegrationPoint2D(25 / 96, 3 / 5, 1 / 5),
            IntegrationPoint2D(25 / 96, 1 / 5, 1 / 5),
            IntegrationPoint2D(25 / 96, 1 / 5, 3 / 5),
        )
    end
end

"""
returns integration Points of Tetrahedral regions, given the number of gauss/integration points.
Source:: Concepts and applications of Finite Element Analysis by Robert D. Cook, Fourth Edition, Page Number: 268
"""
function get_tetrahedral_region_integration_points(
    num_ip::Int64,  # number of integration points
)::Tuple{Vararg{IntegrationPoint3D}}
    wt::Vector{Float64} = []
    ξ::Vector{Float64} = []
    η::Vector{Float64} = []
    if num_ip == 1
        return (
            IntegrationPoint3D(1 / 6, 1 / 4, 1 / 4, 1 / 4),
        )
    elseif num_ip == 4
        a::Float64 = (5 + (3 * sqrt(5))) / 20
        b::Float64 = (5 - sqrt(5)) / 20
        return (
            IntegrationPoint3D(1 / 24, a, b, b),
            IntegrationPoint3D(1 / 24, b, b, b),
            IntegrationPoint3D(1 / 24, b, b, a),
            IntegrationPoint3D(1 / 24, b, a, b),
        )
    elseif num_ip == 5
        return (
            IntegrationPoint3D(3 / 40, 1 / 2, 1 / 6, 1 / 6),
            IntegrationPoint3D(3 / 40, 1 / 6, 1 / 6, 1 / 6),
            IntegrationPoint3D(3 / 40, 1 / 6, 1 / 6, 1 / 2),
            IntegrationPoint3D(3 / 40, 1 / 6, 1 / 2, 1 / 6),
            IntegrationPoint3D(-2 / 15, 1 / 4, 1 / 4, 1 / 4),
        )
    end
end

# ===================================================================
#                       MESH
# ===================================================================

abstract type AbstractNode end
abstract type AbstractElement end

abstract type AbstractFENode <: AbstractNode end
abstract type AbstractFElement <:AbstractElement end

abstract type SolidFElement <: AbstractFElement end
abstract type PlaneFElement <: AbstractFElement end
abstract type LineFElement <: AbstractFElement end

struct FENodePair
    tag::Union{String,Int}
    p_nodetag::Int  # parent node
    c_nodetag::Int  # child node
end

abstract type AbstractFElementSet end

"""Structure for storing a particular type of finite elements"""
@with_kw struct FiniteElementSet <: AbstractFElementSet
    tag::Union{Int64,String}
    elements::Vector{AbstractFElement}
    material::Union{Material, Nothing} = nothing
end

function get_num_ele(fele_set::FiniteElementSet)::Dict{DataType, Int}
    fele_type_summary::Dict{DataType, Int} = Dict{DataType, Int}()
    for a_fele in fele_set.elements
        fele_type = typeof(a_fele)
        if fele_type in keys(fele_type_summary)
            fele_type_summary[fele_type] += 1
        else
            fele_type_summary[fele_type] = 1
        end
    end
    return fele_type_summary
end

function get_num_ele(fele_sets::Vector{FiniteElementSet})
    fele_type_summary::Dict{DataType, Int} = Dict{DataType, Int}()
    for a_fele_set in fele_sets
        aset_fele_type_summary = get_num_ele(a_fele_set)
        for (ak, av) in aset_fele_type_summary
            if ak in keys(fele_type_summary)
                fele_type_summary[ak] += av
            else
                fele_type_summary[ak] = av
            end
        end
    end
    return fele_type_summary
end

""" PhaseFiniteElementConnectivity

Phase is the union of different types of elements, with same material
type"""
struct PhaseFiniteElementConnectivity
    tag::String
    ele_conn::Dict{Int, Matrix{Int}}
    material::Material
end

raw"""

    Elements local coordinates

    source: gmsh documentation

                C3D4                                 C3D10

                        v
                      .
                    ,/
                   /
                3                                     3
              ,/|`\                                 ,/|`\
            ,/  |  `\                             ,/  |  `\
          ,/    '.   `\                         ,7    '.   `6
        ,/       |     `\                     ,/       9     `\
      ,/         |       `\                 ,/         |       `\
     1-----------'.--------2 --> u         1--------5--'.--------2
      `\.         |      ,/                 `\.         |      ,/
         `\.      |    ,/                      `\.      |    ,10
            `\.   '. ,/                           `8.   '. ,/
               `\. |/                                `\. |/
                  `4                                    `4
                     `\.
                        ` w

         C3D8                       C3D20          Hexahedron27:

            v
     4----------3            3----13----2           3----13----2
     |\     ^   |\           |\         |\          |\         |\
     | \    |   | \          | 15       | 14        |15    24  | 14
     |  \   |   |  \         9  \       11 \        9  \ 20    11 \
     |   8------+---7        |   7----19+---6       |   7----19+---6
     |   |  +-- |-- | -> u   |   |      |   |       |22 |  26  | 23|
     1---+---\--2   |        0---+-8----1   |       0---+-8----1   |
      \  |    \  \  |         \  17      \  18       \ 17    25 \  18
       \ |     \  \ |         10 |        12|        10 |  21    12|
        \|      w  \|           \|         \|          \|         \|
         5----------6            4----16----5           4----16----5

              C3D6                 C3D15

                w
                ^
                |
                4                       3                      3
              ,/|`\                   ,/|`\                  ,/|`\
            ,/  |  `\               12  |  13              12  |  13
          ,/    |    `\           ,/    |    `\          ,/    |    `\
         5------+------6         4------14-----5        4------14-----5
         |      |      |         |      8      |        |      8      |
         |    ,/|`\    |         |      |      |        |    ,/|`\    |
         |  ,/  |  `\  |         |      |      |        |  15  |  16  |
         |,/    |    `\|         |      |      |        |,/    |    `\|
        ,|      |      |\        10     |      11       10-----17-----11
      ,/ |      1      | `\      |      0      |        |      0      |
     u   |    ,/ `\    |    v    |    ,/ `\    |        |    ,/ `\    |
         |  ,/     `\  |         |  ,6     `7  |        |  ,6     `7  |
         |,/         `\|         |,/         `\|        |,/         `\|
         2-------------3         1------9------2        1------9------2

"""

@with_kw struct C3D4 <: SolidFElement
    tag::Int64
    node_tags::SVector{4,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

@with_kw struct C3D6 <: SolidFElement
    tag::Int64
    node_tags::SVector{6,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

@with_kw struct C3D8 <: SolidFElement
    tag::Int64
    node_tags::SVector{8,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

raw"""
 Triangle:               Triangle6:          Triangle9/10:          Triangle12/15:

 v
 ^                                                                   2
 |                                                                   | \
 3                       2                    2                      9   8
 |`\                     |`\                  | \                    |     \
 |  `\                   |  `\                7   6                 10 (14)  7
 |    `\                 5    `4              |     \                |         \
 |      `\               |      `\            8  (9)  5             11 (12) (13) 6
 |        `\             |        `\          |         \            |             \
 1----------2 --> u      0-----3----1         0---3---4---1          0---3---4---5---1

 CPS4                         CPS8                  CPS9

       v
       ^
       |
 3-----------2          4-----7-----3           4-----7-----3
 |     |     |          |           |           |           |
 |     |     |          |           |           |           |
 |     +---- | --> u    8           6           8     9     6
 |           |          |           |           |           |
 |           |          |           |           |           |
 0-----------1          1-----5-----2           1-----5-----2

"""

@with_kw struct CPS3 <: PlaneFElement
    tag::Int64
    node_tags::SVector{3,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

@with_kw struct CPS4 <: PlaneFElement
    tag::Int64
    node_tags::SVector{4,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

@with_kw struct CPS8 <: PlaneFElement
    tag::Int64
    node_tags::SVector{8,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

@with_kw struct T2D2 <: LineFElement
    tag::Int64
    node_tags::SVector{2,Int64}
    ip_volumes::Tuple{Vararg{Float64}} = ()
    sfgrad_xyz::Tuple{Vararg{SMatrix}} = ()
end

""" Returns the GEOMETRIC DIMENSION of a finite element
"""
dim(fe_dtype::DataType)::Int64 = begin
    if fe_dtype in (C3D4, C3D6, C3D8)
        return 3
    elseif fe_dtype in (CPS3, CPS4, CPS8,)
        return 2
    elseif fe_dtype in (T2D2, )
        return 1
    else
        throw(error("Encountered incorrect FE data type $fe "))
    end
end

"""
It returns finite elemenet data type (with ABAQUS convention) for a given gmsh elemenet naming convention.
"""
function get_felement_data_type(
    gmsh_ele_type::Int64
)::DataType
    if gmsh_ele_type == 4
        return C3D4
    elseif gmsh_ele_type == 5
        return C3D8
    elseif gmsh_ele_type == 6
        return C3D6
    elseif gmsh_ele_type == 2
        return CPS3
    elseif gmsh_ele_type == 3
        return CPS4
    elseif gmsh_ele_type == 16
        return CPS8
    elseif gmsh_ele_type == 1
        return T2D2
    end
end

function get_felement_properties(
    element_type::Int64,
    gq_order::Int64 = -1,
    num_ip::Int64 = -1,
)
    ele_data_type::DataType = get_felement_data_type(element_type)
    return (
        dim(ele_data_type),
        ele_data_type,
        get_integration_points(ele_data_type; num_ip = num_ip, gq_order = gq_order),
        )
end

function get_num_nodes_per_ele(e::DataType)
    num_nodes_per_ele::Int = if (e===C3D8) || (e===CPS8)
        8
    elseif e===C3D6
        6
    elseif (e===C3D4) || (e===CPS4)
        4
    elseif e===CPS3
        3
    else
        @error "Element datatype $e is illegal, please provide correct element data type!."
    end
    return num_nodes_per_ele
end

# =============================================================================
#                                 SHAPE FUNCTIONS
# =============================================================================

"""
Returns `SMatrix` of `2x2` shape with each column for a Node/Shape function. First row is shape function evaluations and last row is shape function gradient evaluations at a given integration point.
"""
function shf_T2D2(
    intPoint::IntegrationPoint1D,
)::SMatrix{2,2,Float64}
    ξ::Float64 = intPoint.r
    @assert -1 ≤ ξ ≤ 1 "parent coordinate ξ must be between 0 & 1 but $(ξ) is found."
    N1::Float64 = 0.5 * (1 - ξ)
    N2::Float64 = 0.5 * (1 + ξ)
    dN1::Float64 = -0.5
    dN2::Float64 = 0.5
    return SMatrix{2,2,Float64}([N1 N2; dN1 dN2])
end

"""
Returns `SMatrix` of `3x3` shape with each column for a Node/Shape function. First row is shape function evaluations and remaining two rows are shape function gradient evaluations at a given integration point.
"""
function shf_CPS3(
    intPoint::IntegrationPoint2D,
)::SMatrix{3,3,Float64}
    ξ::Float64 = intPoint.r
    η::Float64 = intPoint.s
    @assert all(0 .≤ (ξ, η) .≤ 1) "parent coordinates have to be between 0 & 1"
    @assert ξ + η ≤ 1.0 "Sum of parent coordinates must be ≤ 1 for Triangular element"
    #
    N1::Float64 = 1 - ξ - η
    N2::Float64 = ξ
    N3::Float64 = η
    dN1::Vector{Float64} = [-1.0; -1.0]
    dN2::Vector{Float64} = [1.0; 0.0]
    dN3::Vector{Float64} = [0.0; 1.0]
    #
    return SMatrix{3,3,Float64}([[N1 N2 N3]; [dN1 dN2 dN3]])
end

"""
Returns `SMatrix` of `3x4` shape with each column for a Node/Shape function. First row is shape function evaluations and remaining two rows are shape function gradient evaluations at a given integration point.
"""
function shf_CPS4(
    intPoint::IntegrationPoint2D,
)::SMatrix{3,4,Float64}
    ξ::Float64 = intPoint.r
    η::Float64 = intPoint.s
    @assert all(-1 .≤ (ξ, η) .≤ 1) "parent coordinates for quad FElement have to be between 0 & 1"
    #
    N1::Float64 = (1 - ξ) * (1 - η)
    N2::Float64 = (1 + ξ) * (1 - η)
    N3::Float64 = (1 + ξ) * (1 + η)
    N4::Float64 = (1 - ξ) * (1 + η)
    dN1::Vector{Float64} = [η - 1.0; ξ - 1.0]
    dN2::Vector{Float64} = [1.0 - η; -1.0 - ξ]
    dN3::Vector{Float64} = [1.0 + η; 1.0 + ξ]
    dN4::Vector{Float64} = [-1.0 - η; 1.0 - ξ]
    #
    return 0.25 * SMatrix{3,4,Float64}([[N1 N2 N3 N4]; [dN1 dN2 dN3 dN4]])
end

"""
Returns `SMatrix` of `4x4` shape with each column for a Node/Shape function. First row is shape function evaluations and remaining three rows are shape function gradient evaluations at a given integration point.
"""
function shf_C3D4(
    intPoint::IntegrationPoint3D,
)::SMatrix{4,4,Float64}
    ξ::Float64 = intPoint.r
    η::Float64 = intPoint.s
    ζ::Float64 = intPoint.t
    @assert all(0 .≤ (ξ, η, ζ) .≤ 1) "parent coordinates have to be between 0 & 1 but $ξ, $η, $ζ are supplied."
    @assert ξ + η + ζ ≤ 1.0 "Sum of parent coordinates must be ≤ 1"
    #
    N1::Float64 = 1 - ξ - η - ζ
    N2::Float64 = ξ
    N3::Float64 = η
    N4::Float64 = ζ
    dN1::Vector{Float64} = [-1.0; -1.0; -1.0]
    dN2::Vector{Float64} = [1.0; 0.0; 0.0]
    dN3::Vector{Float64} = [0.0; 1.0; 0.0]
    dN4::Vector{Float64} = [0.0; 0.0; 1.0]
    #
    return SMatrix{4,4,Float64}([[N1 N2 N3 N4]; [dN1 dN2 dN3 dN4]])
end

"""
Returns `SMatrix` of `4x6` shape with each column for a Node/Shape function. First row is shape function evaluations and remaining three rows are shape function gradient evaluations at a given integration point.
"""
function shf_C3D6(
    intPoint::IntegrationPoint3D,
)::SMatrix{4,6,Float64}
    ξ::Float64 = intPoint.r
    η::Float64 = intPoint.s
    ζ::Float64 = intPoint.t
    @assert all(0 .≤ (ξ, η) .≤ 1) "parent coordinates ξ, η have to be between 0 & 1"
    @assert -1 ≤ ζ ≤ 1 "Parent coordinate ζ along the axis of the wedge element must be between -1 and 1"
    @assert ξ + η ≤ 1.0 "Sum of parent coordinates must be ≤ 1"
    #
    kn::Float64, kp::Float64  = (1.0 - ζ) * 0.5, (1.0 + ζ) * 0.5
    N1::Float64 = (1.0 - ξ - η) * kn
    N2::Float64 = ξ * kn
    N3::Float64 = η * kn
    N4::Float64 = (1.0 - ξ - η) * kp
    N5::Float64 = ξ * kp
    N6::Float64 = η * kp
    dN1::Vector{Float64} = [-kn; -kn; (-1.0 + ξ + η) * 0.5]
    dN2::Vector{Float64} = [kn; 0.0; -ξ * 0.5]
    dN3::Vector{Float64} = [0.0; kn; -η * 0.5]
    dN4::Vector{Float64} = [-kp; -kp; (1.0 - ξ - η) * 0.5]
    dN5::Vector{Float64} = [kp; 0.0; ξ * 0.5]
    dN6::Vector{Float64} = [0.0; kp; η * 0.5]
    #
    return SMatrix{4,6,Float64}([[N1 N2 N3 N4 N5 N6]; [dN1 dN2 dN3 dN4 dN5 dN6]])
end

"""
Returns `SMatrix` of `4x8` shape with each column for a Node/Shape function. First row is shape function evaluations and remaining three rows are shape function gradient evaluations at a given integration point.

"""
function shf_C3D8(
    intPoint::IntegrationPoint3D,
)::SMatrix{4,8,Float64}
    ξ::Float64 = intPoint.r
    η::Float64 = intPoint.s
    ζ::Float64 = intPoint.t
    @assert all(-1 .≤ (ξ, η, ζ) .≤ 1) " All parent coordinates have to be between 0 & 1"
    N1::Float64 = (1.0 - ξ) * (1.0 - η) * (1.0 - ζ)
    N2::Float64 = (1.0 + ξ) * (1.0 - η) * (1.0 - ζ)
    N3::Float64 = (1.0 + ξ) * (1.0 + η) * (1.0 - ζ)
    N4::Float64 = (1.0 - ξ) * (1.0 + η) * (1.0 - ζ)
    N5::Float64 = (1.0 - ξ) * (1.0 - η) * (1.0 + ζ)
    N6::Float64 = (1.0 + ξ) * (1.0 - η) * (1.0 + ζ)
    N7::Float64 = (1.0 + ξ) * (1.0 + η) * (1.0 + ζ)
    N8::Float64 = (1.0 - ξ) * (1.0 + η) * (1.0 + ζ)
    dN1::Vector{Float64} = [(η - 1.0) * (1.0 - ζ); (ξ - 1.0) * (1.0 - ζ); (η - 1.0) * (1.0 - ξ)]
    dN2::Vector{Float64} = [(1.0 - η) * (1.0 - ζ); (ξ + 1.0) * (ζ - 1.0); (η - 1.0) * (1.0 + ξ)]
    dN3::Vector{Float64} = [(1.0 + η) * (1.0 - ζ); (ξ + 1.0) * (1.0 - ζ); (-η - 1.0) * (1.0 + ξ)]
    dN4::Vector{Float64} = [(1.0 + η) * (ζ - 1.0); (1.0 - ξ) * (1.0 - ζ); (η + 1.0) * (ξ - 1.0)]
    dN5::Vector{Float64} = [(η - 1.0) * (1.0 + ζ); (ξ - 1.0) * (1.0 + ζ); (1.0 - ξ) * (1.0 - η)]
    dN6::Vector{Float64} = [(1.0 - η) * (1.0 + ζ); (-ξ - 1.0) * (1.0 + ζ); (1.0 - η) * (1.0 + ξ)]
    dN7::Vector{Float64} = [(1.0 + η) * (1.0 + ζ); (ξ + 1.0) * (1.0 + ζ); (η + 1.0) * (1.0 + ξ)]
    dN8::Vector{Float64} = [(-1.0 - η) * (1.0 + ζ); (1.0 - ξ) * (1.0 + ζ); (1.0 - ξ) * (1.0 + η)]
    return 0.125 .* SMatrix{4,8,Float64}([[N1 N2 N3 N4 N5 N6 N7 N8]; [dN1 dN2 dN3 dN4 dN5 dN6 dN7 dN8]])
end

# ===================================>

"""
    get_shape_functions_info() -> Dict{AbstractIntegrationPoint, SMatrix{1+dim, num_shf, Float64}}

# Args
- `**element_type**: DataType`
- `**integration_points**:Tuple{Vararg{AbstractIntegrationPoint}}`

"""
function get_shape_functions_info(
    element_type::DataType,
    integration_points::Tuple{Vararg{AbstractIntegrationPoint}};
)::Dict{AbstractIntegrationPoint,SMatrix}
    #
    shape_functions_info = begin
        if element_type == CPS3
            shf_CPS3
        elseif element_type == CPS4
            shf_CPS4
        elseif element_type == C3D4
            shf_C3D4
        elseif element_type == C3D6
            shf_C3D6
        elseif element_type == C3D8
            shf_C3D8
        elseif element_type == T2D2
            shf_T2D2
        else
            throw(FEPrepError("Unable to get the shape functions for the element type $element_type"))
        end
    end
    return Dict((
        a_int_pnt => shape_functions_info(a_int_pnt)
    ) for a_int_pnt in integration_points)
end

# ===================================================================
#                       ANALYSIS
# ===================================================================

abstract type AbstractFEAnalysis end

abstract type FEA_3Danalysis <: AbstractFEAnalysis end
abstract type FEA_2Danalysis <: AbstractFEAnalysis end

#
struct Elastic_3DFEA <: FEA_3Danalysis end
struct Piezo_Electric_3DFEA <: FEA_3Danalysis end
struct Piezo_Magnetic_3DFEA <: FEA_3Danalysis end
struct Magneto_Electro_Elastic_3DFEA <: FEA_3Danalysis end
struct Thermal_Conductivity_3DFEA <: FEA_3Danalysis end
#
struct Thermo_Elastic_3DFEA <: FEA_3Danalysis end
struct Thermo_Piezo_Electric_3DFEA <: FEA_3Danalysis end
struct Thermo_Piezo_Magnetic_3DFEA <: FEA_3Danalysis end
struct Thermo_Magneto_Electro_Elastic_3DFEA <: FEA_3Danalysis end
#
abstract type Elastic_2DFEA <: FEA_2Danalysis end
struct PlaneStress_2DFEA <: Elastic_2DFEA end
struct PlaneStrain_2DFEA <: Elastic_2DFEA end

"""
Returns dof associated with each node for a given analysis and dimension
"""
function get_nodal_dof(
    analysis_type::DataType,
)::Tuple{Vararg{Symbol}}
    if analysis_type in (Elastic_3DFEA, Thermo_Elastic_3DFEA)
        return (:ux, :uy, :uz,)
    elseif analysis_type === Elastic_2DFEA
        return (:ux, :uy,)
    elseif analysis_type === Thermal_Conductivity_3DFEA
        return (:T,)
    elseif analysis_type === Piezo_Electric_3DFEA
        return (:ux, :uy, :uz, :ϕe,)
    elseif analysis_type === Magneto_Electro_Elastic_3DFEA
        return (:ux, :uy, :uz, :ϕe, :ϕm,)
    end
end

"""
    get_material_matrix_size(analysis_type) -> Tuple(Int, Int)

"""
function get_material_matrix_size(
    analysis_type::DataType,
)
    if analysis_type in (Elastic_3DFEA, Thermo_Elastic_3DFEA, )
        return (6, 6)
    elseif analysis_type in (Elastic_2DFEA, Thermal_Conductivity_3DFEA, )
        return (3, 3)
    elseif analysis_type == Piezo_Electric_3DFEA
        return (9, 9)
    elseif analysis_type == Magneto_Electro_Elastic_3DFEA
        return (12, 12)
    else
        @error "Encountered an unknown analysis type for returning material matrix size"
    end
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#                         Archives/Dev code
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# abstract type AbstractMesh end
# abstract type AbstractFEMesh <: AbstractMesh end

# struct Mesh <: AbstractMesh
#     nodes::StructVector{AbstractNode}
#     elements::StructVector{AbstractElement}
# end

# struct FEMesh <: AbstractFEMesh
#     nodes::StructVector{AbstractFENode}
#     elements::StructVector{AbstractFElement}
# end

# @with_kw struct FENode <: AbstractFENode
#     tag::Int
#     x::Float64
#     y::Float64 = 0.0
#     z::Float64 = 0.0
# end

# struct C3D{N} <: SolidFElement
#     tag::Int
#     nodes::SVector{N, FENode}
# end

# struct FiniteElement{N} <: AbstractFElement
#     tag::Int64
#     nodes::SVector{N, FENode}
#     ip_volumes::SVector{N, Float64}
#     shf_grad_xyz::SVector{N, SMatrix}
# end

# """
#     get_felement_data_type() -> DataType

#  Returns finite elemenet data type converted to a common data type used in the
#  present code.

#  # Args
#  - `element_ID::Union{Int, String, Symbol}`
#  - `convention::AbstractString`, can be either "GMSH" or "ABAQUS"

#  # KWArgs
#  - `sos::String="PS"`

# """
# function get_felement_data_type(
#     element_ID::Union{Int, String, Symbol},
#     convention::String;  #
#     sos::String="PS",
# )::DataType
#     convention = uppercase(convention)

#     if convention=="GMSH"
#         if element_ID == 4
#             return C3D{4}
#         elseif element_ID == 5
#             return C3D{8}
#         elseif element_ID == 6
#             return C3D{6}
#         elseif element_ID == 2
#             if uppercase(sos)=="PS"
#                 return CPS{3}
#             elseif uppercase(sos)=="PE"
#                 return CPE{3}
#             end
#         elseif element_ID == 3
#             if uppercase(sos)=="PS"
#                 return CPS{4}
#             elseif uppercase(sos)=="PE"
#                 return CPE{4}
#             end
#         elseif element_ID == 16
#             if uppercase(sos)=="PS"
#                 return CPS{8}
#             elseif uppercase(sos)=="PE"
#                 return CPE{8}
#             end
#         elseif element_ID == 1
#             return T2D{2}
#         end
#     elseif convention=="ABAQUS"
#         if element_ID == "C3D4"
#             return C3D{4}
#         elseif element_ID == "C3D8"
#             return C3D{8}
#         elseif element_ID == "C3D6"
#             return C3D{6}
#         elseif element_ID == "CPS3"
#             return CPS{3}
#         elseif element_ID == "CPE3"
#             return CPE{3}
#         elseif element_ID == "CPS4"
#             return CPS{4}
#         elseif element_ID == "CPE4"
#             return CPE{4}
#         elseif element_ID == "CPS8"
#             return CPS{8}
#         elseif element_ID == "CPE8"
#             return CPE{8}
#         elseif element_ID == "T2D2"
#             return T2D{2}
#         end
#     else
#         throw(FEPrepError("At present, only gmsh and ABAQUS based element codes are handled "))
#     end
# end

# #
# gmsh_ele_type(et::T2D2)::Int = 1
# #
# gmsh_ele_type(et::CPS{3})::Int = 2
# gmsh_ele_type(et::CPS{4})::Int = 3
# gmsh_ele_type(et::CPS{8})::Int = 16
# #
# gmsh_ele_type(et::C3D{4})::Int = 4
# gmsh_ele_type(et::C3D{6})::Int = 6
# gmsh_ele_type(et::C3D{8})::Int = 5

# """
# Returns SMatrix of adjacent local node labels for each local node of the
# given element type.

# Each column holds a local node neighbours in order

# """
# function adjacent_local_node_indices(
#     el_dtype::DataType
# )::SMatrix
#     if el_dtype == T2D{2}
#         return SMatrix{1,2,Int64}([2 1])
#     elseif el_dtype in (CPS{3}, CPE{3}, )
#         return SMatrix{2,3,Int64}([[2, 3] [3, 1] [1, 2]])
#     elseif el_dtype in (CPS{4}, CPE{4}, )
#         return SMatrix{2,4,Int64}([[2, 4] [3, 1] [4, 2] [1, 3]])
#     elseif el_dtype in (CPS{8}, CPE{8}, )
#         return SMatrix{2,8,Int64}(
#             [[5, 8] [5, 6] [6, 7] [7, 8] [1, 2] [2, 3] [3, 4] [4, 1]]
#         )
#     elseif el_dtype == C3D{4}
#         return SMatrix{3,4,Int64}([[2, 3, 4] [3, 1, 4] [1, 2, 4] [1, 2, 3]])
#     elseif el_dtype == C3D{6}
#         return SMatrix{3,6,Int64}(
#             [[2, 3, 4] [1, 3, 5] [1, 2, 6] [1, 5, 6] [2, 4, 6] [3, 4, 5]]
#         )
#     elseif el_dtype == C3D{8}
#         return SMatrix{3,8,Int64}(
#             [[2, 4, 5] [3, 1, 6] [7, 2, 4] [1, 8, 3] [1, 8, 6] [2, 5, 7] [3, 6, 8] [4, 5, 7]]
#         )
#     end
# end


text:
using StatsBase
function unique(x; count = false)
    s = [v for v in Set(x)]
    if count
        c = zeros(1, size(s)[1])
        for i in 1:length(x)
            for j in 1:length(s)
                if x[i] == s[j]
                    c[j] += 1
                    break
                end
            end
        end
        k = zeros(Int, size(c))
        for i in 1:length(c)
            k[i] = Int(c[i])
        end
        return s, k
    end
    return s
end
function DecisionTree(x, t; alpha = 0.01)
    function get_branching_shold(xs)
        u = unique(xs)
        return (u[2:length(u)] + u[1:length(u)-1]) / 2
    end
    function delta_gini_index(left, right, alpha, branch_count)
        n_left = length(left)
        n_right = length(right)
        n_total = n_left + n_right
        _, counts = unique(left, count = true)
        ratio_classes = counts / n_left
        left_gain = (n_left / n_total) * (1 .- sum(ratio_classes .^ 2))
        _, counts = unique(right, count = true)
        ratio_classes = counts / n_right
        right_gain = (n_right / n_total) * (1 .- sum(ratio_classes .^ 2))
        return left_gain + right_gain + alpha * branch_count
    end
    function branch(x, t, alpha, branch_count)
        gains = []
        rules = []
        for feature_id in 1 : size(x)[2]
            xs = x[:, feature_id]
            thresholds = get_branching_shold(xs)
            for th in thresholds
                f = findall(x -> x < th, xs)
                left_y = [t[i] for i in f]
                f = findall(x -> x >= th, xs)
                right_y = [t[i] for i in f]
                gain = delta_gini_index(left_y, right_y, alpha, branch_count)
                push!(gains, gain)
                if length(rules) == 0
                    rules = [feature_id th]
                end
                rules = vcat(rules, [feature_id th])
            end
        end
        best_rule = rules[argmin(gains), :]
        best_gini = minimum(gains)
        feature_id = Int(best_rule[1])
        threshold = best_rule[2]
        split = x[:, feature_id] .< threshold
        newaxis = [CartesianIndex()]
        left_x, left_t, right_x, right_t = [], [], [], []
        global left_x, left_t, right_x, right_t
        for i in 1 : length(split)
            if split[i]
                if length(left_x) == 0
                    left_x = x[i, :]
                    left_x = left_x[newaxis, :]
                    left_t = t[:, i]
                    left_t = left_t[newaxis, :]
                else
                    kx = x[i, :]
                    kx = kx[newaxis, :]
                    left_x = vcat(left_x, kx)
                    kx = t[:, i]
                    kx = kx[newaxis, :]
                    left_t = vcat(left_t, kx)
                end
            else
                if length(right_x) == 0
                    right_x = x[i, :]
                    right_x = right_x[newaxis, :]
                    right_t = t[:, i]
                    right_t = right_t[newaxis, :]
                else
                    kx = x[i, :]
                    kx = kx[newaxis, :]
                    right_x = vcat(right_x, kx)
                    kx = t[:, i]
                    kx = kx[newaxis, :]
                    right_t = vcat(right_t, kx)
                end
            end
        end
        return left_x, left_t, right_x, right_t, feature_id, threshold, best_gini
    end
    function grow(x, t, classes, gini, alpha, branch_count)
        if size(x)[1] == size(t)[1]
            t = t'
        end
        uniques, counts = unique(t, count = true)
        counter = Dict(zip(uniques, counts))
        class_count = [(c in keys(counter)) ? counter[c] : 0 for c in classes]
        this = Dict("class_count" => class_count, "feature_id" => nothing, "threshold" => nothing, "left" => nothing, "right" => nothing)
        if length(t) == 1
            return this
        elseif length(uniques) == 1
            return this
        end
        a = true
        for i in 1:size(x)[2]
            a = a && length(unique(x[:, i])) == 1
        end
        if a
            return this
        end
        left_x, left_t, right_x, right_t, feature_id, threshold, best_gini = branch(x, t, alpha, branch_count)
        if best_gini < gini
            branch_count += 1
            left = grow(left_x, left_t, classes, best_gini, alpha, branch_count)
            right = grow(right_x, right_t, classes, best_gini, alpha, branch_count)
            return Dict("class_count" => class_count, "feature_id" => feature_id, "threshold" => threshold, "left" => left, "right" => right)
        else
            return this
        end
    end
    classes = unique(t)
    return grow(x, t, classes, Inf, alpha, 1), classes
end
function tree_predict(xs, tree, classes)
    function predict_one(x, tree)
        while tree["feature_id"] != nothing
            tree = (x[tree["feature_id"]] < tree["threshold"]) ? tree["left"] : tree["right"]
        end
        return argmax(tree["class_count"])
    end
    predicts = [predict_one(xs[i, :], tree) for i in 1:size(xs)[1]]
    return predicts
end
function bootstrap(x, t, n_trees)
    n_features = size(x)[2]
    n_features_forest = Int(floor(sqrt(n_features)))
    bootstrapped_x = zeros((length(t), n_features_forest, n_trees))
    bootstrapped_t = zeros(Int, n_trees, length(t))
    using_feature = zeros(Int, n_trees, n_features_forest)
    newaxis = [CartesianIndex()]
    global bootstrapped_x
    global bootstrapped_t
    global using_feature
    for i in 1 : n_trees
        ind = rand(1:length(t), length(t))
        col = sample(1 : n_features, n_features_forest, replace = false)
        k = x[ind, col]
        bootstrapped_x[:, :, i] = k[:, :, newaxis]
        bootstrapped_t[i, :] = t[ind]
        using_feature[i, :] = col
    end
    return bootstrapped_x, bootstrapped_t, using_feature
end
function fit(x, t; n_trees = 10, alpha = 0.01)
    classes = unique(t)
    forest = Array{Any}(undef, n_trees)
    classes = Array{Any}(undef, n_trees)
    bootstrapped_x, bootstrapped_t, using_feature = bootstrap(x, t, n_trees)
    for i in 1:n_trees
        forest[i], classes[i] = DecisionTree(bootstrapped_x[:, :, i], bootstrapped_t[i, :], alpha = alpha)
    end
    return forest, classes, using_feature
end
function predict(x, forest, classes, using_feature)
    solution = []
    for (tree, classes, feature) in zip(forest, classes, using_feature)
        s = tree_predict(x[:, feature], tree, classes)
        if solution == []
            solution = s
        else
            solution = hcat(solution, s)
        end
    end
    predicts = []
    for i in 1:length(solution[:, 1])
        class, counts = unique(solution[i, :], count = true)
        push!(predicts, class[argmax(counts)[2]])
    end
    d = zeros(Int, size(classes[1]))
    for i in 1 : length(d)
        d[i] = i
    end
    d = Dict(zip(d, classes[1]))
    for i in 1 : length(predicts)
        predicts[i] = d[predicts[i]]
    end
    return predicts
end


text:
abstract EnKF_iter
typealias IntPairs Tuple{Int64, Int64}

immutable EAKF <: EnKF_iter end

function update( iiz :: Array{IntPairs}, loc :: Localization, kf :: EnKF_iter,
                 rraw :: Array{Float64,2},
                 obs:: Array{Float64,1},error_obs :: Array{Float64, 1})
  res = copy(rraw)
  for ii = iiz
    (ob, err_ob, ra, y) = iter_select(ii, loc, res, obs, error_obs)
    upyy = iter_update(ii, loc, kf, ra, y, ob, err_ob)
    iter_refill!(ii, loc, res, upyy)
  end
  return res
end

function iter_update( ii :: IntPairs, loc :: NearPoints, kf :: EAKF,
                      raw0 :: Array{Float64,2}, y0 :: Array{Float64,1},
                      obs :: Float64, error_obs :: Float64,
                    )
  y1 = y0 - mean(y0)
  raw1 = raw0 .- mean(raw0,2)
  covr = (raw1 * y1) / sum(y1.^2)

  ## Updating the mean
  sigma_ratio = err_obs^2 / (error_obs^2 + var(y0))
  delta_ym = (obs - mean(y0) ) * (1-sigma_ratio)
  delta_yd :: Array{Float64,1} = y1 (sqrt(sigma_ratio)-1)
  delta_upyy0 :: Array{Float64,1} = covr* delta_y
  delta_upyy1 :: Array{Float64,2} = covr * delta_yd'
  return ((delta_upyy0 .+ delta_upyy1) + raw0)
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct LogQuantUniform <: AbstractDistDelayed
    low::NestedFloat
    high::NestedFloat
    q::NestedFloat
end

function logquantuniform(low::Float64, high::Float64, q::Float64)::Float64
    return exp(quantuniform(low, high, q))
end


text:
# -----------------------------------------------------------------------
# Experiments pre processing
# -----------------------------------------------------------------------
# Required packages
using JLD, DataFrames, Dates, Statistics

"""
    function op2column(A::AbstractArray, v::AbstractVector, op=+)

## Description
Apply operation for the columns of a matrix
## Arguments
- A: AbstractArray
- v: vector to perform elementwise operation
- op: operation
## Example
```jldoctest
julia> A = [[1,2,3] [4,5,6] [7,8,9]]
3×3 Array{Int64,2}:
 1  4  7
 2  5  8
 3  6  9

julia> v = [1,1,1]
 3-element Array{Int64,1}:
 1
 1
 1

julia> op2column(A, v, op=-)
3×3 Array{Int64,2}:
 0  3  6
 1  4  7
 2  5  8
```
"""
function op2column(A::AbstractArray, v::AbstractVector; op=+)
    return op.(A, reshape(repeat(v, inner=size(A)[1]), size(A)))
end

"""
    correct_zero(fname::String, fzero_init::String, fzero_end::String)

## Decription
Correct zero by interpolating zero measurement at the begining and at the end.
## Arguments
- fname: string (or array) with the name of the test file
- fzero_init: name of the file with the zero read before the test
- fzero_end: name of the file with the zero read after the test
"""
function correct_zero(fname::String, fzero_init::String, fzero_end::String)

    p,h = load(fname, "pressure", "h")
    pz1,hz1 = load(fzero_init, "pressure", "h")
    pz2,hz2 = load(fzero_end, "pressure", "h")

    Δhz = (hz2 - hz1)/(Millisecond(1.0)*1000) # dt in seconds
    Δh = (h .- hz1) ./ (Millisecond(1.0)*1000)

    pz1_m = [mean(pz1[:,i]) for i in 1:16]
    pz2_m = [mean(pz2[:,i]) for i in 1:16]
    pz_m = pz2_m .- pz1_m

    Δp = [pz1_m + (x/Δhz)*pz_m for x in Δh]

    p_corr = [op2column(p[i], Δp[i], op=-) for i in 1:length(p)]

    return p_corr
end

function correct_zero(fname::AbstractArray, fzero_init::String, fzero_end::String)

    nf = length(fname)
    pz1,hz1 = load(fzero_init, "pressure", "h")
    pz2,hz2 = load(fzero_end, "pressure", "h")

    pz1_m = [mean(pz1[:,i]) for i in 1:16]
    pz2_m = [mean(pz2[:,i]) for i in 1:16]
    pz_m = pz2_m .- pz1_m

    Δhz = (hz2 - hz1)/(Millisecond(1.0)*1000) # dt in seconds

    p = []
    h = []
    for i in 1:nf
        ptemp,htemp = load(fname[i], "pressure", "h")
        push!(p, ptemp)
        push!(h, htemp)
    end

    Δh = (h .- hz1) ./ (Millisecond(1.0)*1000)
    Δp = [pz1_m + (x/Δhz)*pz_m for x in Δh]

    p_corr = [op2column(p[i], Δp[i], op=-) for i in 1:length(p)]

    return p_corr
end

"""
    calc_meanpressure(p::AbstractArray, nch)

## Description
Compute the mean pressure values given an array of scanivalve measurements
## Arguments
- p: pressure measured data (already corrected if desired)
- nch: scanivalve channels to consider
"""
function calc_meanpressure(p::AbstractArray, nch)
    pm = [mean(i[:,nch]) for i in p]
end

"""
    calc_meancp(p::AbstractArray, nch)

## Description
Compute the mean pressure coefficient values given an array of scanivalve measurements
## Arguments
- p: pressure measured data (already corrected if desired)
- nch: scanivalve channels to consider
- chref: reference channel (dinamic pressure measurement)
"""
function calc_meancp(p::AbstractArray, nch, chref=16)
    pref = [mean(i[:,chref]) for i in p]
    cpm = calc_meanpressure(p, nch) ./ pref
    return [cpm, pref]
end

"""
    calc_stdpressure(p::AbstractArray, nch)

## Description
Compute the standard deviation pressure values given an array of scanivalve measurements
## Arguments
- p: pressure measured data (already corrected if desired)
- nch: scanivalve channels to consider
"""
function calc_stdpressure(p::AbstractArray, nch)
    pstd = [std(i[:,k]) for i in p for k in nch]
    pstd = reshape(pstd, (length(nch), length(p)))
    return pstd
end

"""
    calc_stdcp(p::AbstractArray, nch)

## Description
Compute the standard deviation pressure coefficient values given an array of scanivalve measurements
## Arguments
- p: pressure measured data (already corrected if desired)
- nch: scanivalve channels to consider
- chref: reference channel (dinamic pressure measurements)
"""
function calc_stdcp(p::AbstractArray, nch, chref=16)
    cpstd = [std(i[:,k] ./ i[:,chref]) for i in p for k in nch]
    cpstd = reshape(cpstd, (length(nch), length(p)))
    return cpstd
end