text:
struct PoincareInvariant2ndTrapezoidal{ET,DT,TT,ΩT}
    equ::ET
    ω::ΩT
    Δt::TT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::Vector{DT}
    J::Vector{DT}
    L::Vector{DT}
end

function PoincareInvariant2ndTrapezoidal(f_equ::Function, f_surface::Function, ω::ΩT, Δt::TT, d::Int, nx::Int, ny::Int, ntime::Int, nsave::Int, DT=Float64) where {TT,ΩT}

    println()
    println("Second Euler-Poincaré Integral Invariant")
    println("========================================")
    println()
    println(" nx    = ", nx)
    println(" ny    = ", ny)
    println(" ntime = ", ntime)
    println(" nsave = ", nsave)
    println(" Δt    = ", Δt)
    println()

    # compute initial conditions
    q₀ = zeros(d, nx*ny)

    for j in 1:ny
        for i in 1:nx
            q₀[:,nx*(j-1)+i] = f_surface(i/nx, j/nx)
        end
    end

    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = zeros(DT, nt+1)
    J = zeros(DT, nt+1)
    L = zeros(DT, nt+1)

    PoincareInvariant2ndTrapezoidal{typeof(equ),DT,TT,ΩT}(equ, ω, Δt, nx, ny, ntime, nsave, nt, I, J, L)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant2ndTrapezoidal, sol::Solution)
    for i in 1:size(sol.q.d,2)
        pinv.I[i] = surface_integral(sol.t.t[i], sol.q.d[:,i,:], pinv.ω, pinv.nx, pinv.ny)
        pinv.J[i] = surface_integral_canonical(sol.q.d[:,i,:], sol.p.d[:,i,:], pinv.nx, pinv.ny)
    end

    return (pinv.I, pinv.J, pinv.L)
end

function interpolate_trajectory(x, i1, j1, λ, μ, γ, nx, ny)
    @assert length(γ) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in 1:length(γ)
        γ[k] = x[k, nx*(j1-1)+i1] * (1-λ) * (1-μ) +
               x[k, nx*(j1-1)+i2] *    λ  * (1-μ) +
               x[k, nx*(j2-1)+i1] * (1-λ) *    μ  +
               x[k, nx*(j2-1)+i2] *    λ  *    μ
    end
end

function interpolate_derivative_i(x, i1, j1, λ, μ, γ̇, nx, ny)
    @assert length(γ̇) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in 1:length(γ̇)
        γ̇[k] = (x[k, nx*(j1-1)+i2] - x[k, nx*(j1-1)+i1]) * (1-μ) +
               (x[k, nx*(j2-1)+i2] - x[k, nx*(j2-1)+i1]) *    μ
    end
end

function interpolate_derivative_j(x, i1, j1, λ, μ, γ̇, nx, ny)
    @assert length(γ̇) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in 1:length(γ̇)
        γ̇[k] = (x[k, nx*(j2-1)+i1] - x[k, nx*(j1-1)+i1]) * (1-λ) +
               (x[k, nx*(j2-1)+i2] - x[k, nx*(j1-1)+i2]) *    λ
    end
end

function integrate(t, γ, γ̇ᵢ, γ̇ⱼ, ω, b::Vector{TT}, c::Vector{TT}, q::Vector{DT}, vᵢ::Vector{DT}, vⱼ::Vector{DT}, B::Matrix{DT}) where {DT,TT}
    @assert length(b) == length(c)

    local result = zero(DT)

    for i in 1:length(b)
        for j in 1:length(b)
            γ(c[i], c[j], q)
            γ̇ᵢ(c[i], c[j], vᵢ)
            γ̇ⱼ(c[i], c[j], vⱼ)
            ω(t, q, B)

            result += b[i] * b[j] * vector_matrix_vector_product(vⱼ, B, vᵢ)
        end
    end

    return result
end

function surface_integral(t, x::Matrix{DT}, ω, nx, ny) where {DT}
    local b = [0.5, 0.5]
    local c = [0.0, 1.0]

    local q  = zeros(DT, size(x,1))
    local vᵢ = zeros(DT, size(x,1))
    local vⱼ = zeros(DT, size(x,1))
    local B  = zeros(DT, size(x,1), size(x,1))
    local I  = zero(DT)

    integrate_trapezoidal = (γ, γ̇ᵢ, γ̇ⱼ) -> integrate(t, γ, γ̇ᵢ, γ̇ⱼ, ω, b, c, q, vᵢ, vⱼ, B)

    for j in 1:ny-1
        for i in 1:nx-1
            γ  = (λ, μ, y) -> interpolate_trajectory(x, i, j, λ, μ, y, nx, ny)
            γ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(x, i, j, λ, μ, y, nx, ny)
            γ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(x, i, j, λ, μ, y, nx, ny)
            I += integrate_trapezoidal(γ, γ̇ᵢ, γ̇ⱼ)
        end
    end

    return I
end

function integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b::Vector{TT}, c::Vector{TT}, vᵢ::Vector{DT}, vⱼ::Vector{DT}) where {DT,TT}
    @assert length(b) == length(c)

    local result = zero(DT)

    for i in 1:length(b)
        for j in 1:length(b)
            Θ̇ᵢ(c[i], c[j], vᵢ)
            γ̇ⱼ(c[i], c[j], vⱼ)
            result += b[i] * b[j] * dot(vᵢ,vⱼ)

            γ̇ᵢ(c[i], c[j], vᵢ)
            Θ̇ⱼ(c[i], c[j], vⱼ)
            result -= b[i] * b[j] * dot(vᵢ,vⱼ)
        end
    end

    return result
end

function surface_integral_canonical(q::Matrix{DT}, p::Matrix{DT}, nx, ny) where {DT}
    local b = [0.5, 0.5]
    local c = [0.0, 1.0]

    local vᵢ = zeros(DT, size(q,1))
    local vⱼ = zeros(DT, size(q,1))
    local I  = zero(DT)

    integrate_trapezoidal = (γ, γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ) -> integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b, c, vᵢ, vⱼ)

    for j in 1:ny-1
        for i in 1:nx-1
            γ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(q, i, j, λ, μ, y, nx, ny)
            γ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(q, i, j, λ, μ, y, nx, ny)
            Θ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(p, i, j, λ, μ, y, nx, ny)
            Θ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(p, i, j, λ, μ, y, nx, ny)
            I += integrate_trapezoidal(γ, γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ)
        end
    end

    return I
end

function CommonFunctions.write_to_hdf5(pinv::PoincareInvariant2ndTrapezoidal, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing
        isdefined(sol, :λ) ? write(h5, "L", pinv.L) : nothing

    end
end


text:
module SparseCntr

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOSmoothMcCormickGrad

nx = 3
A = spzeros(Float64,nx,nx)
c = [3.0, 0.5, 1.0]
P = spzeros(Float64,nx,nx)
A[1,1] = P[1,1] = 0.5
A[1,2] = P[1,2] = 0.2
A[1,3] = P[1,3] = 0.1
A[2,1] = P[2,1] = 0.3
A[2,2] = P[2,2] = 0.5
A[2,3] = P[2,3] = 0.5
A[3,1] = P[3,1] = 0.4
A[3,2] = P[3,2] = 0.5
A[3,3] = P[3,3] = 0.3

lufct = lufact(P)
At = copy(A)
ct = copy(c)
Pt = copy(P)
preA = inv(lufct)*A
prec = inv(lufct)*c
SSto = EAGOSmoothMcCormickGrad.SparseInSto()
SSto.Xh = zeros(Float64,nx)
SSto.Yh = zeros(Float64,nx)
SSto.Zh = zeros(Float64,nx)
SSto.Xj = spzeros(Float64,nx,nx)
SSto.Yj = spzeros(Float64,nx,nx)
SSto.Zj = spzeros(Float64,nx,nx)
SSto.nx = nx
EAGOSmoothMcCormickGrad.Sparse_Precondition!(c,A,P,SSto)

@test (7.7272-10E-4 <= c[1] <= 7.7272+10E-4)
@test (-5.0-10E-4 <= c[2] <= -5.0+10E-4)
@test (1.3636-10E-4 <= c[3] <= 1.3636+10E-4)

end

ext: jl


text:
using RequiredKeywords
using Base.Test

# Assignment form
@required_keywords f(a; x, y::Int, z::Int=2) = a * x * y * z

# Standard form
@required_keywords function g(a; x, y::Int, z::Int=2)
    a * x * y * z
end


text:
module QuantumStateDistributions
    using ClassicalOrthogonalPolynomials
    using ProgressMeter
    using QuantumStateBase
    using Distributions
    using LinearAlgebra
    using Random

    include("gaussian_state_bhd.jl")
    include("pdf.jl")
    include("mle.jl")
end


text:
function writemaster(master,filename)

ndims = zeros(20,1);
ndims[1] = master.nd;
ndims[2] = master.elemtype;
ndims[3] = master.nodetype;
ndims[4] = master.porder;
ndims[5] = master.pgauss;
ndims[6] = master.npe;
ndims[7] = master.npf;
ndims[8] = master.nge;
ndims[9] = master.ngf;
ndims[10] = length(master.xp1d);
ndims[11] = length(master.gp1d);

nsize = zeros(22,1);
nsize[1] = length(ndims[:]);
nsize[2] = length(master.shapegt[:]);
nsize[3] = length(master.shapegw[:]);
nsize[4] = length(master.shapfgt[:]);
nsize[5] = length(master.shapfgw[:]);
nsize[6] = length(master.shapent[:]);
nsize[7] = length(master.shapen[:]);
nsize[8] = length(master.shapfnt[:]);
nsize[9] = length(master.shapfn[:]);
nsize[10] = length(master.xpe[:]);
nsize[11] = length(master.gpe[:]);
nsize[12] = length(master.gwe[:]);
nsize[13] = length(master.xpf[:]);
nsize[14] = length(master.gpf[:]);
nsize[15] = length(master.gwf[:]);
nsize[16] = length(master.shap1dgt[:]);
nsize[17] = length(master.shap1dgw[:]);
nsize[18] = length(master.shap1dnt[:]);
nsize[19] = length(master.shap1dn[:]);
nsize[20] = length(master.xp1d[:]);
nsize[21] = length(master.gp1d[:]);
nsize[22] = length(master.gw1d[:]);

print("Writing master into file...\n");
fileID = open(filename,"w");

# write master structure to files
write(fileID,Float64(length(nsize[:])));
write(fileID,Float64.(nsize[:]));
write(fileID,Float64.(ndims[:]));
write(fileID,master.shapegt[:]);
write(fileID,master.shapegw[:]);
write(fileID,master.shapfgt[:]);
write(fileID,master.shapfgw[:]);
write(fileID,master.shapent[:]);
write(fileID,master.shapen[:]);
write(fileID,master.shapfnt[:]);
write(fileID,master.shapfn[:]);
write(fileID,master.xpe[:]);
write(fileID,master.gpe[:]);
write(fileID,master.gwe[:]);
write(fileID,master.xpf[:]);
write(fileID,master.gpf[:]);
write(fileID,master.gwf[:]);
write(fileID,master.shap1dgt[:]);
write(fileID,master.shap1dgw[:]);
write(fileID,master.shap1dnt[:]);
write(fileID,master.shap1dn[:]);
write(fileID,master.xp1d[:]);
write(fileID,master.gp1d[:]);
write(fileID,master.gw1d[:]);

close(fileID);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
export  Obstacles,
        Particle,
        Event,
        Metadata

mutable struct Obstacles # We will use this as a struct of arrays
    x::Array{Float64,1}
    y::Array{Float64,1}
    r::Array{Float64,1}
    n::Int64
end

mutable struct Metadata
    nx::Int64
    ny::Int64
    edgesx::Array{Float64,1}
    edgesy::Array{Float64,1}
end

mutable struct Particle
    x::Float64
    y::Float64
    r::Float64
    vx::Float64
    vy::Float64
end

mutable struct Event
    oldx::Float64
    oldy::Float64
    oldvx::Float64
    oldvy::Float64
    oldi::Int64
    oldj::Int64
    newx::Float64
    newy::Float64
    newvx::Float64
    newvy::Float64
    newi::Int64
    newj::Int64
    distance::Float64
    collision::Bool
    crossinterfaceperiodic::Bool
end


text:
function compute_derivatives(obj::FieldSecondaryODE, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement}, arg3::Vector{RealFieldElement})
    return jcall(obj, "computeDerivatives", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}, Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function get_dimension(obj::FieldSecondaryODE)
    return jcall(obj, "getDimension", jint, ())
end

function init(obj::FieldSecondaryODE, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement}, arg3::RealFieldElement)
    return jcall(obj, "init", void, (RealFieldElement, Vector{RealFieldElement}, Vector{RealFieldElement}, RealFieldElement), arg0, arg1, arg2, arg3)
end


text:
if R123_USE_AESNI

import Random: seed!
using Test: @test

using RandomNumbers
using Random123

import RandomNumbers: split_uint
import Random123: __m128i, AESNIKey

x = zero(__m128i)
ctr = __m128i(0x9799b5d54f7b9227, 0xb47607190d0dfefb)
key = 0x07b8e4b6aa98ec245a7da274d3b8146a
aesni_key = AESNIKey(key)
@test rand(AESNI1x(x, ctr, aesni_key), UInt128) ≡ 0x60f4c27fe48fe1b8c5f4568a585b0dc0

r = AESNI1x(key)
r1 = AESNI4x(split_uint(key, UInt32))
@test seed_type(r) ≡ UInt128
@test seed_type(r1) ≡ NTuple{4, UInt32}
@test copyto!(copy(r), r) == r
@test copyto!(copy(r1), r1) == r1
@test UInt128(r.x) ≡ rand(r1, UInt128)
@test rand(r, UInt128) ≡ rand(r1, UInt128)
set_counter!(r, 0)
set_counter!(r1, 1)
@test rand(r, Tuple{UInt128})[1] ≡ rand(r1, UInt128)

end


text:
module BlockSparseMatrices

using BlockArrays

import BlockArrays: nblocks, getblock, getblock!, setblock!, BlockIndex

import Base: Order.Forward, A_mul_B!, SparseMatrixCSC

export SparseMatrixBSC, nnzblocks, nblocks

include("SparseMatrixBSC.jl")

end # module


text:
using MPIReco, JLD, SensorSelection, PyPlot

# estimate measurement noise from the training data
function varMeas(x::Matrix{T},rs::Int64) where T
  roi = vec(x[1:rs,1:rs].-mean(x[1:rs,1:rs]))
  append!(roi,x[end-rs+1:end,1:rs].-mean(x[end-rs+1:end,1:rs]))
  append!(roi,x[1:rs,end-rs+1:end].-mean(x[1:rs,end-rs+1:end]))
  append!(roi,x[end-rs+1:end,end-rs+1:end].-mean(x[end-rs+1:end,end-rs+1:end]))

  return sum(abs.(vec(roi)).^2)/(length(roi)-1)
end

####################
# load system matrix
####################
# Dataset store handling
# uncomment the next two lines if running this script separately
# datadir = artifact"MDFStore"
# store = MDFDatasetStore(datadir)

# system matrix
datadirSF1 = joinpath(calibdir(store), "1.mdf")
fSF1 = MPIFile(datadirSF1)
nx,ny,nz = calibSize(fSF1)
S1 = reshape( getSystemMatrix(fSF1,bgCorrection=true), nx, ny, :)
snr1 = MPIReco.getSNRAllFrequencies(fSF1)

# select a slice and apply SNR thresholding
snr_idx1 = findall(x->x>3.0, vec(snr1))[2:end]
S1 = S1[:,:,snr_idx1]

# normalize the frequency components
for k=1:size(S1,3)
  S1[:,:,k] ./= norm(S1[:,:,k])
end

# chose frequency components for training
Random.seed!(1234)
p = shuffle(collect(4:2:size(S1,3)))
St = S1[:,:,p[1:30]]  # training data

########################################################
# setup up measurement operators and covariance matrices
########################################################

# mean amplitudes of measurements
st_amp = zeros(nx,ny,size(St,3))
for k=1:size(St,3)
  st_amp[:,:,k] .= sqrt.( 1/size(St,3)*dropdims(sum(abs.(St[:,:,:]).^2,dims=3),dims=3) )
end
st_amp = reshape(st_amp,nx*ny,:)

# measurement matrix
dctOp = DCTOp(ComplexF64,(nx,ny))
H = Matrix(adjoint(dctOp))
Ht = collect(transpose(real.(H))) # build transpose of measurement matrix

# DCT representation of training data
St_dct = zeros(ComplexF64,nx*ny,size(St,3));
for k=1:size(St,3)
    St_dct[:,k] .= dctOp*vec(St[:,:,k])
end

# obtain sparse measurement matrix
thresh = 1.e-1  # portion of the total number of sparse coefficients to use for optimization
numCoeff = floor(Int64,thresh*nx*ny)
Ht_sub = zeros(numCoeff,size(Ht,2),size(St,3))
for k=1:size(St,3)
  Sk_dct = vec(abs.(St_dct[:,k]))
  idx_sig = maxIndices(Sk_dct,numCoeff)
  Ht_sub[:,:,k] .= Ht[idx_sig,:]
end

# estimate measurement noise for each frequency component
measNoise = [varMeas(St[:,:,k],3) for k=1:size(St,3)]
r = zeros(nx*ny,size(St,3))
for k=1:size(St,3)
  r[:,k] .= fill(measNoise[k],nx*ny)
end

######################
# Forward optimization
######################
λ=3.0                             # proportionality factor for the prior covariance matrices
numSamples = floor(Int64, nx*ny)  # number of samples to add to the sampling pattern
w = zeros(Bool,nx*ny)             # initial sampling pattern (empty)
m = SensorSelection.Experiment(Ht_sub, r, Ht_sub, λ.*st_amp.^2,w)
@time idx_opt = optSamplingFW!(m,numSamples,4,180)

#############################
# plot some sampling patterns
#############################
# precomputed Poisson disk patterns
msk2_pd = zeros(nx,ny); msk2_pd[load(datadir*"/samplingPatterns/pd_r2.jld","idx")] .= 1
msk6_pd = zeros(nx,ny); msk6_pd[load(datadir*"/samplingPatterns/pd_r6.jld","idx")] .= 1
msk10_pd = zeros(nx,ny); msk10_pd[load(datadir*"/samplingPatterns/pd_r10.jld","idx")] .= 1
# optimized patterns
msk2_opt = zeros(nx,ny); msk2_opt[idx_opt[1:div(nx*ny,2)]] .= 1
msk6_opt = zeros(nx,ny); msk6_opt[idx_opt[1:div(nx*ny,6)]] .= 1
msk10_opt = zeros(nx,ny); msk10_opt[idx_opt[1:div(nx*ny,10)]] .= 1

figure("sampling Patterns", figsize=(8,6))
subplot(2,3,1); imshow(msk2_pd); title("r=2"); ylabel("PD"); xticks([]); yticks([])
subplot(2,3,2); imshow(msk6_pd); title("r=6"); xticks([]); yticks([])
subplot(2,3,3); imshow(msk10_pd); title("r=10"); xticks([]); yticks([])
subplot(2,3,4); imshow(msk2_opt); ylabel("Opt"); xticks([]); yticks([])
subplot(2,3,5); imshow(msk6_opt); xticks([]); yticks([])
subplot(2,3,6); imshow(msk10_opt); xticks([]); yticks([])
tight_layout()


text:
using Test
using LaTeXStrings: latexstring
using Printf: @printf, @sprintf
using QuantumLattices.Essentials.QuantumOperators
using QuantumLattices.Essentials: dtype
using QuantumLattices.Interfaces: id, value, rank, add!, sub!, mul!, div!
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Prerequisites.Combinatorics: Combinations
using QuantumLattices.Prerequisites.Traits: contentnames, parameternames, parametertype, isparameterbound

import QuantumLattices.Interfaces: ⊗, ⋅, permute
import QuantumLattices.Essentials.QuantumOperators: script

struct AID{O<:Real, S<:Real} <: OperatorUnit
    orbital::O
    nambu::S
end
@inline Base.adjoint(id::AID) = AID(id.orbital, 3-id.nambu)
@inline ⊗(m₁::Operator{<:Number, <:ID{AID}}, m₂::Operator{<:Number, <:ID{AID}}) = m₁ * m₂
@inline ⋅(m₁::Operator{<:Number, <:ID{AID}}, m₂::Operator{<:Number, <:ID{AID}}) = m₁ * m₂
@inline script(::Val{:orbital}, id::AID; kwargs...) = id.orbital
@inline script(::Val{:nambu}, id::AID; kwargs...) = id.nambu==2 ? "\\dagger" : ""
latexformat(AID, LaTeX{(:nambu,), (:orbital,)}('c'))

function permute(u₁::AID, u₂::AID)
    @assert u₁ ≠ u₂ "permute error: permuted operator units should not be equal to each other."
    if (u₁.nambu == 3-u₂.nambu) && (u₁.orbital == u₂.orbital)
        if u₁.nambu == 2
            return (Operator(1), Operator(1, ID(u₂, u₁)))
        else
            return (Operator(-1), Operator(1, ID(u₂, u₁)))
        end
    else
        return (Operator(1, ID(u₂, u₁)),)
    end
end

@testset "ID" begin
    @test ID|>rank == Any
    @test promote_type(ID{AID{Int, Int}, 1}, ID{AID{Int, Int}}) == ID{AID{Int, Int}}
    @test promote_type(ID{AID{Int, Int}}, ID{AID{Int, Int}, 1}) == ID{AID{Int, Int}}
    @test promote_type(ID{AID{Int, Int}}, ID{AID{Int, Float}}) == ID{AID{Int, <:Real}}
    @test promote_type(Tuple{}, ID{AID{Int, Int}}) == ID{AID{Int, Int}}
    @test promote_type(Tuple{}, ID{AID{Int, Int}, 2}) == ID{AID{Int, Int}}

    sid = AID(1, 1)
    @test string(sid) == "AID(1, 1)"
    @test sid'==AID(1, 2)

    cid = ID(AID(2, 1), AID(Inf, 1))
    @test cid == ID((AID(2, 1), AID(Inf, 1)))
    @test ID(sid, cid) == ID(AID(1, 1), AID(2, 1), AID(Inf, 1))
    @test ID(cid, sid) == ID(AID(2, 1), AID(Inf, 1), AID(1, 1))
    @test ID(cid, cid) == ID(AID(2, 1), AID(Inf, 1), AID(2, 1), AID(Inf, 1))
    @test cid == ID(AID, (2, Inf), (1, 1))
    @test cid|>propertynames == (:orbitals, :nambus)
    @test cid.orbitals == (2, Inf)
    @test cid.nambus == (1, 1)
    @test cid|>string == "ID(AID(2, 1), AID(Inf, 1))"
    @test cid|>eltype == AID{<:Real, Int}
    @test cid|>rank == cid|>typeof|>rank == 2
    @test cid'==ID(AID(Inf, 2), AID(2, 2))
    @test ishermitian(cid)==false
    @test ishermitian(ID(cid, cid'))
end

@testset "QuantumOperator" begin
    @test contentnames(OperatorProd) == (:value, :id)
    @test parameternames(OperatorProd) == (:value, :id)
    @test isparameterbound(OperatorProd, :value, Any) == false
    @test isparameterbound(OperatorProd, :id, ID{AID{Int, Int}, 2}) == false

    @test valtype(OperatorProd) == parametertype(OperatorProd, :value) == parametertype(OperatorProd, 1) == Any
    @test valtype(OperatorProd{Int}) == parametertype(OperatorProd{Int}, :value) == parametertype(OperatorProd{Int}, 1) == Int
    @test idtype(OperatorProd{<:Number}) == parametertype(OperatorProd{<:Number}, :id) == parametertype(OperatorProd{<:Number}, 2) == ID{OperatorUnit}
    @test idtype(OperatorProd{<:Number, ID{AID}}) == parametertype(OperatorProd{<:Number, ID{AID}}, :id) == parametertype(OperatorProd{<:Number, ID{AID}}, 2) == ID{AID}

    @test promote_type(Operator{Int}, Operator) == Operator
    @test promote_type(Operator, Operator{Int}) == Operator
    @test promote_type(Operator{Int}, Operator{Float}) == Operator{Float}
    @test promote_type(Operator{Int, ID{AID{Int, Int}, 2}}, Operator{Float, ID{AID{Int, Int}, 2}}) == Operator{Float, ID{AID{Int, Int}, 2}}
    @test promote_type(Operator{Int, ID{AID}}, Operator{Float, ID{AID}}) == Operator{Float, <:ID{AID}}
    @test promote_type(Operator{Int}, Float) == Operator{Float}

    opt = Operator(2.0, AID(1, 1))
    @test deepcopy(opt)==opt && isequal(deepcopy(opt), opt)
    @test eltype(opt) == eltype(typeof(opt)) == AID{Int, Int}
    @test collect(opt) == [AID(1, 1)]
    @test opt|>rank == opt|>typeof|>rank == 1
    @test opt|>valtype == opt|>typeof|>valtype == parametertype(opt|>typeof, :value) == Float
    @test opt|>idtype == opt|>typeof|>idtype == parametertype(opt|>typeof, :id) == ID{AID{Int, Int}, 1}
    @test opt|>dtype == opt|>typeof|>dtype == Float
    @test value(opt) == 2.0
    @test id(opt) == ID(AID(1, 1))
    @test replace(opt, 3) == replace(opt, value=3) == Operator(3, AID(1, 1))
    @test isapprox(opt, replace(opt, value=opt.value+10^-6); atol=10^-5)
    @test length(opt) == 1 && firstindex(opt) == 1 && lastindex(opt) ==1
    @test opt[1]==AID(1, 1) && opt[1:1]==Operator(1.0, AID(1, 1))
    @test split(opt) == (2.0, AID(1, 1))
    @test opt|>typeof|>one == opt|>one == Operator(1.0)
    @test convert(Operator{Float, Tuple{AID{Int, Int}}}, Operator(2, AID(1, 1))) == Operator(2.0, AID(1, 1))
    @test Operator(2.0) == Operator(2.0)
    @test string(opt) == "Operator(2.0, AID(1, 1))"
    @test opt' == Operator(2.0, AID(1, 2))
    @test ishermitian(opt) == false
    @test ishermitian(Operator(2.0, AID(1, 1), AID(1, 2)))
    @test sequence(opt, Dict(AID(1, 1)=>1, AID(1, 2)=>2)) == (1,)
    @test convert(typeof(opt), AID(1, 1)) == Operator(1.0, AID(1, 1))
    @test convert(Operator{Float, <:ID{AID}}, AID(1, 1)) == Operator(1.0, AID(1, 1))
    @test convert(Operator{Float, Tuple{}}, 2) == Operator(2.0)
    @test convert(Operator{Float, <:ID{AID}}, 2) == Operator(2.0)

    opt₁ = Operator(1.0, AID(1, 1))
    opt₂ = Operator(2.0, AID(1, 2))
    opts = Operators(opt₁, opt₂)
    @test opts == Operators{eltype(opts)}(opt₁, opt₂)
    @test opts == OperatorSum(opt₁, opt₂) == OperatorSum((opt₁, opt₂)) == OperatorSum{eltype(opts)}(opt₁, opt₂) == OperatorSum{eltype(opts)}((opt₁, opt₂))
    @test eltype(opts) == eltype(typeof(opts)) == Operator{Float, ID{AID{Int, Int}, 1}}
    @test collect(opts) == collect(values(opts.contents))
    @test length(opts) == 2
    @test summary(opts) == "Operators"
    @test string(opts) == @sprintf "Operators with 2 Operator\n  %s\n" join(opts, "\n  ")
    @test haskey(opts, id(opt₁)) && haskey(opts, id(opt₂)) && !haskey(opts, ID(AID(3, 1)))
    @test opts[id(opt₁)]==opt₁ && opts[id(opt₂)]==opt₂
    @test empty(opts) == empty!(deepcopy(opts)) == zero(opts)
    @test merge(empty(opts), opts) == merge!(empty(opts), opts) == opts
    optp₁ = promote_type(typeof(opts), OperatorSum{Operator{Complex{Int}, NTuple{2, AID{Float, Float}}}, NTuple{2, AID{Float, Float}}})
    optp₂ = OperatorSum{Operator{Complex{Float}, <:Tuple{AID, Vararg{AID{Float, Float}}}}, Tuple{AID, Vararg{AID{Float, Float}}}}
    @test optp₁ == optp₂
    @test isapprox(opts, deepcopy(opts)) && isapprox(opts, opts+10^-6; atol=10^-5) && isapprox(opts+10^-6, opts; atol=10^-5)
    @test !isapprox(opts, opts+10^-6) && !isapprox(opts+10^-6, opts)
    @test opts' == Operators(opt₁', opt₂')
    @test !ishermitian(opts) && ishermitian(Operators(opt₁, opt₂, opt₁', opt₂'))
    @test zero(opts) == zero(typeof(opts)) == OperatorSum{eltype(opts)}()
    @test add!(deepcopy(opts))==opts && add!(zero(opts), opts)==opts
    @test sub!(deepcopy(opts))==opts && sub!(deepcopy(opts), opts)==zero(opts)
    @test mul!(deepcopy(opts), 2.0) == opts*2
    @test div!(deepcopy(opts), 2.0) == opts/2

    @test optype(AID(1, 1)) == optype(AID{Int, Int}) == Operator{Int, Tuple{AID{Int, Int}}}
    @test optype(AID) == Operator{Int, <:Tuple{AID}}
    @test optype(opt) == optype(typeof(opt)) == Operator{Float, Tuple{AID{Int, Int}}}
    @test optype(opts) == optype(typeof(opts)) == eltype(opts)

    @test +opt == opt
    @test -opt == Operator(-2.0, AID(1, 1))
    @test opt*2 == 2*opt == Operator(4.0, AID(1, 1))
    @test opt/2 == Operator(1.0, AID(1, 1))
    @test opt^2 == opt*opt == Operator(4.0, AID(1, 1), AID(1, 1))
    @test opt+1 == 1+opt == Operators(Operator(1), opt)
    @test 1-opt == Operators(Operator(1), -opt)
    @test opt-1 == Operators(opt, Operator(-1))
    @test 2*AID(1, 1) == AID(1, 1)*2 == Operator(2, AID(1, 1))
    @test AID(1, 1)*AID(1, 2) == Operator(1, AID(1, 1), AID(1, 2))
    @test Operator(2, AID(1, 1))//3 == Operator(2//3, AID(1, 1))

    @test +opts == opts
    @test -opts == Operators(-opt₁, -opt₂)
    @test opts*2 == 2*opts == Operators(2opt₁, 2opt₂)
    @test opts/2 == Operators(opt₁/2, opt₂/2)
    @test opts^2 == opts*opts == opts⊗opts == opts⋅opts == Operators(opt₁*opt₁, opt₁*opt₂, opt₂*opt₁, opt₂*opt₂)
    @test opts+1 == 1+opts == Operators(Operator(1), opt₁, opt₂)
    @test 1-opts == Operators(Operator(1), -opt₁, -opt₂)
    @test opts-1 == Operators(opt₁, opt₂, Operator(-1))
    @test opt₁+opt₂ == opts
    @test opt₁+opts == opts+opt₁ == Operators(opt₁*2, opt₂)
    @test opts+opts == Operators(opt₁*2, opt₂*2)
    @test opts-opt₁ == Operators(opt₂)
    @test opt₁-opts == Operators(-opt₂)
    @test opts-opts == zero(opts) == zero(typeof(opts))
    @test opts*opt == Operators(opt₁*opt, opt₂*opt) == opts⊗opt == opts⋅opt
    @test opt*opts == Operators(opt*opt₁, opt*opt₂) == opt⊗opts == opt⋅opts
end

@testset "LaTeX" begin
    @test latexname(AID) == :AID
    @test latexformat(AID) == LaTeX{(:nambu,), (:orbital,)}('c')

    latex = LaTeX{(:nambu,), (:orbital,)}('d')
    @test superscript(latex|>typeof) == (:nambu,)
    @test subscript(latex|>typeof) == (:orbital,)
    latexformat(AID, latex)

    aid = AID(1, 2)
    @test script(Val(:BD), aid, latex) == 'd'
    @test script(Val(:SP), aid, latex) == ("\\dagger",)
    @test script(Val(:SB), aid, latex) == (1,)
    @test latexstring(aid) == "d^{\\dagger}_{1}"

    opt = Operator(1.0, AID(2, 2), AID(1, 1))
    @test latexstring(opt) == "d^{\\dagger}_{2}d^{}_{1}"
    io = IOBuffer()
    show(io, MIME"text/latex"(), opt)
    @test String(take!(io)) == "\$d^{\\dagger}_{2}d^{}_{1}\$"
    @test latexstring(Operator(1.0, aid, aid)) == "(d^{\\dagger}_{1})^2"

    latexformat(AID, LaTeX{(:nambu,), (:orbital,)}('c'))
    opts = Operators(
            Operator(1.0-1.0im, AID(2, 2), AID(1, 1)),
            Operator(-1.0, AID(1, 2), AID(1, 1))
            )
    candidate₁ = "(1.0-1.0im)c^{\\dagger}_{2}c^{}_{1}-c^{\\dagger}_{1}c^{}_{1}"
    candidate₂ = "-c^{\\dagger}_{1}c^{}_{1}+(1.0-1.0im)c^{\\dagger}_{2}c^{}_{1}"
    @test latexstring(opts) ∈ (candidate₁, candidate₂)
end

struct UnitCoeffAddition <: Transformation end
Base.valtype(::Type{UnitCoeffAddition}, M::Type{<:Union{Operator, Operators}}) = M
@inline (unitcoeffaddition::UnitCoeffAddition)(m::Operator) = replace(m, value=one(valtype(m))+value(m))

@testset "Transformation" begin
    m = Operator(1, ID(AID(1, 1)))
    s = Operators(m)

    i = Identity()
    @test i==deepcopy(i) && isequal(i, deepcopy(i))
    @test valtype(i, m) == valtype(typeof(i), typeof(m)) == typeof(m)
    @test valtype(i, s) == valtype(typeof(i), typeof(s)) == typeof(s)
    @test i(m)==m && i(s)==s

    n = Numericalization{Float64}()
    @test valtype(n) == valtype(typeof(n)) == Float64
    @test valtype(n, m) == valtype(typeof(n), typeof(m)) == Operator{Float, Tuple{AID{Int, Int}}}
    @test valtype(n, s) == valtype(typeof(n), typeof(s)) == OperatorSum{Operator{Float, Tuple{AID{Int, Int}}}, Tuple{AID{Int, Int}}}
    @test n(m) == replace(m, value=1.0)
    @test n(s) == Operators(replace(m, value=1.0))

    unitcoeffaddition = UnitCoeffAddition()
    @test map!(unitcoeffaddition, s) == s == Operators(Operator(2, ID(AID(1, 1))))
    @test map!(unitcoeffaddition, zero(s), s) == Operators(Operator(3, ID(AID(1, 1))))
end

@testset "Permutation" begin
    id₁, id₂ = AID(1, 1), AID(1, 2)
    opt₀ = Operator(1.5)
    opt₁ = Operator(1.5, id₁, id₂)
    opt₂ = Operator(1.5, id₂, id₁)

    p = Permutation(Dict(id₁=>1, id₂=>2))
    @test p(opt₁) == opt₂-opt₀
    @test p(opt₁+opt₂) == 2*opt₂-opt₀

    p = Permutation(Dict(id₁=>2, id₂=>1))
    @test p(opt₂) == opt₁+opt₀
    @test p(opt₁+opt₂) == 2*opt₁+opt₀
end

@testset "Substitution" begin
    id₁, id₂ = AID(1, 1), AID(1, 2)
    opt = Operator(1.5, id₁, id₂)

    id₃, id₄ = AID(2, 1), AID(2, 2)
    ops₁ = Operator(2.0, id₃) + Operator(3.0, id₄) + 1
    ops₂ = Operator(2.0, id₃) - Operator(3.0, id₄)

    M = promote_type(eltype(ops₁), eltype(ops₂))
    table = Dict{eltype(opt), Operators{M, idtype(M)}}(id₁=>ops₁, id₂=>ops₂)
    substitution = UnitSubstitution(table)
    @test substitution(opt) == 1.5*ops₁*ops₂
    @test substitution(Operators(opt)) == 1.5*ops₁*ops₂
end

@testset "RankFilter" begin
    op₀, op₁, op₂ = Operator(1), Operator(2, AID(1, 1)), Operator(3, AID(1, 1), AID(2, 1))
    ops = Operators(op₀, op₁, op₂)
    @test RankFilter(0)(ops) == Operators(op₀)
    @test RankFilter(1)(ops) == Operators(op₁)
    @test RankFilter(2)(ops) == Operators(op₂)
end


text:
function change_point_test_prep(N=10000, λ=N/10)
    L = @SMatrix [1. 0.;
                  0. 1.]
    Σ = @SMatrix [10^(-5) 0.;
         0. 10^(-5)]

    θ₀ = [10.0, -8.0, 15.0, 0.0, 3.0]

    param = :complexConjug
    # Target law
    P˟ = FitzhughDiffusion(param, θ₀...)

    # Auxiliary law
    t₀ = 1.0
    T = 2.0
    x0 = ℝ{2}(-0.5, 2.25)
    xT = ℝ{2}(1.0, 0.0)
    P̃ = FitzhughDiffusionAux(param, θ₀..., t₀, L*x0, T, L*xT)

    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    dt = (T-t₀)/N
    tt = τ(t₀,T).(t₀:dt:T)

    P₁ = GuidPropBridge(eltype(x0), tt, P˟, P̃, L, L*x0, Σ;
                        changePt=NoChangePt(), solver=Vern7())

    P₂ = GuidPropBridge(eltype(x0), tt, P˟, P̃, L, L*x0, Σ;
                        changePt=SimpleChangePt(λ), solver=Vern7())
    P₁, P₂
end

@testset "change point between ODE solvers" begin

    parametrisation = POSSIBLE_PARAMS[5]
    include(joinpath(SRC_DIR, "fitzHughNagumo.jl"))
    include(joinpath(SRC_DIR, "types.jl"))
    include(joinpath(SRC_DIR, "vern7.jl"))
    include(joinpath(SRC_DIR, "guid_prop_bridge.jl"))
    N = 10000
    P₁, P₂ = change_point_test_prep(N)
    @testset "comparing H[$i]" for i in 1:div(N,20):N
        @test P₁.H[i] ≈ P₂.H[i]
    end
    @testset "comparing Hν[$i]" for i in 1:div(N,20):N
        @test P₁.Hν[i] ≈ P₂.Hν[i]
    end
    @testset "comparing c[$i]" for i in 1:div(N,20):N
        @test P₁.c[i] ≈ P₂.c[i]
    end
end


text:
using DiffEqSensitivity, SafeTestsets
using Test, Pkg

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

@time begin

if GROUP == "All" || GROUP == "Core1" || GROUP == "Downstream"
    @time @safetestset "Forward Sensitivity" begin include("forward.jl") end
    @time @safetestset "Adjoint Sensitivity" begin include("adjoint.jl") end
    @time @safetestset "Sparse Adjoint Sensitivity" begin include("sparse_adjoint.jl") end
    @time @safetestset "Second Order Sensitivity" begin include("second_order.jl") end
    @time @safetestset "Concrete Solve Derivatives" begin include("concrete_solve_derivatives.jl") end
    @time @safetestset "Branching Derivatives" begin include("branching_derivatives.jl") end
    @time @safetestset "Derivative Shapes" begin include("derivative_shapes.jl") end
    @time @safetestset "save_idxs" begin include("save_idxs.jl") end
    @time @safetestset "ArrayPartitions" begin include("array_partitions.jl") end
    @time @safetestset "Complex Adjoints" begin include("complex_adjoints.jl") end
    @time @safetestset "Forward Remake" begin include("forward_remake.jl") end
    @time @safetestset "DiscreteProblem Adjoints" begin include("discrete.jl") end
    @time @safetestset "Time Type Mixing Adjoints" begin include("time_type_mixing.jl") end
    @time @safetestset "GSA tests" begin include("gsa.jl") end
end

if GROUP == "All" || GROUP == "Core2"
    @time @safetestset "hasbranching" begin include("hasbranching.jl") end
    @time @safetestset "Literal Adjoint" begin include("literal_adjoint.jl") end
    @time @safetestset "ForwardDiff Chunking Adjoints" begin include("forward_chunking.jl") end
    @time @safetestset "Stiff Adjoints" begin include("stiff_adjoints.jl") end
    @time @safetestset "Autodiff Events" begin include("autodiff_events.jl") end
    @time @safetestset "Null Parameters" begin include("null_parameters.jl") end
    @time @safetestset "Forward Mode Prob Kwargs" begin include("forward_prob_kwargs.jl") end
    @time @safetestset "Steady State Adjoint" begin include("steady_state.jl") end
    @time @safetestset "Concrete Solve Derivatives of Second Order ODEs" begin include("second_order_odes.jl") end
end

if GROUP == "All" || GROUP == "Core3"
    @time @safetestset "Shadowing Tests" begin include("shadowing.jl") end
end

if GROUP == "All" || GROUP == "SDE1"
    @time @safetestset "SDE Adjoint" begin include("sde_stratonovich.jl") end
    @time @safetestset "SDE Scalar Noise" begin include("sde_scalar_stratonovich.jl") end
    @time @safetestset "SDE Checkpointing" begin include("sde_checkpointing.jl") end
end

if GROUP == "All" || GROUP == "SDE2"
    @time @safetestset "SDE Non-Diagonal Noise" begin include("sde_nondiag_stratonovich.jl") end
end

if GROUP == "All" || GROUP == "SDE3"
    @time @safetestset "RODE Tests" begin include("rode.jl") end
    @time @safetestset "SDE Ito Conversion Tests" begin include("sde_transformation_test.jl") end
    @time @safetestset "SDE Ito Scalar Noise" begin include("sde_scalar_ito.jl") end
end

if GROUP == "DiffEqFlux"
    activate_downstream_env()
    @time @safetestset "Callback - ReverseDiff" begin include("downstream/callback_reversediff.jl") end
    @time @safetestset "ForwardDiff Sparsity Components" begin include("downstream/forwarddiffsensitivity_sparsity_components.jl") end
    @time @safetestset "SDE - Neural" begin include("downstream/sde_neural.jl") end
    @time @safetestset "Complex No u" begin include("downstream/complex_no_u.jl") end
    @time @safetestset "HybridNODE" begin include("downstream/HybridNODE.jl") end
end

if GROUP == "Callbacks"
    @time @safetestset "Callbacks with Adjoints" begin include("callbacks.jl") end
end

if GROUP == "GPU"
    activate_downstream_env()
    @time @safetestset "Mixed GPU/CPU" begin include("gpu/mixed_gpu_cpu_adjoint.jl") end
end
end


text:
module BBVI

"""TODO: Remove this function."""
greet(name::String) = "Hello, $(name)!"

end # module


text:
################################################################################
#
#   METHODS FOR CONSTRUCTION OF FERMI SURFACES (2D & 3D)
#
#   STRUCTURE OF THE FILE
#
#   1) CALCULATION OF FERMI SURFACE
#
#   2) PLOTTING OF FERMI SURFACE
#       - plotting from points
#       - plotting from unitcell
#
################################################################################

# Function for calculating Fermi surface in 2D (not exported)
function getFermiSurface2D(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(2),
            bounds_upper::Array{Float64,1}=2*pi.*ones(2),
            refold_to_first_BZ::Bool=true
        )

    # list of k values that contribute to the BZ
    k_values = zeros(Float64, N_points, 2)

    # funtion of energy
    function energy(kvector)
        eigenvalues = eigvals(getInteractionMatrixKSpace(unitcell, kvector, enforce_hermitian=enforce_hermitian)) .- fermi_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energy(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:100
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energy(k .+ [epsilon_k, 0]),
                energy(k .+ [0, epsilon_k])
            ]
            dH = (H_eps - H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energy(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = unitcell.lattice_vectors[1]
        a2 = unitcell.lattice_vectors[2]
        # get the reciprocal lattice vectors
        b1 = [a2[2], -a2[1]]
        b2 = [a1[2], -a1[1]]
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)

        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in [b1, b2, -b1, -b2, b1+b2, -b2-b1, b1-b2, -b2+b1]
                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # return the list
    return k_values

end

# Function for calculating Fermi surface in 3D (not exported)
function getFermiSurface3D(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(3),
            bounds_upper::Array{Float64,1}=2*pi.*ones(3),
            refold_to_first_BZ::Bool=true
        )

    # list of k values that contribute to the BZ
    k_values = zeros(Float64, N_points, 3)

    # funtion of energy
    function energy(kvector)
        eigenvalues = eigvals(getInteractionMatrixKSpace(unitcell, kvector, enforce_hermitian=enforce_hermitian)) .- fermi_energy
        eigenvalues = eigenvalues .* eigenvalues
        return minimum(eigenvalues)
    end

    # the current search index
    index = 1
    # search until there are enough points
    while index <= N_points

        # find a suitable starting point for the Newton algorithm
        k = Float64[rand(), rand(), rand()]
        for j in 1:length(k)
            k[j] = k[j] * bounds_lower[j] + (1-k[j]) * bounds_upper[j]
        end

        # start with the initial energy
        e0 = energy(k)
        # iterate i over 100 newton steps (maximum)
        for i in 1:100
            # check if the energy is already converged
            if e0 < epsilon
                # save the k vector
                k_values[index,:] = k
                # increment the index
                index = index+1
                # break the newton loop
                break
            end
            # the current energy
            H_0 = e0
            # the gradient of the energy
            H_eps =
            [
                energy(k .+ [epsilon_k, 0, 0]),
                energy(k .+ [0, epsilon_k, 0]),
                energy(k .+ [0, 0, epsilon_k])
            ]
            dH = (H_eps - H_0) ./ epsilon_k
            # absolute value of the gradient
            dHdH = dot(dH, dH)
            # break the Newton loop if the gradient diverges or is flattened too much
            if abs(dHdH) < 1e-20 || abs(dHdH) > 1e20
                break
            end
            # increment k
            dk = dH * (H_0 / dHdH)
            k -= dk
            # calculate a new energy
            e0 = energy(k)
        end
    end

    # if refolding is enabled, refold all points
    if refold_to_first_BZ
        # get the lattice vectors
        a1 = unitcell.lattice_vectors[1]
        a2 = unitcell.lattice_vectors[2]
        a3 = unitcell.lattice_vectors[3]
        # get the reciprocal lattice vectors
        b1 = cross(a2, a3)
        b2 = cross(a3, a1)
        b3 = cross(a1, a2)
        # normalize the vectors
        b1 .*= 2*pi/sum(b1.*a1)
        b2 .*= 2*pi/sum(b2.*a2)
        b3 .*= 2*pi/sum(b3.*a3)
        # put together refolding
        refoldings = Array{Float64,1}[]
        for i in [-1,0,1]
        for j in [-1,0,1]
        for l in [-1,0,1]
            refolding = i*b1 + j*b2 + l*b3
            if sum(abs.(refolding))>1e-8
                push!(refoldings, refolding)
            end
        end
        end
        end
        # iterate over all k points
        for index in 1:N_points
            # assumend to be refolded to begin with
            refolded = true
            # iterate while still refolding possible
            while refolded
                # new iteration, not refolded
                refolded = false
                # try all possible refoldings
                for r in refoldings
                    if sum((k_values[index,:].-r).*(k_values[index,:].-r)) < sum((k_values[index,:]).*(k_values[index,:]))
                        k_values[index,:] .-= r
                        refolded=true
                    end
                end
            end
        end
    end

    # return the list
    return k_values

end

# obtain Fermi surface
"""
    getFermiSurface(
                unitcell::Unitcell,
                N_points::Int64
             [; fermi_energy::Float64=0.0,
                enforce_hermitian::Bool=false,
                epsilon::Float64=1e-10,
                epsilon_k::Float64=1e-10,
                bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
                bounds_upper::Array{Float64,1}=2*pi.*ones(4),
                refold_to_first_BZ::Bool=true ]
            )

calculates `N_points` points which belong to the Fermi surface of the tight-binding model given by
the `Unitcell` object `unitcell`. The Fermi energy can be adjusted by using `fermi_energy=...`.

The points are located using Newton's method with random starting positions. The procedure can be
modified by using some of the optional keywords:
- `enforce_hermitian` determines if the matrix used for energy calculation is made hermitian by construction
- `epsilon` energy threshold for values to be considered close to the Fermi energy
- `epsilon_k` small k increment that is used for calculation of derivatives
- `bounds_upper` / `bounds_lower` the bounds between which the random starting location for the Newton method is searched

After the search has been completed, the values will optionally be refolded into the first brillouin zone by choosing `refold_to_first_BZ=true`.

# Examples

```julia-repl
julia> fermi_surace = getFermiSurface(unitcell, 100)    # 2D unitcell
100x2 Array{Float64}
...

julia> fermi_surace = getFermiSurface(unitcell, 100)    # 3D unitcell
100x3 Array{Float64}
...

julia> fermi_surace = getFermiSurface(unitcell, 100, fermi_energy=1.0)
100x2 Array{Float64}
...

julia> fermi_surace = getFermiSurface(unitcell, 100, refold_to_first_BZ=false)
100x2 Array{Float64}
...

```
"""
function getFermiSurface(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(3),
            bounds_upper::Array{Float64,1}=2*pi.*ones(3),
            refold_to_first_BZ::Bool=true
        )
    # check if the unitcell can be put in either method
    if length(unitcell.basis[1]) != length(unitcell.lattice_vectors)
        println("Unitcell has not the same number of lattice vectors as dimensions")
        return zeros(1,1)
    end
    # check which function to pass to
    if length(unitcell.lattice_vectors) == 2
        # 2D case
        return getFermiSurface2D(
                unitcell,
                N_points,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )
    else length(unitcell.lattice_vectors) == 3
        # 3D case
        return getFermiSurface3D(
                unitcell,
                N_points,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ
            )
    end

end
export getFermiSurface

# plot the Fermi surface (2D)
function plotFermiSurface2D(
            k_values::Array{Float64, 2};
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT FERMI SURFACE
    ###########################

    # if brillouin zone not empty, plot it as well
    if length(brillouin_zone.points) > 0
        plotBrillouinZone(brillouin_zone, new_figure=false)
    end

    # scatter the points
    scatter(k_values[:,1], k_values[:,2], color=plot_color)

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("Fermi surface")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig

end
function plotFermiSurface2D(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
            bounds_upper::Array{Float64,1}=2*pi.*ones(4),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end

    # calculate the Fermi surface first
    fermi_surface = getFermiSurface2D(
            unitcell,
            N_points,
            fermi_energy=fermi_energy,
            enforce_hermitian=enforce_hermitian,
            epsilon=epsilon,
            epsilon_k=epsilon_k,
            bounds_lower=bounds_lower,
            bounds_upper=bounds_upper,
            refold_to_first_BZ=refold_to_first_BZ
        )

    # plot the fermi surface
    plotFermiSurface2D(
            fermi_surface,
            brillouin_zone=brillouin_zone,
            plot_title=plot_title,
            plot_color=plot_color,
            figsize=figsize,
            showPlot=showPlot,
            save_filename=save_filename
        )

end

# plot the Fermi surface (3D)
function plotFermiSurface3D(
            k_values::Array{Float64, 2};
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT FERMI SURFACE
    ###########################

    # if brillouin zone not empty, plot it as well
    if length(brillouin_zone.points) > 0
        plotBrillouinZone(brillouin_zone, new_figure=false)
    end

    # scatter the points
    scatter3D(k_values[:,1], k_values[:,2], k_values[:,3], color=plot_color)

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("Fermi surface")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig

end
function plotFermiSurface3D(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(3),
            bounds_upper::Array{Float64,1}=2*pi.*ones(3),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end

    # calculate the Fermi surface first
    fermi_surface = getFermiSurface3D(
            unitcell,
            N_points,
            fermi_energy=fermi_energy,
            enforce_hermitian=enforce_hermitian,
            epsilon=epsilon,
            epsilon_k=epsilon_k,
            bounds_lower=bounds_lower,
            bounds_upper=bounds_upper,
            refold_to_first_BZ=refold_to_first_BZ
        )

    # plot the fermi surface
    plotFermiSurface3D(
            fermi_surface,
            brillouin_zone=brillouin_zone,
            plot_title=plot_title,
            plot_color=plot_color,
            figsize=figsize,
            showPlot=showPlot,
            save_filename=save_filename
        )

end

"""
    plotFermiSurface(
                unitcell::Unitcell,
                N_points::Int64
             [; fermi_energy::Float64=0.0,
                enforce_hermitian::Bool=false,
                epsilon::Float64=1e-10,
                epsilon_k::Float64=1e-10,
                bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
                bounds_upper::Array{Float64,1}=2*pi.*ones(4),
                refold_to_first_BZ::Bool=true,
                brillouin_zone::BrillouinZone=BrillouinZone(),
                plot_title::String="",
                plot_color="b",
                figsize::Tuple=(6,6),
                showPlot::Bool=true,
                save_filename::String="NONE" ]
            )

calculates `N_points` points which belong to the Fermi surface of the tight-binding model given by
the `Unitcell` object `unitcell` and plots them using `PyPlot`. The Fermi energy can be adjusted by using `fermi_energy=...`.

The points are located using Newton's method with random starting positions. The procedure can be
modified by using some of the optional keywords:
- `enforce_hermitian` determines if the matrix used for energy calculation is made hermitian by construction
- `epsilon` energy threshold for values to be considered close to the Fermi energy
- `epsilon_k` small k increment that is used for calculation of derivatives
- `bounds_upper` / `bounds_lower` the bounds between which the random starting location for the Newton method is searched

After the search has been completed, the values will optionally be refolded into the first brillouin zone by choosing `refold_to_first_BZ=true`.

For plotting, several more options can be used.

# Examples

```julia-repl
julia> plotFermiSurface(unitcell, 100)
PyPlot.Figure(...)

julia> plotFermiSurface(unitcell, 100, fermi_energy=1.0, showPlot=false)
PyPlot.Figure(...)

```
"""
function plotFermiSurface(
            k_values::Array{Float64, 2};
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    # check which function to pass to
    if size(k_values,2) == 2
        # 2D case
        return plotFermiSurface2D(
                k_values,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
    else size(k_values,2) == 3
        # 3D case
        return plotFermiSurface3D(
                k_values,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
    end

end
function plotFermiSurface(
            unitcell::Unitcell,
            N_points::Int64;
            fermi_energy::Float64=0.0,
            enforce_hermitian::Bool=false,
            epsilon::Float64=1e-10,
            epsilon_k::Float64=1e-10,
            bounds_lower::Array{Float64,1}=-2*pi.*ones(4),
            bounds_upper::Array{Float64,1}=2*pi.*ones(4),
            refold_to_first_BZ::Bool=true,
            brillouin_zone::BrillouinZone=BrillouinZone(),
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,6),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    # check if the unitcell can be put in either method
    if length(unitcell.basis[1]) != length(unitcell.lattice_vectors)
        println("Unitcell has not the same number of lattice vectors as dimensions")
        return zeros(1,1)
    end
    # check if BZ should be calculated
    if length(brillouin_zone.points) == 0
        brillouin_zone = createBrillouinZone(unitcell)
    end
    # check which function to pass to
    if length(unitcell.lattice_vectors) == 2
        # 2D case
        return plotFermiSurface2D(
                unitcell,
                N_points,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
    else length(unitcell.lattice_vectors) == 3
        # 3D case
        return plotFermiSurface3D(
                unitcell,
                N_points,
                fermi_energy=fermi_energy,
                enforce_hermitian=enforce_hermitian,
                epsilon=epsilon,
                epsilon_k=epsilon_k,
                bounds_lower=bounds_lower,
                bounds_upper=bounds_upper,
                refold_to_first_BZ=refold_to_first_BZ,
                brillouin_zone=brillouin_zone,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
    end

end
export plotFermiSurface


text:
using MomentClosure
using MomentClosure: define_M, define_μ
using Symbolics: value, expand
using Test
using Catalyst

@parameters t, k_on, k_off, k_p, γ_p, b
@variables p(t), g(t)

vars = [g, p]
ps = [k_on, k_off, k_p, γ_p, b]
S = [1 -1 0 0;
   0 0 b -1]
as = [k_on*(1-g),  # 0 -> g
    k_off*g*(p^2), # g -> 0
    k_p*g,         # 0 -> mP, m ~ Geometric(mean=b)
    γ_p*p]         # p -> 0
binary_vars = [1]
rn = ReactionSystemMod(t, vars, ps, as, S)

μ = define_μ(2,5)
M = define_M(2,5)

sys = generate_raw_moment_eqs(rn, 3)
expr1 = sys.odes.eqs[4].rhs
expr2 = k_on*μ[0,1] + b*k_p*μ[2,0] - k_off*μ[1,3] - k_on*μ[1,1]-γ_p*μ[1,1]
@test isequal(expr1, expr2)
@test length(sys.odes.eqs) == 9

sys_clean = bernoulli_moment_eqs(sys, binary_vars)
expr1 = sys_clean.odes.eqs[3].rhs
expr2 = k_on*μ[0,1] + b*k_p*μ[1,0] - k_off*μ[1,3] - k_on*μ[1,1]-γ_p*μ[1,1]
expr2 = simplify(value.(expr2))
@test isequal(expr1, expr2)
@test length(sys_clean.odes.eqs) == 6

@test_throws ErrorException moment_closure(sys, "conditional derivative matching")

closed_eqs = moment_closure(sys, "conditional gaussian", binary_vars)
expr1 = closed_eqs.closure[μ[1,4]]
expr2 = 4*μ[1,1]*μ[1,3]*μ[1,0]^-1 + 3*μ[1,2]^2*μ[1,0]^-1 - 12*μ[1,2]*μ[1,1]^2*μ[1,0]^-2 + 6*μ[1,1]^4*μ[1,0]^-3
@test isequal(expr1, expr2)
expr1 = closed_eqs.closure[μ[1,3]]
expr2 = 3*μ[1,2]*μ[1,1]*μ[1,0]^-1 - 2*μ[1,1]^3*μ[1,0]^-2
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars)
expr1 = closed_eqs.closure[μ[1,4]]
expr2 = μ[1,3]^4*μ[1,1]^4*μ[1,0]^-1*μ[1,2]^-6
@test isequal(expr1, expr2)
expr1 = closed_eqs.closure[μ[1,3]]
expr2 = μ[1,2]^3*μ[1,0]*μ[1,1]^-3
@test isequal(expr1, expr2)

sys = generate_central_moment_eqs(rn, 3, 5)
expr1 = sys.odes.eqs[1].rhs
expr2 = -k_off*M[1,2] - k_on*μ[1,0] - k_off*M[0,2]*μ[1,0]-2*k_off*M[1,1]*μ[0,1] - k_off*μ[1,0]*μ[0,1]^2 + k_on
@test isequal(expand(expr1), expr2)
@test length(sys.odes.eqs) == 9

sys_clean = bernoulli_moment_eqs(sys, binary_vars)
expr1 = expand(sys_clean.odes.eqs[3].rhs)
expr2 = b*k_p*μ[1,0] - k_off*M[1,3] - k_on*M[1,1] - M[1,1]*γ_p - b*k_p*μ[1,0]^2 - k_off*M[0,3]*μ[1,0] -
      k_off*M[1,1]*μ[0,1]^2 - 2*k_off*M[1,2]*μ[0,1] - 2*k_off*M[0,2]*μ[0,1]*μ[1,0]
@test isequal(expr1, expr2)
@test length(sys_clean.odes.eqs) == 6

closed_eqs = moment_closure(sys, "conditional gaussian", binary_vars)
expr1 = expand(closed_eqs.closure[M[0,5]])
expr2 = 10*M[0,2]*M[0,3] + 5*M[0,4]*μ[0,1] - 15*μ[0,1]*M[0,2]^2
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars)
expr1 = closed_eqs.closure[M[1,3]]
expr2 = μ[1,0]*(M[1,1]+μ[0,1]*μ[1,0])^-3*(M[1,2]+M[0,2]*μ[1,0]+μ[1,0]*μ[0,1]^2+2*M[1,1]*μ[0,1])^3 -
    M[0,3]*μ[1,0] - 3*M[1,1]*μ[0,1]^2 - 3*M[1,2]*μ[0,1] - μ[1,0]*μ[0,1]^3- 3*M[0,2]*μ[0,1]*μ[1,0]
@test isequal(expand(expr1), expand(expr2))


text:
const TEST = get(ENV, "TEST", "UnitTests")
const CURVE = get(ENV, "RELIC_TOOLKIT_CURVE", "BLS381")

@info TEST, CURVE

eval(Meta.parse("using RelicToolkit.$CURVE"))
include(joinpath(@__DIR__, "$TEST.jl"))


text:
using Photon

λ = 1.55;
Lx = 10 * λ;
sim = Simulation(λ₀ = λ ;  grid = Grid(extent = Lx , spacing =  λ/40 ))
sim(setpml!,λ)
sim(setsrc!, PlaneWave(k̂ = (2pi/λ,), ê = (0,0,1)))
sim(setTFSF!,3λ)
material = Material(ϵᵣ = 7.4)
inBarrier(x) = abs(x[1] - 0.5Lx[1]) < 0.5λ
sim(setmaterial!,material,inBarrier)
sim(solve!)


text:
function forcevec20!(Out, In)

    Out[1+0] = -0.1*In[22] + 0.15958366*^(In[23], 2)*sin(In[2]) + 0.151604477*^(In[24], 2)*sin(In[3]) + 0.143625294*^(In[25], 2)*sin(In[4]) + 0.135646111*^(In[26], 2)*sin(In[5]) + 0.127666928*^(In[27], 2)*sin(In[6]) + 0.119687745*^(In[28], 2)*sin(In[7]) + 0.111708562*^(In[29], 2)*sin(In[8]) + 0.103729379*^(In[30], 2)*sin(In[9]) + 0.095750196*^(In[31], 2)*sin(In[10]) + 0.087771013*^(In[32], 2)*sin(In[11]) + 0.07979183*^(In[33], 2)*sin(In[12]) + 0.071812647*^(In[34], 2)*sin(In[13]) + 0.063833464*^(In[35], 2)*sin(In[14]) + 0.055854281*^(In[36], 2)*sin(In[15]) + 0.047875098*^(In[37], 2)*sin(In[16]) + 0.039895915*^(In[38], 2)*sin(In[17]) + 0.031916732*^(In[39], 2)*sin(In[18]) + 0.023937549*^(In[40], 2)*sin(In[19]) + 0.015958366*^(In[41], 2)*sin(In[20]) + 0.007979183*^(In[42], 2)*sin(In[21]) + In[43];
    Out[1+1] = -0.01*In[23] - 0.151604477*^(In[24], 2)*(0.0589*sin(In[2])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[2])) - 0.143625294*^(In[25], 2)*(0.0589*sin(In[2])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[2])) - 0.135646111*^(In[26], 2)*(0.0589*sin(In[2])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[2])) - 0.127666928*^(In[27], 2)*(0.0589*sin(In[2])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[2])) - 0.119687745*^(In[28], 2)*(0.0589*sin(In[2])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[2])) - 0.111708562*^(In[29], 2)*(0.0589*sin(In[2])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[2])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[2])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[2])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[2]) + 0.0589*sin(In[2])*cos(In[21])) - 1.5655157046*sin(In[2]);
    Out[1+2] = -0.151604477*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[2])) + 0.01*In[23] - 0.01*In[24] - 0.143625294*^(In[25], 2)*(0.0589*sin(In[3])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[3])) - 0.135646111*^(In[26], 2)*(0.0589*sin(In[3])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[3])) - 0.127666928*^(In[27], 2)*(0.0589*sin(In[3])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[3])) - 0.119687745*^(In[28], 2)*(0.0589*sin(In[3])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[3])) - 0.111708562*^(In[29], 2)*(0.0589*sin(In[3])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[3])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[3])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[3])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[3]) + 0.0589*sin(In[3])*cos(In[21])) - 1.48723991937*sin(In[3]);
    Out[1+3] = -0.143625294*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[2])) - 0.143625294*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[3])) + 0.01*In[24] - 0.01*In[25] - 0.135646111*^(In[26], 2)*(0.0589*sin(In[4])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[4])) - 0.127666928*^(In[27], 2)*(0.0589*sin(In[4])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[4])) - 0.119687745*^(In[28], 2)*(0.0589*sin(In[4])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[4])) - 0.111708562*^(In[29], 2)*(0.0589*sin(In[4])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[4])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[4])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[4])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[4]) + 0.0589*sin(In[4])*cos(In[21])) - 1.40896413414*sin(In[4]);
    Out[1+4] = -0.135646111*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[2])) - 0.135646111*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[3])) - 0.135646111*^(In[25], 2)*(-0.0589*sin(In[4])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[4])) + 0.01*In[25] - 0.01*In[26] - 0.127666928*^(In[27], 2)*(0.0589*sin(In[5])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[5])) - 0.119687745*^(In[28], 2)*(0.0589*sin(In[5])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[5])) - 0.111708562*^(In[29], 2)*(0.0589*sin(In[5])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[5])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[5])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[5])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[5]) + 0.0589*sin(In[5])*cos(In[21])) - 1.33068834891*sin(In[5]);
    Out[1+5] = -0.127666928*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[2])) - 0.127666928*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[3])) - 0.127666928*^(In[25], 2)*(-0.0589*sin(In[4])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[4])) - 0.127666928*^(In[26], 2)*(-0.0589*sin(In[5])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[5])) + 0.01*In[26] - 0.01*In[27] - 0.119687745*^(In[28], 2)*(0.0589*sin(In[6])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[6])) - 0.111708562*^(In[29], 2)*(0.0589*sin(In[6])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[6])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[6])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[6])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[6]) + 0.0589*sin(In[6])*cos(In[21])) - 1.25241256368*sin(In[6]);
    Out[1+6] = -0.119687745*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[2])) - 0.119687745*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[3])) - 0.119687745*^(In[25], 2)*(-0.0589*sin(In[4])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[4])) - 0.119687745*^(In[26], 2)*(-0.0589*sin(In[5])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[5])) - 0.119687745*^(In[27], 2)*(-0.0589*sin(In[6])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[6])) + 0.01*In[27] - 0.01*In[28] - 0.111708562*^(In[29], 2)*(0.0589*sin(In[7])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[7])) - 0.103729379*^(In[30], 2)*(0.0589*sin(In[7])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[7])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[7]) + 0.0589*sin(In[7])*cos(In[21])) - 1.17413677845*sin(In[7]);
    Out[1+7] = -0.111708562*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[2])) - 0.111708562*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[3])) - 0.111708562*^(In[25], 2)*(-0.0589*sin(In[4])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[4])) - 0.111708562*^(In[26], 2)*(-0.0589*sin(In[5])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[5])) - 0.111708562*^(In[27], 2)*(-0.0589*sin(In[6])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[6])) - 0.111708562*^(In[28], 2)*(-0.0589*sin(In[7])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[7])) + 0.01*In[28] - 0.01*In[29] - 0.103729379*^(In[30], 2)*(0.0589*sin(In[8])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[8])) - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[8]) + 0.0589*sin(In[8])*cos(In[21])) - 1.09586099322*sin(In[8]);
    Out[1+8] = -0.103729379*^(In[23], 2)*(-0.0589*sin(In[2])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[2])) - 0.103729379*^(In[24], 2)*(-0.0589*sin(In[3])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[3])) - 0.103729379*^(In[25], 2)*(-0.0589*sin(In[4])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[4])) - 0.103729379*^(In[26], 2)*(-0.0589*sin(In[5])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[5])) - 0.103729379*^(In[27], 2)*(-0.0589*sin(In[6])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[6])) - 0.103729379*^(In[28], 2)*(-0.0589*sin(In[7])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[7])) - 0.103729379*^(In[29], 2)*(-0.0589*sin(In[8])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[8])) + 0.01*In[29] - 0.01*In[30] - 0.095750196*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[10])) - 0.087771013*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[11])) - 0.07979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[12])) - 0.071812647*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[13])) - 0.063833464*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[14])) - 0.055854281*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[15])) - 0.047875098*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[16])) - 0.039895915*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[17])) - 0.031916732*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[18])) - 0.023937549*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[19])) - 0.015958366*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[20])) - 0.007979183*^(In[42], 2)*(-0.0589*sin(In[21])*cos(In[9]) + 0.0589*sin(In[9])*cos(In[21])) - 1.01758520799*sin(In[9]);
    Out[1+9] = -0.095750196*^(In[23], 2)*(0.0589*sin(In[10])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[10])) - 0.095750196*^(In[24], 2)*(0.0589*sin(In[10])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[10])) - 0.095750196*^(In[25], 2)*(0.0589*sin(In[10])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[10])) - 0.095750196*^(In[26], 2)*(0.0589*sin(In[10])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[10])) - 0.095750196*^(In[27], 2)*(0.0589*sin(In[10])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[10])) - 0.095750196*^(In[28], 2)*(0.0589*sin(In[10])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[10])) - 0.095750196*^(In[29], 2)*(0.0589*sin(In[10])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[10])) - 0.095750196*^(In[30], 2)*(0.0589*sin(In[10])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[10])) + 0.01*In[30] - 0.01*In[31] - 0.087771013*^(In[32], 2)*(0.0589*sin(In[10])*cos(In[11]) - 0.0589*sin(In[11])*cos(In[10])) - 0.07979183*^(In[33], 2)*(0.0589*sin(In[10])*cos(In[12]) - 0.0589*sin(In[12])*cos(In[10])) - 0.071812647*^(In[34], 2)*(0.0589*sin(In[10])*cos(In[13]) - 0.0589*sin(In[13])*cos(In[10])) - 0.063833464*^(In[35], 2)*(0.0589*sin(In[10])*cos(In[14]) - 0.0589*sin(In[14])*cos(In[10])) - 0.055854281*^(In[36], 2)*(0.0589*sin(In[10])*cos(In[15]) - 0.0589*sin(In[15])*cos(In[10])) - 0.047875098*^(In[37], 2)*(0.0589*sin(In[10])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[10])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[10])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[10])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[10])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[10])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[10])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[10])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[10])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[10])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[10])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[10])) - 0.93930942276*sin(In[10]);
    Out[1+10] = -0.087771013*^(In[23], 2)*(0.0589*sin(In[11])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[11])) - 0.087771013*^(In[24], 2)*(0.0589*sin(In[11])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[11])) - 0.087771013*^(In[25], 2)*(0.0589*sin(In[11])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[11])) - 0.087771013*^(In[26], 2)*(0.0589*sin(In[11])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[11])) - 0.087771013*^(In[27], 2)*(0.0589*sin(In[11])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[11])) - 0.087771013*^(In[28], 2)*(0.0589*sin(In[11])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[11])) - 0.087771013*^(In[29], 2)*(0.0589*sin(In[11])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[11])) - 0.087771013*^(In[30], 2)*(0.0589*sin(In[11])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[11])) - 0.087771013*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[11]) + 0.0589*sin(In[11])*cos(In[10])) + 0.01*In[31] - 0.01*In[32] - 0.07979183*^(In[33], 2)*(0.0589*sin(In[11])*cos(In[12]) - 0.0589*sin(In[12])*cos(In[11])) - 0.071812647*^(In[34], 2)*(0.0589*sin(In[11])*cos(In[13]) - 0.0589*sin(In[13])*cos(In[11])) - 0.063833464*^(In[35], 2)*(0.0589*sin(In[11])*cos(In[14]) - 0.0589*sin(In[14])*cos(In[11])) - 0.055854281*^(In[36], 2)*(0.0589*sin(In[11])*cos(In[15]) - 0.0589*sin(In[15])*cos(In[11])) - 0.047875098*^(In[37], 2)*(0.0589*sin(In[11])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[11])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[11])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[11])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[11])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[11])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[11])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[11])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[11])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[11])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[11])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[11])) - 0.86103363753*sin(In[11]);
    Out[1+11] = -0.07979183*^(In[23], 2)*(0.0589*sin(In[12])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[12])) - 0.07979183*^(In[24], 2)*(0.0589*sin(In[12])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[12])) - 0.07979183*^(In[25], 2)*(0.0589*sin(In[12])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[12])) - 0.07979183*^(In[26], 2)*(0.0589*sin(In[12])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[12])) - 0.07979183*^(In[27], 2)*(0.0589*sin(In[12])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[12])) - 0.07979183*^(In[28], 2)*(0.0589*sin(In[12])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[12])) - 0.07979183*^(In[29], 2)*(0.0589*sin(In[12])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[12])) - 0.07979183*^(In[30], 2)*(0.0589*sin(In[12])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[12])) - 0.07979183*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[12]) + 0.0589*sin(In[12])*cos(In[10])) - 0.07979183*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[12]) + 0.0589*sin(In[12])*cos(In[11])) + 0.01*In[32] - 0.01*In[33] - 0.071812647*^(In[34], 2)*(0.0589*sin(In[12])*cos(In[13]) - 0.0589*sin(In[13])*cos(In[12])) - 0.063833464*^(In[35], 2)*(0.0589*sin(In[12])*cos(In[14]) - 0.0589*sin(In[14])*cos(In[12])) - 0.055854281*^(In[36], 2)*(0.0589*sin(In[12])*cos(In[15]) - 0.0589*sin(In[15])*cos(In[12])) - 0.047875098*^(In[37], 2)*(0.0589*sin(In[12])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[12])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[12])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[12])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[12])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[12])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[12])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[12])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[12])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[12])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[12])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[12])) - 0.7827578523*sin(In[12]);
    Out[1+12] = -0.071812647*^(In[23], 2)*(0.0589*sin(In[13])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[13])) - 0.071812647*^(In[24], 2)*(0.0589*sin(In[13])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[13])) - 0.071812647*^(In[25], 2)*(0.0589*sin(In[13])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[13])) - 0.071812647*^(In[26], 2)*(0.0589*sin(In[13])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[13])) - 0.071812647*^(In[27], 2)*(0.0589*sin(In[13])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[13])) - 0.071812647*^(In[28], 2)*(0.0589*sin(In[13])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[13])) - 0.071812647*^(In[29], 2)*(0.0589*sin(In[13])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[13])) - 0.071812647*^(In[30], 2)*(0.0589*sin(In[13])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[13])) - 0.071812647*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[13]) + 0.0589*sin(In[13])*cos(In[10])) - 0.071812647*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[13]) + 0.0589*sin(In[13])*cos(In[11])) - 0.071812647*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[13]) + 0.0589*sin(In[13])*cos(In[12])) + 0.01*In[33] - 0.01*In[34] - 0.063833464*^(In[35], 2)*(0.0589*sin(In[13])*cos(In[14]) - 0.0589*sin(In[14])*cos(In[13])) - 0.055854281*^(In[36], 2)*(0.0589*sin(In[13])*cos(In[15]) - 0.0589*sin(In[15])*cos(In[13])) - 0.047875098*^(In[37], 2)*(0.0589*sin(In[13])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[13])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[13])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[13])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[13])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[13])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[13])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[13])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[13])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[13])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[13])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[13])) - 0.70448206707*sin(In[13]);
    Out[1+13] = -0.063833464*^(In[23], 2)*(0.0589*sin(In[14])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[14])) - 0.063833464*^(In[24], 2)*(0.0589*sin(In[14])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[14])) - 0.063833464*^(In[25], 2)*(0.0589*sin(In[14])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[14])) - 0.063833464*^(In[26], 2)*(0.0589*sin(In[14])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[14])) - 0.063833464*^(In[27], 2)*(0.0589*sin(In[14])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[14])) - 0.063833464*^(In[28], 2)*(0.0589*sin(In[14])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[14])) - 0.063833464*^(In[29], 2)*(0.0589*sin(In[14])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[14])) - 0.063833464*^(In[30], 2)*(0.0589*sin(In[14])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[14])) - 0.063833464*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[14]) + 0.0589*sin(In[14])*cos(In[10])) - 0.063833464*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[14]) + 0.0589*sin(In[14])*cos(In[11])) - 0.063833464*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[14]) + 0.0589*sin(In[14])*cos(In[12])) - 0.063833464*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[14]) + 0.0589*sin(In[14])*cos(In[13])) + 0.01*In[34] - 0.01*In[35] - 0.055854281*^(In[36], 2)*(0.0589*sin(In[14])*cos(In[15]) - 0.0589*sin(In[15])*cos(In[14])) - 0.047875098*^(In[37], 2)*(0.0589*sin(In[14])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[14])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[14])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[14])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[14])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[14])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[14])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[14])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[14])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[14])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[14])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[14])) - 0.62620628184*sin(In[14]);
    Out[1+14] = -0.055854281*^(In[23], 2)*(0.0589*sin(In[15])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[15])) - 0.055854281*^(In[24], 2)*(0.0589*sin(In[15])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[15])) - 0.055854281*^(In[25], 2)*(0.0589*sin(In[15])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[15])) - 0.055854281*^(In[26], 2)*(0.0589*sin(In[15])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[15])) - 0.055854281*^(In[27], 2)*(0.0589*sin(In[15])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[15])) - 0.055854281*^(In[28], 2)*(0.0589*sin(In[15])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[15])) - 0.055854281*^(In[29], 2)*(0.0589*sin(In[15])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[15])) - 0.055854281*^(In[30], 2)*(0.0589*sin(In[15])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[15])) - 0.055854281*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[15]) + 0.0589*sin(In[15])*cos(In[10])) - 0.055854281*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[15]) + 0.0589*sin(In[15])*cos(In[11])) - 0.055854281*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[15]) + 0.0589*sin(In[15])*cos(In[12])) - 0.055854281*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[15]) + 0.0589*sin(In[15])*cos(In[13])) - 0.055854281*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[15]) + 0.0589*sin(In[15])*cos(In[14])) + 0.01*In[35] - 0.01*In[36] - 0.047875098*^(In[37], 2)*(0.0589*sin(In[15])*cos(In[16]) - 0.0589*sin(In[16])*cos(In[15])) - 0.039895915*^(In[38], 2)*(0.0589*sin(In[15])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[15])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[15])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[15])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[15])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[15])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[15])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[15])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[15])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[15])) - 0.54793049661*sin(In[15]);
    Out[1+15] = -0.047875098*^(In[23], 2)*(0.0589*sin(In[16])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[16])) - 0.047875098*^(In[24], 2)*(0.0589*sin(In[16])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[16])) - 0.047875098*^(In[25], 2)*(0.0589*sin(In[16])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[16])) - 0.047875098*^(In[26], 2)*(0.0589*sin(In[16])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[16])) - 0.047875098*^(In[27], 2)*(0.0589*sin(In[16])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[16])) - 0.047875098*^(In[28], 2)*(0.0589*sin(In[16])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[16])) - 0.047875098*^(In[29], 2)*(0.0589*sin(In[16])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[16])) - 0.047875098*^(In[30], 2)*(0.0589*sin(In[16])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[16])) - 0.047875098*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[10])) - 0.047875098*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[11])) - 0.047875098*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[12])) - 0.047875098*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[13])) - 0.047875098*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[14])) - 0.047875098*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[16]) + 0.0589*sin(In[16])*cos(In[15])) + 0.01*In[36] - 0.01*In[37] - 0.039895915*^(In[38], 2)*(0.0589*sin(In[16])*cos(In[17]) - 0.0589*sin(In[17])*cos(In[16])) - 0.031916732*^(In[39], 2)*(0.0589*sin(In[16])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[16])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[16])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[16])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[16])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[16])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[16])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[16])) - 0.46965471138*sin(In[16]);
    Out[1+16] = -0.039895915*^(In[23], 2)*(0.0589*sin(In[17])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[17])) - 0.039895915*^(In[24], 2)*(0.0589*sin(In[17])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[17])) - 0.039895915*^(In[25], 2)*(0.0589*sin(In[17])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[17])) - 0.039895915*^(In[26], 2)*(0.0589*sin(In[17])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[17])) - 0.039895915*^(In[27], 2)*(0.0589*sin(In[17])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[17])) - 0.039895915*^(In[28], 2)*(0.0589*sin(In[17])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[17])) - 0.039895915*^(In[29], 2)*(0.0589*sin(In[17])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[17])) - 0.039895915*^(In[30], 2)*(0.0589*sin(In[17])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[17])) - 0.039895915*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[10])) - 0.039895915*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[11])) - 0.039895915*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[12])) - 0.039895915*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[13])) - 0.039895915*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[14])) - 0.039895915*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[15])) - 0.039895915*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[17]) + 0.0589*sin(In[17])*cos(In[16])) + 0.01*In[37] - 0.01*In[38] - 0.031916732*^(In[39], 2)*(0.0589*sin(In[17])*cos(In[18]) - 0.0589*sin(In[18])*cos(In[17])) - 0.023937549*^(In[40], 2)*(0.0589*sin(In[17])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[17])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[17])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[17])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[17])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[17])) - 0.39137892615*sin(In[17]);
    Out[1+17] = -0.031916732*^(In[23], 2)*(0.0589*sin(In[18])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[18])) - 0.031916732*^(In[24], 2)*(0.0589*sin(In[18])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[18])) - 0.031916732*^(In[25], 2)*(0.0589*sin(In[18])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[18])) - 0.031916732*^(In[26], 2)*(0.0589*sin(In[18])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[18])) - 0.031916732*^(In[27], 2)*(0.0589*sin(In[18])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[18])) - 0.031916732*^(In[28], 2)*(0.0589*sin(In[18])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[18])) - 0.031916732*^(In[29], 2)*(0.0589*sin(In[18])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[18])) - 0.031916732*^(In[30], 2)*(0.0589*sin(In[18])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[18])) - 0.031916732*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[10])) - 0.031916732*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[11])) - 0.031916732*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[12])) - 0.031916732*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[13])) - 0.031916732*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[14])) - 0.031916732*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[15])) - 0.031916732*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[16])) - 0.031916732*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[18]) + 0.0589*sin(In[18])*cos(In[17])) + 0.01*In[38] - 0.01*In[39] - 0.023937549*^(In[40], 2)*(0.0589*sin(In[18])*cos(In[19]) - 0.0589*sin(In[19])*cos(In[18])) - 0.015958366*^(In[41], 2)*(0.0589*sin(In[18])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[18])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[18])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[18])) - 0.31310314092*sin(In[18]);
    Out[1+18] = -0.023937549*^(In[23], 2)*(0.0589*sin(In[19])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[19])) - 0.023937549*^(In[24], 2)*(0.0589*sin(In[19])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[19])) - 0.023937549*^(In[25], 2)*(0.0589*sin(In[19])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[19])) - 0.023937549*^(In[26], 2)*(0.0589*sin(In[19])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[19])) - 0.023937549*^(In[27], 2)*(0.0589*sin(In[19])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[19])) - 0.023937549*^(In[28], 2)*(0.0589*sin(In[19])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[19])) - 0.023937549*^(In[29], 2)*(0.0589*sin(In[19])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[19])) - 0.023937549*^(In[30], 2)*(0.0589*sin(In[19])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[19])) - 0.023937549*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[10])) - 0.023937549*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[11])) - 0.023937549*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[12])) - 0.023937549*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[13])) - 0.023937549*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[14])) - 0.023937549*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[15])) - 0.023937549*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[16])) - 0.023937549*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[17])) - 0.023937549*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[19]) + 0.0589*sin(In[19])*cos(In[18])) + 0.01*In[39] - 0.01*In[40] - 0.015958366*^(In[41], 2)*(0.0589*sin(In[19])*cos(In[20]) - 0.0589*sin(In[20])*cos(In[19])) - 0.007979183*^(In[42], 2)*(0.0589*sin(In[19])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[19])) - 0.23482735569*sin(In[19]);
    Out[1+19] = -0.015958366*^(In[23], 2)*(0.0589*sin(In[20])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[20])) - 0.015958366*^(In[24], 2)*(0.0589*sin(In[20])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[20])) - 0.015958366*^(In[25], 2)*(0.0589*sin(In[20])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[20])) - 0.015958366*^(In[26], 2)*(0.0589*sin(In[20])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[20])) - 0.015958366*^(In[27], 2)*(0.0589*sin(In[20])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[20])) - 0.015958366*^(In[28], 2)*(0.0589*sin(In[20])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[20])) - 0.015958366*^(In[29], 2)*(0.0589*sin(In[20])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[20])) - 0.015958366*^(In[30], 2)*(0.0589*sin(In[20])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[20])) - 0.015958366*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[10])) - 0.015958366*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[11])) - 0.015958366*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[12])) - 0.015958366*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[13])) - 0.015958366*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[14])) - 0.015958366*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[15])) - 0.015958366*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[16])) - 0.015958366*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[17])) - 0.015958366*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[18])) - 0.015958366*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[20]) + 0.0589*sin(In[20])*cos(In[19])) + 0.01*In[40] - 0.01*In[41] - 0.007979183*^(In[42], 2)*(0.0589*sin(In[20])*cos(In[21]) - 0.0589*sin(In[21])*cos(In[20])) - 0.15655157046*sin(In[20]);
    Out[1+20] = -0.007979183*^(In[23], 2)*(0.0589*sin(In[21])*cos(In[2]) - 0.0589*sin(In[2])*cos(In[21])) - 0.007979183*^(In[24], 2)*(0.0589*sin(In[21])*cos(In[3]) - 0.0589*sin(In[3])*cos(In[21])) - 0.007979183*^(In[25], 2)*(0.0589*sin(In[21])*cos(In[4]) - 0.0589*sin(In[4])*cos(In[21])) - 0.007979183*^(In[26], 2)*(0.0589*sin(In[21])*cos(In[5]) - 0.0589*sin(In[5])*cos(In[21])) - 0.007979183*^(In[27], 2)*(0.0589*sin(In[21])*cos(In[6]) - 0.0589*sin(In[6])*cos(In[21])) - 0.007979183*^(In[28], 2)*(0.0589*sin(In[21])*cos(In[7]) - 0.0589*sin(In[7])*cos(In[21])) - 0.007979183*^(In[29], 2)*(0.0589*sin(In[21])*cos(In[8]) - 0.0589*sin(In[8])*cos(In[21])) - 0.007979183*^(In[30], 2)*(0.0589*sin(In[21])*cos(In[9]) - 0.0589*sin(In[9])*cos(In[21])) - 0.007979183*^(In[31], 2)*(-0.0589*sin(In[10])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[10])) - 0.007979183*^(In[32], 2)*(-0.0589*sin(In[11])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[11])) - 0.007979183*^(In[33], 2)*(-0.0589*sin(In[12])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[12])) - 0.007979183*^(In[34], 2)*(-0.0589*sin(In[13])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[13])) - 0.007979183*^(In[35], 2)*(-0.0589*sin(In[14])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[14])) - 0.007979183*^(In[36], 2)*(-0.0589*sin(In[15])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[15])) - 0.007979183*^(In[37], 2)*(-0.0589*sin(In[16])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[16])) - 0.007979183*^(In[38], 2)*(-0.0589*sin(In[17])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[17])) - 0.007979183*^(In[39], 2)*(-0.0589*sin(In[18])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[18])) - 0.007979183*^(In[40], 2)*(-0.0589*sin(In[19])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[19])) - 0.007979183*^(In[41], 2)*(-0.0589*sin(In[20])*cos(In[21]) + 0.0589*sin(In[21])*cos(In[20])) + 0.01*In[41] - 0.01*In[42] - 0.07827578523*sin(In[21]);

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
@inline blockIdx_x() = (ccall("llvm.nvvm.read.ptx.sreg.ctaid.x", llvmcall, UInt32, ()))+UInt32(1)

function kernel_vadd()
    blockIdx_x()-UInt32(1)
end

code_llvm(stdout, kernel_vadd, Tuple{})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx,
                typeof(kernel_vadd)})