using Plots, BenchmarkTools

#NMRパラメータの設定
H0 = 3.0; #(T)外部磁場
gamma0 = 17.235; #(MHz/T)核磁気回転比

#等方的ナイトシフト
Kiso = 1.5; #(%)
#異方的ナイトシフト。Kanx+Kany+Kanz=0, |Kanz|>|Kanx|>|Kany|
Kanz = 2.4;#(%)
Kanx = -1.7;#(%)
Kany = -0.7;#(%)

#計算周波数範囲の設定
fmin0 = 51;#MHz
fmax0 = 54;#MHz
df0   = 0.01;#MHz

nu0 = H0*gamma0; #(MHz) K=0の共鳴周波数
KX = (Kiso+Kanx)/100;
KY = (Kiso+Kany)/100;
KZ = (Kiso+Kanz)/100;

calfreq = collect(fmin0:df0:fmax0);#周波数の列を作る
nagasa = length(calfreq);
answave = zeros(nagasa, 2); #周波数分の列で2行の行列を作る
answave[:,1] = calfreq;           #一列目に周波数を代入

function main()
    freq1 = [0.0]
    int1 = [0.0]
    for θ = 0:0.002:π/2
        for φ = 0:0.02:π/2
            freq0 = nu0*(1.0 + KZ*cos(θ)^2+sin(θ)^2*(KX*cos(φ)^2+KY*sin(φ)^2))#共鳴周波数
            int0 = abs(sin(θ)) #遷移確率はθに比例
            freq1 = push!(freq1, freq0)
            int1 = push!(int1, int0)
        end
    end
    for jjj = 1:nagasa
        answave[jjj,2] = answave[jjj, 2] + sum(int1 .* (freq1 .<= (answave[jjj, 1] + df0/2)) .* (freq1 .> (answave[jjj,1] - df0/2)))
    end
end

@time main()
plot(answave[:,1], answave[:,2], fill = (0, :skyblue))
plot!(xlabel="Frequency (MHz)", ylabel="Intensity", legend=false, size=(400,300))

#NMRパラメータの設定
H0 = 3.0; #(T)外部磁場
gamma0 = 17.235; #(MHz/T)核磁気回転比

#等方的ナイトシフト
Kiso = 1.0; #(%)
#異方的ナイトシフト
Kan = 1.0;#(%)

#計算周波数範囲の設定
fmin0 = 51;#MHz
fmax0 = 54;#MHz
df0   = 0.01;#MHz

nu0 = H0*gamma0; #(MHz) K=0の共鳴周波数
Kiso = Kiso/100;#%の入力を処理
Kan = Kan/100;#%の入力を処理

calfreq = collect(fmin0:df0:fmax0);#周波数の列を作る
nagasa = length(calfreq);
answave = zeros(nagasa, 2); #周波数分の列で2行の行列を作る
answave[:,1] = calfreq;           #一列目に周波数を代入\

function main()
    freq1 = [0.0]
    int1 = [0.0]
    for θ = 0:0.001:π/2
        freq0 = nu0*(1.0 + Kiso + Kan*(3*cos(θ)^2 - 1))#共鳴周波数
        int0 = abs(sin(θ)) #遷移確率はθに比例
        freq1 = push!(freq1, freq0)
        int1 = push!(int1, int0)
    end
    for jjj = 1:nagasa
        answave[jjj,2] = answave[jjj, 2] + sum(int1 .* (freq1 .<= (answave[jjj, 1] + df0/2)) .* (freq1 .> (answave[jjj,1] - df0/2)))
    end
end

@time main()
plot(answave[:,1], answave[:,2], fill = (0, :skyblue))
plot!(xlabel="Frequency (MHz)", ylabel="Intensity", legend=false, size=(400,300))