#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>

double l_band_HRH(double log_l_bol, double nu);
double return_ratio_to_hard_xray(double nu);
double return_ratio_to_b_band(double nu);
double l_band_GTR(double log_L_bol, double nu);
double l_band_MAR(double log_l_bol, double nu);
double ratio_of_vandenberk_to_continuum(double nu);

/* 
	Script to return the model intrinsic (un-reddened/un-obscured) quasar spectrum
		as a function of frequency and bolometric quasar luminosity, constructed from
		observations from ~300 microns to ~100 keV in Hopkins, Richards, 
		and Hernquist 2006, ApJ (HRH06)

	  These model spectra are compiled in the paper above, which details the 
	    observational data sets, but in short, we draw heavily from the 
        observations in a number of papers, including :
        Richards et al. 2006, Steffen et al. 2006, Hatziminaoglou et al. 2005, 
        Tozzi et al. 2006, Strateva et al. 2005, Telfer et al. 2002, Vanden Berk et al. 2001, 
        George et al. 1998, Elvis et al. 1994
		also, Marconi et al. 2004 (we follow a similar methodology to their derivation)
		and note that we generate the X-ray reflection component 
          following Ueda et al. 2003 in the PEXRAV code, Magdziarz & Zdziarski 1995

	number of input parameters : 1-2
		1 : log_{10}(L_bol/L_sun) -- the log of the bolometric quasar luminosity 
				for which the mean spectrum is desired. This uses the model 
				constructed from the observations above in HRH06, in which the mean
				intrinsic quasar SED is a function only of bolometric luminosity. 
				units are :
					L_sun = the bolometric solar luminosity = 3.9 x 10^33 erg/s
		2 : SDSS_KEY -- (OPTIONAL) keyword. By default, the template spectrum of 
				vanden Berk et al. 2001 from SDSS is overlaid on the optical portion 
				of the spectrum, as described in HRH06. If for some reason you 
				want to turn this feature *off* (for example, you need to smooth 
				over or otherwise don't want to include line details), then 
				add a second input parameter, set to 0

	results are simply printed to stdout, in three columns. The header 
	  printed out will describe them, in the general format : 
	  
			;;//         
			;;// Intrinsic (un-reddened) AGN spectrum for 
			;;//   log_{10}(L_bol/L_sun)   = xx.xx
			;;//   log_{10}(L_bol/[erg/s]) = xx.xx 
			;;//     (throughout, L_sun is the bolometric solar L = 3.9x10^33 erg/s)  
			;;// columns are frequency, AGN spectrum (nu*L_nu) in solar and CGS units :  
			;;//   log_{10}(nu/Hz)   log_{10}(nu*L_nu/L_sun)  log_{10}(nu*L_nu/[erg/s])  
			;;// the output frequency list spans ~300 microns to ~100 keV  
			;;// the last four columns compute bands frequently adopted :  
			;;//   nu = -1.00  ::  B-band (4400 Angstroms)  
			;;//   nu = -2.00  ::  mid-IR (15 microns)  
			;;//   nu = -3.00  ::  soft X-ray (integrated luminosity from 0.5-2 keV)  
			;;//   nu = -4.00  ::  hard X-ray (integrated luminosity from 2-10  keV)  
			;;//         


	Example use : 
        compile       ::    > gcc agn_spectrum.c -o agn_spec

        standard call ::    > agn_spec 14.0
             (print the template spectrum for a L_bol = 10^14 L_sun quasar to stdout)

        dump to file  ::    > agn_spec 14.0 > filename
             (the output results above get dumped to 'filename')

        turn off the SDSS template overlay  ::    > agn_spec 14.0 0
             (as the above, but adds the second input set to 0, to turn off the 
                overlay of the SDSS model spectrum in the optical)

*/


int main(int argc, char** argv) 
{
	double nu, log_L_bol, l_band, *l_band_vec;
	char *MODEL_SPECTRUM_SET;
	char *band_check;
	char *SLOAN_SET;
	int BAND_KEY=0;	int MODEL_SPECTRUM_KEY; int SLOAN_KEY=0; int i_nu;
		// sets if simply calling the bolometric correction for a prefab band : 
		//		BB = B-band, IR = 15 micron, SX = 0.5-2 keV, HX = 2-10 keV	
  	if (argc < 2) {
  	  fprintf(stderr, "Expected >=1 arguments (found %d)\n",argc-1); 
  	  exit(0); 
  	}

int N_nu = 375;	
double nu_vec[375] = {
12.00, 12.02, 12.04, 12.06, 12.08, 12.10, 12.12, 12.14, 12.16, 12.18, 12.20, 12.22, 12.24,
12.26, 12.28, 12.30, 12.32, 12.34, 12.36, 12.38, 12.40, 12.42, 12.44, 12.46, 12.48, 
12.50, 12.52, 12.54, 12.56, 12.58, 12.60, 12.62, 12.64, 12.66, 12.68, 12.70, 12.72, 12.74,
12.76, 12.78, 12.80, 12.82, 12.84, 12.86, 12.88, 12.90, 12.92, 12.94, 12.96, 12.98, 13.00,
13.02, 13.04, 13.06, 13.08, 13.10, 13.12, 13.14, 13.16, 13.18, 13.20, 13.22, 13.24, 13.26,
13.28, 13.30, 13.32, 13.34, 13.36, 13.38, 13.40, 13.42, 13.44, 13.46, 13.48, 13.50, 13.52,
13.54, 13.56, 13.58, 13.60, 13.62, 13.64, 13.66, 13.68, 13.70, 13.72, 13.74, 13.76, 13.78,
13.80, 13.82, 13.84, 13.86, 13.88, 13.90, 13.92, 13.94, 13.96, 13.98, 14.00, 14.02, 14.04,
14.06, 14.08, 14.10, 14.12, 14.14, 14.16, 14.18, 14.20, 14.22, 14.24, 14.26, 14.28, 14.30,
14.32, 14.34, 14.36, 14.38, 14.40, 14.42, 14.44, 14.46, 14.48, 14.50, 14.52, 14.54, 14.56,
14.58, 14.60, 14.62, 14.64, 14.66, 14.68, 14.70, 14.72, 14.74, 14.76, 14.78, 14.80, 14.82,
14.84, 14.86, 14.88, 14.90, 14.92, 14.94, 14.96, 14.98, 15.00, 15.02, 15.04, 15.06, 15.08,
15.10, 15.12, 15.14, 15.16, 15.18, 15.20, 15.22, 15.24, 15.26, 15.28, 15.30, 15.32, 15.34,
15.36, 15.38, 15.40, 15.42, 15.44, 15.46, 15.48, 15.50, 15.52, 15.54, 15.56, 15.58, 15.60,
15.62, 15.64, 15.66, 15.68, 15.70, 15.72, 15.74, 15.76, 15.78, 15.80, 15.82, 15.84, 15.86,
15.88, 15.90, 15.92, 15.94, 15.96, 15.98, 16.00, 16.02, 16.04, 16.06, 16.08, 16.10, 16.12,
16.14, 16.16, 16.18, 16.20, 16.22, 16.24, 16.26, 16.28, 16.30, 16.32, 16.34, 16.36, 16.38,
16.40, 16.42, 16.44, 16.46, 16.48, 16.50, 16.52, 16.54, 16.56, 16.58, 16.60, 16.62, 16.64,
16.66, 16.68, 16.70, 16.72, 16.74, 16.76, 16.78, 16.80, 16.82, 16.84, 16.86, 16.88, 16.90,
16.92, 16.94, 16.96, 16.98, 17.00,
17.02, 17.04, 17.06, 17.08, 17.10, 17.12, 17.14, 17.16, 17.18, 17.20, 17.22, 17.24, 17.26,
17.28, 17.30, 17.32, 17.34, 17.36, 17.38, 17.40, 17.42, 17.44, 17.46, 17.48, 17.50, 17.52,
17.54, 17.56, 17.58, 17.60, 17.62, 17.64, 17.66, 17.68, 17.70, 17.72, 17.74, 17.76, 17.78,
17.80, 17.82, 17.84, 17.86, 17.88, 17.90, 17.92, 17.94, 17.96, 17.98, 18.00, 18.02, 18.04,
18.06, 18.08, 18.10, 18.12, 18.14, 18.16, 18.18, 18.20, 18.22, 18.24, 18.26, 18.28, 18.30,
18.32, 18.34, 18.36, 18.38, 18.40, 18.42, 18.44, 18.46, 18.48, 18.50, 18.52, 18.54, 18.56,
18.58, 18.60, 18.62, 18.64, 18.66, 18.68, 18.70, 18.72, 18.74, 18.76, 18.78, 18.80, 18.82,
18.84, 18.86, 18.88, 18.90, 18.92, 18.94, 18.96, 18.98, 19.00, 19.02, 19.04, 19.06, 19.08,
19.10, 19.12, 19.14, 19.16, 19.18, 19.20, 19.22, 19.24, 19.26, 19.28, 19.30, 19.32, 19.34,
19.36, 19.38, 19.40, -1.00, -2.00, -3.00, -4.00 };
l_band_vec = calloc(N_nu,sizeof(double));

log_L_bol = atof(argv[1]);
SLOAN_KEY = 1;
if (argc >= 3) {SLOAN_KEY  = atoi(argv[2]);}	
   if ((SLOAN_KEY != 0)&&(SLOAN_KEY != 1)) {SLOAN_KEY=0;}
  // if specified, can turn on/off the superposition of the SDSS mean spectrum (default = ON)


	MODEL_SPECTRUM_KEY = 0;  // switch to 1 for the Richards et al. 2006 luminosity-independent spectrum
	for (i_nu=0;i_nu<N_nu;i_nu++) {
	nu = pow(10.,nu_vec[i_nu]);
	if (nu_vec[i_nu]<=0.) {nu=nu_vec[i_nu];}
	if (MODEL_SPECTRUM_KEY==0) l_band = log10(l_band_HRH(log_L_bol,nu));
	if (MODEL_SPECTRUM_KEY==1) l_band = l_band_GTR(log_L_bol,nu);
	double log_nu;
		if (nu <= 0.) {
			if (nu==-1.0) log_nu = log10(2.998e8/(4.4e-7));
			if (nu==-2.0) log_nu = log10(2.998e8/(1.5e-5));
			if (nu==-3.0) log_nu = log10(1.0 * 2.418e17);
			if (nu==-4.0) log_nu = log10(5.0 * 2.418e17);
		} else {log_nu = log10(nu);}
	// check if inside range of Vanden Berk spectrum (note, this is not the 
	//    exact range over which that spectrum is compiled, since there is 
	//    serious instrumental cutoff at the extreme edges)
	if ((SLOAN_KEY==1)&&(log_nu>14.575)&&(log_nu<15.500)) 
		l_band += log10(ratio_of_vandenberk_to_continuum(nu));
	l_band_vec[i_nu] = l_band;
	}


	// now print the resulting spectrum
	double to_cgs = log10(3.9) + 33.0;
 	 printf(";;//         \n");
 	 printf(";;// Intrinsic (un-reddened) AGN spectrum for \n");
 	 printf(";;//   log_{10}(L_bol/L_sun)   = %5.2f \n",log_L_bol);
 	 printf(";;//   log_{10}(L_bol/[erg/s]) = %5.2f \n",log_L_bol+to_cgs);
 	 printf(";;//     (throughout, L_sun is the bolometric solar L = 3.9x10^33 erg/s)  \n",log_L_bol);
 	 printf(";;// columns are frequency, AGN spectrum (nu*L_nu) in solar and CGS units :  \n");
 	 printf(";;//   log_{10}(nu/Hz)   log_{10}(nu*L_nu/L_sun)  log_{10}(nu*L_nu/[erg/s])  \n");
 	 printf(";;// the output frequency list spans ~300 microns to ~100 keV  \n");
 	 printf(";;// the last four rows compute bands frequently adopted :  \n");
 	 printf(";;//   nu = -1.00  ::  B-band (4400 Angstroms)  \n");
 	 printf(";;//   nu = -2.00  ::  mid-IR (15 microns)  \n");
 	 printf(";;//   nu = -3.00  ::  soft X-ray (integrated luminosity from 0.5-2 keV)  \n");
 	 printf(";;//   nu = -4.00  ::  hard X-ray (integrated luminosity from 2-10  keV)  \n");
 	 printf(";;//         \n");	
    for (i_nu=0;i_nu<N_nu;i_nu++) printf(" %8.2f %10.5f %10.5f \n",nu_vec[i_nu],l_band_vec[i_nu],l_band_vec[i_nu]+to_cgs);

	return 0;
}


// code to return the vanden Berk mean spectrum -- has been carefully gap-repaired 
//  and the low-nu continuum recalibrated following Haziminaglou et al. 2006 and 
//  Richards et al. 2006 to give the correct continuum slope where the geometric 
//  averaging routine was flattening things out -- then clipped at the edges for 
//  instrumental/IGM absorption cutoff at the extreme high-nu/low-nu ends, 
//  and renormalized to the boxcar continuum -- so here is the final ratio to 
//  continuum over the appropriate nu range, normalized s.t. total luminosity in 
//  this range is conserved -- can then be overlayed on any template spectrum
double ratio_of_vandenberk_to_continuum(double nu)
{
if (nu<= 0.) {
	if (nu==-1.0) return 1.0147898;
	// the other prefab bands are clearly outside the appopriate range for this overlay
	return 1.0;
}
double log_nu[926] = {
14.575, 14.576, 14.577, 14.578, 14.579, 14.580, 14.581, 14.582, 14.583, 14.584, 14.585,
14.586, 14.587, 14.588, 14.589, 14.590, 14.591, 14.592, 14.593, 14.594, 14.595, 14.596,
14.597, 14.598, 14.599, 14.600, 14.601, 14.602, 14.603, 14.604, 14.605, 14.606, 14.607,
14.608, 14.609, 14.610, 14.611, 14.612, 14.613, 14.614, 14.615, 14.616, 14.617, 14.618,
14.619, 14.620, 14.621, 14.622, 14.623, 14.624, 14.625, 14.626, 14.627, 14.628, 14.629,
14.630, 14.631, 14.632, 14.633, 14.634, 14.635, 14.636, 14.637, 14.638, 14.639, 14.640,
14.641, 14.642, 14.643, 14.644, 14.645, 14.646, 14.647, 14.648, 14.649, 14.650, 14.651,
14.652, 14.653, 14.654, 14.655, 14.656, 14.657, 14.658, 14.659, 14.660, 14.661, 14.662,
14.663, 14.664, 14.665, 14.666, 14.667, 14.668, 14.669, 14.670, 14.671, 14.672, 14.673,
14.674, 14.675, 14.676, 14.677, 14.678, 14.679, 14.680, 14.681, 14.682, 14.683, 14.684,
14.685, 14.686, 14.687, 14.688, 14.689, 14.690, 14.691, 14.692, 14.693, 14.694, 14.695,
14.696, 14.697, 14.698, 14.699, 14.700, 14.701, 14.702, 14.703, 14.704, 14.705, 14.706,
14.707, 14.708, 14.709, 14.710, 14.711, 14.712, 14.713, 14.714, 14.715, 14.716, 14.717,
14.718, 14.719, 14.720, 14.721, 14.722, 14.723, 14.724, 14.725, 14.726, 14.727, 14.728,
14.729, 14.730, 14.731, 14.732, 14.733, 14.734, 14.735, 14.736, 14.737, 14.738, 14.739,
14.740, 14.741, 14.742, 14.743, 14.744, 14.745, 14.746, 14.747, 14.748, 14.749, 14.750,
14.751, 14.752, 14.753, 14.754, 14.755, 14.756, 14.757, 14.758, 14.759, 14.760, 14.761,
14.762, 14.763, 14.764, 14.765, 14.766, 14.767, 14.768, 14.769, 14.770, 14.771, 14.772,
14.773, 14.774, 14.775, 14.776, 14.777, 14.778, 14.779, 14.780, 14.781, 14.782, 14.783,
14.784, 14.785, 14.786, 14.787, 14.788, 14.789, 14.790, 14.791, 14.792, 14.793, 14.794,
14.795, 14.796, 14.797, 14.798, 14.799, 14.800, 14.801, 14.802, 14.803, 14.804, 14.805,
14.806, 14.807, 14.808, 14.809, 14.810, 14.811, 14.812, 14.813, 14.814, 14.815, 14.816,
14.817, 14.818, 14.819, 14.820, 14.821, 14.822, 14.823, 14.824, 14.825, 14.826, 14.827,
14.828, 14.829, 14.830, 14.831, 14.832, 14.833, 14.834, 14.835, 14.836, 14.837, 14.838,
14.839, 14.840, 14.841, 14.842, 14.843, 14.844, 14.845, 14.846, 14.847, 14.848, 14.849,
14.850, 14.851, 14.852, 14.853, 14.854, 14.855, 14.856, 14.857, 14.858, 14.859, 14.860,
14.861, 14.862, 14.863, 14.864, 14.865, 14.866, 14.867, 14.868, 14.869, 14.870, 14.871,
14.872, 14.873, 14.874, 14.875, 14.876, 14.877, 14.878, 14.879, 14.880, 14.881, 14.882,
14.883, 14.884, 14.885, 14.886, 14.887, 14.888, 14.889, 14.890, 14.891, 14.892, 14.893,
14.894, 14.895, 14.896, 14.897, 14.898, 14.899, 14.900, 14.901, 14.902, 14.903, 14.904,
14.905, 14.906, 14.907, 14.908, 14.909, 14.910, 14.911, 14.912, 14.913, 14.914, 14.915,
14.916, 14.917, 14.918, 14.919, 14.920, 14.921, 14.922, 14.923, 14.924, 14.925, 14.926,
14.927, 14.928, 14.929, 14.930, 14.931, 14.932, 14.933, 14.934, 14.935, 14.936, 14.937,
14.938, 14.939, 14.940, 14.941, 14.942, 14.943, 14.944, 14.945, 14.946, 14.947, 14.948,
14.949, 14.950, 14.951, 14.952, 14.953, 14.954, 14.955, 14.956, 14.957, 14.958, 14.959,
14.960, 14.961, 14.962, 14.963, 14.964, 14.965, 14.966, 14.967, 14.968, 14.969, 14.970,
14.971, 14.972, 14.973, 14.974, 14.975, 14.976, 14.977, 14.978, 14.979, 14.980, 14.981,
14.982, 14.983, 14.984, 14.985, 14.986, 14.987, 14.988, 14.989, 14.990, 14.991, 14.992,
14.993, 14.994, 14.995, 14.996, 14.997, 14.998, 14.999, 15.000, 15.001, 15.002, 15.003,
15.004, 15.005, 15.006, 15.007, 15.008, 15.009, 15.010, 15.011, 15.012, 15.013, 15.014,
15.015, 15.016, 15.017, 15.018, 15.019, 15.020, 15.021, 15.022, 15.023, 15.024, 15.025,
15.026, 15.027, 15.028, 15.029, 15.030, 15.031, 15.032, 15.033, 15.034, 15.035, 15.036,
15.037, 15.038, 15.039, 15.040, 15.041, 15.042, 15.043, 15.044, 15.045, 15.046, 15.047,
15.048, 15.049, 15.050, 15.051, 15.052, 15.053, 15.054, 15.055, 15.056, 15.057, 15.058,
15.059, 15.060, 15.061, 15.062, 15.063, 15.064, 15.065, 15.066, 15.067, 15.068, 15.069,
15.070, 15.071, 15.072, 15.073, 15.074, 15.075, 15.076, 15.077, 15.078, 15.079, 15.080,
15.081, 15.082, 15.083, 15.084, 15.085, 15.086, 15.087, 15.088, 15.089, 15.090, 15.091,
15.092, 15.093, 15.094, 15.095, 15.096, 15.097, 15.098, 15.099, 15.100, 15.101, 15.102,
15.103, 15.104, 15.105, 15.106, 15.107, 15.108, 15.109, 15.110, 15.111, 15.112, 15.113,
15.114, 15.115, 15.116, 15.117, 15.118, 15.119, 15.120, 15.121, 15.122, 15.123, 15.124,
15.125, 15.126, 15.127, 15.128, 15.129, 15.130, 15.131, 15.132, 15.133, 15.134, 15.135,
15.136, 15.137, 15.138, 15.139, 15.140, 15.141, 15.142, 15.143, 15.144, 15.145, 15.146,
15.147, 15.148, 15.149, 15.150, 15.151, 15.152, 15.153, 15.154, 15.155, 15.156, 15.157,
15.158, 15.159, 15.160, 15.161, 15.162, 15.163, 15.164, 15.165, 15.166, 15.167, 15.168,
15.169, 15.170, 15.171, 15.172, 15.173, 15.174, 15.175, 15.176, 15.177, 15.178, 15.179,
15.180, 15.181, 15.182, 15.183, 15.184, 15.185, 15.186, 15.187, 15.188, 15.189, 15.190,
15.191, 15.192, 15.193, 15.194, 15.195, 15.196, 15.197, 15.198, 15.199, 15.200, 15.201,
15.202, 15.203, 15.204, 15.205, 15.206, 15.207, 15.208, 15.209, 15.210, 15.211, 15.212,
15.213, 15.214, 15.215, 15.216, 15.217, 15.218, 15.219, 15.220, 15.221, 15.222, 15.223,
15.224, 15.225, 15.226, 15.227, 15.228, 15.229, 15.230, 15.231, 15.232, 15.233, 15.234,
15.235, 15.236, 15.237, 15.238, 15.239, 15.240, 15.241, 15.242, 15.243, 15.244, 15.245,
15.246, 15.247, 15.248, 15.249, 15.250, 15.251, 15.252, 15.253, 15.254, 15.255, 15.256,
15.257, 15.258, 15.259, 15.260, 15.261, 15.262, 15.263, 15.264, 15.265, 15.266, 15.267,
15.268, 15.269, 15.270, 15.271, 15.272, 15.273, 15.274, 15.275, 15.276, 15.277, 15.278,
15.279, 15.280, 15.281, 15.282, 15.283, 15.284, 15.285, 15.286, 15.287, 15.288, 15.289,
15.290, 15.291, 15.292, 15.293, 15.294, 15.295, 15.296, 15.297, 15.298, 15.299, 15.300,
15.301, 15.302, 15.303, 15.304, 15.305, 15.306, 15.307, 15.308, 15.309, 15.310, 15.311,
15.312, 15.313, 15.314, 15.315, 15.316, 15.317, 15.318, 15.319, 15.320, 15.321, 15.322,
15.323, 15.324, 15.325, 15.326, 15.327, 15.328, 15.329, 15.330, 15.331, 15.332, 15.333,
15.334, 15.335, 15.336, 15.337, 15.338, 15.339, 15.340, 15.341, 15.342, 15.343, 15.344,
15.345, 15.346, 15.347, 15.348, 15.349, 15.350, 15.351, 15.352, 15.353, 15.354, 15.355,
15.356, 15.357, 15.358, 15.359, 15.360, 15.361, 15.362, 15.363, 15.364, 15.365, 15.366,
15.367, 15.368, 15.369, 15.370, 15.371, 15.372, 15.373, 15.374, 15.375, 15.376, 15.377,
15.378, 15.379, 15.380, 15.381, 15.382, 15.383, 15.384, 15.385, 15.386, 15.387, 15.388,
15.389, 15.390, 15.391, 15.392, 15.393, 15.394, 15.395, 15.396, 15.397, 15.398, 15.399,
15.400, 15.401, 15.402, 15.403, 15.404, 15.405, 15.406, 15.407, 15.408, 15.409, 15.410,
15.411, 15.412, 15.413, 15.414, 15.415, 15.416, 15.417, 15.418, 15.419, 15.420, 15.421,
15.422, 15.423, 15.424, 15.425, 15.426, 15.427, 15.428, 15.429, 15.430, 15.431, 15.432,
15.433, 15.434, 15.435, 15.436, 15.437, 15.438, 15.439, 15.440, 15.441, 15.442, 15.443,
15.444, 15.445, 15.446, 15.447, 15.448, 15.449, 15.450, 15.451, 15.452, 15.453, 15.454,
15.455, 15.456, 15.457, 15.458, 15.459, 15.460, 15.461, 15.462, 15.463, 15.464, 15.465,
15.466, 15.467, 15.468, 15.469, 15.470, 15.471, 15.472, 15.473, 15.474, 15.475, 15.476,
15.477, 15.478, 15.479, 15.480, 15.481, 15.482, 15.483, 15.484, 15.485, 15.486, 15.487,
15.488, 15.489, 15.490, 15.491, 15.492, 15.493, 15.494, 15.495, 15.496, 15.497, 15.498,
15.499, 15.500
};
double nuLnu[926] = {
0.915,  0.956,  0.960,  0.938,  0.969,  0.969,  0.999,  0.915,  0.944, 0.930,  0.928, 
0.935,  0.896,  0.985,  0.932,  0.931,  1.004,  0.957,  0.927,  0.972, 0.939,  1.018, 
0.949,  0.947,  0.981,  0.964,  0.921,  0.972,  0.946,  1.014,  0.918, 0.940,  0.965, 
0.992,  0.946,  0.926,  0.967,  0.958,  1.002,  0.925,  0.945,  0.937, 0.933,  0.932, 
0.945,  0.908,  0.940,  0.917,  0.939,  0.929,  0.953,  0.984,  1.008, 0.983,  0.992, 
0.975,  0.954,  0.973,  0.927,  0.955,  0.931,  0.943,  0.952,  0.941, 0.942,  0.932, 
0.965,  0.984,  0.952,  0.959,  0.945,  0.962,  0.954,  0.977,  1.078, 1.040,  1.059, 
1.073,  1.136,  1.198,  1.272,  1.457,  1.742,  2.781,  3.065,  3.514, 2.354,  1.792, 
1.435,  1.195,  1.191,  1.114,  1.080,  1.076,  1.025,  1.012,  0.990, 0.997,  1.019, 
0.998,  0.985,  0.945,  0.982,  0.957,  0.955,  0.929,  0.952,  0.935, 0.901,  0.914, 
0.932,  0.921,  0.923,  0.941,  0.932,  0.941,  0.932,  0.973,  0.974, 0.955,  0.965, 
0.960,  0.935,  0.888,  0.953,  0.944,  0.960,  0.957,  0.934,  0.966, 0.994,  0.977, 
1.038,  1.059,  0.975,  0.985,  0.969,  0.953,  0.962,  0.918,  0.925, 0.920,  0.927, 
0.923,  0.956,  0.918,  0.931,  0.925,  0.927,  0.904,  0.896,  0.925, 0.906,  0.908, 
0.892,  0.916,  0.914,  0.935,  0.934,  0.955,  0.925,  0.917,  0.908, 0.898,  0.919, 
0.912,  0.941,  0.932,  0.942,  0.919,  0.941,  0.953,  0.979,  0.979, 0.976,  0.996, 
1.042,  1.033,  1.047,  1.021,  0.977,  1.006,  1.019,  0.988,  0.997, 0.984,  1.023, 
0.973,  0.982,  1.031,  0.980,  0.966,  0.944,  0.939,  0.927,  0.951, 0.930,  0.960, 
0.981,  0.981,  1.035,  1.102,  2.217,  1.243,  1.047,  1.048,  1.270, 1.192,  1.111, 
1.163,  1.168,  1.228,  1.276,  1.389,  1.529,  1.844,  1.527,  1.356, 1.222,  1.148, 
1.083,  1.059,  1.011,  1.001,  0.996,  1.002,  1.010,  0.996,  1.010, 1.000,  1.027, 
1.084,  1.058,  1.063,  1.067,  1.097,  1.122,  1.113,  1.099,  1.082, 1.128,  1.133, 
1.132,  1.132,  1.112,  1.076,  1.106,  1.102,  1.102,  1.058,  1.042, 1.036,  1.003, 
0.995,  0.972,  0.973,  0.992,  1.016,  1.011,  1.021,  1.038,  1.089, 1.202,  1.179, 
1.253,  1.168,  1.102,  1.035,  0.984,  0.979,  0.960,  0.926,  0.931, 0.930,  0.943, 
0.928,  0.901,  0.897,  0.881,  0.899,  0.916,  0.905,  0.904,  0.893, 0.896,  0.918, 
0.943,  0.957,  0.969,  0.986,  0.971,  0.931,  0.929,  0.902,  0.889, 0.886,  0.873, 
0.887,  0.888,  0.883,  0.884,  0.903,  0.912,  0.930,  0.931,  0.915, 0.897,  0.858, 
0.903,  0.909,  0.914,  0.928,  0.941,  0.939,  1.049,  0.955,  0.937, 0.909,  0.891, 
0.902,  0.896,  0.894,  0.872,  0.898,  0.917,  0.915,  0.932,  0.949, 0.918,  0.921, 
1.089,  0.944,  0.916,  0.915,  0.914,  0.912,  0.904,  0.909,  0.911, 0.910,  0.909, 
0.915,  0.922,  0.912,  0.911,  0.908,  0.911,  0.912,  0.908,  0.904, 0.897,  0.893, 
0.889,  0.892,  0.898,  0.913,  0.922,  0.933,  0.942,  0.941,  0.942, 0.946,  0.944, 
0.937,  0.934,  0.936,  0.942,  1.000,  0.955,  0.938,  0.935,  0.919, 0.917,  0.919, 
0.906,  0.901,  0.900,  0.917,  0.910,  0.905,  0.896,  0.898,  0.906, 0.918,  0.913, 
0.919,  0.929,  0.922,  0.924,  0.922,  0.915,  0.919,  0.918,  0.925, 0.924,  0.933, 
0.932,  0.939,  0.943,  0.926,  0.916,  0.911,  0.899,  0.894,  0.892, 0.892,  0.905, 
0.895,  0.877,  0.863,  0.849,  0.839,  0.841,  0.843,  0.839,  0.836, 0.831,  0.829, 
0.828,  0.825,  0.819,  0.817,  0.825,  0.831,  0.850,  0.865,  0.880, 0.891,  0.897, 
0.899,  0.905,  0.911,  0.907,  0.907,  0.902,  0.901,  0.885,  0.878, 0.870,  0.876, 
0.882,  0.895,  0.899,  0.913,  0.927,  0.934,  0.947,  0.965,  0.976, 1.004,  1.014, 
1.049,  1.123,  1.239,  1.348,  1.368,  1.265,  1.140,  1.056,  1.013, 0.983,  0.973, 
0.967,  0.957,  0.936,  0.927,  0.914,  0.904,  0.894,  0.884,  0.875, 0.870,  0.863, 
0.854,  0.860,  0.861,  0.851,  0.844,  0.844,  0.855,  0.866,  0.882, 0.888,  0.892, 
0.895,  0.893,  0.891,  0.887,  0.892,  0.894,  0.900,  0.901,  0.905, 0.909,  0.911, 
0.914,  0.920,  0.923,  0.925,  0.925,  0.928,  0.931,  0.932,  0.933, 0.936,  0.938, 
0.937,  0.936,  0.938,  0.943,  0.937,  0.925,  0.925,  0.932,  0.939, 0.944,  0.944, 
0.958,  0.953,  0.931,  0.925,  0.913,  0.911,  0.909,  0.907,  0.903, 0.900,  0.902, 
0.904,  0.903,  0.900,  0.904,  0.911,  0.918,  0.929,  0.937,  0.923, 0.911,  0.898, 
0.897,  0.890,  0.883,  0.878,  0.866,  0.860,  0.850,  0.842,  0.839, 0.831,  0.832, 
0.825,  0.820,  0.818,  0.826,  0.833,  0.832,  0.828,  0.829,  0.819, 0.820,  0.823, 
0.827,  0.831,  0.836,  0.842,  0.842,  0.838,  0.835,  0.835,  0.839, 0.845,  0.849, 
0.849,  0.848,  0.851,  0.855,  0.863,  0.867,  0.876,  0.883,  0.888, 0.885,  0.890, 
0.895,  0.901,  0.897,  0.901,  0.894,  0.884,  0.870,  0.862,  0.858, 0.850,  0.845, 
0.844,  0.843,  0.845,  0.837,  0.844,  0.845,  0.845,  0.848,  0.848, 0.844,  0.846, 
0.839,  0.832,  0.839,  0.847,  0.857,  0.872,  0.881,  0.890,  0.898, 0.907,  0.933, 
0.971,  1.007,  1.053,  1.119,  1.197,  1.245,  1.226,  1.178,  1.155, 1.121,  1.068, 
1.019,  0.986,  0.965,  0.954,  0.955,  0.949,  0.949,  0.935,  0.914, 0.895,  0.881, 
0.871,  0.863,  0.857,  0.862,  0.862,  0.860,  0.857,  0.857,  0.851, 0.854,  0.857, 
0.866,  0.869,  0.862,  0.860,  0.859,  0.855,  0.857,  0.864,  0.872, 0.881,  0.894, 
0.889,  0.876,  0.866,  0.867,  0.871,  0.866,  0.870,  0.866,  0.862, 0.856,  0.850, 
0.835,  0.843,  0.837,  0.845,  0.840,  0.856,  0.866,  0.883,  0.905, 0.915,  0.922, 
0.918,  0.919,  0.910,  0.916,  0.935,  0.955,  0.957,  0.951,  0.942, 0.928,  0.929, 
0.922,  0.926,  0.919,  0.921,  0.915,  0.916,  0.918,  0.923,  0.928, 0.934,  0.947, 
0.964,  0.987,  1.016,  1.067,  1.125,  1.196,  1.330,  1.489,  1.595, 1.580,  1.468, 
1.351,  1.257,  1.175,  1.109,  1.038,  0.993,  0.947,  0.922,  0.906, 0.877,  0.867, 
0.866,  0.852,  0.846,  0.846,  0.826,  0.820,  0.817,  0.808,  0.808, 0.792,  0.792, 
0.797,  0.782,  0.781,  0.779,  0.780,  0.776,  0.783,  0.791,  0.793, 0.802,  0.818, 
0.828,  0.840,  0.866,  0.890,  0.916,  0.960,  1.016,  1.069,  1.069, 1.070,  1.061, 
1.035,  0.992,  0.959,  0.938,  0.908,  0.904,  0.898,  0.913,  0.905, 0.893,  0.901, 
0.892,  0.906,  0.916,  0.932,  0.943,  0.973,  0.978,  0.972,  0.959, 0.935,  0.954, 
0.951,  0.958,  1.009,  1.052,  1.076,  1.098,  1.097,  1.067,  1.044, 1.030,  1.008, 
1.010,  1.017,  1.030,  1.057,  1.073,  1.118,  1.180,  1.245,  1.299, 1.310,  1.355, 
1.387,  1.463,  1.550,  1.680,  1.900,  2.100,  2.187,  2.188,  2.164, 2.233,  2.371, 
2.536,  2.856,  3.141,  3.665,  3.235,  2.608,  2.279,  1.928,  1.795, 1.521,  1.477, 
1.307,  1.250,  1.183,  1.067,  0.985,  1.028,  1.045,  1.056,  1.000, 0.931,  0.886, 
0.877,  0.856,  0.856,  0.882,  0.876,  0.816,  0.842,  0.833,  0.804, 0.823,  0.865, 
0.867,  0.851,  0.937,  0.855,  0.905,  0.975,  0.945,  0.962,  0.881, 0.975,  0.910, 
0.920,  0.913,  0.827,  0.819,  0.869,  0.873,  0.849,  0.910,  0.905, 0.906,  0.920, 
1.090,  0.975,  1.019,  1.020,  1.097,  1.075,  1.149,  1.073,  1.163, 1.085,  1.059, 
1.186,  1.140,  1.147,  1.147,  1.200,  1.402,  1.505,  1.587,  1.458, 1.534,  1.490, 
1.384,  1.105,  1.240,  1.008,  1.110,  0.984,  0.900,  0.948,  0.980, 0.929,  0.845, 
0.868,  0.962,  1.033,  1.008,  1.076,  1.240,  1.048,  0.942,  1.102, 0.935,  1.051, 
1.093,  1.039,  0.927,  1.081,  0.881,  0.905,  0.868,  0.842,  0.716, 0.832,  0.972, 
0.921,  0.952 
};
double log_nu_obs = log10(nu);
if (log_nu_obs < log_nu[0])   return 1.0;
if (log_nu_obs > log_nu[924]) return 1.0;
if ((log_nu_obs>=log_nu[0])&&(log_nu_obs<=log_nu[924])) {
	int n0 = (int )((log_nu_obs-log_nu[0])/0.001);
	// interpolate over the observed range
	return nuLnu[n0] + (nuLnu[n0+1]-nuLnu[n0]) * 
							((log_nu_obs-log_nu[n0])/(log_nu[n0+1]-log_nu[n0]));
}
return 1.0;
}




// load the spectrum based on the all-quasar template of Richards et al. 2006
double l_band_GTR(double log_L_bol, double nu)
{
if (nu==0.)  return log_L_bol;
if (nu==-1.) return log_L_bol - 0.93863501;
if (nu==-2.) return log_L_bol - 0.99164163;
if (nu==-3.) return log_L_bol - 1.96992370;
if (nu==-4.) return log_L_bol - 1.90528590;

double log_nu[226]={
12.50, 12.52, 12.54, 12.56, 12.58, 12.60, 12.62, 12.64, 12.66, 12.68, 12.70, 12.72, 12.74,
12.76, 12.78, 12.80, 12.82, 12.84, 12.86, 12.88, 12.90, 12.92, 12.94, 12.96, 12.98, 13.00,
13.02, 13.04, 13.06, 13.08, 13.10, 13.12, 13.14, 13.16, 13.18, 13.20, 13.22, 13.24, 13.26,
13.28, 13.30, 13.32, 13.34, 13.36, 13.38, 13.40, 13.42, 13.44, 13.46, 13.48, 13.50, 13.52,
13.54, 13.56, 13.58, 13.60, 13.62, 13.64, 13.66, 13.68, 13.70, 13.72, 13.74, 13.76, 13.78,
13.80, 13.82, 13.84, 13.86, 13.88, 13.90, 13.92, 13.94, 13.96, 13.98, 14.00, 14.02, 14.04,
14.06, 14.08, 14.10, 14.12, 14.14, 14.16, 14.18, 14.20, 14.22, 14.24, 14.26, 14.28, 14.30,
14.32, 14.34, 14.36, 14.38, 14.40, 14.42, 14.44, 14.46, 14.48, 14.50, 14.52, 14.54, 14.56,
14.58, 14.60, 14.62, 14.64, 14.66, 14.68, 14.70, 14.72, 14.74, 14.76, 14.78, 14.80, 14.82,
14.84, 14.86, 14.88, 14.90, 14.92, 14.94, 14.96, 14.98, 15.00, 15.02, 15.04, 15.06, 15.08,
15.10, 15.12, 15.14, 15.16, 15.18, 15.20, 15.22, 15.24, 15.26, 15.28, 15.30, 15.32, 15.34,
15.36, 15.38, 15.40, 15.42, 15.44, 15.46, 15.48, 15.50, 15.52, 15.54, 15.56, 15.58, 15.60,
15.62, 15.64, 15.66, 15.68, 15.70, 15.72, 15.74, 15.76, 15.78, 15.80, 15.82, 15.84, 15.86,
15.88, 15.90, 15.92, 15.94, 15.96, 15.98, 16.00, 16.02, 16.04, 16.06, 16.08, 16.10, 16.12,
16.14, 16.16, 16.18, 16.20, 16.22, 16.24, 16.26, 16.28, 16.30, 16.32, 16.34, 16.36, 16.38,
16.40, 16.42, 16.44, 16.46, 16.48, 16.50, 16.52, 16.54, 16.56, 16.58, 16.60, 16.62, 16.64,
16.66, 16.68, 16.70, 16.72, 16.74, 16.76, 16.78, 16.80, 16.82, 16.84, 16.86, 16.88, 16.90,
16.92, 16.94, 16.96, 16.98, 17.00
};
double log_nuLnu[226]={
44.43, 44.49, 44.54, 44.59, 44.65, 44.70, 44.74, 44.79, 44.83, 44.87, 44.90, 44.93, 44.96,
44.99, 45.02, 45.04, 45.06, 45.09, 45.11, 45.13, 45.15, 45.16, 45.18, 45.20, 45.22, 45.23,
45.25, 45.26, 45.27, 45.28, 45.30, 45.31, 45.32, 45.33, 45.33, 45.34, 45.35, 45.36, 45.37,
45.37, 45.38, 45.39, 45.39, 45.40, 45.40, 45.41, 45.41, 45.42, 45.42, 45.43, 45.43, 45.44,
45.44, 45.44, 45.44, 45.45, 45.45, 45.45, 45.45, 45.45, 45.45, 45.45, 45.45, 45.45, 45.45,
45.45, 45.45, 45.45, 45.45, 45.44, 45.44, 45.44, 45.44, 45.43, 45.43, 45.42, 45.42, 45.41,
45.40, 45.39, 45.37, 45.36, 45.34, 45.32, 45.30, 45.28, 45.26, 45.24, 45.21, 45.20, 45.18,
45.17, 45.16, 45.15, 45.15, 45.16, 45.16, 45.17, 45.17, 45.19, 45.20, 45.21, 45.22, 45.23,
45.25, 45.26, 45.27, 45.29, 45.30, 45.31, 45.32, 45.33, 45.35, 45.37, 45.38, 45.40, 45.42,
45.44, 45.46, 45.48, 45.51, 45.54, 45.56, 45.59, 45.61, 45.63, 45.64, 45.65, 45.65, 45.65,
45.65, 45.66, 45.66, 45.67, 45.69, 45.70, 45.71, 45.72, 45.73, 45.74, 45.75, 45.76, 45.76,
45.76, 45.76, 45.74, 45.73, 45.70, 45.67, 45.64, 45.60, 45.57, 45.54, 45.52, 45.49, 45.48,
45.46, 45.45, 45.43, 45.41, 45.40, 45.38, 45.36, 45.34, 45.32, 45.30, 45.27, 45.25, 45.23,
45.20, 45.18, 45.16, 45.14, 45.12, 45.09, 45.07, 45.05, 45.03, 45.01, 44.98, 44.96, 44.94,
44.92, 44.90, 44.87, 44.85, 44.83, 44.81, 44.79, 44.77, 44.74, 44.72, 44.70, 44.68, 44.66,
44.63, 44.61, 44.59, 44.57, 44.55, 44.52, 44.50, 44.48, 44.46, 44.44, 44.42, 44.41, 44.39,
44.37, 44.36, 44.34, 44.33, 44.31, 44.30, 44.29, 44.29, 44.28, 44.27, 44.27, 44.27, 44.27,
44.27, 44.27, 44.26, 44.26, 44.26
};

double GTR_lbol = 46.372012;
double log_nu_obs = log10(nu);
double nuLnu_obs;

if (log_nu_obs < log_nu[0])   nuLnu_obs = log_nuLnu[0] + 2.0*(log_nu_obs - log_nu[0]);
if (log_nu_obs > log_nu[224]) nuLnu_obs = log_nuLnu[224];	// assumes Gamma=2.0
if ((log_nu_obs>=log_nu[0])&&(log_nu_obs<=log_nu[224])) {
	int n0 = (int )((log_nu_obs-log_nu[0])/0.02);
	nuLnu_obs = log_nuLnu[n0] + (log_nuLnu[n0+1]-log_nuLnu[n0]) * 
									((log_nu_obs-log_nu[n0])/(log_nu[n0+1]-log_nu[n0]));
}
return nuLnu_obs-GTR_lbol + log_L_bol;
}



// return the intrinsic band luminosity for some bolometric luminosity and frequency; 
//   nu = 0 (l_bol), -1 (B-band), -2 (15 microns), -3 (0.5-2 keV), -4 (2-10 keV), 
//   otherwise nu is the observed frequency and the return is nu*L_nu
double l_band_HRH(double log_l_bol, double nu)
{
	double x = log_l_bol - 10.;
	double lband = 0.;
	double P0,P1,P2,P3;
	if (nu==(0.))  return pow(10.,log_l_bol);
	if (nu < 0.) {
		if (nu==(-1.)) {P0=8.99833; P1=6.24800; P2=-0.370587; P3=-0.0115970;}
		if (nu==(-2.)) {P0=10.6615; P1=7.40282; P2=-0.370587; P3=-0.0115970;}
		if (nu==(-3.)) {P0=10.0287; P1=17.8653; P2=0.276804;  P3=-0.0199558;}
		if (nu==(-4.)) {P0=6.08087; P1=10.8331; P2=0.276802;  P3=-0.0199597;}
		lband = P0*pow(10.,P3*x) + P1*pow(10.,P2*x);
		return pow(10.,log_l_bol)/lband;
	}
	// if not one of the specified bands, then take advantage of the fact that our model 
	//   spectrum is not l-dependent below 500 angstroms or above 50 angstroms, so 
	//   just take the appropriate ratios to renormalize to those luminosities
	double nu_angstrom = (2.998e8)/(1.0e-10);
	double nu500 = nu_angstrom/500.;
	double nu50  = nu_angstrom/50.;
	if (nu <= nu500) {
		// just take the ratio relative to B-band
		double P[4] = {8.99833   ,   6.24800  ,  -0.370587  , -0.0115970};
		lband = P[0]*pow(10.,P[3]*x) + P[1]*pow(10.,P[2]*x);
		return return_ratio_to_b_band(nu)*pow(10.,log_l_bol)/lband;
	}
	if (nu >= nu50) {
		// just take the ratio relative to the hard X-rays
		double P[4] = {6.08087   ,   10.8331  ,   0.276802  , -0.0199597};
		lband = P[0]*pow(10.,P[3]*x) + P[1]*pow(10.,P[2]*x);
		return return_ratio_to_hard_xray(nu)*pow(10.,log_l_bol)/lband;
	}
	if ((nu>nu500)&&(nu<nu50)) {
		// interpolate between both regimes
		double P[4] = {8.99833   ,   6.24800  ,  -0.370587  , -0.0115970};
		double L500 = return_ratio_to_b_band(nu500)/(P[0]*pow(10.,P[3]*x)+P[1]*pow(10.,P[2]*x));
		double Q[4] = {6.08087   ,   10.8331  ,   0.276802  , -0.0199597};
		double L50  = return_ratio_to_hard_xray(nu50)/(Q[0]*pow(10.,Q[3]*x)+Q[1]*pow(10.,Q[2]*x));
		double L00  = log10(L500) + log10(L50/L500) * (log10(nu/nu500)/log10(nu50/nu500));	
		return pow(10.,L00)*pow(10.,log_l_bol);
	}
}

// load the x-ray template, based on the observations in text and 
//     specifically the Magdziarz & Zdziarski 1995 PEXRAV model with Gamma=1.8 
//     (Tozzi et al., George et al.), theta=2pi, solar abundances
double return_ratio_to_hard_xray(double nu)
{
double log_nu[275]={
16.00, 16.02, 16.04, 16.06, 16.08, 16.10, 16.12, 16.14, 16.16, 16.18, 16.20, 16.22, 16.24,
16.26, 16.28, 16.30, 16.32, 16.34, 16.36, 16.38, 16.40, 16.42, 16.44, 16.46, 16.48, 16.50,
16.52, 16.54, 16.56, 16.58, 16.60, 16.62, 16.64, 16.66, 16.68, 16.70, 16.72, 16.74, 16.76,
16.78, 16.80, 16.82, 16.84, 16.86, 16.88, 16.90, 16.92, 16.94, 16.96, 16.98, 17.00, 17.02,
17.04, 17.06, 17.08, 17.10, 17.12, 17.14, 17.16, 17.18, 17.20, 17.22, 17.24, 17.26, 17.28,
17.30, 17.32, 17.34, 17.36, 17.38, 17.40, 17.42, 17.44, 17.46, 17.48, 17.50, 17.52, 17.54,
17.56, 17.58, 17.60, 17.62, 17.64, 17.66, 17.68, 17.70, 17.72, 17.74, 17.76, 17.78, 17.80,
17.82, 17.84, 17.86, 17.88, 17.90, 17.92, 17.94, 17.96, 17.98, 18.00, 18.02, 18.04, 18.06,
18.08, 18.10, 18.12, 18.14, 18.16, 18.18, 18.20, 18.22, 18.24, 18.26, 18.28, 18.30, 18.32,
18.34, 18.36, 18.38, 18.40, 18.42, 18.44, 18.46, 18.48, 18.50, 18.52, 18.54, 18.56, 18.58,
18.60, 18.62, 18.64, 18.66, 18.68, 18.70, 18.72, 18.74, 18.76, 18.78, 18.80, 18.82, 18.84,
18.86, 18.88, 18.90, 18.92, 18.94, 18.96, 18.98, 19.00, 19.02, 19.04, 19.06, 19.08, 19.10,
19.12, 19.14, 19.16, 19.18, 19.20, 19.22, 19.24, 19.26, 19.28, 19.30, 19.32, 19.34, 19.36,
19.38, 19.40, 19.42, 19.44, 19.46, 19.48, 19.50, 19.52, 19.54, 19.56, 19.58, 19.60, 19.62,
19.64, 19.66, 19.68, 19.70, 19.72, 19.74, 19.76, 19.78, 19.80, 19.82, 19.84, 19.86, 19.88,
19.90, 19.92, 19.94, 19.96, 19.98, 20.00, 20.02, 20.04, 20.06, 20.08, 20.10, 20.12, 20.14,
20.16, 20.18, 20.20, 20.22, 20.24, 20.26, 20.28, 20.30, 20.32, 20.34, 20.36, 20.38, 20.40,
20.42, 20.44, 20.46, 20.48, 20.50, 20.52, 20.54, 20.56, 20.58, 20.60, 20.62, 20.64, 20.66,
20.68, 20.70, 20.72, 20.74, 20.76, 20.78, 20.80, 20.82, 20.84, 20.86, 20.88, 20.90, 20.92,
20.94, 20.96, 20.98, 21.00, 21.02, 21.04, 21.06, 21.08, 21.10, 21.12, 21.14, 21.16, 21.18,
21.20, 21.22, 21.24, 21.26, 21.28, 21.30, 21.32, 21.34, 21.36, 21.38, 21.40, 21.42, 21.44,
21.46, 21.48 
};
double log_nuLnu[275]={
-2.1132, -2.1092, -2.1052, -2.1012, -2.0972, -2.0932, -2.0892, -2.0852, -2.0812, -2.0772, -2.0732, 
-2.0692, -2.0652, -2.0612, -2.0572, -2.0532, -2.0492, -2.0452, -2.0412, -2.0372, -2.0332, -2.0292, 
-2.0252, -2.0212, -2.0172, -2.0132, -2.0092, -2.0052, -2.0012, -1.9972, -1.9932, -1.9892, -1.9852, 
-1.9812, -1.9772, -1.9732, -1.9692, -1.9652, -1.9611, -1.9571, -1.9531, -1.9491, -1.9452, -1.9412, 
-1.9372, -1.9332, -1.9292, -1.9252, -1.9212, -1.9172, -1.9132, -1.9092, -1.9052, -1.9012, -1.8971, 
-1.8931, -1.8894, -1.8854, -1.8814, -1.8774, -1.8734, -1.8694, -1.8654, -1.8614, -1.8574, -1.8534, 
-1.8495, -1.8455, -1.8415, -1.8374, -1.8334, -1.8294, -1.8253, -1.8213, -1.8172, -1.8132, -1.8091, 
-1.8050, -1.8009, -1.7968, -1.7927, -1.7885, -1.7843, -1.7803, -1.7761, -1.7718, -1.7675, -1.7631, 
-1.7587, -1.7547, -1.7502, -1.7457, -1.7410, -1.7363, -1.7314, -1.7266, -1.7215, -1.7163, -1.7109, 
-1.7053, -1.6999, -1.6940, -1.6878, -1.6814, -1.6747, -1.6677, -1.6604, -1.6527, -1.6451, -1.6369, 
-1.6283, -1.6194, -1.6342, -1.6266, -1.6188, -1.6113, -1.6036, -1.5950, -1.5861, -1.5763, -1.5655, 
-1.5530, -1.5391, -1.5253, -1.5116, -1.4968, -1.4808, -1.4635, -1.4451, -1.4259, -1.4062, -1.3867, 
-1.3676, -1.3494, -1.3323, -1.3165, -1.3021, -1.2891, -1.2775, -1.2671, -1.2580, -1.2501, -1.2432, 
-1.2372, -1.2320, -1.2277, -1.2240, -1.2210, -1.2185, -1.2166, -1.2152, -1.2143, -1.2141, -1.2144, 
-1.2150, -1.2160, -1.2172, -1.2188, -1.2207, -1.2229, -1.2253, -1.2280, -1.2310, -1.2342, -1.2376, 
-1.2413, -1.2453, -1.2495, -1.2540, -1.2588, -1.2640, -1.2694, -1.2752, -1.2814, -1.2881, -1.2952, 
-1.3028, -1.3109, -1.3196, -1.3289, -1.3389, -1.3496, -1.3609, -1.3731, -1.3861, -1.4000, -1.4148, 
-1.4302, -1.4456, -1.4602, -1.4741, -1.4877, -1.5012, -1.5147, -1.5284, -1.5424, -1.5567, -1.5714, 
-1.5866, -1.6023, -1.6186, -1.6354, -1.6529, -1.6711, -1.6901, -1.7098, -1.7304, -1.7518, -1.7742, 
-1.7975, -1.8219, -1.8474, -1.8740, -1.9018, -1.9308, -1.9612, -1.9930, -2.0263, -2.0611, -2.0975, 
-2.1356, -2.1755, -2.2173, -2.2611, -2.3070, -2.3551, -2.4055, -2.4584, -2.5138, -2.5720, -2.6331, 
-2.6971, -2.7642, -2.8347, -2.9087, -2.9864, -3.0678, -3.1533, -3.2430, -3.3371, -3.4358, -3.5394, 
-3.6481, -3.7620, -3.8814, -4.0067, -4.1382, -4.2760, -4.4205, -4.5719, -4.7306, -4.8971, -5.0716, 
-5.2545, -5.4462, -5.6471, -5.8576, -6.0783, -6.3096, -6.5520, -6.8059, -7.0719, -7.3506, -7.6428, 
-7.9489, -8.2696, -8.6055, -8.9574, -9.3260, -9.7124, -10.1171, -10.5410, -10.9850, -11.4499, -11.9371};

// want the ratio with respect to the intrinsic 2-10 keV:
double L_HX  = -1.47406;
double log_nu_obs = log10(nu);
double nuLnu_obs = 0.;

if (log_nu_obs < log_nu[0])   nuLnu_obs = log_nuLnu[0];
if (log_nu_obs > log_nu[273]) nuLnu_obs = -2.0 - pow(10.0,log_nu_obs-20.1204)*(0.43429448);
if ((log_nu_obs>=log_nu[0])&&(log_nu_obs<=log_nu[273])) {
	int n0 = (int )((log_nu_obs-log_nu[0])/0.02);
	nuLnu_obs = log_nuLnu[n0] + (log_nuLnu[n0+1]-log_nuLnu[n0]) * 
									((log_nu_obs-log_nu[n0])/(log_nu[n0+1]-log_nu[n0]));
}
return pow(10.0,nuLnu_obs-L_HX);
}



// load the optical-IR template, based on the observations in text and specifically 
//		 the Richards et al. 2006 mean blue SED 
double return_ratio_to_b_band(double nu)
{
double log_nu[226]={
12.50, 12.52, 12.54, 12.56, 12.58, 12.60, 12.62, 12.64, 12.66, 12.68, 12.70, 12.72, 12.74,
12.76, 12.78, 12.80, 12.82, 12.84, 12.86, 12.88, 12.90, 12.92, 12.94, 12.96, 12.98, 13.00,
13.02, 13.04, 13.06, 13.08, 13.10, 13.12, 13.14, 13.16, 13.18, 13.20, 13.22, 13.24, 13.26,
13.28, 13.30, 13.32, 13.34, 13.36, 13.38, 13.40, 13.42, 13.44, 13.46, 13.48, 13.50, 13.52,
13.54, 13.56, 13.58, 13.60, 13.62, 13.64, 13.66, 13.68, 13.70, 13.72, 13.74, 13.76, 13.78,
13.80, 13.82, 13.84, 13.86, 13.88, 13.90, 13.92, 13.94, 13.96, 13.98, 14.00, 14.02, 14.04,
14.06, 14.08, 14.10, 14.12, 14.14, 14.16, 14.18, 14.20, 14.22, 14.24, 14.26, 14.28, 14.30,
14.32, 14.34, 14.36, 14.38, 14.40, 14.42, 14.44, 14.46, 14.48, 14.50, 14.52, 14.54, 14.56,
14.58, 14.60, 14.62, 14.64, 14.66, 14.68, 14.70, 14.72, 14.74, 14.76, 14.78, 14.80, 14.82,
14.84, 14.86, 14.88, 14.90, 14.92, 14.94, 14.96, 14.98, 15.00, 15.02, 15.04, 15.06, 15.08,
15.10, 15.12, 15.14, 15.16, 15.18, 15.20, 15.22, 15.24, 15.26, 15.28, 15.30, 15.32, 15.34,
15.36, 15.38, 15.40, 15.42, 15.44, 15.46, 15.48, 15.50, 15.52, 15.54, 15.56, 15.58, 15.60,
15.62, 15.64, 15.66, 15.68, 15.70, 15.72, 15.74, 15.76, 15.78, 15.80, 15.82, 15.84, 15.86,
15.88, 15.90, 15.92, 15.94, 15.96, 15.98, 16.00, 16.02, 16.04, 16.06, 16.08, 16.10, 16.12,
16.14, 16.16, 16.18, 16.20, 16.22, 16.24, 16.26, 16.28, 16.30, 16.32, 16.34, 16.36, 16.38,
16.40, 16.42, 16.44, 16.46, 16.48, 16.50, 16.52, 16.54, 16.56, 16.58, 16.60, 16.62, 16.64,
16.66, 16.68, 16.70, 16.72, 16.74, 16.76, 16.78, 16.80, 16.82, 16.84, 16.86, 16.88, 16.90,
16.92, 16.94, 16.96, 16.98, 17.00
};
double log_nuLnu[226]={
44.39, 44.44, 44.50, 44.55, 44.60, 44.65, 44.70, 44.74, 44.78, 44.82, 44.86, 44.89, 44.92,
44.95, 44.97, 45.00, 45.02, 45.04, 45.06, 45.08, 45.10, 45.12, 45.14, 45.16, 45.17, 45.19,
45.20, 45.22, 45.23, 45.24, 45.25, 45.26, 45.27, 45.28, 45.29, 45.30, 45.31, 45.31, 45.32,
45.33, 45.34, 45.34, 45.35, 45.35, 45.36, 45.36, 45.37, 45.38, 45.38, 45.38, 45.39, 45.39,
45.39, 45.40, 45.40, 45.40, 45.40, 45.41, 45.41, 45.41, 45.41, 45.41, 45.41, 45.41, 45.41,
45.41, 45.40, 45.40, 45.40, 45.40, 45.40, 45.40, 45.39, 45.39, 45.38, 45.38, 45.37, 45.36,
45.35, 45.34, 45.33, 45.31, 45.30, 45.28, 45.26, 45.24, 45.22, 45.19, 45.17, 45.15, 45.13,
45.12, 45.11, 45.11, 45.11, 45.11, 45.12, 45.12, 45.13, 45.14, 45.15, 45.16, 45.18, 45.19,
45.20, 45.22, 45.23, 45.25, 45.26, 45.27, 45.28, 45.30, 45.32, 45.34, 45.36, 45.38, 45.40,
45.42, 45.44, 45.47, 45.49, 45.52, 45.55, 45.58, 45.60, 45.62, 45.64, 45.65, 45.65, 45.66,
45.66, 45.67, 45.68, 45.69, 45.71, 45.72, 45.74, 45.75, 45.77, 45.78, 45.79, 45.80, 45.80,
45.80, 45.80, 45.79, 45.77, 45.74, 45.71, 45.68, 45.64, 45.60, 45.57, 45.55, 45.52, 45.51,
45.49, 45.47, 45.46, 45.44, 45.42, 45.40, 45.38, 45.36, 45.34, 45.32, 45.29, 45.27, 45.25,
45.22, 45.20, 45.18, 45.15, 45.13, 45.11, 45.08, 45.06, 45.04, 45.02, 44.99, 44.97, 44.95,
44.92, 44.90, 44.88, 44.86, 44.83, 44.81, 44.79, 44.77, 44.74, 44.72, 44.70, 44.67, 44.65,
44.63, 44.61, 44.58, 44.56, 44.54, 44.52, 44.50, 44.47, 44.45, 44.43, 44.41, 44.40, 44.38,
44.36, 44.35, 44.33, 44.32, 44.30, 44.29, 44.28, 44.27, 44.27, 44.26, 44.26, 44.26, 44.26,
44.26, 44.25, 44.25, 44.25, 44.25
};

// want the ratio with respect to the intrinsic B-band:
double nu_BB = 14.833657;
double L_BB  = 45.413656;
double log_nu_obs = log10(nu);
double nuLnu_obs = 0.;

if (log_nu_obs < log_nu[0])   nuLnu_obs = log_nuLnu[0] + 2.0*(log_nu_obs - log_nu[0]);
if (log_nu_obs > log_nu[224]) nuLnu_obs = log_nuLnu[224];
	// assumes Gamma=2.0; the calling code will actually use X-ray template for this case
if ((log_nu_obs>=log_nu[0])&&(log_nu_obs<=log_nu[224])) {
	int n0 = (int )((log_nu_obs-log_nu[0])/0.02);
	nuLnu_obs = log_nuLnu[n0] + (log_nuLnu[n0+1]-log_nuLnu[n0]) * 
									((log_nu_obs-log_nu[n0])/(log_nu[n0+1]-log_nu[n0]));
}
return pow(10.0,nuLnu_obs-L_BB);
}