Merge some changes that had been made in main only

README.md

@@ -11,7 +11,7 @@ Documentation is hosted on [ReadTheDocs](http://sirocco-rt.readthedocs.io/en/dev
 
 [![C/C++ CI](https://github.com/sirocco-rt/sirocco/actions/workflows/build.yml/badge.svg)](https://github.com/sirocco-rt/sirocco/actions/workflows/build.yml)
 
-[![Documentation Status](https://readthedocs.org/projects/agnwinds/badge/?version=latest)](https://agnwinds.readthedocs.io/en/latest/?badge=latest)
+[![Documentation Status](https://readthedocs.org/projects/sirocco-rt/badge/?version=latest)](https://sirocco-rt.readthedocs.io/en/latest/?badge=latest)
 
 ## Installation
 

docs/pydocs/README.md

@@ -2,3 +2,5 @@ This directory contains help files created from doc strings from
 python scripts in the py_progs directory.  
 
 To create the documentation yourself, simply run `write_docs.py`.
+
+This procedure should be largely superseded by https://agnwinds.readthedocs.io/en/dev/py_progs.html but that page is not currently complete.

docs/sphinx/requirements.txt

@@ -9,7 +9,7 @@ numpy
 matplotlib
 scipy
 sqlalchemy
-jinja2==3.0.3
+jinja2==3.1.4
 sphinx_gallery
 pyhdf
 sphinx-autoapi==3.1.1

examples/basic/tde_standard.pf

@@ -1,16 +1,16 @@
 System_type(star,cv,bh,agn,previous)                  agn
 
 ### Parameters for the Central Object
-Central_object.mass(msol)                  3e7
-Central_object.radius(cm)                  2.65e13
+Central_object.mass(msol)                  3e6
+Central_object.radius(cm)                  2.65e12
 
 ### Parameters for the Disk (if there is one)
-Disk.type(none,flat,vertically.extended)                 flat
+Disk.type(none,flat,vertically.extended,rmin>central.obj.rad)                 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)                         0.01
-Disk.radmax(cm)                            1e15
+Disk.mdot(msol/yr)                         9.98e-03
+Disk.radmax(cm)                            1.00885e+14
 
 ### Parameters for Boundary Layer or the compact object in an X-ray Binary or AGN
 Central_object.radiation(yes,no)                   no
@@ -24,52 +24,53 @@ Wind.dim.in.z_or_theta.direction           100
 
 ### Parameters associated with photon number, cycles,ionization and radiative transfer options
 Photons_per_cycle                          5e7
-Ionization_cycles                          15
+Ionization_cycles                          25
 Spectrum_cycles                            5
 Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est)           matrix_pow
 Line_transfer(pure_abs,pure_scat,sing_scat,escape_prob,thermal_trapping,macro_atoms_escape_prob,macro_atoms_thermal_trapping) macro_atoms_thermal_trapping
-Atomic_data                                data/h20_hetop_standard80.dat
+Matom_transition_mode(mc_jumps,matrix)               matrix
 Surface.reflection.or.absorption(reflect,absorb,thermalized.rerad)              reflect
 Wind_heating.extra_processes(none,adiabatic,nonthermal,both)            adiabatic
+Atomic_data                                data/h20_hetop_standard80.dat
 
 ### Parameters for Domain 0
-Wind.mdot(msol/yr)                         0.01
+Wind.mdot(msol/yr)                         9.98e-03
 SV.diskmin(units_of_rstar)                 1
-SV.diskmax(units_of_rstar)                 30
+SV.diskmax(units_of_rstar)                 3.80699e+01
 SV.thetamin(deg)                           20
 SV.thetamax(deg)                           65
-SV.mdot_r_exponent                         0
-SV.v_infinity(in_units_of_vescape          1
-SV.acceleration_length(cm)                 5e16
-SV.acceleration_exponent                   1.5
+SV.mdot_r_exponent                         2.0
+SV.v_infinity(in_units_of_vescape          0.3
+SV.acceleration_length(cm)                 8.856670e+14
+SV.acceleration_exponent                   4.000000
 SV.gamma(streamline_skew;1=usually)        1.0
 SV.v_zero_mode(fixed,sound_speed)          sound_speed
 SV.v_zero(multiple_of_sound_speed)         1
 Wind.radmax(cm)                            5e17
-Wind.t.init                                40000
+Wind.t.init                                35000
 Wind.filling_factor(1=smooth,<1=clumped)   0.1
 
 ### Parameters defining the spectra seen by observers
 
-Disk.rad_type_in_final_spectrum(bb,models,uniform)                   bb
+Disk.rad_type_in_final_spectrum(bb,models,uniform,mono,mod_bb)                   bb
 
-### The minimum and maximum wavelengths in the final spectra
-Spectrum.nwave                    10000
-Spectrum.wavemin(Angstroms)                100
-Spectrum.wavemax(Angstroms)                10000
+### The minimum and maximum wavelengths in the final spectra and the number of wavelength bins
+Spectrum.nwave                                10000
+Spectrum.wavemin(Angstroms)                    100
+Spectrum.wavemax(Angstroms)                  10000
 
 ### The observers and their location relative to the system
-Spectrum.no_observers                      5
-Spectrum.angle(0=pole)                     10
-Spectrum.angle(0=pole)                     35
-Spectrum.angle(0=pole)                     60
-Spectrum.angle(0=pole)                     75
-Spectrum.angle(0=pole)                     85
+Spectrum.no_observers                            5
+Spectrum.angle(0=pole)                           10
+Spectrum.angle(0=pole)                          35
+Spectrum.angle(0=pole)                          60
+Spectrum.angle(0=pole)                          75
+Spectrum.angle(0=pole)                          85
 Spectrum.live_or_die(live.or.die,extract)              extract
 Spectrum.type(flambda,fnu,basic)              flambda
 
 ### Parameters for Reverberation Modeling (if needed)
 Reverb.type(none,photon,wind,matom)                 none
 
 ### Other parameters
-Photon_sampling.approach(T_star,cv,yso,AGN,tde_bb,min_max_freq,user_bands,cloudy_test,wide,logarithmic)                    agn
+Photon_sampling.approach(T_star,cv,yso,AGN,tde_bb,min_max_freq,user_bands,cloudy_test,wide,logarithmic)               tde_bb

source/rad_hydro_files.c

@@ -42,6 +42,8 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
+#include <gsl/gsl_errno.h>
+
 
 #include "atomic.h"
 #include "sirocco.h"
@@ -143,7 +145,7 @@ main (argc, argv)
 
   struct photon ptest;          //We need a test photon structure in order to compute t
 
-  FILE *fptr_hc, *fptr_drive, *fptr_ion, *fptr_spec, *fptr_pcon, *fptr_debug, *fptr_flux, *fptr_flux_theta, *fptr_flux_phi, *fptr_flux_r, *fopen ();  /*This is the file to communicate with zeus */
+  FILE *fptr_hc, *fptr_drive, *fptr_ion, *fptr_spec, *fptr_pcon, *fptr_debug, *fptr_flux, *fptr_flux_theta, *fptr_flux_phi, *fptr_flux_r, *fopen ();    /*This is the file to communicate with zeus */
   domain = geo.hydro_domain_number;
 
   /* Initialize  MPI, which is needed because some of the routines are MPI enabled */
@@ -166,7 +168,8 @@ main (argc, argv)
 
   /* MPI intialiazation is complete */
 
-
+  /* Disable GSL error handling, so we can handle errors ourselves */
+  gsl_set_error_handler_off ();
 
   strcpy (parameter_file, "NONE");
 
@@ -181,9 +184,13 @@ main (argc, argv)
   strcat (windsavefile, ".wind_save");
   strcat (outputfile, ".txt");
 
-
-/* Read in the wind file */
-
+  /* Allocate memory required for domain, and read in the wind file */
+  zdom = calloc (MAX_DOM, sizeof (domain_dummy));
+  if (zdom == NULL)
+  {
+    Error ("Failed to allocate memory for domain\n");
+    return (EXIT_FAILURE);
+  }
   if (wind_read (windsavefile) < 0)
   {
     Error ("swind: Could not open %s", windsavefile);
@@ -269,7 +276,7 @@ main (argc, argv)
   }
   else if (zdom[domain].coord_type == RTHETA)
   {
-    fprintf (fptr_drive, "i j rcen thetacen vol rho ne F_vis_x F_vis_y F_vis_z F_vis_mod F_UV_theta F_UV_phi F_UV_r F_UV_mod F_Xray_x F_Xray_y F_Xray_z F_Xray_mod es_f_x es_f_y es_f_z es_f_mod bf_f_x bf_f_y bf_f_z bf_f_mod\n");   //directional flux by band
+    fprintf (fptr_drive, "i j rcen thetacen vol rho ne F_vis_x F_vis_y F_vis_z F_vis_mod F_UV_theta F_UV_phi F_UV_r F_UV_mod F_Xray_x F_Xray_y F_Xray_z F_Xray_mod es_f_x es_f_y es_f_z es_f_mod bf_f_x bf_f_y bf_f_z bf_f_mod\n");     //directional flux by band
     fprintf (fptr_flux, "i j rcen thetacen F_vis_x F_vis_y F_vis_z F_vis_mod F_UV_x F_UV_y F_UV_z F_UV_mod F_Xray_x F_Xray_y F_Xray_z F_Xray_mod\n");   //directional flux by band
   }
 
@@ -305,9 +312,18 @@ main (argc, argv)
       nplasma = wmain[nwind].nplasma;
       wind_n_to_ij (domain, plasmamain[nplasma].nwind, &i, &j);
 
-      fprintf (fptr_flux_theta, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);  //output geometric things
-      fprintf (fptr_flux_phi, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);  //output geometric things
-      fprintf (fptr_flux_r, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);  //output geometric things
+      if (zdom[domain].coord_type == SPHERICAL || zdom[domain].coord_type == RTHETA)
+      {
+        fprintf (fptr_flux_theta, "%3d %3d      %3d %10.8e %10.8e ", i, j, wmain[nwind].inwind, wmain[nwind].rcen, wmain[nwind].thetacen / RADIAN);     //output geometric things
+        fprintf (fptr_flux_phi, "%3d %3d      %3d %10.8e %10.8e ", i, j, wmain[nwind].inwind, wmain[nwind].rcen, wmain[nwind].thetacen / RADIAN);       //output geometric things
+        fprintf (fptr_flux_r, "%3d %3d      %3d %10.8e %10.8e ", i, j, wmain[nwind].inwind, wmain[nwind].rcen, wmain[nwind].thetacen / RADIAN); //output geometric things
+      }
+      else if (zdom[domain].coord_type == CYLIND)
+      {
+        fprintf (fptr_flux_theta, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);    //output geometric things
+        fprintf (fptr_flux_phi, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);      //output geometric things
+        fprintf (fptr_flux_r, "%3d %3d      %3d %10.3e %10.3e ", i, j, wmain[nwind].inwind, wmain[nwind].xcen[0], wmain[nwind].xcen[2]);        //output geometric things
+      }
       for (ii = 0; ii < NFLUX_ANGLES; ii++)
       {
         fprintf (fptr_flux_theta, "%10.3e ", plasmamain[nplasma].F_UV_ang_theta_persist[ii]);