diff --git a/.travis.yml b/.travis.yml index e44143447..d2242c57f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -49,11 +49,14 @@ script: - cd ../examples/travis - Setup_Py_Dir - py balmer_test > balmer_test.out - - cd .. + - cd ../core/ - Setup_Py_Dir - py -i cv_macro_benchmark - py -i cv_standard - py -i -d sv_detailedmode - py -i fiducial_agn - py -i 1d_sn + - cd ../beta/ - py -i agn_ss_2010_modela.pf + - py -i ulx1 + - py -i ngc5548 diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 000000000..0af5f8b93 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,26 @@ +### examples + +This folder contains a number of standard parameter files, used for testing and porting to the user. + +### core + +Core files which are standard tests + +* cv_standard.pf -- Standard CV model (Shlosman & Vitello 1993; Long & Knigge 2002). +* fiducial_agn.pf -- QSO model from Higginbottom et al. (2013). +* 1d_sn.pf -- simple Supernovae model as presented in Kerzendorf & Sim (2013). +* star.pf -- simple spherical stellar wind model + +More performance intensive core matom models: + +* cv_macro_benchmark.pf -- Macro-atom CV run from Matthews et al. (2015). +* m16_agn.pf -- Macro-atom AGN run from Matthews et al. (2016). + + +### beta + +Files of astrophysical interest not currently fully tested. + +* ngc5548.pf -- designed to represent rough params for NGC 5548 +* ulx1.pf -- First attempt at a ULX model +* agn_ss_2010_modela.pf -- AGN model A from Sim et al. (2010). diff --git a/examples/agn_ss_2010_modela.pf b/examples/beta/agn_ss_2010_modela.pf similarity index 100% rename from examples/agn_ss_2010_modela.pf rename to examples/beta/agn_ss_2010_modela.pf diff --git a/examples/beta/ngc5548.pf b/examples/beta/ngc5548.pf new file mode 100644 index 000000000..83e3abfbc --- /dev/null +++ b/examples/beta/ngc5548.pf @@ -0,0 +1,59 @@ +Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell) 0 +Atomic_data data/standard78 +photons_per_cycle 20000000 +Ionization_cycles 30 +spectrum_cycles 20 +Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 1 +Wind.dim.in.x_or_r.direction 100 +Wind.dim.in.z_or_theta.direction 100 +Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow,8=matrix_bb,9=matrix_pow) 7 +Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering) 3 +Thermal_balance_options(0=everything.on,1=no.adiabatic) 0 +System_type(0=star,1=binary,2=agn) 2 +disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk) 1 +Disk_radiation(y=1) 1 +Wind_radiation(y=1) 1 +QSO_BH_radiation(y=1) 1 +Rad_type_for_disk(0=bb,1=models)_to_make_wind 0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_to_make_wind 3 +mstar(msol) 1e7 +rstar(cm) 8.85667e+12 +disk.mdot(msol/yr) 0.02 +Disk.illumination.treatment(0=no.rerad,1=high.albedo,2=thermalized.rerad,3=analytic) 0 +Disk.temperature.profile(0=standard;1=readin) 0 +disk.radmax(cm) 1e15 +lum_agn(ergs/s) 1e43 +agn_power_law_index -0.9 +wind.radmax(cm) 1e17 +wind.t.init 1e5 +wind.mdot(msol/yr) 5 +sv.diskmin(wd_rad) 50 +sv.diskmax(wd_rad) 100 +sv.thetamin(deg) 70 +sv.thetamax(deg) 82 +sv.mdot_r_exponent 0 +sv.v_infinity(in_units_of_vescape 1 +sv.acceleration_length(cm) 1e16 +sv.acceleration_exponent 1.0 +wind.filling_factor(1=smooth,<1=clumpted) 1 +Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum 0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_in_final_spectrum 3 +spectrum_wavemin 200 +spectrum_wavemax 2600 +no_observers 8 +angle(0=pole) 20 +angle(0=pole) 40 +angle(0=pole) 60 +angle(0=pole) 70 +angle(0=pole) 75 +angle(0=pole) 80 +angle(0=pole) 85 +angle(0=pole) 89 +live.or.die(0).or.extract(anything_else) 1 +Select_specific_no_of_scatters_in_spectra(y/n) n +Select_photons_by_position(y/n) n +spec.type(flambda(1),fnu(2),basic(other) 1 +Extra.diagnostics(0=no) 0 +Use.standard.care.factors(1=yes) 1 +reverb.type(0=None,1=Photon,2=Wind) 0 +Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic) 7 diff --git a/examples/beta/ulx1.pf b/examples/beta/ulx1.pf new file mode 100644 index 000000000..437aaa407 --- /dev/null +++ b/examples/beta/ulx1.pf @@ -0,0 +1,59 @@ +Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell) 0 +Atomic_data std_79 +photons_per_cycle 20000000 +Ionization_cycles 30 +spectrum_cycles 20 +Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 1 +Wind.dim.in.x_or_r.direction 100 +Wind.dim.in.z_or_theta.direction 100 +Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow,8=matrix_bb,9=matrix_pow) 9 +Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering) 5 +Thermal_balance_options(0=everything.on,1=no.adiabatic) 0 +System_type(0=star,1=binary,2=agn) 2 +disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk) 0 +Disk_radiation(y=1) 0 +Wind_radiation(y=1) 1 +QSO_BH_radiation(y=1) 1 +Rad_type_for_disk(0=bb,1=models)_to_make_wind 0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_to_make_wind 3 +mstar(msol) 10.0 +rstar(cm) 8.85667e+5 +disk.mdot(msol/yr) 1.25e-8 +Disk.illumination.treatment(0=no.rerad,1=high.albedo,2=thermalized.rerad,3=analytic) 0 +Disk.temperature.profile(0=standard;1=readin) 0 +disk.radmax(cm) 1e7 +lum_agn(ergs/s) 1e39 +agn_power_law_index -0.9 +wind.radmax(cm) 1e10 +wind.t.init 1e6 +wind.mdot(msol/yr) 5e-6 +sv.diskmin(wd_rad) 100 +sv.diskmax(wd_rad) 200 +sv.thetamin(deg) 70 +sv.thetamax(deg) 82 +sv.mdot_r_exponent 0 +sv.v_infinity(in_units_of_vescape 2 +sv.acceleration_length(cm) 1e8 +sv.acceleration_exponent 0.2 +wind.filling_factor(1=smooth,<1=clumpted) 1 +Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum 0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_in_final_spectrum 3 +spectrum_wavemin 2 +spectrum_wavemax 50 +no_observers 8 +angle(0=pole) 20 +angle(0=pole) 40 +angle(0=pole) 60 +angle(0=pole) 70 +angle(0=pole) 75 +angle(0=pole) 80 +angle(0=pole) 85 +angle(0=pole) 89 +live.or.die(0).or.extract(anything_else) 1 +Select_specific_no_of_scatters_in_spectra(y/n) n +Select_photons_by_position(y/n) n +spec.type(flambda(1),fnu(2),basic(other) 2 +Extra.diagnostics(0=no) 0 +Use.standard.care.factors(1=yes) 1 +reverb.type(0=None,1=Photon,2=Wind) 0 +Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic) 7 diff --git a/examples/1d_sn.pf b/examples/core/1d_sn.pf similarity index 97% rename from examples/1d_sn.pf rename to examples/core/1d_sn.pf index 67e029a70..4c2c54b8a 100644 --- a/examples/1d_sn.pf +++ b/examples/core/1d_sn.pf @@ -2,7 +2,7 @@ Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=homologous,7=yso, Atomic_data data/standard_sn_kurucz photons_per_cycle 80000 Ionization_cycles 20 -spectrum_cycles 1 +spectrum_cycles 10 Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 0 Wind.dim.in.x_or_r.direction 30 Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow,8=matrix_bb,9=matrix_pow) 3 diff --git a/examples/cv_macro_benchmark.pf b/examples/core/cv_macro_benchmark.pf similarity index 98% rename from examples/cv_macro_benchmark.pf rename to examples/core/cv_macro_benchmark.pf index 68c08e808..a386bc305 100644 --- a/examples/cv_macro_benchmark.pf +++ b/examples/core/cv_macro_benchmark.pf @@ -1,8 +1,8 @@ Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell) 0 Atomic_data data/h20_hetop_standard78 photons_per_cycle 100000 -Ionization_cycles 1 -spectrum_cycles 0 +Ionization_cycles 30 +spectrum_cycles 10 Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 1 Wind.dim.in.x_or_r.direction 30 Wind.dim.in.z_or_theta.direction 30 diff --git a/examples/cv_standard.pf b/examples/core/cv_standard.pf similarity index 100% rename from examples/cv_standard.pf rename to examples/core/cv_standard.pf diff --git a/examples/fiducial_agn.pf b/examples/core/fiducial_agn.pf similarity index 98% rename from examples/fiducial_agn.pf rename to examples/core/fiducial_agn.pf index dda1cfb7f..c25a275e9 100644 --- a/examples/fiducial_agn.pf +++ b/examples/core/fiducial_agn.pf @@ -1,6 +1,6 @@ Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell) 0 Atomic_data data/standard78 -photons_per_cycle 200000 +photons_per_cycle 2000000 Ionization_cycles 30 spectrum_cycles 20 Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 1 diff --git a/examples/core/m16_agn.pf b/examples/core/m16_agn.pf new file mode 100644 index 000000000..c5997ea04 --- /dev/null +++ b/examples/core/m16_agn.pf @@ -0,0 +1,72 @@ +Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell) 0.0 +Atomic_data data/h10_hetop_lohe1_standard78 +photons_per_cycle 20000000.0 +Ionization_cycles 20.0 +spectrum_cycles 5.0 +Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 1.0 +Wind.dim.in.x_or_r.direction 100.0 +Wind.dim.in.z_or_theta.direction 200.0 +Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow) 9.0 +Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering) 7.0 +Thermal_balance_options(0=everything.on,1=no.adiabatic) 0.0 +System_type(0=star,1=binary,2=agn) 2.0 +disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk) 1.0 +Disk_radiation(y=1) 1.0 +Wind_radiation(y=1) 1.0 +QSO_BH_radiation(y=1) 1.0 +Rad_type_for_disk(0=bb,1=models)_to_make_wind 0.0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_to_make_wind 3.0 +mstar(msol) 1000000000.0 +rstar(cm) 8.85667e+14 +disk.mdot(msol/yr) 5.0 +Disk.illumination.treatment(0=no.rerad,1=high.albedo,2=thermalized.rerad,3=analytic) 0.0 +Disk.temperature.profile(0=standard;1=readin) 0.0 +disk.radmax(cm) 1e+17 +lum_agn(ergs/s) 1e+45 +agn_power_law_index -0.9 +Torus(0=no,1=yes) 0.0 +wind.radmax(cm) 1e+20 +wind.t.init 100000.0 +wind.mdot(msol/yr) 5.0 +sv.diskmin(wd_rad) 50 +sv.diskmax(wd_rad) 100.0 +sv.thetamin(deg) 70 +sv.thetamax(deg) 82.0 +sv.mdot_r_exponent 0.0 +sv.v_infinity(in_units_of_vescape 1.0 +sv.acceleration_length(cm) 1e19 +sv.acceleration_exponent 0.5 +wind.filling_factor(1=smooth,<1=clumped) 0.01 +Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum 0.0 +Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_in_final_spectrum 3.0 +spectrum_wavemin 500.0 +spectrum_wavemax 7500.0 +no_observers 20.0 +angle(0=pole) 20.0 +angle(0=pole) 40.0 +angle(0=pole) 60.0 +angle(0=pole) 70.0 +angle(0=pole) 71.0 +angle(0=pole) 72.0 +angle(0=pole) 73.0 +angle(0=pole) 74.0 +angle(0=pole) 75.0 +angle(0=pole) 76.0 +angle(0=pole) 77.0 +angle(0=pole) 78.0 +angle(0=pole) 79.0 +angle(0=pole) 80.0 +angle(0=pole) 81.0 +angle(0=pole) 82.0 +angle(0=pole) 83.0 +angle(0=pole) 84.0 +angle(0=pole) 85.0 +angle(0=pole) 89.0 +live.or.die(0).or.extract(anything_else) 1.0 +Select_specific_no_of_scatters_in_spectra(y/n) n +Select_photons_by_position(y/n) n +spec.type(flambda(1),fnu(2),basic(other) 1.0 +reverb.type 0.0 +Extra.diagnostics(0=no) 0.0 +Use.standard.care.factors(1=yes) 1.0 +Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined), 7.0 diff --git a/examples/core/star.pf b/examples/core/star.pf new file mode 100644 index 000000000..fdb20ad83 --- /dev/null +++ b/examples/core/star.pf @@ -0,0 +1,46 @@ +Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=homologous,7=yso,8=elvis,9=shell) 1 +Atomic_data data/standard78 +photons_per_cycle 100000 +Ionization_cycles 20 +spectrum_cycles 20 +Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var) 0 +Wind.dim.in.x_or_r.direction 30 +Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,6=pairwise_bb,7=pairwise_pow) 3 +Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering) 3 +Thermal_balance_options(0=everything.on,1=no.adiabatic) 0 +System_type(0=star,1=binary,2=agn) 0 +disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk) 0 +Star_radiation(y=1) 1 +Disk_radiation(y=1) 0 +Boundary_layer_radiation(y=1) 0 +Wind_radiation(y=1) 1 +Rad_type_for_star(0=bb,1=models)_to_make_wind 0 +Model_file kurucz91/kurucz91.ls +mstar(msol) 52.5 +rstar(cm) 1.32e12 +tstar 42000 +disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk) 0 +Torus(0=no,1=yes) 0 +wind.radmax(cm) 1e13 +wind.t.init 40000 +stellar_wind_mdot(msol/yr) 5.1e-6 +stellar.wind.radmin(cm) 1.32e12 +stellar.wind_vbase(cm) 2e+07 +stellar.wind.v_infinity(cm) 2.25e+08 +stellar.wind.acceleration_exponent 1 +wind.filling_factor(1=smooth,<1=clumped) 1 +Rad_type_for_star(0=bb,1=models,2=uniform)_in_final_spectrum 0 +Model_file kurucz91/kurucz91.ls +spectrum_wavemin 850 +spectrum_wavemax 1850 +no_observers 1 +angle(0=pole) 45 +phase(0=inferior_conjunction) 0.5 +live.or.die(0).or.extract(anything_else) 1 +Select_specific_no_of_scatters_in_spectra(y/n) n +Select_photons_by_position(y/n) n +spec.type(flambda(1),fnu(2),basic(other) 1 +Extra.diagnostics(0=no) 0 +Use.standard.care.factors(1=yes) 1 +reverb.type 0 +Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined,5=cloudy_test,6=wide,7=AGN,8=logarithmic) 2 diff --git a/examples/sv_detailedmode.pf b/examples/core/sv_detailedmode.pf similarity index 100% rename from examples/sv_detailedmode.pf rename to examples/core/sv_detailedmode.pf diff --git a/examples/core/template_ionloop.pf b/examples/core/template_ionloop.pf new file mode 100644 index 000000000..52168b975 --- /dev/null +++ b/examples/core/template_ionloop.pf @@ -0,0 +1,47 @@ +Wind_type() 9 +Atomic_data data/standard78 +photons_per_cycle 100000 +Ionization_cycles 1 +spectrum_cycles 0 +Coord.system() 0 +Wind.dim.in.x_or_r.direction 4 +Wind_ionization() 7 +Line_transfer() 5 +Thermal_balance_options 0 +System_type() 2 +disk.type() 0 +Star_radiation(y=1) 0 +Disk_radiation(y=1) 0 +Wind_radiation(y=1) 0 +QSO_BH_radiation(y=1) 1 +Rad_type_for_star(0=bb,1=models)_to_make_wind 0 +Rad_type_for_agn()_to_make_wind 4 +mstar(msol) 0.8 +rstar(cm) 1e10 +tstar 1000000 +lum_agn(ergs/s) 1e+21 +agn_power_law_index -0.9 +low_energy_break(ev) 0.136 +high_energy_break(ev) 20000 +Torus(0=no,1=yes) 0 +wind.radmax(cm) 1.00000000001e11 +wind.t.init 10000 +shell_wind_mdot(msol/yr) 0.00472e-17 +shell.wind.radmin(cm) 1e11 +shell.wind_v_at_rmin(cm) 1.00000 +shell.wind.v_at_rmax(cm) 1.000010 +shell.wind.acceleration_exponent 1 +wind.filling_factor(1=smooth,<1=clumped) (1) 1 +spec.type(flambda(1),fnu(2),basic(other) 2 +reverb.type 0 +Use.standard.care.factors(1=yes) 1 +Photon.sampling.approach() 5 +Num.of.frequency.bands 3 +Lowest_energy_to_be_considered(eV) 0.0001 +Highest_energy_to_be_considered(eV) 100000000 +Band.boundary(eV) 2000 +Band.boundary(eV) 10000 +Band.minimum_fraction) 0.3 +Band.minimum_fraction) 0.4 +Band.minimum_fraction) 0.3 + diff --git a/py_progs/hydro_2_python.py b/py_progs/hydro_2_python.py index 653f23366..a149f7577 100755 --- a/py_progs/hydro_2_python.py +++ b/py_progs/hydro_2_python.py @@ -88,13 +88,18 @@ def get_hdf_data(fname): if name[0:4]=="Data": sds=hdf.select(name) long_name=sds.attributes()["long_name"] - short_name=long_name[:-17] + for i in range(len(sds.attributes()["long_name"])): + if long_name[i:i+2]=="AT": + junk=i-1 + short_name=long_name[:junk] data_sets.append([name,long_name,short_name]) #Get the time from the last long name - - time=float(long_name[-8:]) - + if long_name[junk+9:] != "********": + time=float(long_name[junk+9:]) + else: + time=0.0 + #Get the dimensions from the last data set dims=len((sds.info()[2])) @@ -122,10 +127,6 @@ def get_hdf_data(fname): #Loop over all the data sets in the hdf file - name each of the resulting dictionaries with the short name - - - - for i in range (len(data_sets)): print data_sets[i][2] sds=hdf.select(data_sets[i][0]) @@ -165,6 +166,9 @@ def get_hdf_data(fname): for i in range(len(x2)-1): theta_edge.append(theta_edge[-1]+dtheta) dtheta=dtheta*theta_ratio + if (theta_edge[-1]+(x2[-1]-theta_edge[-1])*2.0)>(np.pi/2.0): + x2[-1]=(theta_edge[-1]+(np.pi/2.0))/2.0 + alldat["r_edge"]=r_edge alldat["theta_edge"]=theta_edge @@ -296,7 +300,7 @@ def get_ndf_data(fname,r_file,theta_file): #We need to compute the temperature in order to supply it to python. -temp=(3.0/2.0)*data["TOTAL ENERGY"] +temp=(2.0/3.0)*data["TOTAL ENERGY"] data["TEMPERATURE"]=temp/((data["DENSITY"]/(consts.m_p.cgs*MU))*consts.k_B.cgs) # Open an output file diff --git a/py_progs/py_plot_output.py b/py_progs/py_plot_output.py index cab9fba7e..cce4c5a69 100755 --- a/py_progs/py_plot_output.py +++ b/py_progs/py_plot_output.py @@ -38,6 +38,7 @@ import numpy as np import os, sys import py_plot_util as util +import brewer2mpl has_astropy = True try: @@ -45,7 +46,47 @@ except ImportError: has_astropy = False +use_pretty = True +try: + import brewer2mpl +except ImportError: + use_pretty = False + + +colormaps = ["Set1", "Set2", "Dark2", "Paired", "Paired2", "Accent"] + +def set_pretty(cmap="jm"): + + if type(cmap) is str: + if cmap != "jm": + color = cmap + else: + color = "Set1" + elif type(cmap) is int: + color = colormaps[cmap] + else: + raise TypeError("cmap must be string or integer.") + + + + #r.setpars() + + # Get "Set2" colors from ColorBrewer (all colorbrewer scales: http://bl.ocks.org/mbostock/5577023) + set2 = brewer2mpl.get_map(color, 'qualitative', 8).mpl_colors + + if cmap == "jm": + g = set2[2] + b = set2[1] + red = set2[0] + set2[0] = b + set2[1] = g + set2[2] = red + + p.rc("axes", color_cycle=set2) + p.rcParams["legend.frameon"] = "False" + + return set2 def make_spec_plot(s, fname, smooth_factor = 10, angles = True, components = False, with_composite=False): @@ -86,6 +127,7 @@ def make_spec_plot(s, fname, smooth_factor = 10, angles = True, components = Fal ncomponents = 9 + if angles: # first make viewing angle plot @@ -123,7 +165,7 @@ def make_spec_plot(s, fname, smooth_factor = 10, angles = True, components = Fal p.xlabel("Wavelength") p.ylabel("Flux") - p.savefig("spectrum_%s.png" % (fname)) + p.savefig("spectrum_%s.png" % (fname), dpi=300) p.clf() if components: @@ -142,7 +184,7 @@ def make_spec_plot(s, fname, smooth_factor = 10, angles = True, components = Fal p.ylabel("Flux") p.legend() - p.savefig("spec_components_%s.png" % (fname)) + p.savefig("spec_components_%s.png" % (fname), dpi=300) p.clf() return 0 @@ -361,6 +403,8 @@ def make_spec_comparison_plot (s_array, labels, fname="comparison", smooth_facto ncomponents = 9 + if use_pretty: set_pretty() + if angles: @@ -392,7 +436,7 @@ def make_spec_comparison_plot (s_array, labels, fname="comparison", smooth_facto if i == 0 and j == (len(s_array) - 1): p.legend() - p.savefig("spectrum_%s.png" % (fname)) + p.savefig("spectrum_%s.png" % (fname), dpi=300) p.clf() if components: @@ -418,7 +462,7 @@ def make_spec_comparison_plot (s_array, labels, fname="comparison", smooth_facto p.ylabel("Flux") p.legend() - p.savefig("spec_components_%s.png" % (fname)) + p.savefig("spec_components_%s.png" % (fname), dpi=300) p.clf() return 0 @@ -438,6 +482,9 @@ def make_spec_comparison_plot (s_array, labels, fname="comparison", smooth_facto print __doc__ sys.exit(1) + if use_pretty: + set_pretty() + # try to read a parms file in the directory io_print = True try: diff --git a/source/Makefile b/source/Makefile index b7f338519..2a34486a3 100755 --- a/source/Makefile +++ b/source/Makefile @@ -113,7 +113,7 @@ LDFLAGS= -L$(LIB) -lm -lcfitsio -lgsl -lgslcblas #Note that version should be a single string without spaces. -VERSION = 80 +VERSION = 81 CHOICE=1 // Compress plasma as much as possible diff --git a/source/atomic.h b/source/atomic.h index af9e86e64..ffb1c3de3 100755 --- a/source/atomic.h +++ b/source/atomic.h @@ -30,7 +30,7 @@ #define ANGSTROM 1.e-8 /*Definition of an Angstrom in units of this code, e.g. cm */ #define EV2ERGS 1.602192e-12 -#define RADIAN 57.29578 +#define RADIAN 57.29577951308232 #define RYD2ERGS 2.1798741e-11 /* Rydberg in units of ergs */ diff --git a/source/brem.c b/source/brem.c index 51d8d4d81..0d0e832c9 100644 --- a/source/brem.c +++ b/source/brem.c @@ -35,7 +35,7 @@ integ_brem (freq) double freq; { double answer; - answer = geo.const_agn * pow(freq,-0.2) * exp ((-1.0 * H * freq) / (BOLTZMANN * geo.brem_temp)); + answer = geo.const_agn * pow(freq,geo.brem_alpha) * exp ((-1.0 * H * freq) / (BOLTZMANN * geo.brem_temp)); return (answer); } @@ -45,7 +45,7 @@ double double alpha; { double answer; - answer=pow(alpha,-0.2)*exp(-1.0*alpha); + answer=pow(alpha,geo.brem_alpha)*exp(-1.0*alpha); return(answer); } @@ -90,7 +90,7 @@ get_rand_brem (freqmin, freqmax) - /*First time through create the array containing the proper boundaries for the integral of the BB function, + /*First time through create the array containing the proper boundaries for the integral of the brem function, Note calling pdf_gen_from func also defines ylo and yhi */ if (ninit_brem == 0) @@ -102,7 +102,7 @@ get_rand_brem (freqmin, freqmax) Error ("get_rand_brem: on return from pdf_gen_from_func %d\n", echeck); } - /* We need the integral of the bb function outside of the regions of interest as well */ + /* We need the integral of the brem function outside of the regions of interest as well */ cdf_brem_tot = (pow(BREM_ALPHAMIN/100.0,0.8))/0.8+qromb (brem_d, BREM_ALPHAMIN/100.0, BREM_ALPHABIG, 1e-8); cdf_brem_lo = (pow(BREM_ALPHAMIN/100.0,0.8))/0.8+qromb (brem_d, BREM_ALPHAMIN/100.0, BREM_ALPHAMIN, 1e-8) / cdf_brem_tot; //position in the full cdf of low frequcny boundary @@ -191,7 +191,7 @@ reset. A careful review of them is warranted. if (y <= cdf_brem_lo || brem_alphamax < BREM_ALPHAMIN) { - alpha = get_rand_pow (brem_lo_freq_alphamin, brem_lo_freq_alphamax, -0.2); + alpha = get_rand_pow (brem_lo_freq_alphamin, brem_lo_freq_alphamax,geo.brem_alpha); } else if (y >= cdf_brem_hi || brem_alphamin > BREM_ALPHAMAX) { diff --git a/source/diag.c b/source/diag.c index bbb8c92c1..c588f8c0c 100755 --- a/source/diag.c +++ b/source/diag.c @@ -250,10 +250,10 @@ int get_extra_diagnostics() int -save_photon_stats (one, p, ds) +save_photon_stats (one, p, ds, w_ave) WindPtr one; PhotPtr p; - double ds; + double ds,w_ave; { int i; @@ -267,8 +267,8 @@ save_photon_stats (one, p, ds) if (one->nplasma == ncell_stats[i]) { fprintf (pstatptr, - "PHOTON_DETAILS %3d %8.3e %8.3e %8.3e cell%3d wind cell%3d\n", - geo.wcycle, p->freq, p->w, ds, one->nplasma, one->nwind); + "PHOTON_DETAILS cycle %3d n_photon %d freq %8.3e w %8.3e ave_w %8.3e ds %8.3e nscat %d plasma cell %3d wind cell %3d\n", + geo.wcycle, p->np, p->freq, p->w, w_ave,ds, p->nscat, one->nplasma, one->nwind); } } return (0); diff --git a/source/direct_ion.c b/source/direct_ion.c index 38696b7cd..11cdab2b1 100644 --- a/source/direct_ion.c +++ b/source/direct_ion.c @@ -7,8 +7,11 @@ #include "python.h" #include "recipes.h" -#include //We need this gsl library to evaluate the first exponential integral +#include //We need this gsl library to evaluate the first exponential integral +/* Ratio of hnu / kT beyond which we don't bother calculating + see #197 */ +#define ALPHABIG_DIRECT_ION 100. /* direct_ion contains the routines relating to direct (collisional) ionization and threebody recombination */ @@ -43,16 +46,15 @@ compute_di_coeffs (T) { if (ion[n].dere_di_flag == 0) { - //printf ("NO COEFFS\n"); di_coeffs[n] = 0.0; } else - { - di_coeffs[n] = q_ioniz_dere(n, T); - } + { + di_coeffs[n] = q_ioniz_dere(n, T); + } - } //End of loop over ions + } //End of loop over ions return (0); } @@ -96,7 +98,7 @@ compute_qrecomb_coeffs(T) if (ion[n].istate > 1) //There is only any point doing this is we are not a neutral ion { - if (ion[n].dere_di_flag == 0) + if (ion[n-1].dere_di_flag == 0) { //printf ("NO COEFFS\n"); qrecomb_coeffs[n] = 0.0; @@ -129,6 +131,10 @@ compute_qrecomb_coeffs(T) } } //End of if statement for neutral ions + else //We are a neutral ion - so there can be no recombination + { + qrecomb_coeffs[n] = 0.0; + } } //End of loop over ions return (0); @@ -159,41 +165,41 @@ compute_qrecomb_coeffs(T) double total_di (one, t_e) - WindPtr one; // Pointer to the current wind cell - we need the cell volume, this is not in the plasma structure - double t_e; //Current electron temperature of the cell + WindPtr one; // Pointer to the current wind cell - we need the cell volume, this is not in the plasma structure + double t_e; //Current electron temperature of the cell { - double x; //The returned variable - int nplasma; //The cell number in the plasma array + double x; //The returned variable + int nplasma; //The cell number in the plasma array PlasmaPtr xplasma; //pointer to the relevant cell in the plasma structure - int n; //loop pointers + int n; //loop pointers - nplasma = one->nplasma; //Get the correct plasma cell related to this wind cell - xplasma = &plasmamain[nplasma]; //copy the plasma structure for that cell to local variable - x = 0; //zero the luminosity + nplasma = one->nplasma; //Get the correct plasma cell related to this wind cell + xplasma = &plasmamain[nplasma]; //copy the plasma structure for that cell to local variable + x = 0; //zero the luminosity - compute_di_coeffs (t_e); //Calculate the DR coefficients for this cell + compute_di_coeffs (t_e); //Calculate the DR coefficients for this cell for (n = 0; n < nions; n++) { //We have no DI data for this ion - if (ion[n].dere_di_flag == 0) - { - x += 0.0; //Add nothing to the sum of coefficients - } + if (ion[n].dere_di_flag == 0) + { + x += 0.0; //Add nothing to the sum of coefficients + } else - { + { - x += xplasma->vol * xplasma->ne * xplasma->density[n] * di_coeffs[n] * - dere_di_rate[ion[n].nxderedi].xi * EV2ERGS; + x += xplasma->vol * xplasma->ne * xplasma->density[n] * di_coeffs[n] * + dere_di_rate[ion[n].nxderedi].xi * EV2ERGS; //printf ("n=%i V=%e ne=%e rho=%e coeff=%e xi=%e cooling=%e\n",n, V , //xplasma->ne , xplasma->density[n] , di_coeffs[n] , //dere_di_rate[ion[n].nxderedi].xi*EV2ERGS,x); - } + } } return (x); } @@ -228,6 +234,10 @@ q_ioniz_dere (nion, t_e) t = (BOLTZMANN * t_e) / (dere_di_rate[nrec].xi * EV2ERGS); scaled_t = 1.0 - ((log (2.0)) / (log (2.0 + t))); + /* 1/t is (hnu) / (k t_e). when this ratio is >100 we return 0, see #197 */ + if ( (1.0 / t) > ALPHABIG_DIRECT_ION) + return 0.0; + if (scaled_t < dere_di_rate[nrec].temps[0]) //we are below the range of DI data data { @@ -307,6 +317,10 @@ q_ioniz (cont_ptr, electron_temperature) nion = cont_ptr->nion; gaunt = 0.1 * ion[nion].z; //for now - from Mihalas for hydrogen and Helium + /* when hnu / kT ratio is >100 we return 0, see #197 */ + if (u0 > ALPHABIG_DIRECT_ION) + return (0.0); + /* if ion[n].dere_di_flag == 1 then we have direct ionization data for this ion only do this if it is the ground state */ @@ -353,6 +367,11 @@ q_recomb_dere (cont_ptr, electron_temperature) gaunt = 0.1; //for now - from Mihalas for hydrogen u0 = cont_ptr->freq[0] * H_OVER_K / electron_temperature; nion = cont_ptr->nion; + + /* when hnu / kT ratio is >100 we return 0, see #197 */ + if (u0 > ALPHABIG_DIRECT_ION) + return (0.0); + /* if ion[n].dere_di_flag == 1 then we have direct ionization data for this ion only do this if it is the ground state */ if (ion[nion].dere_di_flag == 1 && config[cont_ptr->nlev].ilv == 1) @@ -415,6 +434,9 @@ q_recomb (cont_ptr, electron_temperature) u0 = cont_ptr->freq[0] * H_OVER_K / electron_temperature; gaunt = 0.1 * ion[nion].z; //for now - from Mihalas for hydrogen and Helium + /* when hnu / kT ratio is >100 we return 0, see #197 */ + if (u0 > ALPHABIG_DIRECT_ION) + return (0.0); /* if ion[n].dere_di_flag == 1 then we have direct ionization data for this ion only do this if it is the ground state */ diff --git a/source/emission.c b/source/emission.c index 4890521fa..80bc04c76 100755 --- a/source/emission.c +++ b/source/emission.c @@ -315,14 +315,25 @@ adiabatic_cooling (one, t) double t; { double cooling; - int nplasma; + int nplasma, nion; + double nparticles; PlasmaPtr xplasma; nplasma = one->nplasma; xplasma = &plasmamain[nplasma]; //JM 1401 -- here was an old factor of 3/2 which KSL and JM believe to be incorrect. - cooling = xplasma->ne * BOLTZMANN * t * xplasma->vol * one->div_v; + //JM 1601 -- this should include the pressure from all particles, rather than just ne + //cooling = xplasma->ne * BOLTZMANN * t * xplasma->vol * one->div_v; + nparticles = xplasma->ne; + + for (nion = 0; nion < nions; nion++) + { + /* loop over all ions as they all contribute to the pressure */ + nparticles += xplasma->density[nion]; + } + + cooling = nparticles * BOLTZMANN * t * xplasma->vol * one->div_v; return (cooling); } @@ -692,7 +703,10 @@ total_free (one, t_e, f1, f2) } x = BREMS_CONSTANT * xplasma->ne * (sum) / H_OVER_K; } - + + /* JM 1604 -- The reason why this is proportional to t_e**1/2, + rather than t_e**(-1/2) as in equation 40 of LK02 is because + one gets an extra factor of (k*t_e/h) when one does the integral */ x *= sqrt (t_e) * xplasma->vol; x *= (exp (-H_OVER_K * f1 / t_e) - exp (-H_OVER_K * f2 / t_e)); diff --git a/source/estimators.c b/source/estimators.c index 8706719ac..89b166051 100755 --- a/source/estimators.c +++ b/source/estimators.c @@ -87,7 +87,7 @@ bf_estimators_increment (one, p, ds) */ if (modes.save_cell_stats && ncstat > 0) { - save_photon_stats(one, p, ds); // save photon statistics (extra diagnostics) + save_photon_stats(one, p, ds,p->w); // save photon statistics (extra diagnostics) } diff --git a/source/foo.h b/source/foo.h index d660256d5..896ba6594 100755 --- a/source/foo.h +++ b/source/foo.h @@ -193,6 +193,7 @@ int rtheta_make_hydro_grid(WindPtr w); int rtheta_hydro_volumes(WindPtr w); int hydro_frac(double coord, double coord_array[], int imax, int *cell1, int *cell2, double *frac); double hydro_interp_value(double array[], int im, int ii, int jm, int jj, double f1, double f2); +int hydro_restart(void); /* corona.c */ int get_corona_params(void); double corona_velocity(double x[], double v[]); @@ -257,7 +258,7 @@ double gs_rrate(int nion, double T); /* diag.c */ int open_diagfile(void); int get_extra_diagnostics(void); -int save_photon_stats(WindPtr one, PhotPtr p, double ds); +int save_photon_stats(WindPtr one, PhotPtr p, double ds, double w_ave); /* sv.c */ int get_sv_wind_params(void); double sv_velocity(double x[], double v[]); @@ -559,6 +560,7 @@ int convergence_all(WindPtr w, char rootname[], int ochoice); int model_bands(WindPtr w, char rootname[], int ochoice); int heatcool_summary(WindPtr w, char rootname[], int ochoice); int complete_physical_summary(WindPtr w, char rootname[], int ochoice); +int complete_ion_summary(WindPtr w, char rootname[], int ochoice); double get_density_or_frac(PlasmaPtr xplasma, int element, int istate, int frac_choice); int find_ion(int element, int istate); int find_element(int element); @@ -590,3 +592,4 @@ int main(int argc, char *argv[]); int one_choice(int choice, char *root, int ochoice); int py_wind_help(void); /* test_saha.c */ +int main(int argc, char *argv[]); diff --git a/source/get_atomicdata.c b/source/get_atomicdata.c index 4cfbe18a0..2408476e9 100755 --- a/source/get_atomicdata.c +++ b/source/get_atomicdata.c @@ -179,8 +179,6 @@ get_atomic_data (masterfile) int levl, levu; int in, il; //The levels used in inner shell data double q; - //int nelectrons; - //double et, emax, e0, sigma, ya, p, yw, y0, y1; double freq, f, exx, lambda, alpha, beta, tm, et, p; double the_ground_frac[20]; char choice; @@ -966,6 +964,7 @@ a level type has not been established } else if (ion[n].lev_type != lev_type) { + //XXX Why are these lines commented out? ksl 160212 //OLD Error //OLD ("Get_atomic_data: file %s Reading lev_type (%d) for ion %d with lev_type (%d). Not allowed\n", //OLD file, lev_type, n, ion[n].lev_type); @@ -1116,7 +1115,7 @@ is already incremented { exx *= EV2ERGS; qqnum = ilv = qnum; - lev_type = -2; // It's an old stylle record, one which is only here for backward compatibility + lev_type = -2; // It's an old style record, one which is only here for backward compatibility } else { @@ -2490,7 +2489,8 @@ BAD_T_RR 5 0 1 1 4.647E-10 0.7484 6.142E+01 1.753E+07*/ } //end of loop over ions break; -/* NSH 120921 The following are lines to read in temperature averaged gaunt factors from the data of Sutherland (1997). The atomic file is basically unchaned from the data on the website, just with the top few lines commented out, and a label prepended to each line */ +/* NSH 120921 The following are lines to read in temperature averaged gaunt factors from the data of Sutherland (1997). The atomic file is basically unchanged + * from the data on the website, just with the top few lines commented out, and a label prepended to each line */ case 'g': nparam = sscanf (aline, "%*s %le %le %le %le %le", &gsqrdtemp, &gfftemp, &s1temp, &s2temp, &s3temp); //split and assign the line @@ -2609,7 +2609,7 @@ ion[n].dere_di_flag=1; if (inner_cross[n].n_fluor_yield==-1) /*This is the first yield data for this vacancy */ { inner_fluor_yield[n_fluor_yield_tot].nion=n; /*This yield refers to this ion*/ - inner_cross[n].n_elec_yield=n_fluor_yield_tot; + inner_cross[n].n_fluor_yield=n_fluor_yield_tot; inner_fluor_yield[n_fluor_yield_tot].z=z; inner_fluor_yield[n_fluor_yield_tot].istate=istate; inner_fluor_yield[n_fluor_yield_tot].n=in; diff --git a/source/hydro_import.c b/source/hydro_import.c index c4713d388..e89ffab3c 100755 --- a/source/hydro_import.c +++ b/source/hydro_import.c @@ -508,8 +508,10 @@ hydro_temp (x) temp=hydro_interp_value(temp_input,im,ii,jm,jj,f1,f2); - if (temp < 1e4) //Set a lower limit. - temp = 1e4; + /* NSH 16/2/29 - removed lower limit - set this in the hydro translation software or hydro model */ + +// if (temp < 1e4) //Set a lower limit. +// temp = 1e4; @@ -597,7 +599,7 @@ rtheta_make_hydro_grid (w) hydro_theta_cent[j],theta,thetacen,(theta+2*(thetacen-theta))); } - else if (j >= ihydro_theta) //We are setting up a cell past where there is any data - this is our own ghost cell + else if (j > ihydro_theta) //We are setting up a cell past where there is any data - this is our own ghost cell { Log ("we are past ihydro_theta , %i, %i \n",j,ihydro_theta); thetacen = hydro_thetamax; //Set the center and the edge to the maximum extent of the data/interest @@ -624,7 +626,7 @@ rtheta_make_hydro_grid (w) } } - +/* for (i = 0; i < NDIM; i++) { wind_ij_to_n (i, 0, &n); @@ -636,7 +638,7 @@ for (i = 0; i < MDIM; i++) wind_ij_to_n (0, i, &n); // printf ("hydro_grid j=%i, ihydrotheta=%i, n=%i, theta=%f, thetacen=%f\n",i,ihydro_theta,n,w[n].theta,w[n].thetacen); } - +*/ /* Now set up the wind cones that are needed for calclating ds in a cell */ @@ -822,4 +824,97 @@ double // (1. - f1) * ((1. - f2) * array[im * MAXHYDRO + jm] + f2 * array[im * MAXHYDRO + jj]) + // f1 * ((1. - f2) * array[ii * MAXHYDRO + jm] + f2 * array[ii * MAXHYDRO + jj]); return(value); + } + + + +/*********************************************************** + Southampton + + Synopsis: + hydro_restart is a subrotine which permits a previous wind save + file to be used in a hydro simulation. The density and temperature + for each cell are those from the hydro simulation. The ion fractions + are taken from the windsave file +Arguments: + none + +Returns: + +Description: + + + This sets up a restarted hydro model + + +History: + 16feb nsh 80 -- Coded and debugged. + + +**************************************************************/ + + + + + + int + hydro_restart() + { + int n,nion; + int nwind; + double x[3]; + double old_density; + WindPtr w; + + w = wmain; + geo.wind_type=3; //Temporarily set the wind type to hydro, so we can use the normal routines + NDIM = ndim = geo.ndim; + MDIM = mdim = geo.mdim; + NDIM2 = NDIM * MDIM; + for (n = 0; n < NDIM2; n++) + { + /* 04aug -- ksl -52 -- The next couple of lines are part of the changes + * made in the program to allow more that one coordinate system in python + */ + + model_velocity (w[n].x, w[n].v); + model_vgrad (w[n].x, w[n].v_grad); + } + for (n = 0; n < NPLASMA; n++) + { + nwind = plasmamain[n].nwind; + stuff_v (w[nwind].xcen, x); + old_density=plasmamain[n].rho; + plasmamain[n].rho = model_rho (x)/geo.fill; + plasmamain[n].t_r = plasmamain[n].t_e =hydro_temp (x); + for (nion=0;nionnrad); Log ("xyz %8.2e %8.2e %8.2e vel %8.2e %8.2e %8.2e\n", w[n].x[0], w[n].x[1], w[n].x[2], w[n].v[0], w[n].v[1], w[n].v[2]); + Log ("r theta %12.6e %12.6e \n", w[n].rcen, w[n].thetacen/RADIAN); + Log ("nh %8.2e ne %8.2e t_r %8.2e t_e %8.2e w %8.2e vol %8.2e\n", xplasma->rho * rho2nh, xplasma->ne, xplasma->t_r, xplasma->t_e, xplasma->w, w[n].vol); @@ -1515,9 +1517,9 @@ rdint ("Wind.array.element", &n); ("t_e %8.2e cool_tot %8.2e lum_lines %8.2e lum_ff %8.2e lum_fb %8.2e cool_comp %8.2e cool_adiab %8.2e cool_DR %8.2e \n", xplasma->t_e, xplasma->lum_rad_ioniz+xplasma->lum_comp_ioniz+xplasma->lum_adiabatic_ioniz+xplasma->lum_dr_ioniz, xplasma->lum_lines_ioniz, xplasma->lum_ff_ioniz, xplasma->lum_fb_ioniz, xplasma->lum_comp_ioniz, xplasma->lum_adiabatic_ioniz, xplasma->lum_dr_ioniz); Log - ("t_r %8.2e heat_tot %8.2e heat_lines %8.2e heat_ff %8.2e heat_photo %8.2e heat_comp %8.2e heat_icomp %8.2e\n", + ("t_r %8.2e heat_tot %8.2e heat_lines %8.2e heat_ff %8.2e heat_photo %8.2e heat_auger %8.2e heat_comp %8.2e heat_icomp %8.2e\n", xplasma->t_r, xplasma->heat_tot, xplasma->heat_lines, xplasma->heat_ff, - xplasma->heat_photo, xplasma->heat_comp,xplasma->heat_ind_comp); + xplasma->heat_photo, xplasma->heat_auger,xplasma->heat_comp,xplasma->heat_ind_comp); @@ -3235,14 +3237,14 @@ complete_physical_summary (w, rootname, ochoice) printf("n\tnplasma\tinwind\ti\tj\tx\tz\tv\tvx\tvy\tvz\tdvds_ave\tvol\t \ rho\tne\tte\ttr\tnphot\tw\tave_freq\tIP\tconv\tconv_tr\tconv_te\tconv_hc\t \ lum_tot\tlum_rad\tlum_fb\tlum_ff\tlum_lines\tlum_adiabatic\tlum_comp\tlum_dr\t \ -heat_tot\theat_photo\theat_lines\theat_ff\theat_comp\theat_ind_comp\t \ +heat_tot\theat_photo\theat_auger\theat_lines\theat_ff\theat_comp\theat_ind_comp\t \ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\n"); if (ochoice) - fprintf(fptr, "n\tnplasma\tinwind\ti\tj\tx\tz\tv\tvx\tvy\tvz\tdvds_ave\tvol\t \ -rho\tne\tte\ttr\tnphot\tw\tave_freq\tIP\tconv\tconv_tr\tconv_te\tconv_hc\t \ + fprintf(fptr, "n\tnplasma\tinwind\ti\tj\tx\tz\tr\ttheta\tv\tvx\tvy\tvz\tdvds_ave\tvol\t \ +rho\tne\tte\ttr\tnphot\tw\tave_freq\tIP\tXi\tconv\tconv_tr\tconv_te\tconv_hc\t \ lum_tot\tlum_rad\tlum_fb\tlum_ff\tlum_lines\tlum_adiabatic\tlum_comp\tlum_dr\t \ -heat_tot\theat_photo\theat_lines\theat_ff\theat_comp\theat_ind_comp\t \ +heat_tot\theat_photo\theat_auger\theat_lines\theat_ff\theat_comp\theat_ind_comp\t \ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\n"); @@ -3272,7 +3274,7 @@ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\ o6den = get_density_or_frac(xplasma,8,6, frac_choice); si4den = get_density_or_frac(xplasma,14,4, frac_choice); - /* printf("%i %i %i %i %i %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ + /* printf("%i %i %i %i %i %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %8.4e \ %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ @@ -3289,22 +3291,23 @@ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\ */ if (ochoice) - fprintf(fptr, "%i %i %i %i %i %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ - %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %8.4e \ - %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ + fprintf(fptr, "%i %i %i %i %i %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ + %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %8.4e %i %8.4e %8.4e %8.4e %8.4e \ + %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e\ %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e \ %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e %8.4e\n", - n, np, w[n].inwind, ii, jj, w[n].x[0], w[n].x[2], vtot, w[n].v[0], w[n].v[1], w[n].v[2], w[n].dvds_ave, w[n].vol, + n, np, w[n].inwind, ii, jj, w[n].x[0], w[n].x[2],w[n].rcen,w[n].thetacen/RADIAN, vtot, w[n].v[0], w[n].v[1], w[n].v[2], w[n].dvds_ave, w[n].vol, plasmamain[np].rho, plasmamain[np].ne, plasmamain[np].t_e, plasmamain[np].t_r, plasmamain[np].ntot, - plasmamain[np].w, plasmamain[np].ave_freq, plasmamain[np].ip, plasmamain[np].converge_whole, + plasmamain[np].w, plasmamain[np].ave_freq, plasmamain[np].ip, plasmamain[np].xi, plasmamain[np].converge_whole, plasmamain[np].converge_t_r, plasmamain[np].converge_t_e, plasmamain[np].converge_hc, - plasmamain[np].lum_ioniz, plasmamain[np].lum_rad, plasmamain[np].lum_fb, + plasmamain[np].lum_rad+plasmamain[np].lum_comp+plasmamain[np].lum_adiabatic+plasmamain[np].lum_dr, + plasmamain[np].lum_rad, plasmamain[np].lum_fb, plasmamain[np].lum_ff, plasmamain[np].lum_lines, plasmamain[np].lum_adiabatic, - plasmamain[np].lum_comp, plasmamain[np].lum_dr, plasmamain[np].heat_tot, plasmamain[np].heat_photo, + plasmamain[np].lum_comp, plasmamain[np].lum_dr, plasmamain[np].heat_tot, plasmamain[np].heat_photo, plasmamain[np].heat_auger, plasmamain[np].heat_lines , plasmamain[np].heat_ff , plasmamain[np].heat_comp, plasmamain[np].heat_ind_comp, h1den, h2den, he1den, he2den, he3den, c3den, c4den, c5den, n5den, o6den, si4den); } - else + else { /* if we aren't inwind then print out a load of zeroes */ @@ -3315,30 +3318,109 @@ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", n, np, w[n].inwind, ii, jj, w[n].x[0], w[n].x[2]); */ - + if (ochoice) - fprintf(fptr, "%i %i %i %i %i %8.4e %8.4e 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ + fprintf(fptr, "%i %i %i %i %i %8.4e %8.4e 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ - 0.0 0.0 0.0 0.0 0.0 0.0 \ + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", - n, np, w[n].inwind, ii, jj, w[n].x[0], w[n].x[2]); + n, np, w[n].inwind, ii, jj, w[n].x[0], w[n].x[2]); } } if (ochoice) { - fclose (fptr); printf("\nSaved summary of physical quantites in %s, use py_read_output.py to read\n", filename); } return (0); + + +} + +/************************************************************************** + + + Synopsis: + complete_ion_summary outputs a file with all of the ion fractions for a given cell + + + History: + 1602 NSH coded + +************************************************************************/ + + + + int + complete_ion_summary (w, rootname, ochoice) + WindPtr w; + char rootname[]; + int ochoice; + + { + char cell[5]; + PlasmaPtr xplasma; + FILE *fptr, *fopen (); + char filename[LINELENGTH]; + + + int n, mm; + n = 50; + a: printf("There are %i wind elements in this model\n",NDIM2); + rdint ("Wind.array.element", &n); + + if (n < 0) + goto b; + else if (n > NDIM2) + { + printf("No, there are %i wind elements, not %i\n",NDIM2,n); + goto a; + } + printf ("OK, using cell %i\n",n); + xplasma = &plasmamain[w[n].nplasma]; + + if (ochoice) + { + strcpy (filename, rootname); + strcat (filename,"_cell_"); + sprintf(cell,"%04d",n); + strcat (filename,cell); + printf ("opening file %s\n",filename); + + strcat (filename, ".ions"); + printf ("opening file %s\n",filename); + fptr = fopen (filename, "w"); + } + + printf ("ion z n(ion) n(ion)/n(H)\n"); + + + if (ochoice) + fprintf(fptr, "ion z n(ion) n(ion)/n(H)\n"); + + + for (mm = 0; mm < nions; mm++) + { + + printf ("%i %i %e %e\n",mm,ion[mm].z,xplasma->density[mm],xplasma->density[mm]/(xplasma->rho*rho2nh)); + if (ochoice) + { + fprintf (fptr,"%i %i %e %e\n",mm,ion[mm].z,xplasma->density[mm],xplasma->density[mm]/(xplasma->rho*rho2nh)); + } + + } + goto a; + + b:return (0); } + /************************************************************************** diff --git a/source/python.c b/source/python.c index 965f2ca6b..3e57de8fd 100755 --- a/source/python.c +++ b/source/python.c @@ -38,6 +38,7 @@ not attempt to seek a new temperature, but it does output heating and cooling rates --version print out python version, commit hash and if there were files with uncommitted changes + --seed set the random number seed to be time based, rather than fixed. if one simply types py or pyZZ where ZZ is the version number one is queried for a name @@ -217,6 +218,7 @@ #include #include #include "atomic.h" +#include //To allow the used of the clock command without errors!! #include "python.h" @@ -503,6 +505,15 @@ main (argc, argv) ("There is a problem in allocating memory for the photon structure\n"); exit (0); } + else + { + /* JM 1605 -- large photon numbers can cause problems / runs to crash. Report to use (see #209) */ + Log("Allocated %10d bytes for each of %5d elements of photon structure totaling %10.1f Mb \n", + sizeof (p_dummy), NPHOT, 1.e-6 * NPHOT * sizeof (p_dummy)); + if ( (NPHOT * sizeof (p_dummy)) > 1e9) + Error("Over 1 GIGABYTE of photon structure allocated. Could cause serious problems.\n"); + } + /* Define the coordinate system for the grid and allocate memory for the wind structure by reading from user */ @@ -674,6 +685,11 @@ main (argc, argv) } } } + if (modes.zeus_connect==1) /* We are in rad-hydro mode, we want the new density and temperature*/ + { + Log ("We are going to read in the density and temperature from a zeus file\n"); + get_hydro (); //This line just populates the hydro structures + } } /* 121219 NSH Set up DFUDGE to be a value that makes some kind of sense @@ -998,6 +1014,11 @@ main (argc, argv) } // Do not reinit if you want to use old windfile +else if (modes.zeus_connect==1) //We have restarted, but are in zeus connect mode, so we want to update density, temp and velocities +{ + hydro_restart(); +} + w = wmain; if (modes.save_cell_stats) @@ -1007,8 +1028,17 @@ main (argc, argv) } /* initialize the random number generator */ - // srand( (n=(unsigned int) clock())); - srand (1084515760+(13*rank_global)); + /* JM 1503 -- Sometimes it is useful to vary the random number seed. + Default is fixed, but will vary with different processor numbers */ + /* We don't want to run the same photons each cycle in zeus mode, so + everytime we are using zeus we also set to use the clock */ + if ( (modes.rand_seed_usetime == 1) || (modes.zeus_connect == 1) ) + { + n=(unsigned int) clock()*(rank_global+1); + srand(n); + } + else + srand (1084515760+(13*rank_global)); /* 68b - 0902 - ksl - Start with photon history off */ @@ -1253,7 +1283,17 @@ main (argc, argv) /* Calculate and store the amount of heating of the disk due to radiation impinging on the disk */ + /* We only want one process to write to the file */ +#ifdef MPI_ON + if (rank_global == 0) + { +#endif qdisk_save (files.disk, ztot); + +#ifdef MPI_ON + } + MPI_Barrier(MPI_COMM_WORLD); +#endif /* Completed writing file describing disk heating */ @@ -1265,25 +1305,6 @@ main (argc, argv) wind_update (w); -/* In a diagnostic mode save the wind file for each cycle (from thread 0) */ - - if (modes.keep_ioncycle_windsaves) - { - strcpy (dummy, ""); - sprintf (dummy, "python%02d.wind_save", geo.wcycle); - -#ifdef MPI_ON - if (rank_global == 0) - { -#endif - wind_save (dummy); -#ifdef MPI_ON - } -#endif - Log ("Saved wind structure in %s\n", dummy); - } - - Log ("Completed ionization cycle %d : The elapsed TIME was %f\n", geo.wcycle, timer ()); @@ -1305,13 +1326,13 @@ main (argc, argv) #endif spectrum_summary (files.wspec, "w", 0, 6, 0, 1., 0); spectrum_summary (files.lspec, "w", 0, 6, 0, 1., 1); /* output the log spectrum */ - + phot_gen_sum (files.phot, "w"); /* Save info about the way photons are created and absorbed + by the disk */ #ifdef MPI_ON } MPI_Barrier(MPI_COMM_WORLD); #endif - phot_gen_sum (files.phot, "w"); /* Save info about the way photons are created and absorbed - by the disk */ + /* Save everything after each cycle and prepare for the next cycle JM1304: moved geo.wcycle++ after xsignal to record cycles correctly. First cycle is cycle 0. */ @@ -1332,7 +1353,18 @@ main (argc, argv) #endif wind_save (files.windsave); Log_silent ("Saved wind structure in %s after cycle %d\n", files.windsave, - geo.wcycle); + geo.wcycle); + /* In a diagnostic mode save the wind file for each cycle (from thread 0) */ + + if (modes.keep_ioncycle_windsaves) + { + strcpy (dummy, ""); + sprintf (dummy, "python%02d.wind_save", geo.wcycle); + wind_save (dummy); + Log ("Saved wind structure in %s\n", dummy); + } + + #ifdef MPI_ON } MPI_Barrier(MPI_COMM_WORLD); @@ -1537,9 +1569,14 @@ main (argc, argv) /* XXXX -- END CYCLE TO CALCULATE DETAILED SPECTRUM */ - +#ifdef MPI_ON + if (rank_global == 0) + { +#endif phot_gen_sum (files.phot, "a"); - +#ifdef MPI_ON + } +#endif /* 57h - 07jul -- ksl -- Write out the freebound information */ #ifdef MPI_ON @@ -2135,6 +2172,7 @@ int init_advanced_modes() modes.quit_after_inputs = 0; // testing mode which quits after reading in inputs modes.fixed_temp = 0; // do not attempt to change temperature - used for testing modes.zeus_connect = 0; // connect with zeus + modes.rand_seed_usetime = 0; // default random number seed is fixed, not based on time //note this is defined in atomic.h, rather than the modes structure diff --git a/source/python.h b/source/python.h index cea316a17..927ac1c83 100755 --- a/source/python.h +++ b/source/python.h @@ -343,6 +343,8 @@ struct geometry // The next set of parameters relate to the central source of an AGN double brem_temp; /*The temperature of a bremsstrahlung source */ + double brem_alpha; /*The exponent of the nu term for a bremstrahlung source */ + double pl_low_cutoff; /* accessible only in advanced mode- see #34. default to zero */ double alpha_agn; /*The power law index of a BH at the center of an AGN. Note that the luminosity @@ -1174,6 +1176,7 @@ struct advanced_modes int quit_after_inputs; // quit after inputs read in, testing mode int fixed_temp; // do not alter temperature from that set in the parameter file int zeus_connect; // We are connecting to zeus, do not seek new temp and output a heating and cooling file + int rand_seed_usetime; // default random number seed is fixed, not based on time } modes; diff --git a/source/radiation.c b/source/radiation.c index c8859fad8..c9661e15d 100755 --- a/source/radiation.c +++ b/source/radiation.c @@ -462,7 +462,7 @@ if (freq > phot_freq_min) */ if (modes.save_cell_stats && ncstat > 0) { - save_photon_stats (one, p, ds); // save photon statistics (extra diagnostics) + save_photon_stats (one, p, ds, w_ave); // save photon statistics (extra diagnostics) } diff --git a/source/recomb.c b/source/recomb.c index d4ac29441..7d39a86b7 100644 --- a/source/recomb.c +++ b/source/recomb.c @@ -165,10 +165,23 @@ fb_topbase_partial (freq) return (0.0); // No recombination at frequencies lower than the threshold freq occur nion = fb_xtop->nion; - gn = config[fb_xtop->nlev].g; + + /* JM -- below lines to address bug #195 */ + if (ion[nion].phot_info > 0) // it's a topbase record + gn = config[fb_xtop->nlev].g; + else if (ion[nion].phot_info == 0) // it's a VFKY record, so shouldn't really use levels + gn = ion[nion].g; + else + { + Error("fb_topbase_partial: Did not understand cross-section type %i for ion %i. Setting multiplicity to zero!\n", + ion[nion].phot_info, nion); + gn = 0.0; + } + + gion = ion[nion + 1].g; // Want the g factor of the next ion up x = sigma_phot (fb_xtop, freq); -// Now calculate emission using Ferland's expression + // Now calculate emission using Ferland's expression partial = @@ -178,7 +191,7 @@ fb_topbase_partial (freq) -// 0=emissivity, 1=heat loss from electrons, 2=photons emissivity + // 0=emissivity, 1=heat loss from electrons, 2=photons emissivity if (fbfr == 1) partial *= (freq - fthresh) / freq; else if (fbfr == 2) @@ -1345,7 +1358,7 @@ total_rrate (nion, T) else if (total_rr[ion[nion].nxtotalrr].type == RRTYPE_SHULL) { rate = - total_rr[ion[nion].nxtotalrr].params[0] * pow ((T / 1.0e4), + total_rr[ion[nion].nxtotalrr].params[0] * pow ((T / 1.0e4),-1.0* total_rr[ion [nion]. nxtotalrr]. diff --git a/source/setup.c b/source/setup.c index 0c7b37297..5c2bfe816 100644 --- a/source/setup.c +++ b/source/setup.c @@ -104,6 +104,12 @@ int parse_command_line(argc, argv) { modes.fixed_temp = 1; } + + /* JM 1503 -- Sometimes it is useful to vary the random number seed. Set a mode for that */ + else if (strcmp (argv[i], "--rseed") == 0) + { + modes.rand_seed_usetime = 1; + } else if (strcmp (argv[i], "-z") == 0) { modes.zeus_connect = 1; @@ -481,7 +487,6 @@ int get_radiation_sources() "Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table,5=bremsstrahlung)_to_make_wind", &geo.agn_ion_spectype); - printf ("NSH %i\n",geo.agn_ion_spectype); /* 130621 - ksl - This is a kluge to add a power law to stellar systems. What id done is to remove the bl emission, which we always assume to some kind of temperature driven source, and replace it with a power law source @@ -845,8 +850,10 @@ int get_bl_and_agn_params (lstar) else if (geo.agn_ion_spectype == SPECTYPE_BREM) { geo.brem_temp=1.16e8; //10kev + geo.brem_alpha=-0.2; //This is the cloudy form of bremstrahlung geo.const_agn=1.0; - rddoub ("agn_bremsstrahung_temp(K)",&geo.brem_temp); + rddoub ("agn_bremsstrahlung_temp(K)",&geo.brem_temp); + rddoub ("agn_bremsstrahlung_alpha",&geo.brem_alpha); temp_const_agn = geo.lum_agn / qromb(integ_brem,4.84e17,2.42e18,1e-4); geo.const_agn=temp_const_agn; Log ("AGN Input parameters give a Bremsstrahlung constant of %e\n", temp_const_agn); diff --git a/source/templates.h b/source/templates.h index d660256d5..896ba6594 100755 --- a/source/templates.h +++ b/source/templates.h @@ -193,6 +193,7 @@ int rtheta_make_hydro_grid(WindPtr w); int rtheta_hydro_volumes(WindPtr w); int hydro_frac(double coord, double coord_array[], int imax, int *cell1, int *cell2, double *frac); double hydro_interp_value(double array[], int im, int ii, int jm, int jj, double f1, double f2); +int hydro_restart(void); /* corona.c */ int get_corona_params(void); double corona_velocity(double x[], double v[]); @@ -257,7 +258,7 @@ double gs_rrate(int nion, double T); /* diag.c */ int open_diagfile(void); int get_extra_diagnostics(void); -int save_photon_stats(WindPtr one, PhotPtr p, double ds); +int save_photon_stats(WindPtr one, PhotPtr p, double ds, double w_ave); /* sv.c */ int get_sv_wind_params(void); double sv_velocity(double x[], double v[]); @@ -559,6 +560,7 @@ int convergence_all(WindPtr w, char rootname[], int ochoice); int model_bands(WindPtr w, char rootname[], int ochoice); int heatcool_summary(WindPtr w, char rootname[], int ochoice); int complete_physical_summary(WindPtr w, char rootname[], int ochoice); +int complete_ion_summary(WindPtr w, char rootname[], int ochoice); double get_density_or_frac(PlasmaPtr xplasma, int element, int istate, int frac_choice); int find_ion(int element, int istate); int find_element(int element); @@ -590,3 +592,4 @@ int main(int argc, char *argv[]); int one_choice(int choice, char *root, int ochoice); int py_wind_help(void); /* test_saha.c */ +int main(int argc, char *argv[]); diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c index 156986dd1..714a85630 100644 --- a/source/wind_updates2d.c +++ b/source/wind_updates2d.c @@ -665,7 +665,7 @@ free (commbuffer); { Log("Outputting heatcool file for connecting to zeus\n"); fptr = fopen ("py_heatcool.dat", "w"); - fprintf(fptr,"i j rcen thetacen vol temp xi ne heat_xray heat_comp heat_lines heat_ff cool_comp cool_lines cool_ff\n"); + fprintf(fptr,"i j rcen thetacen vol temp xi ne heat_xray heat_comp heat_lines heat_ff cool_comp cool_lines cool_ff rho n_h\n"); } @@ -729,24 +729,10 @@ free (commbuffer); - if (modes.zeus_connect==1) //If we are running in zeus connect mode, we output heating and cooling rates. - { - wind_n_to_ij (plasmamain[nplasma].nwind, &i, &j); - vol=w[plasmamain[nplasma].nwind].vol; - fprintf(fptr,"%d %d %e %e %e ",i,j,w[plasmamain[nplasma].nwind].rcen,w[plasmamain[nplasma].nwind].thetacen/RADIAN,vol); //output geometric things - fprintf(fptr,"%e %e %e ",plasmamain[nplasma].t_e,plasmamain[nplasma].xi,plasmamain[nplasma].ne); //output temp, xi and ne to ease plotting of heating rates - fprintf(fptr,"%e ",(plasmamain[nplasma].heat_photo+plasmamain[nplasma].heat_auger)/vol); //Xray heating - or photoionization - fprintf(fptr,"%e ",(plasmamain[nplasma].heat_comp)/vol); //Compton heating - fprintf(fptr,"%e ",(plasmamain[nplasma].heat_lines)/vol); //Line heating 28/10/15 - not currently used in zeus - fprintf(fptr,"%e ",(plasmamain[nplasma].heat_ff)/vol); //FF heating 28/10/15 - not currently used in zeus - fprintf(fptr,"%e ",(plasmamain[nplasma].lum_comp)/vol); //Compton cooling - fprintf(fptr,"%e ",(plasmamain[nplasma].lum_lines+plasmamain[nplasma].lum_fb+plasmamain[nplasma].lum_dr)/vol); //Line cooling must include all recombinatiobs cooling - fprintf(fptr,"%e\n",(plasmamain[nplasma].lum_ff)/vol); //ff cooling - } + } - if (modes.zeus_connect==1) - fclose(fptr); + /* JM130621- bugfix for windsave bug- needed so that we have the luminosities from ionization cycles in the windsavefile even if the spectral cycles are run */ @@ -782,6 +768,28 @@ free (commbuffer); asum = wind_luminosity (0.0, VERY_BIG); /*We call wind_luminosity here to obtain an up to date set of cooling rates*/ + if (modes.zeus_connect==1) //If we are running in zeus connect mode, we output heating and cooling rates. + { + for (nplasma = 0; nplasma < NPLASMA; nplasma++) + { + + wind_n_to_ij (plasmamain[nplasma].nwind, &i, &j); + vol=w[plasmamain[nplasma].nwind].vol; + fprintf(fptr,"%d %d %e %e %e ",i,j,w[plasmamain[nplasma].nwind].rcen,w[plasmamain[nplasma].nwind].thetacen/RADIAN,vol); //output geometric things + fprintf(fptr,"%e %e %e ",plasmamain[nplasma].t_e,plasmamain[nplasma].xi,plasmamain[nplasma].ne); //output temp, xi and ne to ease plotting of heating rates + fprintf(fptr,"%e ",(plasmamain[nplasma].heat_photo+plasmamain[nplasma].heat_auger)/vol); //Xray heating - or photoionization + fprintf(fptr,"%e ",(plasmamain[nplasma].heat_comp)/vol); //Compton heating + fprintf(fptr,"%e ",(plasmamain[nplasma].heat_lines)/vol); //Line heating 28/10/15 - not currently used in zeus + fprintf(fptr,"%e ",(plasmamain[nplasma].heat_ff)/vol); //FF heating 28/10/15 - not currently used in zeus + fprintf(fptr,"%e ",(plasmamain[nplasma].lum_comp)/vol); //Compton cooling + fprintf(fptr,"%e ",(plasmamain[nplasma].lum_lines+plasmamain[nplasma].lum_fb+plasmamain[nplasma].lum_dr)/vol); //Line cooling must include all recombinatiobs cooling + fprintf(fptr,"%e ",(plasmamain[nplasma].lum_ff)/vol); //ff cooling + fprintf(fptr,"%e ",plasmamain[nplasma].rho); //density + fprintf(fptr,"%e\n",plasmamain[nplasma].rho*rho2nh); //hydrogen number density + } + fclose(fptr); +} + /* 1108 NSH Added commands to report compton heating */ Log ("!!wind_update: Absorbed flux %8.2e (photo %8.2e ff %8.2e compton %8.2e auger %8.2e induced_compton %8.2e lines %8.2e)\n", xsum, psum, fsum, csum, ausum, icsum, lsum); diff --git a/source/windsave2table.c b/source/windsave2table.c index a7d017a9e..79ea12707 100644 --- a/source/windsave2table.c +++ b/source/windsave2table.c @@ -227,7 +227,10 @@ create_master_table (rootname) c[10]=get_one("dmo_dt_z"); strcpy (column_name[10], "dmo_dt_z"); - ncols = 11; + c[11] = get_one ("ntot"); + strcpy (column_name[11], "ntot"); + + ncols = 12; converge=get_one("converge"); @@ -670,6 +673,11 @@ get_one (variable_name) { x[n]=plasmamain[nplasma].dmo_dt[2]; } + else if (strcmp (variable_name, "ntot")==0) + { + x[n]=plasmamain[nplasma].ntot; + } +