using Plots
using Plots.PlotMeasures
using LaTeXStrings

B(x) = ((2J+1)/2J)*coth(((2J+1)/2J)*x) - (1/2J)*coth(x/2J)
J = 7/2;
x = collect(0:0.1:30);
y = B.(x);

plot(x, y)
plot!(size=(500, 300), xlabel=L"x", ylabel=L"B(x)", legend = false)

μB = 9.27410*10^(-21); #ボーア磁子
kB = 1.3866*10^(-16); #ボルツマン定数
NA = 6.02217*10^(23); #アボガドロ数
J = 7/2; #全角運動量
L = 0;   #軌道角運動量
S = 7/2; #スピン
T = 2;   #温度
g = 3/2 + (S*(S+1)-L*(L+1))/(2J*(J+1)) #g因子
B(x) = ((2J+1)/2J)*coth(((2J+1)/2J)*x) - (1/2J)*coth(x/2J) #ブリュアン関数
MH(H) = g*J*B((J*g*μB*H)/(kB*T)) #磁化の式(磁場依存性として設定)

h = collect(1:100:70000); #磁場（Oe）
m = MH.(h) #磁化（μB/ion）

plot(h,m)
plot!(size = (500, 300), xlabel = L"H\quad \textrm{(Oe)}", ylabel = L"M\quad (\mu_\textrm{B}/\textrm{ion})", ylims = (0,7), legend = false)

J = 7/2; #全角運動量
t1 = collect(0.1:0.2:0.9); #計算する温度（転移温度で規格化した温度 : T/T_C）
x = collect(0:0.1:30);
m1 =  ((2J+1)/2J) * coth.(((2J+1)/2J) * x) - (1/2J) * coth.(x/2J)
m2 = (t1' * (J+1)/3J) .* x

plot(x, [m1 m2])
plot!(size=(500, 300), xlabel = L"x", ylabel = L"M/M_0", legend = false, xlims = (0, 20), ylims = (0, 1.2))

using Roots

function FindRootCalc(J, t)
    m1(x) =  ((2J+1)/2J) * coth(((2J+1)/2J) * x) - (1/2J) * coth(x/2J)
    m2(x) = (t * (J+1)/3J) * x
    m3(x) = m1(x)-m2(x)
    sss = find_zero(m3, 100, Order16())
end

J = 7/2; #全角運動量
calcT = collect(0.01:0.01:1); #T/Tc
calcM = FindRootCalc.(J, calcT) .* calcT * (J+1)/3J; #M/M0
    
plot(calcT, calcM)
plot!(size=(500, 300), xlabel = L"T/T_c", ylabel = L"M/M_0", legend = false, xlims = (0, 1.1), ylims = (0, 1.1))

using CSV, DataFrames
function saveData2(xwave, y1wave, filename::String)
    df = DataFrame(TdivTc = xwave, MdivM0 = y1wave);
    df |> CSV.write(filename)
end

saveData2(calcT, calcM, "molecularfield.csv")