Positrons implemented.

model/iharm/model.c

@@ -47,6 +47,9 @@ static double Mdot_dump;
 static double MdotEdd_dump;
 static double Ladv_dump;
 
+//ARR: positron to electron ratio
+double positronRatio = 0.0;
+
 // MAYBES
 //static double t0;
 
@@ -99,6 +102,9 @@ void try_set_model_parameter(const char *word, const char *value)
   set_by_word_val(word, value, "dump", (void *)fnam, TYPE_STR);
   set_by_word_val(word, value, "counterjet", &counterjet, TYPE_INT);
 
+  //ARR:  positron to electron ratio
+  set_by_word_val(word, value, "positronRatio", &positronRatio, TYPE_DBL);
+
   set_by_word_val(word, value, "tp_over_te", &tp_over_te, TYPE_DBL);
   set_by_word_val(word, value, "trat_small", &trat_small, TYPE_DBL);
   set_by_word_val(word, value, "trat_large", &trat_large, TYPE_DBL);

model/iharm/model_params.h

@@ -18,6 +18,7 @@
 extern double DTd;
 extern double sigma_cut;
 extern int counterjet;
+extern double positronRatio;
 extern double rmax_geo;
 
 #endif // MODEL_PARAMS_H

src/ipolarray.c

@@ -125,6 +125,19 @@ void evolve_N(double Xi[NDIM], double Kconi[NDIM],
       }
     }
 
+    //ARR:  Modify according to positron ratio
+    jI *= (1 + positronRatio);
+    jQ *= (1 + positronRatio);
+    jU *= (1 + positronRatio);
+    jV *= (1 - positronRatio);
+    aI *= (1 + positronRatio);
+    aQ *= (1 + positronRatio);
+    aU *= (1 + positronRatio);
+    aV *= (1 - positronRatio);
+    rQ *= (1 + positronRatio);
+    rU *= (1 + positronRatio);
+    rV *= (1 - positronRatio);
+
     /* make plasma tetrad */
     B = get_model_b(Xf); /* field in G */
     if (B < 0.) {

src/par.c

@@ -67,6 +67,7 @@ void load_par_from_argv(int argc, char *argv[], Params *params) {
   params->trace_stride = 1;
   params->trace_i = -1;
   params->trace_j = -1;
+
   // This is what the par infra does.
   // I'm not sure there's still any advantage to "const" if we do this,
   // but hey, no warnings

src/par.h

@@ -24,6 +24,8 @@ typedef struct params_t {
   int quench_output;    // Quench output, i.e. "quench" argument
   int only_unpolarized; // Unpolarized transport only
 
+  double positronRatio; //Ratio of positrons to electrons.
+
   const char dump[STRLEN];
   const char outf[STRLEN];