import ADerrors
const beta_to_idx = Dict{Float64,Int64}(3.40 =>1, 3.46=>2, 3.55=>3, 3.70=>4, 3.85=>5)

let 
  C = zeros(8, 8)
  M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916, 5279.65, 5366.3] #MD, MD*, MDs, MDs*, \eta_c, J/\psi , MB, MBs(MeV)
  M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011, 0.12, 0.6]
  M_ = Ref{Vector{ADerrors.uwreal}}()
  for i = 1:8
    C[i, i] = M_error[i] ^ 2
  end
  global function M()
    if !isassigned(M_) 
      M_.x = ADerrors.cobs(M_values, C, "charmed meson masses")
    end
    return M_.x;
  end
end

let 
  #1802.05243 taking into account correlations in zm_tm
  C = [0.375369e-6 0.429197e-6 -0.186896e-5; 0.429197e-6 0.268393e-4 0.686776e-4; -0.186896e-5 0.686776e-4 0.212386e-3]
  z = Ref{Vector{ADerrors.uwreal}}();
  global function ZP(beta::Float64)
    if !isassigned(z)
      z.x =  ADerrors.cobs([0.348629, 0.020921, 0.070613], C, "z")
    end
    return z.x[1] + z.x[2] *(beta-3.79) + z.x[3]*(beta-3.79)^2
  end
  CC = [0.164635e-4 0.215658e-4 -0.754203e-4; 0.215658e-4 0.121072e-2 0.308890e-2; -0.754203e-4 0.308890e-2 0.953843e-2]
  zm = Ref{Vector{ADerrors.uwreal}}()
  Mrat_ = Ref{ADerrors.uwreal}()

  global function Mrat()
    if !isassigned(Mrat_)
      Mrat_.x =ADerrors.uwreal([.9148, 0.0088], "M/mhad")
    end
    return Mrat_.x
  end

  global function z_covar(beta::Float64)
    if !isassigned(zz)
      zm.x = ADerrors.cobs([2.270073, 0.121644, -0.464575], CC, "zm")
    end
    return Mrat().mean*(zz.x[1] + zz.x[2]*(beta-3.79)+zz.x[3]*(b-3.79)^2)
  end
  
  global  zm_tm_covar(beta::Float64) = Mrat() / ZP(beta)
end

## taking into account correlations in r_m
#Crm = [[3.699760, -2.198586, -1.476913e-3] [-2.198586, 1.306512, 8.776569e-4] [-1.476913e-3, 8.776569e-4, 5.895710e-7] ]
#c_i = ADerrors.cobs([-7.86078, 5.49175, -0.54078], Crm, "c_i in rm")
#rm(b::Float64) = 1.0 +  0.004630*(6/b)^4 * ((1. +c_i[1]*(6/b)^2 + c_i[2]*(6/b)^4) / (1. + c_i[3]*(6/b)^2))

let C = zeros(5, 5)
  #1802.05243 
  ZM_error = [35, 27, 33, 38, 45] .* 1e-4
  ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]

  for i in 1:5
    C[i, i] = ZM_error[i] ^ 2
  end
  ZM = Ref{Vector{ADerrors.uwreal}}()
  global function zm(beta::Float64)
    if !isassigned(ZM)
      ZM.x =  ADerrors.cobs(ZM_data, C, "ZM values")
    end
    return ZM[beta_to_idx[beta]]
  end
end

let C = zeros(5, 5)
  #1802.05243 
  ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
  ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
  for i in 1:5
    C[i, i] = ZM_tm_error[i] ^ 2
  end
  ZM_tm = Ref{Vector{ADerrors.uwreal}}()
  
  global function zm_tm(beta::Float64)
    if !isassigned(ZM_tm)
      ZM_tm.x = ADerrors.cobs(ZM_tm_data, C, "ZM_tm values")
    end
    return ZM_tm.x[beta_to_idx[beta]]
  end
end

let C = zeros(5, 5)
  #2101.02694
  t0_data = [2.846, 3.634, 5.140, 8.564, 14.064]
  t0_error = [8, 13, 21, 24, 63] .* 1e-3

  for i in 1:5
    C[i, i] = t0_error[i] ^ 2
  end
  
  t0_ = Ref{Vector{ADerrors.uwreal}}()
  global function t0(beta::Float64)
    if !isassigned(beta)
      t0_.x = ADerrors.cobs(t0_data, C, "t0")
    end
    return t0_.x[beta_to_idx[beta]];
  end

  t0_ph_ = Ref{ADerrors.uwreal}()
  global function t0_ph()
    if !isassigned(t0_ph)
      t0_ph_.x = ADerrors.uwreal([0.4118,0.25e-4],"sqrt(8 t0) (fm)") 
    end
    return t0_ph_.x
  end

  global function a(beta::Float64)
    return t0_ph() / sqrt(8*t0_.x[beta_to_idx[beta]]);
  end
end


let C = zeros(5, 5)
  #1808.09236
  ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
  ZA_err = [72, 93, 43, 47, 47] .*1e-5

  for i in 1:5
    C[i, i] = ZA_err[i] ^ 2
  end
  Za = Ref{Vector{ADerrors.uwreal}}()

  global function za(beta::Float64)
    if !isassigned(Za)
      Za.x = ADerrors.cobs(ZA_data, C, "ZA")
    end
    return za.x[beta_to_idx[beta]];
  end
end

let C = zeros(5, 5)
  #1808.09236
  ZV_data = [0.72259, 0.72845, 0.73886, 0.75574, 0.77074]
  ZV_err = [74, 89, 46, 48, 49] .* 1e-5
  
  for i in 1:5
    C[i, i] = ZV_err[i] ^ 2 
  end
  Zv = Ref{Vector{ADerrors.uwreal}}()

  global function zv(beta::Float64)
    if !isassigned(Zv)
      Zv.x = ADerrors.cobs(ZV_data, C, "ZV")
    end
    return zv.x[beta_to_idx[beta]];
  end

end

let C = zeros(5, 5)
  #2005.01352
  ZS_over_ZP_data = [1.3497, 1.2914, 1.2317, 1.1709, 1.1343]
  ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4

  C = zeros(5,5)
  for i in 1:5
    C[i, i] = ZS_over_ZP_err[i] ^ 2
  end
  ZS_over_ZP = Ref{Vector{ADerrors.uwreal}}()
  global function zs_over_zp(beta::Float64)
    if !isassigned(ZS_over_ZP_data)
      ZS_over_ZP.x = ADerrors.cobs(ZS_over_ZP_data, C, "ZS/ZP")
    end
    return ZS_over_ZP.x[beta_to_idx[beta]];
  end 
end

let C = zeros(5, 5)
  #2101.10969
  RM_err  = [31, 19, 14, 16, 18] .* 1e-3
  RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
  C=zeros(5,5)
  for i in 1:5
    C[i, i] = RM_err[i] ^2 
  end
  RM  = Ref{Vector{ADerrors.uwreal}}()
  global function rm(beta::Float64)
    if !isassigned(rm)
      RM.x = ADerrors.cobs(RM_data, C, "rm")
    end
    return rm.x[beta_to_idx[beta]];
  end  
end


let C = zeros(5, 5)
  #1906.03445 LCP1 for heavy quarks
  BA_MINUS_BP_data = [-0.324, -0.265, -0.196, -0.119, -0.073]
  BA_MINUS_BP_err  = [17, 14, 14, 14, 12] .*1e-3 
  C=zeros(5,5)
  for i in 1:5
    C[i,i]  = BA_MINUS_BP_err[i] ^2
  end
  BA_BP = Ref{Vector{ADerrors.uwreal}}()
  global function ba_bp(beta::Float64)
    if !isassigned(BA_BP)
      BA_BP.x =  ADerrors.cobs(BA_MINUS_BP_data, C, "ba-bp")
    end
    return  BA_BP.x[beta_to_idx[beta]];
  end
end

let C = zeros(5, 5)
  #1906.03445 LCP1 for heavy quarks 
  ZDATA = [0.8981, 0.9468, 1.0015, 1.0612, 1.0971]
  ZERR  = [35, 35, 30, 17, 18 ] .* 1e-4
  for i in 1:5
    C[i,i]  = ZERR[i] ^2
  end
  ZZ = Ref{Vector{ADerrors.uwreal}}();
  global function Z(beta::Float64)
    if !isassinged(ZZ)
      ZZ.x =ADerrors.cobs(ZDATA, C, "Z")
    end
    return ZZ[beta_to_idx[beta]];
  end
end