update

sed/magdziarz1997/magdziarz1997_powlaw_coords.tab

@@ -1,6 +1,6 @@
 0.0004  0.0003
-0.03     0.028
+0.03    0.028
 0.4     0.024
-3       0.013
+3       0.012
 50      0.029
 500     0.00017

sed/mehdipour2013/mehdipour2013_powlaw_coords.tab

@@ -1,6 +1,6 @@
 0.0004  3.2e10
-0.006   4.9e13
+0.006   4.9e12
-.4      2.2e13
+.4      2.2e12
-2       1.4e13
+2       1.4e12
-100     7.7e13
+100     6.58e12
-500     1.5e13
+500     1.5e12

sed/mehdipour2013/mehdipour2013_samples.tab

@@ -1,7 +1,8 @@
-0.0013   9e12
+0.0013  9e11
-0.002   2e13
+0.002   2e12
-0.004   4.5e13
+0.004   4.5e12
-0.005   5e13
+0.005   5e12
-1.3     1.3e13
+1.3     1.3e12
-240     8.2e13
+200     8e12
-300     8e13les
+240     8.2e12
+300     8e12

sed/plots.R

@@ -1,5 +1,8 @@
 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")
@@ -8,10 +11,12 @@ plot (data,log="xy",type="l")
 
 
 data = read.table("mehdipour2013.tab")
-plot (data,log="xy",type="l",xlim=c(.001,300),ylim=c(10e11,10e13))
+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)
 
 data = read.table("mehdipour2013.tab")
 plot (data,log="xy",type="l")
@@ -19,10 +24,15 @@ plot (data,log="xy",type="l")
 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)
@@ -43,9 +53,99 @@ abline(v=0.012)
 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")
+
+
+
+
+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-6,1e7),lwd=2,xlab="eV",ylab="nuFnu")
+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="eV",ylab="nuFnu")
+lines(x2[,],y2[,],type="l",col="red",lwd=2)
+powers = seq(-10,10,by=1)
+coefficients = seq(2,8,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,8,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
+
+
 
 
 

src/sed.hpp

@@ -38,9 +38,11 @@ struct sed_table {
 // To account for the four main powerlaws in a typical
 // AGN SED.
 
-// Hardcoded infrared and gamma ray power laws.
+// Hardcoded infrared and gamma ray power laws and cutoffs.
 const double IR_POWER = 3;
-const double GAMMA_POWER = -2;
+const double GAMMA_POWER = -5;
+const double RADIO_CUTOFF = 1e-4; // IN KEV
+const double GAMMA_CUTOFF = 1e4;    // IN KEV
 
 struct powerlaw_bounds {
     double ir_min;
@@ -68,7 +70,14 @@ public:
         _power(power),
         _normal(exp(log(x0.second)-(_power*log(x0.first))))
         {}
-    double eval(double hnu) { return _normal*pow(hnu,_power); }
+    double eval(double hnu) { 
+        return 
+            _normal
+            * pow(hnu,_power)
+            * exp(-(hnu)/GAMMA_CUTOFF)
+            * exp(-(RADIO_CUTOFF/hnu))
+        ;
+    }
 };
 
 class sed {