Implementing advanced mode with -d flag to access

examples/sv_detailedmode.pf

@@ -0,0 +1,87 @@
+Wind_type()                  0 
+Atomic_data                       data/standard77
+write_atomicdata 					1
+photons_per_cycle                          10000
+Ionization_cycles                               1
+spectrum_cycles                                 100
+Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var)                    1
+Wind.dim.in.x_or_r.direction                     5
+Wind.dim.in.z_or_theta.direction                   5
+adjust_grid(0=no,1=yes)                       	1
+geo.xlog_scale 									1e7
+geo.zlog_scale 									1e7
+Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow)                  3 
+Line_transfer()                    3
+Thermal_balance_options(0=everything.on,1=no.adiabatic) 0
+System_type(0=star+disk,1=star+disk+secondary,2=agn)                   0
+Star_radiation(y=1)				1
+Disk_radiation(y=1)                               1
+Boundary_layer_radiation(y=1)			1
+Wind_radiation(y=1)				0
+Rad_type_for_star(0=bb,1=models)_to_make_wind				0
+Rad_type_for_disk(0=bb,1=models)_to_make_wind                    0
+Rad_type_for_bl(0=bb,1=models)_to_make_wind			0
+mstar(msol)                                   0.8 
+rstar(cm)                              7e8 
+tstar					40000
+disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk)                    1
+disk.mdot(msol/yr)                             1e-8
+Disk.illumination.treatment()                    0
+Disk.temperature.profile(0=standard;1=readin)                    0
+disk.radmax(cm)                         2.38e10
+lum_bl					0
+t_bl					10000
+Torus(0=no,1=yes)                               0 
+wind.radmax(cm)                             1e12
+wind.t.init                                 40000
+wind.mdot(msol/yr)                            1e-20
+sv.diskmin(wd_rad)                              4
+sv.diskmax(wd_rad)                             12
+sv.thetamin(deg)                                20
+sv.thetamax(deg)                                65
+sv.mdot_r_exponent                              0
+sv.v_infinity(in_units_of_vescape                   3
+sv.acceleration_length(cm)                    7e10
+sv.acceleration_exponent                         1.5
+wind.filling_factor(1=smooth,<1=clumpted)		1
+Rad_type_for_star()_in_final_spectrum						0
+Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum                    0
+Model_file			data/disk14/kur_tlusty_hybrid.ls
+spectrum_wavemin                             4000
+spectrum_wavemax                             6900 
+no_observers                                     10
+angle(0=pole)                                   0
+angle(0=pole)                                   5
+angle(0=pole)                                   10
+angle(0=pole)                                   15
+angle(0=pole)                                   20
+angle(0=pole)                                   30
+angle(0=pole)                                   40
+angle(0=pole)                                   50
+angle(0=pole)                                   60
+angle(0=pole)                                   80
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+phase(0=inferior_conjunction)                   0.5
+live.or.die(0).or.extract(anything_else)                    1
+Select_specific_no_of_scatters_in_spectra(y/n)                    n
+Select_photons_by_position(y/n)                    n
+spec.type(flambda(1),fnu(2),basic(other)                    1
+Use.standard.care.factors(1=yes)                    1
+Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined),                   2
+Extra.diagnostics(0=no,1=yes)       1
+ispymode(1=yes)             		1
+keep_ioncycle_windsaves(1=yes)     	1
+track_resonant_scatters(1=yes)      1
+print_windrad_summary(1=yes) 		1
+print_dvds_info(1=yes) 				1
+save_extract_photons(1=yes) 		1
+

source/Makefile

@@ -187,7 +187,7 @@ py_wind_objects = py_wind.o get_atomicdata.o py_wind_sub.o windsave.o py_wind_io
 		matom.o estimators.o yso.o elvis.o photon2d.o cylindrical.o rtheta.o spherical.o  \
 		cylind_var.o bilinear.o gridwind.o py_wind_macro.o partition.o auger_ionization.o\
 		spectral_estimators.o shell_wind.o compton.o torus.o zeta.o dielectronic.o \
-                variable_temperature.o bb.o rdpar.o log.o direct_ion.o
+        variable_temperature.o bb.o rdpar.o log.o direct_ion.o diag.o 
 
 
 

source/atomic.h

@@ -497,3 +497,7 @@ typedef struct gaunt_total
 } Gaunt_total, *Gaunt_totalptr;
 
 Gaunt_total gaunt_total[MAX_GAUNT_N_GSQRD];	//Set up the structure
+
+/* a variable which controls whether to save a summary of atomic data
+   this is defined in atomic.h, rather than the modes structure */
+int write_atomicdata;

source/corona.c

@@ -63,8 +63,12 @@ get_corona_params ()
   geo.wind_thetamin = 0.0;
   geo.wind_thetamax = 0.0;
 
-  geo.xlog_scale = 0.3 * geo.corona_rmin;
-  geo.zlog_scale = 0.3 * geo.corona_scale_height;
+  /* if modes.adjust_grid is 1 then we have already adjusted the grid manually */
+  if (modes.adjust_grid == 0)
+    {
+      geo.xlog_scale = 0.3 * geo.corona_rmin;
+      geo.zlog_scale = 0.3 * geo.corona_scale_height;
+    }
 
 /* Prior to 01dec, windcones were defined here.  But this broke a capability to continue
    a calculation.  To fix this, wind_cone definition was moved backed to python.c.  To

source/diag.c

@@ -1,9 +1,14 @@
+/* diag.c contains a number of routines relating to extra diagnositics
+   and related file i/o. 
+   JM 1410 -- moved a few routines here relating to extra diagnostics
+*/
 
 /***********************************************************
                                        Southampton
 
  Synopsis:
 	 open_diagfile sets up diagnostic files for use in python when the extra.diagnostics flag is 
+
 Arguments:
 
 	PhotPtr p;	the photon
@@ -20,6 +25,7 @@
 
 	12jun 	nsh	72 Added lines to set up a file for outputting photons in given cells. The cells are read in from a file called diag_cells.dat. The existance of this file defined wether the diagnostic takes place.
 	13jul	jm	changed print statements to logs and made more descriptive
+	14oct   jm  changed to reflect move to 'advanced mode', see #111 and #120
 
 **************************************************************/
 
@@ -29,6 +35,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <string.h>
 
 #include "atomic.h"
 #include "python.h"
@@ -42,14 +49,14 @@ open_diagfile ()
   FILE *cellfile;		/*File that may or may not exist, pointing to cells we want to write out photon stats for */
   int cell;			/*Temporary storage of cell to use */
 
-  if (eplinit == 0)
+  if (eplinit == 0 && modes.save_extract_photons)
     {
       epltptr = fopen ("python.ext", "w");
       eplinit = 1;
     }
 
   ncstat = 0;			/*Zero the counter for the number of cells to be tracked */
-  if (pstatinit == 0)		/* Check we havent already done this */
+  if (pstatinit == 0 && modes.save_cell_stats)		/* Check we havent already done this */
     {
       cellfile = fopen ("diag_cells.dat", "r");	/*This is the file containing cells to track */
       if (cellfile != NULL)	/*If there actually *is* a file read it */
@@ -82,6 +89,187 @@ open_diagfile ()
       pstatinit = 1;		/* We have initialised this routine */
     }
 
-//  diag_on_off = 0;            // 0=off everything else is on
+  // diag_on_off = 0;            // 0=off everything else is on
+  return (0);
+}
+
+
+
+
+
+
+/***********************************************************
+				University of Southampton
+
+Synopsis:
+	get_extra_diagnostics reads in extra diagnostics if 
+	the user has asked for them. It uses rd_extra() in rdpar.c
+	to get the actual lines, and then checks them against 
+	hardwired options via a number of if loops.
+
+Arguments:	
+    none	
+
+Returns:
+    modifies the modes structure to turn on whichever modes
+    are asked for.
+
+Description:	
+
+Notes:
+    see #111 and #120
+
+History:
+    1410 -- JM -- Coded
+**************************************************************/
+
+int get_extra_diagnostics()
+{
+  int wordlength;
+  char firstword[LINELENGTH];
+  int noptions, rdstat;
+  double answer;
+
+
+  noptions = 0;
+  wordlength = 0;
+  rdstat = 0;
+
+  if (modes.iadvanced == 0)
+  	Error("Getting extra_diagnostics but advanced mode is off!\n");
+
+  Log("get_extra_diagnostics: Getting extra diagnostics as requested...\n");
+
+  while (noptions < NMAX_OPTIONS && rdstat == 0)
+  {
+
+    rdstat = rd_extra(firstword, &answer, &wordlength);
+
+
+    /* if rdstat is 1 we reached EOF, so exit this routine */
+    if (rdstat)
+      {
+      	if (noptions == 0)
+      	  Error("get_extra_diagnostics: EOF: %i options read, but extra diagnostics on!\n", noptions);
+
+        return(0);
+      }
+
+    Log("%s %8.4e\n", firstword, answer);
+
+    /* would you like to save cell photon statistics */
+	if (strncmp ("save_cell_statistics", firstword, wordlength) == 0)
+	  {
+	    modes.save_cell_stats = answer;
+  	    Log("You are tracking photon statistics by cell\n");
+	  }
+
+    /* would you like to use ispy mode */
+	else if (strncmp ("ispymode", firstword, wordlength) == 0)
+      {
+	    modes.ispy = answer;
+	    Log("ISPY mode is on\n");
+	  }
+
+	/* would you like to track resonant scatters */
+	else if (strncmp ("track_resonant_scatters", firstword, wordlength) == 0)
+      {
+	    modes.track_resonant_scatters = answer;
+	    Log("You are tracking resonant scatters\n");
+	  }
+
+    /* would you like keep windsave files for each cycle */
+	else if (strncmp ("keep_ioncycle_windsaves", firstword, wordlength) == 0)
+      {
+	    modes.keep_ioncycle_windsaves = answer;
+	    Log("You are keeping windsave files for each cycle\n");
+	  }
+
+	/* would you like to save data on extract */
+	else if (strncmp ("save_extract_photons", firstword, wordlength) == 0)
+      {
+	    modes.save_extract_photons = answer;
+	    Log("You are saving data on extract\n");
+	  }
+
+	/* would you like to print wind_rad_summary*/
+	else if (strncmp ("print_windrad_summary", firstword, wordlength) == 0)
+      {
+	    modes.print_windrad_summary = answer;
+	    Log("You are printing wind_rad_summary\n");
+	  }
+
+	/* would you like to print dvds_info */
+	else if (strncmp ("print_dvds_info", firstword, wordlength) == 0)
+      {
+	    modes.print_dvds_info = answer;
+	    Log("You are printing dvds_info\n");
+	  }
+
+    else
+      {
+      	Error("get_extra_diagnostics: didn't understand question %s, continuing!\n", firstword); 
+      	noptions--;		// this isn't a real option, so decrement it before we increment it...!
+      }
+
+    noptions++;	// increment noptions
+  }
+
+  return 0;
+}
+
+
+
+
+/***********************************************************
+                Southampton University
+
+Synopsis: 
+	save_photon_stats prints photon statistics to a file
+
+Arguments:	
+	One 		WindPtr for the cell
+	p 			Photon pointer
+	ds 			ds travelled
+
+Returns:
+
+Description:
+   the loop below is if the user requires extra diagnostics and
+   has provided a file diag_cells.dat to store photons stats for cells they have specified
+
+Notes:
+   Moved here to save duplicating code between bf_estimators_increment and radiation.
+
+History:
+   1410 JM 		Coding began
+   1410 JM      Moved here from python.c	
+
+**************************************************************/
+
+
+
+int
+save_photon_stats (one, p, ds)
+     WindPtr one;
+     PhotPtr p;
+     double ds;
+{
+  int i;
+
+  /* JM -- 1310 -- the loop below is if the user requires extra diagnostics and
+     has provided a file diag_cells.dat to store photons stats for cells they have specified
+   */
+
+  for (i = 0; i < ncstat; i++)
+    {
+      /* check if the cell is in the specified list - ncell_stats is global variable */
+      if (one->nplasma == ncell_stats[i])
+	{
+	  fprintf (pstatptr,
+		   "PHOTON_DETAILS %3d %8.3e %8.3e %8.3e cell%3d wind cell%3d\n",
+		   geo.wcycle, p->freq, p->w, ds, one->nplasma, one->nwind);
+	}
+    }
   return (0);
 }

source/elvis.c

@@ -117,8 +117,13 @@ get_elvis_wind_params ()
   geo.wind_thetamin = geo.sv_thetamin;
   geo.wind_thetamax = geo.sv_thetamax;
 //OLD  geo.xlog_scale = geo.sv_rmin + (geo.elvis_offset * tan (geo.sv_thetamin));
-  geo.xlog_scale = geo.sv_rmin;
-  geo.zlog_scale = 1e15;	/* Big number - for AGN */
+
+  /* if modes.adjust_grid is 1 then we have already adjusted the grid manually */
+  if (modes.adjust_grid == 0)
+    {
+      geo.xlog_scale = geo.sv_rmin;
+      geo.zlog_scale = 1e15;	/* Big number - for AGN */
+    }
 
 /*Now calculate the normalization factor for the wind*/
 

source/estimators.c

@@ -86,7 +86,7 @@ bf_estimators_increment (one, p, ds)
   /* JM -- 1310 -- check if the user requires extra diagnostics and
      has provided a file diag_cells.dat to store photons stats for cells they have specified
    */
-  if (diag_on_off == 1 && ncstat > 0)
+  if (modes.save_cell_stats && ncstat > 0)
     {
       save_photon_stats(one, p, ds);  // save photon statistics (extra diagnostics)
     }

source/extract.c

@@ -156,7 +156,7 @@ one is odd. We do frequency here but weighting is carried out in  extract */
 
 	    }
 
-	  if (diag_on_off && 1545.0 < 2.997925e18 / pp.freq
+	  if (modes.save_extract_photons && 1545.0 < 2.997925e18 / pp.freq
 	      && 2.997925e18 / pp.freq < 1565.0)
 	    {
 	      fprintf (epltptr,