From e48d0d9403722194d4565342270fdb4e3ee409da Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 9 Aug 2024 12:52:04 -0500 Subject: [PATCH 01/23] Dummy routine for windsave2fits --- source/Makefile | 7 + source/windsave2fits.c | 571 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 578 insertions(+) create mode 100644 source/windsave2fits.c diff --git a/source/Makefile b/source/Makefile index 31646e821..7afdde127 100644 --- a/source/Makefile +++ b/source/Makefile @@ -264,6 +264,7 @@ rhf_objects = rad_hydro_files.o $(python_objects) modify_wind_objects = modify_wind.o $(python_objects) inspect_wind_objects = inspect_wind.o $(python_objects) py_optd_obj = py_optd.o py_optd_output.o py_optd_trans.o py_optd_util.o py_optd_extern_init.o $(python_objects) +windsave2fits_objects = windsave2fits.o $(python_objects) run_indent: ../py_progs/run_indent.py -all_no_headers @@ -280,6 +281,12 @@ windsave2table: startup $(table_objects) $(CUDA_OBJECTS) mv $@ $(BIN)/windsave2table$(VERSION) @if [ $(INDENT) = yes ] ; then ../py_progs/run_indent.py -changed ; fi +windsave2fits: startup $(windsave2fits_objects) $(CUDA_OBJECTS) + $(CC) $(CFLAGS) $(windsave2fits_objects) $(CUDA_OBJECTS) $(LDFLAGS) -lcfitsio -o windsave2fits + cp $@ $(BIN) + mv $@ $(BIN)/windsave2fits$(VERSION) + @if [ $(INDENT) = yes ] ; then ../py_progs/run_indent.py -changed ; fi + rad_hydro_files: startup $(rhf_objects) $(CUDA_OBJECTS) $(CC) $(CFLAGS) $(rhf_objects) $(CUDA_OBJECTS) $(LDFLAGS) -o rad_hydro_files cp $@ $(BIN) diff --git a/source/windsave2fits.c b/source/windsave2fits.c new file mode 100644 index 000000000..26b5848f3 --- /dev/null +++ b/source/windsave2fits.c @@ -0,0 +1,571 @@ + +/***********************************************************/ +/** @file inspect_wind.c + * @author ksl + * @date October, 2021 + * + * @brief Routines to inspect variables in a wind structure + * + *###Notes### + * This is intended just as a diagnostic routine + * so that one can print out whatever variables in + * a windstrucutre one wants in order to diagnose + * a problem. It was written so that we could inspect + * some of the macro atom variables in paralell mode + * in diagnosing issue #898 and #910, but anyon + * should change it so that other problems might + * be addressed. + * + * + ***********************************************************/ + +#include +#include +#include +#include + +#include "atomic.h" +#include "python.h" +#include "fitsio.h" + + +char inroot[LINELENGTH], outroot[LINELENGTH], model_file[LINELENGTH], folder[LINELENGTH]; +int model_flag, ksl_flag, cmf2obs_flag, obs2cmf_flag; + +double line_matom_lum_single (double lum[], PlasmaPtr xplasma, int uplvl); +int line_matom_lum (int uplvl); +int create_matom_level_map (); + +/**********************************************************/ +/** + * @brief parses the command line options + * + * @param [in] int argc the number of command line arguments + * @param [in] char * argv[] The command line arguments + * + * + * ###Notes### + * + * The general purpose of each of the command line options + * should be fairly obvious from reading the code. + * + * + * Although this routine uses the standard Log and Error commands + * the diag files have not been created yet and so this information + * is really simply written to the terminal. + * + **********************************************************/ + +int +xparse_command_line (argc, argv) + int argc; + char *argv[]; +{ + int j = 0; + int i; + char dummy[LINELENGTH]; + int mkdir (); + char *fgets_rc; + + + sprintf (outroot, "%s", "new"); + + model_flag = ksl_flag = obs2cmf_flag = cmf2obs_flag = 0; + + if (argc == 1) + { + printf ("Parameter file name (e.g. my_model.pf, or just my_model):"); + fgets_rc = fgets (dummy, LINELENGTH, stdin); + if (!fgets_rc) + { + printf ("Input rootname is NULL or invalid\n"); + exit (1); + } + get_root (inroot, dummy); + } + else + { + + for (i = 1; i < argc; i++) + { + if (strcmp (argv[i], "-out_root") == 0) + { + if (sscanf (argv[i + 1], "%s", dummy) != 1) + { + printf ("python: Expected out_root after -out_root switch\n"); + exit (0); + } + + get_root (outroot, dummy); + i++; + j = i; + + } + if (strcmp (argv[i], "-model_file") == 0) + { + if (sscanf (argv[i + 1], "%s", dummy) != 1) + { + printf ("python: Expected a model file containing density, velocity and temperature after -model_file switch\n"); + exit (0); + } + get_root (model_file, dummy); + i++; + j = i; + printf ("got a model file %s\n", model_file); + model_flag = 1; + } + else if (strcmp (argv[i], "-ksl") == 0) + { + printf ("Carrying out a simple hard wired ion modification\n"); + ksl_flag = 1; + } + else if (strcmp (argv[i], "--dry-run") == 0) + { + modes.quit_after_inputs = 1; + j = i; + } + else if (strcmp (argv[i], "-cmf") == 0) + { + obs2cmf_flag = 1; + } + else if (strcmp (argv[i], "-obs") == 0) + { + cmf2obs_flag = 1; + } + else if (strncmp (argv[i], "-", 1) == 0) + { + printf ("python: Unknown switch %s\n", argv[i]); + exit (0); + } + } + + /* The last command line variable is always the windsave file */ + + if (j + 1 == argc) + { + printf ("All of the command line has been consumed without specifying a parameter file name, so exiting\n"); + exit (1); + } + strcpy (dummy, argv[argc - 1]); + get_root (inroot, dummy); + + } + + return (0); +} + + + +/**********************************************************/ +/** + * @brief the main routine which carries out the effort + * + * @param [in] int argc the number of command line arguments + * @param [in] char * argv[] The command line arguments + * + * + * ###Notes### + * + * This routine oversees the effort. The basic steps are + * + * - parse the command line to get the names of files + * - read the old windsave file + * - read the densities from in this case H + * - modify the densities + * - write out the new windsave file + * + * + **********************************************************/ + + +int +main (argc, argv) + int argc; + char *argv[]; +{ + + char infile[LINELENGTH], outfile[LINELENGTH]; + int n, i; + FILE *fptr, *fopen (); + int ii, jj, ndom, nnwind; + int mkdir (); + + + xparse_command_line (argc, argv); + + sprintf (infile, "%.150s.wind_save", inroot); + sprintf (outfile, "%.150s.txt", inroot); + + printf ("Reading %s and writing to %s\n", infile, outfile); + + zdom = calloc (MAX_DOM, sizeof (domain_dummy)); + if (zdom == NULL) + { + printf ("Unable to allocate memory for domain\n"); + return EXIT_FAILURE; + } + + wind_read (infile); + + if (nlevels_macro == 0) + { + printf ("Currently this routine only looks at macro atom values, and this is a simple atom file\n"); + exit (0); + } + + + + fptr = fopen (outfile, "w"); + + fprintf (fptr, "# Results for %s\n", infile); + + fprintf (fptr, "# Size %d %d %d \n", size_Jbar_est, size_gamma_est, size_alpha_est); + + fprintf (fptr, "%15s %4s %4s %4s ", "Variable", "np", "i", "j"); + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "MacLev%02d ", i); + } + fprintf (fptr, "\n"); + + fprintf (fptr, "%15s %4s %4s %4s ", "---------------", "----", "----", "----"); + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "-------- "); + } + fprintf (fptr, "\n"); + + for (n = 0; n < NPLASMA; n++) + { + + fprintf (fptr, "%15s", "jbar"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].jbar[i]); + } + fprintf (fptr, "\n"); + } + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "jbar_old"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].jbar_old[i]); + } + fprintf (fptr, "\n"); + } + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "gamma"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].gamma[i]); + } + fprintf (fptr, "\n"); + } + + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "alpha_st"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].alpha_st[i]); + } + fprintf (fptr, "\n"); + } + + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "recomb_sp"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].recomb_sp[i]); + } + fprintf (fptr, "\n"); + } + + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "matom_abs"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].matom_abs[i]); + } + fprintf (fptr, "\n"); + } + + + + + for (n = 0; n < NPLASMA; n++) + { + fprintf (fptr, "%15s", "matom_emiss"); + + nnwind = plasmamain[n].nwind; + ndom = wmain[nnwind].ndom; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, " %4d %4d %4d ", n, ii, jj); + + for (i = 0; i < nlevels_macro; i++) + { + fprintf (fptr, "%8.2e ", macromain[n].matom_emiss[i]); + } + fprintf (fptr, "\n"); + } + fclose (fptr); + + + /* calculate all the line luminosities in macro-atom mode */ + /* these are stored in a folder with a separate file for each upper level */ + sprintf (folder, "matom_linelum_%.200s", inroot); + mkdir (folder, 0777); + printf ("Calculating macro-atom line luminosities for all lines in wavelength range %.1f to %.1f Angstroms...\n", + VLIGHT / geo.sfmax / ANGSTROM, VLIGHT / geo.sfmin / ANGSTROM); + + /* create a file that tells the user the wavelengths of every macro-atom line and the + corresponding upper and lower levels */ + create_matom_level_map (); + + /* loop over all macro atom upper levels */ + for (i = 0; i < nlevels_macro; i++) + { + line_matom_lum (i); + } + printf ("Done for %d levels.\n", nlevels_macro); + + exit (0); + +} + +/**********************************************************/ +/** + * @brief print out which line wavelengths correspond to which upper and lower levels + * + **********************************************************/ + +int +create_matom_level_map () +{ + int uplvl, nbbd, n; + char outfile[LINELENGTH]; + FILE *fptr, *fopen (); + + /* open a file in the folder where we store the matom line luminosities */ + sprintf (outfile, "%.200s/line_map.txt", folder); + + /* print some header information to the file */ + fptr = fopen (outfile, "w"); + fprintf (fptr, "# model name %s\n", inroot); + fprintf (fptr, "# which line wavelengths (Angstroms) correspond to which upper and lower levels\n"); + fprintf (fptr, "upper lower wavelength\n"); + + for (uplvl = 0; uplvl < nlevels_macro; uplvl++) + { + nbbd = xconfig[uplvl].n_bbd_jump; + for (n = 0; n < nbbd; n++) + { + fprintf (fptr, "%d %d %12.2f\n", uplvl, line[xconfig[uplvl].bbd_jump[n]].nconfigl, + VLIGHT / line[xconfig[uplvl].bbd_jump[n]].freq / ANGSTROM); + } + } + fclose (fptr); + return (0); +} + + +/**********************************************************/ +/** + * @brief calculate line luminosities for all cells for a given upper level and save to file + * + **********************************************************/ + +int +line_matom_lum (uplvl) + int uplvl; +{ + int n, nbbd, i, ii, jj, nnwind, ndom, inwind; + double lum[NBBJUMPS]; + char outfile[LINELENGTH]; + FILE *fptr, *fopen (); + + nbbd = xconfig[uplvl].n_bbd_jump; + + /* open a file in the folder where we store the matom line luminosities */ + sprintf (outfile, "%.100s/linelums_%.100s_lev%d.txt", folder, inroot, uplvl); + + /* print some header information to the file */ + fptr = fopen (outfile, "w"); + fprintf (fptr, "# model name %s\n", inroot); + fprintf (fptr, "# Line luminosities from upper level %d from %.2f to %.2f Angstroms\n", uplvl, VLIGHT / geo.sfmax / ANGSTROM, + VLIGHT / geo.sfmin / ANGSTROM); + fprintf (fptr, "# Lower Levels: "); + for (n = 0; n < nbbd; n++) + { + fprintf (fptr, "%12d", line[xconfig[uplvl].bbd_jump[n]].nconfigl); + } + fprintf (fptr, "\n# Line Wavelengths:"); + for (n = 0; n < nbbd; n++) + { + fprintf (fptr, " %12.2f", VLIGHT / line[xconfig[uplvl].bbd_jump[n]].freq / ANGSTROM); + } + fprintf (fptr, "\n"); + + /* columns of wind cells */ + fprintf (fptr, "%12s %12s %12s %4s %4s %12s %12s %12s", "nwind", "x", "z", "i", "j", "inwind", "wind_vol", "plasma_vol"); + + /* header info for each lower level */ + for (n = 0; n < nbbd; n++) + { + fprintf (fptr, " LowerLev%03d ", line[xconfig[uplvl].bbd_jump[n]].nconfigl); + } + fprintf (fptr, "\n"); + // fprintf (fptr, "%12s %12s %12s %4s %4s %12s %12s %12s", "------------", "------------", "------------", "----", "----", "------------", + // "------------", "------------"); + // for (n = 0; n < nbbd; n++) + // { + // fprintf (fptr, " ----------- "); + // } + // fprintf (fptr, "\n"); + + + /* loop over each plasma cell and do the calculation */ + for (nnwind = 0; nnwind < NDIM2; nnwind++) + { + ndom = wmain[nnwind].ndom; + inwind = wmain[nnwind].inwind; + wind_n_to_ij (ndom, nnwind, &ii, &jj); + fprintf (fptr, "%12d %12.4e %12.4e %4d %4d %12d %12.4e", nnwind, wmain[nnwind].xcen[0], wmain[nnwind].xcen[2], ii, jj, inwind, + wmain[nnwind].vol); + if (wmain[nnwind].inwind >= 0) + { + n = wmain[nnwind].nplasma; + line_matom_lum_single (lum, &plasmamain[n], uplvl); + /* print the filled volume */ + fprintf (fptr, " %13.4e", plasmamain[n].vol); + for (i = 0; i < nbbd; i++) + fprintf (fptr, " %13.4e", lum[i]); + } + else + { + /* print a 0.0 instead of the filled volume and line lums outside the wind */ + fprintf (fptr, " %13.4e", 0.0); + for (i = 0; i < nbbd; i++) + fprintf (fptr, " %13.4e", 0.0); + } + + fprintf (fptr, "\n"); + } + fclose (fptr); + return (0); +} + +/**********************************************************/ +/** + * @brief calculate line luminosities for a single cell for a given upper level + * + **********************************************************/ + +double +line_matom_lum_single (lum, xplasma, uplvl) + double lum[]; + PlasmaPtr xplasma; + int uplvl; +{ + int n, nbbd, m; + double penorm, bb_cont; + double freq_min, freq_max, lum_tot; + struct lines *line_ptr; + double eprbs[NBBJUMPS]; + freq_min = geo.sfmin; + freq_max = geo.sfmax; + /* identify number of bb downward jumps */ + nbbd = xconfig[uplvl].n_bbd_jump; + penorm = 0.0; + m = 0; + lum_tot = 0.0; + /* work out how often we come out in each of these locations */ + for (n = 0; n < nbbd; n++) + { + line_ptr = &line[xconfig[uplvl].bbd_jump[n]]; + if ((line_ptr->freq > freq_min) && (line_ptr->freq < freq_max)) // correct range + { + // bb_cont = (a21 (line_ptr) * p_escape (line_ptr, xplasma)); + bb_cont = (a21 (line_ptr) * 1.0); + eprbs[m] = bb_cont * (xconfig[uplvl].ex - xconfig[line[xconfig[uplvl].bbd_jump[n]].nconfigl].ex); //energy difference + penorm += eprbs[m]; + } + else + { + eprbs[m] = 0.0; + } + m++; + } + + /* correctly normalise the probabilities */ + for (n = 0; n < nbbd; n++) + { + if (penorm == 0) + { + lum[n] = 0.0; + } + else + { + eprbs[n] = eprbs[n] / penorm; + lum[n] = eprbs[n] * macromain[xplasma->nplasma].matom_emiss[uplvl]; + } + lum_tot += lum[n]; + } + + return (lum_tot); +} From fc1c294963efcf7a94177ac27bb116cce5e26de6 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 9 Aug 2024 17:07:04 -0500 Subject: [PATCH 02/23] First version of windsave2fits that works --- source/windsave2fits.c | 232 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 231 insertions(+), 1 deletion(-) diff --git a/source/windsave2fits.c b/source/windsave2fits.c index 26b5848f3..88d5dfda2 100644 --- a/source/windsave2fits.c +++ b/source/windsave2fits.c @@ -1,6 +1,6 @@ /***********************************************************/ -/** @file inspect_wind.c +/** @file windsave2fits.c * @author ksl * @date October, 2021 * @@ -36,6 +36,233 @@ double line_matom_lum_single (double lum[], PlasmaPtr xplasma, int uplvl); int line_matom_lum (int uplvl); int create_matom_level_map (); +// Define a structure to hold spectral data +typedef struct +{ + int num_wavelengths; + int num_spectra; + float **data; // 2D array to hold spectra data +} Spectra; + +/* Next fucntions are utilities that it should be possitle to call repeatedly*/ + +// Function to prepare 2D image data from a 2D array +float * +prepare_image_data (float **data, int width, int height) +{ + // Allocate memory for a contiguous 1D array to store the 2D image data + float *image_data = calloc (width * height, sizeof (float)); + if (!image_data) + { + fprintf (stderr, "Memory allocation error\n"); + return NULL; + } + + // Flatten the 2D data into a 1D array for FITS writing + for (int i = 0; i < height; i++) + { + for (int j = 0; j < width; j++) + { + image_data[i * width + j] = data[i][j]; + } + } + + return image_data; +} + +// Function to prepare binary table data +void +prepare_binary_table (int num_rows, int num_columns, float **table_data) +{ + // Allocate memory for each row of the table + for (int i = 0; i < num_rows; i++) + { + table_data[i] = malloc (num_columns * sizeof (float)); + if (!table_data[i]) + { + fprintf (stderr, "Memory allocation error\n"); + return; + } + + // Fill each row with some sample data + for (int j = 0; j < num_columns; j++) + { + table_data[i][j] = (float) (i * num_columns + j); // Example data + } + } +} + + + +// Function to write a 2D image to a FITS extension with a given name +int +write_image_extension (fitsfile *fptr, float *image_data, int width, int height, const char *extname) +{ + int status = 0; // CFITSIO status value MUST be initialized to zero + long naxes[2]; // Array to store dimension sizes + + // Set the dimensions of the 2D image: [width, height] + naxes[0] = width; // Number of columns + naxes[1] = height; // Number of rows + + // Create a new image extension of type FLOAT_IMG + if (fits_create_img (fptr, FLOAT_IMG, 2, naxes, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Set the extension name using the EXTNAME keyword + if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Write the image data to the FITS extension + if (fits_write_img (fptr, TFLOAT, 1, width * height, image_data, &status)) + { + fits_report_error (stderr, status); + return status; + } + + return 0; +} + +// Function to write a binary table to a FITS extension with a given name +int +write_binary_table_extension (fitsfile *fptr, int num_rows, int num_columns, float **table_data, const char *extname) +{ + int status = 0; + char *ttype[] = { "COL1" }; // Name for each column + char *tform[] = { "E" }; // Format for each column (E = float) + char *tunit[] = { "unit" }; // Unit for each column + + // Create a new binary table extension + if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_columns, ttype, tform, tunit, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Set the extension name using the EXTNAME keyword + if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Write the table data + for (int i = 0; i < num_rows; i++) + { + if (fits_write_col (fptr, TFLOAT, 1, i + 1, 1, num_columns, table_data[i], &status)) + { + fits_report_error (stderr, status); + return status; + } + } + + return 0; +} + + + + + + +int +make_spec (inroot) + char *inroot; +{ + + fitsfile *fptr; // Pointer to the FITS file + int status = 0; // CFITSIO status value MUST be initialized to zero + + // Create a new FITS file + if (fits_create_file (&fptr, "!multi_extension.fits", &status)) + { + fits_report_error (stderr, status); // Print any error message + return status; + } + +// Create an empty primary HDU + if (fits_create_img (fptr, FLOAT_IMG, 0, NULL, &status)) + { + fits_report_error (stderr, status); + } + + printf ("Hello World %s \n", inroot); + printf ("Plasma %d \n", NPLASMA); + printf ("NBINS in spec %d \n", NBINS_IN_CELL_SPEC); + + + Spectra spectra; + spectra.num_spectra = NPLASMA; + spectra.num_wavelengths = NBINS_IN_CELL_SPEC; + spectra.data = calloc (spectra.num_spectra, sizeof (float *)); + + for (int i = 0; i < spectra.num_spectra; i++) + { + spectra.data[i] = calloc (spectra.num_wavelengths, sizeof (float)); + + // Fill with sample data + for (int j = 0; j < spectra.num_wavelengths; j++) + { + // spectra.data[i][j] = (float) (i + j); + spectra.data[i][j] = (float) plasmamain[i].cell_spec_flux[j]; + } + } + + // Prepare image data for the first extension + float *image_data1 = prepare_image_data (spectra.data, spectra.num_wavelengths, spectra.num_spectra); + if (!image_data1) + { + // Free allocated memory before exiting + for (int i = 0; i < spectra.num_spectra; i++) + { + free (spectra.data[i]); + } + free (spectra.data); + return 1; + } + + + status = write_image_extension (fptr, image_data1, spectra.num_wavelengths, spectra.num_spectra, "Jnu"); + + Spectra freq; + freq.num_spectra = 1; + freq.num_wavelengths = NBINS_IN_CELL_SPEC; + freq.data = calloc (freq.num_spectra, sizeof (float *)); + freq.data[0] = calloc (freq.num_wavelengths, sizeof (float)); + + for (int i = 0; i < NBINS_IN_CELL_SPEC; i++) + { + freq.data[0][i] = (float) pow (10., geo.cell_log_freq_min + i * geo.cell_delta_lfreq); + } + + + float *image_data2 = prepare_image_data (freq.data, NBINS_IN_CELL_SPEC, 1); +/* + double freq[2000][1]; + int i; + + for (i = 0; i < NBINS_IN_CELL_SPEC; i++) + { + freq[i][0] = (float) pow (10., geo.cell_log_freq_min + i * geo.cell_delta_lfreq); + } + + float *image_data2 = prepare_image_data (freq, NBINS_IN_CELL_SPEC, 1); +*/ + status = write_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, 1, "nu"); + + if (fits_close_file (fptr, &status)) + { + fits_report_error (stderr, status); + } + + return 0; +} + /**********************************************************/ /** * @brief parses the command line options @@ -178,6 +405,7 @@ xparse_command_line (argc, argv) **********************************************************/ + int main (argc, argv) int argc; @@ -207,6 +435,8 @@ main (argc, argv) wind_read (infile); + make_spec (inroot); + if (nlevels_macro == 0) { printf ("Currently this routine only looks at macro atom values, and this is a simple atom file\n"); From 9ba8bf79146f5d22288d3fa862cae5d7bda51007 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Sat, 10 Aug 2024 08:14:12 -0500 Subject: [PATCH 03/23] Remove extra code --- source/windsave2fits.c | 364 +---------------------------------------- 1 file changed, 1 insertion(+), 363 deletions(-) diff --git a/source/windsave2fits.c b/source/windsave2fits.c index 88d5dfda2..1e509680e 100644 --- a/source/windsave2fits.c +++ b/source/windsave2fits.c @@ -413,9 +413,7 @@ main (argc, argv) { char infile[LINELENGTH], outfile[LINELENGTH]; - int n, i; - FILE *fptr, *fopen (); - int ii, jj, ndom, nnwind; + FILE *fopen (); int mkdir (); @@ -437,365 +435,5 @@ main (argc, argv) make_spec (inroot); - if (nlevels_macro == 0) - { - printf ("Currently this routine only looks at macro atom values, and this is a simple atom file\n"); - exit (0); - } - - - - fptr = fopen (outfile, "w"); - - fprintf (fptr, "# Results for %s\n", infile); - - fprintf (fptr, "# Size %d %d %d \n", size_Jbar_est, size_gamma_est, size_alpha_est); - - fprintf (fptr, "%15s %4s %4s %4s ", "Variable", "np", "i", "j"); - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "MacLev%02d ", i); - } - fprintf (fptr, "\n"); - - fprintf (fptr, "%15s %4s %4s %4s ", "---------------", "----", "----", "----"); - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "-------- "); - } - fprintf (fptr, "\n"); - - for (n = 0; n < NPLASMA; n++) - { - - fprintf (fptr, "%15s", "jbar"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].jbar[i]); - } - fprintf (fptr, "\n"); - } - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "jbar_old"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].jbar_old[i]); - } - fprintf (fptr, "\n"); - } - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "gamma"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].gamma[i]); - } - fprintf (fptr, "\n"); - } - - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "alpha_st"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].alpha_st[i]); - } - fprintf (fptr, "\n"); - } - - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "recomb_sp"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].recomb_sp[i]); - } - fprintf (fptr, "\n"); - } - - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "matom_abs"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].matom_abs[i]); - } - fprintf (fptr, "\n"); - } - - - - - for (n = 0; n < NPLASMA; n++) - { - fprintf (fptr, "%15s", "matom_emiss"); - - nnwind = plasmamain[n].nwind; - ndom = wmain[nnwind].ndom; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, " %4d %4d %4d ", n, ii, jj); - - for (i = 0; i < nlevels_macro; i++) - { - fprintf (fptr, "%8.2e ", macromain[n].matom_emiss[i]); - } - fprintf (fptr, "\n"); - } - fclose (fptr); - - - /* calculate all the line luminosities in macro-atom mode */ - /* these are stored in a folder with a separate file for each upper level */ - sprintf (folder, "matom_linelum_%.200s", inroot); - mkdir (folder, 0777); - printf ("Calculating macro-atom line luminosities for all lines in wavelength range %.1f to %.1f Angstroms...\n", - VLIGHT / geo.sfmax / ANGSTROM, VLIGHT / geo.sfmin / ANGSTROM); - - /* create a file that tells the user the wavelengths of every macro-atom line and the - corresponding upper and lower levels */ - create_matom_level_map (); - - /* loop over all macro atom upper levels */ - for (i = 0; i < nlevels_macro; i++) - { - line_matom_lum (i); - } - printf ("Done for %d levels.\n", nlevels_macro); - - exit (0); - -} - -/**********************************************************/ -/** - * @brief print out which line wavelengths correspond to which upper and lower levels - * - **********************************************************/ - -int -create_matom_level_map () -{ - int uplvl, nbbd, n; - char outfile[LINELENGTH]; - FILE *fptr, *fopen (); - - /* open a file in the folder where we store the matom line luminosities */ - sprintf (outfile, "%.200s/line_map.txt", folder); - - /* print some header information to the file */ - fptr = fopen (outfile, "w"); - fprintf (fptr, "# model name %s\n", inroot); - fprintf (fptr, "# which line wavelengths (Angstroms) correspond to which upper and lower levels\n"); - fprintf (fptr, "upper lower wavelength\n"); - - for (uplvl = 0; uplvl < nlevels_macro; uplvl++) - { - nbbd = xconfig[uplvl].n_bbd_jump; - for (n = 0; n < nbbd; n++) - { - fprintf (fptr, "%d %d %12.2f\n", uplvl, line[xconfig[uplvl].bbd_jump[n]].nconfigl, - VLIGHT / line[xconfig[uplvl].bbd_jump[n]].freq / ANGSTROM); - } - } - fclose (fptr); return (0); } - - -/**********************************************************/ -/** - * @brief calculate line luminosities for all cells for a given upper level and save to file - * - **********************************************************/ - -int -line_matom_lum (uplvl) - int uplvl; -{ - int n, nbbd, i, ii, jj, nnwind, ndom, inwind; - double lum[NBBJUMPS]; - char outfile[LINELENGTH]; - FILE *fptr, *fopen (); - - nbbd = xconfig[uplvl].n_bbd_jump; - - /* open a file in the folder where we store the matom line luminosities */ - sprintf (outfile, "%.100s/linelums_%.100s_lev%d.txt", folder, inroot, uplvl); - - /* print some header information to the file */ - fptr = fopen (outfile, "w"); - fprintf (fptr, "# model name %s\n", inroot); - fprintf (fptr, "# Line luminosities from upper level %d from %.2f to %.2f Angstroms\n", uplvl, VLIGHT / geo.sfmax / ANGSTROM, - VLIGHT / geo.sfmin / ANGSTROM); - fprintf (fptr, "# Lower Levels: "); - for (n = 0; n < nbbd; n++) - { - fprintf (fptr, "%12d", line[xconfig[uplvl].bbd_jump[n]].nconfigl); - } - fprintf (fptr, "\n# Line Wavelengths:"); - for (n = 0; n < nbbd; n++) - { - fprintf (fptr, " %12.2f", VLIGHT / line[xconfig[uplvl].bbd_jump[n]].freq / ANGSTROM); - } - fprintf (fptr, "\n"); - - /* columns of wind cells */ - fprintf (fptr, "%12s %12s %12s %4s %4s %12s %12s %12s", "nwind", "x", "z", "i", "j", "inwind", "wind_vol", "plasma_vol"); - - /* header info for each lower level */ - for (n = 0; n < nbbd; n++) - { - fprintf (fptr, " LowerLev%03d ", line[xconfig[uplvl].bbd_jump[n]].nconfigl); - } - fprintf (fptr, "\n"); - // fprintf (fptr, "%12s %12s %12s %4s %4s %12s %12s %12s", "------------", "------------", "------------", "----", "----", "------------", - // "------------", "------------"); - // for (n = 0; n < nbbd; n++) - // { - // fprintf (fptr, " ----------- "); - // } - // fprintf (fptr, "\n"); - - - /* loop over each plasma cell and do the calculation */ - for (nnwind = 0; nnwind < NDIM2; nnwind++) - { - ndom = wmain[nnwind].ndom; - inwind = wmain[nnwind].inwind; - wind_n_to_ij (ndom, nnwind, &ii, &jj); - fprintf (fptr, "%12d %12.4e %12.4e %4d %4d %12d %12.4e", nnwind, wmain[nnwind].xcen[0], wmain[nnwind].xcen[2], ii, jj, inwind, - wmain[nnwind].vol); - if (wmain[nnwind].inwind >= 0) - { - n = wmain[nnwind].nplasma; - line_matom_lum_single (lum, &plasmamain[n], uplvl); - /* print the filled volume */ - fprintf (fptr, " %13.4e", plasmamain[n].vol); - for (i = 0; i < nbbd; i++) - fprintf (fptr, " %13.4e", lum[i]); - } - else - { - /* print a 0.0 instead of the filled volume and line lums outside the wind */ - fprintf (fptr, " %13.4e", 0.0); - for (i = 0; i < nbbd; i++) - fprintf (fptr, " %13.4e", 0.0); - } - - fprintf (fptr, "\n"); - } - fclose (fptr); - return (0); -} - -/**********************************************************/ -/** - * @brief calculate line luminosities for a single cell for a given upper level - * - **********************************************************/ - -double -line_matom_lum_single (lum, xplasma, uplvl) - double lum[]; - PlasmaPtr xplasma; - int uplvl; -{ - int n, nbbd, m; - double penorm, bb_cont; - double freq_min, freq_max, lum_tot; - struct lines *line_ptr; - double eprbs[NBBJUMPS]; - freq_min = geo.sfmin; - freq_max = geo.sfmax; - /* identify number of bb downward jumps */ - nbbd = xconfig[uplvl].n_bbd_jump; - penorm = 0.0; - m = 0; - lum_tot = 0.0; - /* work out how often we come out in each of these locations */ - for (n = 0; n < nbbd; n++) - { - line_ptr = &line[xconfig[uplvl].bbd_jump[n]]; - if ((line_ptr->freq > freq_min) && (line_ptr->freq < freq_max)) // correct range - { - // bb_cont = (a21 (line_ptr) * p_escape (line_ptr, xplasma)); - bb_cont = (a21 (line_ptr) * 1.0); - eprbs[m] = bb_cont * (xconfig[uplvl].ex - xconfig[line[xconfig[uplvl].bbd_jump[n]].nconfigl].ex); //energy difference - penorm += eprbs[m]; - } - else - { - eprbs[m] = 0.0; - } - m++; - } - - /* correctly normalise the probabilities */ - for (n = 0; n < nbbd; n++) - { - if (penorm == 0) - { - lum[n] = 0.0; - } - else - { - eprbs[n] = eprbs[n] / penorm; - lum[n] = eprbs[n] * macromain[xplasma->nplasma].matom_emiss[uplvl]; - } - lum_tot += lum[n]; - } - - return (lum_tot); -} From ee5284710d75e0d80bc608cb3a2da255f0adb69e Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Sat, 10 Aug 2024 13:55:53 -0500 Subject: [PATCH 04/23] Info to evaluate cellspec models now included --- source/windsave2fits.c | 324 ++++++++++++++++++++++++++++++++++------- 1 file changed, 274 insertions(+), 50 deletions(-) diff --git a/source/windsave2fits.c b/source/windsave2fits.c index 1e509680e..e98562b87 100644 --- a/source/windsave2fits.c +++ b/source/windsave2fits.c @@ -2,19 +2,18 @@ /***********************************************************/ /** @file windsave2fits.c * @author ksl - * @date October, 2021 + * @date August,2024 * - * @brief Routines to inspect variables in a wind structure + * @brief Routines to write variious portions of in a wind structure + * to a fits file * *###Notes### - * This is intended just as a diagnostic routine - * so that one can print out whatever variables in - * a windstrucutre one wants in order to diagnose - * a problem. It was written so that we could inspect - * some of the macro atom variables in paralell mode - * in diagnosing issue #898 and #910, but anyon - * should change it so that other problems might - * be addressed. + * This routine was written primarily to better understand + * how well the spectra we use for estimating ionization rates + * model the actual cell spectra; but it can be rather + * straightforwardly adatpted to create images or table + * of other portions of a windsave file, that involve + * large amounts of data. * * ***********************************************************/ @@ -29,6 +28,9 @@ #include "fitsio.h" +/* The first portions of the program contain some utility routines + for working with cfitsio */ + char inroot[LINELENGTH], outroot[LINELENGTH], model_file[LINELENGTH], folder[LINELENGTH]; int model_flag, ksl_flag, cmf2obs_flag, obs2cmf_flag; @@ -44,9 +46,18 @@ typedef struct float **data; // 2D array to hold spectra data } Spectra; -/* Next fucntions are utilities that it should be possitle to call repeatedly*/ +/* Next functions are utilities that it should be possitle to call repeatedly*/ + +/**********************************************************/ +/** + * @brief Function to perpare 2d image data from a 2d array + * + * @details The routine simply returns a flattened version of the 2d array + * ### Notes ### + * + * + **********************************************************/ -// Function to prepare 2D image data from a 2D array float * prepare_image_data (float **data, int width, int height) { @@ -70,31 +81,21 @@ prepare_image_data (float **data, int width, int height) return image_data; } -// Function to prepare binary table data -void -prepare_binary_table (int num_rows, int num_columns, float **table_data) -{ - // Allocate memory for each row of the table - for (int i = 0; i < num_rows; i++) - { - table_data[i] = malloc (num_columns * sizeof (float)); - if (!table_data[i]) - { - fprintf (stderr, "Memory allocation error\n"); - return; - } - - // Fill each row with some sample data - for (int j = 0; j < num_columns; j++) - { - table_data[i][j] = (float) (i * num_columns + j); // Example data - } - } -} +/**********************************************************/ +/** + * @brief Function to write a 2d image data to a fits file with + * a ginve extension name + * + * @details The routine normally returns 0 + * ### Notes ### + * This routine is called after prepare_image data, which will + * have flattened and converted the original array to a float + * + * + **********************************************************/ -// Function to write a 2D image to a FITS extension with a given name int write_image_extension (fitsfile *fptr, float *image_data, int width, int height, const char *extname) { @@ -129,14 +130,63 @@ write_image_extension (fitsfile *fptr, float *image_data, int width, int height, return 0; } -// Function to write a binary table to a FITS extension with a given name + + + +/**********************************************************/ +/** + * @brief Function to write a 1d image extension with + * a ginve extension name + * + * @details The routine normally returns 0 + * ### Notes ### + * + * + **********************************************************/ + int -write_binary_table_extension (fitsfile *fptr, int num_rows, int num_columns, float **table_data, const char *extname) +write_1d_image_extension (fitsfile *fptr, float *data, long length, const char *extname) +{ + int status = 0; // CFITSIO status value MUST be initialized to zero + long naxes[1]; // Array to store the dimension size (1D array) + + naxes[0] = length; // Length of the 1D array + + // Create a new image extension of type FLOAT_IMG + if (fits_create_img (fptr, FLOAT_IMG, 1, naxes, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Set the extension name using the EXTNAME keyword + if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) extname, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Write the 1D array data to the FITS image extension + if (fits_write_img (fptr, TFLOAT, 1, length, data, &status)) + { + fits_report_error (stderr, status); + return status; + } + + return 0; +} + + +/* This routine is not longer used but has been kept for reference +// Function to write a binary table with mixed data types to a FITS extension +int +write_mixed_table_extension (fitsfile *fptr, int num_rows, int num_columns, int *int_data, float *float_data, const char *extname) { int status = 0; - char *ttype[] = { "COL1" }; // Name for each column - char *tform[] = { "E" }; // Format for each column (E = float) - char *tunit[] = { "unit" }; // Unit for each column + // Define the names, formats, and units for each column + char *ttype[] = { "INT_COL", "FLOAT_COL" }; // Column names + char *tform[] = { "J", "E" }; // Formats: 'J' for integer, 'E' for float + char *tunit[] = { "", "" }; // Units // Create a new binary table extension if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_columns, ttype, tform, tunit, NULL, &status)) @@ -152,22 +202,173 @@ write_binary_table_extension (fitsfile *fptr, int num_rows, int num_columns, flo return status; } - // Write the table data - for (int i = 0; i < num_rows; i++) + // Write the integer data to the first column + if (fits_write_col (fptr, TINT, 1, 1, 1, num_rows, int_data, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Write the float data to the second column + if (fits_write_col (fptr, TFLOAT, 2, 1, 1, num_rows, float_data, &status)) + { + fits_report_error (stderr, status); + return status; + } + + return 0; +} +*/ + +/* This ends the utility routines */ + +/**********************************************************/ +/** + * @brief write_spectral_model_table extension to the + * fits file + * + * @details The routine normally returns 0 + * ### Notes ### + * + * + **********************************************************/ + +int +write_spectra_model_table (fitsfile *fptr) +{ + + /* We need one row for each cell and band so + our table will be NPLASMA*nbands long + */ + + int i, j, k; + int num_rows, num_cols; + int nbands; + int *ichoice, *iplasma; + float *exp_w, *exp_temp, *pl_log_w, *pl_alpha; + + nbands = geo.nxfreq; + num_rows = nbands * NPLASMA; + printf ("xtest %d %d\n", nbands, NPLASMA); + num_cols = 6; // for now + char *table_name = "spec_model"; + + ichoice = calloc (num_rows, sizeof (int *)); + iplasma = calloc (num_rows, sizeof (int *)); + exp_w = calloc (num_rows, sizeof (float *)); + exp_temp = calloc (num_rows, sizeof (float *)); + pl_log_w = calloc (num_rows, sizeof (float *)); + pl_alpha = calloc (num_rows, sizeof (float *)); + + k = 0; + for (i = 0; i < NPLASMA; i++) { - if (fits_write_col (fptr, TFLOAT, 1, i + 1, 1, num_columns, table_data[i], &status)) + for (j = 0; j < nbands; j++) { - fits_report_error (stderr, status); - return status; + iplasma[k] = i; + ichoice[k] = plasmamain[i].spec_mod_type[j]; + exp_w[k] = plasmamain[i].exp_w[j]; + exp_temp[k] = plasmamain[i].exp_temp[j]; + pl_log_w[k] = plasmamain[i].pl_log_w[j]; + pl_alpha[k] = plasmamain[i].pl_alpha[j]; + k++; } } + + int status = 0; + // Define the names, formats, and units for each column + char *ttype[] = { "nplasma", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha" }; // Column names + char *tform[] = { "J", "J", "E", "E", "E", "E" }; // Formats: 'J' for integer, 'E' for float + char *tunit[] = { "", "", "", "", "", "" }; // Units + + // Create a new binary table extension + if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_cols, ttype, tform, tunit, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Set the extension name using the EXTNAME keyword + if (fits_update_key (fptr, TSTRING, "EXTNAME", (void *) table_name, NULL, &status)) + { + fits_report_error (stderr, status); + return status; + } + + printf ("num_rows: %d\n", num_rows); + for (i = 0; i < num_rows; i++) + { + printf ("%d %d %d\n", i, iplasma[i], ichoice[i]); + } + + // Write the integer data to the first column + if (fits_write_col (fptr, TINT, 1, 1, 1, num_rows, iplasma, &status)) + { + fits_report_error (stderr, status); + return status; + } + + // Write the float data to the second column + if (fits_write_col (fptr, TINT, 2, 1, 1, num_rows, ichoice, &status)) + { + fits_report_error (stderr, status); + return status; + } + + + // Write the float data to the third column + if (fits_write_col (fptr, TFLOAT, 3, 1, 1, num_rows, exp_w, &status)) + { + fits_report_error (stderr, status); + return status; + } + + + + // Write the float data to the fourth column + if (fits_write_col (fptr, TFLOAT, 4, 1, 1, num_rows, exp_temp, &status)) + { + fits_report_error (stderr, status); + return status; + } + + + // Write the float data to the fifth column + if (fits_write_col (fptr, TFLOAT, 5, 1, 1, num_rows, pl_log_w, &status)) + { + fits_report_error (stderr, status); + return status; + } + + + // Write the float data to the sixth column + if (fits_write_col (fptr, TFLOAT, 6, 1, 1, num_rows, pl_alpha, &status)) + { + fits_report_error (stderr, status); + return status; + } + return 0; } +/**********************************************************/ +/** + * @brief write_a fits file that contains information + * releveant to how well the cell spectra are modelled + * fits file + * + * @details The routine normally returns 0 + * ### Notes ### + * + * This is the main routine, in the sense that everyting is + * controlled frrom here + * + **********************************************************/ + int @@ -177,9 +378,16 @@ make_spec (inroot) fitsfile *fptr; // Pointer to the FITS file int status = 0; // CFITSIO status value MUST be initialized to zero + char outfile[LINELENGTH]; + + sprintf (outfile, "%s_cellspec.fits", inroot); + printf ("outfile %s\n", outfile); + + // Create a new FITS file - if (fits_create_file (&fptr, "!multi_extension.fits", &status)) + //if (fits_create_file (&fptr, "!multi_extension.fits", &status)) + if (fits_create_file (&fptr, outfile, &status)) { fits_report_error (stderr, status); // Print any error message return status; @@ -253,7 +461,22 @@ make_spec (inroot) float *image_data2 = prepare_image_data (freq, NBINS_IN_CELL_SPEC, 1); */ - status = write_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, 1, "nu"); + // status = write_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, 1, "nu"); + status = write_1d_image_extension (fptr, image_data2, NBINS_IN_CELL_SPEC, "nu"); + + /* Elimainate this for now + + // Prepare data for the binary table + int num_rows = 5; + int int_data[] = { 1, 2, 3, 4, 5 }; // Example integer data + float float_data[] = { 0.1, 0.2, 0.3, 0.4, 0.5 }; // Example float data + + write_mixed_table_extension (fptr, num_rows, 2, int_data, float_data, "Table_Ext"); + */ + +/* Now write my table */ + + write_spectra_model_table (fptr); if (fits_close_file (fptr, &status)) { @@ -280,6 +503,10 @@ make_spec (inroot) * Although this routine uses the standard Log and Error commands * the diag files have not been created yet and so this information * is really simply written to the terminal. + * + * WARNING: this has not been updeated. Currently there + * routine expects only the rootname of a windsave file + * * **********************************************************/ @@ -396,10 +623,7 @@ xparse_command_line (argc, argv) * This routine oversees the effort. The basic steps are * * - parse the command line to get the names of files - * - read the old windsave file - * - read the densities from in this case H - * - modify the densities - * - write out the new windsave file + * - call the routine to craetate fits file * * **********************************************************/ From dfe7b3f5032f59f5f4e781e5290ad5541f12fb8b Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Sun, 11 Aug 2024 14:10:56 -0500 Subject: [PATCH 05/23] More updtates to windsave2fits --- source/windsave2fits.c | 56 +++++++++++++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/source/windsave2fits.c b/source/windsave2fits.c index e98562b87..9415de647 100644 --- a/source/windsave2fits.c +++ b/source/windsave2fits.c @@ -244,17 +244,18 @@ write_spectra_model_table (fitsfile *fptr) int i, j, k; int num_rows, num_cols; int nbands; - int *ichoice, *iplasma; + int *ichoice, *iplasma, *iband; float *exp_w, *exp_temp, *pl_log_w, *pl_alpha; nbands = geo.nxfreq; num_rows = nbands * NPLASMA; printf ("xtest %d %d\n", nbands, NPLASMA); - num_cols = 6; // for now + num_cols = 7; // for now char *table_name = "spec_model"; ichoice = calloc (num_rows, sizeof (int *)); iplasma = calloc (num_rows, sizeof (int *)); + iband = calloc (num_rows, sizeof (int *)); exp_w = calloc (num_rows, sizeof (float *)); exp_temp = calloc (num_rows, sizeof (float *)); pl_log_w = calloc (num_rows, sizeof (float *)); @@ -266,6 +267,7 @@ write_spectra_model_table (fitsfile *fptr) for (j = 0; j < nbands; j++) { iplasma[k] = i; + iband[k] = j; ichoice[k] = plasmamain[i].spec_mod_type[j]; exp_w[k] = plasmamain[i].exp_w[j]; exp_temp[k] = plasmamain[i].exp_temp[j]; @@ -278,9 +280,9 @@ write_spectra_model_table (fitsfile *fptr) int status = 0; // Define the names, formats, and units for each column - char *ttype[] = { "nplasma", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha" }; // Column names - char *tform[] = { "J", "J", "E", "E", "E", "E" }; // Formats: 'J' for integer, 'E' for float - char *tunit[] = { "", "", "", "", "", "" }; // Units + char *ttype[] = { "nplasma", "band", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha" }; // Column names + char *tform[] = { "J", "J", "J", "E", "E", "E", "E" }; // Formats: 'J' for integer, 'E' for float + char *tunit[] = { "", "", "", "", "", "", "" }; // Units // Create a new binary table extension if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_cols, ttype, tform, tunit, NULL, &status)) @@ -309,33 +311,33 @@ write_spectra_model_table (fitsfile *fptr) return status; } + // Write the float data to the second column - if (fits_write_col (fptr, TINT, 2, 1, 1, num_rows, ichoice, &status)) + if (fits_write_col (fptr, TINT, 2, 1, 1, num_rows, iband, &status)) { fits_report_error (stderr, status); return status; } - // Write the float data to the third column - if (fits_write_col (fptr, TFLOAT, 3, 1, 1, num_rows, exp_w, &status)) + if (fits_write_col (fptr, TINT, 3, 1, 1, num_rows, ichoice, &status)) { fits_report_error (stderr, status); return status; } - - // Write the float data to the fourth column - if (fits_write_col (fptr, TFLOAT, 4, 1, 1, num_rows, exp_temp, &status)) + // Write the float data to the fourth column + if (fits_write_col (fptr, TFLOAT, 4, 1, 1, num_rows, exp_w, &status)) { fits_report_error (stderr, status); return status; } - // Write the float data to the fifth column - if (fits_write_col (fptr, TFLOAT, 5, 1, 1, num_rows, pl_log_w, &status)) + + // Write the float data to the fifth columnmn + if (fits_write_col (fptr, TFLOAT, 5, 1, 1, num_rows, exp_temp, &status)) { fits_report_error (stderr, status); return status; @@ -343,7 +345,15 @@ write_spectra_model_table (fitsfile *fptr) // Write the float data to the sixth column - if (fits_write_col (fptr, TFLOAT, 6, 1, 1, num_rows, pl_alpha, &status)) + if (fits_write_col (fptr, TFLOAT, 6, 1, 1, num_rows, pl_log_w, &status)) + { + fits_report_error (stderr, status); + return status; + } + + + // Write the float data to the seventh column + if (fits_write_col (fptr, TFLOAT, 7, 1, 1, num_rows, pl_alpha, &status)) { fits_report_error (stderr, status); return status; @@ -380,7 +390,10 @@ make_spec (inroot) int status = 0; // CFITSIO status value MUST be initialized to zero char outfile[LINELENGTH]; - sprintf (outfile, "%s_cellspec.fits", inroot); + /* The obscure exclamation mark at the beginning of the name means yu want to overwrie + and existing file, if it exists */ + + sprintf (outfile, "!%s_cellspec.fits", inroot); printf ("outfile %s\n", outfile); @@ -474,6 +487,19 @@ make_spec (inroot) write_mixed_table_extension (fptr, num_rows, 2, int_data, float_data, "Table_Ext"); */ + +/* Now create a 1d extension hold the frequencies of the broad bands */ + + int i; + float xfreq[100]; + + for (i = 0; i <= geo.nxfreq; i++) + { + xfreq[i] = geo.xfreq[i]; + } + + status = write_1d_image_extension (fptr, xfreq, geo.nxfreq + 1, "nu_model"); + /* Now write my table */ write_spectra_model_table (fptr); From 68b4d54dc16ea1c1481c0490ddf690b905299678 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Tue, 3 Sep 2024 15:53:49 -0500 Subject: [PATCH 06/23] Add better rad temp estimate --- source/Makefile | 2 +- source/estimators_simple.c | 84 +++++++++++++++++++++++++++++++++++++- 2 files changed, 84 insertions(+), 2 deletions(-) diff --git a/source/Makefile b/source/Makefile index 7afdde127..26395ad47 100644 --- a/source/Makefile +++ b/source/Makefile @@ -147,7 +147,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas $(CUDA_LIBS) #Note that version should be a single string without spaces. -VERSION = 88c +VERSION = 88d diff --git a/source/estimators_simple.c b/source/estimators_simple.c index 27359ba05..1e309030d 100644 --- a/source/estimators_simple.c +++ b/source/estimators_simple.c @@ -16,6 +16,7 @@ #include "atomic.h" #include "python.h" +#include /**********************************************************/ @@ -471,6 +472,83 @@ update_force_estimators (xplasma, p, phot_mid, ds, w_ave, ndom, z, frac_ff, frac } +typedef struct +{ + double T; // Temperature in Kelvin +} planck_params; + +// Planck spectral radiance function +double +planck_spectral_radiance (double nu, void *params) +{ + planck_params *p = (planck_params *) params; + double T = p->T; + double exponent = PLANCK * nu / (BOLTZMANN * T); + return (2 * PLANCK * pow (nu, 3) / pow (VLIGHT, 2)) / (exp (exponent) - 1); +} + +// Function to calculate nu * J(nu) +double +nu_times_radiance (double nu, void *params) +{ + return nu * planck_spectral_radiance (nu, params); +} + +// Function to calculate mean frequency over a specified range +double +mean_frequency_nu_range (double T, double nu_min, double nu_max) +{ + gsl_integration_workspace *w = gsl_integration_workspace_alloc (1000); + + double numerator, denominator, error; + planck_params params = { T }; + + gsl_function F; + F.function = &nu_times_radiance; + F.params = ¶ms; + + gsl_integration_qags (&F, nu_min, nu_max, 0, 1e-7, 1000, w, &numerator, &error); + + F.function = &planck_spectral_radiance; + + gsl_integration_qags (&F, nu_min, nu_max, 0, 1e-7, 1000, w, &denominator, &error); + + gsl_integration_workspace_free (w); + + return numerator / denominator; +} + +// Function to estimate temperature given a mean frequency and a frequency range +double +estimate_temperature_from_mean_frequency (double mean_nu_target, double nu_min, double nu_max, double initial_guess) +{ + double T_guess = initial_guess; + double tol = 1e-6; // Convergence tolerance + int max_iter = 100; // Maximum number of iterations + double mean_nu; + + + for (int i = 0; i < max_iter; i++) + { + // Calculate the mean frequency for the current temperature guess + mean_nu = mean_frequency_nu_range (T_guess, nu_min, nu_max); + + // Adjust the temperature guess + double T_new = T_guess * (mean_nu_target / mean_nu); + + // Check for convergence + if (fabs (T_new - T_guess) < tol) + { + return T_new; + } + + T_guess = T_new; + } + + return T_guess; // Return the estimated temperature +} + + /**********************************************************/ /** * @brief Normalise simple estimators and cell spectra @@ -527,7 +605,11 @@ normalise_simple_estimators (xplasma) xplasma->j_scatt /= (4. * PI * invariant_volume_time); xplasma->t_r_old = xplasma->t_r; // Store the previous t_r in t_r_old immediately before recalculating - radiation_temperature = xplasma->t_r = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); + + radiation_temperature = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); + radiation_temperature = xplasma->t_r = + estimate_temperature_from_mean_frequency (xplasma->ave_freq, xband.f1[0], xband.f2[xband.nbands - 1], radiation_temperature); + xplasma->w = PI * xplasma->j / (STEFAN_BOLTZMANN * radiation_temperature * radiation_temperature * radiation_temperature * radiation_temperature); From 46136758ae205c4baebd3e67fac1e26c69949e2e Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Thu, 12 Sep 2024 17:14:06 -0500 Subject: [PATCH 07/23] Update to 88e --- source/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/Makefile b/source/Makefile index 26395ad47..2fec4bfbb 100644 --- a/source/Makefile +++ b/source/Makefile @@ -147,7 +147,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas $(CUDA_LIBS) #Note that version should be a single string without spaces. -VERSION = 88d +VERSION = 88e From 348eb438c0887e1c968685ea386caa9de686aadb Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Thu, 12 Sep 2024 17:31:07 -0500 Subject: [PATCH 08/23] Add hook for new mode Only a place holder for a new mode has been added at this point. --- source/python.h | 4 ++++ source/setup.c | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/source/python.h b/source/python.h index 41aac6daf..8c54b88b4 100644 --- a/source/python.h +++ b/source/python.h @@ -1223,6 +1223,10 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; #define IONMODE_MATRIX_BB 8 /**< matrix solver BB model */ #define IONMODE_MATRIX_SPECTRALMODEL 9 /**< matrix solver spectral model based on power laws */ #define IONMODE_MATRIX_ESTIMATORS 10 /**< matrix solver spectral model based on power laws */ +#define IONMODE_MATRIX_MULTISHOT 11 /**< matrix solver spectral model based on power laws which + * updates T_e multiple times before arriving at a final + * solution + */ // and the corresponding modes in nebular_concentrations #define NEBULARMODE_TR 0 /**< LTE using t_r */ diff --git a/source/setup.c b/source/setup.c index 0cd7fec2e..6d786e122 100644 --- a/source/setup.c +++ b/source/setup.c @@ -714,7 +714,8 @@ init_ionization () strcpy (answer, "matrix_bb"); geo.ioniz_mode = - rdchoice ("Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est)", "0,3,1,4,2,8,9,10", answer); + rdchoice ("Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est,matrix_te)", "0,3,1,4,2,8,9,10,11", + answer); if (geo.ioniz_mode == IONMODE_FIXED) { From b23718d6331fdd76475472e8d48f88d9d5b285b5 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Thu, 12 Sep 2024 18:08:00 -0500 Subject: [PATCH 09/23] Added dummy routine for new ionization mode Basically, at this point matrix_populations2 is identical to matrix_populations, so this is just a routine to be used to develop a new ionization mode --- source/Makefile | 2 +- source/matrix_ion2.c | 416 +++++++++++++++++++++++++++++++++++++++++++ source/python.h | 5 + source/saha.c | 4 + source/templates.h | 6 + 5 files changed, 432 insertions(+), 1 deletion(-) create mode 100644 source/matrix_ion2.c diff --git a/source/Makefile b/source/Makefile index 2fec4bfbb..7fa3ce367 100644 --- a/source/Makefile +++ b/source/Makefile @@ -201,7 +201,7 @@ python_source = agn.c anisowind.c atomic_extern_init.c atomicdata.c atomicdata_i extract.c frame.c gradv.c gridwind.c homologous.c hydro_import.c import.c \ import_calloc.c import_cylindrical.c import_rtheta.c import_spherical.c ionization.c \ janitor.c knigge.c levels.c lines.c macro_accelerate.c macro_gen_f.c macro_gov.c \ - matom.c matom_diag.c matrix_cpu.c matrix_ion.c models_extern_init.c para_update.c \ + matom.c matom_diag.c matrix_cpu.c matrix_ion.c matrix_ion2.c models_extern_init.c para_update.c \ parse.c partition.c paths.c phot_util.c photon2d.c photon_gen.c photon_gen_matom.c \ pi_rates.c python_extern_init.c radiation.c random.c rdpar.c rdpar_init.c recipes.c \ recomb.c resonate.c reverb.c roche.c rtheta.c run.c saha.c setup.c setup_disk.c setup_domains.c \ diff --git a/source/matrix_ion2.c b/source/matrix_ion2.c new file mode 100644 index 000000000..2caf0c09f --- /dev/null +++ b/source/matrix_ion2.c @@ -0,0 +1,416 @@ + +/***********************************************************/ +/** @file matrix_ion.c + * @author ksl + * @date January, 2018 + * + * @brief Contains routines which compute the ionization state using a matrix inversion technique + * + ***********************************************************/ + +#include +#include +#include +#include + +#include +#include "atomic.h" +#include "python.h" + + +/**********************************************************/ +/** + * @brief A matrix solver for the ionization state in a cell + * + * @param [in,out] PlasmaPtr xplasma The plasma cell we are working on + * @param [in] int mode What type of model to use for J_nu - 1=fitted model 2=blackbody + * @return 0 if successful + * + * The abundances contained in xplasma are updated witht he results of + * the calculation. + * + * @details + * modes: + * 1 - use a power law and or exponential model for J + * 2 - use a dilute BB model for J defined by t_r and w + * The general structure of matrix_ion_populations is as follows: + * First we compute all of the rates that we have data for linking all ionization stages + * We make an initial guess at the electron density + * We attempt to solve the ionization rates matrix, and produce a new guess at the electron density + * We then calculate new rates, re-solve and proceed until the electron density converges. + * We solve the matrix equation A x = b for the vector x, where A is a square rate matrix linking + * each state with each other state (often containing zeros if the ions arent linked) b is a vector + * containing the total density of each element arranged in a certain way (see later) and + * x is a vector containing the relative ion populations. We invert A to get x=bA^-1 + * + * ### Notes ### + * Uses a relative abundance scheme, in order to reduce large number issues + * + * Various parameters for the calculation, and in particular the t_e are passed + * via the PlasmaPtr + * + **********************************************************/ + +int +matrix_ion_populations2 (xplasma, mode) + PlasmaPtr xplasma; + int mode; + +{ + double elem_dens[NELEMENTS]; //The density of each element + int nn, mm, nrows; + double rate_matrix[nions][nions]; //The rate matrix that we are going to try and solve + double newden[NIONS]; //A temporary array to hold our intermediate solutions + double nh, nh1, nh2, t_e; + double xne, xxne, xxxne; //Various stores for intermediate guesses at electron density + double b_temp[nions]; //The b matrix + double *b_data, *a_data; //These arrays are allocated later and sent to the matrix solver + double *populations; //This array is allocated later and is retrieved from the matrix solver + int matrix_err, niterate; //counters for errors and the number of iterations we have tried to get a converged electron density + double xnew; + int xion[nions]; // This array keeps track of what ion is in each line +//OLD int xelem[nions]; // This array keeps track of the element for each ion + double pi_rates[nions]; //photoionization rate coefficients + double rr_rates[nions]; //radiative recombination rate coefficients + double inner_rates[n_inner_tot]; //This array contains the rates for each of the inner shells. Where they go to requires the electron yield array + + nh1 = nh2 = 0; + + /* Copy some quantities from the cell into local variables */ + + nh = xplasma->rho * rho2nh; // The number density of hydrogen ions - computed from density + t_e = xplasma->t_e; // The electron temperature in the cell - used for collisional processes + + + /* We now calculate the total abundances for each element to allow us to use fractional abundances */ + + /* Belt and braces, we zero our array */ + for (mm = 0; mm < ion[nions - 1].z + 1; mm++) + { + elem_dens[mm] = 0.0; + } + + /* Now we populate the elemental abundance array */ + for (mm = 0; mm < nions; mm++) + { + elem_dens[ion[mm].z] = elem_dens[ion[mm].z] + xplasma->density[mm]; + } + + /* Dielectronic recombination, collisional ionization coefficients, three body recombination and + charge_exchange rate coefficients depend only on electron temperature, calculate them now - + they will not change they are all stored in global arrays */ + + compute_dr_coeffs (t_e); + compute_di_coeffs (t_e); + compute_qrecomb_coeffs (t_e); + compute_ch_ex_coeffs (t_e); + + /* In the following loop, over all ions in the simulation, we compute the radiative recombination rates, and photionization + rates OUT OF each ionization stage. The PI rates are calculated either using the modelled mean intensity in a cell, or + using the dilute blackbody approximation, depending on which mode we are in. At the same time, we copy the ion densities + from the plasma structure into a local array. We will only overwrite the numbers in the structure if we believe the + results are an improvement on what is there. */ + + for (mm = 0; mm < nions; mm++) + { + newden[mm] = xplasma->density[mm] / elem_dens[ion[mm].z]; // newden is our local fractional density array + xion[mm] = mm; // xion is an array we use to track which ion is in which row of the matrix + if (mm != ele[ion[mm].nelem].firstion) // We can recombine since we are not in the first ionization stage + { + rr_rates[mm] = total_rrate (mm, xplasma->t_e); // radiative recombination rates + } + if (ion[mm].istate != ele[ion[mm].nelem].istate_max) // we can photoionize, since we are not in the highest ionization state + { + if (mode == NEBULARMODE_MATRIX_BB) + { + pi_rates[mm] = calc_pi_rate (mm, xplasma, 2, 1); // PI rate for the BB model + } + else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL) + { + pi_rates[mm] = calc_pi_rate (mm, xplasma, 1, 1); // PI rate for an explicit spectral model + } + else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) + { + pi_rates[mm] = xplasma->ioniz[mm] / xplasma->density[mm]; // PI rate logged during the photon passage + } + else + { + Error ("matrix_ion_populations: Unknown mode %d\n", mode); + Exit (0); + } + } + +/* ksl: This whole lope is not actually used, and produces warnings in gcc11 */ +//OLD for (nn = 0; nn < nelements; nn++) +//OLD { +//OLD if (ion[mm].z == ele[nn].z) +//OLD { +//OLD xelem[mm] = nn; /* xelem logs which element each row in the arrays refers to. This is important because we need +//OLD to know what the total density will be for a group of rows all representing the same +//OLD element. */ +//OLD } +//OLD } + } + + + /* The next loop generates the inner shell ionization rates, if they are present in the atomic data and + we wish to compute auger ionizaion rates. This only computes the rates out of each ion, we also need to + consult the electron yield array if we are to compute the change in electron number */ + + if (geo.auger_ionization == 1) + { + for (mm = 0; mm < n_inner_tot; mm++) + { + if (mode == NEBULARMODE_MATRIX_BB) + { + inner_rates[mm] = calc_pi_rate (mm, xplasma, 2, 2); + } + else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL) + { + inner_rates[mm] = calc_pi_rate (mm, xplasma, 1, 2); + } + else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) //We are using estimators - so we will need to reorder the rates from freq order to cross section order to match the electron yields. This takes time, so we only do it if we need to. + { + for (nn = 0; nn < n_inner_tot; nn++) + { + if (inner_cross[mm].nion == inner_cross_ptr[nn]->nion && inner_cross[mm].freq[0] == inner_cross_ptr[nn]->freq[0]) //Check for a match + { + inner_rates[mm] = xplasma->inner_ioniz[nn] / xplasma->density[inner_cross_ptr[nn]->nion]; + } + } + } + } + } + + /* This next line sets the partition function for each ion. This has always been the place here python calculates the + partition functions and sets the level densities for each ion. It needs to be done, or other parts of the code which rely + on sensible level populations don't work properly. In the case of the dilute blackbody, the code works well, however we do + not currently (v78) have a procedure to calucate the levels for a spectral model case. We therefore call partition + functions with mode 4 - this is a special mode which forces the ions to be in the ground state. This is reasonable for a + radiation dominated plasma, since any excitations into higher states will quickly be followed by radative de-excitation. + We should really do better though... */ + + if (mode == NEBULARMODE_MATRIX_BB) + { + partition_functions (xplasma, NEBULARMODE_ML93); // We use t_r and the radiative weight + } + else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL || mode == NEBULARMODE_MATRIX_ESTIMATORS) + { + partition_functions (xplasma, NEBULARMODE_LTE_GROUND); // Set to ground state + } + + /* Next we need to obtain an initial guess for the electron density. In the past this has been done by calculating the + hydrogen density directly from the Saha equation at the current electron temperature. In initial testing of this mode - + this seemed to be a little unstable. At the moment, we use the last ionization state to compute n_e. For the first + iteraion, this will just be calculated in LTE from the initial temperature structure of the wind, which will give almost + the same result as the original procedure, or for successive calculations, it should be a better guess. I've leftin the + original code, commented out... */ + + xne = xxne = xxxne = get_ne (xplasma->density); //Even though the abundances are fractional, we need the real electron density + + + /* xne is the current working number xxne */ + + + /* We are now going to iterate on the electron density - MAXITERATIONS is set in python.h and is currently (78) set to 200. + We would normally expect to converge much fater than this */ + + niterate = 0; + while (niterate < MAXITERATIONS) + { + + /* In order to compute charge exchange reactions, we need the number density of neutral and ionized hydrogen + the current ethod of matrix solving does not allow is to link ions multiplicatively so we need to compute + these outside the loop. If hydrogen were not dominant - this could cause iddues since the ionization(recombination) + of a metal or helium via this process would cause an identical recombination(ionization) of hydrogen. The loop below + is a bit belt and braces, since one would almost always expect hydrogen to be ion=0 and 1 */ + + for (nn = 0; nn < nions; nn++) + { + if (ion[nn].z == 1 && ion[nn].istate == 1) + { + nh1 = newden[nn] * elem_dens[ion[nn].z]; + } + if (ion[nn].z == 1 && ion[nn].istate == 2) + { + nh2 = newden[nn] * elem_dens[ion[nn].z]; + } + } + + + populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp, xne, nh1, nh2); + + + /* The array is now fully populated, and we can begin the process of solving it */ + + nrows = nions; /* This is a placeholder, we may end up removing rows and columns that have no rates (or very + low rates) connecting them to other rows. This may improve stability but will need to be + done carefully */ + + /* Here we solve the matrix equation M x = b, where x is our vector containing level populations as a fraction w.r.t the + whole element. The actual LU decomposition- the process of obtaining a solution - is done by the routine solve_matrix() */ + + a_data = (double *) calloc (sizeof (double), nrows * nrows); + populations = (double *) calloc (nrows, sizeof (double)); + + /* This b_data column matrix is the total number density for each element, placed into the row which relates to the neutral + ion. This matches the row in the rate matrix which is just 1 1 1 1 for all stages. NB, we could have chosen any line for + this. */ + + b_data = (double *) calloc (nrows, sizeof (double)); + + /* We now copy our rate matrix into the prepared matrix */ + for (mm = 0; mm < nrows; mm++) + { + for (nn = 0; nn < nrows; nn++) + { + a_data[mm * nrows + nn] = rate_matrix[mm][nn]; /* row-major */ + } + } + + for (nn = 0; nn < nrows; nn++) + { + b_data[nn] = b_temp[nn]; + } + + matrix_err = solve_matrix (a_data, b_data, nrows, populations, xplasma->nplasma); + + if (matrix_err) + { + Error ("matrix_ion_populations: %s\n", get_matrix_error_string (matrix_err)); + } + + /* free memory */ + free (a_data); + free (b_data); + + if (matrix_err == 4) + { + free (populations); + return (-1); + } + + /* Calculate level populations for macro-atoms */ + if (geo.macro_ioniz_mode == MACRO_IONIZ_MODE_ESTIMATORS) + { + int mp_err = macro_pops (xplasma, xne); + if (mp_err != EXIT_SUCCESS) + { + return -1; + } + } + + /* We now have the populations of all the ions stored in the matrix populations. We copy this data into the newden array + which will temperarily store all the populations. We wont copy this to the plasma structure until we are sure thatwe + have made things better. We just loop over all ions, and copy. The complexity in the loop is to future proof us against + the possibility that there are some ions that are not included in the matrix scheme becuse there is no route into or + out of it. */ + + for (nn = 0; nn < nions; nn++) + { + newden[nn] = 0.0; + + /* if the ion is being treated by macro_pops then use the populations just computed */ + if ((ion[nn].macro_info == TRUE) && (geo.macro_simple == FALSE) && (geo.macro_ioniz_mode == MACRO_IONIZ_MODE_ESTIMATORS) + && (modes.no_macro_pops_for_ions == FALSE)) + { + newden[nn] = xplasma->density[nn] / elem_dens[ion[nn].z]; + } + + /* if the ion is "simple" then find it's calculated ionization state in populations array */ + else + { + + for (mm = 0; mm < nrows; mm++) // inner loop over the elements of the population array + { + if (xion[mm] == nn) // if this element contains the population of the ion is question + { + newden[nn] = populations[mm]; // get the population + } + } + } + + if (newden[nn] < DENSITY_MIN) // this wil also capture the case where population doesnt have a value for this ion + newden[nn] = DENSITY_MIN; + } + free (populations); + + +/* We need to get the 'true' new electron density so we need to do a little loop here to compute it */ + + + xnew = 0.0; + for (nn = 0; nn < nions; nn++) + { + xnew += newden[nn] * (ion[nn].istate - 1) * elem_dens[ion[nn].z]; + } + + + + if (xnew < DENSITY_MIN) + xnew = DENSITY_MIN; /* fudge to keep a floor on ne */ + + + if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR || xnew < 1.e-6) /* We have converged, or have the situation where we + have a neutral plasma */ + { + break; /* Break out of the while loop - we have finished our iterations */ + } + xne = xxxne = (xnew + xne) / 2.; /* New value of ne */ + + niterate++; + + + if (niterate == MAXITERATIONS) + { + Error ("matrix_ion_populations: failed to converge for cell %i t %e nh %e xnew %e\n", xplasma->nplasma, t_e, nh, xnew); + + + for (nn = 0; nn < geo.nxfreq; nn++) + { + Log + ("numin= %e (%e) numax= %e (%e) Model= %2d PL_log_w= %e PL_alpha= %e Exp_w= %e EXP_temp= %e\n", + xplasma->fmin_mod[nn], geo.xfreq[nn], xplasma->fmax_mod[nn], + geo.xfreq[nn + 1], xplasma->spec_mod_type[nn], + xplasma->pl_log_w[nn], xplasma->pl_alpha[nn], xplasma->exp_w[nn], xplasma->exp_temp[nn]); + } + + Error ("matrix_ion_populations: xxne %e theta %e\n", xxne); + + return (-1); /* If we get to MAXITERATIONS, we return without copying the new populations into plasma */ + } + } /* This is the end of the iteration loop */ + + + xplasma->ne = xnew; + for (nn = 0; nn < nions; nn++) + { + /* If statement added here to suppress interference with macro populations */ + if (ion[nn].macro_info == FALSE || geo.macro_ioniz_mode == MACRO_IONIZ_MODE_NO_ESTIMATORS || geo.macro_simple == TRUE) + { + xplasma->density[nn] = newden[nn] * elem_dens[ion[nn].z]; //We return to absolute densities here + } + if ((sane_check (xplasma->density[nn])) || (xplasma->density[nn] < 0.0)) + Error ("matrix_ion_populations: ion %i has population %8.4e in cell %i\n", nn, xplasma->density[nn], xplasma->nplasma); + } + + xplasma->ne = get_ne (xplasma->density); + + if (n_charge_exchange > 0) + { + xplasma->heat_ch_ex = ch_ex_heat (&wmain[xplasma->nwind], xplasma->t_e); //Compute the charge exchange heating + xplasma->heat_tot += xplasma->heat_ch_ex; + } + + /*We now need to populate level densities in order to later calculate line emission (for example). + We call partition functions to do this. At present, we do not have a method for accurately computing + level populations for modelled mean intensities. We get around this by setting all level populations + other than ground state to zero - this is a pretty good situation for very dilute radiation fields, + which is what we are normally looking at, however one should really do better in the long term. + */ + + partition_functions (xplasma, NEBULARMODE_LTE_GROUND); + + return (0); +} + diff --git a/source/python.h b/source/python.h index 8c54b88b4..a52204bcd 100644 --- a/source/python.h +++ b/source/python.h @@ -1240,6 +1240,11 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; #define NEBULARMODE_MATRIX_SPECTRALMODEL 9 /**< matrix solver spectral model */ #define NEBULARMODE_MATRIX_ESTIMATORS 10 /**< matrix solver spectral model */ +#define NEBULAR_MATRIX_MULTISHOT 11 /**< matrix solver spectral model based on power laws which + * updates T_e multiple times before arriving at a final + * solution + */ + // modes for the wind_luminosity routine #define MODE_OBSERVER_FRAME_TIME 0 #define MODE_CMF_TIME 1 diff --git a/source/saha.c b/source/saha.c index 2da06f038..5c1d3c145 100644 --- a/source/saha.c +++ b/source/saha.c @@ -80,6 +80,10 @@ nebular_concentrations (xplasma, mode) m = matrix_ion_populations (xplasma, mode); } + else if (mode == NEBULAR_MATRIX_MULTISHOT) + { + m = matrix_ion_populations2 (xplasma, mode); + } else { Error ("nebular_concentrations: Unknown mode %d\n", mode); diff --git a/source/templates.h b/source/templates.h index 14892e70c..4b7f4a7f0 100644 --- a/source/templates.h +++ b/source/templates.h @@ -200,6 +200,10 @@ double alpha_st_e_integrand(double freq, void *params); int update_banded_estimators(PlasmaPtr xplasma, PhotPtr p, double ds, double w_ave, int ndom); int update_flux_estimators(PlasmaPtr xplasma, PhotPtr phot_mid, double ds_obs, double w_ave, int ndom); int update_force_estimators(PlasmaPtr xplasma, PhotPtr p, PhotPtr phot_mid, double ds, double w_ave, int ndom, double z, double frac_ff, double frac_auger, double frac_tot); +double planck_spectral_radiance(double nu, void *params); +double nu_times_radiance(double nu, void *params); +double mean_frequency_nu_range(double T, double nu_min, double nu_max); +double estimate_temperature_from_mean_frequency(double mean_nu_target, double nu_min, double nu_max, double initial_guess); int normalise_simple_estimators(PlasmaPtr xplasma); void update_persistent_directional_flux_estimators(int nplasma, double flux_persist_scale); /* extract.c */ @@ -355,6 +359,8 @@ int invert_matrix(double *matrix, double *inverted_matrix, int num_rows); /* matrix_ion.c */ int matrix_ion_populations(PlasmaPtr xplasma, int mode); int populate_ion_rate_matrix(double rate_matrix[nions][nions], double pi_rates[nions], double inner_rates[n_inner_tot], double rr_rates[nions], double b_temp[nions], double xne, double nh1, double nh2); +/* matrix_ion2.c */ +int matrix_ion_populations2(PlasmaPtr xplasma, int mode); /* models_extern_init.c */ /* para_update.c */ int get_parallel_nrange(int rank, int ntotal, int nproc, int *my_nmin, int *my_nmax); From 81475186558cc2a1e09cf3d098ff61f678f5d968 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 13 Sep 2024 15:56:58 -0500 Subject: [PATCH 10/23] Additinal diagnostics The normal modes are unchanged but various changes have been made in matrix_ion2 to slow down how rapidly the temperature can change. This still is not fundamentally different from the previous approach. --- source/ionization.c | 13 ++++++++++++ source/matrix_ion2.c | 49 ++++++++++++++++++++++++++------------------ source/photon_gen.c | 4 ++++ source/python.h | 2 +- source/saha.c | 2 +- 5 files changed, 48 insertions(+), 22 deletions(-) diff --git a/source/ionization.c b/source/ionization.c index 11d41c838..d8efc0e3c 100644 --- a/source/ionization.c +++ b/source/ionization.c @@ -125,6 +125,19 @@ ion_abundances (PlasmaPtr xplasma, int mode) convergence (xplasma); } + else if (mode == IONMODE_MATRIX_MULTISHOT) + { +/* New mode to carry more extensive search for best fit + Note that it unclear that we could not simplify a lot of this + routine but what is done here is most analagous to the other modes + excetp we avoid one_shot*/ + spectral_estimators (xplasma); + update_old_plasma_variables (xplasma); + // ireturn = one_shot (xplasma, mode); + ireturn = nebular_concentrations (xplasma, mode); + + convergence (xplasma); + } else { Error ("ion_abundances: Could not calculate abundances for mode %d\n", mode); diff --git a/source/matrix_ion2.c b/source/matrix_ion2.c index 2caf0c09f..2dd30ca1e 100644 --- a/source/matrix_ion2.c +++ b/source/matrix_ion2.c @@ -69,7 +69,6 @@ matrix_ion_populations2 (xplasma, mode) int matrix_err, niterate; //counters for errors and the number of iterations we have tried to get a converged electron density double xnew; int xion[nions]; // This array keeps track of what ion is in each line -//OLD int xelem[nions]; // This array keeps track of the element for each ion double pi_rates[nions]; //photoionization rate coefficients double rr_rates[nions]; //radiative recombination rate coefficients double inner_rates[n_inner_tot]; //This array contains the rates for each of the inner shells. Where they go to requires the electron yield array @@ -79,8 +78,13 @@ matrix_ion_populations2 (xplasma, mode) /* Copy some quantities from the cell into local variables */ nh = xplasma->rho * rho2nh; // The number density of hydrogen ions - computed from density - t_e = xplasma->t_e; // The electron temperature in the cell - used for collisional processes + // t_e = xplasma->t_e; // The electron temperature in the cell - used for collisional processes + xplasma->t_e_old = xplasma->t_e; + // t_e = calc_te (xplasma, 0.7 * xplasma->t_e, 1.3 * xplasma->t_e); + t_e = calc_te (xplasma, 0.95 * xplasma->t_e, 1.05 * xplasma->t_e); + xplasma->t_e = t_e; + zero_emit (xplasma->t_e); //Get the heating and cooling rates correctly for the new temperature /* We now calculate the total abundances for each element to allow us to use fractional abundances */ @@ -129,27 +133,21 @@ matrix_ion_populations2 (xplasma, mode) { pi_rates[mm] = calc_pi_rate (mm, xplasma, 1, 1); // PI rate for an explicit spectral model } + else if (mode == NEBULARMODE_MATRIX_MULTISHOT) + { + pi_rates[mm] = calc_pi_rate (mm, xplasma, 1, 1); // PI rate for an explicit spectral model + } else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) { pi_rates[mm] = xplasma->ioniz[mm] / xplasma->density[mm]; // PI rate logged during the photon passage } else { - Error ("matrix_ion_populations: Unknown mode %d\n", mode); + Error ("matrix_ion_populations2: Unknown mode %d\n", mode); Exit (0); } } -/* ksl: This whole lope is not actually used, and produces warnings in gcc11 */ -//OLD for (nn = 0; nn < nelements; nn++) -//OLD { -//OLD if (ion[mm].z == ele[nn].z) -//OLD { -//OLD xelem[mm] = nn; /* xelem logs which element each row in the arrays refers to. This is important because we need -//OLD to know what the total density will be for a group of rows all representing the same -//OLD element. */ -//OLD } -//OLD } } @@ -169,6 +167,10 @@ matrix_ion_populations2 (xplasma, mode) { inner_rates[mm] = calc_pi_rate (mm, xplasma, 1, 2); } + else if (mode == NEBULARMODE_MATRIX_MULTISHOT) + { + inner_rates[mm] = calc_pi_rate (mm, xplasma, 1, 2); + } else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) //We are using estimators - so we will need to reorder the rates from freq order to cross section order to match the electron yields. This takes time, so we only do it if we need to. { for (nn = 0; nn < n_inner_tot; nn++) @@ -198,6 +200,10 @@ matrix_ion_populations2 (xplasma, mode) { partition_functions (xplasma, NEBULARMODE_LTE_GROUND); // Set to ground state } + else if (mode == NEBULARMODE_MATRIX_MULTISHOT) + { + partition_functions (xplasma, NEBULARMODE_LTE_GROUND); // Set to ground state + } /* Next we need to obtain an initial guess for the electron density. In the past this has been done by calculating the hydrogen density directly from the Saha equation at the current electron temperature. In initial testing of this mode - @@ -277,7 +283,7 @@ matrix_ion_populations2 (xplasma, mode) if (matrix_err) { - Error ("matrix_ion_populations: %s\n", get_matrix_error_string (matrix_err)); + Error ("matrix_ion_populations2: %s\n", get_matrix_error_string (matrix_err)); } /* free memory */ @@ -351,11 +357,15 @@ matrix_ion_populations2 (xplasma, mode) xnew = DENSITY_MIN; /* fudge to keep a floor on ne */ - if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR || xnew < 1.e-6) /* We have converged, or have the situation where we - have a neutral plasma */ +// if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR || xnew < 1.e-6) /* We have converged, or have the situation where we +// have a neutral plasma */ + if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR) /* We have converged, or have the situation where we + have a neutral plasma */ { break; /* Break out of the while loop - we have finished our iterations */ } + + Log ("xxxxxx %3d xnew %.2e xne %.2e -> %.2e \n", niterate, xnew, xne, (xnew + xne) / 2.); xne = xxxne = (xnew + xne) / 2.; /* New value of ne */ niterate++; @@ -363,7 +373,7 @@ matrix_ion_populations2 (xplasma, mode) if (niterate == MAXITERATIONS) { - Error ("matrix_ion_populations: failed to converge for cell %i t %e nh %e xnew %e\n", xplasma->nplasma, t_e, nh, xnew); + Error ("matrix_ion_populations2: failed to converge for cell %i t %e nh %e xnew %e\n", xplasma->nplasma, t_e, nh, xnew); for (nn = 0; nn < geo.nxfreq; nn++) @@ -375,7 +385,7 @@ matrix_ion_populations2 (xplasma, mode) xplasma->pl_log_w[nn], xplasma->pl_alpha[nn], xplasma->exp_w[nn], xplasma->exp_temp[nn]); } - Error ("matrix_ion_populations: xxne %e theta %e\n", xxne); + Error ("matrix_ion_populations2: xxne %e theta %e\n", xxne); return (-1); /* If we get to MAXITERATIONS, we return without copying the new populations into plasma */ } @@ -391,7 +401,7 @@ matrix_ion_populations2 (xplasma, mode) xplasma->density[nn] = newden[nn] * elem_dens[ion[nn].z]; //We return to absolute densities here } if ((sane_check (xplasma->density[nn])) || (xplasma->density[nn] < 0.0)) - Error ("matrix_ion_populations: ion %i has population %8.4e in cell %i\n", nn, xplasma->density[nn], xplasma->nplasma); + Error ("matrix_ion_populations2: ion %i has population %8.4e in cell %i\n", nn, xplasma->density[nn], xplasma->nplasma); } xplasma->ne = get_ne (xplasma->density); @@ -413,4 +423,3 @@ matrix_ion_populations2 (xplasma, mode) return (0); } - diff --git a/source/photon_gen.c b/source/photon_gen.c index ae6e8cb12..2ebb4f228 100644 --- a/source/photon_gen.c +++ b/source/photon_gen.c @@ -159,6 +159,10 @@ define_phot (p, f1, f2, nphot_tot, ioniz_or_extract, iwind, freq_sampling) iphot_start += xband.nphot[n]; } + else + { + Error ("photon_gen: No photons for band %d\n", n); + } } diff --git a/source/python.h b/source/python.h index a52204bcd..015e9820c 100644 --- a/source/python.h +++ b/source/python.h @@ -1240,7 +1240,7 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; #define NEBULARMODE_MATRIX_SPECTRALMODEL 9 /**< matrix solver spectral model */ #define NEBULARMODE_MATRIX_ESTIMATORS 10 /**< matrix solver spectral model */ -#define NEBULAR_MATRIX_MULTISHOT 11 /**< matrix solver spectral model based on power laws which +#define NEBULARMODE_MATRIX_MULTISHOT 11 /**< matrix solver spectral model based on power laws which * updates T_e multiple times before arriving at a final * solution */ diff --git a/source/saha.c b/source/saha.c index 5c1d3c145..21f5b77c8 100644 --- a/source/saha.c +++ b/source/saha.c @@ -80,7 +80,7 @@ nebular_concentrations (xplasma, mode) m = matrix_ion_populations (xplasma, mode); } - else if (mode == NEBULAR_MATRIX_MULTISHOT) + else if (mode == NEBULARMODE_MATRIX_MULTISHOT) { m = matrix_ion_populations2 (xplasma, mode); } From ebd80ed762876de2aa8fff921fdef166e4d1a949 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 13 Sep 2024 16:54:47 -0500 Subject: [PATCH 11/23] Add nxtot to windsave2fits --- source/windsave2fits.c | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/source/windsave2fits.c b/source/windsave2fits.c index 9415de647..f56e2e186 100644 --- a/source/windsave2fits.c +++ b/source/windsave2fits.c @@ -244,13 +244,13 @@ write_spectra_model_table (fitsfile *fptr) int i, j, k; int num_rows, num_cols; int nbands; - int *ichoice, *iplasma, *iband; + int *ichoice, *iplasma, *iband, *nxtot; float *exp_w, *exp_temp, *pl_log_w, *pl_alpha; nbands = geo.nxfreq; num_rows = nbands * NPLASMA; printf ("xtest %d %d\n", nbands, NPLASMA); - num_cols = 7; // for now + num_cols = 8; // for now char *table_name = "spec_model"; ichoice = calloc (num_rows, sizeof (int *)); @@ -260,6 +260,7 @@ write_spectra_model_table (fitsfile *fptr) exp_temp = calloc (num_rows, sizeof (float *)); pl_log_w = calloc (num_rows, sizeof (float *)); pl_alpha = calloc (num_rows, sizeof (float *)); + nxtot = calloc (num_rows, sizeof (int *)); k = 0; for (i = 0; i < NPLASMA; i++) @@ -273,6 +274,8 @@ write_spectra_model_table (fitsfile *fptr) exp_temp[k] = plasmamain[i].exp_temp[j]; pl_log_w[k] = plasmamain[i].pl_log_w[j]; pl_alpha[k] = plasmamain[i].pl_alpha[j]; + nxtot[k] = plasmamain[i].nxtot[j]; + printf ("nxtot %d\n", plasmamain[i].nxtot[j]); k++; } } @@ -280,9 +283,9 @@ write_spectra_model_table (fitsfile *fptr) int status = 0; // Define the names, formats, and units for each column - char *ttype[] = { "nplasma", "band", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha" }; // Column names - char *tform[] = { "J", "J", "J", "E", "E", "E", "E" }; // Formats: 'J' for integer, 'E' for float - char *tunit[] = { "", "", "", "", "", "", "" }; // Units + char *ttype[] = { "nplasma", "band", "spec_mod_type", "exp_w", "exp_temp", "pl_log_w", "pl_alpha", "nxtot" }; // Column names + char *tform[] = { "J", "J", "J", "E", "E", "E", "E", "J" }; // Formats: 'J' for integer, 'E' for float + char *tunit[] = { "", "", "", "", "", "", "", "" }; // Units // Create a new binary table extension if (fits_create_tbl (fptr, BINARY_TBL, num_rows, num_cols, ttype, tform, tunit, NULL, &status)) @@ -359,6 +362,14 @@ write_spectra_model_table (fitsfile *fptr) return status; } + // Write the integer data to the eightth column + if (fits_write_col (fptr, TINT, 8, 1, 1, num_rows, nxtot, &status)) + { + fits_report_error (stderr, status); + return status; + } + + return 0; } From a0ccda859cb15d84d2313eaf4db5c42f051a2756 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 20 Sep 2024 10:30:26 -0500 Subject: [PATCH 12/23] Small updates --- source/matrix_ion2.c | 16 ++++++++++++++-- source/python.h | 2 +- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/source/matrix_ion2.c b/source/matrix_ion2.c index 2dd30ca1e..feae57fad 100644 --- a/source/matrix_ion2.c +++ b/source/matrix_ion2.c @@ -81,8 +81,14 @@ matrix_ion_populations2 (xplasma, mode) // t_e = xplasma->t_e; // The electron temperature in the cell - used for collisional processes xplasma->t_e_old = xplasma->t_e; + + // t_e = calc_te (xplasma, 0.7 * xplasma->t_e, 1.3 * xplasma->t_e); + Log ("xxxxxx before calc_te %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, + xplasma->cool_tot); t_e = calc_te (xplasma, 0.95 * xplasma->t_e, 1.05 * xplasma->t_e); + Log ("xxxxxx after calc_te %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, + xplasma->cool_tot); xplasma->t_e = t_e; zero_emit (xplasma->t_e); //Get the heating and cooling rates correctly for the new temperature @@ -221,6 +227,8 @@ matrix_ion_populations2 (xplasma, mode) /* We are now going to iterate on the electron density - MAXITERATIONS is set in python.h and is currently (78) set to 200. We would normally expect to converge much fater than this */ + Log ("xxxxxx before_loop %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, + xplasma->cool_tot); niterate = 0; while (niterate < MAXITERATIONS) { @@ -359,8 +367,10 @@ matrix_ion_populations2 (xplasma, mode) // if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR || xnew < 1.e-6) /* We have converged, or have the situation where we // have a neutral plasma */ - if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR) /* We have converged, or have the situation where we - have a neutral plasma */ +// if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR) /* We have converged, or have the situation where we +// have a neutral plasma */ + if (fabs ((xne - xnew) / (xnew)) < 0.01) /* We have converged, or have the situation where we + have a neutral plasma */ { break; /* Break out of the while loop - we have finished our iterations */ } @@ -420,6 +430,8 @@ matrix_ion_populations2 (xplasma, mode) */ partition_functions (xplasma, NEBULARMODE_LTE_GROUND); + Log ("xxxxxx after At_end %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, + xplasma->cool_tot); return (0); } diff --git a/source/python.h b/source/python.h index 015e9820c..bb1959dc7 100644 --- a/source/python.h +++ b/source/python.h @@ -1208,7 +1208,7 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; //These constants are used in the various routines which compute ionization state #define SAHA 4.82907e15 /**< 2* (2.*PI*MELEC*k)**1.5 / h**3 (Calculated in constants) */ #define MAXITERATIONS 200 /**< /The number of loops to do to try to converge in ne */ -#define FRACTIONAL_ERROR 0.03 /**< The change in n_e which causes a break out of the loop for ne */ +#define FRACTIONAL_ERROR 0.001 /**< The change in n_e which causes a break out of the loop for ne */ #define THETAMAX 1e4 /**< Used in initial calculation of n_e */ #define MIN_TEMP 100. /**< ??? this is another minimum temperature, which is * used in saha.c and variable_temperature.c ) From 8164f01f8b5c035652453f29bf2e92aebb246092 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 20 Sep 2024 15:42:10 -0500 Subject: [PATCH 13/23] Prevent duplicate column name in heat.txt --- source/windsave2table_sub.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/windsave2table_sub.c b/source/windsave2table_sub.c index ba3e72f94..9a5c3b9f7 100644 --- a/source/windsave2table_sub.c +++ b/source/windsave2table_sub.c @@ -498,10 +498,10 @@ create_heat_table (ndom, rootname) strcpy (column_name[25], "heat_shock"); c[26] = get_one (ndom, "heat_lines_macro"); - strcpy (column_name[26], "heat_lines_macro"); + strcpy (column_name[26], "ht_ln_macro"); c[27] = get_one (ndom, "heat_photo_macro"); - strcpy (column_name[27], "heat_photo_macro"); + strcpy (column_name[27], "ht_ph_macro"); /* This should be the maximum number above +1 */ ncols = 28; From 4468c85539d6ad871e4abba9e63716f25931e881 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Tue, 8 Oct 2024 15:44:59 -0400 Subject: [PATCH 14/23] Print disk_diag each cycle if also saving tables Previously setting make_tables to yes, just ran windsave totable on each cycle, now the disk_diag file is also saved. Cleaned up writing of the disk_diag file when no photons hit disk. --- source/disk_init.c | 20 ++++++++++++-------- source/run.c | 7 +++++++ 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/source/disk_init.c b/source/disk_init.c index b87121e39..d79093f59 100644 --- a/source/disk_init.c +++ b/source/disk_init.c @@ -161,8 +161,8 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot) } else if (spectype == SPECTYPE_BB_FCOL) { - fcol = disk_colour_correction(t); - emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0); + fcol = disk_colour_correction (t); + emit = emittance_bb (freqmin, freqmax, fcol * t) / pow (fcol, 4.0); } else { @@ -220,8 +220,8 @@ disk_init (rmin, rmax, m, mdot, freqmin, freqmax, ioniz_or_extract, ftot) } else if (spectype == SPECTYPE_BB_FCOL) { - fcol = disk_colour_correction(t); - emit = emittance_bb (freqmin, freqmax, fcol * t) / pow(fcol, 4.0); + fcol = disk_colour_correction (t); + emit = emittance_bb (freqmin, freqmax, fcol * t) / pow (fcol, 4.0); } else { @@ -386,9 +386,7 @@ qdisk_save (diskfile, ztot) double area, theat, ttot; qptr = fopen (diskfile, "w"); fprintf (qptr, - " r v zdisk t_disk g heat nhit nhit/emit t_heat t_irrad W_irrad t_tot\n"); -// fprintf (qptr, -// "---------- ---------- ---------- --------- --------- --------- ----- --------- --------- --------- --------- ---------\n"); + " r v zdisk t_disk g heat nhit nhit/emit t_heat t_irrad W_irrad t_tot\n"); for (n = 0; n < NRINGS; n++) { @@ -411,10 +409,16 @@ qdisk_save (diskfile, ztot) //Basic conversion from freq to T qdisk.w[n] = qdisk.heat[n] / (4. * PI * STEFAN_BOLTZMANN * area * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n]); } + else + { + qdisk.ave_freq[n] = 0.0; + qdisk.t_hit[n] = 0.0; + qdisk.w[n] = 0.0; + } ttot = pow (qdisk.t[n], 4) + pow (theat, 4); ttot = pow (ttot, 0.25); fprintf (qptr, - "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %5d %8.3e %8.3e %8.3e %8.3e %8.3e\n", + "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %5d %8.3e %8.3e %8.3e %10.3e %8.3e\n", qdisk.r[n], qdisk.v[n], zdisk (qdisk.r[n]), qdisk.t[n], qdisk.g[n], qdisk.heat[n], qdisk.nhit[n], qdisk.heat[n] * NRINGS / ztot, theat, qdisk.t_hit[n], qdisk.w[n], ttot); } diff --git a/source/run.c b/source/run.c index cd4956072..2a9c7791e 100644 --- a/source/run.c +++ b/source/run.c @@ -324,6 +324,13 @@ calculate_ionization (restart_stat) if (geo.disk_type != DISK_NONE) { qdisk_save (files.disk, ztot); + if (modes.make_tables) + { + strcpy (dummy, ""); + sprintf (dummy, "diag_%.100s/%.100s.%02d.disk.diag", files.root, files.root, geo.wcycle + 1); + qdisk_save (dummy, ztot); + + } } #ifdef MPI_ON } From 69d70f9310f19985505d3b90392793bab4c8791a Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Wed, 9 Oct 2024 10:52:23 -0400 Subject: [PATCH 15/23] Version updatre and documentation --- source/Makefile | 2 +- source/ionization.c | 10 +++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/source/Makefile b/source/Makefile index 7fa3ce367..c0d0b1e9d 100644 --- a/source/Makefile +++ b/source/Makefile @@ -147,7 +147,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas $(CUDA_LIBS) #Note that version should be a single string without spaces. -VERSION = 88e +VERSION = 88f diff --git a/source/ionization.c b/source/ionization.c index d8efc0e3c..4d1697899 100644 --- a/source/ionization.c +++ b/source/ionization.c @@ -585,7 +585,8 @@ calc_te (PlasmaPtr xplasma, double tmin, double tmax) * at a specifc temperature. * * @details - * This routine is used to estiamate a new temperature for a cell given the + * This routine is used in the process of estimating a new temperature for a cell + * given the * heating of the cell in the previous ionization cycle. When 0 the heating * and cooling are matched. * @@ -596,6 +597,13 @@ calc_te (PlasmaPtr xplasma, double tmin, double tmax) * * * ### Notes ### + * + * This routine is poorly named; what it does is simply calculate the + * sum the heating that occurred during the current cycle, and calculate + * the cooling for the input temperature given the current abundances + * it returns the difference between these two numbers. It does update + * t_e in the plasma cell. + * * The abundances of ions in the cell are not modified. Results are stored * in the cell of interest. This routine is used in connection with a zero * finding routine From 5399916bc4834ae8994cf4226be1b79ef551a88b Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Fri, 11 Oct 2024 15:21:39 -0500 Subject: [PATCH 16/23] Improvments to the disk.diag file Variables have been added to the disk structure to allow for better tracking of where photons are emitted from the disk, and where they hit the disk. --- docs/sphinx/source/output/diagnostics.rst | 5 + .../py_progs/MakeMacro/MacroCombine.rst | 3 + source/communicate_plasma.c | 12 +- source/disk_init.c | 117 ++++++++++++++++-- source/disk_photon_gen.c | 4 +- source/python.h | 9 +- source/run.c | 16 +-- source/templates.h | 3 +- 8 files changed, 134 insertions(+), 35 deletions(-) diff --git a/docs/sphinx/source/output/diagnostics.rst b/docs/sphinx/source/output/diagnostics.rst index 48874b679..485f0debb 100644 --- a/docs/sphinx/source/output/diagnostics.rst +++ b/docs/sphinx/source/output/diagnostics.rst @@ -10,3 +10,8 @@ where the in a multiprocessor run, the number refers to the log of a specific th Inspecting these logs is important for understanding whether a Python run is successful, and for understanding why if failed if that is the case. + +Additionally in this directory one will find a file **whatever.disk.diag** if the system contained a disk. This file contains information +about the structure of the disk as a function of disk radius. It also indicates how many photons hit the disk with how much total energy +as a function of radius. It provides several estimates of the effective temperature of the radiation that hits the disk, and determines +how the disk temperature of the disk would so that it would reradiate the sum of the viscous energy and the energy associated with illumination. diff --git a/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst b/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst index 5b3744935..34c24f0be 100644 --- a/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst +++ b/docs/sphinx/source/py_progs/MakeMacro/MacroCombine.rst @@ -14,8 +14,11 @@ .. autosummary:: add_gtot + combine_collisions plot_xsec + read_collisions read_phot + redo_collisions redo_lines redo_phot reweight diff --git a/source/communicate_plasma.c b/source/communicate_plasma.c index 768e46a56..1e84baa98 100644 --- a/source/communicate_plasma.c +++ b/source/communicate_plasma.c @@ -986,11 +986,11 @@ reduce_simple_estimators (void) ion_helper2 = calloc (sizeof (double), NPLASMA * nions); inner_ion_helper = calloc (sizeof (double), NPLASMA * n_inner_tot); inner_ion_helper2 = calloc (sizeof (double), NPLASMA * n_inner_tot); - /* JM -- added routine to average the qdisk quantities. The 2 is because - we only have two doubles to worry about (heat and ave_freq) and + /* Routine to average the qdisk quantities. The 3 is because + we have three doubles to worry about (emit, heat and ave_freq) and two integers (nhit and nphot) */ - qdisk_helper = calloc (sizeof (double), NRINGS * 2); - qdisk_helper2 = calloc (sizeof (double), NRINGS * 2); + qdisk_helper = calloc (sizeof (double), NRINGS * 3); + qdisk_helper2 = calloc (sizeof (double), NRINGS * 3); flux_helper = calloc (sizeof (double), NPLASMA * NFLUX_ANGLES * 3); flux_helper2 = calloc (sizeof (double), NPLASMA * NFLUX_ANGLES * 3); @@ -1068,6 +1068,7 @@ reduce_simple_estimators (void) { qdisk_helper[mpi_i] = qdisk.heat[mpi_i] / np_mpi_global; qdisk_helper[mpi_i + NRINGS] = qdisk.ave_freq[mpi_i] / np_mpi_global; + qdisk_helper[mpi_i + 2 * NRINGS] = qdisk.emit[mpi_i] / np_mpi_global; } @@ -1080,7 +1081,7 @@ reduce_simple_estimators (void) MPI_Allreduce (flux_helper, flux_helper2, NPLASMA * 3 * NFLUX_ANGLES, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce (ion_helper, ion_helper2, NPLASMA * nions, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce (inner_ion_helper, inner_ion_helper2, NPLASMA * n_inner_tot, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce (qdisk_helper, qdisk_helper2, 2 * NRINGS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce (qdisk_helper, qdisk_helper2, 3 * NRINGS, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); /* Unpacking stuff */ for (mpi_i = 0; mpi_i < NPLASMA; mpi_i++) @@ -1157,6 +1158,7 @@ reduce_simple_estimators (void) { qdisk.heat[mpi_i] = qdisk_helper2[mpi_i]; qdisk.ave_freq[mpi_i] = qdisk_helper2[mpi_i + NRINGS]; + qdisk.emit[mpi_i] = qdisk_helper2[mpi_i + 2 * NRINGS]; } /* now we've done all the doubles so we can free their helper arrays */ diff --git a/source/disk_init.c b/source/disk_init.c index d79093f59..2411e0571 100644 --- a/source/disk_init.c +++ b/source/disk_init.c @@ -354,13 +354,72 @@ qdisk_init (rmin, rmax, m, mdot) } +/**********************************************************/ +/** + * @brief Reinitialize portions a structure (qdisk) for recording information about + * photons/energy impinging + * on the disk. + * + * @return Always return zero + * + * + * ###Notes### + * + * The structure qdisk records information about how many photons hit the disk + * and what the impact of these photons would be if these + * photons are used to modify the temperatures structure + * of the disk + * + * This routine records information about the + * photons which were emitted from the disk + * and zeros various arrays that will record information + * about photons that will hit the disk + * + * + **********************************************************/ + +int +qdisk_reinit (p) + PhotPtr p; +{ + int nphot, i, n; + double rho; + struct photon pp; + for (n = 0; n < NRINGS; n++) + { + qdisk.emit[n] = qdisk.nhit[n] = qdisk.heat[n] = qdisk.nphot[n] = qdisk.w[n] = qdisk.ave_freq[n] = 0; + } + + for (nphot = 0; nphot < NPHOT; nphot++) + { + stuff_phot (&p[nphot], &pp); + if (pp.origin_orig == PTYPE_DISK) + { + rho = sqrt (pp.x[0] * pp.x[0] + pp.x[1] * pp.x[1]); + + i = 0; + while (rho > qdisk.r[i] && i < NRINGS - 1) + i++; + i--; /* So that the heating refers to the heating between i and i+1 */ + + qdisk.emit[i] += pp.w; + qdisk.nphot[i] += 1; + + + } + } + + return (0); +} + + /**********************************************************/ /** * @brief Save the information in the qdisk structure to a file * * @param [in] char * diskfile Name of the file whihc is writteen - * @param [in] double ztot The total luminosity of the disk as - * calculate over multiple cycles + * @param [in] int ichoice A switch that just says whether the weighted + * frequency has been converted to a real frequency (0 = No, 1=Yes) * @return Always returns 0 * * The routine reformats the data about disk heating which has @@ -368,59 +427,91 @@ qdisk_init (rmin, rmax, m, mdot) * * ###Notes### * - * The data concerning heating by the disk is built up during - * a the ionization cycles * * The file that is produced should be readable as an astropy * table * + * Here is a definition of the columnts + * + * * r the radius of the wirng + * * v the velocity of the ring + * * zdisk the zheight ot the rign + * * t_disk the temperature derived from viscous heating or the input grid + * * g the log of the gravity + * * emit the luminosity of photons emitted by the ring in this cycle + * * nemit the number photons emitted by the ring in this cycle + * * heat the total luminosity of photons that hit the ring + * * nhit the number of photon bundles hitting a ring + * * ehit/emit the ratio of the energy of photons that hit the ring/to those + * to those that are emitted by the ring + * * t_heat the temperature calculated from the energy of the + * photons hitting a ring + * * t_freq The temperature calculated from the weighted frequency of + * photons hitting the ring + * * W_freq Assuming t_freq for the temperature, this gives ther + * ratio of the energy flux to LTE + * * t_rerad The temperature required to radiate away both the + * the viscous heating and the heating from photons + * in the next cycle +* **********************************************************/ int -qdisk_save (diskfile, ztot) +qdisk_save (diskfile, ichoice) char *diskfile; - double ztot; + int ichoice; { FILE *qptr; int n; double area, theat, ttot; + double ratio[NRINGS]; qptr = fopen (diskfile, "w"); fprintf (qptr, - " r v zdisk t_disk g heat nhit nhit/emit t_heat t_irrad W_irrad t_tot\n"); + " r v zdisk t_disk g emit nemit heat nhit ehit/emit t_heat t_freq W_freq t_rerad\n"); for (n = 0; n < NRINGS; n++) { if (n < NRINGS - 1) { + // The factor of 2 comes from the fact that the disk has two sides area = (2. * PI * (qdisk.r[n + 1] * qdisk.r[n + 1] - qdisk.r[n] * qdisk.r[n])); } else { area = 0; } - theat = qdisk.heat[n] / area; - theat = pow (theat / STEFAN_BOLTZMANN, 0.25); - //theat is temperature if no internal energy production if (qdisk.nhit[n] > 0) { - qdisk.ave_freq[n] /= qdisk.heat[n]; + /* During photon transport ave_freq is actually the w x frequency */ + if (ichoice == 0) + qdisk.ave_freq[n] /= qdisk.heat[n]; + qdisk.t_hit[n] = PLANCK * qdisk.ave_freq[n] / (BOLTZMANN * 3.832); //Basic conversion from freq to T qdisk.w[n] = qdisk.heat[n] / (4. * PI * STEFAN_BOLTZMANN * area * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n] * qdisk.t_hit[n]); + ratio[n] = 99.; + if (qdisk.emit[n] > 0.0) + ratio[n] = qdisk.heat[n] / qdisk.emit[n]; + theat = qdisk.heat[n] / area; + theat = pow (theat / STEFAN_BOLTZMANN, 0.25); + //theat is temperature if no internal energy production + } else { qdisk.ave_freq[n] = 0.0; qdisk.t_hit[n] = 0.0; qdisk.w[n] = 0.0; + ratio[n] = 0.0; + theat = 0; } ttot = pow (qdisk.t[n], 4) + pow (theat, 4); ttot = pow (ttot, 0.25); fprintf (qptr, - "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %5d %8.3e %8.3e %8.3e %10.3e %8.3e\n", + "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %10lld %8.3e %10lld %8.3e %8.3e %8.3e %10.3e %8.3e\n", qdisk.r[n], qdisk.v[n], zdisk (qdisk.r[n]), qdisk.t[n], qdisk.g[n], - qdisk.heat[n], qdisk.nhit[n], qdisk.heat[n] * NRINGS / ztot, theat, qdisk.t_hit[n], qdisk.w[n], ttot); + qdisk.emit[n], qdisk.nphot[n], qdisk.heat[n], qdisk.nhit[n], ratio[n], theat, qdisk.t_hit[n], qdisk.w[n], ttot); } fclose (qptr); diff --git a/source/disk_photon_gen.c b/source/disk_photon_gen.c index 9d5f12cd1..225d13add 100644 --- a/source/disk_photon_gen.c +++ b/source/disk_photon_gen.c @@ -189,7 +189,7 @@ photo_gen_disk (p, weight, f1, f2, spectype, istart, nphot) } else if (spectype == SPECTYPE_BB_FCOL) { - fcol = disk_colour_correction(disk.t[nring]); + fcol = disk_colour_correction (disk.t[nring]); p[i].freq = planck (fcol * disk.t[nring], freqmin, freqmax); } else if (spectype == SPECTYPE_UNIFORM) @@ -277,7 +277,7 @@ disk_photon_summary (filename, mode) x = 0; else x = disk.nphot[n] / (disk.r[n] - disk.r[n - 1]); - fprintf (ptr, "%d %8.2e %8.2e %8d %8.2e %8d\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); + fprintf (ptr, "%d %8.2e %8.2e %8lld %8.2e %8lld\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); } fclose (ptr); diff --git a/source/python.h b/source/python.h index bb1959dc7..3f4598060 100644 --- a/source/python.h +++ b/source/python.h @@ -755,12 +755,13 @@ struct xdisk double t[NRINGS]; /**< The temperature at the middle of the annulus */ double g[NRINGS]; /**< The gravity at the middle of the annulus */ double v[NRINGS]; /**< The velocity at the middle of the annulus */ + double emit[NRINGS]; /**< The radiative energy of the photons emitted from the disk */ double heat[NRINGS]; /**< The total energy flux of photons hitting each annulus */ double ave_freq[NRINGS]; /**< The flux weighted average of frequency of photons hitting each annulus */ - double w[NRINGS]; /**< The radiative weight of the photons that hit the disk */ - double t_hit[NRINGS]; /**< The effective T of photons hitting the disk */ - int nphot[NRINGS]; /**< The number of photons created in each annulus */ - int nhit[NRINGS]; /**< The number of photons which hit each annulus */ + double w[NRINGS]; /**< The weight relative to a BB of the photons that hit the disk */ + double t_hit[NRINGS]; /**< The effective T of photons hitting the disk, based on the average frequency */ + int64_t nphot[NRINGS]; /**< The number of disk photons created in each annulus */ + int64_t nhit[NRINGS]; /**< The number of photons which hit each annulus */ }; extern struct xdisk disk, qdisk; /**< disk defines zones in the disk which in a specified frequency band emit equal amounts of radiation. disk gets reinitialized whenever the frequency interval of interest diff --git a/source/run.c b/source/run.c index 2a9c7791e..93d8c92a6 100644 --- a/source/run.c +++ b/source/run.c @@ -48,8 +48,8 @@ int calculate_ionization (restart_stat) int restart_stat; { - int n, nn; - double zz, z_abs_all, z_abs_all_orig, z_orig[N_ISTAT], z_abs[N_ISTAT], z_else, z_else_orig, ztot; + int nn; + double zz, z_abs_all, z_abs_all_orig, z_orig[N_ISTAT], z_abs[N_ISTAT], z_else, z_else_orig; double radiated[20], radiated_orig[20]; int nphot_istat[N_ISTAT]; WindPtr w; @@ -146,7 +146,6 @@ calculate_ionization (restart_stat) geo.n_ioniz = 0.0; geo.cool_tot_ioniz = 0.0; - ztot = 0.0; /* ztot is the luminosity of the disk multipled by the number of cycles, which is used by save_disk_heating */ if (!geo.wind_radiation || (geo.wcycle == 0 && geo.run_type != RUN_TYPE_PREVIOUS)) iwind = -1; /* Do not generate photons from wind */ @@ -190,11 +189,9 @@ calculate_ionization (restart_stat) geo.lum_star_back = 0; geo.lum_disk_back = 0; + /* Prepare qdisk for recording photon pages */ + qdisk_reinit (p); - for (n = 0; n < NRINGS; n++) - { - qdisk.heat[n] = qdisk.nphot[n] = qdisk.w[n] = qdisk.ave_freq[n] = 0; - } zz = 0.0; @@ -205,7 +202,6 @@ calculate_ionization (restart_stat) Log ("!!python: Total photon luminosity before transphot %18.12e\n", zz); Log_flush (); - ztot += zz; /* Total luminosity in all cycles, used for calculating disk heating */ /* kbf_need determines how many & which bf processes one needs to considere. It was introduced * as a way to speed up the program. It has to be recalculated evey time one changes @@ -323,12 +319,12 @@ calculate_ionization (restart_stat) #endif if (geo.disk_type != DISK_NONE) { - qdisk_save (files.disk, ztot); + qdisk_save (files.disk, 1); if (modes.make_tables) { strcpy (dummy, ""); sprintf (dummy, "diag_%.100s/%.100s.%02d.disk.diag", files.root, files.root, geo.wcycle + 1); - qdisk_save (dummy, ztot); + qdisk_save (dummy, 0); } } diff --git a/source/templates.h b/source/templates.h index 4b7f4a7f0..e90e95a94 100644 --- a/source/templates.h +++ b/source/templates.h @@ -163,7 +163,8 @@ double disk_colour_correction(double t); /* disk_init.c */ double disk_init(double rmin, double rmax, double m, double mdot, double freqmin, double freqmax, int ioniz_or_extract, double *ftot); int qdisk_init(double rmin, double rmax, double m, double mdot); -int qdisk_save(char *diskfile, double ztot); +int qdisk_reinit(PhotPtr p); +int qdisk_save(char *diskfile, int ichoice); int read_non_standard_disk_profile(char *tprofile); /* disk_photon_gen.c */ int photo_gen_disk(PhotPtr p, double weight, double f1, double f2, int spectype, int istart, int nphot); From 167d2cc7d94e8ed7ae6fdecd9415363ddc964cd2 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Mon, 14 Oct 2024 09:32:22 -0500 Subject: [PATCH 17/23] Small formatting change --- source/disk_init.c | 2 +- source/disk_photon_gen.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/disk_init.c b/source/disk_init.c index 2411e0571..643ab9001 100644 --- a/source/disk_init.c +++ b/source/disk_init.c @@ -509,7 +509,7 @@ qdisk_save (diskfile, ichoice) ttot = pow (qdisk.t[n], 4) + pow (theat, 4); ttot = pow (ttot, 0.25); fprintf (qptr, - "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %10lld %8.3e %10lld %8.3e %8.3e %8.3e %10.3e %8.3e\n", + "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %10ld %8.3e %10ld %8.3e %8.3e %8.3e %10.3e %8.3e\n", qdisk.r[n], qdisk.v[n], zdisk (qdisk.r[n]), qdisk.t[n], qdisk.g[n], qdisk.emit[n], qdisk.nphot[n], qdisk.heat[n], qdisk.nhit[n], ratio[n], theat, qdisk.t_hit[n], qdisk.w[n], ttot); } diff --git a/source/disk_photon_gen.c b/source/disk_photon_gen.c index 225d13add..5c9d88291 100644 --- a/source/disk_photon_gen.c +++ b/source/disk_photon_gen.c @@ -277,7 +277,7 @@ disk_photon_summary (filename, mode) x = 0; else x = disk.nphot[n] / (disk.r[n] - disk.r[n - 1]); - fprintf (ptr, "%d %8.2e %8.2e %8lld %8.2e %8lld\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); + fprintf (ptr, "%d %8.2e %8.2e %8ld %8.2e %8ld\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); } fclose (ptr); From 8e6278190801468b629c81335e96ec2cfe711d95 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Mon, 14 Oct 2024 10:55:23 -0500 Subject: [PATCH 18/23] Removed int64_t not needed once other errors corrected. --- source/disk_init.c | 2 +- source/disk_photon_gen.c | 2 +- source/python.h | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/source/disk_init.c b/source/disk_init.c index 643ab9001..4b41a29ad 100644 --- a/source/disk_init.c +++ b/source/disk_init.c @@ -509,7 +509,7 @@ qdisk_save (diskfile, ichoice) ttot = pow (qdisk.t[n], 4) + pow (theat, 4); ttot = pow (ttot, 0.25); fprintf (qptr, - "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %10ld %8.3e %10ld %8.3e %8.3e %8.3e %10.3e %8.3e\n", + "%9.4e %9.4e %0.4e %8.3e %8.3e %8.3e %10d %8.3e %10d %8.3e %8.3e %8.3e %10.3e %8.3e\n", qdisk.r[n], qdisk.v[n], zdisk (qdisk.r[n]), qdisk.t[n], qdisk.g[n], qdisk.emit[n], qdisk.nphot[n], qdisk.heat[n], qdisk.nhit[n], ratio[n], theat, qdisk.t_hit[n], qdisk.w[n], ttot); } diff --git a/source/disk_photon_gen.c b/source/disk_photon_gen.c index 5c9d88291..826f3d51d 100644 --- a/source/disk_photon_gen.c +++ b/source/disk_photon_gen.c @@ -277,7 +277,7 @@ disk_photon_summary (filename, mode) x = 0; else x = disk.nphot[n] / (disk.r[n] - disk.r[n - 1]); - fprintf (ptr, "%d %8.2e %8.2e %8ld %8.2e %8ld\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); + fprintf (ptr, "%d %8.2e %8.2e %8d %8.2e %8d\n", n, disk.r[n], disk.t[n], disk.nphot[n], x, disk.nhit[n]); } fclose (ptr); diff --git a/source/python.h b/source/python.h index 3f4598060..84bdb9bef 100644 --- a/source/python.h +++ b/source/python.h @@ -760,8 +760,8 @@ struct xdisk double ave_freq[NRINGS]; /**< The flux weighted average of frequency of photons hitting each annulus */ double w[NRINGS]; /**< The weight relative to a BB of the photons that hit the disk */ double t_hit[NRINGS]; /**< The effective T of photons hitting the disk, based on the average frequency */ - int64_t nphot[NRINGS]; /**< The number of disk photons created in each annulus */ - int64_t nhit[NRINGS]; /**< The number of photons which hit each annulus */ + int nphot[NRINGS]; /**< The number of disk photons created in each annulus */ + int nhit[NRINGS]; /**< The number of photons which hit each annulus */ }; extern struct xdisk disk, qdisk; /**< disk defines zones in the disk which in a specified frequency band emit equal amounts of radiation. disk gets reinitialized whenever the frequency interval of interest From f1c2fe7fc7e7f8ec6d2d844633c432a1c5bffeeb Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Tue, 15 Oct 2024 07:24:58 -0500 Subject: [PATCH 19/23] Begin removal of matrix_ion2.c --- source/Makefile | 2 +- source/matrix_ion2.c | 437 ------------------------------------------- 2 files changed, 1 insertion(+), 438 deletions(-) delete mode 100644 source/matrix_ion2.c diff --git a/source/Makefile b/source/Makefile index c0d0b1e9d..33efedaf7 100644 --- a/source/Makefile +++ b/source/Makefile @@ -201,7 +201,7 @@ python_source = agn.c anisowind.c atomic_extern_init.c atomicdata.c atomicdata_i extract.c frame.c gradv.c gridwind.c homologous.c hydro_import.c import.c \ import_calloc.c import_cylindrical.c import_rtheta.c import_spherical.c ionization.c \ janitor.c knigge.c levels.c lines.c macro_accelerate.c macro_gen_f.c macro_gov.c \ - matom.c matom_diag.c matrix_cpu.c matrix_ion.c matrix_ion2.c models_extern_init.c para_update.c \ + matom.c matom_diag.c matrix_cpu.c matrix_ion.c models_extern_init.c para_update.c \ parse.c partition.c paths.c phot_util.c photon2d.c photon_gen.c photon_gen_matom.c \ pi_rates.c python_extern_init.c radiation.c random.c rdpar.c rdpar_init.c recipes.c \ recomb.c resonate.c reverb.c roche.c rtheta.c run.c saha.c setup.c setup_disk.c setup_domains.c \ diff --git a/source/matrix_ion2.c b/source/matrix_ion2.c deleted file mode 100644 index feae57fad..000000000 --- a/source/matrix_ion2.c +++ /dev/null @@ -1,437 +0,0 @@ - -/***********************************************************/ -/** @file matrix_ion.c - * @author ksl - * @date January, 2018 - * - * @brief Contains routines which compute the ionization state using a matrix inversion technique - * - ***********************************************************/ - -#include -#include -#include -#include - -#include -#include "atomic.h" -#include "python.h" - - -/**********************************************************/ -/** - * @brief A matrix solver for the ionization state in a cell - * - * @param [in,out] PlasmaPtr xplasma The plasma cell we are working on - * @param [in] int mode What type of model to use for J_nu - 1=fitted model 2=blackbody - * @return 0 if successful - * - * The abundances contained in xplasma are updated witht he results of - * the calculation. - * - * @details - * modes: - * 1 - use a power law and or exponential model for J - * 2 - use a dilute BB model for J defined by t_r and w - * The general structure of matrix_ion_populations is as follows: - * First we compute all of the rates that we have data for linking all ionization stages - * We make an initial guess at the electron density - * We attempt to solve the ionization rates matrix, and produce a new guess at the electron density - * We then calculate new rates, re-solve and proceed until the electron density converges. - * We solve the matrix equation A x = b for the vector x, where A is a square rate matrix linking - * each state with each other state (often containing zeros if the ions arent linked) b is a vector - * containing the total density of each element arranged in a certain way (see later) and - * x is a vector containing the relative ion populations. We invert A to get x=bA^-1 - * - * ### Notes ### - * Uses a relative abundance scheme, in order to reduce large number issues - * - * Various parameters for the calculation, and in particular the t_e are passed - * via the PlasmaPtr - * - **********************************************************/ - -int -matrix_ion_populations2 (xplasma, mode) - PlasmaPtr xplasma; - int mode; - -{ - double elem_dens[NELEMENTS]; //The density of each element - int nn, mm, nrows; - double rate_matrix[nions][nions]; //The rate matrix that we are going to try and solve - double newden[NIONS]; //A temporary array to hold our intermediate solutions - double nh, nh1, nh2, t_e; - double xne, xxne, xxxne; //Various stores for intermediate guesses at electron density - double b_temp[nions]; //The b matrix - double *b_data, *a_data; //These arrays are allocated later and sent to the matrix solver - double *populations; //This array is allocated later and is retrieved from the matrix solver - int matrix_err, niterate; //counters for errors and the number of iterations we have tried to get a converged electron density - double xnew; - int xion[nions]; // This array keeps track of what ion is in each line - double pi_rates[nions]; //photoionization rate coefficients - double rr_rates[nions]; //radiative recombination rate coefficients - double inner_rates[n_inner_tot]; //This array contains the rates for each of the inner shells. Where they go to requires the electron yield array - - nh1 = nh2 = 0; - - /* Copy some quantities from the cell into local variables */ - - nh = xplasma->rho * rho2nh; // The number density of hydrogen ions - computed from density - - // t_e = xplasma->t_e; // The electron temperature in the cell - used for collisional processes - xplasma->t_e_old = xplasma->t_e; - - - // t_e = calc_te (xplasma, 0.7 * xplasma->t_e, 1.3 * xplasma->t_e); - Log ("xxxxxx before calc_te %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, - xplasma->cool_tot); - t_e = calc_te (xplasma, 0.95 * xplasma->t_e, 1.05 * xplasma->t_e); - Log ("xxxxxx after calc_te %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, - xplasma->cool_tot); - xplasma->t_e = t_e; - zero_emit (xplasma->t_e); //Get the heating and cooling rates correctly for the new temperature - - /* We now calculate the total abundances for each element to allow us to use fractional abundances */ - - /* Belt and braces, we zero our array */ - for (mm = 0; mm < ion[nions - 1].z + 1; mm++) - { - elem_dens[mm] = 0.0; - } - - /* Now we populate the elemental abundance array */ - for (mm = 0; mm < nions; mm++) - { - elem_dens[ion[mm].z] = elem_dens[ion[mm].z] + xplasma->density[mm]; - } - - /* Dielectronic recombination, collisional ionization coefficients, three body recombination and - charge_exchange rate coefficients depend only on electron temperature, calculate them now - - they will not change they are all stored in global arrays */ - - compute_dr_coeffs (t_e); - compute_di_coeffs (t_e); - compute_qrecomb_coeffs (t_e); - compute_ch_ex_coeffs (t_e); - - /* In the following loop, over all ions in the simulation, we compute the radiative recombination rates, and photionization - rates OUT OF each ionization stage. The PI rates are calculated either using the modelled mean intensity in a cell, or - using the dilute blackbody approximation, depending on which mode we are in. At the same time, we copy the ion densities - from the plasma structure into a local array. We will only overwrite the numbers in the structure if we believe the - results are an improvement on what is there. */ - - for (mm = 0; mm < nions; mm++) - { - newden[mm] = xplasma->density[mm] / elem_dens[ion[mm].z]; // newden is our local fractional density array - xion[mm] = mm; // xion is an array we use to track which ion is in which row of the matrix - if (mm != ele[ion[mm].nelem].firstion) // We can recombine since we are not in the first ionization stage - { - rr_rates[mm] = total_rrate (mm, xplasma->t_e); // radiative recombination rates - } - if (ion[mm].istate != ele[ion[mm].nelem].istate_max) // we can photoionize, since we are not in the highest ionization state - { - if (mode == NEBULARMODE_MATRIX_BB) - { - pi_rates[mm] = calc_pi_rate (mm, xplasma, 2, 1); // PI rate for the BB model - } - else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL) - { - pi_rates[mm] = calc_pi_rate (mm, xplasma, 1, 1); // PI rate for an explicit spectral model - } - else if (mode == NEBULARMODE_MATRIX_MULTISHOT) - { - pi_rates[mm] = calc_pi_rate (mm, xplasma, 1, 1); // PI rate for an explicit spectral model - } - else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) - { - pi_rates[mm] = xplasma->ioniz[mm] / xplasma->density[mm]; // PI rate logged during the photon passage - } - else - { - Error ("matrix_ion_populations2: Unknown mode %d\n", mode); - Exit (0); - } - } - - } - - - /* The next loop generates the inner shell ionization rates, if they are present in the atomic data and - we wish to compute auger ionizaion rates. This only computes the rates out of each ion, we also need to - consult the electron yield array if we are to compute the change in electron number */ - - if (geo.auger_ionization == 1) - { - for (mm = 0; mm < n_inner_tot; mm++) - { - if (mode == NEBULARMODE_MATRIX_BB) - { - inner_rates[mm] = calc_pi_rate (mm, xplasma, 2, 2); - } - else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL) - { - inner_rates[mm] = calc_pi_rate (mm, xplasma, 1, 2); - } - else if (mode == NEBULARMODE_MATRIX_MULTISHOT) - { - inner_rates[mm] = calc_pi_rate (mm, xplasma, 1, 2); - } - else if (mode == NEBULARMODE_MATRIX_ESTIMATORS) //We are using estimators - so we will need to reorder the rates from freq order to cross section order to match the electron yields. This takes time, so we only do it if we need to. - { - for (nn = 0; nn < n_inner_tot; nn++) - { - if (inner_cross[mm].nion == inner_cross_ptr[nn]->nion && inner_cross[mm].freq[0] == inner_cross_ptr[nn]->freq[0]) //Check for a match - { - inner_rates[mm] = xplasma->inner_ioniz[nn] / xplasma->density[inner_cross_ptr[nn]->nion]; - } - } - } - } - } - - /* This next line sets the partition function for each ion. This has always been the place here python calculates the - partition functions and sets the level densities for each ion. It needs to be done, or other parts of the code which rely - on sensible level populations don't work properly. In the case of the dilute blackbody, the code works well, however we do - not currently (v78) have a procedure to calucate the levels for a spectral model case. We therefore call partition - functions with mode 4 - this is a special mode which forces the ions to be in the ground state. This is reasonable for a - radiation dominated plasma, since any excitations into higher states will quickly be followed by radative de-excitation. - We should really do better though... */ - - if (mode == NEBULARMODE_MATRIX_BB) - { - partition_functions (xplasma, NEBULARMODE_ML93); // We use t_r and the radiative weight - } - else if (mode == NEBULARMODE_MATRIX_SPECTRALMODEL || mode == NEBULARMODE_MATRIX_ESTIMATORS) - { - partition_functions (xplasma, NEBULARMODE_LTE_GROUND); // Set to ground state - } - else if (mode == NEBULARMODE_MATRIX_MULTISHOT) - { - partition_functions (xplasma, NEBULARMODE_LTE_GROUND); // Set to ground state - } - - /* Next we need to obtain an initial guess for the electron density. In the past this has been done by calculating the - hydrogen density directly from the Saha equation at the current electron temperature. In initial testing of this mode - - this seemed to be a little unstable. At the moment, we use the last ionization state to compute n_e. For the first - iteraion, this will just be calculated in LTE from the initial temperature structure of the wind, which will give almost - the same result as the original procedure, or for successive calculations, it should be a better guess. I've leftin the - original code, commented out... */ - - xne = xxne = xxxne = get_ne (xplasma->density); //Even though the abundances are fractional, we need the real electron density - - - /* xne is the current working number xxne */ - - - /* We are now going to iterate on the electron density - MAXITERATIONS is set in python.h and is currently (78) set to 200. - We would normally expect to converge much fater than this */ - - Log ("xxxxxx before_loop %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, - xplasma->cool_tot); - niterate = 0; - while (niterate < MAXITERATIONS) - { - - /* In order to compute charge exchange reactions, we need the number density of neutral and ionized hydrogen - the current ethod of matrix solving does not allow is to link ions multiplicatively so we need to compute - these outside the loop. If hydrogen were not dominant - this could cause iddues since the ionization(recombination) - of a metal or helium via this process would cause an identical recombination(ionization) of hydrogen. The loop below - is a bit belt and braces, since one would almost always expect hydrogen to be ion=0 and 1 */ - - for (nn = 0; nn < nions; nn++) - { - if (ion[nn].z == 1 && ion[nn].istate == 1) - { - nh1 = newden[nn] * elem_dens[ion[nn].z]; - } - if (ion[nn].z == 1 && ion[nn].istate == 2) - { - nh2 = newden[nn] * elem_dens[ion[nn].z]; - } - } - - - populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp, xne, nh1, nh2); - - - /* The array is now fully populated, and we can begin the process of solving it */ - - nrows = nions; /* This is a placeholder, we may end up removing rows and columns that have no rates (or very - low rates) connecting them to other rows. This may improve stability but will need to be - done carefully */ - - /* Here we solve the matrix equation M x = b, where x is our vector containing level populations as a fraction w.r.t the - whole element. The actual LU decomposition- the process of obtaining a solution - is done by the routine solve_matrix() */ - - a_data = (double *) calloc (sizeof (double), nrows * nrows); - populations = (double *) calloc (nrows, sizeof (double)); - - /* This b_data column matrix is the total number density for each element, placed into the row which relates to the neutral - ion. This matches the row in the rate matrix which is just 1 1 1 1 for all stages. NB, we could have chosen any line for - this. */ - - b_data = (double *) calloc (nrows, sizeof (double)); - - /* We now copy our rate matrix into the prepared matrix */ - for (mm = 0; mm < nrows; mm++) - { - for (nn = 0; nn < nrows; nn++) - { - a_data[mm * nrows + nn] = rate_matrix[mm][nn]; /* row-major */ - } - } - - for (nn = 0; nn < nrows; nn++) - { - b_data[nn] = b_temp[nn]; - } - - matrix_err = solve_matrix (a_data, b_data, nrows, populations, xplasma->nplasma); - - if (matrix_err) - { - Error ("matrix_ion_populations2: %s\n", get_matrix_error_string (matrix_err)); - } - - /* free memory */ - free (a_data); - free (b_data); - - if (matrix_err == 4) - { - free (populations); - return (-1); - } - - /* Calculate level populations for macro-atoms */ - if (geo.macro_ioniz_mode == MACRO_IONIZ_MODE_ESTIMATORS) - { - int mp_err = macro_pops (xplasma, xne); - if (mp_err != EXIT_SUCCESS) - { - return -1; - } - } - - /* We now have the populations of all the ions stored in the matrix populations. We copy this data into the newden array - which will temperarily store all the populations. We wont copy this to the plasma structure until we are sure thatwe - have made things better. We just loop over all ions, and copy. The complexity in the loop is to future proof us against - the possibility that there are some ions that are not included in the matrix scheme becuse there is no route into or - out of it. */ - - for (nn = 0; nn < nions; nn++) - { - newden[nn] = 0.0; - - /* if the ion is being treated by macro_pops then use the populations just computed */ - if ((ion[nn].macro_info == TRUE) && (geo.macro_simple == FALSE) && (geo.macro_ioniz_mode == MACRO_IONIZ_MODE_ESTIMATORS) - && (modes.no_macro_pops_for_ions == FALSE)) - { - newden[nn] = xplasma->density[nn] / elem_dens[ion[nn].z]; - } - - /* if the ion is "simple" then find it's calculated ionization state in populations array */ - else - { - - for (mm = 0; mm < nrows; mm++) // inner loop over the elements of the population array - { - if (xion[mm] == nn) // if this element contains the population of the ion is question - { - newden[nn] = populations[mm]; // get the population - } - } - } - - if (newden[nn] < DENSITY_MIN) // this wil also capture the case where population doesnt have a value for this ion - newden[nn] = DENSITY_MIN; - } - free (populations); - - -/* We need to get the 'true' new electron density so we need to do a little loop here to compute it */ - - - xnew = 0.0; - for (nn = 0; nn < nions; nn++) - { - xnew += newden[nn] * (ion[nn].istate - 1) * elem_dens[ion[nn].z]; - } - - - - if (xnew < DENSITY_MIN) - xnew = DENSITY_MIN; /* fudge to keep a floor on ne */ - - -// if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR || xnew < 1.e-6) /* We have converged, or have the situation where we -// have a neutral plasma */ -// if (fabs ((xne - xnew) / (xnew)) < FRACTIONAL_ERROR) /* We have converged, or have the situation where we -// have a neutral plasma */ - if (fabs ((xne - xnew) / (xnew)) < 0.01) /* We have converged, or have the situation where we - have a neutral plasma */ - { - break; /* Break out of the while loop - we have finished our iterations */ - } - - Log ("xxxxxx %3d xnew %.2e xne %.2e -> %.2e \n", niterate, xnew, xne, (xnew + xne) / 2.); - xne = xxxne = (xnew + xne) / 2.; /* New value of ne */ - - niterate++; - - - if (niterate == MAXITERATIONS) - { - Error ("matrix_ion_populations2: failed to converge for cell %i t %e nh %e xnew %e\n", xplasma->nplasma, t_e, nh, xnew); - - - for (nn = 0; nn < geo.nxfreq; nn++) - { - Log - ("numin= %e (%e) numax= %e (%e) Model= %2d PL_log_w= %e PL_alpha= %e Exp_w= %e EXP_temp= %e\n", - xplasma->fmin_mod[nn], geo.xfreq[nn], xplasma->fmax_mod[nn], - geo.xfreq[nn + 1], xplasma->spec_mod_type[nn], - xplasma->pl_log_w[nn], xplasma->pl_alpha[nn], xplasma->exp_w[nn], xplasma->exp_temp[nn]); - } - - Error ("matrix_ion_populations2: xxne %e theta %e\n", xxne); - - return (-1); /* If we get to MAXITERATIONS, we return without copying the new populations into plasma */ - } - } /* This is the end of the iteration loop */ - - - xplasma->ne = xnew; - for (nn = 0; nn < nions; nn++) - { - /* If statement added here to suppress interference with macro populations */ - if (ion[nn].macro_info == FALSE || geo.macro_ioniz_mode == MACRO_IONIZ_MODE_NO_ESTIMATORS || geo.macro_simple == TRUE) - { - xplasma->density[nn] = newden[nn] * elem_dens[ion[nn].z]; //We return to absolute densities here - } - if ((sane_check (xplasma->density[nn])) || (xplasma->density[nn] < 0.0)) - Error ("matrix_ion_populations2: ion %i has population %8.4e in cell %i\n", nn, xplasma->density[nn], xplasma->nplasma); - } - - xplasma->ne = get_ne (xplasma->density); - - if (n_charge_exchange > 0) - { - xplasma->heat_ch_ex = ch_ex_heat (&wmain[xplasma->nwind], xplasma->t_e); //Compute the charge exchange heating - xplasma->heat_tot += xplasma->heat_ch_ex; - } - - /*We now need to populate level densities in order to later calculate line emission (for example). - We call partition functions to do this. At present, we do not have a method for accurately computing - level populations for modelled mean intensities. We get around this by setting all level populations - other than ground state to zero - this is a pretty good situation for very dilute radiation fields, - which is what we are normally looking at, however one should really do better in the long term. - */ - - partition_functions (xplasma, NEBULARMODE_LTE_GROUND); - Log ("xxxxxx after At_end %d t_e %.1e n_e %.1e heat_tot %.1e %.1e\n", geo.wcycle, xplasma->t_e, xplasma->ne, xplasma->heat_tot, - xplasma->cool_tot); - - return (0); -} From 79efbc46fa95970f5f853f754efeb28f66d5c07f Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Tue, 15 Oct 2024 07:31:45 -0500 Subject: [PATCH 20/23] Remove uncompleted mulitshot mode --- source/ionization.c | 13 ------------- source/python.h | 4 ---- source/saha.c | 4 ---- source/setup.c | 3 +-- 4 files changed, 1 insertion(+), 23 deletions(-) diff --git a/source/ionization.c b/source/ionization.c index 4d1697899..127cc3e92 100644 --- a/source/ionization.c +++ b/source/ionization.c @@ -125,19 +125,6 @@ ion_abundances (PlasmaPtr xplasma, int mode) convergence (xplasma); } - else if (mode == IONMODE_MATRIX_MULTISHOT) - { -/* New mode to carry more extensive search for best fit - Note that it unclear that we could not simplify a lot of this - routine but what is done here is most analagous to the other modes - excetp we avoid one_shot*/ - spectral_estimators (xplasma); - update_old_plasma_variables (xplasma); - // ireturn = one_shot (xplasma, mode); - ireturn = nebular_concentrations (xplasma, mode); - - convergence (xplasma); - } else { Error ("ion_abundances: Could not calculate abundances for mode %d\n", mode); diff --git a/source/python.h b/source/python.h index 84bdb9bef..f6c0a4d4e 100644 --- a/source/python.h +++ b/source/python.h @@ -1224,10 +1224,6 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; #define IONMODE_MATRIX_BB 8 /**< matrix solver BB model */ #define IONMODE_MATRIX_SPECTRALMODEL 9 /**< matrix solver spectral model based on power laws */ #define IONMODE_MATRIX_ESTIMATORS 10 /**< matrix solver spectral model based on power laws */ -#define IONMODE_MATRIX_MULTISHOT 11 /**< matrix solver spectral model based on power laws which - * updates T_e multiple times before arriving at a final - * solution - */ // and the corresponding modes in nebular_concentrations #define NEBULARMODE_TR 0 /**< LTE using t_r */ diff --git a/source/saha.c b/source/saha.c index 21f5b77c8..2da06f038 100644 --- a/source/saha.c +++ b/source/saha.c @@ -80,10 +80,6 @@ nebular_concentrations (xplasma, mode) m = matrix_ion_populations (xplasma, mode); } - else if (mode == NEBULARMODE_MATRIX_MULTISHOT) - { - m = matrix_ion_populations2 (xplasma, mode); - } else { Error ("nebular_concentrations: Unknown mode %d\n", mode); diff --git a/source/setup.c b/source/setup.c index 6d786e122..0cd7fec2e 100644 --- a/source/setup.c +++ b/source/setup.c @@ -714,8 +714,7 @@ init_ionization () strcpy (answer, "matrix_bb"); geo.ioniz_mode = - rdchoice ("Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est,matrix_te)", "0,3,1,4,2,8,9,10,11", - answer); + rdchoice ("Wind.ionization(on.the.spot,ML93,LTE_tr,LTE_te,fixed,matrix_bb,matrix_pow,matrix_est)", "0,3,1,4,2,8,9,10", answer); if (geo.ioniz_mode == IONMODE_FIXED) { From 236b2601efe2231e993ea4a937def2b2e4c59748 Mon Sep 17 00:00:00 2001 From: "Knox S. Long" Date: Wed, 16 Oct 2024 09:26:34 -0500 Subject: [PATCH 21/23] Restore fractional error to 0.03 --- source/python.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/python.h b/source/python.h index f6c0a4d4e..4b7424e9a 100644 --- a/source/python.h +++ b/source/python.h @@ -1209,7 +1209,7 @@ extern int size_Jbar_est, size_gamma_est, size_alpha_est; //These constants are used in the various routines which compute ionization state #define SAHA 4.82907e15 /**< 2* (2.*PI*MELEC*k)**1.5 / h**3 (Calculated in constants) */ #define MAXITERATIONS 200 /**< /The number of loops to do to try to converge in ne */ -#define FRACTIONAL_ERROR 0.001 /**< The change in n_e which causes a break out of the loop for ne */ +#define FRACTIONAL_ERROR 0.03 /**< The change in n_e which causes a break out of the loop for ne */ #define THETAMAX 1e4 /**< Used in initial calculation of n_e */ #define MIN_TEMP 100. /**< ??? this is another minimum temperature, which is * used in saha.c and variable_temperature.c ) From 1062c87c564e59f8a079c0c9cbc39c5b93b5dc6d Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Thu, 17 Oct 2024 15:09:49 +0100 Subject: [PATCH 22/23] make radiation temperature calculation depend on BAND_CORRECTED_TRAD --- source/estimators_simple.c | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/source/estimators_simple.c b/source/estimators_simple.c index 1e309030d..2718da15b 100644 --- a/source/estimators_simple.c +++ b/source/estimators_simple.c @@ -574,6 +574,10 @@ estimate_temperature_from_mean_frequency (double mean_nu_target, double nu_min, * **********************************************************/ +/* switch to choose how to calculate T_r. This could be made to depend on the ionization mode + if desired, or it may be that TRUE is the desired behaviour long-term */ +#define BAND_CORRECTED_TRAD FALSE + int normalise_simple_estimators (xplasma) PlasmaPtr xplasma; @@ -606,9 +610,17 @@ normalise_simple_estimators (xplasma) xplasma->t_r_old = xplasma->t_r; // Store the previous t_r in t_r_old immediately before recalculating - radiation_temperature = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); - radiation_temperature = xplasma->t_r = - estimate_temperature_from_mean_frequency (xplasma->ave_freq, xband.f1[0], xband.f2[xband.nbands - 1], radiation_temperature); + /* the method of calculation of the band corrected radiation temperature depends on + the flag BAND_CORRECTED_TRAD -- see issue #1097 */ + if (BAND_CORRECTED_TRAD == FALSE) + { + radiation_temperature = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); + } + else + { + radiation_temperature = xplasma->t_r = + estimate_temperature_from_mean_frequency (xplasma->ave_freq, xband.f1[0], xband.f2[xband.nbands - 1], radiation_temperature); + } xplasma->w = PI * xplasma->j / (STEFAN_BOLTZMANN * radiation_temperature * radiation_temperature * radiation_temperature * radiation_temperature); From 8d3e35acc4407886044da25f55d218b1ce63fbbc Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Thu, 17 Oct 2024 15:17:19 +0100 Subject: [PATCH 23/23] corrected --- source/estimators_simple.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/estimators_simple.c b/source/estimators_simple.c index 2718da15b..c85f55f98 100644 --- a/source/estimators_simple.c +++ b/source/estimators_simple.c @@ -614,7 +614,7 @@ normalise_simple_estimators (xplasma) the flag BAND_CORRECTED_TRAD -- see issue #1097 */ if (BAND_CORRECTED_TRAD == FALSE) { - radiation_temperature = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); + radiation_temperature = xplasma->t_r = PLANCK * xplasma->ave_freq / (BOLTZMANN * 3.832); } else {