data = read.table("mehdipour2013.tab")
plot (data,log="xy",type="l")
coords=read.table("mehdipour2013_samples.tab")
points(coords)


data=read.table("magdziarz1997.tab")
plot (data,log="xy",type="l")



# mehdipour zoomed with log lines
data = read.table("mehdipour2013.tab")
plot (data,log="xy",type="l",xlim=c(.001,300),ylim=c(1e11,1e13))
# mehdipour axes
abline(v=0.001)
abline(h=40e11)
coords=read.table("mehdipour2013_samples.tab")
points(coords,pch="|")
powers = seq(-10,10,by=1)
coefficients = seq(2,9,by=1)
for (i in powers) {
    abline(v=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(v=j*(10^i),col="black",lty=2)
    }
}
powers = seq(-45,45,by=1)
coefficients = seq(2,9,by=1)
for (i in powers) {
    abline(h=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(h=j*(10^i),col="black",lty=2)
    }
}






data = read.table("mehdipour2013.tab")
plot (data,log="xy",type="l")
# mehdipour axes
abline(v=0.001)
abline(h=40e11)

minorticks=



data=read.table("magdziarz1997.tab")
plot (data,log="xy",type="l",xlim=c(.001,300),ylim=c(1e-3,1e-1))
plot (data,log="xy",type="l")
coords=read.table("magdziarz1997_samples.tab")
points(coords)
#magdziarz boundaries
abline(v=.0004)
abline(v=.03)
abline(v=.4)
abline(v=3)
abline(v=50)
abline(v=500)

#low uv peak
abline(h=0.04)
abline(v=0.008)

#high uv peak
abline(h=0.024)
abline(v=0.012)

#xray peak
abline(h=0.03)
abline(v=60)

data=read.table("magdziarz1997.tab")
plot (data,log="xy",type="l",xlim=c(.001,300),ylim=c(1e-3,1e-1))
coords=read.table("magdziarz1997_samples.tab")
points(coords)

data=read.table("magdziarz1997.tab")
plot (data,log="xy",type="l")
coords=read.table("magdziarz1997_samples.tab")
points(coords)

data=read.table("mehdipour2013.tab")
plot (data,log="xy",type="l",xlim=c(.001,300),ylim=c(1e11,1e13))
coords=read.table("mehdipour2013_samples.tab")
points(coords)

data=read.table("mehdipour2013.tab")
plot (data,log="xy",type="l")
coords=read.table("mehdipour2013_samples.tab")
points(coords)


# model boundaries, if needed
abline(v=.0004,col="black")
abline(v=.03,col="black")
abline(v=.4,col="black")
abline(v=3,col="black")
abline(v=50,col="black")
abline(v=500,col="black")


abline(v=0.0004,col="red")
abline(v=0.006,col="red")
abline(v=.4,col="red")
abline(v=2,col="red")
abline(v=100,col="red")
abline(v=500,col="red")




# print pdfs of normalized SEDs
data1 = read.table("magdziarz_incident_continuum")
data2 = read.table("mehdipour_incident_continuum")
x1 = data1["V1"]
y1 = data1["V2"]
x2 = data2["V1"]
y2 = data2["V2"]
pdf("sed_overlay_with_boundaries.pdf")
plot(x1[,],y1[,],log="xy",type="l",xlim=c(1e-8,1e5),lwd=2,xlab=expression(paste("h",nu," [keV]")),ylab=expression(paste(nu,"F",nu," [erg cm^-2 s^-1]")))
lines(x2[,],y2[,],type="l",col="red",lwd=2)
powers = seq(-10,10,by=1)
coefficients = c(2,5)
for (i in powers) {
    abline(v=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(v=j*(10^i),col="black",lty=2)
    }
}

powers = seq(-45,45,by=1)
coefficients = c(2,5)
for (i in powers) {
    abline(h=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(h=j*(10^i),col="black",lty=2)
    }
}
dev.off()
pdf("sed_overlay_with_boundaries_zoomed.pdf")
plot(x1[,],y1[,],log="xy",type="l",lwd=2,xlim=c(0.001,1000),ylim=c(5e7,5e9),xlab=expression(paste("h",nu," [keV]")),ylab=expression(paste(nu,"F",nu," [erg cm^-2 s^-1]")))
lines(x2[,],y2[,],type="l",col="red",lwd=2)
powers = seq(-10,10,by=1)
coefficients = seq(2,9,by=1)
for (i in powers) {
    abline(v=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(v=j*(10^i),col="black",lty=2)
    }
}

powers = seq(-45,45,by=1)
coefficients = seq(2,9,by=1)
for (i in powers) {
    abline(h=10^i,col="black",lty=1)
    for (j in coefficients) {
        abline(h=j*(10^i),col="black",lty=2)
    }
}
dev.off()


# Logspace minor ticks






hc=9.113e2

inci = read.table("cont_combined_inci",header=TRUE)
diffuse = read.table("cont_combined_diffuse",header=TRUE)




# To plot continuum with bins, both tables already in angstroms

# prep tables
cut -f1,4 magdziarz_emitted_continuum_angstroms > cont 
cut -f1,3 magdziarz_continuum_bins_angstroms > bins

#R code
cont = read.table("cont",header=TRUE)
bins = read.table("bins",header=TRUE)
plot(cont,log="xy",type="b",xlim=c(1000,1200),ylim=c(1e7,1e9));abline(v=1420)


# old stuff
plot(Enr,nFn,log="x",xlim=c(20,25),type="l")
arrows(Enr-d.anu.,nFn,Enr+d.anu.,nFn,length=0.05,angle=90,code=3)
arrows(Enr-d.anu.,nFn+5,Enr+d.anu.,nFn-5,length=0.05,angle=90,code=3)
arrows(Enr-d.anu.,nFn+5,Enr+d.anu.,nFn+5,length=0.05,angle=90,code=3)
arrows(Enr-d.anu.,nFn+50,Enr+d.anu.,nFn+50,length=0.05,angle=90,code=3)
arrows(Enr-d.anu.,nFn+1e9,Enr+d.anu.,nFn+1e9,length=0.05,angle=90,code=3)



plot(energy.Ryd,Total,log="x",type="l",xlim=c(1e-2,5))
arrows(energy.Ryd-d.anu.,Total,energy.Ryd+d.anu.,Total,length=0.05,angle=90,code=3)
plot(energy.Ryd,ConEmitLocal,log="x",type="l",xlim=c(1e-2,5))

lambda=hc/energy.Ryd
width_lambda=lambda*d.anu./energy.Ryd
plot(lambda,Total,log="x",type="l",ylim=c(0,0.2),xlim=c(2000,9400))
arrows(lambda-width_lambda/2,Total,lambda+width_lambda/2,Total,length=0.05,angle=90,code=3)