text:
function _test_normal_form(G, D, B, p, prec = 0)
  return _test_normal_form_isometry(G, D, p) &&
         _test_normal_form_congruence(G, D, B, p, prec) &&
         Hecke._ispadic_normal_form(D, p)
end

function _test_normal_form_isometry(M, D, p)
  return islocally_isometric(Zlattice(gram = M), Zlattice(gram = D), p)
end

function _test_normal_form_congruence(G, D, B, p, prec = 0)
  X = denominator(G) * B * G * transpose(B) - denominator(G) * D
  v = valuation(reduce(gcd, Hecke._eltseq(X)), p)
  if prec == 0
    return v >= 1
  else
    return v >= prec
  end
end

@testset "helpers" begin
  R = ResidueRing(ZZ, 9)
  @test Hecke._issquare(R(1), ZZ(3))
  R = ResidueRing(ZZ, ZZ(9))
  @test Hecke._issquare(R(1), ZZ(3))
end

@testset "NormalForm" begin
  for R in [ZZ, GF(3)]
    W = R[1;]
    V = R[0 1; 1 0]

    L = diagonal_matrix([W, V])
    @test (@inferred diagonal_matrix(Hecke.collect_small_blocks(L))) == L

    L = diagonal_matrix([W, W])
    @test diagonal_matrix(Hecke.collect_small_blocks(L)) == L

    L = diagonal_matrix([W, V, W])
    @test diagonal_matrix(Hecke.collect_small_blocks(L)) == L

    L = diagonal_matrix([W])
    @test diagonal_matrix(Hecke.collect_small_blocks(L)) == L

    L = diagonal_matrix([V])
    @test diagonal_matrix(Hecke.collect_small_blocks(L)) == L

    L = diagonal_matrix([V, W, W, V, V, W, W])
    @test diagonal_matrix(Hecke.collect_small_blocks(L)) == L

    L = zero_matrix(ZZ, 0, 0)
    @test Hecke.collect_small_blocks(L) == typeof(L)[]
  end

  G1 = matrix(FlintQQ, 4, 4, [2, 0, 0, 1,
                              0, 2, 0, 1,
                              0, 0, 4, 2,
                              1, 1, 2, 6])
  G2 = matrix(FlintQQ, 4, 4, [2, 1, 1, 0,
                              1, 2, 0, 0,
                              1, 0, 2, 0,
                              0, 0, 0, 16])

  D1, U1 = @inferred Hecke.padic_normal_form(G1, 2, prec = 30)
  @test _test_normal_form(G1, D1, U1, 2, 30)

  D2, U2 = @inferred Hecke.padic_normal_form(G2, 2, prec = 30)
  @test _test_normal_form(G2, D2, U2, 2, 30)

  DD = matrix(base_ring(D1), 4, 4, [2, 1, 0, 0,
                                    1, 2, 0, 0,
                                    0, 0, 4 + 8, 0,
                                    0, 0, 0, 16])
  @test DD == D1
  DD = matrix(base_ring(D2), 4, 4, [2, 1, 0, 0,
                                    1, 2, 0, 0,
                                    0, 0, 4 + 8, 0,
                                    0, 0, 0, 16])
  @test DD == D2

  D4 = matrix(FlintQQ,  4,  4,  [2, -1, -1, -1,
                                 -1, 2, 0, 0,
                                 -1, 0, 2, 0,
                                 -1, 0, 0, 2])

  D, B = Hecke.padic_normal_form(D4, 2)
  @test _test_normal_form(D4, D, B, 2)

  @test D == matrix(base_ring(D), 4, 4, [2, 1, 0, 0,
                                         1, 2, 0, 0,
                                         0, 0, 4, 2,
                                         0, 0, 2, 4])
  A4 = matrix(FlintQQ, 4, 4, [2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2])
  D, B = Hecke.padic_normal_form(A4, 2)
  K = base_ring(D)
  @test D == matrix(K, 4, 4, [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 2, 1, 0, 0, 1, 2])
  @test _test_normal_form(A4, D, B, 2)

  A4_extended = matrix(FlintQQ, 5, 5, [2, -1, 0, 0, -1,
                                       -1, 2, -1, 0, 0,
                                       0, -1, 2, -1, 0,
                                       0, 0, -1, 2, -1,
                                       -1, 0, 0, -1, 2])

  D, B = Hecke.padic_normal_form(A4_extended, 5)
  K = base_ring(D)
  @test D == matrix(K, 5, 5, [1, 0, 0, 0, 0,
                              0, 1, 0, 0, 0,
                              0, 0, 1, 0, 0,
                              0, 0, 0, 5, 0,
                              0, 0, 0, 0, 0])
  # isometry only works for non-degenerate lattices
  #@test _test_normal_form_isometry(A4_extended, D)
  @test  _test_normal_form_congruence(A4_extended, D, B, 5)
  #@test Hecke._ispadic_normal_form(D, 5)

  A4dual = inv(A4)
  D, B = Hecke.padic_normal_form(A4dual, 5)
  K = base_ring(D)
  @test D == matrix(K, 4, 4, [inv(K(5)), 0, 0, 0,
                              0, 1, 0, 0,
                              0, 0, 1, 0,
                              0, 0, 0, 1])
  @test _test_normal_form(A4dual, D, B, 5)

  # Corner cases
  Z = matrix(FlintQQ, 0, 0, [])
  B, D = Hecke.padic_normal_form(Z, 2)
  @test iszero(B)
  @test iszero(D)

  Z = zero_matrix(FlintQQ, 10, 10)
  B, D = Hecke.padic_normal_form(Z, 3)
  @test iszero(B)
  @test iszero(D)

  # 2-adic cases

  M = matrix(FlintQQ, 3, 3, fmpq[2, 1, 0, 1, 3, 1, 0, 1, 2])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 3, 3, fmpq[400, 0, 0, 0, 160, 0, 0, 0, 80])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 6, 6, fmpq[12544, 3584, 0, 896, 8064, 0, 3584, 1152, 896, 512, 2304, 1152, 0, 896, 18816, 5376, 3584, 8064, 896, 512, 5376, 1728, 1600, 2816, 8064, 2304, 3584, 1600, 6272, 0, 0, 1152, 8064, 2816, 0, 12544])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 4, 4, fmpq[4, 0, 2, 1, 0, 4, 1, 2, 2, 1, 5, 1, 1, 2, 1, 5])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 3, 3, fmpq[1, 0, 0, 0, 1, 0, 0, 0, 1])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 6, 6, fmpq[2, 1, 0, 1, 1, 1, 1, 2, 1, 0, 0, 0, 0, 1, 2, 0, 0, 0, 1, 0, 0, 3, 0, 0, 1, 0, 0, 0, 3, 0, 1, 0, 0, 0, 0, 3])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 10, 10, fmpq[3, 1, 0, 0, 0, -1, 1, -1, 0, 1, 1, 3, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 3, 1, 0, 1, -1, -1, -1, -1, 0, 1, 1, 3, -1, 0, -1, -1, 1, -1, 0, 0, 0, -1, 3, 1, 1, 2, 1, 2, -1, 0, 1, 0, 1, 3, 0, 1, 0, -1, 1, 1, -1, -1, 1, 0, 4, 1, -1, 2, -1, 0, -1, -1, 2, 1, 1, 4, 1, 2, 0, 0, -1, 1, 1, 0, -1, 1, 4, 1, 1, 1, -1, -1, 2, -1, 2, 2, 1, 5])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 6, 6, fmpq[4, 1, -1, -1, 1, -1, 1, 4, 0, 1, 2, 1, -1, 0, 4, -1, 2, -1, -1, 1, -1, 4, -1, 0, 1, 2, 2, -1, 4, -1, -1, 1, -1, 0, -1, 4])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 5, 5, fmpq[5, -1, -1, -1, -1, -1, 5, -1, -1, -1, -1, -1, 5, -1, -1, -1, -1, -1, 5, -1, -1, -1, -1, -1, 5])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 3, 3, fmpq[2, 1, 0, 1, 4, 1, 0, 1, 4])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 7, 7, fmpq[8, -4, 3, 4, 0, 1, 1, -4, 8, 1, 0, 4, 1, 1, 3, 1, 8, 4, 0, 1, 1, 4, 0, 4, 8, 3, 0, 4, 0, 4, 0, 3, 8, 4, 0, 1, 1, 1, 0, 4, 8, -4, 1, 1, 1, 4, 0, -4, 8])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  M = matrix(FlintQQ, 24, 24, fmpq[2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18014398509481985//18014398509481984, -1//18014398509481984, -1//18014398509481984, -1//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 1, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1//36028797018963968, -1//18014398509481984, -1//18014398509481984, -1//18014398509481984, 0, 0, 0, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -36028797018963969//36028797018963968, -36028797018963969//36028797018963968, -36028797018963969//36028797018963968, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18014398509481985//18014398509481984, 0, 0, -18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 2, 0, -1, -1, 1, 0, 0, 0, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, -1, 0, 0, 0, 0, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, 2, 1, -1, 0, 0, 0, 0, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, 2, 0, 0, 0, 0, 0, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 0, 2, 0, 0, 0, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, -1, -1, 18014398509481983//18014398509481984, 0, -36028797018963969//36028797018963968, -18014398509481983//18014398509481984, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 18014398509481983//18014398509481984, 0, 0, -18014398509481983//18014398509481984, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2, 1, -18014398509481983//18014398509481984, 0, 36028797018963969//36028797018963968, 18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 2, 0, 0, 36028797018963969//36028797018963968, 0, 0, 0, 0, 0, 1, 18014398509481985//18014398509481984, 1//36028797018963968, 0, 0, 0, 0, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 5192296858534827340300120177508381//1298074214633706907132624082305024, 649037107316853399523116512706557//324518553658426726783156020576256, -162259276829213336369980246065155//162259276829213363391578010288128, 18014398509481979//324518553658426726783156020576256, 54043195528445949//649037107316853453566312041152512, -1298074214633706708974240478003207//649037107316853453566312041152512, 324518553658426690754359001612289//162259276829213363391578010288128, -324518553658426690754359001612289//162259276829213363391578010288128, 1, -1//18014398509481984, -1//18014398509481984, -36028797018963969//36028797018963968, 18014398509481985//18014398509481984, 0, 0, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 0, 0, 0, 649037107316853399523116512706557//324518553658426726783156020576256, 2596148429267413814265248164610069//649037107316853453566312041152512, 288230376151711757//1298074214633706907132624082305024, 2596148429267413814265248164610069//1298074214633706907132624082305024, -649037107316853435551913531670527//324518553658426726783156020576256, -1298074214633706817060631534895105//649037107316853453566312041152512, 324518553658426690754359001612289//324518553658426726783156020576256, 54043195528445949//649037107316853453566312041152512, 1, -1//18014398509481984, -1//18014398509481984, -36028797018963969//36028797018963968, 0, 0, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -36028797018963969//36028797018963968, 0, 36028797018963969//36028797018963968, 36028797018963969//36028797018963968, -162259276829213336369980246065155//162259276829213363391578010288128, 288230376151711757//1298074214633706907132624082305024, 2596148429267413814265248164610069//649037107316853453566312041152512, 1298074214633707051247812158160913//1298074214633706907132624082305024, -54043195528445949//649037107316853453566312041152512, 649037107316853327465522474778629//649037107316853453566312041152512, -324518553658426690754359001612289//324518553658426726783156020576256, 324518553658426690754359001612289//324518553658426726783156020576256, 1, -1//18014398509481984, -1//18014398509481984, -36028797018963969//36028797018963968, 0, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, 18014398509481979//324518553658426726783156020576256, 2596148429267413814265248164610069//1298074214633706907132624082305024, 1298074214633707051247812158160913//1298074214633706907132624082305024, 5192296858534827340300120177508381//1298074214633706907132624082305024, -1298074214633706817060631534895105//649037107316853453566312041152512, -54043195528445949//649037107316853453566312041152512, 0, 0, 0, 18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 0, 0, 0, 54043195528445949//649037107316853453566312041152512, -649037107316853435551913531670527//324518553658426726783156020576256, -54043195528445949//649037107316853453566312041152512, -1298074214633706817060631534895105//649037107316853453566312041152512, 324518553658426690754359001612289//81129638414606681695789005144064, 324518553658426690754359001612289//162259276829213363391578010288128, -324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//324518553658426726783156020576256, 0, 18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, 0, 0, 0, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, -1298074214633706708974240478003207//649037107316853453566312041152512, -1298074214633706817060631534895105//649037107316853453566312041152512, 649037107316853327465522474778629//649037107316853453566312041152512, -54043195528445949//649037107316853453566312041152512, 324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//81129638414606681695789005144064, -324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//162259276829213363391578010288128, 0, 0, 0, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//324518553658426726783156020576256, -324518553658426690754359001612289//324518553658426726783156020576256, 0, -324518553658426690754359001612289//162259276829213363391578010288128, -324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//81129638414606681695789005144064, -324518553658426690754359001612289//162259276829213363391578010288128, 0, 0, 0, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 0, 18014398509481983//18014398509481984, 18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, -18014398509481983//18014398509481984, 18014398509481983//18014398509481984, 0, -324518553658426690754359001612289//162259276829213363391578010288128, 54043195528445949//649037107316853453566312041152512, 324518553658426690754359001612289//324518553658426726783156020576256, 0, 324518553658426690754359001612289//324518553658426726783156020576256, 324518553658426690754359001612289//162259276829213363391578010288128, -324518553658426690754359001612289//162259276829213363391578010288128, 324518553658426690754359001612289//81129638414606681695789005144064])
  D, B = Hecke.padic_normal_form(M, 2)
  @test _test_normal_form(M, D, B, 2)

  # non-dyadic tests
  M = matrix(FlintQQ, 2, 2, fmpq[2, -1, -1, 2])
  D, B = Hecke.padic_normal_form(M, 3)
  @test _test_normal_form(M, D, B, 3)

  M = matrix(FlintQQ, 2, 2, fmpq[2, -1, -1, 2])
  D, B = Hecke.padic_normal_form(M, 3)
  @test _test_normal_form(M, D, B, 3)

  M = matrix(FlintQQ, 3, 3, fmpq[2, -1, 0, -1, 2, -1, 0, -1, 2])
  D, B = Hecke.padic_normal_form(M, 3)
  @test _test_normal_form(M, D, B, 3)

  M = matrix(FlintQQ, 3, 3, fmpq[3, -1, -1, -1, 3, -1, -1, -1, 3])
  D, B = Hecke.padic_normal_form(M, 3)
  @test _test_normal_form(M, D, B, 3)

  R = ResidueRing(FlintZZ, fmpz(2)^30)
  U = matrix(R,2,2,[0,1,1,0])
  V = matrix(R,2,2,[2,1,1,2])
  W1 = matrix(R,1,1,[1])
  W3 = matrix(R,1,1,[3])
  W5 = matrix(R,1,1,[5])
  W7 = matrix(R,1,1,[7])

  p = fmpz(2)
   G = diagonal_matrix(W1,W1)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [5, 0,
                                         0, 5])

   G = diagonal_matrix(W1,W3)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [5, 0,
                                         0, 7])

   G = diagonal_matrix(W1,W5)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [5, 0,
                                        0, 1])
   G = diagonal_matrix(W1,W7)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [5, 0,
                                        0, 3])

   G = diagonal_matrix(W3,W3)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 7])

   G = diagonal_matrix(W3,W5)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 1])

   G = diagonal_matrix(W3,W7)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 3])

   G = diagonal_matrix(W5,W5)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [1, 0,
                                        0, 1])
   G = diagonal_matrix(W5,W7)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [1, 0,
                                        0, 3])
   G = diagonal_matrix(W7,W7)
   b = Hecke._relations(G,1, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [3, 0,
                                        0, 3])

   G = diagonal_matrix(V,V)
   b = Hecke._relations(G,3, p)
   @test b * G * transpose(b) == matrix(R, 4, 4, [0, 1, 0, 0,
                                        1, 0, 0, 0,
                                        0, 0, 0, 1,
                                        0, 0, 1, 0])
   G = diagonal_matrix(V,W1, W1)
   b = Hecke._relations(G,5, p)
   @test b * G * transpose(b) == matrix(R, 4, 4, [0, 1, 0, 0,
                                        1, 0, 0, 0,
                                        0, 0, 7, 0,
                                        0, 0, 0, 3])

   G = diagonal_matrix(V,W1, W5)
   b = Hecke._relations(G,5, p)
   @test b * G * transpose(b) == matrix(R, 4, 4, [0, 1, 0, 0,
                                        1, 0, 0, 0,
                                        0, 0, 3, 0,
                                        0, 0, 0, 3])
   G = diagonal_matrix(V,W3, W7)
   b = Hecke._relations(G,5, p)
   @test b * G * transpose(b) == matrix(R, 4, 4, [0, 1, 0, 0,
                                        1, 0, 0, 0,
                                        0, 0, 5, 0,
                                        0, 0, 0, 5])

   G = diagonal_matrix(W1,2*W1)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [3, 0,
                                        0, 6])

   G = diagonal_matrix(W1,2*W3)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 10])

   G = diagonal_matrix(W1,2*W5)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [3, 0,
                                        0, 14])

   G = diagonal_matrix(W1,2*W7)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 2])

   G = diagonal_matrix(W3,2*W5)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [5, 0,
                                        0, 6])

   G = diagonal_matrix(W3,2*W3)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [1, 0,
                                        0, 2])

   G = diagonal_matrix(2*W5, 4*W7)
   b = Hecke._relations(G,6, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [6, 0,
                                        0, 4])

   G = diagonal_matrix(W3, 2*V)
   b = Hecke._relations(G, 7, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [7, 0, 0,
                                        0, 0, 2,
                                        0, 2, 0])

   G = diagonal_matrix(W7, 2*V)
   b = Hecke._relations(G, 7, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [3, 0, 0,
                                        0, 0, 2,
                                        0, 2, 0])

   G = diagonal_matrix(U, 2*W1)
   b = Hecke._relations(G, 8, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [2, 1, 0,
                                        1, 2, 0,
                                        0, 0, 10])
   G = diagonal_matrix(U, 2*W5)
   b = Hecke._relations(G, 8, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [2, 1, 0,
                                        1, 2, 0,
                                        0, 0, 2])
   G = diagonal_matrix(V, 2*W1)
   b = Hecke._relations(G, 8, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [0, 1, 0,
                                        1, 0, 0,
                                        0, 0, 10])
   G = diagonal_matrix(V, 2*W7)
   b = Hecke._relations(G, 8, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [0, 1, 0,
                                        1, 0, 0,
                                        0, 0, 6])

   G = diagonal_matrix(W1, W5, 2*W5)
   b = Hecke._relations(G, 9, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [3, 0, 0,
                                        0, 3, 0,
                                        0, 0, 2])

   G = diagonal_matrix(W3, W3, 2*W5)
   b = Hecke._relations(G, 9, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [5, 0, 0,
                                        0, 1, 0,
                                        0, 0, 2])
   G = diagonal_matrix(W3, W3, 2*W1)
   b = Hecke._relations(G, 9, p)
   @test b * G * transpose(b) == matrix(R, 3, 3, [5, 0, 0,
                                        0, 1, 0,
                                        0, 0, 10])

   G = diagonal_matrix(W3, 4*W1)
   b = Hecke._relations(G, 10, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [7, 0,
                                        0, 20])

   G = diagonal_matrix(W5, 4*W5)
   b = Hecke._relations(G, 10, p)
   @test b * G * transpose(b) == matrix(R, 2, 2, [1, 0,
                                        0, 4])

   # _relations errors

   @test_throws ErrorException Hecke._relations(zero_matrix(R, 0, 0), 4, p)
   @test_throws ErrorException Hecke._relations(zero_matrix(R, 0, 0), 11, p)
   @test_throws ErrorException Hecke._relations(matrix(R, 4, 4, [1, 0, 0, 0,
                                                                 0, 1, 0, 0,
                                                                 0, 0, 1, 0,
                                                                 0, 0, 0, 3]), 5, p)
   @test_throws ErrorException Hecke._relations(matrix(R, 2, 2, [1, 0,
                                                                 0, 1]), 6, p)

   @test_throws ErrorException Hecke._normalize_twobytwo(matrix(R, 2, 2, [1, 2, 2, 1]), p)

   # _two_adic_normal_forms
   G = diagonal_matrix(2*W1,2*W1,4*V)
   D, B = Hecke._two_adic_normal_forms(G, p)
   @test D == matrix(R, 4, 4, [2, 0, 0, 0,
                               0, 10, 0, 0,
                               0, 0, 0, 4,
                               0, 0, 4, 0])
   @test B * G * transpose(B) == D

   G = diagonal_matrix(W1,2*V,2*W3,2*W5)
   D, B = Hecke._two_adic_normal_forms(G, p)
   @test D == matrix(R, 5, 5, [3, 0, 0, 0, 0,
                              0, 0, 2, 0, 0,
                              0, 2, 0, 0, 0,
                              0, 0, 0, 2, 0,
                              0, 0, 0, 0, 2])
   @test D == B * G * transpose(B)

   G = diagonal_matrix(U,2*V,2*W3,2*W5)
   D, B = Hecke._two_adic_normal_forms(G, p)
   @test D == matrix(R, 6, 6, [2, 1, 0, 0, 0, 0,
                               1, 2, 0, 0, 0, 0,
                               0, 0, 4, 2, 0, 0,
                               0, 0, 2, 4, 0, 0,
                               0, 0, 0, 0, 2, 0,
                               0, 0, 0, 0, 0, 6])
   @test D == B * G * transpose(B)

   @inferred Hecke._two_adic_normal_forms(G, p, partial = true)

   # test _ispadic_normal_form

   # not the right block structure
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 3, 3, [1, 1, 1, 1, 1, 1, 1, 1, 1]), 2)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 3, 3, [1, 1, 1, 1, 1, 1, 1, 1, 1]), 3)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 2, 2, [1, 1, 1, 1]), 3)

   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 2, 2, [2, 0, 0, 1]), 3)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 2, 2, [1, 0, 0, 3]), 5)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 2, 2, [1, 0, 0, 4]), 5)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 2, 2, [1, 1, 1, 1]), 2)
   @test !Hecke._ispadic_normal_form(matrix(FlintQQ, 1, 1, [9]), 2)

 end


text:
using SpecialFunctions
using DistributedArrays
using LinearAlgebra

mutable struct BoundaryPrecomputedCalcs
  b::Boundary
  # Number of discretized steps (values go from 0:(2N-1))
  N::Int32
  k::Complex{Float64}
  j0NiMatrix::Array{Complex{Float64},2}
  j0NeMatrix::Array{Complex{Float64},2}
  j1NiMatrix::Array{Complex{Float64},2}
  j1NeMatrix::Array{Complex{Float64},2}
  h0NiMatrix::Array{Complex{Float64},2}
  h0NeMatrix::Array{Complex{Float64},2}
  h1NiMatrix::Array{Complex{Float64},2}
  h1NeMatrix::Array{Complex{Float64},2}
  absDiff::Array{Float64,2}
end

function absDiff(b::Boundary, t::Float64, τ::Float64)
  a = boundaryPoint(b, t)
  b = boundaryPoint(b, τ)
  norm(a-b)
end

function absDiff(b::Boundary, tτMatrix::Array{Array{Float64,1},2})
  map(pair -> norm(boundaryPoint(b, pair[1])-boundaryPoint(b, pair[2])), tτMatrix)
end

# @param l index (from 1:(2N-1) )
function γ(b::Boundary, N::Int64, l::Int64)
  norm(boundaryPointDeriv(b, l*π/N))+0im
end

# Quadrature weight factor
function R(j::Int64, N::Int64)
  cosAry = Array{Float64}(undef, N-1)
  for m in 1:N-1
    cosAry[m] = 1/m * cos(m*j*π/N)
  end
  -2π/N * sum(cosAry) - (-1)^j*π/N^2
end

################
j0(x::Complex{Float64}) = besselj(0, x)
j1(x::Complex{Float64}) = besselj(1, x)
h0(x::Complex{Float64}) = hankelh1(0, x)
# this is more accurate for some small values but it didn't seem to help me
# h0(x::Complex{Float64}) = j0(x)+bessely(0, x)*im
h1(x::Complex{Float64}) = hankelh1(1, x)
# h1(x::Complex{Float64}) = j1(x)+bessely(1, x)*im
nVals(b::Boundary, t::Float64) = (nInt(b), nExt(b))
xVals(b::Boundary, t::Float64, τ::Float64) = (convert(Array{Float64,1}, boundaryPoint(b, t)), convert(Array{Float64,1}, boundaryPoint(b, τ)))
log4Term(t::Float64, τ::Float64)::Float64 = log(4*sin((t-τ)/2)^2)

################
# H functions
function H1(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    # Note: Need to check that we can use the k values at t, not τ, this is a deviation from the Heider paper
    (ni, ne) = nVals(b, t)
    (xt, xτ) = xVals(b, t, τ)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j1ni = bpc.j1NiMatrix[tIdx, τIdx]
    j1ne = bpc.j1NeMatrix[tIdx, τIdx]
    -1/(2π)*dot(normVector(b, τ), (xt-xτ)/ad)*(k*ne*j1ne-k*ni*j1ni)
  end
end

function H2(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    H(b, t, τ, k, bpc)-H1(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function H(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (xt, xτ) = xVals(b, t, τ)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h1ni = bpc.h1NiMatrix[tIdx, τIdx]
  h1ne = bpc.h1NeMatrix[tIdx, τIdx]
  im/2*dot(normVector(b, τ), (xt-xτ)/ad)*(k*ne*h1ne-k*ni*h1ni)
end

H1star(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H1(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))
H2star(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H2(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))
Hstar(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))

###################
# L functions
function L1(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  if (t == τ)
    ((k*ni)^2-(k*ne)^2) / (2π)
  else
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    -1/(2π)*((k*ne)^2*j0ne-(k*ni)^2*j0ni)
  end
end

function L2(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    (ni, ne) = nVals(b, t)
    (im/2 - Base.MathConstants.eulergamma/π - 1/π*log(norm(boundaryPointDeriv(b, t)) / 2)) * ((k*ne)^2-(k*ni)^2) - 1/π*((k*ne)^2*log(k*ne) - (k*ni)^2*log(k*ni))
  else
    L(b, t, τ, k, bpc)-L1(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function L(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  im/2*((k*ne)^2*h0ne-(k*ni)^2*h0ni)
end

###################
# M functions
function M1(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    (ni, ne) = nVals(b, t)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    -1/(2π)*(j0ne-j0ni)
  end
end

function M2(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    (ni, ne) = nVals(b, t)
    1/π*log(ni/ne)
  else
    M(b, t, τ, k, bpc) - M1(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function M(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  im/2 * (h0ne-h0ni)
end

###################
# N functions
function N(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  xtP = boundaryPointDeriv(b, t)
  xτP = boundaryPointDeriv(b, τ)
  xpdot = dot(xtP, xτP)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  h1ni = bpc.h1NiMatrix[tIdx, τIdx]
  h1ne = bpc.h1NeMatrix[tIdx, τIdx]
  nbar = Nbar(b, t, τ, bpc)
  im/2*nbar*((k*ne)^2*h0ne-2*k*ne*h1ne/ad) +
    im/2*k*ne*xpdot/ad*h1ne - im/2*nbar*((k*ni)^2*h0ni-2*k*ni*h1ni/ad) -
    im/2*k*ni*xpdot/ad*h1ni
end

function Nbar(b::Boundary, t::Float64, τ::Float64, bpc::BoundaryPrecomputedCalcs)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  (xt, xτ) = xVals(b, t, τ)
  xtP = boundaryPointDeriv(b, t)
  xτP = boundaryPointDeriv(b, τ)

  xtm = (xt-xτ)
  dot(xtm, xtP)*dot(xtm, xτP)/ad^2
end

function N1(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  xtP = boundaryPointDeriv(b, t)

  if (t == τ)
    norm(xtP)^2*((k*ni)^2-(k*ne)^2)/(4π)
  else
    xτP = boundaryPointDeriv(b, τ)
    xpdot = dot(xtP, xτP)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    j1ni = bpc.j1NiMatrix[tIdx, τIdx]
    j1ne = bpc.j1NeMatrix[tIdx, τIdx]

    -1/(2π)*Nbar(b, t, τ, bpc)*((k*ne)^2*j0ne-2*k*ne*j1ne/ad) -
      k*ne*xpdot/(2π*ad)*j1ne +
      1/(2π)*Nbar(b, t, τ, bpc)*((k*ni)^2*j0ni-2*k*ni*j1ni/ad) +
      k*ni*xpdot/(2π*ad)*j1ni
  end
end

function N2(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    (xt, xτ) = xVals(b, t, τ)
    (ni, ne) = nVals(b, t)
    xtP = boundaryPointDeriv(b, t)
    normXtP = norm(xtP)
    normXtPSq = norm(xtP)^2

    normXtPSq/(4π) * (
      ((k*ne)^2-(k*ni)^2)*(π*im-1-2*Base.MathConstants.eulergamma) -
        2*(k*ne)^2*log(k*ne*normXtP/2)+2*(k*ni)^2*log(k*ni*normXtP/2)
    )
  else
    N(b, t, τ, k, bpc)-N1(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

###################
# A matrix blocks
#
A11(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ -(R(abs(j-l), N)*H1(b, j*π/N, l*π/N, k, bpc) + π/N*H2(b, j*π/N, l*π/N, k, bpc)) for j in 0:(2N-1), l in 0:(2N-1) ]

function A12(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  @DArray [ -(R(abs(j-l), N)*M1(b, j*π/N, l*π/N, k, bpc)*γ(b, N, l) +
    π/N*M2(b, j*π/N, l*π/N, k, bpc)*γ(b,N,l))
    for j in 0:(2N-1), l in 0:(2N-1) ]
end

function A21entry(b::Boundary, N::Int64, k::Complex{Float64}, j::Int64, l::Int64, bpc::BoundaryPrecomputedCalcs)
  xjP = boundaryPointDeriv(b, j*π/N)
  xlP = boundaryPointDeriv(b, l*π/N)
  R(abs(j-l), N) * (L1(b, j*π/N, l*π/N, k, bpc)*dot(xjP, xlP) - N1(b, j*π/N, l*π/N, k, bpc)) +
    π/N*(L2(b, j*π/N, l*π/N, k, bpc)*dot(xjP, xlP) - N2(b, j*π/N, l*π/N, k, bpc))
end
A21(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ A21entry(b, N, k, j, l, bpc) for j in 0:(2N-1), l in 0:(2N-1) ]

A22(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ R(abs(j-l), N)*H1star(b, j*π/N, l*π/N, k, bpc) + π/N*H2star(b, j*π/N, l*π/N, k, bpc) for j in 0:(2N-1), l in 0:(2N-1) ]

###################
# Other matrix blocks
γMatrix(b::Boundary, N::Int64) = Diagonal([γ(b, N, l) for l in 0:(2N-1)])

###################
# Full matrix block
function A(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  a11 = A11(b, N, k, bpc)
  a12 = A12(b, N, k, bpc)
  a21 = A21(b, N, k, bpc)
  a22 = A22(b, N, k, bpc)
  2*[ Matrix(1.0I, 2N, 2N) zeros(2N,2N); zeros(2N,2N) γMatrix(b, N)] - [convertAndClose(a11) convertAndClose(a12); convertAndClose(a21) convertAndClose(a22)]
end

function A(b::Boundary, N::Int64, k::Complex{Float64})
  bpc = doPrecalculation(b, N, k)
  A(b, N, k, bpc)
end

###################
# A'(k) derivative functions

################
# H functions
function H1′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    # Note: Need to check that we can use the k values at t, not τ, this is a deviation from the Heider paper
    (ni, ne) = nVals(b, t)
    (xt, xτ) = xVals(b, t, τ)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    k*ad/(2π)*dot(normVector(b, τ), (xt-xτ)/ad)*(ni^2*j0ni-ne^2*j0ne)
  end
end

function H2′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    H′(b, t, τ, k, bpc)-H1′(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function H′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (xt, xτ) = xVals(b, t, τ)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  im*k*ad/2*dot(normVector(b, τ), (xt-xτ)/ad)*(ne^2*h0ne-ni^2*h0ni)
end

H1star′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H1′(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))
H2star′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H2′(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))
Hstar′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = H′(b, τ, t, k, bpc)*norm(boundaryPointDeriv(b, τ))

###################
# L functions
function L1′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  if (t == τ)
    k*(ni^2-ne^2) / π
  else
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    j1ni = bpc.j1NiMatrix[tIdx, τIdx]
    j1ne = bpc.j1NeMatrix[tIdx, τIdx]
    1/(2π)*k*(k*ad*(ne^3*j1ne-ni^3*j1ni)-2*ne^2*j0ne+2*ni^2*j0ni)
  end
end

function L2′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    (ni, ne) = nVals(b, t)
    k/π*(-2*ne^2*log(k*ne)+2*ni^2*log(k*ni)-(ne-ni)*(ne+ni)*(2*log(norm(boundaryPointDeriv(b, t)) / 2)-π*im+2*Base.MathConstants.eulergamma+1))
  else
    L′(b, t, τ, k, bpc)-L1′(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function L′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  h1ni = bpc.h1NiMatrix[tIdx, τIdx]
  h1ne = bpc.h1NeMatrix[tIdx, τIdx]
  -k*im/2*(k*ad*(ne^3*h1ne-ni^3*h1ni)-2*ne^2*h0ne+2*ni^2*h0ni)
end

###################
# M functions
function M1′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    (ni, ne) = nVals(b, t)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j1ni = bpc.j1NiMatrix[tIdx, τIdx]
    j1ne = bpc.j1NeMatrix[tIdx, τIdx]
    ad/(2π)*(ne*j1ne-ni*j1ni)
  end
end

function M2′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    0.0+0im
  else
    M′(b, t, τ, k, bpc) - M1′(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

function M′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h1ni = bpc.h1NiMatrix[tIdx, τIdx]
  h1ne = bpc.h1NeMatrix[tIdx, τIdx]
  -im*ad/2 * (ne*h1ne-ni*h1ni)
end

###################
# N functions
function N′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  xtP = boundaryPointDeriv(b, t)
  xτP = boundaryPointDeriv(b, τ)
  xpdot = dot(xtP, xτP)
  (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
  ad = bpc.absDiff[tIdx, τIdx]
  h0ni = bpc.h0NiMatrix[tIdx, τIdx]
  h0ne = bpc.h0NeMatrix[tIdx, τIdx]
  h1ni = bpc.h1NiMatrix[tIdx, τIdx]
  h1ne = bpc.h1NeMatrix[tIdx, τIdx]
  k*im/2 * (
    k*ad*Nbar(b, t, τ, bpc) * (ni^3*h1ni-ne^3*h1ne) +
    xpdot*(ne^2*h0ne-ni^2*h0ni)
  )
end

function N1′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  (ni, ne) = nVals(b, t)
  xtP = boundaryPointDeriv(b, t)
  xτP = boundaryPointDeriv(b, τ)

  if (t == τ)
    norm(xtP)^2*k*(ni^2-ne^2)/(2π)
  else
    xτP = boundaryPointDeriv(b, τ)
    xpdot = dot(xtP, xτP)
    (tIdx, τIdx) = tAndτIndices(bpc, t, τ)
    ad = bpc.absDiff[tIdx, τIdx]
    j0ni = bpc.j0NiMatrix[tIdx, τIdx]
    j0ne = bpc.j0NeMatrix[tIdx, τIdx]
    j1ni = bpc.j1NiMatrix[tIdx, τIdx]
    j1ne = bpc.j1NeMatrix[tIdx, τIdx]
    k/(2π)*(
      k*ad*Nbar(b, t, τ, bpc)*(ne^3*j1ne-ni^3*j1ni) +
      xpdot*(ni^2*j0ni-ne^2*j0ne)
    )
  end
end

function N2′(b::Boundary, t::Float64, τ::Float64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  if (t == τ)
    (ni, ne) = nVals(b, t)
    xtP = boundaryPointDeriv(b, t)
    normXtP = norm(xtP)
    normXtPSq = norm(xtP)^2

    k*normXtPSq/(2π) * (
      -2*ne^2*log(k/2*ne*normXtP) +
      2*ni^2*log(k/2*ni*normXtP) -
      (2+2*Base.MathConstants.eulergamma-π*im)*(ne-ni)*(ne+ni)
    )
  else
    N′(b, t, τ, k, bpc)-N1′(b, t, τ, k, bpc)*log4Term(t, τ)
  end
end

###################
# A matrix blocks′
#
A11′(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ -(R(abs(j-l), N)*H1′(b, j*π/N, l*π/N, k, bpc) + π/N*H2′(b, j*π/N, l*π/N, k, bpc)) for j in 0:(2N-1), l in 0:(2N-1) ]

function A12′(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  @DArray [ -(R(abs(j-l), N)*M1′(b, j*π/N, l*π/N, k, bpc)*γ(b, N, l) +
    π/N*M2′(b, j*π/N, l*π/N, k, bpc)*γ(b,N,l))
    for j in 0:(2N-1), l in 0:(2N-1) ]
end

function A21entry′(b::Boundary, N::Int64, k::Complex{Float64}, j::Int64, l::Int64, bpc::BoundaryPrecomputedCalcs)
  xjP = boundaryPointDeriv(b, j*π/N)
  xlP = boundaryPointDeriv(b, l*π/N)
  R(abs(j-l), N) * (L1′(b, j*π/N, l*π/N, k, bpc)*dot(xjP, xlP) - N1′(b, j*π/N, l*π/N, k, bpc)) +
    π/N*(L2′(b, j*π/N, l*π/N, k, bpc)*dot(xjP, xlP) - N2′(b, j*π/N, l*π/N, k, bpc))
end
A21′(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ A21entry′(b, N, k, j, l, bpc) for j in 0:(2N-1), l in 0:(2N-1) ]

A22′(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs) = @DArray [ R(abs(j-l), N)*H1star′(b, j*π/N, l*π/N, k, bpc) + π/N*H2star′(b, j*π/N, l*π/N, k, bpc) for j in 0:(2N-1), l in 0:(2N-1) ]

###################
# Full matrix block
function A′(b::Boundary, N::Int64, k::Complex{Float64}, bpc::BoundaryPrecomputedCalcs)
  a11′ = A11′(b, N, k, bpc)
  a12′ = A12′(b, N, k, bpc)
  a21′ = A21′(b, N, k, bpc)
  a22′ = A22′(b, N, k, bpc)
  - [convertAndClose(a11′) convertAndClose(a12′); convertAndClose(a21′) convertAndClose(a22′)]
end

function A′(b::Boundary, N::Int64, k::Complex{Float64})
  bpc = doPrecalculation(b, N, k)
  A′(b, N, k, bpc)
end

###################
# Utility functions
function tAndτIndices(bpc::BoundaryPrecomputedCalcs, t::Float64, τ::Float64)
  (Base.round(Int, t*bpc.N/π)+1, Base.round(Int, τ*bpc.N/π)+1)
end

function convertAndClose(a::DistributedArrays.DArray)
  tmp = Base.convert(Array, a)
  close(a)
  tmp
end

function doPrecalculation(b::Boundary, N::Int64, k::Complex{Float64})
  range = 0:(2N-1)
  nis = [nVals(b, l*π/N)[1] for l in range]
  nes = [nVals(b, l*π/N)[2] for l in range]

  tτMatrix = [[l*π/N,j*π/N] for l in range, j in range]
  absDiffM = absDiff(b, tτMatrix)

  localNiArgs = k*nis.*absDiffM
  localNeArgs = k*nes.*absDiffM
  niArgs = distribute(localNiArgs)
  neArgs = distribute(localNeArgs)
  j0NiM = besselj.(0, niArgs)
  j0NeM = besselj.(0, neArgs)
  j1NiM = besselj.(1, niArgs)
  j1NeM = besselj.(1, neArgs)

  # Set diagonal to NaN to prevent exception for hankel functions
  [localNiArgs[i, i] = NaN for i in 1:2N]
  [localNeArgs[i, i] = NaN for i in 1:2N]
  close(niArgs)
  close(neArgs)
  niArgs = distribute(localNiArgs)
  neArgs = distribute(localNeArgs)
  h0NiM = h0.(niArgs)
  h0NeM = h0.(neArgs)
  h1NiM = h1.(niArgs)
  h1NeM = h1.(neArgs)
  close(niArgs)
  close(neArgs)
  bpc = BoundaryPrecomputedCalcs(b, N, k, convertAndClose(j0NiM), convertAndClose(j0NeM), convertAndClose(j1NiM), convertAndClose(j1NeM),
      convertAndClose(h0NiM), convertAndClose(h0NeM), convertAndClose(h1NiM), convertAndClose(h1NeM), absDiffM)
  # @everywhere gc()
  bpc
end

ext: jl


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(z; timeaxis=0.0:0.001:3.0)
png(path*"IS_ex0.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
plot(𝑆; timeaxis=0.0:0.001:3.0)
png(path*"IS_ex1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
plot(𝑆,timeaxis=0.0:0.001:3.0,view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_ex2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->10*t^2+50,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(2t),t->25 + 10.5*sin(t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
plot(𝑆)

using ISA, Plots
𝐶₀ = AMFMtriplet(t->0.2+0.8cos(11t), t->200.0, 0.0)
𝐶₁ = AMFMtriplet(t->exp(-t^2),t->100.0,0.0)
𝑆 = compSet([𝐶₀,𝐶₁])
plot(𝑆,realProj=true)
png(path*"realProject.png")

using ISA, Plots
ψ₀ = AMFMcomp(t->t,t->25cos(t)+50,0.0)
plot(ψ₀, timeaxis=0.0:0.001:3.0)
png(path*"defaultColor.png")

plot(ψ₀, timeaxis=0.0:0.001:3.0, colorMap="cubeYF")
png(path*"cubeYFcolor.png")

myColorMap  = reverse(colormap("Blues",256))
plot(ψ₀, timeaxis=0.0:0.001:3.0, colorMap = myColorMap)
png(path*"myColorMap.png")

#---------------------------------------------------

𝐶₀ = AMFMtriplet(t->0.2+0.8cos(11t), t->200.0, 0.0)
𝐶₁ = AMFMtriplet(t->exp(-abs(t/3)), t->100t, 0.1)
𝐶₂ = AMFMtriplet(t->𝒩ᵤ(t; μ=1.5, σ=1.0), t->150 + 125sin(5t), π)
𝐶₃ = AMFMtriplet(t->u(t-1.725)-u(t-2.475), t->50, 0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂,𝐶₃])
plot(𝑆; timeaxis=0.0:0.001:3.0)
png(path*"defaultIS.png")

plot(𝑆,view="TF", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"TFviewIS.png")

plot(𝑆,view="TR", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"TRviewIS.png")

plot(𝑆,view="FR", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"FRviewIS.png")


text:
#--
# Threads.nthreads() displays number of available threads
# if ==1, start julia with JULIA_NUM_THREADS=4 julia
#--

#--
# If no graphics (X-windows) required, then ENV["GKSwstype"]="nul"
#--

using FDModBiotCyl
#
using CPUTime
using ProgressMeter
using Plots
using Printf

#--Suppressing graphics
#ENV["GKSwstype"]="nul"

@__DIR__
@__FILE__

#================
Drawing a model
================#
function drawmodel(Field_Var,nz,dz,nr,dr,LPML_r,LPML_z)
   println("Drawing a model...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
#   plt1=heatmap(rvec,zvec,Vp,yflip=true,xlabel="R (m)",ylabel="Z (m)",ratio=1)
   plt1=heatmap(rvec,zvec,Field_Var,yflip=true,xlabel="R (m)",ylabel="Z (m)")
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)
   display(plot(plt1))
   tmp_filename=@sprintf "tmp_FDmodel.png"
   tmp_filename_full=string(@__DIR__,"/",tmp_filename)
   png(tmp_filename_full)
   println("--------")
   println("Model config figure is saved in:", tmp_filename_full)
   println("--------")
end

#================
Drawing a snapshot
================#
function drawsnap(Data1,nz,dz,nr,dr,LPML_r,LPML_z,title_str)
   println("Drawing a snapshot...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
   plt1=heatmap(rvec,zvec,Data1,yflip=true,xlabel="R (m)",ylabel="Z (m)",legend=:none)
   heatmap!(title=title_str)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)
   display(plot(plt1))
end
# Two snapshots(1,2)
function drawsnap_1by2(Data1,Data2,nz,dz,nr,dr,LPML_r,LPML_z,title_str1,title_str2)
   println("Drawing a snapshot...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
   plt1=heatmap(rvec,zvec,Data1,yflip=true,xlabel="R (m)",ylabel="Z (m)",legend=:none)
   heatmap!(title=title_str1)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   plt2=heatmap(rvec,zvec,Data2,yflip=true,xlabel="R (m)",ylabel="Z (m)",legend=:none)
   heatmap!(title=title_str2)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)

   display(plot(plt1,plt2,layout=(1,2)))
end

#=========================
Creating a homogeneous model
=========================#
function makemodel_homogeneous_PoroElastic()
   #Necessary input matrices to the main loop function are:
   #Rho: Bulk density
   #Rhof: Fluid density
   #M,C,H: Poroelastic moduli
   #G: Formation shear modulus
   #D1,D2: Poroelastic moduli in Ou's formulation
   #Flag_AC,Flag_E: flag specifying acoustic region or elastic region
   #Other input parameters to the main loop function are:
   #
   # This function creates the above input parameters assuming a
   # homogeneous poroelastic media

   # Model size
   nr=251 #samples
   nz=351 #samples
   dr=0.5 #meter
   dz=0.5 #meter

   #==============================================
   Initializing grain/frame parameters
   ==============================================#
   #First creating grain/frame parameters and then converting them
   #to necessary poroelastic parameters for our FD
   #----solid phase-----
   Km=zeros(nz,nr) #frame
   Ks=zeros(nz,nr) #grain
   G=zeros(nz,nr) #bulk
   Rho=zeros(nz,nr) #bulk (weighted average of Rhos and Rhof)
   Rhos=zeros(nz,nr) #grain
   #---fluid phase----
   Kf=zeros(nz,nr)
   Rhof=zeros(nz,nr)
   Kappa0=zeros(nz,nr) #static permeability (m^2, 1D=9.869E-13m^2)
   Eta=zeros(nz,nr) #dynamic viscosity Pa.s=kg/(m.s)
   #---other parameters---
   Phi=zeros(nz,nr)
   Tot=zeros(nz,nr) #Tortuosity factor (see Sidler 2014)
   #End initializing poroelastic parameters

   #==============================================
   Fluid parameters
   ==============================================#
   #Assuming that fluid properties are identical everywhere
   Kf[:,:]=ones(nz,nr)*2.25*10^9
   Rhof[:,:]=ones(nz,nr)*1000
   Eta[:,:]=ones(nz,nr)*1.0*10^(-3) #water 1E-3 Pa.s

   Flag_AC=zeros(nz,nr)

   #==========================
   Homogeneous PoroElastic model
   ==========================#
   Kappa0[:,:]=ones(nz,nr)*1*9.869*10^(-13) #m^2
   Phi[:,:]=ones(nz,nr)*0.3
   #Tortuosity factor from Ou
   m_Ou=8.0 #
   alpha_inf_Ou=3.0 #tortuosity
   Tot=ones(nz,nr)*(1+2/m_Ou)*alpha_inf_Ou
   #Solid parameters from dry_Vp, dry_Vs, grain density, and grain bulk modulus
   Vp_dry=2170.0
   Vs_dry=1360.0
   Rhos[:,:]=ones(nz,nr)*3143.0
   Ks[:,:]=ones(nz,nr)*5*1E10
   #
   Rho_dry=(ones(nz,nr)-Phi).*Rhos
   G[:,:]=Vs_dry^2*Rho_dry[:,:]
   Km[:,:]=Vp_dry^2*Rho_dry[:,:]-4/3*G[:,:]
   Rho=(ones(nz,nr)-Phi).*Rhos+Phi.*Rhof

   Flag_E=zeros(nz,nr)

   #=====================================================
   Converting grain/frame parameters into poroelastic parameters
   =====================================================#
   #Rho, Rhof, M, C, H, G, D1, D2
   #Rho as is
   #Rhof as is
   Alpha=ones(nz,nr)-Km./Ks
   M=( (Alpha-Phi)./Ks+Phi./Kf ).^(-1)
   C=M.*Alpha
   H=Km+4/3*G+M.*Alpha.^2
   #G as is
   D1=Tot.*Rhof./Phi #(1+2/m)*alpha_inf*rhof/phi
   D2=Eta./Kappa0

   #Correction of H and C for elastic media
   for ir=1:nr
      for iz=1:nz
         if(Flag_E[iz,ir]==1)
            H[iz,ir]=Ks[iz,ir]+4/3*G[iz,ir]
            C[iz,ir]=Ks[iz,ir] #Ou's
         end

         if(Flag_AC[iz,ir]==1)
#            D1[iz,ir]=1.0 #Guan's
            D1[iz,ir]=1000. #Ou's D1=rhof
            D2[iz,ir]=0. #Ou's D2=0
         end
      end
   end

   return nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E
end

#=========================
Computing Vp_sat, Vs_sat
from poroelastic properties
=========================#
function Gassmann(vp_dry,vs_dry,rhos,phi,Ks)
#   vp_dry=5170; #dry vp
#   vs_dry=3198; #dry vs
#   rhos=3143; #grain density
#   phi=0.3; #porosity
#   Ks=10E10; #grain bulk modulus
   Kf=1000*1500^2; #fluid bulk modulus
   Rhof=1000;

   rho_dry=(1-phi)*rhos;
   rho_bulk=(1-phi)*rhos+phi*Rhof;
   G=vs_dry^2*rho_dry;
   K_dry=vp_dry^2*rho_dry-4/3*G;
   alpha=1-K_dry./Ks;
   M=( (alpha-phi)./Ks+phi./Kf ).^(-1);
   C=M.*alpha;
   H=K_dry+4/3*G+M.*alpha.^2; #H=K_undrained+4/3*G (see 2.6c, Guan2011)

   #Gassmann K_undrained=alpha^2*M+K_dry (see Guan2011)
   #Gassmann K_sat=K_undrained
   Ksat=K_dry+(1-K_dry/Ks)^2/(
       phi/Kf+(1-phi)/Ks-K_dry/Ks^2);
   K_undrained=K_dry+alpha^2*M;
   Vp0=sqrt(H/rho_bulk); #->sqrt((K_undrained+4/3G)/rho_bulk)
   Vs0=sqrt(G/rho_bulk); #->sqrt((K_undrained+4/3G)/rho_bulk)
   #Skempton Coefficient
   B=(1/K_dry-1/Ks)/(1/K_dry-1/Ks+phi*(1/Kf-1/Ks));

   return Vp0,Vs0,rho_bulk
end

#----MAIN loop function with src----
function main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   src_index,src_dn,src_func,
   Flag_AC,Flag_E,
   nrec,rec_vr,rec_vz,rec_pf,index_allrec_vr,index_allrec_vz,index_allrec_pf,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

#println("start")

#--
Flag_vf_zero=get_Flag_vf_zero(Flag_AC,Flag_E,nr,nz)

#initializing matrices of PML
Prz_T,Pzz_T,Rz_T,Srz_T,PzzPE_T,RzPE_T,
Prz_B,Pzz_B,Rz_B,Srz_B,PzzPE_B,RzPE_B,
Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,PrrPE_R,RrPE_R,RpPE_R,
Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,PrrPE_TR,PzzPE_TR,RrPE_TR,RpPE_TR,RzPE_TR,
Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,PrrPE_BR,PzzPE_BR,RrPE_BR,RpPE_BR,RzPE_BR=
init_PML_Por(LPML_r,LPML_z,nr,nz)

memT_vr,memT_vz,memT_trr,memT_tpp,memT_tzz,memT_trz,memT_vfr,memT_vfz,memT_pf,
memB_vr,memB_vz,memB_trr,memB_tpp,memB_tzz,memB_trz,memB_vfr,memB_vfz,memB_pf,
memR_vr,memR_vz,memR_trr,memR_tpp,memR_tzz,memR_trz,memR_vfr,memR_vfz,memR_pf=
init_memory_variables_Por(LPML_r,LPML_z,nr,nz)
#return

@showprogress for ii=1:nt

#----Time at (ii-1)*dt---(updating velocities)-----

#--PML: save velocity at previous step
PML_save_vel!(memT_vr,memT_vz,memT_vfr,memT_vfz,
              memB_vr,memB_vz,memB_vfr,memB_vfz,
              memR_vr,memR_vz,memR_vfr,memR_vfz,
              vr,vz,vfr,vfz,nr,nz,LPML_z,LPML_r)

# Main velocity update!
update_velocity_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,LPML_z,LPML_r)

# PML: updating velocity
PML_update_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
              Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,
              Prz_T,Pzz_T,PzzPE_T,
              Prz_B,Pzz_B,PzzPE_B,
              Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,
              PrrPE_R,RrPE_R,RpPE_R,
              Prz_TR,Pzz_TR,PzzPE_TR,
              Prr_TR,Qrp_TR,Pzr_TR,Qzp_TR,Rr_TR,Rp_TR,Rrz_TR,
              PrrPE_TR,RrPE_TR,RpPE_TR,
              Prz_BR,Pzz_BR,PzzPE_BR,
              Prr_BR,Qrp_BR,Pzr_BR,Qzp_BR,Rr_BR,Rp_BR,Rrz_BR,
              PrrPE_BR,RrPE_BR,RpPE_BR,
              LPML_z,LPML_r)

# Making sure RightBC for velocity (zero values -> vr)
ApplyBCRight_vel!(vr,vz,
      trr,tpp,tzz,trz,
      vfr,vfz,pf,
      Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt)

#---velocity B.Cs-----
#Periodic Left Edge
ApplyBCLeft_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                          Rho,Rhof,D1,D2,Flag_vf_zero,
                          nr,nz,dr,dz,dt,
                          Prz_T,Pzz_T,PzzPE_T,
                          Prz_B,Pzz_B,PzzPE_B,
                          LPML_z)

#PML: updating memory variables for stress (Rx and Sxx) using velocity at two-time steps
PML_update_memRS!(Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  memT_vr,memT_vz,memT_vfr,memT_vfz,
                  memB_vr,memB_vz,memB_vfr,memB_vfz,
                  memR_vr,memR_vz,memR_vfr,memR_vfz,
                  vr,vz,vfr,vfz,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#--src injection (velocity)
srcamp=src_func[ii]
srcamp=srcamp*dt/Rhof[1,1]
srcapply!(vz,src_index,src_dn,srcamp)
#srcapply!(vr,src_index,src_dn,srcamp)
#srcapply!(vfz,src_index,src_dn,srcamp)
#srcapply!(vfr,src_index,src_dn,srcamp)

#----Time at (ii-1)*dt+dt/2----(updating stress)---

#--PML: saving stress at the previous step
PML_save_stress!(memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                 memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                 memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                 trr,tpp,tzz,trz,pf,nr,nz,LPML_z,LPML_r)

# Main stress update!
update_stress_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
     M,C,H,G,nr,nz,dr,dz,dt,LPML_z,LPML_r)

# PML: updating stress
PML_update_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                  M,C,H,G,nr,nz,dr,dz,dt,
                  Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  LPML_z,LPML_r)

# Making sure RightBC for stress (zero values -> trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
     trz,
     G,nr,nz,dr,dz,dt)

#---stress B.Cs
ApplyBCLeft_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                    M,C,H,G,nr,nz,dr,dz,dt,
                    Rz_T,Srz_T,RzPE_T,
                    Rz_B,Srz_B,RzPE_B,
                    LPML_z)

#--src injection (stress:monopole src)
#srcamp=src_func[ii]
#srcamp=srcamp*dt/Rhof[1,1]
#srcapply!(pf,src_index,src_dn,srcamp)

#--PML: updating memory variables for velocity (Pxx and Qxx) using stress at two-time steps
PML_update_memPQ!(Prz_T,Pzz_T,PzzPE_T,
                  Prz_B,Pzz_B,PzzPE_B,
                  Prr_R,Qrp_R,Pzr_R,Qzp_R,PrrPE_R,
                  Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,PrrPE_TR,PzzPE_TR,
                  Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,PrrPE_BR,PzzPE_BR,
                  memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                  memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                  memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                  trr,tpp,tzz,trz,pf,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#-----End of updating field----

#Receiver field
getRecData_from_index!(vr,rec_vr,index_allrec_vr,nrec,ii)
getRecData_from_index!(vz,rec_vz,index_allrec_vz,nrec,ii)
getRecData_from_index!(pf,rec_pf,index_allrec_pf,nrec,ii)

#--Snapshots
check_snap=findall(x -> x==ii,itvec_snap)
if (length(check_snap)!=0)
   println("ii=",ii)
   cnt_snap=Int(check_snap[1])
   get_snapshots!(snapshots_vz,cnt_snap,vz,nr,nz)
   get_snapshots!(snapshots_vr,cnt_snap,vr,nr,nz)
   get_snapshots!(snapshots_pf,cnt_snap,pf,nr,nz)

#   drawsnap(snapshots_vz[:,:,cnt_snap],nz,dz,nr,dr,
#            LPML_r,LPML_z,"Snapshots: Vz")

   drawsnap_1by2(snapshots_vz[:,:,cnt_snap],snapshots_pf[:,:,cnt_snap],nz,dz,nr,dr,
            LPML_r,LPML_z,"Snapshots: Vz","Snapshots: Pf")

end

end #i Time

end #function
#------------------------------

#-----ENTRY POINT HERE----
CPUtic()
start=time()

#======================
General variables
======================#

#time samples
dt=1E-4
nt=751

T=(nt-1)*dt
tvec=range(0.0,T,length=nt) # range object (no memory allocation)
tvec=collect(tvec) # a vector

#PML thickness in samples
LPML_r=30
LPML_z=30

#======================
Creating a model
======================#
nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E=makemodel_homogeneous_PoroElastic()
drawmodel(H,nz,dz,nr,dr,LPML_r,LPML_z)

#======================
Stability check
======================#
Vmax=maximum((H./Rho).^(0.5))
check_stability01(dt,dr,Vmax,0)
check_stability01(dt,dz,Vmax,0)

#======================
Creating a src wavelet
======================#
f0=100 #src Freq
delay=1/f0*1.0
#src_func=myricker2(tvec,f0,delay,2) #when using 2nd derivative Gaussian (Ricker with a negative sign)
src_func=myricker2(tvec,f0,delay,1) #when using 1st derivative Gaussian
#src_func=myricker2(tvec,f0,delay,3) #when using 3rd derivative Gaussian
tmp_maxamp=maximum(map(abs,src_func))
src_func=src_func/tmp_maxamp*1E+9 #this scaling is just for visualization

display(plot(tvec,src_func[:],title="src function"))
tmpfilename_src=string(@__DIR__,"/src_function.png")
png(tmpfilename_src)
println("Source signature figure saved in ",tmpfilename_src)

#================================
Point src geometry
================================#
srcgeom=zeros(1,2) #(z,r)
srcgeom[1,1]=(nz-1)*dz/2 #z meter
#srcgeom[1,2]=(nr-LPML_r+1)*dr #r meter
srcgeom[1,2]=0*dr #r meter
#Gaussian point src
wsize=7 #window size (odd number)
wsigma=1 #std
src_index,src_dn=get_srcindex_pGauss(srcgeom,dr,dz,wsize,wsigma)
# if you want to see the amplitude distribution,
# plot(src_index[1,:],src_index[2,:],src_dn[:],marker=:circle,seriestype=:scatter,zcolor=src_dn[:],camera=(90,90))

#==============================
Initializing field variables
==============================#
vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf=init_fields_Por(nz,nr) #

# Creating PML profile
LPML_r,LPML_z,PML_Wr,PML_Wz,
PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2=init_PML_profile(LPML_r,LPML_z,Vmax,dr,dz,nr,f0)
#PML_check(LPML_r,LPML_z,Vmax,dr,dz,f0)

#===================
Receiver geometry
====================#
#--Fluid pressure at borehole center
nrec=nz
recgeom=zeros(nrec,2) #(z,r)
for irec=1:nrec
   recgeom[irec,1]=dz*(irec-1)
   recgeom[irec,2]=25.0
end
rec_vr,rec_vz,index_allrec_vr,index_allrec_vz=init_receiver_geophone(recgeom,nrec,dr,dz,nr,nz,nt)
rec_pf,index_allrec_pf=init_receiver_hydrophone(recgeom,nrec,dr,dz,nr,nz,nt)

#==============================
Snapshot settings
==============================#
nskip=30 #change here to adjust time sampling (1~nt)
snapshots_vz,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_vr,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_pf,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
#error()

#==============================
Start main FD Loop
==============================#

main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   src_index,src_dn,src_func,
   Flag_AC,Flag_E,
   nrec,rec_vr,rec_vz,rec_pf,index_allrec_vr,index_allrec_vz,index_allrec_pf,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

CPUtoc()
println("elapsed real time: ", round(time() - start;digits=3)," seconds")

println("FD finished without a problem.")

println("Plotting receiver responses...")
offset=recgeom[:,1]-srcgeom[1]*ones(size(recgeom[:,1]))
plt1=heatmap(offset[:],tvec[:],rec_vz,xlabel="Offset (m)",ylabel="time (s)",title="Vz")
heatmap!(yflip=true)
plt2=heatmap(offset[:],tvec[:],rec_pf,xlabel="Offset (m)",ylabel="time (s)",title="Pf")
heatmap!(yflip=true)
display(plot(plt1,plt2,layout=(1,2)))


text:
# This file contains the low-level process wrapper for the solidity compiler

"""
    function solcWrapper()
Wrapper function for calling to solc.

# Arguments
- solc_binary: Location of the `solc` binary. If not given, the current default binary is used.
- stdin: Input to pass to `solc` via stdin.
- source_files: Path, or list of paths, of sources to compole.

"""
function solcWrapper(solc_binary::String,
                    stdin::String,
                    source_files::Union{Vector, String})

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"""Node is a worker node in Kubernetes. Each node will have a unique identifier in the cache (i.e. in etcd).

    IoK8sApiCoreV1Node(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - spec::IoK8sApiCoreV1NodeSpec : Spec defines the behavior of a node. https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#spec-and-status
    - status::IoK8sApiCoreV1NodeStatus : Most recently observed status of the node. Populated by the system. Read-only. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#spec-and-status
"""
mutable struct IoK8sApiCoreV1Node <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeStatus } # spec name: status

    function IoK8sApiCoreV1Node(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Node, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1Node, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1Node, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1Node, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiCoreV1Node, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiCoreV1Node

const _property_map_IoK8sApiCoreV1Node = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiCoreV1Node = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiCoreV1NodeSpec", Symbol("status")=>"IoK8sApiCoreV1NodeStatus")
Base.propertynames(::Type{ IoK8sApiCoreV1Node }) = collect(keys(_property_map_IoK8sApiCoreV1Node))
Swagger.property_type(::Type{ IoK8sApiCoreV1Node }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Node[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Node }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Node[property_name]

function check_required(o::IoK8sApiCoreV1Node)
    true
end

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


text:
using Plots
using ContinuousPOMDPTreeSearchExperiments
using POMDPs
using DiscreteValueIteration
using POMCPOW
using ParticleFilters
using POMDPToolbox

pyplot()

@show max_time = 1.0
@show max_depth = 20

pomdp = SimpleLightDark()

solvers = Dict{String, Union{Solver,Policy}}(
    "pomcpow" => begin
        rng = MersenneTwister(13)
        ro = ValueIterationSolver()
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(90.0),
                               # final_criterion=MaxTries(),
                               max_depth=max_depth,
                               max_time=max_time,
                               enable_action_pw=false,
                               k_observation=5.0,
                               alpha_observation=1/15.0,
                               estimate_value=FOValue(ro),
                               check_repeat_obs=false,
                               # default_action=ReportWhenUsed(-1),
                               rng=rng
                              )
    end,

    "pomcpdpw" => begin
        rng = MersenneTwister(15)
        ro = ValueIterationSolver()
        solver = PDPWSolver(tree_queries=10_000_000,
                            c=100.0,
                            max_depth=max_depth,
                            max_time=max_time,
                            enable_action_pw=false,
                            k_observation=4.0,
                            alpha_observation=1/10,
                            estimate_value=FOValue(ro),
                            check_repeat_obs=false,
                            # default_action=ReportWhenUsed(-1),
                            rng=rng
                           )
    end
   )

names = Dict("pomcpow"=>"POMCPOW", "pomcpdpw"=>"POMCP-DPW (≈QMDP)")
if !isdefined(:results) || isempty(results)
    results = Dict{String, Any}()
end

for (k, solver) in solvers
    if !haskey(results, k)
        planner = solve(solver, pomdp)
        rng = MersenneTwister(7)
        sim = HistoryRecorder(show_progress=true, max_steps=100, rng=rng, initial_state=3)
        filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(1_000),
                                           0.05, MersenneTwister(90_000)
                                          )
        hist = simulate(sim, pomdp, planner, filter)
        results[k] = hist
    end
end

smin = -10
smax = 20

tmax = maximum(length(state_hist(h)) for h in values(results))
plots = []
p = nothing
for (k, hist) in results
    vsh = collect(filter(s->!isterminal(pomdp,s), state_hist(hist)))
    bh = belief_hist(hist)[1:end-1]

    pts = Int[]
    pss = Int[]
    pws = Float64[]

    for t in 0:length(bh)-1
        b = bh[t+1]
        for s in smin:smax
            w = 10.0*sqrt(pdf(b, s))
            if 0.0<w<1.0
                w = 1.0
            end
            push!(pts, t)
            push!(pss, s)
            push!(pws, w)
        end
    end

    T = linspace(0.0, tmax)
    S = linspace(-1.0, 21.0)
    inv_grays = cgrad([RGB(1.0, 1.0, 1.0),RGB(0.0,0.0,0.0)])
    p = contour(T, S, (t,s)->abs(s-pomdp.light_loc),
            bg_inside=:black,
            fill=true,
            xlim=(0, tmax),
            ylim=(smin, smax),
            color=inv_grays,
            ylabel="State",
            cbar=false,
            legend=:topright
           )
    plot!(p, [0, tmax], [0, 0], linewidth=2, color="green", label="Goal", line=:dash)
    scatter!(p, pts, pss, color="lightblue", label="Belief Particles", markersize=pws, marker=stroke(0.1, 0.3))
    plot!(p, 0:length(vsh)-1, vsh, linewidth=3, color="orangered", label="Trajectory", title=names[k])
    ffamily = "Times"
    smfont = Plots.font(ffamily, 12)
    mfont = Plots.font(ffamily, 12)
    tfont = Plots.font(ffamily, 16)
    plot!(p, titlefont=tfont, legendfont=mfont, tickfont=smfont, guidefont=mfont)
    push!(plots, p)
end

xlabel!(p, "Time")
plot!(p, legend=false)

plot(plots...; layout=(length(plots),1))
gui()


text:
# https://www.hackerrank.com/challenges/30-hello-world/problem

input = readline()

function greet(input_string::String)
    println("Hello, World.")
    println(input_string)
end

greet(input)

max_stars_repo_path: Practice/Tutorials/30 Days of Code/Day 00 Hello, World.jl

max_stars_repo_name: mtaufanr/hackerrank-code-julia

max_issues_repo_name: mtaufanr/hackerrank-code-julia

max_forks_repo_name: mtaufanr/hackerrank-code-julia


text:
"""
    null1(N::BinaryNetwork)

Given a matrix `A`, `null1(A)` returns a matrix with the same dimensions,
where every interaction happens with a probability equal to the connectance of
`A`.
"""
function null1(N::BinaryNetwork)
    return linearfilter(N; α=[0.0, 0.0, 0.0, 1.0])
end

"""
    null3out(N::BinaryNetwork)

Given a matrix `A`, `null3out(A)` returns a matrix with the same dimensions,
where every interaction happens with a probability equal to the out-degree
(number of successors) of each species, divided by the total number of possible
successors.
"""
function null3out(N::BinaryNetwork)
  return linearfilter(N; α=[0.0, 1.0, 0.0, 0.0])
end

"""
    null3in(N::BinaryNetwork)

Given a matrix `A`, `null3in(A)` returns a matrix with the same dimensions,
where every interaction happens with a probability equal to the in-degree
(number of predecessors) of each species, divided by the total number of
possible predecessors.
"""
function null3in(N::BinaryNetwork)
  return linearfilter(N; α=[0.0, 0.0, 1.0, 0.0])
end

"""
    null2(N::BinaryNetwork)

Given a matrix `A`, `null2(A)` returns a matrix with the same dimensions, where
every interaction happens with a probability equal to the degree of each
species.
"""
function null2(N::BinaryNetwork)
  return linearfilter(N; α=[0.0, 0.5, 0.5, 0.0])
end


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c8e11800-1874-11eb-2e58-53acb96689fb
# Add several packages to our environment

begin
	using Compose
	using PlutoUI
end

# ╔═╡ a18e4c30-1880-11eb-0d80-e34d6fb8f5fe
md"# Math for FUN!!!"

# ╔═╡ a9d1b190-187e-11eb-02c5-2909b144bc33
md"# FRACTALS!!!  😊"

# ╔═╡ 73d06d30-187d-11eb-1b07-c16f07d51e2f
md"# What is a Fractal? 🤔
### A fractal is a shape where you can zoom in infinitely on a shape and it will still look similar to itself. "

# ╔═╡ a95adf70-1874-11eb-20f6-11f9b697522b
md"# _The Sierpinksi's triangle_"

# ╔═╡ 2884fbe0-187b-11eb-0be4-4d03f526a50c
md"Complexity = $(@bind number Slider(0:7, show_value=true))"

# ╔═╡ 98768e90-187c-11eb-19e5-17c8b90fe324
md"## Sierpinski's triangle is defined _recursively_:
- Sierpinski's triangle of complexity N is a figure in the form of a triangle which is made of 3 triangular figures which are themselves Sierpinski's triangles of complexity N-1. In other words, the Sierspinki Triangle is made up of three copies of itself which are around an upside down equilateral triangle with the same size.

- A Sierpinski's triangle of complexity 0 is a simple solid equilateral triangle."

# ╔═╡ 7838c9a0-21ef-11eb-1672-b3aca3f79b32
md"# How can I draw it?"

# ╔═╡ fe972230-21ef-11eb-2c2f-35e42e64de5f
step = 1

# ╔═╡ 468628d2-1880-11eb-27d7-71598afaf445
md"# See you in the next video."

# ╔═╡ 52c3bf82-21f0-11eb-1895-fd0e7807793a
md"# Bye!!! 👋"

# ╔═╡ 595093e0-21f1-11eb-2581-0db61d8c2262
md"# How did I program the Sierspinki Triangle with Julia?"

# ╔═╡ b3f0d570-1874-11eb-002c-fd0abf2c1d9d
md"To draw Sierpinski's triangle, we are going to use two external packages, [_Compose.jl_](https://giovineitalia.github.io/Compose.jl/latest/tutorial) and [_PlutoUI.jl_](https://github.com/fonsp/PlutoUI.jl). Let's import the packages into our notebook. (You can download them in the terminal or Atom REPL.)
A package contains a coherent set of functionality that you cannot do with Julia."

# ╔═╡ 673a7362-187b-11eb-00ec-4983d864eb34
triangle() = compose(context(), polygon([(1, 1), (0, 1), (1 / 2, 0)]))

# ╔═╡ aa6b60e0-187b-11eb-1d9f-3dfbe0bdf48c
# It does not matter which order you define the building blocks (functions) of the
# program in. The best way to organize code is the one that promotes understanding.

function place_in_3_corners(t)
	# Uses the Compose library to place 3 copies of t
	# in the 3 corners of a triangle.
	# treat this function as a black box,
	# or learn how it works from the Compose documentation here https://giovineitalia.github.io/Compose.jl/latest/tutorial/#Compose-is-declarative-1
	compose(context(),
			(context(1 / 4,   0, 1 / 2, 1 / 2), t),
			(context(0, 1 / 2, 1 / 2, 1 / 2), t),
			(context(1 / 2, 1 / 2, 1 / 2, 1 / 2), t))
end

# ╔═╡ c4eaac50-187b-11eb-3c92-29f462b7e2d4
function sierpinski(n)
	if n == 0
		triangle()
	else
		t = sierpinski(n - 1) # recursively construct a smaller sierpinski's triangle
		place_in_3_corners(t) # place it in the 3 corners of a triangle
	end
end

# ╔═╡ d0bf15c0-187b-11eb-3a7a-c50bb793cdd8
sierpinski(number)

# ╔═╡ b9bc0ab0-1879-11eb-3b22-c558830517c3
sierpinski.(0:6)

# ╔═╡ 15a8be20-21f0-11eb-2938-0f4033a36423
sierpinski(step - 1)

# ╔═╡ Cell order:
# ╟─a18e4c30-1880-11eb-0d80-e34d6fb8f5fe
# ╟─a9d1b190-187e-11eb-02c5-2909b144bc33
# ╟─73d06d30-187d-11eb-1b07-c16f07d51e2f
# ╟─a95adf70-1874-11eb-20f6-11f9b697522b
# ╟─2884fbe0-187b-11eb-0be4-4d03f526a50c
# ╟─d0bf15c0-187b-11eb-3a7a-c50bb793cdd8
# ╠═b9bc0ab0-1879-11eb-3b22-c558830517c3
# ╟─98768e90-187c-11eb-19e5-17c8b90fe324
# ╟─7838c9a0-21ef-11eb-1672-b3aca3f79b32
# ╠═fe972230-21ef-11eb-2c2f-35e42e64de5f
# ╟─15a8be20-21f0-11eb-2938-0f4033a36423
# ╟─468628d2-1880-11eb-27d7-71598afaf445
# ╟─52c3bf82-21f0-11eb-1895-fd0e7807793a
# ╟─595093e0-21f1-11eb-2581-0db61d8c2262
# ╟─b3f0d570-1874-11eb-002c-fd0abf2c1d9d
# ╠═c8e11800-1874-11eb-2e58-53acb96689fb
# ╠═673a7362-187b-11eb-00ec-4983d864eb34
# ╠═aa6b60e0-187b-11eb-1d9f-3dfbe0bdf48c
# ╠═c4eaac50-187b-11eb-3c92-29f462b7e2d4

max_stars_repo_path: Fractals/Sierspinki's Triangle.jl

max_stars_repo_name: tinfungster/FRACTALS

max_issues_repo_path: Fractals/Sierspinki's Triangle.jl

max_issues_repo_name: tinfungster/FRACTALS

max_forks_repo_path: Fractals/Sierspinki's Triangle.jl

max_forks_repo_name: tinfungster/FRACTALS


text:
import FinanceLib
import Dates

@testset "FinanceLib                                                  " begin

  @testset "tv" begin
    @test FinanceLib.yearFrac(Dates.Date(2027,2,12), Dates.Date(2018,2,12)) ≈ -8.999315537303216
    @test FinanceLib.invYearFrac(Dates.Date(2027,2,12), -8.999315537303216) == Dates.Date(2018,2,12)

    @test FinanceLib.disFactAnnual(0.07) == 0.9345794392523364
    @test FinanceLib.disFact(0.09, 3) == 0.7721834800610642
    @test FinanceLib.fwdDisFact((0.07, 1), (0.09, 3)) == 0.8262363236653387

    @test FinanceLib.xdisFact(0.09, Dates.Date(2015,3,15), Dates.Date(2018,10,8)) == 0.7353328680759499

    @test FinanceLib.tMul(0.06/12, -120.0) == 0.5496327333641637
    @test FinanceLib.tMul(0.06, -10.0, 12.0) == 0.5496327333641637
    @test FinanceLib.rateGwth(7.35, 8.52, 5.0) == -0.029111071029244595
    @test FinanceLib.periodGwth(100.0,50.0,0.07) == 10.244768351058712

  end

  @testset "pv" begin
    @test FinanceLib.pv(10_000_000., 0.09, 5.0) == 6_499_313.862983453
    @test FinanceLib.pv(12_704_891.6109538, 0.06, 4.0, 12.0) ≈ 10_000_000.
    @test FinanceLib.pvr(10_000_000., 1.09, 5.0) == 6_499_313.862983453
    @test FinanceLib.pvc(11_735.108709918102, 0.08, 2.0) == 10_000
  end

  @testset "fv" begin
    @test FinanceLib.fv(6_499_313.862983453, 0.09, 5.0) == 10_000_000.0
    @test FinanceLib.fv(10_000_000.0, 0.06, 4.0, 12.0) ≈ 12_704_891.6109538
    @test FinanceLib.fvc(10_000., 0.08, 2.0) == 11_735.108709918102
  end

  @testset "annuity" begin
    @test FinanceLib.pvAnnuity(1000.0, 0.12, 5.0) == 3_604.776202345007
    @test FinanceLib.pvAnnuity(7.33764573879378, 0.08, 30.0, 12.0) == 1000

    @test FinanceLib.pvAnnuity(100.0, 0.05) == 2000.0

    @test FinanceLib.fvAnnuity(1000.0, 0.05, 5.0) == 5_525.631250000007
    @test FinanceLib.fvAnnuity(2000.0, 0.24, 5.0, 3.0) == 54_304.2278549568

    @test FinanceLib.pmt(3_604.776202345007, 0.12, 5.0) == 1000.0
    @test FinanceLib.pmt(1000.0, 0.08, 30.0, 12.0) == 7.33764573879378

    @test FinanceLib.fmt(5_525.631250000007, 0.05, 5.0) == 1000
    @test FinanceLib.fmt(54_304.2278549568, 0.24, 5.0, 3.0) == 2000

    @test FinanceLib.pv(FinanceLib.pvAnnuity(10.0^6,.05,30.0),0.05,9.0) == 9_909_218.99605011
  end

  @testset "effective rates" begin
    @test FinanceLib.effR(0.08, 2.0) ≈ 0.0816
    @test FinanceLib.expR(0.08, 2.0) == 0.07844142630656266
    @test FinanceLib.expR(0.08) == 0.0769610411361284
    @test FinanceLib.nomR(FinanceLib.effR(0.08, 4), 4) ≈ 0.08

    FinanceLib.pvc(20,FinanceLib.expR(0.07,4),4.25) == FinanceLib.pvr(20,1+FinanceLib.effR(0.07,4),4.25)

    eT = FinanceLib.RateCurve{FinanceLib.NomRate}([0.0016, 0.0021, 0.0027, 0.0033, 0.0037, 0.0041], 2)

    eR = FinanceLib.effR(eT)

    @test eR.rate[1] ≈ 0.0016006400
    @test eR.rate[2] ≈ 0.0021011025
    @test eR.rate[3] ≈ 0.0027018225
    @test eR.rate[4] ≈ 0.0033027225
    @test eR.rate[5] ≈ 0.0037034225
    @test eR.rate[6] ≈ 0.0041042025

    eN = FinanceLib.nomR(eR)

    @test eN.rate[1] ≈ 0.0016
    @test eN.rate[2] ≈ 0.0021
    @test eN.rate[3] ≈ 0.0027
    @test eN.rate[4] ≈ 0.0033
    @test eN.rate[5] ≈ 0.0037
    @test eN.rate[6] ≈ 0.0041

    eZ = FinanceLib.nomR(FinanceLib.effR(FinanceLib.expR(eT))) # N - X - E - N
    eY = FinanceLib.nomR(FinanceLib.expR(FinanceLib.effR(eT))) # N - E - X - N
    eW = FinanceLib.nomR(FinanceLib.expR(eT)) # N - X - N

    @test eZ.rate[3] == eW.rate[3]
    @test eY.rate[4] == eW.rate[4]
    @test eZ.rate[2] == eY.rate[2]
    @test eZ.rate[5] == eY.rate[5]
    @test eZ.rate[1] == eW.rate[1]

  end

  @testset "npv" begin
    @test FinanceLib.npv(0.05, [0.0:1.0:4.0;],[1000.,2000.0,4000.0,5000.0,6000.0],-1.45)==14709.923338335731
    @test FinanceLib.npv(0.08, [0.25,6.25,3.5,4.5,1.25], [-6.25,1.2,1.25,3.6,2.5], -0.45) == 0.36962283798505946
    @test FinanceLib.npv(0.08, zip([0.25,6.25,3.5,4.5,1.25],[-6.25,1.2,1.25,3.6,2.5]), -0.45) == 0.36962283798505946
    @test FinanceLib.npv(0.08, [0.25,6.25,3.5,4.5,1.25], [-6.25,1.2,1.25,3.6,2.5], 6.25) == 0.619010419015909

    @test FinanceLib.irr([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25]) ≈ 0.31813386476788824

    ts = collect(zip([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25])) :: FinanceLib.PeriodSeries

    @test FinanceLib.irr(ts) ≈ 0.31813386476788824

    @test FinanceLib.irr(zip([0.125,0.29760274,0.49760274,0.55239726,0.812671233], [-10.25,-2.5,3.5,9.5,1.25])) ≈ 0.31813386476788824

    @test FinanceLib.xnpv(0.08, [Dates.Date(2012,2,25), Dates.Date(2012,6,28),
      Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
      [-15, 5, 25, -10, 50], Dates.Date(2012,1,10) ) ==  44.15557928534869

    @test FinanceLib.xnpv(0.08, zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
      Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
      [-15, 5, 25, -10, 50.]), Dates.Date(2012,1,10) ) ==  44.15557928534869

    @test FinanceLib.xirr([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
      Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
      [-115, 5, 25, -10, 200] ) ==  0.2783166029306355

    td = collect(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
    Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
    [-115, 5, 25, -10, 200]))

    td1 = FinanceLib.dateToPeriodSeries(Dates.Date(2010,05,12), td)
    @test td1[3][1] == 2.7652292950034223

    @test FinanceLib.irr(td1) ==  0.2783166029306353
    @test FinanceLib.xirr(td) ==  0.2783166029306355

    @test FinanceLib.xirr(zip([Dates.Date(2012,2,25), Dates.Date(2012,6,28),
      Dates.Date(2013,2,15), Dates.Date(2014,9,18), Dates.Date(2015,2,20)],
      [-115, 5, 25, -10, 200]) ) ==  0.2783166029306355

    @test FinanceLib.npv(FinanceLib.PeriodSeries([(0.5,0.05), (1.25, 0.0575), (2, 0.0485), (3.5, 0.0625), (4.25, 0.055)]), [-150, 20, 15, 80, 100], 0.3) == 31.530253870718543

    @test FinanceLib.xnpv(FinanceLib.DateSeries([(Dates.Date(2014,9,20),0.05), (Dates.Date(2015,2,1), 0.0575), (Dates.Date(2016,10,5), 0.0485), (Dates.Date(2017,12,5), 0.0625), (Dates.Date(2019,1,5), 0.055)]), [-150, 20, 15, 80, 100], Dates.Date(2014,2,15)) ==  29.323165765999597

  end

  @testset "Sharpe" begin
    @test FinanceLib.sharpe(1.58,9.26,22.36) ≈ 0.3434704830053667
  end

  @testset "Rates" begin

    @test FinanceLib.discFactorToNominalRate(FinanceLib.DiscountFactor([0.9524, 0.89, 0.8163, 0.735],1)).rate[3] == 0.0699990723472752

    dsc = FinanceLib.discFactorToNominalRate(FinanceLib.DiscountFactor([ 0.99920063949, 0.99790330288, 0.99596091045, 0.99342713542, 0.99080111671, 0.98778777227 ],2)).rate

    @test dsc[1] ≈ 0.0016
    @test dsc[2] ≈ 0.0021
    @test dsc[3] ≈ 0.0027
    @test dsc[4] ≈ 0.0033
    @test dsc[5] ≈ 0.0037
    @test dsc[6] ≈ 0.0041

    @test FinanceLib.estimR(FinanceLib.RateCurve{FinanceLib.NomRate}([0.05, 0.06, 0.07, 0.08], 2), 1.5) == 0.07
    @test FinanceLib.estimR(FinanceLib.RateCurve{FinanceLib.NomRate}([0.05, 0.06, 0.07, 0.08], 2), 1.2) == 0.064

  end

end

include("FixedIncomes/runtests.jl")
include("Derivatives/runtests.jl")
include("Statements/runtests.jl")


text:
module moreUnits

using Unitful

const radEarth = 6378e3

@unit Re "Re" EarthRadii (radEarth)*Unitful.m false

const localpromotion = Unitful.promotion
function __init__()
    merge!(Unitful.promotion, localpromotion)
    Unitful.register(moreUnits)
end

end

using .moreUnits

functiona(v1::Vector, v2::Vector) = cross(v1, v2)
functiona(c::Real, A::Array) = A*c
functiona(v::Vector, M::Matrix) = M*v

function main()

    c = 2
    v1 = [1, 2, 3]
    v2 = [4, 5, 6]
    M = ones(3,3)

    a = 2*u"Re"
    a = uconvert(u"m", a)
    b = a / (radEarth*u"m")

    println(b)
    println(a)

end

main()

#     a = uconvert(u"m", a)
#     b = a / radEarth
#     println(b)

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
module IO

using DocStringExtensions
using StaticArrays
using Printf
using RecipesBase
using WriteVTK
using OrderedCollections

using WaveProp.Utils
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh
using WaveProp.Nystrom

import WaveProp.Mesh: reference_nodes

export read_geo
export read_msh
export vtk_mesh_file

include("gmshIO.jl")
include("plotsIO.jl")
include("vtkIO.jl")

end # module

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using MassInstallAction
using Test

@testset "MassInstallAction.jl" begin
    @testset "types.jl" begin
        @test MassInstallAction.Workflow("a", Dict("file"=>"stuff"), Set(["bad"])) isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", Dict("file"=>"stuff"), "bad") isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", Dict("file"=>"stuff")) isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", "file"=>"stuff", Set(["bad"])) isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", "file"=>"stuff", "bad") isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", "file"=>"stuff") isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", ["file1"=>"stuff1", "file2"=>"stuff2"], Set(["bad"])) isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", ["file1"=>"stuff1", "file2"=>"stuff2"], "bad") isa MassInstallAction.Workflow
        @test MassInstallAction.Workflow("a", ["file1"=>"stuff1", "file2"=>"stuff2"]) isa MassInstallAction.Workflow
    end

    @testset "public.jl" begin
    end

    @testset "repo_url.jl" begin
        @test MassInstallAction._get_repo_url_with_auth("MYORG", "MYPKG"; token = "MYTOKEN") == "https://x-access-token:MYTOKEN@github.com/MYORG/MYPKG.jl"
        @test MassInstallAction._get_repo_url_with_auth("MYORG", "MYPKG"; token = nothing) == "https://github.com/MYORG/MYPKG.jl"
    end
end