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

docs/sphinx/source/input/parameters/Disk/Disk.colour_correction.rst

@@ -0,0 +1,22 @@
+Disk.colour_correction
+==========================
+Type of colour correction to use
+
+Type
+  Enumerator
+
+Values
+  Done12
+    Temperature dependent form of colour correction from Done 2012 (see :doc:`Disk <../radiation/disk>`)
+
+
+File
+  `setup_disk.c <https://github.com/agnwinds/python/blob/master/source/setup_disk.c>`_
+
+
+Parent(s)
+  * :ref:`Disk.radiation`: ``True``
+
+  * :ref:`Disk.type`: ``flat``, ``vertically.extended``
+
+  * :ref:`Disk.rad_type_to_make_wind`: ``bb``, ``models``, ``mod_bb``

docs/sphinx/source/input/parameters/Disk/Disk.rad_type_in_final_spectrum.rst

@@ -7,13 +7,13 @@ Type
 
 Values
   bb
-    Multi-line description, must keep indentation.
+    Blackbody from each annulus
 
   models
-    Multi-line description, must keep indentation.
+    Use model files such as stellar atmosphers
 
-  uniform
-    Multi-line description, must keep indentation.
+  mod_bb
+    modified blackbody (colour correction)
 
 
 File

docs/sphinx/source/input/parameters/Disk/Disk.rad_type_to_make_wind.rst

@@ -7,10 +7,13 @@ Type
 
 Values
   bb
-    Multi-line description, must keep indentation.
+    Blackbody from each annulus
 
   models
-    Multi-line description, must keep indentation.
+    Use model files such as stellar atmosphers
+
+  mod_bb
+  	modified blackbody (colour correction)
 
 
 File
@@ -26,3 +29,5 @@ Parent(s)
 Child(ren)
   * :ref:`Input_spectra.model_file`
 
+  * :ref:`Disk.colour_correction`
+

docs/sphinx/source/radiation/disk.rst

@@ -1,24 +1,33 @@
 The Disk
 ########
 
-The disk is normally treated as infinitely thin and defined by an inner boundary and an outer boundary.  It assumed to be in  Keplerian  rotation about 
-the central object in the system.   The temperature distribution of the disk
-is normally assumed to be that of a standard Shakura-Sunyaev disk, with a hard
-boundary at its inner edge.   Options are provided for reading in a non-standard
-temperature distribution.
+The disk is normally treated as infinitely thin and defined by an inner boundary and an outer boundary.  It assumed to be in  Keplerian  rotation about the central object in the system.   The temperature distribution of the disk is normally assumed to be that of a standard Shakura-Sunyaev disk, with a hard boundary at its inner edge. Options are provided for reading in a non-standard temperature distribution.
 
 An option is provide for a vertically extended disk, whose thickness increases
-as with distance from the central object object.   
+as with distance from the  central object object.   
 
 The parameters involved in describing a flat disk are::
 
     Disk.type(none,flat,vertically.extended)                 flat
     Disk.radiation(yes,no)                          yes
-    Disk.rad_type_to_make_wind(bb,models)                   bb
+    Disk.rad_type_to_make_wind(bb,models,mod_bb)                   bb
     Disk.temperature.profile(standard,readin)             standard
     Disk.mdot(msol/yr)                         5
     Disk.radmax(cm)                            1e17
 
+Colour Correction (mod_bb)
+=============================
+
+A simple form of the disc colour correction is available in the code, accessible via the 
+:ref:`Disk.rad_type_to_make_wind(bb,models,mod_bb)` keyword. The colour correction factor, :math:`f_{\rm col}`, is defined such that 
+
+.. math::
+    B_\nu (\nu, T) \to f_{\rm col}^{-4} B_\nu (\nu, f_{\rm col} T). 
+
+This correction is designed to approximate the effect of radiative transfer in the disc atmosphere. We adopt the form given by `Done et al. 2012 <https://academic.oup.com/mnras/article/420/3/1848/977649>`_, in which :math:`f_{\rm col}=1` for :math:`T<3\times10^4~{\rm K}`, and for :math:`T>3\times10^4~{\rm K}`
+
+.. math::
+    f_{\rm col}(T)=\left(\frac{T}{3\times10^4~{\rm K}} \right)^{0.82}.
 
 
 Vertically Extended disk (Details)
@@ -52,16 +61,9 @@ not take into account the fact that the disk area of a vertically extended disk
 Non-Standard Temperature Profile
 ================================================
 
-If desired the user can read the temperature profile for the disk from a file. Each
-line in the file should consist of a radius and a temperature (and optionally a value of log g)
-separated by whitespace (in the
-first two columns)  The values are assumed to be entered in a logical order, that is in 
-ascending values of radius. Lines, such as comments or header names of an astropy table, will be ignored.
+If desired the user can read the temperature profile for the disk from a file. Each line in the file should consist of a radius (in cm) and a temperature (in K), and optionally a value of :math:`\log g`. The values separated by whitespace (in the first two columns). The values are assumed to be entered in a logical order, that is in  ascending values of radius. Lines such as comments or header names of an astropy table, will be ignored.
 
-The log g value is not required to generate BB spectra, but is required if the spectrum from the disk is to be generated from a 
-two-dimensional grid of models, usually a set of spectra generated to represent the spectra from a set of stellar 
-atmospheres calculations.  
+The :math:`\log g` value is not required to generate BB spectra, but is required if the spectrum from the disk is to be generated from a two-dimensional grid of models, usually a set of spectra generated to represent the spectra from a set of stellar atmospheres calculations.  
 
-With this option, the radius of the disk will be set to the maximum radius (the last value of r) in 
-the file.  
+With this option, the radius of the disk will be set to the maximum radius (the last value of r) in  the file.  
 

source/disk.c

@@ -957,3 +957,47 @@ disk_height (double s, void *params)
   return (z1);
 
 }
+
+
+/**********************************************************/
+/**
+ * @brief      Calculate the colour correction 
+ *
+ * @param [in] double  t temperature in the disk
+
+ * @return     
+ *
+ **********************************************************/
+
+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");
+    Exit (0);
+  }
+  /* apply the Done et al. 2012 colour correction */
+  else if (geo.colour_correction == FCOL_DONE)
+  {
+    if (t < 3.0e4)
+      fcol = 1.0;
+    else if (t > 1e5)
+    {
+      fcol = 2.7;
+    }
+    else 
+    {
+      fcol = pow(t / 3.0e4, 0.82);
+    }
+  }
+  else 
+  {
+    Error ("Did not understand colour correction mode. Setting fcol to 1");
+    fcol = 1.0; 
+  }
+
+  return (fcol);
+}

source/disk_init.c

@@ -74,7 +74,7 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot)
   int nrings, i;
   int spectype;
   double emit;
-  double factor;
+  double factor, fcol;
 
 
 
@@ -159,10 +159,14 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot)
     {
       emit = emittance_continuum (spectype, freqmin, freqmax, t, log_g);
     }
+    else if (spectype == SPECTYPE_BB_FCOL)
+    {
+      fcol = disk_colour_correction(t);
+      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0);
+    }
     else
     {
       emit = emittance_bb (freqmin, freqmax, t);
-
     }
     (*ftot) += emit * (2. * r + dr) * dr * factor;
   }
@@ -215,6 +219,11 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot)
     {
       emit = emittance_continuum (spectype, freqmin, freqmax, t, log_g);
     }
+    else if (spectype == SPECTYPE_BB_FCOL)
+    {
+      fcol = disk_colour_correction(t);
+      emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0);
+    }
     else
     {
       emit = emittance_bb (freqmin, freqmax, t);

source/disk_photon_gen.c

@@ -60,6 +60,7 @@ photo_gen_disk (p, weight, f1, f2, spectype, istart, nphot)
   double r, z, theta, phi;
   int nring = 0;
   double north[3];
+  double fcol;
 
   if ((iend = istart + nphot) > NPHOT)
   {
@@ -186,6 +187,11 @@ photo_gen_disk (p, weight, f1, f2, spectype, istart, nphot)
     {
       p[i].freq = planck (disk.t[nring], freqmin, freqmax);
     }
+    else if (spectype == SPECTYPE_BB_FCOL)
+    {
+      fcol = disk_colour_correction(disk.t[nring]);
+      p[i].freq = planck (fcol * disk.t[nring], freqmin, freqmax);
+    }
     else if (spectype == SPECTYPE_UNIFORM)
     {
       p[i].freq = random_number (freqmin, freqmax);

source/python.c

@@ -447,7 +447,7 @@ main (argc, argv)
     if (geo.disk_radiation)
     {
       geo.disk_spectype = geo.disk_ion_spectype;
-      get_spectype (geo.disk_radiation, "Disk.rad_type_in_final_spectrum(bb,models,uniform,mono)", &geo.disk_spectype);
+      get_spectype (geo.disk_radiation, "Disk.rad_type_in_final_spectrum(bb,models,uniform,mono,mod_bb)", &geo.disk_spectype);
     }
 
     if (geo.bl_radiation)

source/python.h

@@ -153,9 +153,13 @@ extern int NWAVE_NOW;         /**< Either NWAVE_IONIZ or NWAVE_EXTRACT depending
 #define SPECTYPE_CL_TAB  -5   // This is to emulate cloudy
 #define SPECTYPE_BREM    -6
 #define SPECTYPE_MONO    -7
-#define SPECTYPE_NONE	 -3
+#define SPECTYPE_NONE	   -3
+#define SPECTYPE_BB_FCOL -8     // colour corrected blackbody 
 #define SPECTYPE_MODEL	 -99    // This is just used briefly, before a model number is assigned
 
+/* definitions of types of colour correction */
+#define FCOL_OFF  0
+#define FCOL_DONE 1 
 
 /* Number of model_lists that one can have, should be the same as NCOMPS in models.h */
 #define NCOMPS 	10
@@ -580,6 +584,8 @@ struct geometry
   int bl_ion_spectype, bl_spectype;     /**<  Same as above but for the boundary layer */
   int agn_ion_spectype, agn_spectype;   /**<  Same as above but for the AGN */
 
+  int colour_correction;    /** colour correction mode */
+
   /* Searchlight mode is very experimental.  See notes in diag.c */
   double searchlight_x[3], searchlight_lmn[3];/**< location and direction of all photons in the spectral
                                                 * cycles when searchlight mode (an advanced option) is

source/rdpar_init.c

@@ -77,10 +77,10 @@ init_choices ()
 {
   /* Initialize the structure that contains all of the types of possible radiation types */
 
-  char *xchoices[] = { "bb", "uniform", "power", "cloudy", "brems", "none", "models", "mono" };
+  char *xchoices[] = { "bb", "uniform", "power", "cloudy", "brems", "none", "models",  "mod_bb", "mono" };
   int xvals[] =
-    { SPECTYPE_BB, SPECTYPE_UNIFORM, SPECTYPE_POW, SPECTYPE_CL_TAB, SPECTYPE_BREM, SPECTYPE_NONE, SPECTYPE_MODEL, SPECTYPE_MONO };
-  int num_choices = 8;
+    { SPECTYPE_BB, SPECTYPE_UNIFORM, SPECTYPE_POW, SPECTYPE_CL_TAB, SPECTYPE_BREM, SPECTYPE_NONE, SPECTYPE_MODEL, SPECTYPE_BB_FCOL, SPECTYPE_MONO};
+  int num_choices = 9;
 
   if (xinit_choices)
     return (0);

source/setup.c

@@ -115,6 +115,7 @@ init_geo ()
   geo.disk_type = DISK_FLAT;    /*1 implies existence of a disk for purposes of absorption */
   geo.disk_rad_max = 2.4e10;
   geo.disk_mdot = 1.e-8 * MSOL / YR;
+  geo.colour_correction = FCOL_OFF;
 
   geo.t_bl = 100000.;
 
@@ -323,7 +324,7 @@ init_advanced_modes ()
   modes.use_debug = FALSE;
   modes.print_dvds_info = FALSE;        // print information on velocity gradients
   modes.quit_after_inputs = FALSE;      // check inputs and quit
-  modes.quit_after_wind_defined = FALSE;        // define wind and quit
+  modes.quit_after_wind_defined = FALSE;      // define wind and quit
   modes.fixed_temp = FALSE;     // do not attempt to change temperature 
   modes.zeus_connect = FALSE;   // connect with zeus
 

source/setup_disk.c

@@ -67,9 +67,18 @@ get_disk_params ()
   geo.disk_radiation = rdchoice ("Disk.radiation(yes,no)", "1,0", answer);
 
   if (geo.disk_radiation)
+  {
     get_spectype (geo.disk_radiation,
-                  //"Disk.rad_type_to_make_wind(0=bb,1=models)", &geo.disk_ion_spectype);
-                  "Disk.rad_type_to_make_wind(bb,models)", &geo.disk_ion_spectype);
+                  "Disk.rad_type_to_make_wind(bb,models,mod_bb)", &geo.disk_ion_spectype);
+
+    /* if we have chosen a colour correction, ask the user which form to use */
+    if (geo.disk_ion_spectype == SPECTYPE_BB_FCOL)
+    {
+      strcpy (answer, "Done12");
+      sprintf (values, "%d", FCOL_DONE);
+      geo.colour_correction = rdchoice ("Disk.colour_correction(Done12)", values, answer);
+    }
+  }
 
 
   geo.disk_tprofile = DISK_TPROFILE_STANDARD;

source/setup_star_bh.c

@@ -109,6 +109,12 @@ get_stellar_params ()
     strcpy (answer, "yes");
     geo.star_radiation = rdchoice ("Central_object.radiation(yes,no)", "1,0", answer);
     get_spectype (geo.star_radiation, "Central_object.rad_type_to_make_wind(bb,models)", &geo.star_ion_spectype);
+
+    if (geo.star_ion_spectype == SPECTYPE_BB_FCOL)
+    {
+      Error ("Colour corrected BB not implemented for star at this stage. Exiting.");
+      Exit (0);
+    }
 
     if (geo.star_radiation)
       rddoub ("Central_object.temp", &geo.tstar_init);
@@ -195,6 +201,11 @@ get_bl_and_agn_params (lstar)
     if (geo.agn_radiation)
       get_spectype (geo.agn_radiation, "Central_object.rad_type_to_make_wind(bb,models,power,cloudy,brems,mono)", &geo.agn_ion_spectype);
 
+    if (geo.agn_ion_spectype == SPECTYPE_BB_FCOL)
+    {
+      Error ("Colour corrected BB not implemented for AGN at this stage. Exiting.");
+      Exit (0);
+    }
   }
   else
   {
@@ -204,6 +215,12 @@ get_bl_and_agn_params (lstar)
 
     if (geo.bl_radiation)
       get_spectype (geo.bl_radiation, "Boundary_layer.rad_type_to_make_wind(bb,models,power)", &geo.bl_ion_spectype);
+
+    if (geo.bl_ion_spectype == SPECTYPE_BB_FCOL)
+    {
+      Error ("Colour corrected BB not implemented for boundary layer. Exiting.");
+      Exit (0);
+    }
   }
 
 

source/templates.h

@@ -226,6 +226,7 @@ double vdisk(double x[], double v[]);
 double zdisk(double r);
 double ds_to_disk(struct photon *p, int allow_negative, int *hit);
 double disk_height(double s, void *params);
+double disk_colour_correction(double t);
 /* lines.c */
 double total_line_emission(PlasmaPtr xplasma, double f1, double f2);
 double lum_lines(PlasmaPtr xplasma, int nmin, int nmax);

source/wind_updates2d.c

@@ -992,6 +992,13 @@ WindPtr (w);
       report_bf_simple_ionpool ();
   }
 
+  /* report a warning if the induced Compton heating is greater than 10% of the heating, see #1016 */
+  if (icsum >= (0.1 * xsum))
+  {
+    Error ("!!wind_update: Induced Compton is responsible for %3.1f percent of radiative heating. Could cause problems.\n",
+           icsum / xsum * 100.0);
+  }
+
 
 
   if (modes.zeus_connect == 1 || modes.fixed_temp == 1) //There is no point in computing temperature changes, because we have fixed them!