% % file: i3d_sphxcyl.sl % % This file adds an ISIS fit function, sphxcyl_line, to the models % available for fitting. % % This also demonstrates spatial-spectral 3D model fitting in ISIS. % % For more information see the comments in routine below. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % It can be useful to skip re-evaluating v3d and MC things if the % geometry and velocity parameters are the same from the previous % evaluation. % Note that two parameters can be part of the model and can be % changed without requiring redoing the geom/MC % work, these are: % - an overall spatial scaling factor, xyz_scale, % and an overall Doppler factor, vel_scale. % These might usefully be included in all models. % % Save these variables for re-use of the Geom/MC simulation: static variable GeomMC_Xs, GeomMC_Ys, GeomMC_Enfs; % % Signal that the Geometry and MC calcs have been calculated for % the present parameter set (=1); set to 0 to force re-calculation: static variable GeomMC_Calculated = 0; % % Identify and save the last parameter values used for parameters % whose change will require a re-calculation: static variable Last_sph_rin, Last_sph_rout, Last_cyl_rin, Last_cyl_rout, Last_cyl_len, Last_theta, Last_phi, Last_vel_type, Last_vel_val; % initialize them to unlikely values... Last_sph_rin=-0.5; Last_sph_rout=-0.5; Last_cyl_rin=-0.5; Last_cyl_rout=-0.5; Last_cyl_len=-0.5; Last_theta=-0.5; Last_phi=-0.5; Last_vel_type=-0.5; Last_vel_val=-0.5; % Show the 3d model when it is evaluated ? public variable I3D_SHOW_MODEL = 1; % % In anycase, make the 3d geometry array available for viewing: public variable I3D_GEOM_ARRAY; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% define sphxcyl_line_fit (lo, hi, pars) { % Return the flux vs dispersed (observed) wavelength % for an object defined by the intersection of a % sphere with a cylinder and emitting in a single line. % Both the spatial (projected) and spectral (Doppler) % effects on the line shape are included. % Assume that the full cross-dispersion range is included. % Assume the source size is not too large so that the usual % arf can be applied (after this "flux" calculation): want the % d_arf of d_lambda-across-the-source to be small. % The rmf is applied as usual to provide the point-source % spectral bluring. % Isis_Active_Dataset gives the dataset being evaluated % and can be used to find the grating, order, etc. which % effects the spatial-spectral mapping. % get_data_info() returns structure with useful tags: % order, part (1,2,3 = H,M,L) % % Other data specific to the spectrum are provided by using % get[set]_dataset_metadata (hist_index [, meta_struct]); % and for this model these data should include field(s): % .roll = roll_nom % This can be set at the time of reading in the data via % something like: % % Set the required meta data for hist index 1 : % variable data_meta = struct { roll }; % data_meta.roll = 37.123; % ROLL_NOM in degrees % set_dataset_metadata ( 1, data_meta); % Note the default roll, if no metadata is available, is 0.0. % variable norm, lam_center, vel_scale, xyz_scale; variable zs, flux_out; variable this_info, this_meta, roll_nom; % parameters that can be adjusted w/o re-evaluating Geom/MC: norm = pars[0]; lam_center = pars[1]; vel_scale = pars[2]; xyz_scale = pars[3]; % parameters that will require re-doing the Geom/MC: % Check if any of them have changed: variable par_vals = [ pars[4], pars[5], pars[6], pars[7], pars[8], pars[9], pars[10], int(0.01+pars[11]), pars[12] ]; variable last_vals = [ Last_sph_rin, Last_sph_rout, Last_cyl_rin, Last_cyl_rout, Last_cyl_len, Last_theta, Last_phi, Last_vel_type, Last_vel_val]; if ( sum( int( par_vals != last_vals ) ) != 0 ) { % one or more Geom/MC params have changed - % save these values and flag recalculation: Last_sph_rin = pars[4]; Last_sph_rout = pars[5]; Last_cyl_rin = pars[6]; Last_cyl_rout = pars[7]; Last_cyl_len = pars[8]; Last_theta = pars[9]; Last_phi = pars[10]; Last_vel_type = int(0.01+pars[11]); Last_vel_val = pars[12]; % GeomMC_Calculated = 0; } % OK, here's the Geometry/MC calculation stuff... % to generate the v3d volume for this model % and from that a set of MC points and their "energy factors" % (Doppler factors) % if ( GeomMC_Calculated == 0) { % any other Geom/MC parameters hard wired for now? variable mc_n_points, mc_v3d_res; mc_n_points = 50000L; mc_v3d_res = 37; % The spatial unit for the MC geometry is "arc seconds" % % setup the 3d cube - some factor times the outer size: v3d_setup(max([Last_sph_rout,Last_cyl_rout])*1.20, mc_v3d_res); % Sphere I3D_GEOM_ARRAY = v3d_sphere(Last_sph_rin, Last_sph_rout, [0.,0.,0.]); % Cylinder I3D_GEOM_ARRAY = I3D_GEOM_ARRAY * v3d_cylinder (Last_cyl_rin, Last_cyl_rout, -0.5*Last_cyl_len, 0.5*Last_cyl_len, [0.,0.,0.], [Last_theta, Last_phi]); % Show the model when calculated if desired :-) if ( I3D_SHOW_MODEL == 1 ) { % Show the 3D solid-ish view of the object v3d_view( I3D_GEOM_ARRAY ); % Could be useful to slice off the W-side of the object to show % the inside... %% v3d_view( I3D_GEOM_ARRAY*v3d_cube(v3d_cube_radius*[0.6,1.2,1.2], %% v3d_cube_radius*[-0.6,0.,0.]) ); % Show the projection of emission on the sky: v3d_project( I3D_GEOM_ARRAY ); } % Generate Monte Carlo points (GeomMC_Xs, GeomMC_Ys, zs) = v3d_mc_points(I3D_GEOM_ARRAY, mc_n_points); % Assign Doppler energy factors to the MC points depending on the % velocity type and value GeomMC_Enfs = v3d_doppler_points(GeomMC_Xs, GeomMC_Ys, zs, [0.,0.,0.], Last_vel_val, Last_vel_type, [0.,0.,0.], [Last_theta, Last_phi]); % OK, let's not do all that again... unless we need to. GeomMC_Calculated = 1; } % At this point we have photon locations in sky coord.s and Doppler factors: % GeomMC_Xs, GeomMC_Ys and GeomMC_Enfs --- bascially equivalent to the % input to a MARX simulation (for an otherwise monochromatic source). % The next part of the code essentially does the equivalent % of a MARX simulation, a marx2fits, and a tg_extract (or fil_extract). % i.e., instrument/observation simulation. % - - - - - Instrument simulation - - - % the output dispersed counts histogram variable pseudo_lams; % Angstrom per arcsecond for the parts=1(HEG),2(MEG),3(LEG) % for Chandra with ACIS-S variable A_per_arcsec=[0.0,0.0055671274, 0.011135172, 0.027566069]/0.492; % Grating dispersion angle offsets: roll_nom + Grat_angle gives disp roll % angle in sky coord.s (in degrees), (from fil_get_data.pro) variable Grat_angle=[0.0, -5.235, 4.725, 0.010-0.069]; % OK, now convert the xs, ys and enfs into pseudo-wavelengths % as if they were from a dispersed image. % The grating type, order, and roll matter so find those out: % Get info on this data set: this_info = get_data_info(Isis_Active_Dataset); % Use the useful meta data for the roll: this_meta = get_dataset_metadata(Isis_Active_Dataset); if (this_meta == NULL) { roll_nom=0.0; message(" sphxcyl_line : Roll meta-data not provided for the dataset "+ string(Isis_Active_Dataset) ); message(" *** Will use roll = 0.0 degrees. "); message(" See comments in i3d_sphxcyl.sl for more information."); } else { roll_nom = this_meta.roll; } % Spectral component = Scaled MC wavelengths % Scale the doppler factors by vel_scale and apply to lam_center wavelength: % (a little messy because enfs is 1/(1+z) ) pseudo_lams = lam_center * ( 1.0 + (vel_scale*(1.0/GeomMC_Enfs -1.0)) ); % Spatial component % Now include the spatial location effect on the observed wavelength % for each MC point. The mapping of spatial location to wavelength % depends on the grating order, the A_per_arcsec appropriate for the % grating (part), and the roll. % % Calculate for each MC point its coordinates in the dispsersion-crossdisp % system using the roll_nom and grating angles. % For roll_nom = 0, the +dispersion axis is in the +W direction which % is also +MC_Xs. As roll_nom increases, the +disp vetor moves toward % -S (= -MC_Ys). variable ang_rads, disp_hat, cdisp_hat, along_disp, cross_disp, spatial_lams; %% variable DTOR = (PI/180.0); % assume we have this from v3d.sl % Unit vectors in disp and cdisp in sky = MC coord.s. % Add grating angle and 0.030 deg (for Sky to ACIS angle offset.) ang_rads = DTOR*(roll_nom + 0.030 + Grat_angle[this_info.part]); disp_hat = [cos(ang_rads), -1.0*sin(ang_rads)]; % not sure which sign to use for cross-disp. Start with this... % (cdisp is +N for roll_nom=0; is +E for roll_nom=270) cdisp_hat = [sin(ang_rads), cos(ang_rads)]; % and the distance along the + dispersion direction includeing scale effect: along_disp = xyz_scale * (disp_hat[0]*GeomMC_Xs + disp_hat[1]*GeomMC_Ys); spatial_lams = this_info.order * A_per_arcsec[this_info.part] * along_disp; pseudo_lams = pseudo_lams + spatial_lams; % use histogram routine to bin the lams flux_out = hist1d (pseudo_lams, lo); % normalize it, et voila. flux_out = norm * flux_out/sum(flux_out); % - - - - - % % For reference: % lo, hi = positive histogram grid, with non-overlapping bins, % in ascending order, in Angstrom units. % return value = array of photons/cm^2/s/bin return flux_out; } % Now add this function: add_slang_function ("sphxcyl_line", ["norm [ph/cm^2/s]", "lam_center [A]", "vel_scale [factor]", "xyz_scale [factor]", "sph_rin [arcsec]", "sph_rout [arcsec]", "cyl_rin [arcsec]", "cyl_rout [arcsec]", "cyl_len [arcsec]", "theta [deg,W-to-Away]", "phi [deg,Equat-to-Npole]", "vel_type [0-r^; 1-r; 2-rot, 3-Kep]", "vel_val [v; v/r; v/r; v*sqrt(r)]"], [0]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Set some defaults for its parameters too % define sphxcyl_line_pdfs (i) { variable values=[0.001, 12.132, 1.0, 1.0, 7.0, 10.0, 0.0, 11.0, 6.0, 270.0, -45.0, 2, 2000.0/10.0]; variable freezes=[0, 0, 1, 1, 1,1,1,1,1, 1,1,1,1]; variable pmins=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -360.0, -90.0, 0.0, 0.0]; variable pmaxs=[1.e6, 150.0, 20.0, 20.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 360.0, 90.0, 3.5, 1.e5]; return (values[i], freezes[i], pmins[i], pmaxs[i]); } set_param_default_hook ("sphxcyl_line", "sphxcyl_line_pdfs"); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % and we're done.