-- Defining single-valued response functions

a = {0, -1,-0.5, 0,   0.5,  1,  1.5}
b = {  0.2, 0.1, 0.1, 0.1, 0.2, 0.3}
GA_L = {a,b,mu=0,type="ListOfPoles", name="A"}
setmetatable(GA_L, ResponseFunctionMeta)

a = {0, 1,   1,   1,   1,   1,  1}
b = {   1, 0.5, 0.5, 0.5, 0.5, 0.5}
GB_T = {a,b,mu=0,type="Tri", name="B"}
setmetatable(GB_T, ResponseFunctionMeta)

a = {0, 1, 1.5,   2, 2.5,   3, 3.5}
b = {   1, 0.5, 0.5, 0.5, 0.5, 0.5}
GC_A = {a,b,mu=0,type="And", name="C"}
setmetatable(GC_A, ResponseFunctionMeta)

acon = {0,         1,   1,   1,   1,   1,  1}
bcon = {   sqrt(1/2), 0.5, 0.5, 0.5, 0.5, 0.5}
Gcon = {acon,bcon,mu=0,type="Tri"}
setmetatable(Gcon, ResponseFunctionMeta)
aval = {0,        -1,  -1,  -1,  -1,  -1, -1}
bval = {   sqrt(1/2), 0.5, 0.5, 0.5, 0.5, 0.5}
Gval = {aval,bval,mu=0,type="Tri"}
setmetatable(Gval, ResponseFunctionMeta)
a0=0
b0=1
GD_N = {{a0,b0},val=Gval,con=Gcon,mu=0,type="Nat", name="D"}
setmetatable(GD_N, ResponseFunctionMeta)

-- For a more efficient storage, one can convert from Table to Userdata
GA_l = ResponseFunction.ToUserdata(GA_L)
GB_t = ResponseFunction.ToUserdata(GB_T)
GC_a = ResponseFunction.ToUserdata(GC_A)
GD_n = ResponseFunction.ToUserdata(GD_N)

print("GA_L:")
print(GA_L)
print("GA_l:")
print(GA_l)

print("GB_T:")
print(GB_T)
print("GB_t:")
print(GB_t)

print("GC_A:")
print(GC_A)
print("GC_a:")
print(GC_a)

print("GD_N:")
print(GD_N)
print("GD_n:")
print(GD_n)

-- response functions in Quanty are functions that, given 
-- a complex number as input w + I gamma/2 return a complex
-- number (single valued functions) or a matrix (matrix functions)
omega = 0.764
gamma = 0.1
print("omega = ", omega)
print("gamma = ", gamma)

print("")
print("GA_L(omega, gamma) = ", GA_L(omega, gamma))
print("GA_l(omega, gamma) = ", GA_l(omega, gamma))

print("GB_T(omega, gamma) = ", GB_T(omega, gamma))
print("GB_t(omega, gamma) = ", GB_t(omega, gamma))

print("GC_A(omega, gamma) = ", GC_A(omega, gamma))
print("GC_a(omega, gamma) = ", GC_a(omega, gamma))

print("GD_N(omega, gamma) = ", GD_N(omega, gamma))
print("GD_n(omega, gamma) = ", GD_n(omega, gamma))


-- Defining matrix response functions

A0 = {{0,0,0},{0,0,0},{0,0,0}}
setmetatable(A0, MatrixMeta)
A1 = {{1,2,3},{2,5,6},{3,6,9}}
setmetatable(A1, MatrixMeta)
A2 = {{2,2,3},{2,5,6},{3,6,9}}
setmetatable(A2, MatrixMeta)
A3 = {{3,2,3},{2,5,6},{3,6,9}}
setmetatable(A3, MatrixMeta)
a1 = -1
a2 = 0
a3 = 1
av1 = -0.5
av2 = -1.0
av3 = -1.5
ac1 = 0.5
ac2 = 1.0
ac3 = 1.5
B0s = {{1,0,0},{0,1,0},{0,0,1}}
setmetatable(B0s, MatrixMeta)
t=sqrt(0.5)
B0vs = {{t,0,0},{0,t,0},{0,0,t}}
setmetatable(B0vs, MatrixMeta)
B0cs = {{t,0,0},{0,t,0},{0,0,t}}
setmetatable(B0cs, MatrixMeta)
B0 = B0s * B0s
B0v = B0vs * B0vs
B0c = B0cs * B0cs
B1s = {{1,I,3},{-I,5,6},{3,6,9}}
setmetatable(B1s, MatrixMeta)
B1 = B1s * B1s
B2s = {{2,0,3},{0,5,6},{3,6,9}}
setmetatable(B2s, MatrixMeta)
B2 = B2s * B2s
B3s = {{3,0,3},{0,5,6},{3,6,9}}
setmetatable(B3s, MatrixMeta)
B3 = B3s * B3s

MA_L = { {A0,a1,a2,a3}, {B1,B2,B3}, mu=0, type="ListOfPoles", name="A"}
setmetatable(MA_L, ResponseFunctionMeta)

MB_T = { {A0,A1,A2,A3}, {B0,B1,B2}, mu=0, type="Tri", name="B"}
setmetatable(MB_T, ResponseFunctionMeta)

MC_A = { {A0,A1,A2,A3}, {B0,B1,B2}, mu=0, type="And", name="C"}
setmetatable(MC_A, ResponseFunctionMeta)

MD_Nv = { {A0,av1,av2,av3}, {B1,B2,B3}, mu=0, type="ListOfPoles", name="Dv"}
setmetatable(MD_Nv, ResponseFunctionMeta)

MD_Nc = { {A0,ac1,ac2,ac3}, {B1,B2,B3}, mu=0, type="ListOfPoles", name="Dc"}
setmetatable(MD_Nc, ResponseFunctionMeta)

MD_N = {{A0,B0},val=MD_Nv,con=MD_Nc,mu=0,type="Nat", name="D"}
setmetatable(MD_N, ResponseFunctionMeta)

-- For a more efficient storage, one can convert from Table to Userdata
MA_l = ResponseFunction.ToUserdata(MA_L)
MB_t = ResponseFunction.ToUserdata(MB_T)
MC_a = ResponseFunction.ToUserdata(MC_A)
MD_n = ResponseFunction.ToUserdata(MD_N)


print("MA_L:")
print(MA_L)
print("MA_l:")
print(MA_l)

print("MB_T:")
print(MB_T)
print("MB_t:")
print(MB_t)

print("MC_A:")
print(MC_A)
print("MC_a:")
print(MC_a)

print("MD_N:")
print(MD_N)
print("MD_n:")
print(MD_n)


-- response functions in Quanty are functions that, given 
-- a complex number as input w + I gamma/2 return a complex
-- number (single valued functions) or a matrix (matrix functions)
omega = 0.764
gamma = 0.1
print("omega = ", omega)
print("gamma = ", gamma)

print("")
print("MA_L(omega, gamma) = ")
print(MA_L(omega, gamma))
print("MA_l(omega, gamma) = ")
print(MA_l(omega, gamma))

print("")
print("MB_T(omega, gamma) = ")
print(MB_T(omega, gamma))
print("MB_t(omega, gamma) = ")
print(MB_t(omega, gamma))

print("")
print("MC_A(omega, gamma) = ")
print(MC_A(omega, gamma))
print("MC_a(omega, gamma) = ")
print(MC_a(omega, gamma))

print("")
print("MD_N(omega, gamma) = ")
print(MD_N(omega, gamma))
print("MD_n(omega, gamma) = ")
print(MD_n(omega, gamma))