Ricardo

.gitignore

@@ -0,0 +1,12 @@
+model.c
+model.h
+*.p
+*.png
+*.dat
+*.pyc
+grmonty
+interpolate.py
+iterate.py
+iterror.py
+plot2.py
+plot3.py

build.py

@@ -136,5 +136,16 @@ def print_config(key, var):
 
   sys.exit()
 
+if len(sys.argv) != 2:
+  print 'ERROR: Format is'
+  print '  python build.py [model]'
+  sys.exit()
+if (not os.path.isfile('model/' + sys.argv[1] + '.c')):
+  print 'ERROR Model %s does not exist' % sys.argv[1]
+  sys.exit()
+
+import shutil
+shutil.copyfile('model/' + sys.argv[1] + '.c', 'model.c')
+shutil.copyfile('model/' + sys.argv[1] + '.h', 'model.h')
 build()
 

compton.c

@@ -17,7 +17,7 @@ Canfield, Howard, and Liang, 1987, ApJ 323, 565.
 
 */
 
-#define OLD_E_SAMP (0)
+#define OLD_E_SAMP (1)
 
 void sample_scattered_photon(double k[4], double p[4], double kp[4])
 {
@@ -43,9 +43,17 @@ void sample_scattered_photon(double k[4], double p[4], double kp[4])
 
 	/* unit vector 1 for scattering coordinate system is
 	   oriented along initial photon wavevector */
-	v0x = ke[1] / ke[0];
-	v0y = ke[2] / ke[0];
-	v0z = ke[3] / ke[0];
+	//v0x = ke[1] / ke[0];
+	//v0y = ke[2] / ke[0];
+	//v0z = ke[3] / ke[0];
+
+  // Explicitly compute kemag instead of using ke[0] to ensure that photon
+  // is created normalized and doesn't inherit light cone errors from the
+  // original superphoton
+  double kemag = sqrt(ke[1]*ke[1] + ke[2]*ke[2] + ke[3]*ke[3]);
+  v0x = ke[1]/kemag;
+  v0y = ke[2]/kemag;
+  v0z = ke[3]/kemag;
 
 	/* randomly pick zero-angle for scattering coordinate system.
 	   There's undoubtedly a better way to do this. */

decs.h

@@ -1,6 +1,7 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <stdarg.h>
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -14,30 +15,32 @@
 #include <omp.h>
 #include <time.h>
 #include "constants.h"
+#include "model.h"
 
-#define NDIM  4
-#define NPRIM 8
+//#define NDIM  4
+//#define NPRIM 8
 
 /* Range of initial superphoton frequencies */
-#define NUMIN 1.e9
-#define NUMAX 1.e20
+/*#define NUMIN 1.e9
+#define NUMAX 1.e28
 #define LNUMIN log(NUMIN)
 #define LNUMAX log(NUMAX)
-#define DLNU ((LNUMAX-LNUMIN)/N_ESAMP)
+#define DLNU ((LNUMAX-LNUMIN)/N_ESAMP)*/
 
-#define THETAE_MAX  1000.
+/*#define THETAE_MAX  1000.
 #define THETAE_MIN  0.3
 //#define TP_OVER_TE  (3.)
 #define WEIGHT_MIN  (1.e28)
 
 #define SYNCHROTRON (1)
-#define BREMSSTRAHLUNG (1)
-#define COMPTON (1)
+#define BREMSSTRAHLUNG (0)
+#define COMPTON (0)
 #define KAPPA (5.)
-#define DIST_KAPPA (1)
+<<<<<<< HEAD
+#define DIST_KAPPA (0)*/
 
 /* mnemonics for primitive vars; conserved vars */
-#define KRHO     0
+/*#define KRHO     0
 #define UU      1
 #define U1      2
 #define U2      3
@@ -46,13 +49,13 @@
 #define B2      6
 #define B3      7
 #define KEL     8
-#define KTOT    9
+#define KTOT    9*/
 
 /* numerical convenience */
-#define SMALL 1.e-40
+//#define SMALL 1.e-40
 
 /* physical parameters */
-#define MMW 0.5   /* mean molecular weight, in units of mp */
+//#define MMW 0.5   /* mean molecular weight, in units of mp */
 
 /** data structures **/
 struct of_photon {
@@ -96,9 +99,9 @@ struct of_spectrum {
   double E0;
 };
 
-#define N_ESAMP   200
-#define N_EBINS   200
-#define N_THBINS  6
+//#define N_ESAMP   200
+//#define N_EBINS   200
+//#define N_THBINS  6
 
 extern struct of_spectrum spect[N_THBINS][N_EBINS];
 
@@ -142,6 +145,11 @@ extern double gam;
 extern double dMsim;
 
 extern double M_unit;
+
+extern double Mdot;
+extern double mdot;
+extern double MdotEdd;
+
 extern double L_unit;
 extern double T_unit;
 extern double RHO_unit;
@@ -172,7 +180,7 @@ void scatter_super_photon(struct of_photon *ph, struct of_photon *php,
   double Ne, double Thetae, double B, double Ucon[NDIM], double Bcon[NDIM],
   double Gcov[NDIM][NDIM]);
 
-void report_bad_input();
+void report_bad_input(int argc);
 
 /* OpenMP specific functions */
 void omp_reduce_spect(void);
@@ -274,7 +282,8 @@ void init_geometry(void);
 void init_data(int argc, char *argv[]);
 void init_nint_table(void);
 void init_storage(void);
-double dOmega_func(double Xi[NDIM], double Xf[NDIM]);
+//double dOmega_func(double Xi[NDIM], double Xf[NDIM]);
+double dOmega_func(int j);
 
 double linear_interp_weight(double nu);