pdf("planck.pdf")

plot(lambda,lfl_T1,
    main="Planck's Law for a Blackbody",
    log="xy",
    xlab=expression(paste("Wavelength (",mu,"m)")),
    ylab=expression(paste(lambda,'F'[lambda]," (erg cm"^-2," s"^-1,")")),
    type="l",
    col="blue",
    xlim=c(0.01,10),
    ylim=c(1e6,1e16),
    xaxt='n',
    yaxt='n'
    )
lines(lambda,lfl_T2,col="green")
lines(lambda,lfl_T3,col="red")
legend(.5,1e8,legend = c("T = 150000 K","T = 40000 K","T = 5800 K"),lty=c(1,1,1),col=c("blue","green","red"))

axis(1,at=0.01,label="0.01")
axis(1,at=0.02,label="0.02")
axis(1,at=0.05,label="0.05")
axis(1,at=0.1,label="0.1")
axis(1,at=0.2,label="0.2")
axis(1,at=0.5,label="0.5")
axis(1,at=1,label="1")
axis(1,at=2,label="2")
axis(1,at=5,label="5")
axis(1,at=10,label="10")


axis(2,at=10^16,label=expression("10"^16))
axis(2,at=10^15,label=expression("10"^15))
axis(2,at=10^14,label=expression("10"^14))
axis(2,at=10^13,label=expression("10"^13))
axis(2,at=10^12,label=expression("10"^12))
axis(2,at=10^11,label=expression("10"^11))
axis(2,at=10^10,label=expression("10"^10))
axis(2,at=10^9,label=expression("10"^9))
axis(2,at=10^8,label=expression("10"^8))
axis(2,at=10^7,label=expression("10"^7))
axis(2,at=10^6,label=expression("10"^6))
axis(2,at=10^5,label=expression("10"^5))
axis(2,at=10^4,label=expression("10"^4))
axis(2,at=10^4,label=expression("10"^4))
axis(2,at=10^3,label=expression("10"^3))
axis(2,at=10^2,label=expression("10"^2))
axis(2,at=10^1,label=expression("10"^1))
axis(2,at=1,label=expression("1"))



dev.off()