Merge branch 'dev' of https://github.com/agnwinds/python into dev

.gitignore

@@ -8,6 +8,9 @@ bin/modify_wind*
 
 # Run files
 *.spec*
+*.pf.old
+*.out.pf
+mod.pf
 
 # Object files
 *.o

source/atomicdata.c

@@ -126,7 +126,7 @@ get_atomic_data (masterfile)
   double auger_branches[NAUGER_ELECTRONS];      /* array to hold branching ratios for number of Auger electrons */
   double Avalue_auger;          /* Auger A value for macro-atom data */
   int ne_records;               /* number of Auger electron entries to read in for each Auger record (normally 4) */
-  int target_istate;
+  int target_istate, n_augertarget;
   char choice;
   int lineno;                   /* the line number in the file beginning with 1 */
   int simple_line_ignore[NIONS], cstren_no_line;
@@ -1429,16 +1429,16 @@ described as macro-levels. */
                    target ion */
                 target_istate = istate + 1 + m;
 
-                n = 0;
-                while ((xconfig[n].z != z || xconfig[n].istate != target_istate || xconfig[n].ilv != 1) && n < nlevels)
-                  n++;
-                if (n == nlevels)
+                n_augertarget = 0;
+                while ((xconfig[n_augertarget].z != z || xconfig[n_augertarget].istate != target_istate || xconfig[n_augertarget].ilv != 1) && n_augertarget < nlevels)
+                  n_augertarget++;
+                if (n_augertarget == nlevels)
                 {
                   Error ("Get_atomic_data: No target ion configuration found to match Auger macro record %d\n", lineno);
                   exit (0);
                 }
 
-                auger_macro[nauger_macro].nconfig_target[m] = n;
+                auger_macro[nauger_macro].nconfig_target[m] = n_augertarget;
               }
             }
 

source/disk.c

@@ -405,7 +405,7 @@ ds_to_disk (p, allow_negative, hit)
 
   /* Begin the calculation for the case of a FLAT disk */
 
-  s_diskplane = ds_to_plane (&diskplane, p);
+  s_diskplane = ds_to_plane (&diskplane, p, FALSE);
   stuff_phot (p, &phit);
   move_phot (&phit, s_diskplane);
   r_diskplane = sqrt (phit.x[0] * phit.x[0] + phit.x[1] * phit.x[1]);
@@ -480,11 +480,11 @@ ds_to_disk (p, allow_negative, hit)
 
   r_phot = sqrt (p->x[0] * p->x[0] + p->x[1] * p->x[1]);
 
-  s_top = ds_to_plane (&disktop, p);
+  s_top = ds_to_plane (&disktop, p, FALSE);
   stuff_phot (p, &phit);
   move_phot (&phit, s_top);
 
-  s_bot = ds_to_plane (&diskbottom, p);
+  s_bot = ds_to_plane (&diskbottom, p, FALSE);
   stuff_phot (p, &phit);
   move_phot (&phit, s_bot);
 
@@ -969,14 +969,15 @@ disk_height (double s, void *params)
  *
  **********************************************************/
 
-double disk_colour_correction(double t)
+double
+disk_colour_correction (double t)
 {
   double fcol;
   fcol = 1.0;
 
   if (geo.colour_correction == FCOL_OFF)
   {
-    Error("trying to apply colour correction when it is turned off! Exiting.\n");
+    Error ("trying to apply colour correction when it is turned off! Exiting.\n");
     Exit (0);
   }
   /* apply the Done et al. 2012 colour correction */
@@ -988,15 +989,15 @@ double disk_colour_correction(double t)
     {
       fcol = 2.7;
     }
-    else 
+    else
     {
-      fcol = pow(t / 3.0e4, 0.82);
+      fcol = pow (t / 3.0e4, 0.82);
     }
   }
-  else 
+  else
   {
     Error ("Did not understand colour correction mode. Setting fcol to 1");
-    fcol = 1.0; 
+    fcol = 1.0;
   }
 
   return (fcol);

source/estimators_simple.c

@@ -614,7 +614,7 @@ normalise_simple_estimators (xplasma)
   electron_density_obs = xplasma->ne * wmain[nwind].xgamma_cen; // Mihalas & Mihalas p146
   volume_obs = wmain[nwind].vol / wmain[nwind].xgamma_cen;
 
-  for (i = 0; i < 4; i++)
+  for (i = 0; i < NFORCE_DIRECTIONS; i++)
   {
     xplasma->rad_force_es[i] *= (volume_obs * electron_density_obs) / (volume_obs * VLIGHT);
     xplasma->F_vis[i] /= volume_obs;

source/matom.c

@@ -369,7 +369,7 @@ matom (p, nres, escape)
     }
     else if (n < (nbbd + nbfd + nbbu + nbfu + nauger))
     {                           /* auger ionization jump */
-      uplvl = auger_macro[iauger].nconfig_target[n - nbbd - nbfd - nbbu - nauger];
+      uplvl = auger_macro[iauger].nconfig_target[n - nbbd - nbfd - nbbu - nbfu];
 //OLD      icheck = 4;
     }
     else

source/phot_util.c

@@ -453,6 +453,8 @@ quadratic (a, b, c, r)
  *
  * @param [in] struct plane *  pl   A structure which descibes a plane in terms of a position and vector normal to the surface
  * @param [in] struct photon *  p   A photon
+ * @param [in] int force_postive_z  Whether to force the photon to be moving towards the positve z quadrant (appropriate for planes that 
+               are symmetric about z=0, but only defined once, such as the windplanes.)
  * @return     The distance to the plane, or VERY_BIG if the photon does not hit the plane (in the positive direction)
  *
  * @details
@@ -469,25 +471,29 @@ quadratic (a, b, c, r)
  **********************************************************/
 
 double
-ds_to_plane (pl, p)
+ds_to_plane (pl, p, force_positive_z)
      struct plane *pl;
      struct photon *p;
+     int force_positive_z;
 {
   double denom, diff[3], numer;
-  double dot ();
+  struct photon ptest;
 
+  stuff_phot (p, &ptest);
+  if (ptest.lmn[2] < 0 && (force_positive_z == TRUE))
+  {
+    ptest.x[2] = -ptest.x[2];   /* force the photon to be in the positive x,z quadrant */
+    ptest.lmn[2] = -ptest.lmn[2];       /* force the photon to moving in the positive z direction */
+  }
 
-  if ((denom = dot (p->lmn, pl->lmn)) == 0)
+  if ((denom = dot (ptest.lmn, pl->lmn)) == 0)
     return (VERY_BIG);
 
-  vsub (pl->x, p->x, diff);
+  vsub (pl->x, ptest.x, diff);
 
   numer = dot (diff, pl->lmn);
 
   return (numer / denom);
-
-
-
 }
 
 

source/photon2d.c

@@ -142,7 +142,7 @@ int
 translate_in_space (pp)
      PhotPtr pp;
 {
-  double ds, delta, s, smax;
+  double ds, delta, s, smax, prhosq;
   int ndom, ndom_next;
   struct photon ptest;
 
@@ -160,16 +160,23 @@ translate_in_space (pp)
     ds += DFUDGE;               //Fix for Bug #592 - we need to keep track of the little DFUDGE we moved the test photon
 
 
-
     /* Note there is a possibility that we reach the other side 
      * of the grid without actually encountering a
      * wind cell
      */
 
+    prhosq = (pp->x[0] * pp->x[0]) + (pp->x[1] * pp->x[1]);
 
-    if (where_in_wind (ptest.x, &ndom_next) < 0)
+    if ((prhosq > (zdom[ndom].wind_rhomin_at_disk * zdom[ndom].wind_rhomin_at_disk)) &&
+        (prhosq < (zdom[ndom].wind_rhomax_at_disk * zdom[ndom].wind_rhomax_at_disk)))
     {
+      stuff_phot (pp, &ptest);
+      ds = 0.0;
+    }
+
 
+    if (where_in_wind (ptest.x, &ndom_next) < 0)
+    {
       smax = ds_to_wind (&ptest, &ndom_next);   // This is the maximum distance can go in this domain
       s = 0;
       while (s < smax && where_in_wind (ptest.x, &ndom_next) < 0)
@@ -197,7 +204,6 @@ translate_in_space (pp)
 
   move_phot (pp, ds + DFUDGE);
 
-
   return (pp->istat);
 }
 
@@ -311,7 +317,7 @@ ds_to_wind (pp, ndom_current)
 
     else if (zdom[ndom].wind_type == CORONA || (zdom[ndom].wind_type == IMPORT && zdom[ndom].coord_type == CYLIND))
     {
-      x = ds_to_plane (&zdom[ndom].windplane[0], &ptest);
+      x = ds_to_plane (&zdom[ndom].windplane[0], &ptest, TRUE);
       if (x > 0 && x < ds)
       {
         stuff_phot (pp, &qtest);
@@ -324,7 +330,7 @@ ds_to_wind (pp, ndom_current)
           xxxbound = BOUND_ZMIN;
         }
       }
-      x = ds_to_plane (&zdom[ndom].windplane[1], &ptest);
+      x = ds_to_plane (&zdom[ndom].windplane[1], &ptest, TRUE);
       if (x > 0 && x < ds)
       {
         stuff_phot (pp, &qtest);

source/python.h

@@ -991,15 +991,6 @@ typedef struct plasma
 
   double cell_spec_flux[NBINS_IN_CELL_SPEC];    /**< The array where the cell spectra are accumulated. */
 
-  /* directional fluxes (in observer frame) in 3 wavebands. - last element contains the  magnitude of flux) */
-  double F_vis[4];
-  double F_UV[4];
-  double F_Xray[4];
-
-  double F_vis_persistent[4];
-  double F_UV_persistent[4];
-  double F_Xray_persistent[4];
-
 #define NFLUX_ANGLES 36 /**< The number of bins into which the directional flux is calculated */
 
 
@@ -1057,10 +1048,18 @@ typedef struct plasma
                                    compute compton cooling-  only needs computing once per cycle
                                  */
 
-#define NUM_RAD_FORCE_DIRECTIONS 3
+#define N_DMO_DT_DIRECTIONS 3
 #define NFORCE_DIRECTIONS 4
+  /* directional fluxes (in observer frame) in 3 wavebands. - last element contains the  magnitude of flux) */
+  double F_vis[NFORCE_DIRECTIONS];
+  double F_UV[NFORCE_DIRECTIONS];
+  double F_Xray[NFORCE_DIRECTIONS];
+
+  double F_vis_persistent[NFORCE_DIRECTIONS];
+  double F_UV_persistent[NFORCE_DIRECTIONS];
+  double F_Xray_persistent[NFORCE_DIRECTIONS];
 
-  double dmo_dt[NUM_RAD_FORCE_DIRECTIONS];             /**< Radiative force of wind */
+  double dmo_dt[N_DMO_DT_DIRECTIONS];             /**< Radiative force of wind */
   double rad_force_es[NFORCE_DIRECTIONS];       /**< Radiative force of wind - 4th element is sum of magnitudes */
   double rad_force_ff[NFORCE_DIRECTIONS];       /**< Radiative force of wind - 4th element is sum of magnitudes */
   double rad_force_bf[NFORCE_DIRECTIONS];       /**< Radiative force of wind - 4th element is sum of magnitudes */

source/reverb.c

@@ -46,7 +46,7 @@ delay_to_observer (PhotPtr pp)
   observer.lmn[0] = pp->lmn[0];
   observer.lmn[1] = pp->lmn[1];
   observer.lmn[2] = pp->lmn[2];
-  return (pp->path + ds_to_plane (&observer, pp));
+  return (pp->path + ds_to_plane (&observer, pp, FALSE));
 }
 
 /**********************************************************/

source/roche.c

@@ -286,7 +286,6 @@ pillbox (p, smin, smax)
   double root[2], ss[2];
   int i, n;
   struct photon pp;
-  double ds_to_plane ();
 
   n = 0;
 
@@ -343,7 +342,7 @@ root is really only possible if the photon is already in the pillbox  */
 planes of the two end caps and then checking the cylindrical radius. Negative
 distances are valid only if the photon is in the pillbox already */
 
-  x1 = ds_to_plane (&plane_m2_near, p); /* Calculate the distance to the L1 plane */
+  x1 = ds_to_plane (&plane_m2_near, p, FALSE);  /* Calculate the distance to the L1 plane */
   stuff_phot (p, &pp);
   move_phot (&pp, x1);          /* So pp is now located at the l1 plane */
   if ((pp.x[1] * pp.x[1] + pp.x[2] * pp.x[2]) <= geo.r2_width * geo.r2_width)
@@ -352,7 +351,7 @@ distances are valid only if the photon is in the pillbox already */
     n++;
   }
 
-  x2 = ds_to_plane (&plane_m2_far, p);  /* Calculate the distance to the plane behind the the star.  */
+  x2 = ds_to_plane (&plane_m2_far, p, FALSE);   /* Calculate the distance to the plane behind the the star.  */
   stuff_phot (p, &pp);
   move_phot (&pp, x2);          /* So pp is now located at the r2_far plane */
   if ((pp.x[1] * pp.x[1] + pp.x[2] * pp.x[2]) <= geo.r2_width * geo.r2_width)

source/templates.h

@@ -392,7 +392,7 @@ double ds_to_cone(ConePtr cc, struct photon *p);
 double ds_to_sphere(double r, struct photon *p);
 double ds_to_sphere2(double x[], double r, struct photon *p);
 int quadratic(double a, double b, double c, double r[]);
-double ds_to_plane(struct plane *pl, struct photon *p);
+double ds_to_plane(struct plane *pl, struct photon *p, int force_positive_z);
 double ds_to_closest_approach(double x[], struct photon *p, double *impact_parameter);
 double ds_to_cylinder(double rho, struct photon *p);
 /* photon2d.c */

source/tests/Makefile

@@ -47,7 +47,8 @@ TEST_SOURCES = \
 	tests/test_matrix.c \
 	tests/test_compton.c \
 	tests/test_define_wind.c \
-	tests/test_run_mode.c
+	tests/test_run_mode.c \
+	tests/test_translate.c
 
 # Using absolute paths
 PYTHON_SOURCES := $(patsubst %,$(PYTHON)/source/%, $(PYTHON_SOURCES))

source/tests/tests/suites.h

@@ -23,4 +23,7 @@ void create_matrix_test_suite (void);
 /* test_run_mode.c */
 void create_run_mode_test_suite (void);
 
+/* test_matrix.c */
+void create_translate_test_suite (void);
+
 #endif