% file: vh1xray.sl % % Routines etc to read, plot, and generate X-ray emission % from the output of a VH-1 hydrodynamic simulation, assumed % to be 1-D and spherical. % % See also: http://space.mit.edu/home/dd/VH1Xray/ % % - - - % OTHER FILES % These files provide useful functions; however, the code is written so % that they are not required (comment out the evalfile line(s) below.) % % mxs_abund_pkg.sl -- This is needed if non-solar (1.0) abundances % are to be used for ejecta, CSM. % The file is available on the web page above. % % tetp_equilib_bmin.sl -- Function to calculate beta (=Te/Tp) as a % function of v_shock and tau: % tetp_equilib_bmin(vsh,tau,bmin[or NULL]); % Includes a beta-min of 0.06, Adelsberg+ 2008 % % oplot_common_HHe.sl -- Simple overplotting of H,He-like line % locations on spectral plot. % % - - - % USAGE % % Start ISIS from linux prompt % [linux] isis % % Load this file to define variables and routines - done once a session. % isis> .load vh1xray % % Set values for user-provided variables either by editing % them in this file in the section below, or % by making a file, e.g., mysim_values.sl: % see/copy the section below as an example of such a file. % isis> .source mysim_values % skip if the values in vh1xray are correct. % % Setup the history file, etc. % isis> vhx_setup; % % This routine takes the abundance values and calculates useful % plasma properties based on them: mu, Nion/Namu, Ne/Namu, NH/Ne. % It is called as needed by other routines but is available to % the user as well, e.g., to see effect of abundances being changed. % isis> vhx_calc_muetc; % isis> print(vhx_cmpmu); % also: 'cmpninamu, 'cmpnenamu, 'cmpnhne % % Read in a time step and plot it % isis> vhx_read_step(250); % isis> vhx_plot_step; % % Determine the location of the CD at a reference time step. % The step is assigned to vhx_cd0step and the CD radius * in pc * % is put in vhx_cd0rad. These "given-CD" values can be put in the % mysim_values file once determined, e.g.: % isis> vhx_cd0step = 165; % isis> vhx_cd0rad = 2.26989; % % Alternately, the user can point and click to measure the CD radius % (at a given time step) by doing: % isis> vhx_cd0step=NULL; % or a desired step number. % isis> vhx_cursor_cd0; % When the user clicks on the CD, the suggested values will be % typed to the screen and can then be put in the mysim_values % file, and/or they can be copy/pasted to the command line, e.g.: % isis> vhx_cd0step = 165; % isis> vhx_cd0rad = 2.26989; % % Provided the user values of 'cd0step and 'cd0rad have been % accurately set, evaluate the CD (and RS, FS) locations for % all time steps. % isis> vhx_calc_cds; % ... takes some time ... % After this is run, the 'plot_steps ouptut will show the CD too. % % Plot the RS, CD, FS vs time % isis> vhx_plot_cds; % or use vhx_psplot_cds; to make .ps plot % and the amount of shocked mass vs time: % isis> vhx_plot_mass; % or use vhx_psplot_mass; to make .ps plot % % Make gif's from the steps ('beg to 'end, every 'nth) % and save them 'outdir. When done it also shows the % user the convert command to make a movie from the images. % isis> vhx_plot_allsteps; % % Setup the mass-cuts 'automatically': % isis> vhx_ncutsej = 15; % defaults can be in mysim_values % isis> vhx_ncutscsm = 20; % isis> vhx_make_cuts; % % If some of the mass-cut shells will be "clumped", then put the % logic, etc here. (Otherwise the default is no clumping via NULL % or by 'make_cuts which sets all of shl_clumped to zeros) % For exmaple, clump all CSM (component = 2) shells: % isis> shl_clumped = 0*shl_icmp; % define and reset clumped array % isis> shl_clumped[ where(shl_icmp == 2) ] = 1; % % Calculate and ouput information (data files) about the mass-cut shells % for all desired time steps. % isis> vhx_calc_shells; % writes shells_1NNN.txt's % % Read in the shells file for a step and show the kT-tau values % in it that will generate X-rays: % isis> vhx_plot_kTtau(tid); or use vhx_psplot_kTtau(tid); to make .ps plot % % Generate an isis.par spectral model file for a time step -- % uses/requires the shells_1NNN.txt data file as input. % This routine may be useful to make a single par file, % and/or it also fills the shl_ variables, but it is generally not % used directly, instead it is called when vhx_calc_allspec is run. % isis> vhx_make_par(tid); % writes model_1.par"); % % With the shells'.txt and model'.par files made for a step, % can calculate the fluxes from each mass-cut shell, % leaving them in shl_fluxes[iEnergyBand, ishl]: % isis> vhx_calc_shellfluxes(tid); % % Make and evaluate all the desired steps' par files, % write spectral plots for each, and write a fluxlocs-vs-time data file. % (Obs or NoNH suffix depending in vhx_nh22 value.) % isis> vhx_calc_allspec; % writes spectrum_1NNN.gif's % and fluxlocs_vs_time_[Obs,NoNH].txt"); % % Plot the fluxes-in-bands vs time (from values in the fluxlocs' file) % isis> vhx_plot_lc;"); % or use vhx_psplot_lc; to make .ps plot % - - - % 12/20-30/12 - started by dd@space.mit.edu; tracking the CD. % 1/ 8-10/13 - improved plotting options, added 'plot_mass. % --2/4/13 - add mass-cuts, abundances; calc fluxes. % --2/25/13 - mostly small changes. Adjusted mass-cut making. % --2/27/13 - lc plotting changes, added options. % -- 4/9/13 - re-number models in par file to be same as shell number; % kTtau plot-symbol size depends on norm. Add shl_fluxes and % 'calc_shellfluxes(tid); output .ps and .txt flux files. % Include shl_rhodr values in shells'txt file, in amu/cm^2. % - - - % - - - % This is needed for the abundances to be set to other than % 1.0 times the "angr" or "wilm" solar values. ()=evalfile("mxs_abund_pkg"); % Routine to estimate Te/Tp = beta = function(v_shock, Tau); % If not available then beta = vhx_betamin . ()=evalfile("tetp_equilib_bmin"); % To overplot H,He-like line locations on spectral plot: ()=evalfile("oplot_common_HHe"); % - - - % User-set variables to control the input, operations, etc. % See following section for some description and example % values. public variable vhx_abovedir, vhx_simdir, vhx_simpre, vhx_simhst; public variable vhx_simbeg, vhx_simend, vhx_simnth, vhx_outdir; public variable vhx_dkpc, vhx_systvel, vhx_nh22, vhx_frac4pi; public variable vhx_cd0step, vhx_cd0rad; public variable vhx_Tshock; public variable shl_mass, shl_icmp, shl_clumped, vhx_clmpfill, vhx_clmpdens; public variable vhx_ncutsej, vhx_ncutscsm, vhx_ncutsnear; public variable vhx_abundstr, vhx_abundzs, vhx_abundej, vhx_abundcsm; public variable vhx_xneiverstr, vhx_xxsectstr, vhx_betamin; public variable vhx_massrelcd, vhx_refstep; public variable vhx_pltxmin, vhx_pltxmax, vhx_pltymin, vhx_pltymax; public variable vhx_fluxmin, vhx_fluxmax, vhx_fluxelos, vhx_fluxehis; public variable vhx_lc_xlog, vhx_lc_showejcsm; % - - - - - - - % - - - - - - - % Example values for the user-set public/global variables. % Can either edit the values here or, preferably, % copy this section to a file, e.g., mysim_values.sl, % and .source it to set the values after vh1xray is loaded % (see USAGE above and output of vhx_help;). % - - - Input and output files % Point to the VH-1 simulation dir, the dirs above it, % the output filenames prefix, the history file vhx_abovedir = "/nfs/spectra/d1/Science/Vikram/Tycho_2013"; vhx_simdir = "tycho6_110216"; vhx_simpre = "tycho"; vhx_simhst = "tychohst"; % Output products go into: vhx_outdir = "./Tycho_out"; % - - - Set start and ending time steps % 1000 will be added to these to generate the step's filename vhx_simbeg = NULL; vhx_simend = NULL; % the RS moves out of range later % output only every nth simulation step (can be NULL to be s/w set) vhx_simnth = NULL; % - - - Information about the source % Distance, systemic velocity, and foreground NH vhx_dkpc = NULL; % ! in kpc vhx_systvel = 0.0; % ! in km/s vhx_nh22 = 0.0; % ! x10^22 /cm^2 % Fraction of the simulation that the source fills, % e.g., if it fills +/-1 half_ang degrees about the equatorial plane % then this fraction would be: sin(half_ang*(PI/180.0)); vhx_frac4pi = 1.0; % 1.0 for full-4pi % - - - CD and shocked region % User define the CD radius at the given-CD time step vhx_cd0step = NULL; vhx_cd0rad = NULL; % ! in pc % Temperature above-which material is considered shocked: vhx_Tshock = 1.e6; % ! in K % - - - Mass-cut related, incl. Clumping % The mass-cuts are defined by these two arrays: % (if 'mass=NULL; then they can be automatically set by vhx_make_cuts) shl_mass=NULL, shl_icmp=NULL; % Set mass-cut numbers for the vhx_make_cuts routine: vhx_ncutsej=12; vhx_ncutscsm=12; vhx_ncutsnear=3; % number of cuts close to the CD, 1 or more. % Indicate shells that will be clumped in their spectral model, % this has to be set after the mass-cuts are setup ('mass, 'icmp). % shl_clumped[ishl] = 0 for a not clumped shell; NULL for no clumping. shl_clumped=NULL; % The volume fraction and density factor for the clumping % (The values are ignored if no clumped shells are specified.) vhx_clmpfill=0.2; % volume filling fraction of the clumps vhx_clmpdens=5.0; % density factor for the clumps % - - - Abundances, Modeling, etc. % Reference abundance set, will be passed to XSPEC models too. % mxs_abund_pkg implements either: "angr" or "wilm". vhx_abundstr="angr"; % The nominal XSPEC abundance Zs: % H He C N O Ne Mg Si S Ar Ca Fe Ni % (The first abundzs must be 1=H, but the rest can % be just a subset of the zs, e.g. [1,8,26] for H,O,Fe only.) vhx_abundzs= [1,2, 6,7,8,10,12,14,16,18,20,26,28]; % Define the Ejecta and CSM relative abundances: % (an overall scale factor does not matter) % Set them to all 1.0: vhx_abundej= 1.0 + 0.0*vhx_abundzs; vhx_abundcsm= 1.0 + 0.0*vhx_abundzs; % Other model-related values; if they are NULL then the % vhx routines will not change them from what the user/isis % have set them to. % The XSPEC "NEIVERS" string for NEI models: vhx_xneiverstr = "1.1"; % "1.0" or "1.1" or "2.0" (or NULL) % The XSPEC cross-sections version: vhx_xxsectstr = NULL; % NULL or "vern" or "bcmc" % Beta-minimum value for tetp_equilib_bmin function % (will also use this beta if tetp_equilib_bmin is not available.) vhx_betamin = 0.06; % nominal value is 0.06 % - - - Plot options for the Locations plots % Plot the cummulative mass from outside-in (0) % or offset it to be zero at the CD (1): vhx_massrelcd = 0; % Overplot the density of a reference step - set step number % or NULL for no overplotting. vhx_refstep=NULL; % Can change the (log) Y-axis values if desired, % (though the axis is used for density, velocity, mass, % and temp, so usually a range of at least 10 to % 10^8 would be used) vhx_pltymin=0.02; vhx_pltymax=2.1e10; % Can set the X-axis min value (or NULL for auto-scaling) vhx_pltxmin=NULL; % ! in pc % Fix the X-axis max value (or NULL for auto-scaling) vhx_pltxmax=NULL; % ! in pc % - - - Energy bands and Spectral Plot % Arrays of E_lo and E_hi, the first one covers a wide range, % others set as desired. (Hard-coded to lengths of 3 for now.) vhx_fluxelos = [0.5, 0.5, 2.0]; vhx_fluxehis = [8.0, 2.0, 8.0]; % Min,Max for the Spectral plots (erg/s/cm^2/keV) vhx_fluxmin = 1.e-16; vhx_fluxmax = 1.e-8; % - - - Lightcurve plot option vhx_lc_xlog = 0; % linear years (0) or log days (1) vhx_lc_showejcsm = 0; % show (1) or not (0) separate ej,CSM contrib.s % End of User defined parameters % - - - - - - - % - - - - - - - % - - - % Variables defined below are filled/used by the routines, etc. % (generally not set directly by the user, although they are % available to the user.) % Starting time and time at each step, % these values are from the 'hst file. % Filled by vhx_setup() public variable vhx_starttime=NULL, vhx_steptime; % s % Values from a single time-step file % (all are in cgs units; except 'mass which is in solar-masses). % Filled by vhx_read_step(tid); public variable vhx_rad, vhx_dens, vhx_pres, vhx_vel, vhx_temp; public variable vhx_rad1; public variable vhx_tid=NULL, vhx_mass; % Variables to save the radius and density of a reference timestep % for overplotting; not valid/plotted if 'refstep == NULL. % Filled by vhx_read_step(tid) when tid==vhx_refstep. public variable vhx_refrad=NULL, vhx_refdens=NULL; % radial location of the CD, FS, RS for each step (in cm as is rad) % Filled by vhx_calc_cds(); public variable vhx_stepcd=NULL; public variable vhx_stepfs=NULL; public variable vhx_steprs=NULL; % shock mass from RS-to-CD and from CD-to-FS at each step - in M_solar - % Filled by vhx_calc_cds(); public variable vhx_stepmassrs=NULL; public variable vhx_stepmassfs=NULL; % Related to the abundance calc.s % Tcie in K is generally not changed so put it here rather than in user section above. public variable vhx_abundTcieK=1.e7; % Define values that are calculated from the abundances, % index of 1 is ejecta and 2 is csm (0 is unused, set to solar). public variable vhx_cmpninamu = Double_Type [3]; public variable vhx_cmpnenamu = Double_Type [3]; public variable vhx_cmpmu = Double_Type [3]; public variable vhx_cmpnhne = Double_Type [3]; % These variables apply to the mass-cut shells within a single time step. % (All shl_ variables are arrays.) See vhx_calc_shells for details. % Properties of the mass-cuts are: % Num_of_VH-1_bins, (outer)Radius[cm], v_kms, n_e, n_H, public variable shl_nvhbins, shl_rad, shl_vkms, shl_ne, shl_nh, shl_rhodr; % T[e]_keV, beta=Te/Ti public variable shl_TkeV, shl_beta; % Processing the time steps in order, the following values can be used to % carry/accumulate mass-cut history from step i-1 to step i for NEI purposes: % partial_shell, tau, public variable shl_partial, shl_tau, shl_Tdtau; % These values are set when vhx_make_par is run: public variable shl_norm; % These are filled by vhx_calc_shellfluxes: public variable shl_fluxes; % - - - % Routines are defined below % public define vhx_help () % Routine to present a simple list of vhx routines, etc. { message(" - - -"); message(" vh1xray.sl VH-1 X-ray Package version TBD, 2013-04-09"); message(""); message(" Summary of usage: > = isis command line prompt, [ ... ] are optional"); message(""); message(" > .load vh1xray"); message(" > .source mysim_values"); message(" > vhx_setup;"); message(" > [ vhx_calc_muetc; % print(vhx_cmpmu); ]"); message(" > [ vhx_list; ]"); message(" > [ vhx_read_step(tid); % tid is step number, e.g., 25 ]"); message(" > [ vhx_plot_step; ]"); message(" > [ vhx_cursor_cd0; % interactively set CD_0 location ]"); message(" > vhx_calc_cds;"); message(" > [ vhx_[ps]plot_cds; ]"); message(" > [ vhx_[ps]plot_mass; ]"); message(" > [ vhx_plot_allsteps; % writes locs_1NNN.gif's"+ ", can then make movie ]"); message(" > vhx_make_cuts;"); message(" > vhx_calc_shells; % writes shells_1NNN.txt's"); message(" > [ vhx_[ps]plot_kTtau(tid); % plot kT-tau from shells_1.txt ]"); message(" > [ vhx_make_par(tid); % writes model_1.par ]"); message(" > [ vhx_calc_shellfluxes(tid); % fills shl_fluxes[iEB,isl]"); message(" > vhx_calc_allspec; % writes '.par's, spectra_1NNN.gif's"+ ", and fluxlocs_vs_time_[Obs,NoNH].txt"); message(" > [ vhx_[ps]plot_lc; ]"); message(""); message(" See also: http://space.mit.edu/home/dd/VH1Xray/"); message(" - - -"); } % end of vhx_help public define vhx_calc_muetc () % Calculate plasma properites of the components from the abunds. % The abundances of the two components set the values of % electrons-per-amu, ions-per-amu, and H-per-ion % and let mu be calculated: % mu = number of amu per particle = 1 / (N_electrons/amu + N_ions/amu) % % Update/calculate these values from the provided values of: % vhx_abundstr, vhx_abundzs, vhx_abundej, vhx_abundcsm % and also: vhx_abundTcieK { % Start by setting the abundance-based values to angr or wilm solar values % default angr: vhx_cmpninamu = 0.775102 + [0.,0.,0.]; vhx_cmpnenamu = 0.852249 + [0.,0.,0.]; vhx_cmpnhne = 0.827312 + [0.,0.,0.]; % wilm? if (vhx_abundstr == "wilm") { vhx_cmpninamu = 0.780512 + [0.,0.,0.]; vhx_cmpnenamu = 0.855025 + [0.,0.,0.]; vhx_cmpnhne = 0.830845 + [0.,0.,0.]; } % % Now, do the real abunds calcs if mxs_abund_pkg is available: #ifexists mxs_abund_set mxs_abund_set(vhx_abundstr); vhx_cmpninamu[0]= 1.0/calc_amu_per_ion(vhx_abundzs, 1.0+0.0*vhx_abundzs); vhx_cmpninamu[1]= 1.0/calc_amu_per_ion(vhx_abundzs, vhx_abundej); vhx_cmpninamu[2]= 1.0/calc_amu_per_ion(vhx_abundzs, vhx_abundcsm); vhx_cmpnenamu[0]=calc_es_per_ion(vhx_abundzs, 1.0+0.0*vhx_abundzs, vhx_abundTcieK)* vhx_cmpninamu[0]; vhx_cmpnenamu[1]=calc_es_per_ion(vhx_abundzs, vhx_abundej, vhx_abundTcieK)* vhx_cmpninamu[1]; vhx_cmpnenamu[2]=calc_es_per_ion(vhx_abundzs, vhx_abundcsm, vhx_abundTcieK)* vhx_cmpninamu[2]; vhx_cmpnhne[0]= 1.0/calc_es_per_H(vhx_abundzs, 1.0+0.0*vhx_abundzs, vhx_abundTcieK); vhx_cmpnhne[1]= 1.0/calc_es_per_H(vhx_abundzs, vhx_abundej, vhx_abundTcieK); vhx_cmpnhne[2]= 1.0/calc_es_per_H(vhx_abundzs, vhx_abundcsm, vhx_abundTcieK); #endif % in any case, calculate mu (amu per particle): vhx_cmpmu = 1.0/(vhx_cmpnenamu + vhx_cmpninamu); } % end of vhx_calc_muetc public define vhx_list () % Routine to list the values of various vhx parameters, % for now it has most of the non-error messages from vhx_setup: { message(""); message(" Output will be in "+vhx_outdir+"/"); message(" History file: .../"+vhx_simdir+"/"+vhx_simhst); message(" Simulation steps (0-"+string(length(vhx_steptime)-1)+ ") go from "+string(vhx_steptime[0]/Const_yr)+ " to "+string(vhx_steptime[-1]/Const_yr)+ " years"); message(" Used steps: beginning--ending = "+string(vhx_simbeg)+"--"+ string(vhx_simend)+", every-nth = "+string(vhx_simnth)); message(" Using "+vhx_abundstr+"-referenced abundances;" +" mu_ejecta = "+string(vhx_cmpmu[1])+", mu_csm = "+string(vhx_cmpmu[2])); message(""); } % end of vhx_list public define vhx_setup () % Routine to setup various parameters, % check the simulation files, etc. % --> Because the user parameters are set before running this % none of those values should be changed here % (except for NULL replacements.) { variable fptr, fline, iline, charloc, tid; % Setup the output directory ()=system("mkdir -pv "+vhx_outdir); message(" Output will be in "+vhx_outdir+"/"); % Read the starting time and the times of all the steps % (variables are saved in units of seconds) % open the history file: fptr = fopen( vhx_abovedir + "/" + vhx_simdir + "/" + vhx_simhst, "r"); if (fptr == NULL) throw OpenError, " *** Could not open file: .../"+vhx_simdir+"/"+vhx_simhst; message(" History file: .../"+vhx_simdir+"/"+vhx_simhst); % go through the lines of the file iline=0; % line we're on for diagnostic purposes if needed. % Find the starting time vhx_starttime=NULL; while ( (-1 != fgets(&fline,fptr)) and vhx_starttime == NULL ) { iline += 1; %%message(string(iline)+": "+fline); % this line has the starting time? charloc=is_substr(fline, "Starting from time"); if (charloc != 0) { % found it %sscanf( substr(fline, charloc+20, -1), fmt, &vhx_starttime ); vhx_starttime = atof( substr(fline, charloc+20, -1) ); message(" Start time: "+string(vhx_starttime)+" yrs"); % convert to seconds: vhx_starttime *= Const_yr; } } % Continue through the file looking for the steps times, % set all times to invalid value (-1), can be up to 2000 of them: vhx_steptime = -1 + Double_Type [2000]; % load the values found into vhx_steptime[tid-1000] while ( -1 != fgets(&fline,fptr) ) { iline += 1; charloc=is_substr(fline, "Wrote "+vhx_simpre+"."); if (charloc != 0) { tid = int(atof( substr(fline, charloc+strlen("Wrote "+vhx_simpre+"."), 4) )) -1000; % get the time following the step number: vhx_steptime[tid] = atof( substr(fline, charloc+strlen("Wrote "+vhx_simpre+".nnnn,"), -1) ); % there could be a "time=" present as well, % if so then use that instead charloc=is_substr(fline, "time="); if (charloc != 0) { % get the time following "time=": vhx_steptime[tid] = atof( substr(fline, charloc+5, -1) ); } % add in the start time too: vhx_steptime[tid] += vhx_starttime; %%message(" tid = "+string(tid)+", time = "+ %% string(vhx_steptime[tid])); } } % close the file ()=fclose(fptr); % Reduce the step times to cover just 0 to the largest valid one vhx_steptime = vhx_steptime[[0:max(where(vhx_steptime > 0)):1]]; % check if any of these values are missing: if (min(vhx_steptime) < 0) { message(" *** vhx_setup: missing time steps in the file "+vhx_simhst); } message(" Simulation steps (0-"+string(length(vhx_steptime)-1)+ ") go from "+string(vhx_steptime[0]/Const_yr)+ " to "+string(vhx_steptime[-1]/Const_yr)+ " years"); % Setup the beg, end, nth step values - if they are not already set if (vhx_simbeg == NULL) vhx_simbeg = min(where(vhx_steptime > 0)); if (vhx_simend == NULL) vhx_simend = length(vhx_steptime)-1; % Keep number of plotted-etc steps below 100 if (vhx_simnth == NULL) vhx_simnth = 1 + int((vhx_simend-vhx_simbeg)/100.0); message(" Used steps: beginning--ending = "+string(vhx_simbeg)+"--"+ string(vhx_simend)+", every-nth = "+string(vhx_simnth)); % Calculate and show the mu values of ejecta and csm: vhx_calc_muetc; message(" Using "+vhx_abundstr+"-referenced abundances;" +" mu_ejecta = "+string(vhx_cmpmu[1])+", mu_csm = "+string(vhx_cmpmu[2])); % "Clear" the currently read-in time step vhx_tid=NULL; % "Clear" the reference rad and dens (not valid now) vhx_refrad=NULL; vhx_refdens=NULL; % "Clear" the CD, etc locations: vhx_stepcd=NULL; vhx_stepfs=NULL; vhx_steprs=NULL; vhx_stepmassrs=NULL; vhx_stepmassfs=NULL; } % end of vhx_setup public define vhx_read_step (tid) % Routine to read the values from the 5-column ouptut file into % the public varaibles (all in cgs units): % vhx_tid, vhx_rad, vhx_dens, vhx_pres, vhx_vel, vhx_temp; % For convenience, also setup an array vhx_rad1 which is the 'next' radius % i.e., zone j goes from vhx_rad[j] to vhx_rad1[j]. { variable vhxfile; % if tid is NULL, set vhx_tid to NULL and return if (tid == NULL) { vhx_tid=NULL; return; } % set the time step of these data: vhx_tid = @tid; vhxfile = vhx_abovedir+"/"+vhx_simdir+"/"+vhx_simpre+"."+string(1000+tid); (vhx_rad, vhx_dens, vhx_pres, vhx_vel, vhx_temp) = readcol(vhxfile,[1,2,3,4,5]); % Format of the VH-1 output data files: % The columns are: % Radius (cms), Density (g/cc), Pressure, Velocity (cms/sec), Temperature % 0.86983E+18 0.74111E-23 0.61268E-11 0.38385E+08 0.10000E+05 % 0.87374E+18 0.74111E-23 0.61268E-11 0.38558E+08 0.10000E+05 % % NOTE: The values of Dens, Press, Vel, Temp in a given row % are their average values over the zone from the row's Radius % to the next row's Radius. % For convenience, setup an array vhx_rad1 which is the 'next' radius % i.e., zone j goes from vhx_rad[j] to vhx_rad1[j]. variable jrad, rad_outer; % assume dr is ~ the same for the outer bin: jrad = length(vhx_rad)-1; rad_outer = 2.0*vhx_rad[jrad] - vhx_rad[jrad-1]; vhx_rad1 = [ vhx_rad[[1:length(vhx_rad)-1:1]], rad_outer ]; % if this step is 'simref then save the radius and density values % ('simref can be NULL to not use this feature) if (tid == vhx_refstep) { vhx_refrad=@vhx_rad; vhx_refdens=@vhx_dens; } % Form the cummulative mass from outside in; leave it in solar masses. % Note that 'dens[i] applies to the region from 'rad[i] to 'rad1[i]. % Start by evaluating dM for each dr bin. vhx_mass = 0.0*vhx_rad; % To avoid differencing cubed values use the approximation: % dM = dens * 4 pi r_ave^2 dr _for jrad (0,length(vhx_rad)-1,1) { vhx_mass[jrad] = vhx_dens[jrad]*4.*PI* (0.5*(vhx_rad1[jrad]+vhx_rad[jrad]))^2*(vhx_rad1[jrad]-vhx_rad[jrad]); } % Now convert to solar-masses and sum from out to in: vhx_mass = reverse(cumsum(reverse(vhx_mass/Const_SM))); } % end of vhx_read_step public define vhx_calc_cds () { % Fill the values of vhx_stepcd[tid] from 'beg to 'end. % (and also fill: 'steprs, 'stepfs, 'stepmassrs, 'stepmassfs) % Use the given-CD information that the CD is at vhx_cd0rad at % the given-CD time step vhx_cd0step. % % Note: if the vhx_stepcd values are later manually changed after % running vhx_calc_cds, then we'd want to re-do the calculation of % the 'stepfs,rs,massrs,massfs values. In that case % separating their calculation from this routine would be useful. % For now we'll assume that stepcd's are set here only. % variable tid, vcd, savetid, lastvel, sel, massatcd; message(""); message(" --- vhx_calc_cds ---"); % Are the 'cd0 step and rad set? if ((vhx_cd0step == NULL) or (vhx_cd0rad == NULL)) { message(""); message(" *** The values of vhx_cd0step and vhx_cd0rad are not both set;"); message(" to interactively set them, run: vhx_cursor_cd0; ***"); message(""); return; } % OK, continue... % Save the current tid value (to restore when done) savetid = @vhx_tid; % Define the arrays and reset the values to undefined (<0): vhx_stepcd = -1.0 + 0.0*vhx_steptime; vhx_stepfs = -1.0 + 0.0*vhx_steptime; vhx_steprs = -1.0 + 0.0*vhx_steptime; vhx_stepmassrs = -1.0 + 0.0*vhx_steptime; vhx_stepmassfs = -1.0 + 0.0*vhx_steptime; % Read the given-CD step tid = @vhx_cd0step; vhx_read_step(tid); % Set the given-CD location (user value is in pc, 'stepcd is in cm) % at the given-CD step: vhx_stepcd[vhx_cd0step] = vhx_cd0rad * Const_pc; % update the velocity at this CD for use in next step % Note: when interpolating in radius, % use the average of 'rad and 'rad1 since vel is the average of the zone. lastvel = interpol (vhx_stepcd[tid], 0.5*(vhx_rad+vhx_rad1), vhx_vel); % Also, set the RS and FS: the min and max locs with T > Tshock sel = where(vhx_temp > vhx_Tshock); if (length(sel) > 0) { % use rad1 when doing max: vhx_stepfs[tid]=vhx_rad1[max(sel)]; vhx_steprs[tid]=vhx_rad[min(sel)]; % and the RS-CD and CD-FS masses (both positive): % Note: when interpolating mass use rad since mass[i] % is cummulative from outside to rad[i]. massatcd = interpol (vhx_stepcd[tid], vhx_rad, vhx_mass); vhx_stepmassrs[tid] = interpol (vhx_steprs[tid], vhx_rad, vhx_mass) - massatcd; vhx_stepmassfs[tid] = massatcd - interpol (vhx_stepfs[tid], vhx_rad, vhx_mass); } % With the given-CD defined, use velocities and delta-t to % find the location of the CD at other time steps. % First, working to later time steps _for tid (vhx_cd0step+1, vhx_simend, 1) { % Read this new step vhx_read_step(tid); % Get the expected velocity of this step's CD % (i.e., at expected location given last vel and last CD) vcd = interpol ( (vhx_stepcd[tid-1] + lastvel*(vhx_steptime[tid]-vhx_steptime[tid-1])), 0.5*(vhx_rad+vhx_rad1), vhx_vel); % Set the new CD at previous plus v_ave times dt: vhx_stepcd[tid] = vhx_stepcd[tid-1] + 0.5*(lastvel+vcd)*(vhx_steptime[tid]-vhx_steptime[tid-1]); % update the new lastvel: lastvel = interpol (vhx_stepcd[tid], 0.5*(vhx_rad+vhx_rad1), vhx_vel); % Also, set the RS and FS: the min and max locs with T > Tshock sel = where(vhx_temp > vhx_Tshock); if (length(sel) > 0) { % use rad1 when doing max: vhx_stepfs[tid]=vhx_rad1[max(sel)]; vhx_steprs[tid]=vhx_rad[min(sel)]; % and the RS-CD and CD-FS masses (both positive): % Note: when interpolating mass use rad since mass[i] % is cummulative from outside to rad[i]. massatcd = interpol (vhx_stepcd[tid], vhx_rad, vhx_mass); vhx_stepmassrs[tid] = interpol (vhx_steprs[tid], vhx_rad, vhx_mass) - massatcd; vhx_stepmassfs[tid] = massatcd - interpol (vhx_stepfs[tid], vhx_rad, vhx_mass); } % show progress... if (tid == 20*int(tid/20.0)) message(" ... done through step "+string(tid)); } % Now, go back to the given-CD step and work to earlier time steps tid = @vhx_cd0step; vhx_read_step(tid); lastvel = interpol (vhx_stepcd[tid], 0.5*(vhx_rad+vhx_rad1), vhx_vel); _for tid (vhx_cd0step-1, vhx_simbeg, -1) { % Read the new (lower) step vhx_read_step(tid); % Get the expected velocity of this step's CD % (i.e., at expected location given last vel and last CD) vcd = interpol ( (vhx_stepcd[tid+1] - lastvel*(vhx_steptime[tid+1]-vhx_steptime[tid])), 0.5*(vhx_rad+vhx_rad1), vhx_vel); % Set the new CD at previous (higher) minus v_ave dt: vhx_stepcd[tid] = vhx_stepcd[tid+1] - 0.5*(vcd+lastvel)*(vhx_steptime[tid+1]-vhx_steptime[tid]); % update the new lastvel: lastvel = interpol (vhx_stepcd[tid], 0.5*(vhx_rad+vhx_rad1), vhx_vel); % Also, set the RS and FS: the min and max locs with T > Tshock sel = where(vhx_temp > vhx_Tshock); if (length(sel) > 0) { % use rad1 when doing max: vhx_stepfs[tid]=vhx_rad1[max(sel)]; vhx_steprs[tid]=vhx_rad[min(sel)]; % and the RS-CD and CD-FS masses (both positive): % Note: when interpolating mass use rad since mass[i] % is cummulative from outside to rad[i]. massatcd = interpol (vhx_stepcd[tid], vhx_rad, vhx_mass); vhx_stepmassrs[tid] = interpol (vhx_steprs[tid], vhx_rad, vhx_mass) - massatcd; vhx_stepmassfs[tid] = massatcd - interpol (vhx_stepfs[tid], vhx_rad, vhx_mass); } % show progress... if (tid == 20*int(tid/20.0)) message(" ... done through step "+string(tid)); } % restore the previous tid and its data (savetid can be NULL) vhx_read_step(savetid); } % end of vhx_calc_cds public define vhx_plot_step () % Make a plot of the values vs radius for the last read in % time step. { variable sel, pltxmin, pltxmax, pltymin, pltymax; variable massatcd=0, ishl, shlradcm, shlradpc, masscutclr; if (vhx_tid == NULL) { message(" *** vhx_plot_step: no valid step has been read in."); return; } % xrange can be set by user or auto-scaled; % it is linear. xlin; % nominal auto scaling, ! x-axis is in pc pltxmin = 0.98*min(vhx_rad/Const_pc); pltxmax = 1.02*max(vhx_rad1/Const_pc); % override that if values are given: if (vhx_pltxmin != NULL) pltxmin=@vhx_pltxmin; if (vhx_pltxmax != NULL) pltxmax=@vhx_pltxmax; % xrange(pltxmin, pltxmax); % yrange is log and covers from densities (<1) % to temps (>10^8). pltymin=@vhx_pltymin; % e.g., 0.02; pltymax=@vhx_pltymax; % e.g., 2.1e10; yrange(pltymin,pltymax); ylog; label("Radius (pc)", "Dens. [Blk,amu/cc]; vel [Grn,km/s]; M [Purp,uSM]; T [Red,K]", "Data from file: .../"+vhx_simdir+"/"+vhx_simpre+"."+string(1000+vhx_tid) + " Time = "+string(vhx_steptime[vhx_tid]/Const_yr)+" yr.s"); set_frame_line_width(3); charsize(1.1); % Density, amu/cm^3 set_line_width(9); plot(0.5*(vhx_rad+vhx_rad1)/Const_pc, vhx_dens/Const_amu, 1); % Underplot a "reference density" (if desired and defined and in-range) if ((vhx_refstep != NULL) and (vhx_refrad != NULL)) { if ((max(vhx_refrad/Const_pc) > pltxmin) and (min(vhx_refrad/Const_pc) < pltxmax)) { % include the reference density xlabel("Radius (pc) [Ref. density (gray) is step "+ string(vhx_refstep)+"]"); oplot(vhx_refrad/Const_pc, vhx_refdens/Const_amu, 15); oplot(0.5*(vhx_rad+vhx_rad1)/Const_pc, vhx_dens/Const_amu, 1); } } % Temp. K, emphasize the greater-than-shocked-level temps set_line_width(3); oplot(0.5*(vhx_rad+vhx_rad1)/Const_pc, vhx_temp,2); sel = where(vhx_temp > vhx_Tshock); % don't have any gaps in this selected region: sel = [min(sel):max(sel):1]; set_line_width(7); oplot(vhx_rad[sel]/Const_pc, vhx_temp[sel],2); % v km/s set_line_width(2); oplot(0.5*(vhx_rad+vhx_rad1)/Const_pc, vhx_vel*1.e-5,3); % Plot the cummulative Mass, in micro-M_solar units. % If the CD is defined and vhx_massrelcd is 1, % then subtract the mass at the CD % and plot the abs value, % i.e. the cummulative mass from the CD in each direction. if ((vhx_stepcd != NULL) and (vhx_massrelcd == 1)) { if (vhx_stepcd[vhx_tid] > 0) { massatcd = interpol ( vhx_stepcd[vhx_tid], vhx_rad, vhx_mass); } } set_line_width(1); oplot(vhx_rad/Const_pc, 1.e6*abs(vhx_mass-massatcd), 6); % Are the CDs defined? Is this one valid and in range? % plot it if so. if (vhx_stepcd != NULL) { % assuming that pltxmin is > 0, then just check this: if( (vhx_stepcd[vhx_tid]/Const_pc > pltxmin) and (vhx_stepcd[vhx_tid]/Const_pc < pltxmax)) { set_line_width(3); oplot(vhx_stepcd[vhx_tid]*[1.,1.]/Const_pc, [pltymin,pltymax], 14); } } % Plot the mass cuts on the plot, if defined if ((vhx_stepcd != NULL) and (shl_mass != NULL)) { % Make sure this CD is assigned/valid if (vhx_stepcd[vhx_tid] > 0) { massatcd = interpol ( vhx_stepcd[vhx_tid], vhx_rad, vhx_mass); shlradcm = interpol ( shl_mass, reverse(vhx_mass - massatcd), reverse(vhx_rad)); % only use the mass-cuts that are within the vhx_rad range: sel = where( (shlradcm > min(vhx_rad)) and (shlradcm < max(vhx_rad1)) ); % add the book-ends to these: if (min(sel) > 0) sel = [min(sel)-1, sel]; if (max(sel) < length(shl_mass)-1) sel = [sel, max(sel)+1]; _for ishl (0,length(sel)-1,1) { % mass-cut color, clumped or not-clumped: masscutclr=15; % not-clumped if (length(shl_clumped) == length (shl_icmp)) { if (shl_clumped[sel[ishl]] != 0) masscutclr=5; % clumped } % plot it if it is within the plot x range (in pc) shlradpc = shlradcm[sel[ishl]]/Const_pc; if ((shlradpc > pltxmin) and (shlradpc < pltxmax)) oplot(shlradpc+[0.,0.], [pltymin,10.0*pltymin],masscutclr); } } } } % end of vhx_plot_step public define vhx_plot_allsteps () % Go through the steps from 'beg to 'end and plot % every 'nth one of them. Write the files to the outdir. % Provide the "convert" command to make a movie from them; % (For mpg's, need ffmpeg, e.g.: apt-get install ffmpeg ) { variable tid, Win, progtid=10; message(""); message(" --- vhx_plot_allsteps ---"); _for tid (vhx_simbeg, vhx_simend, vhx_simnth) { % show progress... if (tid > progtid) { message(" ... done through step "+string(tid)); progtid = tid + 10*vhx_simnth -1; } % vhx_read_step(tid); Win=open_plot(vhx_outdir+"/locs_"+string(1000+tid)+".gif/GIF"); vhx_plot_step; close_plot(Win); } message(""); message(" Can make a movie from these files by doing:"); message("isis> ! convert -delay 10 "+vhx_outdir+"/locs_*.gif "+ vhx_outdir+"/mov_locs.gif \% or .mpg if ffmpeg supported"); message(""); message(" and view it with:"); message("isis> ! eog " + vhx_outdir+"/mov_locs.gif"); } % end of vhx_plot_allsteps public define vhx_plot_cds () % Plot the locations of the RS, CD, FS vs time { variable sel; xlin; xrange(0.0,); ylin; yrange(0.0,); label("Simulation Time (years)", "Radius of RS(red), CD(green), FS(blue) (pc)", "Data from files: .../"+vhx_simdir+"/"+vhx_simpre+".NNNN"); set_frame_line_width(3); charsize(1.1); set_line_width(5); sel=where(vhx_stepcd > 0); plot(vhx_steptime[sel]/Const_yr, vhx_stepfs[sel]/Const_pc, 4); oplot(vhx_steptime[sel]/Const_yr, vhx_steprs[sel]/Const_pc, 2); oplot(vhx_steptime[sel]/Const_yr, vhx_stepcd[sel]/Const_pc, 3); } % end of vhx_plot_cds % public define vhx_psplot_cds () % Do vhx_plot_cds with output to a .ps file in vhx_outdir { variable PSwin; % Do some check that the plot can be made if (vhx_stepfs == NULL) { message(" *** vhx_psplot_cds: No plot; didn't run vhx_calc_cds yet ?"); return; } PSwin=open_plot(vhx_outdir+"/cds_plot.ps/CPS"); vhx_plot_cds; close_plot(PSwin); } % end of vhx_psplot_cds public define vhx_plot_mass () % Plot the amount of shocked mass vs time { variable sel, maxmass, minmass, ishl, iclr; xlin; xrange(0.0,); ylog; maxmass = max([max(vhx_stepmassrs),max(vhx_stepmassfs)]); sel = where(vhx_stepmassrs > 0.0); minmass = min(vhx_stepmassrs[sel]); sel = where(vhx_stepmassfs > 0.0); minmass = min([ minmass, vhx_stepmassfs[sel] ]); yrange(0.5*minmass, 1.5*maxmass); label("Simulation Time (years)", "Shocked-Mass RS-CD(red), CD-FS(blue) (M_solar)", "Data from files: .../"+vhx_simdir+"/"+vhx_simpre+".NNNN"); set_frame_line_width(3); charsize(1.1); set_line_width(5); sel=where(vhx_stepcd > 0); plot(vhx_steptime[sel]/Const_yr, vhx_stepmassfs[sel], 4); oplot(vhx_steptime[sel]/Const_yr, vhx_stepmassrs[sel], 2); % and show the total shocked-mass: oplot(vhx_steptime[sel]/Const_yr, vhx_stepmassrs[sel]+vhx_stepmassfs[sel], 15); % Overplot the mass-cuts, if defined if (shl_mass != NULL) { set_line_width(1); _for ishl (0,length(shl_mass)-1,1) { iclr=7; % 6; % 13; % ejecta color if (shl_mass[ishl] < 0) iclr=5; % CSM ~ light blue oplot( [vhx_steptime[vhx_simbeg], vhx_steptime[vhx_simend]]/Const_yr, abs(shl_mass[ishl])+[0.,0.], iclr); } } } % end of vhx_plot_mass % public define vhx_psplot_mass () % Do vhx_plot_mass with output to a .ps file in vhx_outdir { variable PSwin; % Do some check that the plot can be made if (vhx_stepmassfs == NULL) { message(" *** vhx_psplot_mass: No plot; didn't run vhx_calc_cds yet ?"); return; } PSwin=open_plot(vhx_outdir+"/mass_plot.ps/CPS"); vhx_plot_mass; close_plot(PSwin); } % end of vhx_psplot_mass public define vhx_cursor_cd0 () % Plot the given-CD time step and have the user indicate the CD using the cursor. { variable tid, sel, deltar, xcur, ycur; variable savxmin, savxmax, saveref, savemassrel; % Use the provided cd0step if it is not NULL if (vhx_cd0step == NULL) { message(" vhx_cd0step is NULL, so auto-setting it ..."); % set it to a early time, ~ 20% from beginning to end tid = where(vhx_steptime > (0.80*vhx_steptime[vhx_simbeg] + 0.20*vhx_steptime[vhx_simend]))[0]; } else { message(" vhx_cd0step is defined, so will use that time step."); tid = vhx_cd0step; } message(" Showing time step "+string(tid)); % read it in and plot it around the shocked region vhx_read_step(tid); sel = where(vhx_temp > vhx_Tshock); % No ref density included on plot: saveref=@vhx_refstep; vhx_refstep=NULL; % Show the cummulative mass from outside (not rel to CD): savemassrel = @vhx_massrelcd; vhx_massrelcd = 0; % Save the current xrange settings savxmin = @vhx_pltxmin; savxmax = @vhx_pltxmax; % and zoom-in on the xrange that is around the shocked region deltar=vhx_rad[max(sel)]/Const_pc - vhx_rad[min(sel)]/Const_pc; vhx_pltxmin=vhx_rad[min(sel)]/Const_pc - 0.05*deltar; vhx_pltxmax=vhx_rad[max(sel)]/Const_pc + 0.05*deltar; vhx_plot_step; % activate the cursor to get the x-location of the CD from user message("Move cursor to be at the CD in X (Y doesn't matter), then"); cursor(&xcur,&ycur); oplot(xcur+[0.,0.],[0.1,1.e10],15); % summarize results: message(" This CD can be set by cutting/pasting these lines:"); message(""); message(" vhx_cd0step = "+string(tid)+";"); message(" vhx_cd0rad = "+string(xcur)+";"); message(""); message(" Note: these values are NOT automatically/now set."); message(" *** Copy and paste the 2 lines from above to set the new values. ***"); % put values back: vhx_pltxmin = @savxmin; vhx_pltxmax = @savxmax; vhx_refstep = @saveref; vhx_massrelcd = @savemassrel; } % end of vhx_cursor_cd0 public define vhx_make_cuts () % Setup the mass-cuts based on some logic. % [The values are also typed to the screen in case % the user wants to use them as a starting point % for custom explicit mass cuts.] { variable ncutsej, ncutscsm, mcfrac, csmmass, csmicmp; % number of major mass cuts for auto setup can be user set: ncutsej = @vhx_ncutsej; ncutscsm = @vhx_ncutscsm; % Check if values are already set ... if (shl_mass != NULL) { message(""); message(" *** The mass cuts are already set - replace them?"); message(" If so, then set: shl_mass = NULL;"); message(" and re-run this."); message(""); } % Check that things are ready for doing this step: if (vhx_stepmassrs == NULL) { % ooops - didn't run 'calc_cds ? message(""); message(" *** Need to run: vhx_calc_cds; ***"); message(""); return; } % The typical mass-cut specification, going from small to large % radii which is from postive to negative mass values (in % solar-masses, the CD is at mass=0), looks something like: % ejecta cuts (icmp = 1) : % shl_mass = [0.9, 0.7, 0.6, ... 0.1, 0.01]; % shl_icmp = [ 0, 1, 1, ... 1, 1 ]; % followed by CSM cuts (icmp = 2) : % shl_mass = [shl_mass, -0.01, -0.1, -0.2, ... -1.5, -2.0]; % shl_icmp = [shl_icmp, 0, 2, 2, ... 2, 2 ]; % Auto-generate a similar set of cuts, scaling to the max's of % the 'stepmassrs/fs values. % Ejecta: fill the shl_mass and 'icmp arrays mcfrac=1./(ncutsej+1); shl_mass = [1.0+0.6*mcfrac : mcfrac : -1.0*mcfrac]; shl_icmp = [ 0, 1+Integer_Type[ncutsej] ]; % add more closer to 0 ? shl_mass = [ shl_mass , 1.6*mcfrac/(1.6^[1:vhx_ncutsnear]) ]; shl_icmp = [ shl_icmp, 1+Integer_Type[vhx_ncutsnear] ]; % CSM: assemble in csmmass and csmicmp at first mcfrac=1./(ncutscsm+1); csmmass = [1.0+0.8*mcfrac : mcfrac : -1.0*mcfrac]; csmicmp = [ 0, 2+Integer_Type[ncutscsm] ]; % add more closer to 0 ? csmmass = [ csmmass , 1.6*mcfrac/(1.6^[1:vhx_ncutsnear]) ]; csmicmp = [ csmicmp, 2+Integer_Type[vhx_ncutsnear] ]; % Combine them, scaling for real masses % (reverse and negative the CSM masses, appended to ejecta): shl_mass = [ max(vhx_stepmassrs) * shl_mass, max(vhx_stepmassfs) * reverse( -1.0*csmmass )]; shl_icmp = [shl_icmp, csmicmp]; % Reset the clumped shells to all zeros % (since the number/location of the shells may have changed) shl_clumped = 0 * shl_icmp; } % end of vhx_make_cuts define _vhx_eval_shl (ishl) % This routine is used in the inner loop of vhx_calc_shells (the next % routine below); separated out to reduce the number of indent levels ;-) { variable vhsel, last_partial, shocked_weights, shockmass_weights; variable vshock, davetau, mass_part_last, mass_part_now, tauadj; % Find the vh1 bins that are fully within this mass-cut shell vhsel = where( (vhx_rad > shl_rad[ishl-1]) and (vhx_rad1 < shl_rad[ishl]) ); % number of vh-1 bins in this mass-cut shell: shl_nvhbins[ishl] = length(vhsel); % Do this shell if it has vh bins... if (length(vhsel) > 0) { % Calculate various mass-cut shell values... % Calculate the sum(rho dr) for these shells, amu/cm^2: % (using the 'shockmass_weights' variable, although the name is % not accurate and the variable will be re-used.) shockmass_weights = vhx_dens[vhsel] * (vhx_rad1[vhsel] - vhx_rad[vhsel]); % Save the sum of rho*dr for this mc shell, in amu/cm^2: shl_rhodr[ishl] = sum(shockmass_weights)/Const_amu; % - - - partial shocked fraction % Distinquish between four cases for a shell: % not shocked(0), partially shocked(<1.0), % just fully shocked(0.99), fully shocked(1.0). % and set its shl_partial[ishl] value (above). % Save the partial fraction of this shell from the previous % time step (might be useful) last_partial = shl_partial[ishl]; % Assign shocked/not values to the vh-1 bins, % makes array of 0.0's and 1.0's shocked_weights = double(vhx_temp[vhsel] > vhx_Tshock); % Check for none-shocked case (all 0.0's) if (sum(shocked_weights) < 0.5) shl_partial[ishl]=0.0; % Check for partial shocked case (some 0.0's, some 1.0's) if ( (min(shocked_weights) < 0.5) and (max(shocked_weights) > 0.5) ) { % Put the partially-shocked fraction into shl_parital. % What fraction of the radial distance (e.g. volume) is shocked ? shl_partial[ishl] = sum(shocked_weights* (vhx_rad1[vhsel] - vhx_rad[vhsel])) / sum(vhx_rad1[vhsel] - vhx_rad[vhsel]); } % Check for fully shocked case (all 1.0's) if (min(shocked_weights) > 0.9999) { shl_partial[ishl] = 1.0; % Check if this shell is just-now fully shocked: if (last_partial < 0.9) shl_partial[ishl] = 0.99; } % Is there shocked (X-ray emitting) material in this shell? % And it's not a component 0 shell? % If so, then calc quantities using shock-mass weighted values if ( (shl_partial[ishl] > 0.0) and (shl_icmp[ishl] != 0) ) { % Define the shockmass_weights of the shocked bins: (rho*dr) shockmass_weights = shocked_weights * vhx_dens[vhsel] * (vhx_rad1[vhsel] - vhx_rad[vhsel]) / sum(shocked_weights * vhx_dens[vhsel] * (vhx_rad1[vhsel] - vhx_rad[vhsel])); % Get the shocked-mass-weighted hydro values: shl_vkms[ishl] = 1.e-5*sum(shockmass_weights*vhx_vel[vhsel]); shl_ne[ishl] = vhx_cmpnenamu[shl_icmp[ishl]]* sum(shockmass_weights*vhx_dens[vhsel]/Const_amu); shl_nh[ishl] = vhx_cmpnhne[shl_icmp[ishl]]*shl_ne[ishl]; % Calculate the real temp, tau, etc. % Do the NEI mass-cut history tracking... % The shl_tau is the _average_ tau of the currently-shocked % shell (the distinction is most important in partially-shocked % shells.) % Change in tau is ne * dt: % Assume that step_start is 1 or larger so we have previous time: davetau = shl_ne[ishl] * (vhx_steptime[vhx_tid] - vhx_steptime[vhx_tid-1]); % If this shell is partially shocked % (or was last time so that now partial is 0.99) % then the average tau increase is different. % This is just a crude approximation... if (shl_partial[ishl] < 1.0) { % The newly-shocked material has an average tau ~ 0.5*ne*dt % tau_ave_newly = 0.5*davetau; % The previously shocked material (if any) % increases its ave tau by ne*dt, to have % tau_ave_prev = davetau + tau_previous; % or % tau_ave_prev = 0.5*davetau + 0.5*davetau + tau_previous; % The average tau of the combination of these is then: % tau_ave = f_n tau_ave_new + f_p tau_ave_prev % where f_n f_p are the mass fractions of material[*] % (each >= 0 and their sum is 1). So we can write: % tau_ave = f_n 0.5*davetau + % f_p (0.5*davetau + 0.5*davetau + tau_previous) % = 0.5*davetau + [using f_n + f_p = 1] % f_p (0.5*davetau + tau_previous) % The _change_ in the average Tau is then: % d_ave_tau = 0.5*davetau + % f_p (0.5*davetau + tau_previous) - tau_previous % In principle this says that the average tau can decrease... % % [*] The 'partial' values are fraction of the % shell radius (~volume) and the tau is better % averaged using fraction of mass. Assuming the shocked density % is x4 the unshocked then convert % partial_r to partial_mass via: % partial_mass ~ 4*partial_r / % ( (1-partial_r) + 4*partial_r ); % % Calculate and use the mass partial fractions: mass_part_last = 4.0*last_partial / ((1.0-last_partial)+4.0*last_partial); mass_part_now = 4.0*shl_partial[ishl] / ((1.0-shl_partial[ishl])+4.0*shl_partial[ishl]); % previously, calculated tau change as: %% davetau = ( ( (mass_part_now-mass_part_last)*0.5*davetau + %% mass_part_last*(davetau+shl_tau[ishl]) ) / %% mass_part_now ) - shl_tau[ishl]; % start with the correction to an 0.5*davetau change: tauadj = ( (mass_part_last/mass_part_now) * (0.5*davetau + shl_tau[ishl]) ) - shl_tau[ishl]; % don't reduce davetau below 0.5*davetau: if (tauadj < 0.0) tauadj = 0.0; davetau = 0.5*davetau + tauadj; } % update the tau shl_tau[ishl] = shl_tau[ishl] + davetau; % Start with the hydro(thermo) Temperature: % Derive the thermodynamic temp based on actual mu value: % kT [erg] = mu * amu * p / rho . % Mass-weight and shock-weight the p/rho; % then divide by 1.e3*Const_eV to get kT in keV. shl_TkeV[ishl] = Const_amu * vhx_cmpmu[shl_icmp[ishl]] * sum(shockmass_weights*vhx_pres[vhsel]/vhx_dens[vhsel])/ (1.e3*Const_eV); % Determine the self-consistant shock velocity: % kT = 3/16 mu amu vshock^2 or: % vshock^2 = (16/3) kT / (mu amu) [ convert to km/s ] vshock = 1.e-5 * sqrt( 16.0*shl_TkeV[ishl]*1.e3*Const_eV / (3.0*vhx_cmpmu[shl_icmp[ishl]]*Const_amu) ); shl_beta[ishl] = vhx_betamin; % default if we don't have: #ifexists tetp_equilib_bmin % modify the Temp to be the electron temp based on velocity and tau shl_beta[ishl] = tetp_equilib_bmin(vshock, shl_tau[ishl], vhx_betamin); #endif shl_TkeV[ishl] = shl_TkeV[ishl] * shl_beta[ishl]; % Update the tau-weighted-T sum, Tdtau: shl_Tdtau[ishl] = shl_Tdtau[ishl] + shl_TkeV[ishl]*davetau; } % end if there are shocked bins and not icmp=0 else { % Here if no shocked material or if component type = 0. % There are vh bins in the mass-cut, so do some calculations % and use simple mass weighting. % Define a mass-weighting for the bins: (rho*dr) % (using the 'shockmass_weights' variable, although the name is % not accurate.) shockmass_weights = vhx_dens[vhsel] * (vhx_rad1[vhsel] - vhx_rad[vhsel]); % Normalise the weights: shockmass_weights = shockmass_weights/sum(shockmass_weights); % Save the mass-weighted average over the vh bins of % velocity, in km/s shl_vkms[ishl] = 1.e-5*sum(shockmass_weights*vhx_vel[vhsel]); % Calculate the mass-weighted hydrogen density, amu/cm^3 shl_nh[ishl] = vhx_cmpnhne[shl_icmp[ishl]] * vhx_cmpnenamu[shl_icmp[ishl]] * sum(shockmass_weights*vhx_dens[vhsel])/Const_amu; % Be sure the history values are reset too shl_tau[ishl] = 1.0; shl_Tdtau[ishl] = 0.0; shl_partial[ishl] = 0.0; } } % end doing stuff if there are vh bins in the mc. else { % No bins (fully) in this mass cut % Reset the history-tracking values: shl_tau[ishl] = 1.0; shl_Tdtau[ishl] = 0.0; shl_partial[ishl] = 0.0; } } % end of routine _vhx_eval_shl(ishl) public define vhx_calc_shells () % Go through the desired time steps (one at a time) and % for each, calculate the properties of the set of mass-cut shells at that time. % Specifically, the values of these arrays are determined: % shl_nvhbins, shl_rad, shl_vkms, shl_ne, shl_nh, shl_TkeV, shl_beta; % Num_of_VH-1_bins, (outer)Radius[cm], v_kms, n_e, n_H, T_keV, beta=Te/Ti % Also set the following which are used to carry/accumulate mass-cut history from % mass-cut shell i-1 to shell i for NEI purposes: % shl_partial, shl_tau, shl_Tdtau; % partial_shell, tau, % % The results are saved by writing out to a file-per-time-step, % but only for the the [vhx_sim]nth time steps. % % Consider doing other time-step processing (.par file generation, % flux calculation, etc.) in this routine or in further routines. % { variable writetid, tid, massatcd, ishl, fp; message(""); message(" --- vhx_calc_shells ---"); % Check that the 'stepcd and shl_mass arrays are filled: if ((vhx_stepcd == NULL) or (shl_mass == NULL)) { message(""); message(" *** One or both of vhx_stepcd and shl_mass are not set."); message(" Consider running vhx_calc_cds; and/or vhx_make_cuts; *** "); message(""); return; } % Might also check that the 'stepcd values from 'simbeg to 'simend % are all valid (> 0). Or not: skip that for now. % Check for monotonic shl_mass values if ( min( shl_mass[[0:length(shl_mass)-2]] - shl_mass[[1:length(shl_mass)-1]] ) <= 0.0) { message(""); message(" *** The shl_mass array is not monotonic decreasing. ***"); message(""); return; } % In case shl_clumped is NULL, set it to all 0's: if (shl_clumped == NULL) shl_clumped = 0*shl_icmp; % Calculate/update the plasma properties of the components % based on the abundances vhx_calc_muetc; % Go from steps 'simbeg to 'simend one at a time, % and write out only the 'simnth ones. writetid = Integer_Type [vhx_simend+1]; % all zeroes. writetid[[vhx_simbeg:vhx_simend:vhx_simnth]] = 1; % set to 1 for output tids. % These 'history carrying' variables are initialized here before the tid loop: shl_partial = 0.0*abs(shl_mass); shl_tau = 1.0 + 0.0*abs(shl_mass); shl_Tdtau = 0.0*abs(shl_mass); _for tid (vhx_simbeg, vhx_simend, 1) { % Read this step vhx_read_step(tid); % Reset the shl_ values that will be calculated here; % they are all arrays with length of the mass-cuts. % The shells go from shl_rad[j-1] to shl_rad[j] with % the shell-averaged values at [j], e.g., shl_vkms[j]. shl_rad = 0.0*abs(shl_mass); shl_nvhbins = int(0.0*abs(shl_mass)); shl_vkms = 0.0*abs(shl_mass); % these get defined only when some/all of the shell is shocked: shl_ne = 0.0*abs(shl_mass); shl_nh = 0.0*abs(shl_mass); shl_rhodr = 0.0*abs(shl_mass); shl_TkeV = 0.0*abs(shl_mass); shl_beta = 0.0*abs(shl_mass); % Find the mass at this step's CD massatcd = interpol ( vhx_stepcd[tid], vhx_rad, vhx_mass); % map the mass-cut masses to the mass-cut shell radii shl_rad = interpol ( shl_mass, reverse(vhx_mass - massatcd), reverse(vhx_rad)); % Make a header file to add to the shells ouput values variable tempfp, hdrstr; % Make and add a header to the data file: tempfp=fopen("temp_header.txt", "w"); hdrstr="# Mass-cut shells output file from vhx_calc_shells()," + " defined in vh1xray.sl\n"; ()=fwrite(hdrstr, tempfp); hdrstr="# shell_index Num_of_VH-1_bins i_comp Radius[cm]" + " v_kms n_e n_H Te_keV" + " tau beta=Te/Ti PartialV Clumped\n"; ()=fwrite(hdrstr, tempfp); ()=fclose(tempfp); % Go through the mass-cut shells _for ishl (1,length(shl_mass)-1,1) { % Do the work for this mass-cut shell... _vhx_eval_shl(ishl); } % end loop over shells % The shl_ values are calculated, % output the shell values at this step? if (writetid[tid] == 1) { fp = fopen(vhx_outdir+"/shells_"+string(1000+tid)+".txt","w"); % Include a header ()=fprintf(fp,"# file: shells_"+string(1000+tid)+".txt" +" "+time()+ " \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp, "# shell_index Num_of_VH-1_bins i_comp Radius[cm]" + " v_kms n_e n_H Te_keV" + " tau beta=Te/Ti PartialV" + " Clumped rhodr \n"); % and write the values to the file writecol(fp, [0:length(shl_rad)-1], shl_nvhbins, shl_icmp, shl_rad, shl_vkms, shl_ne, shl_nh, shl_TkeV, shl_tau, shl_Tdtau/shl_tau, shl_beta, shl_partial, shl_clumped*[0:length(shl_rad)-1], shl_rhodr); ()=fclose(fp); } % show progress... if (tid == 20*int(tid/20.0)) message(" ... done through step "+string(tid)); } % end loop over tid } % end of vhx_calc_shells public define _vhx_read_shell (tid) % Read the shell values of a step back into the shl_ variables, % useful for other routines. % Return value is 0 if OK, 1 if a problem. { % Start by reading in the desired time step's shells values, % Does the shells file exist? if (system("ls "+vhx_outdir+"/shells_"+string(1000+tid)+".txt > /dev/null") != 0 ) { message(" *** The shells file is missing for time-step " + string(tid) + ". ***"); % Would the given 'beg, 'end, 'nth include this tid? if ( sum(tid == [vhx_simbeg:vhx_simend:vhx_simnth]) == 1) { % yes, so, suggest: message(" Create the shells files by running: vhx_calc_shells; "); } else { % This tid's shells would not have been made, so: message(" Change the vhx_simbeg, 'end, 'nth to include this step,"); message(" then (re)run the desired routine. "); } return 1; } % Read the values back into the shl_ arrays (the first column is skipped, an index) (shl_nvhbins, shl_icmp, shl_rad, shl_vkms, shl_ne, shl_nh, shl_TkeV, shl_tau, shl_Tdtau, shl_beta, shl_partial, shl_clumped) = readcol(vhx_outdir+"/shells_"+string(1000+tid)+".txt", [ [2:13:1] ]); return 0; } % end of _vhx_read_shell public define _vhx_calc_norms () { variable ishl; % Calculate the XSPEC norms for the shells shl_norm = 0.0*shl_rad; % can calculate them if distance is given, else leave them at 0 if (vhx_dkpc != NULL) { _for ishl (1, length(shl_rad)-1, 1) { % Start with n_e n_H x partial x Volume x 10^-48 : shl_norm[ishl] = shl_partial[ishl] * shl_ne[ishl] * shl_nh[ishl] * 4.0*PI*(0.5*(shl_rad[ishl] + shl_rad[ishl-1])*1.e-16)^2 * 1.e-16*(shl_rad[ishl] - shl_rad[ishl-1]); % Convert from n_e n_H V [10^-48] to % XSPEC norm = 1.e-14 * Volume * n_e n_H /(4 pi d^2) % = 1.e-14 * shl_norm * 1.e48 /(4 pi d^2) % = shl_norm / (4 pi (1.e-17*d)^2) % = shl_norm / (4 pi (1.e-17*1.e3*Const_pc*dkpc)^2) shl_norm[ishl] = shl_norm[ishl] / (4.0*PI*(1.e-14*Const_pc*vhx_dkpc)^2); } } } % end of _vhx_calc_norms public define vhx_plot_kTtau (tid) % Make a plot in kT-tau space of the mass-cut values % for step tid. { variable sel, ipt, ptsize, sumnorm, avenorm; % Read in the shell values if (_vhx_read_shell(tid)) return; % and the norms (could be all 0's if dkpc is not set) _vhx_calc_norms; % limit tau to 5e13 for these plotting purposes: sel = where(shl_tau > 5.e13); shl_tau[sel] = 5.e13; % Criteria to use a shell is as in 'make_par below, % Select shells that are shocked, have at least 1 vh bin, % and have temps above 0.08 keV: sel = where( (shl_partial > 0.0) and (shl_nvhbins > 0) and (shl_TkeV > 0.08) and (shl_Tdtau > 0.08) ); % set the point size for these ptsize = 1.5 + 0.0*sel; % default % adjust the size based on the norms sumnorm = sum(shl_norm[sel]); if (sumnorm > 0.0) { avenorm = sumnorm/length(sel); ptsize = 0.60 + 0.0*sel; % default, less than 0.1 of average ptsize[where(shl_norm[sel] > (0.10*avenorm))] = 0.80; % above 0.1 average ptsize[where(shl_norm[sel] > (0.30*avenorm))] = 1.10; % above 0.3 average ptsize[where(shl_norm[sel] > avenorm)] = 1.50; % above average ptsize[where(shl_norm[sel] > 2.5*avenorm)] = 2.0; % above average } %%message(" Showing plot of kT-Tau values for time step "+string(tid)); label("tau_i [s/cm^3]","kT_i [keV]", "Mass-cut kT-Tau values from " + vhx_outdir+"/shells_"+string(1000+tid)+".txt"); xrange(0.9e8,6.e13); xlog; yrange(0.1,100.0); ylog; set_frame_line_width(3); % show the track in kT-tau space connect_points(1); set_line_width(1); plot(shl_tau[sel], shl_TkeV[sel], 15); % show the discrete points connect_points(0); _for ipt (0,length(sel)-1,1) { % Components ejecta (1) and CSM (2) are colored black(1) and blue(4) color( 1 + 3*(int(0.5+shl_icmp[sel[ipt]])-1) ); point_size(ptsize[ipt]); pointstyle(-4); oplot(min([shl_tau[sel[ipt]],5.e13]), shl_TkeV[sel[ipt]]); % is this point clumped - if so, show the clumped point % fyi: clump parameters are: vhx_clmpfill and vhx_clmpdens . if (shl_clumped[sel[ipt]] != 0) { color(2); point_size(0.8*ptsize[ipt]); oplot(min([vhx_clmpdens*shl_tau[sel[ipt]],5.e13]), shl_TkeV[sel[ipt]]/vhx_clmpdens); } } pointstyle(-1); point_size(1.0); connect_points(1); % key to the color coding set_line_width(5); color(1); xylabel(2.0e9, 50.0,"Ejecta"); color(4); xylabel(2.5e9, 30.0, "CSM"); if (sum(shl_clumped[sel[ipt]]) > 0) { color(2); xylabel(10.0e9, sqrt(30.0*50.0), "Clumped"); } } % end of vhx_plot_kTtau % public define vhx_psplot_kTtau (tid) % Do vhx_plot_kTtau with output to a .ps file in vhx_outdir { variable PSwin; % Do some check that the plot can be made if (_vhx_read_shell(tid)) return; PSwin=open_plot(vhx_outdir+"/kTtau_"+string(1000+tid)+".ps/CPS"); vhx_plot_kTtau(tid); close_plot(PSwin); } % end of vhx_psplot_kTtau % Define two additional (dummy) model functions to provide extra model % parameters for use in the shell par files: define ejecta_fit (l,h,p) { return 1; } add_slang_function ("ejecta", ["norm", "muf", "cldens", "clfill", "H","He","C","N","O","Ne","Mg","Si", "S","Ar","Ca","Fe","Ni"]); define csm_fit (l,h,p) { return 1; } add_slang_function ("csm", ["norm", "muf", "cldens", "clfill", "H","He","C","N","O","Ne","Mg","Si", "S","Ar","Ca","Fe","Ni"]); % end additional isis models public define vhx_make_par (tid) % Generate an isis.par spectral model file for a time step, model_1.par. % This uses/requires the shells_1NNN.txt data file as input. % This routine may be useful but is generally not used directly, % instead it is called for the desired tid's when vhx_calc_allspec is run. % Note that while the shells'.txt file includes partial, clumping, % and NEI information, those could be variously ignored/overridden % at this step, if desired (e.g. no partial, no clumps, no NEI). % % This routine assumes (only) that the vhx values have been filled and % that vhx_setup has been run. ??? Is vhx_setup really required? ??? % % This routine also leaves the model defined and with ancillary % things setup (like xspec_xset command) so the model can be % evaluated if desired. { variable ishl, indx, ejmodstr, csmmodstr; variable sel, ej_vshell, csm_vshell; variable modstr, compstr, clumped; % Read in the shell values (return if a problem) if (_vhx_read_shell(tid)) return; % Make sure some parameters are set: if (vhx_dkpc == NULL) { message(" *** Need to set the distance, vhx_dkpc, to make par file."); return; } if (vhx_nh22 == NULL) { message(" *** The absorption, vhx_nh22, cannot be NULL, set to 0.0 for none."); return; } % Be sure we're using the desired model options (if non-NULL) if (vhx_xneiverstr != NULL) xspec_xset ("NEIVERS", vhx_xneiverstr); if (vhx_xxsectstr != NULL) () = xspec_xsect (vhx_xxsectstr); % Set the XSPEC abundance set to agree with the one used here % (if it is not already equal) if (xspec_abund() != vhx_abundstr) () = xspec_abund(vhx_abundstr); % Clip the TkeV, Tdtau, and tau values: shl_TkeV[where(shl_TkeV > 990.0)]=990.0; shl_Tdtau[where(shl_Tdtau > 990.0)]=990.0; shl_tau[where(shl_tau > 5.0e13)]=5.0e13; shl_tau[where(shl_tau < 6.25e9)]=6.25e9; % Calculate the XSPEC norms for the shells _vhx_calc_norms; % Create an isis model that is separated into ejecta and csm components % with the general form: % % phabs(1)*ejecta(1)*csm(1)*(gsmooth(1,()) + % gsmooth(2,())) % % Define/set the phabs, ejecta, csm values first: fit_fun("phabs(1)*ejecta(1)*csm(1)"); % The overall component norms, mu-factors, and clump info: % The norms are set to the 4pi fraction: set_par("ejecta(1).norm", vhx_frac4pi, 0, 0.0, 1.e2); set_par("csm(1).norm", vhx_frac4pi, 0, 0.0, 1.e2); set_par("ejecta(1).muf", 1.0, 1, 0.0, 1.e2); set_par("csm(1).muf", 1.0, 1, 0.0, 1.e2); % The clumping values - these are ignored except for clumped shells. set_par("ejecta(1).cldens", vhx_clmpdens, 1, 0.0, 1.e2); set_par("csm(1).cldens", vhx_clmpdens, 1, 0.0, 1.e2); set_par("ejecta(1).clfill", vhx_clmpfill, 1, 0.0, 1.0); set_par("csm(1).clfill", vhx_clmpfill, 1, 0.0, 1.0); % Absorption set_par("phabs(1).nH", vhx_nh22, 1); % Assemble the ejecta and CSM model strings ejmodstr=""; csmmodstr=""; % Keep track of the next model component to add: indx=0; _for ishl (1, length(shl_rad)-1, 1) { % Do shells that are shocked, have at least 1 vh bin, % and have temps above 0.08 keV: if ( (shl_partial[ishl] > 0.0) and (shl_nvhbins[ishl] > 0) and (shl_TkeV[ishl] > 0.08) and (shl_Tdtau[ishl] > 0.08) ) { % The next component number to add: %%indx = indx + 1; indx = 1*ishl; % use the shell number for the component index % determine if it is full(vgnei) or partial(vpshock): modstr = "vgnei"; if (shl_partial[ishl] < 0.999) modstr="vpshock"; % is it clumped ? (partial shells cannot be clumped) clumped=0; if ((shl_clumped[ishl] != 0) and (shl_partial[ishl] > 0.999)) clumped=1; % add this model to the ejecta/CSM string: % use the component index to select if it is ejecta or csm: if (shl_icmp[ishl] == 2) { compstr="csm"; if (csmmodstr == "") {csmmodstr = modstr+"("+string(indx)+")";} else {csmmodstr = csmmodstr + "+"+modstr+"("+string(indx)+")";} } else { compstr="ejecta"; if (ejmodstr == "") {ejmodstr = modstr+"("+string(indx)+")";} else {ejmodstr = ejmodstr + "+"+modstr+"("+string(indx)+")";} } % Define the function as of now: fit_fun("phabs(1)*ejecta(1)*csm(1)*(gsmooth(1,("+ejmodstr + "))+gsmooth(2,("+csmmodstr+")))"); % tie the abundances to the appropriate ejecta/csm values: tie(compstr+"(1).H",modstr+"("+string(indx)+").H"); tie(compstr+"(1).He",modstr+"("+string(indx)+").He"); tie(compstr+"(1).C",modstr+"("+string(indx)+").C"); tie(compstr+"(1).N",modstr+"("+string(indx)+").N"); tie(compstr+"(1).O",modstr+"("+string(indx)+").O"); tie(compstr+"(1).Ne",modstr+"("+string(indx)+").Ne"); tie(compstr+"(1).Mg",modstr+"("+string(indx)+").Mg"); tie(compstr+"(1).Si",modstr+"("+string(indx)+").Si"); tie(compstr+"(1).S",modstr+"("+string(indx)+").S"); tie(compstr+"(1).Ar",modstr+"("+string(indx)+").Ar"); tie(compstr+"(1).Ca",modstr+"("+string(indx)+").Ca"); tie(compstr+"(1).Fe",modstr+"("+string(indx)+").Fe"); tie(compstr+"(1).Ni",modstr+"("+string(indx)+").Ni"); % Set the redshift of this shell - depends on the % velocity/opacity properties of the emission model. % Assume optically thin and symmetric for now: set_par(modstr+"("+string(indx)+").Redshift", (vhx_systvel)/(Const_c*1.e-5), 1, -0.2,0.2); % Set the norm parameter value for this model component % including an overall ejecta/csm factor and clumping factor if clumped: if (clumped == 0) { set_par_fun (modstr+"("+string(indx)+").norm", string(shl_norm[ishl])+"*"+compstr+"(1).norm"); } else { set_par_fun (modstr+"("+string(indx)+").norm", string(shl_norm[ishl])+"*"+compstr + "(1).norm*(1.0-"+compstr+"(1).clfill)"); } % Set the kT, , and tau: % set updated limits and then include the appropriate mu-factor: set_par (modstr+"("+string(indx)+").kT", shl_TkeV[ishl], 1, 0.0808, 999.0); set_par_fun (modstr+"("+string(indx)+").kT", string(shl_TkeV[ishl])+"*"+compstr+"(1).muf"); % different parameters for partial shell: if (modstr == "vgnei") { % vgnei: % set the tau-averaged temp set_par (modstr+"("+string(indx)+").", shl_Tdtau[ishl], 1, 0.0808, 999.0); set_par_fun (modstr+"("+string(indx)+").", string(shl_Tdtau[ishl])+"*"+compstr+"(1).muf"); % set the tau set_par(modstr+"("+string(indx)+").Tau", min([shl_tau[ishl], 5.e13]), 1); } else { % vpshock: % set the two taus; the shl_tau is the average tau, % Can set Tau_u/l = 2x 0x shl_tau, % or, to avoid tau=0, use: 1.9x, 0.1x set_par(modstr+"("+string(indx)+").Tau_u", min([1.9*shl_tau[ishl], 5.e13]), 1); set_par(modstr+"("+string(indx)+").Tau_l", min([0.1*shl_tau[ishl], 5.e13]), 1); } % OK, if this shell is clumped then add another vgnei component with the % clumped parameter values... Copy/simplify approriate code from % above with changes/additions indicated at CLUMPED: lines. if (clumped == 1) { % CLUMPED: add clumped componet of type vgnei %%indx=indx+1; indx = indx + 500; modstr = "vgnei"; % add this model to the ejecta/CSM string: % use the component index to select if it is ejecta or csm: if (shl_icmp[ishl] == 2) { compstr="csm"; if (csmmodstr == "") {csmmodstr = modstr+"("+string(indx)+")";} else {csmmodstr = csmmodstr + "+"+modstr+"("+string(indx)+")";} } else { compstr="ejecta"; if (ejmodstr == "") {ejmodstr = modstr+"("+string(indx)+")";} else {ejmodstr = ejmodstr + "+"+modstr+"("+string(indx)+")";} } % Define the function as of now: fit_fun("phabs(1)*ejecta(1)*csm(1)*(gsmooth(1,("+ejmodstr + "))+gsmooth(2,("+csmmodstr+")))"); % tie the abundances to the appropriate ejecta/csm values: tie(compstr+"(1).H",modstr+"("+string(indx)+").H"); tie(compstr+"(1).He",modstr+"("+string(indx)+").He"); tie(compstr+"(1).C",modstr+"("+string(indx)+").C"); tie(compstr+"(1).N",modstr+"("+string(indx)+").N"); tie(compstr+"(1).O",modstr+"("+string(indx)+").O"); tie(compstr+"(1).Ne",modstr+"("+string(indx)+").Ne"); tie(compstr+"(1).Mg",modstr+"("+string(indx)+").Mg"); tie(compstr+"(1).Si",modstr+"("+string(indx)+").Si"); tie(compstr+"(1).S",modstr+"("+string(indx)+").S"); tie(compstr+"(1).Ar",modstr+"("+string(indx)+").Ar"); tie(compstr+"(1).Ca",modstr+"("+string(indx)+").Ca"); tie(compstr+"(1).Fe",modstr+"("+string(indx)+").Fe"); tie(compstr+"(1).Ni",modstr+"("+string(indx)+").Ni"); % CLUMPED: might want to use different redshift for clumps % in future... set_par(modstr+"("+string(indx)+").Redshift", (vhx_systvel)/(Const_c*1.e-5), 1, -0.2,0.2); % CLUMPED: Modify the norm for clumped: set_par_fun (modstr+"("+string(indx)+").norm", string(shl_norm[ishl])+"*"+compstr+"(1).norm*("+ compstr+"(1).clfill*"+compstr+"(1).cldens^2)"); % CLUMPED: Adjust the kT, , and tau for the clumped emission: % T's/dens and tau*dens set_par (modstr+"("+string(indx)+").kT", shl_TkeV[ishl], 1, 0.0808, 999.0); set_par_fun (modstr+"("+string(indx)+").kT", string(shl_TkeV[ishl])+"*"+ compstr+"(1).muf/"+compstr+"(1).cldens"); set_par (modstr+"("+string(indx)+").", shl_Tdtau[ishl], 1, 0.0808, 999.0); set_par_fun (modstr+"("+string(indx)+").", string(shl_Tdtau[ishl])+"*"+ compstr+"(1).muf/"+compstr+"(1).cldens"); set_par_fun (modstr+"("+string(indx)+").Tau", "min(["+string(shl_tau[ishl])+"*"+ compstr+"(1).cldens, 5.e13])"); } % end of adding clumped component } % end of shocked shell } % end of looping over shells for model components % Catch the cases where either or both of the ejecta,CSM have % no model components - set them to a zero-norm frozen model if (ejmodstr == "") { %%indx = indx + 1; indx = 901; ejmodstr = "vgnei("+string(indx)+")"; fit_fun(ejmodstr); % un-define it in case it had been: set_par_fun("vgnei("+string(indx)+").norm", NULL); set_par("vgnei("+string(indx)+").norm", 0.0, 1); freeze([2:18:1]); } if (csmmodstr == "") { %%indx = indx + 1; indx = 902; csmmodstr = "vgnei("+string(indx)+")"; fit_fun(csmmodstr); % un-define it in case it had been: set_par_fun("vgnei("+string(indx)+").norm", NULL); set_par("vgnei("+string(indx)+").norm", 0.0, 1); freeze([2:18:1]); } % Setup the whole model: fit_fun("phabs(1)*ejecta(1)*csm(1)*(gsmooth(1,("+ejmodstr + "))+gsmooth(2,("+csmmodstr+")))"); % Set the remaining global model parameters % Velocities and gsmooth parameters ej_vshell = 0.0; csm_vshell = 0.0; % defaults sel = where((shl_icmp == 1) and (shl_partial > 0.0) and (shl_nvhbins > 0)); if (length(sel) > 0) ej_vshell = sum(shl_norm[sel]*shl_vkms[sel]) / sum(shl_norm[sel]) ; % km/s sel = where((shl_icmp == 2) and (shl_partial > 0.0) and (shl_nvhbins > 0)); if (length(sel) > 0) csm_vshell = sum(shl_norm[sel]*shl_vkms[sel]) / sum(shl_norm[sel]) ; % km/s message(" time-step "+string(tid)+": Ave Ejecta,CSM shell velocities = " + string(ej_vshell)+", "+string(csm_vshell)+" km/s."); % Set the gsmooth values: % The conversion from the (radial) velocity to line-of-sight 1-sigma % blur velocity depends on the geometry of emission. % For a spherical shell the factor is: v_rms = 0.58 * v_radial % (For SN1987A: ring at 45 degrees and +/-15 degrees angle --> 0.505) set_par("gsmooth(1).Sig@6keV", 6.0 * 0.58 * abs(ej_vshell) / (Const_c*1.e-5), 1); set_par("gsmooth(2).Sig@6keV", 6.0 * 0.58 * abs(csm_vshell) / (Const_c*1.e-5), 1); % and the width scales as E^1 (dE/E is constant): set_par("gsmooth(1).Index", 1.0, 1); set_par("gsmooth(2).Index", 1.0, 1); % Abundances: % Set the abundances for ejecta(1) and csm(1), % the rest are already tied to these. % Require that vhx_abundzs array starts with 1(=H) % but it need not have all XSPEC z values in it, % (and could have zs not supported by XSPEC, so skip those.) % Mapping from z value to XSPEC z-string: variable xzs = [1,2,6,7,8,10,12,14,16,18,20,26,28]; variable xstrs = ["H","He","C","N","O","Ne","Mg","Si","S","Ar","Ca","Fe","Ni"]; variable ixz, iwhere; % Set the H abunds to 1 (assumed by the norm definition) set_par("ejecta(1).H", 1.0, 1, 0, 10); set_par("csm(1).H", 1.0, 1, 0, 10); % Go through the rest of the XSPEC zs and set them if preset, % otherwise set them to 0.0: _for ixz (1, length(xzs)-1, 1) { % Is this z in vhx_abundzs ? iwhere = where(vhx_abundzs == xzs[ixz]); if (length(iwhere) > 0) { % This z is here so set the appropriate ej/csm values, % normalized to the H value ('[0]). iwhere = iwhere[0]; set_par("ejecta(1)."+xstrs[ixz], vhx_abundej[iwhere]/vhx_abundej[0], 1, 0,1.e6); set_par("csm(1)."+xstrs[ixz], vhx_abundcsm[iwhere]/vhx_abundcsm[0], 1, 0,1.e6); } else { % Not there, set to zero: set_par("ejecta(1)."+xstrs[ixz], 0.0, 1, 0,1.e6); set_par("csm(1)."+xstrs[ixz], 0.0, 1, 0,1.e6); } } % DONE! Save the parameter file: save_par(vhx_outdir+"/model_"+string(1000+tid)+".par"); } % end of vhx_make_par(tid) public define vhx_calc_shellfluxes (tid) % Calculate the fluxes in the bands (hard-coded to 3 bands) % for each mass-cut shell at time step tid. % The shells'.txt and model'.par files are used. % The fluxes are put in shl_fluxes[iEnergyBand, ishl]. { variable nbins, lamlo, lamhi, fluxbin, hi, lo, iband, sel; variable ishl; variable szfbands = length(vhx_fluxelos); variable szshells, pltxmin, pltxmax, pltymax, Win, fp; % Read in the shl values: if( _vhx_read_shell(tid) ) return; message(" ... loading par file ..."); szshells = length(shl_rad); % Read in the par file for this step, load_par(vhx_outdir+"/model_"+string(1000+tid)+".par"); save_par("temp.par"); % from vh1xray - vhx_calc_allspec: % Setup a grid to evaluate the model on, it's in A with log spacing. % The range is set by the flux bands values and spaced at dE/E ~ 0.001: nbins = int(1000*log(max(vhx_fluxehis)/min(vhx_fluxelos))); (lamlo, lamhi) = linear_grid (log(Const_hc/max(vhx_fluxehis)), log(Const_hc/min(vhx_fluxelos)), nbins); lamlo = exp(lamlo); lamhi = exp(lamhi); % Evaluate the model fluxbin = eval_fun (lamlo, lamhi); % Make lo,hi in keV units: fluxbin = reverse(fluxbin); hi = reverse(Const_hc/lamlo); lo = reverse(Const_hc/lamhi); % Convert flux from photons/cm^2/s/bin to erg/s/cm^2/bin: % (Note that Const_eV is erg per eV) fluxbin = 0.5*(lo+hi) * 1.e3 * Const_eV * fluxbin; % Determine the flux (erg/s/cm^2) in the desired bands _for iband (0,szfbands-1,1) { sel = where( (lo > vhx_fluxelos[iband]) and (hi < vhx_fluxehis[iband]) ); message(" tid="+string(tid)+" : Flux-in-band["+string(iband)+ "] = "+string( sum(fluxbin[sel]) )); } % end loop over iband message(" ... calculating shell fluxes ..."); % Setup the output array shl_fluxes = Double_Type [szfbands, szshells]; variable modstr, clumped; % Do just the flux from each mass-cut shell _for ishl (1,szshells-1,1) { % Does this one have a flux-producing model component? % Same criteria as in make_par: if ( (shl_partial[ishl] > 0.0) and (shl_nvhbins[ishl] > 0) and (shl_TkeV[ishl] > 0.08) and (shl_Tdtau[ishl] > 0.08) ) { % determine if it is full(vgnei) or partial(vpshock): modstr = "vgnei"; if (shl_partial[ishl] < 0.999) modstr="vpshock"; % is it clumped ? (partial shells cannot be clumped) clumped=0; if ((shl_clumped[ishl] != 0) and (shl_partial[ishl] > 0.999)) clumped=1; % Setup the model of this single mass-cut shell if (clumped == 1) { fit_fun("phabs(1)*ejecta(1)*csm(1)*gsmooth("+ string(shl_icmp[ishl])+ ",vgnei("+string(ishl)+")+vgnei("+string(500+ishl)+"))"); } else { fit_fun("phabs(1)*ejecta(1)*csm(1)*gsmooth("+ string(shl_icmp[ishl])+","+modstr+"("+string(ishl)+"))"); } % Evaluate the model fluxbin = eval_fun (lamlo, lamhi); fluxbin = reverse(fluxbin); % Convert flux from photons/cm^2/s/bin to erg/s/cm^2/bin: % (Note that Const_eV is erg per eV) fluxbin = 0.5*(lo+hi) * 1.e3 * Const_eV * fluxbin; % Determine the flux (erg/s/cm^2) in the desired bands _for iband (0,szfbands-1,1) { sel = where( (lo > vhx_fluxelos[iband]) and (hi < vhx_fluxehis[iband]) ); shl_fluxes[iband,ishl]=sum(fluxbin[sel]); %%message(" tid="+string(tid)+" -- shell="+string(ishl)+ %% " : Flux-in-band["+string(iband)+ %% "] = "+string( shl_fluxes[iband,ishl] )); } % end loop over iband } % end of if shell has flux } % end of loop over ishl % Make a simple plot to show the shell fluxes %% xlin; xrange(0.0, % Set the xrange the same as in vhx_plot_step: % xrange can be set by user or auto-scaled; % it is linear. xlin; % nominal auto scaling, ! x-axis is in pc pltxmin = 0.0; pltxmax = 1.05*max(shl_rad[where(shl_fluxes[1,*] > 0.0)]/Const_pc); % override that if values are given: if (vhx_pltxmin != NULL) pltxmin=@vhx_pltxmin; if (vhx_pltxmax != NULL) pltxmax=@vhx_pltxmax; % xrange(pltxmin, pltxmax); % linear y-axis, from 0: %%ylin; yrange(0.0,); % log y-axis, from min to max over all bands: sel = shl_fluxes[where(shl_fluxes > 0.0)]; pltymax = 1.2*max(sel); % at least 2 decades: yrange(min([sel,0.01*pltymax]), pltymax); ylog; label("Radius (pc)","Flux-from-shell (erg/s/cm^2)", "Based on shells par-file: "+vhx_outdir+"/model_"+string(1000+tid)+".par"); plot_bin_integral; hplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[0,[1:szshells-1:1]],1); ohplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[1,[1:szshells-1:1]],2); ohplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[2,[1:szshells-1:1]],4); % Note the energy ranges on the plot: color(1); xylabel(0.9*pltxmin+0.1*pltxmax, 0.3*pltymax, string(vhx_fluxelos[0])+" -- "+string(vhx_fluxehis[0])+" keV"); color(2); xylabel(0.9*pltxmin+0.1*pltxmax, 0.6*0.3*pltymax, string(vhx_fluxelos[1])+" -- "+string(vhx_fluxehis[1])+" keV"); color(4); xylabel(0.9*pltxmin+0.1*pltxmax, 0.6*0.6*0.3*pltymax, string(vhx_fluxelos[2])+" -- "+string(vhx_fluxehis[2])+" keV"); % and write the plot to a file as well Win=open_plot(vhx_outdir+"/shellfluxes_"+string(1000+tid)+".ps/CPS"); plot_bin_integral; hplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[0,[1:szshells-1:1]],1); ohplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[1,[1:szshells-1:1]],2); ohplot(shl_rad[[0:szshells-2:1]]/Const_pc, shl_rad[[1:szshells-1:1]]/Const_pc, shl_fluxes[2,[1:szshells-1:1]],4); % Note the energy ranges on the plot: color(1); xylabel(0.9*pltxmin+0.1*pltxmax, 0.3*pltymax, string(vhx_fluxelos[0])+" -- "+string(vhx_fluxehis[0])+" keV"); color(2); xylabel(0.9*pltxmin+0.1*pltxmax, 0.6*0.3*pltymax, string(vhx_fluxelos[1])+" -- "+string(vhx_fluxehis[1])+" keV"); color(4); xylabel(0.9*pltxmin+0.1*pltxmax, 0.6*0.6*0.3*pltymax, string(vhx_fluxelos[2])+" -- "+string(vhx_fluxehis[2])+" keV"); close_plot(Win); % output the values to a file: % open the output data file fp = fopen(vhx_outdir+"/shellfluxes_"+string(1000+tid)+".txt","w"); % Include a useful header, e.g, ()=fprintf(fp,"# file: shellfluxes_"+string(1000+tid)+".txt" +" "+time()+ " \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# X-ray Fluxes in units of erg/s/cm^2 \n"); ()=fprintf(fp,"# Energy bands (keV): \n"); ()=fprintf(fp,"# Band 0 is "+string(vhx_fluxelos[0])+" to "+ string(vhx_fluxehis[0])+" \n"); ()=fprintf(fp,"# Band 1 is "+string(vhx_fluxelos[1])+" to "+ string(vhx_fluxehis[1])+" \n"); ()=fprintf(fp,"# Band 2 is "+string(vhx_fluxelos[2])+" to "+ string(vhx_fluxehis[2])+" \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# shell_index Flux_0 Flux_1 Flux_2 \n"); % add the data values: writecol(fp, [0:length(shl_rad)-1], shl_fluxes[0,*], shl_fluxes[1,*], shl_fluxes[2,*]); ()=fclose(fp); } % end of vhx_calc_shellfluxes public define vhx_calc_allspec () % Make and evaluate all the desired steps' par files, % write spectral plots for each, and write a fluxlocs-vs-time data file % (spectrum_1NNN.gif's and fluxlocs_vs_time.txt). % This can be (re)run following only vhx_setup provided vhx_calc_shells % had been previously run (to make the shells'.txt files.) { variable sel, weights, fp, modstr, nabs; variable itid, iband, icmp, lamlo, lamhi, lo, hi, nbins, fluxbin; variable Win; % Define arrays over the set of desired timesteps, these % will be filled with flux info, etc. % Sizes of the arrays: variable szcmp = 3; variable step_tids = [vhx_simbeg:vhx_simend:vhx_simnth]; variable szstep = length(step_tids); variable szfbands = length(vhx_fluxelos); % Setup the flux-locs output arrays: variable step_times = Double_Type [szstep], step_vels = Double_Type [szcmp, szstep], step_temps = Double_Type [szcmp, szstep]; variable step_rss = Double_Type [szstep]; variable step_cds = 0.0*step_rss, step_fss = 0.0*step_rss; variable step_fluxes = Double_Type [szcmp, szfbands, szstep]; message(""); message(" --- vhx_calc_allspec ---"); % Make sure the step times (etc) are available if (length(vhx_steptime) < 2) { message(" *** The step times are not defined, do: vhx_setup;"); return; } % and set them step_times = vhx_steptime[step_tids]; % Loop over the desired time steps: _for itid (0, length(step_tids)-1,1) { % Make (and load) the par file: vhx_make_par(step_tids[itid]); % This has also loaded the shl_ variables so we can calculate % average ejecta and CSM values of vel and temp, and locations % of RS, CD, FS. % CD location: sel = where(shl_icmp ==0)[1]; step_cds[itid] = 0.5*(shl_rad[sel] + shl_rad[sel-1]); % Require shocked ejecta shells for these: sel = where((shl_icmp == 1) and (shl_partial > 0.0) and (shl_nvhbins > 0)); if (length(sel) > 0) { weights = shl_norm[sel] / sum(shl_norm[sel]); step_vels[1,itid] = sum( shl_vkms[sel] * weights ); step_temps[1,itid] = sum( shl_TkeV[sel] * weights ); % the RS location step_rss[itid] = shl_rad[min(sel)-1]; } % Require shocked CSM shells for these: sel = where((shl_icmp == 2) and (shl_partial > 0.0) and (shl_nvhbins > 0)); if (length(sel) > 0) { weights = shl_norm[sel] / sum(shl_norm[sel]); step_vels[2,itid] = sum( shl_vkms[sel] * weights ); step_temps[2,itid] = sum( shl_TkeV[sel] * weights ); % the FS location step_fss[itid] = shl_rad[max(sel)]; } % Evaluate the fluxes in the energy bands from ejecta and CSM, % filling the array values: step_fluxes[icmp, iband, itid] % Setup a grid to evaluate the model on, it's in A with log spacing. % The range is set by the flux bands values and spaced at dE/E ~ 0.001: nbins = int(1000*log(max(vhx_fluxehis)/min(vhx_fluxelos))); (lamlo, lamhi) = linear_grid (log(Const_hc/max(vhx_fluxehis)), log(Const_hc/min(vhx_fluxelos)), nbins); lamlo = exp(lamlo); lamhi = exp(lamhi); % Array for the component spectra for plotting variable cmpspec = Double_Type [3, length(lamlo)]; % Loop over the non-0 components (ejecta, CSM) _for icmp (1,2,1) { % Start with the nominal norm values, vhx_frac4pi: set_par("ejecta(1).norm", vhx_frac4pi); set_par("csm(1).norm", vhx_frac4pi); % Set the other component's norm to 0 if (icmp == 1) set_par("csm(1).norm", 0.0); if (icmp == 2) set_par("ejecta(1).norm", 0.0); % Evaluate the model fluxbin = eval_fun (lamlo, lamhi); % Make lo,hi in keV units: fluxbin = reverse(fluxbin); hi = reverse(Const_hc/lamlo); lo = reverse(Const_hc/lamhi); % Convert flux from photons/cm^2/s/bin to erg/s/cm^2/bin: % (Note that Const_eV is erg per eV) fluxbin = 0.5*(lo+hi) * 1.e3 * Const_eV * fluxbin; % Save this component's spectrum (erg/s/cm^2/keV) cmpspec[icmp,*] = fluxbin/(hi-lo); % Determine the flux (erg/s/cm^2) in the desired bands _for iband (0,szfbands-1,1) { sel = where( (lo > vhx_fluxelos[iband]) and (hi < vhx_fluxehis[iband]) ); step_fluxes[icmp, iband, itid] = sum(fluxbin[sel]); } % end loop over iband } % end loop over icmp % Scale to full 4pi and convert to units of 10^30 erg/s: step_fluxes[*,*,itid] *= 4.0 * PI * (vhx_dkpc*1.e3*Const_pc*1.e-15)^2; % Put the sum of ejecta,CSM flux values into the icmp=0 values: step_fluxes[0, *, itid] = step_fluxes[1, *, itid] + step_fluxes[2, *, itid]; message(" --> Total flux ("+string(vhx_fluxelos[0]) + "-"+string(vhx_fluxehis[0])+" keV) = " + sprintf("%.3e", step_fluxes[0,0,itid])+" x10^30 erg/s" ); % Form the ejecta+CSM spectrum and clip the spectra at 'fluxmin cmpspec[0,*] = cmpspec[1,*] + cmpspec[2,*]; sel = where( cmpspec < vhx_fluxmin ); if (length(sel) > 0) cmpspec[sel] = vhx_fluxmin; % Make a plot of the spectrum at this time Win=open_plot(vhx_outdir+"/spectra_"+ string(1000+step_tids[itid])+".gif/GIF"); label("Energy (keV)","Flux [erg/s/cm^2/keV]", "Flux ("+string(vhx_fluxelos[0])+"-"+string(vhx_fluxehis[0])+ " keV) = "+sprintf("%.3e", 1.e30*step_fluxes[0,0,itid])+" erg/s" + " ("+vhx_simdir+") Time["+string(step_tids[itid])+"] = " + string(step_times[itid]/Const_yr)+" yr.s"); xrange(vhx_fluxelos[0],vhx_fluxehis[0]); xlog; yrange(vhx_fluxmin,vhx_fluxmax); ylog; set_frame_line_width(3); charsize(1.0); % Total spectrum plot (will over plot it...) color(1); set_line_width(1); hplot(lo,hi,cmpspec[0,*]); #ifexists oplot_common_HHe % Show the location of reference H,He-like lines set_line_width(1); oplot_common_HHe(vhx_fluxelos[0],vhx_fluxehis[0], 2.0*vhx_fluxmin,0.5*vhx_fluxmax, 1, 0); #endif % comp 1 color(1); set_line_width(5); ohplot(lo,hi,cmpspec[1,*]); % comp 2 color(4); set_line_width(5); ohplot(lo,hi,cmpspec[2,*]); % Total spectrum plot: color(2); set_line_width(3); ohplot(lo,hi,cmpspec[0,*]); close_plot(Win); } % end looping over step_tid values % Write out the Flux-Locations output file % % Modify the filename to indicate different flux calculations? % Note if no (or very low) NH is used: modstr="_Obs"; % default is "observed" fluxes if (vhx_nh22 < 1.e-6) modstr = "_NoNH"; % the model has no NH absorption % open the output data file fp = fopen(vhx_outdir+"/fluxlocs_vs_time"+modstr+".txt","w"); % Include a useful header, e.g, ()=fprintf(fp,"# file: fluxlocs_vs_time"+modstr+".txt" +" "+time()+ " \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# Fluxes and locations vs time for hydro model: "+vhx_simdir+" \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# Distance of "+ string(vhx_dkpc)+" kpc. \n"); ()=fprintf(fp,"# Fluxes include NH22 = "+string(vhx_nh22)+"\n"); ()=fprintf(fp,"# XSPEC NEIVERS = "+string(vhx_xneiverstr)+"\n"); ()=fprintf(fp,"# \n"); % Note: shl_ values will be those of the last vhx_make_par call ()=fprintf(fp,"# Number of mass-cut shells: "+ string(length(where( (shl_icmp == 1) ))) + " (ejecta), "+ string(length(where( (shl_icmp == 2) )))+" (CSM) \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# Number of clumped shells in ejecta, CSM = "+ string(sum( (shl_clumped > 0) and (shl_icmp == 1) )) + ", "+string(sum( (shl_clumped > 0) and (shl_icmp == 2) ))+"\n" ); ()=fprintf(fp,"# Clumping parameters: fill-fraction, density = "+ string(vhx_clmpfill)+", "+string(vhx_clmpdens)+"\n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# Abundances (w.r.t "+vhx_abundstr+"): \n"); ()=fprintf(fp,"# H He C N O Ne Mg Si" + " S Ar Ca Fe Ni \n"); nabs = 1.e-4+int(1000.0*vhx_abundej/max(1.0*vhx_abundej[[0:3:1]]))/1000.0; ()=fprintf(fp,"# Ejecta: ["+sprintf("%S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S", nabs[0],nabs[1],nabs[2],nabs[3],nabs[4],nabs[5], nabs[6],nabs[7],nabs[8],nabs[9], nabs[10],nabs[11],nabs[12])+"]; \n"); nabs = 1.e-4+int(1000.0*vhx_abundcsm/max(1.0*vhx_abundcsm[[0:3:1]]))/1000.0; ()=fprintf(fp,"# CSM: ["+sprintf("%S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S, %S", nabs[0],nabs[1],nabs[2],nabs[3],nabs[4],nabs[5], nabs[6],nabs[7],nabs[8],nabs[9], nabs[10],nabs[11],nabs[12])+"]; \n"); ()=fprintf(fp,"# mu_Ejecta = "+string(vhx_cmpmu[1])+" , mu_CSM = "+string(vhx_cmpmu[2])+" \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# Column Explanations (besides self-explanatory ones): \n"); ()=fprintf(fp,"# F_Ejecta[CSM]_Lo{Hi} : Ejecta[CSM] X-ray Flux in units of 10^30 erg/s \n"); ()=fprintf(fp,"# Lo band is "+string(vhx_fluxelos[1])+" to "+ string(vhx_fluxehis[1])+" keV \n"); ()=fprintf(fp,"# Hi band is "+string(vhx_fluxelos[2])+" to "+ string(vhx_fluxehis[2])+" keV \n"); ()=fprintf(fp,"# \n"); ()=fprintf(fp,"# SimNumber SimTime(s) F_Ejecta_Lo F_Ejecta_Hi F_CSM_Lo" + " F_CSM_Hi r_RS(cm) r_CD(cm) r_FS(cm)" + " v_Ejecta(km/s) v_CSM(km/s) Te_Ejecta(K) Te_CSM(K) \n"); ()=fclose(fp); % % And append/write the data values: fp = fopen(vhx_outdir+"/fluxlocs_vs_time"+modstr+".txt","a"); writecol(fp, 1000+step_tids[*], step_times[*], % Ejecta fluxes step_fluxes[1,1,*], step_fluxes[1,2,*], % CSM fluxes step_fluxes[2,1,*], step_fluxes[2,2,*], % locations step_rss[*], step_cds[*], step_fss[*], % ave velocities and temps step_vels[1,*], step_vels[2,*], step_temps[1,*], step_temps[2,*]); ()=fclose(fp); % Helpful info message(""); message(" Can make a movie from the Spectral files by doing:"); message("isis> ! convert -delay 10 "+vhx_outdir+"/spectra_*.gif "+ vhx_outdir+"/mov_spectra.gif \% or .mpg if ffmpeg supported"); message(""); message(" and view it with:"); message("isis> ! eog " + vhx_outdir+"/mov_spectra.gif"); } % end of vhx_calc_allspec public define vhx_plot_lc () % Plot the fluxes-in-bands vs time (from values in the fluxlocs' file); { variable modstr, yunit, yscale, pltymin, sel; variable xunit, xscale; % setup variables to read into variable step_tids, step_times, % Ejecta fluxes step_flux_ejlo, step_flux_ejhi, % CSM fluxes step_flux_csmlo, step_flux_csmhi, % locations step_rss, step_cds, step_fss, % ave velocities and temps step_vels_ej, step_vels_csm, step_temps_ej, step_temps_csm; % Read the fluxlocs file back in - choose from _Obs or _NoNH % (provided they have been made.) modstr="_Obs"; % default is "observed" fluxes if (vhx_nh22 < 1.e-6) modstr = "_NoNH"; % the model has no NH absorption % do the read (step_tids, step_times, step_flux_ejlo, step_flux_ejhi, step_flux_csmlo, step_flux_csmhi, step_rss, step_cds, step_fss, step_vels_ej, step_vels_csm, step_temps_ej, step_temps_csm) = readcol(vhx_outdir+"/fluxlocs_vs_time"+modstr+".txt", [1:13:1]); step_tids = step_tids - 1000; % Use erg/s or convert to erg/s/cm^2 depending on the file type: % default erg/s: yunit = "10^30 erg/s"; yscale = 1.0; if (modstr == "_Obs") { yunit="erg/s/cm^2"; yscale=1.0e30/(4.0*PI*(vhx_dkpc*1.e3*Const_pc)^2); } % Set the minimum Y-value below the min total flux: pltymin = yscale*(step_flux_ejlo + step_flux_csmlo + step_flux_ejhi + step_flux_csmhi); sel = where(pltymin > 0.0); pltymin = 0.1*min(pltymin[sel]); yrange(pltymin,); ylog; % X-axis % The step_times are in seconds, scale to years or days xunit="years"; xscale=1.0/Const_yr; xrange; xlin; if (vhx_lc_xlog == 1) { xunit="days"; xscale=1.0/(24.0*3600.0); xlog; xrange; % use log if it is days } label("Simulation time ["+xunit+"]","Flux in Bands, low (orange), high (blue)" + " ["+yunit+"]","Lightcurves from simulation "+vhx_simdir); set_frame_line_width(3); charsize(1.1); % Plot "total" flux - this is the sum of ejecta and CSM for the two % bands, because the bands can cover less of more of the total energy % range (e.g.s: 0.5-2, 3-8; or 0.3-2.4, 2-8), this "total" is % not meaningful, so plot it invisible (color=0). set_line_width(1); plot(xscale*step_times, yscale*(step_flux_ejlo + step_flux_csmlo + step_flux_ejhi + step_flux_csmhi ), 0); % Plot lo flux set_line_width(5); oplot(xscale*step_times, yscale*(step_flux_ejlo + step_flux_csmlo), 8); % Plot hi flux oplot(xscale*step_times, yscale*(step_flux_ejhi + step_flux_csmhi), 5); if (vhx_lc_showejcsm == 1) { % Show the individual ejecta (dotted) and CSM (dashed) fluxes too % linestyle: 2 = dashed, 4 = dotted linestyle(2); oplot(xscale*step_times, yscale*(step_flux_csmlo), 8); linestyle(2); oplot(xscale*step_times, yscale*(step_flux_csmhi), 5); linestyle(4); oplot(xscale*step_times, yscale*(step_flux_ejlo), 8); linestyle(4); oplot(xscale*step_times, yscale*(step_flux_ejhi), 5); linestyle(1); } } % end of vhx_plot_lc % public define vhx_psplot_lc () % Do vhx_plot_lc with output to a .ps file in vhx_outdir { variable PSwin, modstr; % Do some check that the plot can be made %%if (...) %% { %% message(" *** vhx_psplot_lc: No plot; didn't run vhx_calc_cds yet ?"); %% return; %% } % Two possible plots: modstr="_Obs"; % default is "observed" fluxes if (vhx_nh22 < 1.e-6) modstr = "_NoNH"; % the model has no NH absorption PSwin=open_plot(vhx_outdir+"/lc"+modstr+"_plot.ps/CPS"); vhx_plot_lc; close_plot(PSwin); } % end of vhx_psplot_lc % End of routines % - - - % - - - Finally, show the vhx_help command and its output message(" > vhx_help;"); vhx_help; % - - - End of vh1xray.sl - - -