using Plots
using LaTeXStrings
pgfplots()
Plots.PGFPlotsBackend()
nbar = 1
tbar = 1
nun = 0.3
nut = 1.2
n(rho) = nbar * ( 1 + nun ) * ( 1 - rho ^ 2 ) ^ nun
t(rho) = tbar * ( 1 + nut ) * ( 1 - rho ^ 2 ) ^ nut
return
jbar = 1
rho_m = 0.6
gamma = 1 / ( 1 - rho_m ^ 2 )
j(rho) = jbar * gamma^2 * ( 1 - rho ^ 2 ) * exp(gamma * rho^2) / ( exp(gamma) - 1 - gamma )
return
lambdat = 1.5
rhoped = 0.9
tped = tbar * 0.45
tsep = tped / 5
ktu = 2 * tped + tsep
ktu *= rhoped
ktu += tped + 2 * tsep
ktu *= ( 1 - rhoped )
ktu /= 3
ktu += tped * rhoped ^ 2
ktd = Base.gamma(1 + nut)
ktd *= Base.gamma(2/lambdat)
ktd /= 1 + nut + 2/lambdat
ktd *= 2
ktd /= lambdat
ktd *= rhoped ^ 2
to = tped + ( tbar - ktu ) / ktd
tpara(rho) = tped + ( to - tped ) * ( 1 - (rho/rhoped) ^lambdat ) ^ nut
tline(rho) = tsep + ( tped - tsep ) * ( 1 - rho ) / ( 1 - rhoped )
tpede(rho) = ( rho <= rhoped ? tpara(rho) : tline(rho) )
return
using QuadGK
pede_scaling = quadgk((rr) -> t(rr)*rr, 0,1)[1] / quadgk((rr) -> tpede(rr)*rr, 0,1)[1]
return
rho = linspace(0,1,251)
nn = n.(rho)
tt = t.(rho)
jj = j.(rho)
pp = pede_scaling*tpede.(rho)
return
temp_plot = plot()
plot!(rho, tt, color=2, label="T -- Parabolic")
plot!(rho, pp, color=5, label="T -- Pedestal")
plot!([0,rhoped], [tped, tped], color=5, label="", style=:dot)
ylims!(0,3)
title!("Temperature Profiles")
xlabel!("Normalized Minor Radius")
yticks!([1,2,3])
profile_plot = plot()
plot!(rho, nn, color=1, label="n -- Density")
plot!(rho, jj, color=3, label="J -- Current")
ylims!(0,3)
title!("Plasma Profiles")
xlabel!("Normalized Minor Radius")
yticks!([1,2,3])
savefig(temp_plot, "../images/temperature_profiles.tex")
savefig(profile_plot, "../images/plasma_profiles.tex")