text:
#####
##### Constant diffusivity
#####

"""
    ConstantIsotropicDiffusivity{T}

    ConstantIsotropicDiffusivity(; ν, κ)

Returns parameters for a constant isotropic diffusivity closure with constant viscosity `ν`
and constant thermal diffusivity `κ`. `ν` and `κ` may represent molecular diffusivities or
turbulent eddy diffusivities.

By default, a molecular viscosity of ``ν = 1.05×10⁻⁶`` m²/s and a molecular thermal
diffusivity of ``κ = 1.46×10⁻⁷`` m²/s is used, corresponding to the use of no turbulent
diffusivity closure at all. These molecular values are the approximate viscosity and
thermal diffusivity for seawater at 20°C and 35 psu, according to Sharqawy et al.,
"Thermophysical properties of seawater: A review of existing correlations and data" (2010).
"""
Base.@kwdef struct ConstantIsotropicDiffusivity{T} <: IsotropicDiffusivity{T}
    ν :: T = ν₀
    κ :: T = κ₀
end

ConstantIsotropicDiffusivity(T; kwargs...) =
    typed_keyword_constructor(T, ConstantIsotropicDiffusivity; kwargs...)

calc_diffusivities!(diffusivities, grid, closure::ConstantIsotropicDiffusivity,
                    args...) = nothing

@inline ∇_κ_∇c(i, j, k, grid, c, closure::ConstantIsotropicDiffusivity, args...) = (
      closure.κ / grid.Δx^2 * δx²_c2f2c(grid, c, i, j, k)
    + closure.κ / grid.Δy^2 * δy²_c2f2c(grid, c, i, j, k)
    + closure.κ / grid.Δz^2 * δz²_c2f2c(grid, c, i, j, k)
)

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, u, v, w, args...) = (
      closure.ν / grid.Δx^2 * δx²_f2c2f(grid, u, i, j, k)
    + closure.ν / grid.Δy^2 * δy²_f2e2f(grid, u, i, j, k)
    + closure.ν / grid.Δz^2 * δz²_f2e2f(grid, u, i, j, k)
)

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, u, v, w, args...) = (
      closure.ν / grid.Δx^2 * δx²_f2e2f(grid, v, i, j, k)
    + closure.ν / grid.Δy^2 * δy²_f2c2f(grid, v, i, j, k)
    + closure.ν / grid.Δz^2 * δz²_f2e2f(grid, v, i, j, k)
)

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, u, v, w, args...) = (
      closure.ν / grid.Δx^2 * δx²_f2e2f(grid, w, i, j, k)
    + closure.ν / grid.Δy^2 * δy²_f2e2f(grid, w, i, j, k)
    + closure.ν / grid.Δz^2 * δz²_f2c2f(grid, w, i, j, k)
)

#####
##### Constant anisotropic diffusivity (tensor diffusivity with heterogeneous
#####                                       diagonal components)
#####

"""
    ConstantAnisotropicDiffusivity{T}

    ConstantAnisotropicDiffusivity(; νh, νv, κh, κv)

Returns parameters for a constant anisotropic diffusivity closure with constant horizontal
and vertical viscosities `νh`, `νv` and constant horizontal and vertical thermal diffusivities
`κh`, `κv`. `ν` and `κ` may represent molecular diffusivities or turbulent eddy diffusivities.

By default, a molecular viscosity of ``ν = 1.05×10⁻⁶`` m²/s and a molecular thermal
diffusivity of ``κ = 1.46×10⁻⁷`` m²/s is used, corresponding to the use of no turbulent
diffusivity closure at all. These molecular values are the approximate viscosity and
thermal diffusivity for seawater at 20°C and 35 psu, according to Sharqawy et al.,
"Thermophysical properties of seawater: A review of existing correlations and data" (2010).
"""
Base.@kwdef struct ConstantAnisotropicDiffusivity{T} <: TensorDiffusivity{T}
    νh :: T = ν₀
    νv :: T = ν₀
    κh :: T = κ₀
    κv :: T = κ₀
end

calc_diffusivities!(diffusivities, grid, closure::ConstantAnisotropicDiffusivity,
                    args...) = nothing

ConstantAnisotropicDiffusivity(T; kwargs...) =
    typed_keyword_constructor(T, ConstantAnisotropicDiffusivity; kwargs...)

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, u, v, w, K) = (
      closure.νh * ∂x²_faa(i, j, k, grid, u)
    + closure.νh * ∂y²_aca(i, j, k, grid, u)
    + closure.νv * ∂z²_aac(i, j, k, grid, u)
    )

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, u, v, w, K) = (
      closure.νh * ∂x²_caa(i, j, k, grid, v)
    + closure.νh * ∂y²_afa(i, j, k, grid, v)
    + closure.νv * ∂z²_aac(i, j, k, grid, v)
    )

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, u, v, w, K) = (
      closure.νh * ∂x²_caa(i, j, k, grid, w)
    + closure.νh * ∂y²_aca(i, j, k, grid, w)
    + closure.νv * ∂z²_aaf(i, j, k, grid, w)
    )

@inline ∇_κ_∇c(i, j, k, grid, c, closure::ConstantAnisotropicDiffusivity, K) = (
      closure.κh * ∂x²_caa(i, j, k, grid, c)
    + closure.κh * ∂y²_aca(i, j, k, grid, c)
    + closure.κv * ∂z²_aac(i, j, k, grid, c)
    )


text:
using Pkg
using Test

mkdir("workspace")
cd("workspace")
isempty(readdir()) || error("workspace not empty")

# Put Pipenv under the current directory.
ENV["PIPENV_VENV_IN_PROJECT"] = "true"
# See: https://pipenv.readthedocs.io/en/latest/advanced/#pipenv.environments.PIPENV_VENV_IN_PROJECT

# Setup Python packages
write("Pipfile", """
[[source]]
url = "https://pypi.org/simple"
verify_ssl = true
name = "pypi"

[packages]

[dev-packages]
""")

run(`pipenv install`)

# Setup Julia packages
write("Project.toml", """
[deps]
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
""")

Pkg.activate(".")
Pkg.add(PackageSpec(url=ENV["PYVIRTUALENV_JL_TEST_TARGET"]))

@debug "Importing PyVirtualenv..."
import PyVirtualenv
PyVirtualenv.activate_pipenv()

@debug "Importing PyCall..."
using PyCall

sys = pyimport("sys")
@testset for key in [:executable, :prefix, :exec_prefix]
    # Note that this works because PIPENV_VENV_IN_PROJECT is set.
    @test startswith(sys[key], pwd())
end

# Cleaning up. Not required, but handy when it's manually invoked:
cd("..")
rm("workspace"; force=true, recursive=true)

ext: jl


text:
#> \brief \b DGBTRF
#
#  =========== DOCUMENTATION ===========
#
# Online html documentation available at
#            http://www.netlib.org/lapack/explore-html/
#
#> \htmlonly
#> Download DGBTRF + dependencies
#>
#> [TGZ]
#>
#> [ZIP]
#>
#> [TXT]
#> \endhtmlonly
#
#  Definition:
#  ===========
#
#       SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO )
#
#       .. Scalar Arguments ..
#       INTEGER            INFO, KL, KU, LDAB, M, N
#       ..
#       .. Array Arguments ..
#       INTEGER            IPIV( * )
#       DOUBLE PRECISION   AB( LDAB, * )
#       ..
#
#
#> \par Purpose:
#  =============
#>
#> \verbatim
#>
#> DGBTRF computes an LU factorization of a real m-by-n band matrix A
#> using partial pivoting with row interchanges.
#>
#> This is the blocked version of the algorithm, calling Level 3 BLAS.
#> \endverbatim
#
#  Arguments:
#  ==========
#
#> \param[in] M
#> \verbatim
#>          M is INTEGER
#>          The number of rows of the matrix A.  M >= 0.
#> \endverbatim
#>
#> \param[in] N
#> \verbatim
#>          N is INTEGER
#>          The number of columns of the matrix A.  N >= 0.
#> \endverbatim
#>
#> \param[in] KL
#> \verbatim
#>          KL is INTEGER
#>          The number of subdiagonals within the band of A.  KL >= 0.
#> \endverbatim
#>
#> \param[in] KU
#> \verbatim
#>          KU is INTEGER
#>          The number of superdiagonals within the band of A.  KU >= 0.
#> \endverbatim
#>
#> \param[in,out] AB
#> \verbatim
#>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
#>          On entry, the matrix A in band storage, in rows KL+1 to
#>          2*KL+KU+1; rows 1 to KL of the array need not be set.
#>          The j-th column of A is stored in the j-th column of the
#>          array AB as follows:
#>          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
#>
#>          On exit, details of the factorization: U is stored as an
#>          upper triangular band matrix with KL+KU superdiagonals in
#>          rows 1 to KL+KU+1, and the multipliers used during the
#>          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
#>          See below for further details.
#> \endverbatim
#>
#> \param[in] LDAB
#> \verbatim
#>          LDAB is INTEGER
#>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
#> \endverbatim
#>
#> \param[out] IPIV
#> \verbatim
#>          IPIV is INTEGER array, dimension (min(M,N))
#>          The pivot indices; for 1 <= i <= min(M,N), row i of the
#>          matrix was interchanged with row IPIV(i).
#> \endverbatim
#>
#> \param[out] INFO
#> \verbatim
#>          INFO is INTEGER
#>          = 0: successful exit
#>          < 0: if INFO = -i, the i-th argument had an illegal value
#>          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
#>               has been completed, but the factor U is exactly
#>               singular, and division by zero will occur if it is used
#>               to solve a system of equations.
#> \endverbatim
#
#  Authors:
#  ========
#
#> \author Univ. of Tennessee
#> \author Univ. of California Berkeley
#> \author Univ. of Colorado Denver
#> \author NAG Ltd.
#
#> \date December 2016
#
#> \ingroup doubleGBcomputational
#
#> \par Further Details:
#  =====================
#>
#> \verbatim
#>
#>  The band storage scheme is illustrated by the following example, when
#>  M = N = 6, KL = 2, KU = 1:
#>
#>  On entry:                       On exit:
#>
#>      *    *    *    +    +    +       *    *    *   u14  u25  u36
#>      *    *    +    +    +    +       *    *   u13  u24  u35  u46
#>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
#>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
#>     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
#>     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
#>
#>  Array elements marked * are not used by the routine; elements marked
#>  + need not be set on entry, but are required by the routine to store
#>  elements of U because of fill-in resulting from the row interchanges.
#> \endverbatim
#>
#  =====================================================================
      function dgbtrf( M::Integer, N::Integer, KL::Integer, KU::Integer, AB::Array, LDAB::Integer, IPIV::Vector{Integer}, INFO::Integer)
#
#  -- LAPACK computational routine (version 3.7.0) --
#  -- LAPACK is a software package provided by Univ. of Tennessee,    --
#  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
#     December 2016
#
#     .. Scalar Arguments ..
#      INTEGER            INFO, KL, KU, LDAB, M, N
#     ..
#     .. Array Arguments ..
#      INTEGER            IPIV( * )
#      DOUBLE PRECISION   AB( ldab, * )
      @assert(size(ab,1) == LDAB)
#     ..
#
#  =====================================================================
#
#     .. Parameters ..

#    DOUBLE PRECISION   ONE, ZERO
#      parameter                ( one = 1.0d+0, zero = 0.0d+0 )
      INTEGER            NBMAX, LDWORK
#      parameter                ( nbmax = 64, ldwork = nbmax+1 )
      nbmax = 64
      ldwork = nbmax+1
#     ..
#     .. Local Scalars ..
#      INTEGER            I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
#     $                   ju, k2, km, kv, nb, nw
#      DOUBLE PRECISION   TEMP
#     ..
#     .. Local Arrays ..
#      DOUBLE PRECISION   WORK13( ldwork, nbmax ),
#     $                   work31( ldwork, nbmax )
      work13 = zeros(ldwork,nbmax)
      work31 = zeros(ldwork,nbmax)
#     ..
#     .. External Functions ..
      INTEGER            IDAMAX, ILAENV
      EXTERNAL           idamax, ilaenv
#     ..
#     .. External Subroutines ..
      EXTERNAL           dcopy, dgbtf2, dgemm, dger, dlaswp, dscal,
     $                   dswap, dtrsm, xerbla
#     ..
#     .. Intrinsic Functions ..
      INTRINSIC          max, min
#     ..
#     .. Executable Statements ..
#
#     KV is the number of superdiagonals in the factor U, allowing for
#     fill-in
#
      kv = ku + kl
#
#     Test the input parameters.
#
      info = 0
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( kl.LT.0 ) THEN
         info = -3
      ELSE IF( ku.LT.0 ) THEN
         info = -4
      ELSE IF( ldab.LT.kl+kv+1 ) THEN
         info = -6
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DGBTRF', -info )
         RETURN
      END IF
#
#     Quick return if possible
#
      IF( m.EQ.0 .OR. n.EQ.0 )
     $   RETURN
#
#     Determine the block size for this environment
#
      nb = ilaenv( 1, 'DGBTRF', ' ', m, n, kl, ku )
#
#     The block size must not exceed the limit set by the size of the
#     local arrays WORK13 and WORK31.
#
      nb = min( nb, nbmax )
#
      IF( nb.LE.1 .OR. nb.GT.kl ) THEN
#
#        Use unblocked code
#
         CALL dgbtf2( m, n, kl, ku, ab, ldab, ipiv, info )
      ELSE
#
#        Use blocked code
#
#        Zero the superdiagonal elements of the work array WORK13
#
         DO 20 j = 1, nb
            DO 10 i = 1, j - 1
               work13( i, j ) = zero
   10       CONTINUE
   20    CONTINUE
#
#        Zero the subdiagonal elements of the work array WORK31
#
         DO 40 j = 1, nb
            DO 30 i = j + 1, nb
               work31( i, j ) = zero
   30       CONTINUE
   40    CONTINUE
252 *
#        Gaussian elimination with partial pivoting
#
#        Set fill-in elements in columns KU+2 to KV to zero
#
         DO 60 j = ku + 2, min( kv, n )
            DO 50 i = kv - j + 2, kl
               ab( i, j ) = zero
   50       CONTINUE
   60    CONTINUE
#
#        JU is the index of the last column affected by the current
#        stage of the factorization
#
         ju = 1
#
         DO 180 j = 1, min( m, n ), nb
            jb = min( nb, min( m, n )-j+1 )
#
#           The active part of the matrix is partitioned
#
#              A11   A12   A13
#              A21   A22   A23
#              A31   A32   A33
#
#           Here A11, A21 and A31 denote the current block of JB columns
#           which is about to be factorized. The number of rows in the
#           partitioning are JB, I2, I3 respectively, and the numbers
#           of columns are JB, J2, J3. The superdiagonal elements of A13
#           and the subdiagonal elements of A31 lie outside the band.
#
            i2 = min( kl-jb, m-j-jb+1 )
            i3 = min( jb, m-j-kl+1 )
#
#           J2 and J3 are computed after JU has been updated.
#
#           Factorize the current block of JB columns
#
            DO 80 jj = j, j + jb - 1
#
#              Set fill-in elements in column JJ+KV to zero
#
               IF( jj+kv.LE.n ) THEN
                  DO 70 i = 1, kl
                     ab( i, jj+kv ) = zero
   70             CONTINUE
               END IF
#
#              Find pivot and test for singularity. KM is the number of
#              subdiagonal elements in the current column.
#
               km = min( kl, m-jj )
               jp = idamax( km+1, ab( kv+1, jj ), 1 )
               ipiv( jj ) = jp + jj - j
               IF( ab( kv+jp, jj ).NE.zero ) THEN
                  ju = max( ju, min( jj+ku+jp-1, n ) )
                  IF( jp.NE.1 ) THEN
#
#                    Apply interchange to columns J to J+JB-1
#
                     IF( jp+jj-1.LT.j+kl ) THEN
#
                        CALL dswap( jb, ab( kv+1+jj-j, j ), ldab-1,
     $                              ab( kv+jp+jj-j, j ), ldab-1 )
                     ELSE
#
#                       The interchange affects columns J to JJ-1 of A31
#                       which are stored in the work array WORK31
#
                        CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
     $                              work31( jp+jj-j-kl, 1 ), ldwork )
                        CALL dswap( j+jb-jj, ab( kv+1, jj ), ldab-1,
     $                              ab( kv+jp, jj ), ldab-1 )
                     END IF
                  END IF
#
#                 Compute multipliers
#
                  CALL dscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
     $                        1 )
#
#                 Update trailing submatrix within the band and within
#                 the current block. JM is the index of the last column
#                 which needs to be updated.
#
                  jm = min( ju, j+jb-1 )
                  IF( jm.GT.jj )
     $               CALL dger( km, jm-jj, -one, ab( kv+2, jj ), 1,
     $                          ab( kv, jj+1 ), ldab-1,
     $                          ab( kv+1, jj+1 ), ldab-1 )
               ELSE
#
#                 If pivot is zero, set INFO to the index of the pivot
#                 unless a zero pivot has already been found.
#
                  IF( info.EQ.0 )
     $               info = jj
               END IF
#
#              Copy current column of A31 into the work array WORK31
#
               nw = min( jj-j+1, i3 )
               IF( nw.GT.0 )
     $            CALL dcopy( nw, ab( kv+kl+1-jj+j, jj ), 1,
     $                        work31( 1, jj-j+1 ), 1 )
   80       CONTINUE
            IF( j+jb.LE.n ) THEN
#
#              Apply the row interchanges to the other blocks.
#
               j2 = min( ju-j+1, kv ) - jb
               j3 = max( 0, ju-j-kv+1 )
#
#              Use DLASWP to apply the row interchanges to A12, A22, and
#              A32.
#
               CALL dlaswp( j2, ab( kv+1-jb, j+jb ), ldab-1, 1, jb,
     $                      ipiv( j ), 1 )
#
#              Adjust the pivot indices.
#
               DO 90 i = j, j + jb - 1
                  ipiv( i ) = ipiv( i ) + j - 1
   90          CONTINUE
#
#              Apply the row interchanges to A13, A23, and A33
#              columnwise.
#
               k2 = j - 1 + jb + j2
               DO 110 i = 1, j3
                  jj = k2 + i
                  DO 100 ii = j + i - 1, j + jb - 1
                     ip = ipiv( ii )
                     IF( ip.NE.ii ) THEN
                        temp = ab( kv+1+ii-jj, jj )
                        ab( kv+1+ii-jj, jj ) = ab( kv+1+ip-jj, jj )
                        ab( kv+1+ip-jj, jj ) = temp
                     END IF
  100             CONTINUE
  110          CONTINUE
#
#              Update the relevant part of the trailing submatrix
#
               IF( j2.GT.0 ) THEN
#
#                 Update A12
#
                  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
     $                        jb, j2, one, ab( kv+1, j ), ldab-1,
     $                        ab( kv+1-jb, j+jb ), ldab-1 )
#
                  IF( i2.GT.0 ) THEN
#
#                    Update A22
#
                     CALL dgemm( 'No transpose', 'No transpose', i2, j2,
     $                           jb, -one, ab( kv+1+jb, j ), ldab-1,
     $                           ab( kv+1-jb, j+jb ), ldab-1, one,
     $                           ab( kv+1, j+jb ), ldab-1 )
                  END IF
#
                  IF( i3.GT.0 ) THEN
#
#                    Update A32
#
                     CALL dgemm( 'No transpose', 'No transpose', i3, j2,
     $                           jb, -one, work31, ldwork,
     $                           ab( kv+1-jb, j+jb ), ldab-1, one,
     $                           ab( kv+kl+1-jb, j+jb ), ldab-1 )
                  END IF
               END IF
#
               IF( j3.GT.0 ) THEN
#
#                 Copy the lower triangle of A13 into the work array
#                 WORK13
#
                  DO 130 jj = 1, j3
                     DO 120 ii = jj, jb
                        work13( ii, jj ) = ab( ii-jj+1, jj+j+kv-1 )
  120                CONTINUE
  130             CONTINUE
#
#                 Update A13 in the work array
#
                  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
     $                        jb, j3, one, ab( kv+1, j ), ldab-1,
     $                        work13, ldwork )
#
                  IF( i2.GT.0 ) THEN
#
#                    Update A23
#
                     CALL dgemm( 'No transpose', 'No transpose', i2, j3,
     $                           jb, -one, ab( kv+1+jb, j ), ldab-1,
     $                           work13, ldwork, one, ab( 1+jb, j+kv ),
     $                           ldab-1 )
                  END IF
#
                  IF( i3.GT.0 ) THEN
#
#                    Update A33
#
                     CALL dgemm( 'No transpose', 'No transpose', i3, j3,
     $                           jb, -one, work31, ldwork, work13,
     $                           ldwork, one, ab( 1+kl, j+kv ), ldab-1 )
                  END IF
#
#                 Copy the lower triangle of A13 back into place
#
                  DO 150 jj = 1, j3
                     DO 140 ii = jj, jb
                        ab( ii-jj+1, jj+j+kv-1 ) = work13( ii, jj )
  140                CONTINUE
  150             CONTINUE
               END IF
            ELSE
#
#              Adjust the pivot indices.
#
               DO 160 i = j, j + jb - 1
                  ipiv( i ) = ipiv( i ) + j - 1
  160          CONTINUE
            END IF
#
#           Partially undo the interchanges in the current block to
#           restore the upper triangular form of A31 and copy the upper
#           triangle of A31 back into place
#
            DO 170 jj = j + jb - 1, j, -1
               jp = ipiv( jj ) - jj + 1
               IF( jp.NE.1 ) THEN
#
#                 Apply interchange to columns J to JJ-1
#
                  IF( jp+jj-1.LT.j+kl ) THEN
#
#                    The interchange does not affect A31
#
                     CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
     $                           ab( kv+jp+jj-j, j ), ldab-1 )
                  ELSE
#
#                    The interchange does affect A31
#
                     CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
     $                           work31( jp+jj-j-kl, 1 ), ldwork )
                  END IF
               END IF
#
#              Copy the current column of A31 back into place
#
               nw = min( i3, jj-j+1 )
               IF( nw.GT.0 )
     $            CALL dcopy( nw, work31( 1, jj-j+1 ), 1,
     $                        ab( kv+kl+1-jj+j, jj ), 1 )
  170       CONTINUE
  180    CONTINUE
      END IF
#
      RETURN
#
#     End of DGBTRF
#
      END


text:
using Distributed
AI_THREADS = haskey(ENV, "AI_THREADS") ? parse(Int64, ENV["AI_THREADS"]) : 1
if AI_THREADS > 1 addprocs(AI_THREADS) end # total: procs + main proc
@info "Start $(nprocs()) processes"

@everywhere begin
  using Distributed
  @info "Start process $(myid())"

  @eval using Pkg
  @eval pkg"activate ." # activate Project
  #@everywhere Pkg.instantiate()

  # add current path
  cd(@__DIR__)
  push!(LOAD_PATH,@__DIR__)
end

################################################################################

const FILE_SCRIPT_ERROR = abspath(joinpath(@__DIR__,"../logs/error_script.log"))
const FILE_OUTPUT = abspath(joinpath(@__DIR__,"../logs/output.log"))
const FILE_REVISELOG = abspath(joinpath(@__DIR__,"../logs/revise.log"))

clearLog() = open(FILE_OUTPUT,"w+") do f write(f,"") end

################################################################################

messages=Channel{String}(Inf)

function sendToPrinter(msg::String)
  global messages
  put!(messages, msg)
end

function Printer()
  global FILE_OUTPUT, messages
  while true
    msgs = ""
    while isready(messages) msgs*=take!(messages) end
    if length(msgs)>0
      open(FILE_OUTPUT,"a+") do f write(f,msgs) end
    end
    sleep(1)
  end
end

################################################################################

logs = nothing

try
    @eval using Revise
    @eval using Base.CoreLogging: Debug

    rlogger = Revise.debug_logger()
    Revise.ReviseLogger(Revise.LogRecord[], Debug)
    logs = filter(r->r.level==Debug, rlogger.logs)
catch ex
  error("Revise Error! Could not load Revise: $ex")
end

################################################################################
# preloads
@everywhere using cywwtaip
@everywhere cywwtaip.init()

################################################################################

print("\nPress enter to start...")
input = readline() # wait

################################################################################

# Green Threads
@async Printer()
#@async Helper()

################################################################################
input = ""

while input != "q" && input != "quit"
  global logs, input

  try

    clearLog()
    println("Start Script...")
    @eval using Client
    Client.main(ARGS;printer=sendToPrinter)

  catch ex
    println("Script Error! For more details look in $FILE_SCRIPT_ERROR")
    open(FILE_SCRIPT_ERROR,"w+") do f Base.showerror(f, ex, catch_backtrace()) end
  end

  if logs != nothing open(FILE_OUTPUT, "w+") do io for log in logs println(io, log) end end end

  print("\nPress enter to restart...")
  input = readline()
  println("Reload Script...")

  try
    Revise.revise()
  catch ex
    println("Revise Error! For more details look in $FILE_REVISELOG")
    open(FILE_REVISELOG,"w+") do f Base.showerror(f, ex, catch_backtrace()) end
  end

  println("\33[2J") # clear console: https://rosettacode.org/wiki/Terminal_control/Clear_the_screen#Julia
end

max_stars_repo_name: Gilga/AI-Cywwtaip

max_issues_repo_name: Gilga/AI-Cywwtaip

max_forks_repo_name: Gilga/AI-Cywwtaip


text:
# Common methods for all `Values` objects

const ScalarValues{ξdim,xdim,T,shape} = Union{CellScalarValues{ξdim,xdim,T,shape},FaceScalarValues{ξdim,xdim,T,shape}}
const VectorValues{ξdim,xdim,T,shape} = Union{CellVectorValues{ξdim,xdim,T,shape},FaceVectorValues{ξdim,xdim,T,shape}}

getnbasefunctions(cv::Values) = size(cv.N, 1)
getngeobasefunctions(cv::Values) = size(cv.M, 1)

getn_scalarbasefunctions(cv::ScalarValues) = size(cv.N, 1)
getn_scalarbasefunctions(cv::VectorValues{ξdim,xdim}) where {ξdim,xdim} = size(cv.N, 1) ÷ xdim

"""
    reinit!(cv::CellValues, x::Vector)
    reinit!(bv::FaceValues, x::Vector, face::Int)

Update the `CellValues`/`FaceValues` object for a cell or face with coordinates `x`.
The derivatives of the shape functions, and the new integration weights are computed.
"""
reinit!

"""
    getnquadpoints(fe_v::Values)

Return the number of quadrature points for the `Values` object.
"""
getnquadpoints(fe::Values) = length(fe.qr_weights)

"""
    getdetJdV(fe_v::Values, q_point::Int)

Return the product between the determinant of the Jacobian and the quadrature
point weight for the given quadrature point: ``\\det(J(\\mathbf{x})) w_q``

This value is typically used when one integrates a function on a
finite element cell or face as

``\\int\\limits_\\Omega f(\\mathbf{x}) d \\Omega \\approx \\sum\\limits_{q = 1}^{n_q} f(\\mathbf{x}_q) \\det(J(\\mathbf{x})) w_q``

``\\int\\limits_\\Gamma f(\\mathbf{x}) d \\Gamma \\approx \\sum\\limits_{q = 1}^{n_q} f(\\mathbf{x}_q) \\det(J(\\mathbf{x})) w_q``

"""
@inline getdetJdV(cv::CellValues, q_point::Int) = cv.detJdV[q_point]
@inline getdetJdV(bv::FaceValues, q_point::Int) = bv.detJdV[q_point, bv.current_face[]]

"""
    shape_value(fe_v::Values, q_point::Int, base_function::Int)

Return the value of shape function `base_function` evaluated in
quadrature point `q_point`.
"""
@inline shape_value(cv::CellValues, q_point::Int, base_func::Int) = cv.N[base_func, q_point]
@inline shape_value(bv::FaceValues, q_point::Int, base_func::Int) = bv.N[base_func, q_point, bv.current_face[]]

@inline geometric_value(cv::CellValues, q_point::Int, base_func::Int) = cv.M[base_func, q_point]
@inline geometric_value(bv::FaceValues, q_point::Int, base_func::Int) = bv.M[base_func, q_point, bv.current_face[]]

"""
    shape_gradient(fe_v::Values, q_point::Int, base_function::Int)

Return the gradient of shape function `base_function` evaluated in
quadrature point `q_point`.
"""
@inline shape_gradient(cv::CellValues, q_point::Int, base_func::Int) = cv.dNdx[base_func, q_point]
@inline shape_gradient(bv::FaceValues, q_point::Int, base_func::Int) = bv.dNdx[base_func, q_point, bv.current_face[]]

"""
    shape_symmetric_gradient(fe_v::Values, q_point::Int, base_function::Int)

Return the symmetric gradient of shape function `base_function` evaluated in
quadrature point `q_point`.
"""
@inline shape_symmetric_gradient(cv::CellVectorValues, q_point::Int, base_func::Int) = symmetric(shape_gradient(cv, q_point, base_func))
const shape_derivative = shape_gradient

"""
    shape_divergence(fe_v::Values, q_point::Int, base_function::Int)

Return the divergence of shape function `base_function` evaluated in
quadrature point `q_point`.
"""
@inline shape_divergence(cv::CellScalarValues, q_point::Int, base_func::Int) = sum(cv.dNdx[base_func, q_point])
@inline shape_divergence(bv::FaceScalarValues, q_point::Int, base_func::Int) = sum(bv.dNdx[base_func, q_point, bv.current_face[]])
@inline shape_divergence(cv::CellVectorValues, q_point::Int, base_func::Int) = tr(cv.dNdx[base_func, q_point])
@inline shape_divergence(bv::FaceVectorValues, q_point::Int, base_func::Int) = tr(bv.dNdx[base_func, q_point, bv.current_face[]])

function shape_curl(cv::VectorValues, q_point::Int, base_func::Int)
    return curl(shape_gradient(cv, q_point, base_func))
end
curl(∇v) = Vec{3}((∇v[3,2] - ∇v[2,3], ∇v[1,3] - ∇v[3,1], ∇v[2,1] - ∇v[1,2]))

"""
    function_value(fe_v::Values, q_point::Int, u::AbstractVector)

Compute the value of the function in a quadrature point. `u` is a vector with values
for the degrees of freedom. For a scalar valued function, `u` contains scalars.
For a vector valued function, `u` can be a vector of scalars (for use of `VectorValues`)
or `u` can be a vector of `Vec`s (for use with ScalarValues).

The value of a scalar valued function is computed as ``u(\\mathbf{x}) = \\sum\\limits_{i = 1}^n N_i (\\mathbf{x}) u_i``
where ``u_i`` are the value of ``u`` in the nodes. For a vector valued function the value is calculated as
``\\mathbf{u}(\\mathbf{x}) = \\sum\\limits_{i = 1}^n N_i (\\mathbf{x}) \\mathbf{u}_i`` where ``\\mathbf{u}_i`` are the
nodal values of ``\\mathbf{u}``.
"""
function function_value(fe_v::Values{ξdim,xdim}, q_point::Int, u::AbstractVector{T}, dof_range::UnitRange = 1:length(u)) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    isa(fe_v, VectorValues) && (n_base_funcs *= xdim)
    @assert length(dof_range) == n_base_funcs
    @boundscheck checkbounds(u, dof_range)
    val = zero(_valuetype(fe_v, u))
    @inbounds for (i, j) in enumerate(dof_range)
        val += shape_value(fe_v, q_point, i) * u[j]
    end
    return val
end

function function_value(fe_v::VectorValues{ξdim,xdim}, q_point::Int, u::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    @assert length(u) == n_base_funcs
    val = zero(Vec{xdim, T})
    basefunc = 1
    @inbounds for i in 1:n_base_funcs
        for j in 1:xdim
            val += shape_value(fe_v, q_point, basefunc) * u[i][j]
            basefunc += 1
        end
    end
    return val
end

Base.@pure _valuetype(::ScalarValues{dim}, ::AbstractVector{T}) where {dim,T} = T
Base.@pure _valuetype(::ScalarValues{dim}, ::AbstractVector{Vec{dim,T}}) where {dim,T} = Vec{dim,T}
Base.@pure _valuetype(::VectorValues{dim}, ::AbstractVector{T}) where {dim,T} = Vec{dim,T}
# Base.@pure _valuetype(::VectorValues{dim}, ::AbstractVector{Vec{dim,T}}) where {dim,T} = Vec{dim,T}

"""
    function_scalar_gradient(fe_v::Values{ξdim,xdim}, q_point::Int, u::AbstractVector)

Compute the gradient of the function in a quadrature point. `u` is a vector with values
for the degrees of freedom. For a scalar valued function, `u` contains scalars.
For a vector valued function, `u` can be a vector of scalars (for use of `VectorValues`)
or `u` can be a vector of `Vec`s (for use with ScalarValues).

The gradient of a scalar function is computed as
``\\mathbf{\\nabla} u(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x}) u_i``
where ``u_i`` are the nodal values of the function.
For a vector valued function the gradient is computed as
``\\mathbf{\\nabla} \\mathbf{u}(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x}) \\otimes \\mathbf{u}_i``
where ``\\mathbf{u}_i`` are the nodal values of ``\\mathbf{u}``.
"""
function function_gradient(fe_v::Values{ξdim,xdim}, q_point::Int, u::AbstractVector{T}, dof_range::UnitRange = 1:length(u)) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    isa(fe_v, VectorValues) && (n_base_funcs *= xdim)
    @assert length(dof_range) == n_base_funcs
    @boundscheck checkbounds(u, dof_range)
    grad = zero(_gradienttype(fe_v, u))
    @inbounds for (i, j) in enumerate(dof_range)
        grad += shape_gradient(fe_v, q_point, i) * u[j]
    end
    return grad
end

Base.@pure _gradienttype(::ScalarValues{ξdim,xdim}, ::AbstractVector{T}) where {ξdim,xdim,T} = Vec{xdim,T}
Base.@pure _gradienttype(::VectorValues{ξdim,xdim}, ::AbstractVector{T}) where {ξdim,xdim,T} = Tensor{2,xdim,T}

function function_gradient(fe_v::ScalarValues{ξdim,xdim}, q_point::Int, u::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    @assert length(u) == n_base_funcs
    grad = zero(Tensor{2,xdim,T})
    @inbounds for i in 1:n_base_funcs
        grad += u[i] ⊗ shape_gradient(fe_v, q_point, i)
    end
    return grad
end

function function_gradient(fe_v::VectorValues{ξdim,xdim}, q_point::Int, u::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    @assert length(u) == n_base_funcs
    grad = zero(Tensor{2,xdim,T})
    basefunc_count = 1
    @inbounds for i in 1:n_base_funcs
        for j in 1:xdim
            grad += u[i][j] * shape_gradient(fe_v, q_point, basefunc_count)
            basefunc_count += 1
        end
    end
    return grad
end

const function_derivative = function_gradient

"""
    function_symmetric_gradient(fe_v::Values, q_point::Int, u::AbstractVector)

Compute the symmetric gradient of the function, see [`function_gradient`](@ref).
Return a `SymmetricTensor`.

The symmetric gradient of a scalar function is computed as
``\\left[ \\mathbf{\\nabla}  \\mathbf{u}(\\mathbf{x_q}) \\right]^\\text{sym} =  \\sum\\limits_{i = 1}^n  \\frac{1}{2} \\left[ \\mathbf{\\nabla} N_i (\\mathbf{x}_q) \\otimes \\mathbf{u}_i + \\mathbf{u}_i  \\otimes  \\mathbf{\\nabla} N_i (\\mathbf{x}_q) \\right]``
where ``\\mathbf{u}_i`` are the nodal values of the function.
"""
function function_symmetric_gradient(fe_v::Values, q_point::Int, u::AbstractVector, dof_range::UnitRange = 1:length(u))
    grad = function_gradient(fe_v, q_point, u, dof_range)
    return symmetric(grad)
end

function function_symmetric_gradient(fe_v::Values, q_point::Int, u::AbstractVector{Vec{xdim, T}}) where {xdim, T}
    grad = function_gradient(fe_v, q_point, u)
    return symmetric(grad)
end

"""
    function_divergence(fe_v::Values, q_point::Int, u::AbstractVector)

Compute the divergence of the vector valued function in a quadrature point.

The divergence of a vector valued functions in the quadrature point ``\\mathbf{x}_q)`` is computed as
``\\mathbf{\\nabla} \\cdot \\mathbf{u}(\\mathbf{x_q}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x_q}) \\cdot \\mathbf{u}_i``
where ``\\mathbf{u}_i`` are the nodal values of the function.
"""
function function_divergence(fe_v::ScalarValues{ξdim,xdim}, q_point::Int, u::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    @assert length(u) == n_base_funcs
    diverg = zero(T)
    @inbounds for i in 1:n_base_funcs
        diverg += shape_gradient(fe_v, q_point, i) ⋅ u[i]
    end
    return diverg
end

function function_divergence(fe_v::VectorValues{ξdim,xdim}, q_point::Int, u::AbstractVector{T}, dof_range::UnitRange = 1:length(u)) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)*xdim
    @assert length(dof_range) == n_base_funcs
    @boundscheck checkbounds(u, dof_range)
    diverg = zero(T)
    @inbounds for (i, j) in enumerate(dof_range)
        grad = shape_gradient(fe_v, q_point, i)
        for k in 1:xdim
            diverg += grad[k, k] * u[j]
        end
    end
    return diverg
end

function function_divergence(fe_v::VectorValues{ξdim,xdim}, q_point::Int, u::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getn_scalarbasefunctions(fe_v)
    @assert length(u) == n_base_funcs
    diverg = zero(T)
    basefunc_count = 1
    @inbounds for i in 1:n_base_funcs
        for j in 1:xdim
            grad = shape_gradient(fe_v, q_point, basefunc_count)
            basefunc_count += 1
            for k in 1:xdim
                diverg += grad[k, k] * u[i][j]
            end
        end
    end
    return diverg
end

function function_curl(fe_v::Values, q_point::Int, u::AbstractVector, dof_range::UnitRange = 1:length(u))
    return curl(function_gradient(fe_v, q_point, u, dof_range))
end

function function_curl(fe_v::Values, q_point::Int, u::AbstractVector{Vec{3, T}}) where T
    return curl(function_gradient(fe_v, q_point, u))
end

"""
    spatial_coordinate(fe_v::Values{ξdim,xdim}, q_point::Int, x::AbstractVector)

Compute the spatial coordinate in a quadrature point. `x` contains the nodal
coordinates of the cell.

The coordinate is computed, using the geometric interpolation, as
``\\mathbf{x} = \\sum\\limits_{i = 1}^n M_i (\\mathbf{x}) \\mathbf{\\hat{x}}_i``
"""
function spatial_coordinate(fe_v::Values{ξdim,xdim}, q_point::Int, x::AbstractVector{Vec{xdim,T}}) where {ξdim,xdim,T}
    n_base_funcs = getngeobasefunctions(fe_v)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{xdim,T})
    @inbounds for i in 1:n_base_funcs
        vec += geometric_value(fe_v, q_point, i) * x[i]
    end
    return vec
end

function Base.show(io::IO, ::MIME"text/plain", fe_v::Values)
    print(io, "$(typeof(fe_v)) with $(getnbasefunctions(fe_v)) shape functions and $(getnquadpoints(fe_v)) quadrature points")
end

# copy
for ValueType in (CellScalarValues, CellVectorValues, FaceScalarValues, FaceVectorValues)
    args = [:(copy(cv.$fname)) for fname in fieldnames(ValueType)]
    @eval begin
        function Base.copy(cv::$ValueType)
            return typeof(cv)($(args...))
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 48c789c8-214e-11eb-0c42-c90206ac9147
# Boilerplate code to pull in required packages in a portable way
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using Plots
	using ERFA
	using LinearAlgebra # for dot product

md"# Phased Array Beam Patterns

This notebook develops a beam pattern analysis for phased arrays of homogenous radio telescopes.
"
end

# ╔═╡ fa66394c-21a1-11eb-1991-fb0fb1181350
md"""## Introduction"""

# ╔═╡ dc0d7392-2175-11eb-0964-5b252639b5fd
md"""For steerable radio telecscope arrays, the array's _control system_ directs the antennas' drive systems to orient the antennas' parabolic reflectors such that the center of each of antenna's primary beam is aligned with the location specified by the observer's instructions. This location is referred to as the _pointing center_. Typically, the pointing center is continuously updated to _track_ the apparent location of a celestial source (i.e. object or location) as it moves across the sky.

Most steerable radio telescopes have _altitude-azimuth_ (or _Alt-Az_) drive systems that are steerable in azimuth and elevation (aka altitude).  Some have equitorial (aka polar) drive systems that are steerable in hour angle and declination.  In either case, the phased array beam analysis is similar.

Non-steerable radio telescopes rely on the apparent motion of sources across the sky (e.g. due to the rotation of Earth) to allow observation of multiple pointing centers.  Many of the phased array concepts presented here are relevant for non-steerable readio telescopes as well, but non-steerable radio telescopes often have non-circular primary beam patterns which can complicate such analysis.

The anaysis presented here is for a radio telescope array whose antennas have a circular primary beam, though the concepts should apply to non-circular primary beams as well.  This analysis also assumes that the radio telescopes are steerable in azimuth and elevation, but the end result is basically the same regardless of the type of drive system.
"""

# ╔═╡ b2c21684-21a5-11eb-0d02-853dfc69164c
md"## Phased Array Theory of Operation"

# ╔═╡ ceed38a2-231f-11eb-1f9a-cf730d080095
md"""The voltage produced by each radio telescope is proportional to the sum of the radio frequency energies emitted by the all the sources within the primary beam (and sidelobes) plus the uncorrelated system noise of the antenna and signal path."""

# ╔═╡ 7f734cbc-219b-11eb-1471-594ab489e084
md"When operating in _phased array_ mode, the voltages from the antennas are aligned in both time and phase relative to a desired sky position and the array's reference position. This desired sky position is known as the _phase center_. Although not strictly required, the phase center is almost always located within the primary beam. When the phase center happens to be the same as the pointing center (i.e. the center of the primary beam), the array is said to be _phased to boresight_.
"

# ╔═╡ fd656238-2176-11eb-3a4e-7dbd106b836b
md"""The time and phase aligned signals from all the antennas are then added together to produce a _phased array sum_ (aka _coherent beam_).  This results in _constructive interference_ at the phase center where the correlated energy (i.e. the energy from celestial sources, not the energy from system noise) from the antennas add together in-phase for maximum gain while the uncorrelated energy (aka system noise) adds together randomly and the correalted energy from elsewhere in the primary beam adds together with sub-optimal phasing that varies across the primary beam. Determining how much gain each primary beam location experiences due to summing is the main objective of the phased array beam pattern analysis.

The time alignemnt depends on the _geometric delays_ from the antennas to the array reference location as the RF wavefront propgates across the array. These geometric delays differ depending on phase center location and array geometry. The signal from each antenna also experiences a finite _fixed delay_ that generally differs between antennas (e.g. due to cable length variations). The geometric delays vary as the phase center moves relative to the array whereas the fixed delays, as their name implies, do not generally vary.  Geometric delays can be calculated when the phase center and array geometry are known, but fixed delays must be measured in situ or given _a priori_.  Other delays, such as those induced by atmosphereic variations over the array, can further complicate things but this analysis will ignore these other sources of delay and delay variation.

Once the geometric and fixed delays are known and compensated for, the signals must still be phase aligned before they will add contructively.  The signal from an antenna can be delay corrected for a single delay value only.  This delay value is the geometric delay for the _delay center_ (generally the same location as the phase center) plus the antenna's fixed delay. Delay corrections are generally performed in discrete steps, which means that there is usually a small residual delay that remains uncorrected.  Phase is the product of frequency and time, ``\phi = \omega \tau``, so the residual delay (a time offset) results in a frequency dependent phase offset.  This phase offset can be removed on a per channel basis after the signal has been channelized into multiple frequency components (e.g. by an FFT or PFB). Locations other than the delay center also have their own residual delays that vary depending on their relationship to the delay center.  Additionally, each antenna has a unique instrumental phase that results from many different factors and must be measured in situ.  Determination of the instrumental phases is beyond the scope of this analysis, so we assume here that the intrumental phases are zero or have already been corrected for.
"""

# ╔═╡ e7853568-21a5-11eb-13de-f17b68612831
md"""Phasing the antennas to a specific phase center requires determination of the geometric delays.  The concept for this is relatively straightforard.  The geometric delay is proportional to the geometric distance that the planar wavefront travels from antenna to array reference location. This distance is the projection of the baseline vector ``\mathbf{b}``, defined to be the vector from the array reference position to the antenna position, onto the unit vector ``\hat{s}``, which points in the direction of the delay center from the array reference position.  This can be computed by taking the dot product of ``\mathbf{b}`` and ``\hat{s}``, which is written as ``\mathbf{b}\cdot\hat{s}``.  The array reference position is usually treated as the origin of the array reference frame so the baseline vector ``\mathbf{b}`` is the same as the antenna position in the array reference frame.  To take the dot product, both ``\mathbf{b}`` and ``\hat{s}`` must be represented in the same frame."""

# ╔═╡ 190282de-21a2-11eb-176e-178fd029cb6d
md"""## Primary Beam Coordinates"""

# ╔═╡ c2295c46-219f-11eb-2d92-cfdeba1d942d
md"""The analysis starts by exploring the coordinates of various points within the primary beam of the relescope. The primary beam is modeled as a circle with an angular diameter equal to the full width half maximum (FWHM) of a single radio telescope. For convenience, this circle is inscribed in a square region of sky which is represented by a two dimensional grid of points with the central point corresponding to the center of the primary beam, i.e. the pointing center. To keep things simple (or at least less complicated), this analysis uses a _flat sky_ approximation.
"""

# ╔═╡ 54a02260-219d-11eb-17b6-17b0e0bc9d8d
md"The axes of this grid correspond to cross-elevation (horizontal, aka `xel`) and elevation (vertical, aka `el`) as seen from the telescope's perspective. The origin of the grid is the center of the primary beam. The unit steps in the cross-elevation and elevations directions are `dxel` and `del`, respectively, though they are sometimes shortened here to `dx` and `de`. The coordinate frame of this grid is referred to here as the _primary beam frame_. A primary beam with FWHM of 2 (arbitrary angular units) and `dxel`, `del` values of 0.25 (arbitrary angular units) can depicted as:
"

# ╔═╡ bc7cb314-2150-11eb-2900-9b1c7ab0a8d9
begin
	plot(cos, sin, 0:0.01:2π, aspect_ratio=:equal, legend=false, xflip=true)
	scatter!(vec([(x,y) for y=-1:0.25:1, x=-1:0.25:1]), ms=2)
	title!("Primary Beam")
	xlabel!("Cross Elevation (dxel=0.25)")
	ylabel!("Relative Elevation (del=0.25)")
end

# ╔═╡ 2f22727e-22ef-11eb-379f-fd85d1a81d55
md"""For reasons that will be clear later, we have opted to orient the cross elevation (horizontal) axis with poitive values to the left of the origin and negative values to the right."""

# ╔═╡ a018417a-2227-11eb-2f6d-b7460d2e140d
md"""After compensating for geometric delay at the delay center, the residual delay at the delay center is zero.  The rest of the antenna's primary beam will have a delay gradient that depends on the geometry of the delay center relative to the antenna location and array reference location.  To compute the delay for a given point, we must create a unit vector, ``\hat{s}``, in the direction of the point and take the dot product of that and the baseline vector ``\mathbf{b}``.  We are free to choose the coordinate frame in which ``\hat{s}`` is represented.  Usually that frame will be different from the frame in which the antenna positions are expressed, so before taking the dot product we need to ensure we convert the antenna positions into the frame of ``\hat{s}`` or vice versa.  For beam pattern analysis we will typically have more ``\hat{s}`` vectors than antennas, so it will be more efficient to transform the antenna positions into the frame of the ``\hat{s}`` vectors."""

# ╔═╡ 9eb14aca-2323-11eb-3543-33f08ff00e2c
md"""As described above, the coordinate frame for ``\hat{s}`` is most naturally expressed in terms of spherical coordintates of ``(\theta, \phi, 1)``, where ``\theta`` is ``\sphericalangle xel``, ``\phi`` is ``\sphericalangle el``, and ``1`` is the unit length of ``\hat{s}``. The ``\hat{s}`` vector needs to be expressed in a rectalinear form so that we can take the dot product with rectalinear baseline vector ``\mathbf{b}``. The rectalinear form of ``\hat{s}`` is also known as _direction cosines_. The ERFA library we are using provides a spherical to Cartesian function, ``\mathtt{s2c}(θ,ϕ)``, that we can use to calculate ``\hat{s}`` for any values of `xel` and `el`.  Here is how to calculate ``\hat{s}`` for the pointing center and two other points (one on each axis):"""

# ╔═╡ 5b49716c-2327-11eb-2df4-37ceb3707973
begin
	dxel = 0.25 * ERFA.DD2R
	del  = 0.25 * ERFA.DD2R
	[ERFA.s2c(θ, ϕ) for (θ, ϕ) in [(0,0), (dxel,0), (0,del)]]
end

# ╔═╡ a85178f6-2327-11eb-2535-c1f23b6a5484
md"""As you can see, `s2c()` takes two angles (NB: in radians) and returns a rectalinear vector that, in our case, corresponds to `(s, xel, el)`, where `s` is the direction of the source (away from the viewer), `xel` is the cross elevation direction, and `el` is the elevation direction. We choose `xel` to be positive to the left of the origin and `el` positive above the origin to make this a right-handed system because the antenna positions are also specified in a right handed system."""

# ╔═╡ a0fd25d6-2335-11eb-223d-fb1a8be4a076
md"""## Antenna Position Coordinates"""

# ╔═╡ ac2e3c86-222b-11eb-0297-f1aa2b01b60d
md"""Antenna positions are often represented in a topocentric (East, North, Up) or ENU frame, which is a right-handed system.  Antenna positions can also be represented in an _Earth Centered Earth Fixed_ (ECEF) frame known as the _International Terrestrial Reference Frame_ (ITRF).  The ECEF frame can be converted to the ENU frame used in the analysis here, so we will not treat it separately.  Antenna positions are usually given in meters, but other units are sometimes used as well (e.g. nanoseconds or wavelengths at a given frequency).  For our purposes, we will ultimately need to convert the units to time and/or a frequency dependent phase."""

# ╔═╡ 5493615c-2336-11eb-2c36-0772eadafb39
md"""## Coordinate Frame Transformations"""

# ╔═╡ 396e99c8-2336-11eb-24a9-fddf80b160b5
md"""The (s, xel, el) primary beam frame described above is beam-centric in that it does not depend on the location of the pointing center on the sky.  On the other hand, the transformation from the topocentric ENU frame to the beam-centric (s, xel, el) frame does depend on the location of the pointing center."""

# ╔═╡ ff98a0d8-2252-11eb-1056-59c804698c62
md"""The transformation from the ENU frame to the (s, xel, el) frame can be represented as a series of two rotations around various axes of the frame iteslf. Transforming from the ENU frame to the (s, xel, el) frame can be performed by:

  1. Rotating the (E,N,U) frame around the U (Up) axis in the clockwise direction (as viewed looking down to the origin) by the topocentric azimuth angle of the pointing center minus 90 degrees.  This results in an intermediate (E',N',U) frame.
  2. Rotating the intermediate (E',N',U) frame around the N' axis in the clockwise direction (as viewed looking at the origin from positive N') by the elevation angle.  This results in a (E", N', U') frame that is equivalent to the (s, xel, el) frame in which our ``\hat{s}`` vectors are expressed.
"""

# ╔═╡ f4f6ab3a-222c-11eb-0564-17ceb3dde858
md"""Each of these rotations can be represented by a 3x3 rotation matrix. These matices can be multiplied together, in the proper order, to produce a single composite 3x3 rotation matrix that transforms the ENU frame into the beam-centric (s, xel, el) frame for the azimuth and elevation of a given the pointing center.  The ENU antenna position(s) represented as a 3 element column vector (or 3xM matrix of M antenna positions) can be left multiplied by this composite rotation matrix to produce the antenna position(s) in the (s, xel, el) frame.  Here we define an `enu2beam(az, el)` function that returns this composite rotation matrix for a pointing center at azimuth `az` and elevation `el` (both given in degrees):"""

# ╔═╡ 57172622-2325-11eb-21f5-87aba95fac55
begin
	"""
	Return a 3x3 rotation matrix for ENU to beam transformation for the
	pointing center at azimuth `az` and elevation `el` (both given in degrees).
	"""
	function enu2beam(az, el)
		ERFA.ry(-el*ERFA.DD2R,
			ERFA.rz(π/2 - az*ERFA.DD2R,
				[1.0 0.0 0.0
				 0.0 1.0 0.0
				 0.0 0.0 1.0]
			)
		)
	end

	# Version that takes a Tuple and auto-splats it
	enu2beam(azel::Tuple{Real, Real}) = enu2beam(azel...)
end

# ╔═╡ 610ca96c-2330-11eb-1869-593e53e74f42
md"""Let's try it out on some basic test cases!"""

# ╔═╡ 6096037a-2330-11eb-1e75-c360e0e0c494
enu2beam(0, 0) # north horizon

# ╔═╡ 76617054-2332-11eb-15f1-91a9f5207081
enu2beam(90, 0) # east horizon

# ╔═╡ 440f67c2-2333-11eb-0975-f362da0edde5
enu2beam(0, 90) # zenith from north

# ╔═╡ b26834ae-2333-11eb-13da-7d924b886cec
enu2beam(180, 90) # zenith from south

# ╔═╡ e902f9fe-2333-11eb-0bb5-9949d2f05225
md"""Notice that in the last example (zenith from south), the sign of `xel` and `el` will be opposite from the example that preceeds it (zenith from north) even though both have the local zenith as the pointing center."""

# ╔═╡ e8db8432-2333-11eb-3195-877803e736ba
md"""## Putting It All Together"""

# ╔═╡ e8b55f50-2333-11eb-0f6c-abc07cdfd3cc
md"""We now have the tools to start building and putting together all the various parts of the analysis:
  - Grid of points covering the primary beam
  - Unit ``\hat{s}`` vectors for each of those points
  - Antenna positions in ENU form
  - Conversion from ENU frame to beam frame for arbitrary pointing center
  - Dot product to compute distance/delay from antenna to array reference poistion
"""

# ╔═╡ d620f9aa-23ba-11eb-11a2-75d99569a381
md"""Here we create a grid of points to cover the primary beam and unit ``\hat{s}`` vectors corresponding to those points:"""

# ╔═╡ 165f5502-23b6-11eb-3824-dde609547349
npoints = 65

# ╔═╡ 5a17d3a8-23b6-11eb-278f-0f399b01b6b1
fwhm = 2.0

# ╔═╡ 11735a72-23b7-11eb-0c75-39f94e10ba04
xels = range(-fwhm/2, +fwhm/2, length=npoints)

# ╔═╡ 4697cad0-23b7-11eb-17fd-d73044fab8bb
els = range(-fwhm/2, +fwhm/2, length=npoints)

# ╔═╡ 7ba15bf6-23b7-11eb-3355-9932eae00a13
grid_coords = collect(Iterators.product(xels, els))

# ╔═╡ a121be52-23b7-11eb-3fa0-f7b7cef14f66
# We could also define an `ERFA.s2c(t::Tuple)` method and
# use broadcasting, but using `map()` is good enough for now.
grid_svecs = map(p->ERFA.s2c((p.*ERFA.DD2R)...), grid_coords)

# ╔═╡ babfe764-23ba-11eb-02a2-714e5c641733
md"""We start with one antenna which you are free to move around (up to 1 km in each direction):

$(@bind user_east html"<input type=range min=-100000 max=100000 value=3000 style='width: 75%'> East
")
$(@bind user_north html"<input type=range min=-100000 max=100000 value=0 style='width: 75%'> North
")
$(@bind user_up html"<input type=range min=-100000 max=100000 value=0 style='width: 75%'> Up")
"""

# ╔═╡ bf0e1a08-23e1-11eb-22e3-2f9770941b9f
user_antpos_enu = [user_east, user_north, user_up] ./ 100

# ╔═╡ 34d02b10-23e1-11eb-1ef6-3540d2475004
md"""And of course you get to point the antenna too:

$(@bind user_az html"<input type=range min=0 max=359 value=90 style='width: 75%'> Azimuth
")
$(@bind user_el html"<input type=range min=0 max=90 value=0 style='width: 75%'> Elevation")
"""

# ╔═╡ ba454464-23ba-11eb-08f7-bb9bd85194e8
user_pointing_azel = (user_az, user_el)

# ╔═╡ 62a0f894-23d7-11eb-347f-c3c1de1f41a8
md"""Now we can transform your specified (E, N, U) antenna position from the ENU frame to the (s, xel, el) frame of your specified (az, el) pointing center."""

# ╔═╡ 18d78b34-23d6-11eb-33f1-8df5a4066d65
user_enu2beam = enu2beam(user_pointing_azel)

# ╔═╡ e899c044-23d6-11eb-2780-a1fe38bec119
user_antpos_pointing = user_enu2beam * user_antpos_enu

# ╔═╡ ece4a6b8-23c3-11eb-38e7-bd248d0f0bc5
md"""With the antenna position in the same frame as the unit ``\hat{s}`` vectors, we can take the dot products to calculate the distances from the antenna to the array reference position (i.e. the origin of the ENU frame) for each point in the primary beam.  We are more interested in the relative distances across the beam rather than the absoulte distances.  To calculate the relative distances we subtract each unit vector ``\hat{s}`` from unit vector ``\hat{d}`` that points in the direction of the delay center (assumed here to be the same as the pointing center).  These relative distances can be plotted as a contour plot.
"""

# ╔═╡ a2bf4b28-23d5-11eb-1ef1-1170b39fcb9e
reldists = map(s->dot(ERFA.s2c(0,0)-s, user_antpos_pointing), grid_svecs);

# ╔═╡ 601ffcd6-23c4-11eb-2012-01bb72855b38
contour(xels, els, reldists',
	aspect_ratio=:equal, fill=true, clabels=true,
	fillcolor=:lightrainbow, linecolor=:black,
	title="Antenna to Origin: Distance\nrelative to (az=$(user_az), el=$(user_el))",
	xlabel="Cross Elevation (degrees)",
	ylabel="Elevation (degrees)"
)

# ╔═╡ a37ac5d2-2483-11eb-3761-c18b249d25c1
md"""For a given frequency, these relative distances can be converted to relative phases which can also be plotted as a heapmap plot:

$(@bind user_mhz html"<input type=range min=500 max=2000 value=1420 style='width: 75%'> Frequency (MHz)")
"""

# ╔═╡ af50e5c2-2487-11eb-1da9-add2857380a8
md"""Frequency: $(user_mhz) MHz"""

# ╔═╡ 13c703ec-2488-11eb-06e2-a93d6fc4ada1
begin
	# ERFA.CMPS is speed of light in meters per second
	λ = ERFA.CMPS / user_mhz / 1e6
	relphases = rem2pi.(2π .* reldists ./ λ, RoundNearest) .* ERFA.DR2D
end;

# ╔═╡ de5be0d6-2489-11eb-21b1-e3218333e30d
begin
	plt=heatmap(xels, els, relphases',
		aspect_ratio=:equal, color=:hsv,
		clims=(-180,180), levels=71,
		title="Antenna to Origin: Degrees of Phase\nrelative to az=$(user_az), el=$(user_el) at $(user_mhz) MHz (λ=$(round(1e5λ)/1e3) cm)\n",
		xlabel="Cross Elevation (degrees)",
		ylabel="Elevation (degrees)"
	)
	# Plot 5 degree contours if pases are within +/- 90 degrees
	# (contouring around phase wraps gets ugly)
	if maximum(abs.(relphases)) <= 90
		contour!(xels, els, relphases',
			aspect_ratio=:equal, fill=false, clabels=true,
			levels=71, linecolor=:black,
			colorbar_entry=false
		)
	end
	plt
end

# ╔═╡ 7324d316-24ae-11eb-2112-a51606afab74
if maximum(abs.(relphases)) <= 90
	md"""(Phases within ±90 degrees: Showing 5 degree contour lines.)"""
else
	md"""(Phases greater than ±90 degrees: Contour lines are not being shown.)"""
end

# ╔═╡ 22338236-24b3-11eb-02aa-0bc76375dfd3
md"""## Adding a Second Antenna"""

# ╔═╡ 4196dd3a-24b3-11eb-2e0d-650a39dc4bb6
md"""Now that we can compute phases across the primary beam of an antenna relative to the array center for a given pointing/delay center at a given frequency, the next thing to do is add another antenna.  With two antennas, we can delay correct them for a common delay/pointing center, calcluate their phases across the primary beam, convert these phases to unit vectors in the complex plane (i.e. complex numbers), and add corresponding complex numbers together to create the complex beam pattern of the coeherent beam created by adding the two delay corrected signals together."""

# ╔═╡ aef63f92-24b3-11eb-1a6c-79b47130a4a9
md"""To keep the UI simple, we will use new antennas, new pointing center, and new frequency.  The new antennas will have an `Up` coordinate of zero, but you can still move them East/West and North/South."""

# ╔═╡ dc9d05c2-24b8-11eb-0ec2-87d8c461845a
md"""
$(@bind ant1_east html"<input type=range min=-100000 max=100000 value=3000 style='width: 75%'> Ant 1 East
")
$(@bind ant1_north html"<input type=range min=-100000 max=100000 value=0 style='width: 75%'> Ant 1 North
")
$(@bind ant2_east html"<input type=range min=-100000 max=100000 value=0 style='width: 75%'> Ant 2 East
")
$(@bind ant2_north html"<input type=range min=-100000 max=100000 value=3000 style='width: 75%'> Ant 2 North
")
$(@bind ant12_az html"<input type=range min=0 max=359 value=90 style='width: 75%'> Ants 1,2 Azimuth
")
$(@bind ant12_el html"<input type=range min=0 max=90 value=0 style='width: 75%'> Ants 1,2 Elevation")
$(@bind ant12_mhz html"<input type=range min=500 max=2000 value=1420 style='width: 75%'> Ants 1,2 Frequency")
"""

# ╔═╡ 15fea64e-24ba-11eb-09cf-4d47980c7fb5
begin
	ant12_enu2beam = enu2beam(ant12_az, ant12_el)

	ant1_antpos_pointing = ant12_enu2beam * [ant1_east, ant1_north, 0] ./ 100
	ant2_antpos_pointing = ant12_enu2beam * [ant2_east, ant2_north, 0] ./ 100

	ant1_reldists = map(s->dot(ERFA.s2c(0,0)-s, ant1_antpos_pointing), grid_svecs)
	ant2_reldists = map(s->dot(ERFA.s2c(0,0)-s, ant2_antpos_pointing), grid_svecs)

	λ12 = ERFA.CMPS / ant12_mhz / 1e6
	ant1_relphases = rem2pi.(2π .* ant1_reldists ./ λ12, RoundNearest)
	ant2_relphases = rem2pi.(2π .* ant2_reldists ./ λ12, RoundNearest)
	#ant1_relphases = 2π .* ant1_reldists ./ λ12
	#ant2_relphases = 2π .* ant2_reldists ./ λ12

	ant1z = exp.(1im*ant1_relphases)
	ant2z = exp.(1im*ant2_relphases)
	ant12z = ant1z .+ ant2z

	p1=heatmap(xels, els, ant1_relphases' .* ERFA.DR2D,
		aspect_ratio=:equal, color=:hsv,
		clims=(-180,180), levels=71,
		title="Ant 1 Phases",
		#xlabel="Cross Elevation (degrees)",
		#ylabel="Elevation (degrees)"
	)

	p2=heatmap(xels, els, ant2_relphases' .* ERFA.DR2D,
		aspect_ratio=:equal, color=:hsv,
		clims=(-180,180), levels=71,
		title="Ant 2 Phases",
		#xlabel="Cross Elevation (degrees)",
		#ylabel="Elevation (degrees)"
	)

	p3=heatmap(xels, els, abs.(ant12z)',
		aspect_ratio=:equal, fill=true, fillcolor=:thermal,
		climes=(0,2), levels=9,
		title="Coherent Beam (amplitude)"
	)

	p4=heatmap(xels, els, angle.(ant12z)' .* ERFA.DR2D,
		aspect_ratio=:equal, color=:hsv,
		clims=(-180,180), levels=71,
		title="Coherent Beam (phase)"
	)

	plot(p1, p2, p3, p4)
end

# ╔═╡ 63963368-24bf-11eb-1237-773f333dcb96
ant1_enu = [ant1_east, ant1_north, 0] ./ 100

# ╔═╡ d853dd02-24bf-11eb-11f2-6d197e32c0c3
ant2_enu = [ant2_east, ant2_north, 0] ./ 100

# ╔═╡ e41097de-24bf-11eb-2fed-a527d5c53dbe
ant12_azel = (ant12_az, ant12_el)

# ╔═╡ 07431902-24c0-11eb-2f8c-07c42874be62
ant12_freq_mhz = ant12_mhz

# ╔═╡ b9e16b1a-23ba-11eb-08ad-8fc21f98dca9
md"""---"""

# ╔═╡ 3beff902-225b-11eb-267d-21ef072bf855
ant1_reldists

# ╔═╡ 428408b8-2151-11eb-0530-d1996a93b332

# ╔═╡ 3e4354d4-2151-11eb-384f-7706d85ce220

# ╔═╡ bccfe344-214c-11eb-2a1a-ff912a978a1e

# ╔═╡ Cell order:
# ╟─48c789c8-214e-11eb-0c42-c90206ac9147
# ╟─fa66394c-21a1-11eb-1991-fb0fb1181350
# ╟─dc0d7392-2175-11eb-0964-5b252639b5fd
# ╟─b2c21684-21a5-11eb-0d02-853dfc69164c
# ╟─ceed38a2-231f-11eb-1f9a-cf730d080095
# ╟─7f734cbc-219b-11eb-1471-594ab489e084
# ╟─fd656238-2176-11eb-3a4e-7dbd106b836b
# ╟─e7853568-21a5-11eb-13de-f17b68612831
# ╟─190282de-21a2-11eb-176e-178fd029cb6d
# ╟─c2295c46-219f-11eb-2d92-cfdeba1d942d
# ╟─54a02260-219d-11eb-17b6-17b0e0bc9d8d
# ╟─bc7cb314-2150-11eb-2900-9b1c7ab0a8d9
# ╟─2f22727e-22ef-11eb-379f-fd85d1a81d55
# ╟─a018417a-2227-11eb-2f6d-b7460d2e140d
# ╟─9eb14aca-2323-11eb-3543-33f08ff00e2c
# ╠═5b49716c-2327-11eb-2df4-37ceb3707973
# ╟─a85178f6-2327-11eb-2535-c1f23b6a5484
# ╟─a0fd25d6-2335-11eb-223d-fb1a8be4a076
# ╟─ac2e3c86-222b-11eb-0297-f1aa2b01b60d
# ╟─5493615c-2336-11eb-2c36-0772eadafb39
# ╟─396e99c8-2336-11eb-24a9-fddf80b160b5
# ╟─ff98a0d8-2252-11eb-1056-59c804698c62
# ╟─f4f6ab3a-222c-11eb-0564-17ceb3dde858
# ╠═57172622-2325-11eb-21f5-87aba95fac55
# ╟─610ca96c-2330-11eb-1869-593e53e74f42
# ╠═6096037a-2330-11eb-1e75-c360e0e0c494
# ╠═76617054-2332-11eb-15f1-91a9f5207081
# ╠═440f67c2-2333-11eb-0975-f362da0edde5
# ╠═b26834ae-2333-11eb-13da-7d924b886cec
# ╟─e902f9fe-2333-11eb-0bb5-9949d2f05225
# ╟─e8db8432-2333-11eb-3195-877803e736ba
# ╟─e8b55f50-2333-11eb-0f6c-abc07cdfd3cc
# ╟─d620f9aa-23ba-11eb-11a2-75d99569a381
# ╟─165f5502-23b6-11eb-3824-dde609547349
# ╟─5a17d3a8-23b6-11eb-278f-0f399b01b6b1
# ╠═11735a72-23b7-11eb-0c75-39f94e10ba04
# ╠═4697cad0-23b7-11eb-17fd-d73044fab8bb
# ╠═7ba15bf6-23b7-11eb-3355-9932eae00a13
# ╠═a121be52-23b7-11eb-3fa0-f7b7cef14f66
# ╟─babfe764-23ba-11eb-02a2-714e5c641733
# ╟─bf0e1a08-23e1-11eb-22e3-2f9770941b9f
# ╟─34d02b10-23e1-11eb-1ef6-3540d2475004
# ╟─ba454464-23ba-11eb-08f7-bb9bd85194e8
# ╟─62a0f894-23d7-11eb-347f-c3c1de1f41a8
# ╠═18d78b34-23d6-11eb-33f1-8df5a4066d65
# ╠═e899c044-23d6-11eb-2780-a1fe38bec119
# ╟─ece4a6b8-23c3-11eb-38e7-bd248d0f0bc5
# ╠═a2bf4b28-23d5-11eb-1ef1-1170b39fcb9e
# ╟─601ffcd6-23c4-11eb-2012-01bb72855b38
# ╟─a37ac5d2-2483-11eb-3761-c18b249d25c1
# ╟─af50e5c2-2487-11eb-1da9-add2857380a8
# ╠═13c703ec-2488-11eb-06e2-a93d6fc4ada1
# ╠═de5be0d6-2489-11eb-21b1-e3218333e30d
# ╟─7324d316-24ae-11eb-2112-a51606afab74
# ╟─22338236-24b3-11eb-02aa-0bc76375dfd3
# ╟─4196dd3a-24b3-11eb-2e0d-650a39dc4bb6
# ╟─aef63f92-24b3-11eb-1a6c-79b47130a4a9
# ╟─dc9d05c2-24b8-11eb-0ec2-87d8c461845a
# ╟─15fea64e-24ba-11eb-09cf-4d47980c7fb5
# ╟─63963368-24bf-11eb-1237-773f333dcb96
# ╟─d853dd02-24bf-11eb-11f2-6d197e32c0c3
# ╟─e41097de-24bf-11eb-2fed-a527d5c53dbe
# ╟─07431902-24c0-11eb-2f8c-07c42874be62
# ╠═b9e16b1a-23ba-11eb-08ad-8fc21f98dca9
# ╠═3beff902-225b-11eb-267d-21ef072bf855
# ╠═428408b8-2151-11eb-0530-d1996a93b332
# ╠═3e4354d4-2151-11eb-384f-7706d85ce220
# ╠═bccfe344-214c-11eb-2a1a-ff912a978a1e


text:
#
# Plot recipes for trajectories!
#

@recipe function f(trajectory::Trajectory; vars=:XYZ, timerange=range(first(trajectory.solution.t); stop=last(trajectory.solution.t), length=1000))
    linewidth    --> 1.75
    title        --> "Trajectory"
    dpi          --> 150
    label        --> :none

    indices = vars isa Tuple ? vars : process_vars(vars)

    labels = Dict(
        0 => "Time" * " ($(timeunit(trajectory)))",
        1 => "X"    * " ($(lengthunit(trajectory)))",
        2 => "Y"    * " ($(lengthunit(trajectory)))",
        3 => "Z"    * " ($(lengthunit(trajectory)))",
        4 => "Ẋ"    * " ($(velocityunit(trajectory)))",
        5 => "Ẏ"    * " ($(velocityunit(trajectory)))",
        6 => "Ż"    * " ($(velocityunit(trajectory)))"
    )

    xguide --> labels[indices[1]]
    yguide --> labels[indices[2]]
    if length(indices) == 3
        zguide --> labels[indices[3]]
    end

    getdata(index) = index == 0 ? timerange : map(t -> trajectory(t; idxs=index), timerange)

    getdata.((indices...,))
end

@recipe function f(trajectory::AbstractVector{<:Orbit}; vars=:XYZ)
    linewidth    --> 1.75
    title        --> "Trajectory"
    dpi          --> 150
    label        --> :none

    indices = vars isa Tuple ? vars : process_vars(vars)

    labels = Dict(
        0 => "Time" * " ($(timeunit(first(trajectory))))",
        1 => "X"    * " ($(lengthunit(first(trajectory))))",
        2 => "Y"    * " ($(lengthunit(first(trajectory))))",
        3 => "Z"    * " ($(lengthunit(first(trajectory))))",
        4 => "Ẋ"    * " ($(velocityunit(first(trajectory))))",
        5 => "Ẏ"    * " ($(velocityunit(first(trajectory))))",
        6 => "Ż"    * " ($(velocityunit(first(trajectory))))"
    )

    xguide --> labels[indices[1]]
    yguide --> labels[indices[2]]
    if length(indices) == 3
        zguide --> labels[indices[3]]
    end

    getdata(index) = index == 0 ? map(epoch, trajectory) : map(orbit -> state(orbit)[index], trajectory)

    getdata.((indices...,))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Handler defines a specific action that should be taken

    IoK8sApiCoreV1Handler(;
        exec=nothing,
        httpGet=nothing,
        tcpSocket=nothing,
    )

    - exec::IoK8sApiCoreV1ExecAction : One and only one of the following should be specified. Exec specifies the action to take.
    - httpGet::IoK8sApiCoreV1HTTPGetAction : HTTPGet specifies the http request to perform.
    - tcpSocket::IoK8sApiCoreV1TCPSocketAction : TCPSocket specifies an action involving a TCP port. TCP hooks not yet supported
"""
mutable struct IoK8sApiCoreV1Handler <: SwaggerModel
    exec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ExecAction } # spec name: exec
    httpGet::Any # spec type: Union{ Nothing, IoK8sApiCoreV1HTTPGetAction } # spec name: httpGet
    tcpSocket::Any # spec type: Union{ Nothing, IoK8sApiCoreV1TCPSocketAction } # spec name: tcpSocket

    function IoK8sApiCoreV1Handler(;exec=nothing, httpGet=nothing, tcpSocket=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Handler, Symbol("exec"), exec)
        setfield!(o, Symbol("exec"), exec)
        validate_property(IoK8sApiCoreV1Handler, Symbol("httpGet"), httpGet)
        setfield!(o, Symbol("httpGet"), httpGet)
        validate_property(IoK8sApiCoreV1Handler, Symbol("tcpSocket"), tcpSocket)
        setfield!(o, Symbol("tcpSocket"), tcpSocket)
        o
    end
end # type IoK8sApiCoreV1Handler

const _property_map_IoK8sApiCoreV1Handler = Dict{Symbol,Symbol}(Symbol("exec")=>Symbol("exec"), Symbol("httpGet")=>Symbol("httpGet"), Symbol("tcpSocket")=>Symbol("tcpSocket"))
const _property_types_IoK8sApiCoreV1Handler = Dict{Symbol,String}(Symbol("exec")=>"IoK8sApiCoreV1ExecAction", Symbol("httpGet")=>"IoK8sApiCoreV1HTTPGetAction", Symbol("tcpSocket")=>"IoK8sApiCoreV1TCPSocketAction")
Base.propertynames(::Type{ IoK8sApiCoreV1Handler }) = collect(keys(_property_map_IoK8sApiCoreV1Handler))
Swagger.property_type(::Type{ IoK8sApiCoreV1Handler }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Handler[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Handler }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Handler[property_name]

function check_required(o::IoK8sApiCoreV1Handler)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Handler }, name::Symbol, val)
end


text:
# A dictionary of speed limits for interpreting the `maxspeedword`
maxspeedword = Dict(
    :no => 130,
    :none => 130,
    :motorway => 130,
    :signals => nothing, # nothing is replased by highway
    :sign => nothing,    # nothing is replased by highway
    :rural => 100,
    :trunk => 100,
    :urban => 50,
    :bicycle_road => 30,
    :walk => 7,
    :living_street => 7
)

# A dictionary of speed limits for `highway` types (in this case for germany)
highwaymaxspeed = Dict(
    :rural => Dict(
        :motorway => 130,
        :link => 80,
        :trunk => 100,
        :primary => 100,
        :secondary => 100,
        :tertiary => 100,
        :unclassified => 100,
        :residential => 100,
        :living_street => 7,
        :service => 30
    ),
    :urban => Dict(
        :motorway => 130,
        :link => 80,
        :trunk => 50,
        :primary => 50,
        :secondary => 50,
        :tertiary => 50,
        :unclassified => 50,
        :residential => 50,
        :living_street => 7,
        :service => 30
    )
)

"""for now we ignore special days
hour per week granularity
"""
function maxspeed(
        way,
        maxspeedword=maxspeedword,
        urban=true,
        highwaymaxspeed=highwaymaxspeed
    )::Dict{Symbol, Union{AbstractArray, Real}}
    if way["type"] != "way"
        throw(ArgumentError("`way` must be an OSM way object."))
    end

    tags = get(way, "tags", "_")
    if tags != "_"
        maxspeed = parsemaxspeed(
            tags,
            maxspeedword,
            "both",
            urban,
            highwaymaxspeed
        )
        maxspeedtagconditional = get(tags, "maxspeed:conditional", "_")
        if maxspeedtagconditional != "_"
            if maxspeed == nothing
                maxspeed = missing
            end
            maxspeed = parsemaxspeedconditional(maxspeedtagconditional, maxspeed)
            if length(unique(maxspeed)) == 1
                maxspeed = maxspeed[1]
            end
        end
        uvmaxspeed = parsemaxspeed(
            tags,
            maxspeedword,
            "forward",
            urban,
            highwaymaxspeed
        )
        vumaxspeed = parsemaxspeed(
            tags,
            maxspeedword,
            "backward",
            urban,
            highwaymaxspeed
        )

        out = Dict{Symbol, Union{AbstractArray, Real}}()
        if uvmaxspeed != nothing
            out[:uvmaxspeed] = uvmaxspeed
        end
        if vumaxspeed != nothing
            out[:vumaxspeed] = vumaxspeed
        end
        if maxspeed != nothing
            out[:maxspeed] = maxspeed
        end
    else
        out = Dict{Symbol, Union{AbstractArray, Real}}()
        @warn "A `way` object without `tags` is discovered it is ignored." ID=way["id"]
    end

    return out
end

""" Parse tags which give the maxspeed of a way
"""
function parsemaxspeed(
        tags,
        maxspeedword=maxspeedword,
        direction="both",
        urban=true,
        highwaymaxspeed=highwaymaxspeed
    )::Union{Real, Nothing}
    # if direction is given check for directional maxspeed
    if direction == "both"
        direction = ""
    else
        direction = ":" * direction
    end

    for tag = ["maxspeed", "source:maxspeed", "zone:maxspeed"]
        maxspeedtag = get(tags, tag * direction, "_")
        if maxspeedtag != "_"
            maxspeed = maxspeednumber(maxspeedtag, maxspeedword)
            if maxspeed != nothing

                return maxspeed
            end
        end
    end

    highwaytag = get(tags, "highway", "_")
    if highwaytag != "_" && length(direction) == 0
        maxspeed = highwaymaxspeenumber(highwaytag, urban, highwaymaxspeed)
        if maxspeed != nothing

            return maxspeed
        end
    end

    return nothing
end

""" Parse conditional tags which give the maxspeed of a way
"""
function parsemaxspeedconditional(
        maxspeedconditional::AbstractString,
        maxspeed::Union{Real,Missing}=missing,
        maxspeedword::Dict=maxspeedword
    )
    conditions = parseconditionaltag(maxspeedconditional)
    condikeys = collect(keys(conditions))

    maxspeed = repeat(Union{Real,Missing}[maxspeed], 168)
    for c = condikeys
        condi = conditions[c]
        value = maxspeednumber(condi[:value], maxspeedword)
        hours = maxspeedconditionhoursofweek(condi[:rules])
        maxspeed[hours .+ 1] .=  value
    end

    return maxspeed
end

""" Maxspeed tag as number
"""
function maxspeednumber(
        tag::AbstractString,
        maxspeedword=maxspeedword
    )::Union{Real, Nothing}
    regexmaxspeednumber = r"(\d+)"
    number = match(regexmaxspeednumber, tag)
    if number != nothing
        number = parse(Int, number.match)
        regexmaxspeedunit = r"(km/h|kmh|kph|mph|knots)$"
        unit = match(regexmaxspeedunit, tag)
        if unit === nothing

            return number
        else
            if unit.match  === nothing

                return number
            elseif unit.match in ["km/h", "kmh", "kph"]

                return number
            elseif unit.match  == "mph"

                return number * 1.609
            elseif unit.match  == "knots"

                return number * 1.852
            end
        end
    else
        wordkeys = collect(keys(maxspeedword))
        wordkeysstr = string.(wordkeys)

        regexmaxspeedword = Regex("(" * join(wordkeysstr, "|") * ")")
        word = match(regexmaxspeedword, tag)

        # if no matching word return nothing
        if word == nothing

            return nothing
        else
            ki = findall(word.match .== wordkeysstr)[1]
            k = wordkeys[ki]

            return maxspeedword[k]
        end
    end
end

""" Use the highway tag to gess the maxspeed
"""
function highwaymaxspeenumber(
        tag::AbstractString,
        urban=true,
        highwaymaxspeed=highwaymaxspeed
    )::Union{Real, Nothing}
    condkey = urban ? :urban : :rural

    wordkeys = collect(keys(highwaymaxspeed[condkey]))
    wordkeysstr = string.(wordkeys)

    regexhighway = Regex("(" * join(wordkeysstr, "|") * ")")
    word = match(regexhighway, tag)

    if word == nothing

        return nothing
    else
        ki = findall(word.match .== wordkeysstr)[1]
        k = wordkeys[ki]

        return highwaymaxspeed[condkey][k]
    end
end

""" Join rules to hours of week
"""
function maxspeedconditionhoursofweek(rules)
    ruleskeys = collect(keys(rules))

    hours = Integer[]
    # one condition with r rules
    for r = ruleskeys
        rule = rules[r]
        if rule[:words] == ["off"]
            setdiff!(hours, rule[:hoursofweek])
        elseif rule[:words] == ["signal"] || rule[:words] == ["wet"]
            # leave `hours` as is
        elseif length(rule[:words]) > 0
            @warn "A unhandel word in a conditional tag was observed." words = rule[:words]
        else
            union!(hours, rule[:hoursofweek])
        end
    end

    return hours
end


text:
module spin

import Base.==

using ..bases, ..states, ..operators, ..operators_sparse

export SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown

"""
Basis for spin-n particles.

The basis can be created for arbitrary spinnumbers by using a rational number,
e.g. SpinBasis(3//2). The Pauli operators are defined for all possible
spinnumbers.

Arguments
---------

spinnumber
    Rational number specifying the spin quantum number.
"""
type SpinBasis <: Basis
    shape::Vector{Int}
    spinnumber::Rational{Int}
    function SpinBasis(spinnumber::Rational{Int})
        @assert den(spinnumber) == 2 || den(spinnumber) == 1
        @assert den(spinnumber) > 0
        new([num(spinnumber*2 + 1)], spinnumber)
    end
end

==(b1::SpinBasis, b2::SpinBasis) = b1.spinnumber==b2.spinnumber

"""
Pauli :math:`\\sigma_x` operator for the given SpinBasis.
"""
function sigmax(b::SpinBasis)
    N = num(b.spinnumber*2 + 1)
    diag = Complex128[complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:num(2*b.spinnumber)]
    data = spdiagm(diag, 1, N, N) + spdiagm(diag, -1, N, N)
    SparseOperator(b, data)
end

"""
Pauli :math:`\\sigma_y` operator for the given SpinBasis.
"""
function sigmay(b::SpinBasis)
    N = num(b.spinnumber*2 + 1)
    diag = Complex128[1im*complex(sqrt(real((b.spinnumber + 1)*2*a - a*(a+1)))) for a=1:num(2*b.spinnumber)]
    data = spdiagm(diag, -1, N, N) - spdiagm(diag, 1, N, N)
    SparseOperator(b, data)
end

"""
Pauli :math:`\\sigma_z` operator for the given SpinBasis.
"""
function sigmaz(b::SpinBasis)
    N = num(b.spinnumber*2 + 1)
    diag = Complex128[complex(2*m) for m=b.spinnumber:-1:-b.spinnumber]
    data = spdiagm(diag, 0, N, N)
    SparseOperator(b, data)
end

"""
Raising operator :math:`\\sigma_+` for the given SpinBasis.
"""
function sigmap(b::SpinBasis)
    N = num(b.spinnumber*2 + 1)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = Complex128[complex(sqrt(float(S - m*(m+1)))) for m=b.spinnumber-1:-1:-b.spinnumber]
    data = spdiagm(diag, 1, N, N)
    SparseOperator(b, data)
end

"""
Lowering operator :math:`\\sigma_-` for the given SpinBasis.
"""
function sigmam(b::SpinBasis)
    N = num(b.spinnumber*2 + 1)
    S = (b.spinnumber + 1)*b.spinnumber
    diag = [complex(sqrt(float(S - m*(m-1)))) for m=b.spinnumber:-1:-b.spinnumber+1]
    data = spdiagm(diag, -1, N, N)
    SparseOperator(b, data)
end

"""
Spin up state for the given SpinBasis.
"""
spinup(b::SpinBasis) = basis_ket(b, 1)

"""
Spin down state for the given SpinBasis.
"""
spindown(b::SpinBasis) = basis_ket(b, b.shape[1])

end #module


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# append x,y,z coordinates and configuration to current figure
function addcoords(x::Coord,y::Coord,Z::Array,conf::CurveConf)
    global gnuplot_state

    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end

    # check coordinates: dimensions, sizes and types
    nex = !isempty(x); ney = !isempty(y); neZ = !isempty(Z)
    # check types
    if nex
        assert(eltype(x)<:Real,"Invalid coordinates")
    end
    if ney
        assert(eltype(y)<:Real,"Invalid coordinates")
    end
    if neZ
        assert(eltype(Z)<:Real,"Invalid coordinates")
        assert(1 < ndims(Z) < 4,"Invalid coordinates")
    end
    # fill missing x, y coordinates when Z is not empty
    if neZ && !nex
        x = 1:size(Z,2)
    end
    if neZ && !ney
        y = 1:size(Z,1)
    end
    # if x,y are matrices, convert to vectors
    if isa(x,Matrix)
        s = size(x)
        if s[1] == 1 || s[2] == 1
            x = squeeze(x)
        else
            error("Invalid abscissa coordinates")
        end
    elseif isa(x,Range1) || isa(x,Range)
        x = [x]
    end
    if isa(y,Matrix)
        s = size(y)
        if s[1] == 1 || s[2] == 1
            y = squeeze(y)
        else
            error("Invalid abscissa coordinates")
        end
    elseif isa(y,Range1) || isa(y,Range)
        y = [y]
    end
    # check number of elements
    if neZ
        if conf.plotstyle == "image" || conf.plotstyle == "rgbimage"
            assert(size(Z,2) == length(x), "Wrong number of columns in Z")
            assert(size(Z,1) == length(y), "Wrong number of rows in Z")
        else
            assert(size(Z,1) == length(x), "Wrong number of columns in Z")
            assert(size(Z,2) == length(y), "Wrong number of rows in Z")
        end
    else
        assert(length(x) == length(y),
        "Abscissa and ordinate must have the same number of elements")
    end

    # check curve configuration: property names
    # TODO: figure out how to check valid color names -- gnuplot supports
    #  112 different color names.
    # check valid values of marker
    assert(validate_marker(conf.marker), "Invalid mark name specified")
    if nex && !neZ ## 2-d plot
        # check valid values of plotstyle
        assert(validate_2d_plotstyle(conf.plotstyle),
            "Invalid plotstyle specified")
    elseif nex && neZ ## 3-d plot or image
        assert(validate_3d_plotstyle(conf.plotstyle),
            "Invalid plotstyle specified")
    end

    conf = copy(conf)       # we need to dereference conf
    # append data to figure
    c = findfigure(gnuplot_state.current)
    fig = gnuplot_state.figs[c]
    if fig.isempty == true
        fig.curves[1] = CurveData(x,y,Z,conf)
        fig.isempty = false
    else
        fig.curves = [fig.curves, CurveData(x,y,Z,conf)]
    end
    gnuplot_state.figs[c] = fig
end
addcoords(y) = addcoords(1:length(y),y,[],CurveConf())
addcoords(y,c::CurveConf) = addcoords(1:length(y),y,[],c)
addcoords(x,y) = addcoords(x,y,[],CurveConf())
addcoords(x,y,c::CurveConf) = addcoords(x,y,[],c)
addcoords(x,y,Z) = addcoords(x,y,Z,CurveConf())
# X, Y data in matrix columns
function addcoords(X::Matrix,Y::Matrix,conf::CurveConf)
    for i = 1:size(X,2)
        addcoords(X[:,i],Y[:,i],[],conf)
    end
end
function addcoords(Y::Matrix,conf::CurveConf)
    x = 1:size(Y,1)
    X = zeros(size(Y))
    for i = 1:size(Y,2)
        X[:,i] = x
    end
    addcoords(X,Y,conf)
end
addcoords(X::Matrix, Y::Matrix) = addcoords(X,Y,CurveConf())
addcoords(Y::Matrix) = addcoords(Y,CurveConf())

# append error data to current set of coordinates
function adderror(yl::Coord,yh::Coord)
    global gnuplot_state
    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end
    c = findfigure(gnuplot_state.current)

    # check arguments and convert to vectors
    if isempty(yl)
        error("Invalid error data")
    else
        assert(eltype(yl)<:Real,"Invalid error data")
        if isa(yl,Matrix)
            s = size(yl)
            if s[1] == 1 || s[2] == 1
                yl = squeeze(yl)
            else
                error("Invalid error data")
            end
        elseif isa(yl,Range1) || isa(yl,Range)
            yl = [yl]
        end
    end
    if !isempty(yh)
        assert(eltype(yh)<:Real,"Invalid error data")
        if isa(yl,Matrix)
            s = size(yh)
            if s[1] == 1 || s[2] == 1
                yh = squeeze(yh)
            else
                error("Invalid error data")
            end
        elseif isa(yh,Range1) || isa(yh,Range)
            yh = [yh]
        end
    end
    # verify vector sizes -- this also implies that x,y coordinates must be
    # added to figure, before error data can be attached to it
    assert(length(gnuplot_state.figs[c].curves[end].x) == length(yl),
        "Error data vector must be of same size as abscissa")
    if !isempty(yh)
        assert(length(yh) == length(yl),
            "Error data vectors must be of same size")
    end

    # set fields in current curve
    gnuplot_state.figs[c].curves[end].ylow = yl
    gnuplot_state.figs[c].curves[end].yhigh = yh

end
adderror(ydelta) = adderror(ydelta,[])

# add axes configuration to current figure
function addconf(conf::AxesConf)
    global gnuplot_state
    # check that at least one figure has been setup
    if gnuplot_state.current == 0
        figure(1)
    end

    # Perform argument validation
    # TODO: find a way to validate the box argument
    # validate axis type
    assert(validate_axis(conf.axis),"Invalid axis type specified")

    conf = copy(conf)
    # select current plot
    c = findfigure(gnuplot_state.current)
    gnuplot_state.figs[c].conf = conf
end

# llplot() is our workhorse plotting function
function llplot()
    global gnuplot_state
    global gaston_config

    # select current plot
    c = findfigure(gnuplot_state.current)
    if c == 0
        println("No current figure")
        return
    end
    figs = gnuplot_state.figs
    config = figs[c].conf

    # setup grid
    ts = gridstring(gaston_config.grid)
    gnuplot_send(ts)
    # Build terminal setup string and send it to gnuplot
    ts = termstring(gaston_config.terminal)
    gnuplot_send(ts)

    # if figure has no data, stop here
    if isempty(figs[c].curves[1].x)
        return
    end

    # datafile filename
    filename = "$(gnuplot_state.tmpdir)figure$(gnuplot_state.current).dat"

    # Send appropriate coordinates and data to gnuplot, depending on
    # whether we are doing 2-d, 3-d or image plots.

    # 2-d plot: Z is empty or plostyle is {,rgb}image
    if isempty(figs[c].curves[1].Z) ||
        figs[c].curves[1].conf.plotstyle == "image" ||
        figs[c].curves[1].conf.plotstyle == "rgbimage"
        # create data file
        f = open(filename,"w")
        for i in figs[c].curves
            ps = i.conf.plotstyle
            if ps == "errorbars" || ps == "errorlines"
                if isempty(i.yhigh)
                    # ydelta (single error coordinate)
                    writedlm(f,[i.x i.y i.ylow],' ')
                else
                    # ylow, yhigh (double error coordinate)
                    writedlm(f,[i.x i.y i.ylow i.yhigh],' ')
                end
            elseif ps == "image"
                # output matrix
                for col = 1:length(i.x)
                    y = length(i.y)
                    for row = 1:length(i.y)
                        writedlm(f,[i.x[col] i.y[row] i.Z[y,col]],' ')
                        y = y-1
                    end
                end
            elseif ps == "rgbimage"
                # output matrix
                for col = 1:length(i.x)
                    y = length(i.y)
                    for row = 1:length(i.y)
                        writedlm(f,
                        [i.x[col] i.y[row] i.Z[y,col,1] i.Z[y,col,2] i.Z[y,col,3]],
                        ' ')
                        y = y-1
                    end
                end
            else
                writedlm(f,[i.x i.y],' ')
            end
            write(f,"\n\n")
        end
        flush(f)
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # send plot command to gnuplot
        gnuplot_send(linestr(figs[c].curves, "plot", filename))

    # 3-d plot: Z is not empty and plotstyle is not {,rgb}image
    elseif !isempty(figs[c].curves[1].Z) &&
            figs[c].curves[1].conf.plotstyle != "image" &&
            figs[c].curves[1].conf.plotstyle != "rgbimage"
        # create data file
        f = open(filename,"w")
        for i in figs[c].curves
            for row in 1:length(i.x)
                for col in 1:length(i.y)
                    writedlm(f,[i.x[row] i.y[col] i.Z[row,col]],' ')
                end
                write(f,"\n")
            end
            write(f,"\n\n")
        end
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # send command to gnuplot
        gnuplot_send(linestr(figs[c].curves, "splot",filename))
    end
    gnuplot_send("reset")
end


text:
export
    load_config!,
    current_dbconfig,

    create_database,
    enable_uuid,
    create_table,
    from_sql,
    to_sql

# #############
# # db config #
# #############

db_config_path() = joinpath(DEPOT_PATH[1], "config", "QuantumStateDB", "dbconfig.json")

const DBCONFIG = Dict{Symbol,Any}()

function load_config!(; auth_file::String=db_config_path())
    global DBCONFIG

    if isfile(auth_file)
        f = open(auth_file, "r")
        for (k, v) in JSON3.read(f)
            DBCONFIG[k] = v
        end
    else
        DBCONFIG[:user] = ENV["PGUSER"]
        DBCONFIG[:password] = ENV["PGPASSWORD"]
        DBCONFIG[:dbname] = ENV["PGDATABASE"]
    end

    return DBCONFIG
end

current_dbconfig() = DBCONFIG

to_config_string(config::Dict) = join(["$k=$v" for (k, v) in config], " ")

# #############
# # create db #
# #############

function create_database(dbname::String; dbconfig=current_dbconfig())
    connection = LibPQ.Connection(to_config_string(dbconfig))
        result = execute(connection, "CREATE DATABASE $(dbname);")
        close(result)
    close(connection)

    return dbname
end

function enable_uuid(dbname::String; dbconfig=current_dbconfig())
    dbconfig = copy(dbconfig)
    dbconfig[:dbname] = dbname
    connection = LibPQ.Connection(to_config_string(dbconfig))
        execute(connection, "CREATE EXTENSION IF NOT EXISTS \"uuid-ossp\";")
    close(connection)

    return 0
end

# ################
# # create table #
# ################

function create_table(table_name::String, sql; dbconfig=current_dbconfig())
    connection = LibPQ.Connection(to_config_string(dbconfig))
        @info "Drop table $table_name if exists!"
        execute(connection, "DROP TABLE IF EXISTS $table_name;")

        result = execute(connection, sql)
        close(result)
    close(connection)

    return table_name
end

# ###############
# # query utils #
# ###############

function from_sql(sql::String; dbconfig=current_dbconfig())
    connection = LibPQ.Connection(to_config_string(dbconfig))
        result = execute(connection, sql)
            df = DataFrame(result)
        close(result)
    close(connection)

    return df
end

# ##### table <-> df #####

function from_sql(table::DataType; dbconfig=current_dbconfig())
    connection = LibPQ.Connection(to_config_string(dbconfig))
        result = execute(connection, "SELECT * FROM $(string(table));")
            df = DataFrame(result)
        close(result)
    close(connection)

    return df
end

function to_sql(df::DataFrame, table::DataType; dbconfig=current_dbconfig())
    col_names = join(lowercase.(names(df)), ", ")
    vals = join(["\$$i" for i in 1:ncol(df)], ", ")
    connection = LibPQ.Connection(to_config_string(dbconfig))
        execute(connection, "BEGIN;")
            LibPQ.load!(
                columntable(df),
                connection,
                "INSERT INTO $(string(table)) ($col_names) VALUES ($vals);",
            )
        execute(connection, "COMMIT;")
    close(connection)

    return 0
end


text:
# Alternative constructors for cudnnTensorDescriptor and cudnnFilterDescriptor

function cudnnTensorDescriptor(   # alternative constructor from array; main one in descriptors.jl
    array;
    format::cudnnTensorFormat_t=CUDNN_TENSOR_NCHW,
    dims::Vector{Cint}=dim4(size(array),Val(format))
)
    @assert length(dims) <= CUDNN_DIM_MAX # length(dims) may not be N
    cudnnTensorDescriptor(format, cudnnDataType(eltype(array)), Cint(length(dims)), dims)
end

function cudnnFilterDescriptor(  # alternative constructor from array; main one in descriptors.jl
    array;
    format::cudnnTensorFormat_t=CUDNN_TENSOR_NCHW,
    dims::Vector{Cint}=dim4(size(array),Val(format))
)
    @assert length(dims) <= CUDNN_DIM_MAX # length(dims) may not be N
    cudnnFilterDescriptor(cudnnDataType(eltype(array)), format, Cint(length(dims)), dims)
end

# From cuDNN docs: Due to historical reasons, the minimum number of dimensions in the filter
# descriptor is three, and at most CUDNN_DIM_MAX dimensions (defined in cudnn.h = 8).
# However many operations only support 4 and 5. So we will pad dims to 4.
# Note also the order of dims reverse from Julia to cuDNN.
# RNN and multiHeadAttn do use 3D descriptors so they do not use dim4.
# Note on formats: even when using the NHWC format the dims are given in NCHW order!

dim4(s::Dims{0}, ::Val{CUDNN_TENSOR_NCHW}) = Cint[1,1,1,1]
dim4(s::Dims{0}, ::Val{CUDNN_TENSOR_NHWC}) = Cint[1,1,1,1]
dim4(s::Dims{1}, ::Val{CUDNN_TENSOR_NCHW}) = Cint[s[1],1,1,1]    # Cy -> Cy,1,1,1
dim4(s::Dims{1}, ::Val{CUDNN_TENSOR_NHWC}) = Cint[s[1],1,1,1]    # Cy -> Cy,1,1,1
dim4(s::Dims{2}, ::Val{CUDNN_TENSOR_NCHW}) = Cint[s[2],s[1],1,1] # Cx,Cy -> Cy,Cx,1,1
dim4(s::Dims{2}, ::Val{CUDNN_TENSOR_NHWC}) = Cint[s[2],s[1],1,1] # Cx,Cy -> Cy,Cx,1,1
dim4(s::Dims{3}, ::Val{CUDNN_TENSOR_NCHW}) = Cint[s[3],s[2],s[1],1] # Xn,Cx,Cy -> Cy,Cx,Xn,1
dim4(s::Dims{3}, ::Val{CUDNN_TENSOR_NHWC}) = Cint[s[3],s[1],s[2],1] # Cx,Xn,Cy -> Cy,Cx,Xn,1
dim4(s::Dims{N}, ::Val{CUDNN_TENSOR_NCHW}) where {N} = Cint[reverse(s)...] # X1,...,Xn,Cx,Cy -> Cy,Cx,Xn,...,X1
dim4(s::Dims{N}, ::Val{CUDNN_TENSOR_NHWC}) where {N} = Cint[s[N],s[1],s[N-1:-1:2]...] # Cx,X1,...,Xn,Cy -> Cy,Cx,Xn,...,X1

# If array is nothing, return nothing for descriptor
cudnnTensorDescriptor(::Nothing; o...) = nothing
cudnnFilterDescriptor(::Nothing; o...) = nothing

# In case we need to get info about a descriptor

function cudnnGetTensorDescriptor(d::cudnnTensorDescriptor)
    nbDimsRequested = CUDNN_DIM_MAX
    dataType = Ref{cudnnDataType_t}(CUDNN_DATA_FLOAT)
    nbDims = Ref{Cint}(0)
    dimA = Array{Cint}(undef, CUDNN_DIM_MAX)
    strideA = Array{Cint}(undef, CUDNN_DIM_MAX)
    cudnnGetTensorNdDescriptor(d, nbDimsRequested, dataType, nbDims, dimA, strideA)
    T = juliaDataType(dataType[])
    D = (dimA[nbDims[]:-1:1]...,)
    S = (strideA[nbDims[]:-1:1]...,)
    return T,D,S
end

function cudnnGetFilterDescriptor(d::cudnnFilterDescriptor)
    nbDimsRequested = CUDNN_DIM_MAX
    dataType = Ref{cudnnDataType_t}(CUDNN_DATA_FLOAT)
    format = Ref{cudnnTensorFormat_t}(CUDNN_TENSOR_NCHW)
    nbDims = Ref{Cint}(0)
    dimA = Array{Cint}(undef, CUDNN_DIM_MAX)
    cudnnGetFilterNdDescriptor(d, nbDimsRequested, dataType, format, nbDims, dimA)
    T = juliaDataType(dataType[])
    D = (dimA[nbDims[]:-1:1]...,)
    return T,D,format[]
end

# Deprecated
function TensorDesc(ptr::cudnnTensorDescriptor_t)
    @warn "TensorDesc is deprecated, use cudnnTensorDescriptor instead." maxlog=1
    cudnnTensorDescriptor(ptr)
end

function TensorDesc(a::DenseCuArray)
    @warn "TensorDesc is deprecated, use cudnnTensorDescriptor instead." maxlog=1
    cudnnTensorDescriptor(a)
end

function TensorDesc(T::Type, size::NTuple{N,Integer}, strides::NTuple{N,Integer} = tuple_strides(size)) where N
    @warn "TensorDesc is deprecated, use cudnnTensorDescriptor instead." maxlog=1
    cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, cudnnDataType(T), Cint(N), dim4(size,Val(CUDNN_TENSOR_NCHW)))
end

function FilterDesc(ptr::cudnnFilterDescriptor_t)
    @warn "FilterDesc is deprecated, use cudnnFilterDescriptor instead." maxlog=1
    cudnnFilterDescriptor(ptr)
end

function FilterDesc(a::DenseCuArray; format = CUDNN_TENSOR_NCHW)
    @warn "FilterDesc is deprecated, use cudnnFilterDescriptor instead." maxlog=1
    cudnnFilterDescriptor(a; format)
end

function FilterDesc(T::Type, size::Tuple; format = CUDNN_TENSOR_NCHW)
    @warn "FilterDesc is deprecated, use cudnnFilterDescriptor instead." maxlog=1
    dims = dim4(size, Val(format))
    cudnnFilterDescriptor(cudnnDataType(T), format, Cint(length(dims)), dims)
end