%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % v3d.sl version: 1.2.1 , 2011-05-03 % % This file is the Volumetric 3D geometry creation library % Copyright (C) 2006-2010 Massachusetss Institute of Technology % Author: Dan Dewey, dd@space.mit.edu % % Support for this work was provided by the National Aeronautics and % Space Administration through the Smithsonian Astrophysical Observatory % contract SV3-73016 to MIT for Support of the Chandra X-Ray Center, % which is operated by the Smithsonian Astrophysical Observatory for and % on behalf of the National Aeronautics Space Administration under contract % NAS8-03060. % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % - - - "Compiles with any S-Lang" - - - % % Any dependencies beyond just slsh-provided infrastructure % are indicated by "#ifexists needed_routine" statements % in the code below and their consequences are summarized here: % % interp_linear : required for v3d_sphere_rlup From: slgsl, gsl % useful for v3d_cyl_azlup % urand OR : required for the MC routines. From: isis OR % ran_flat slgsl % imdisplay : required for v3d_project. From: slgkt, gtk % required for v3d_evtimg. % required for v3d_evtimg % hist2d : required for v3d_evtimg, 'histxy From: histogram module % volview : required for v3d_view. From: ? (volpack) % h5_write : optional for feature of v3d_view. From: slh5, hdf5 % % required: the routine is not available with out it. % useful : some major functionality of the routine is lost. % optional: some secondary feature of the routine is lost. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Public Functions in this file. % General: % % v3d_setup : Sets size and scale of the v3d cube space. % v3d_help : Gives version info and some help. % v3d_list : Lists the few v3d parameters and their values. % Shapes: : Fills a v3d array with values of 1.0 inside to 0.0 % outside of the volumetric shape. % Each shape has its own parameters and % similar optional translation and orientation values: % oxyz = [X,Y,Z] ; rthphi = [theta, phi] (degrees). % % v3d_cone : v3d_cone (angin, angout, axmin, axmax [,oxyz [, rthphi]]) % The half-angles are in degrees. % v3d_cube : v3d_cube (radii, [,oxyz]) % The cube is aligned with the coordinate axes. % v3d_cylinder : v3d_cylinder (rin, rout, axmin, axmax [,oxyz [, rthphi]]) % [v3d_disk] : --> v3d_cylinder % [v3d_jet] : --> v3d_cone % v3d_roche : v3d_roche (r_sep, q_ratio [, thresh_adj [,oxyz [, rthphi]]]) % r_sep = separation between objects % q_ratio = m_lobe_object / M_companion % thresh_adj adjusts the roche threshold, 2.0 default; % smaller values (0, or -3.0) include more of the lobe. % v3d_sphere : v3d_sphere (rin, rout, [,oxyz ]) % v3d_sphere_line : v3d_sphere_line (p1xyz, p2xyz [, rout [, nspheres]]) % pNxyz = [px,py,pz] % v3d_sphere_ring : v3d_sphere_ring (r_ring, r_sph, nspheres [,oxyz [, rthphi]]) % v3d_spheroid : v3d_spheroid (rin, rout, polar_ratio [,oxyz [, rthphi]]) % v3d_torus : v3d_torus (rin, rout [,oxyz [, rthphi]]) % The radii are the inner-most and outer-most % extent of the torus in the equatorial plane. % Advanced 3D shapes % % v3d_sphere_rlup : v3d_sphere (rin, rout, rlups, vlups [,oxyz ]) % Fill the 3D space from a value-of-radius lookup table. % v3d_cyl_azlup : v3d_cyl_azlup (rin, rout, axmin, axmax, oxyz, rthphi, azlups [, vlups]); % Fill cylinder with azimuthal lookup; azlups are in degrees. % v3d_2dto3d : Maps 2d array of quadrant values into the 3d array. % v3d_histxy : Create a histogram in the z=0 plane of x,y input values. % v3d_reproject : Maps an array to a new shape retaining the x,y projection. % % Functions in 3D - fills the v3d array with continuous values... % % v3d_r3dsq : The array cells contain the r_3d^2 from the center. % v3d_f_generic : User provided function of available coordinates. % Copy/edit/rename v3d_f_generic.sl . % Visualization % % v3d_view : v3d_view (v3darr) - wrapper for volview % v3d_project : v3d_project (v3darr) - wrapper for imdisplay % For both of these, values less than 0 are treated % as 0. % v3d_evtimg : v3d_evtimg (xs, ys [, es [, Rrange, Grange, Brange]]); % Rrange = [rmin,rmax] selects on es for the R plane. % % Monte-Carlo % % v3d_mc_cells : (ixs,iys,izs) = v3d_mc_cells (v3darr, nevts); % Return random cells with prob based on v3darr values. % v3d_mc_points : (xs,ys,zs) = v3d_mc_points (v3darr, nevts); % Return random points w/prob based on v3darr values. % v3d_doppler_points : en_factor = v3d_doppler_points (x,y,z, vsxyz [, % velval [,vtype [, oxyz [,rthphi]]]]) % Calculate the line-of-sight velocity for the given % points and velocity specification. % vsxyz=[vsx,vsy,vsz] in km/s; % velval in k/ms [per unit]; % vtype = -1 [no additional v], % 0 [velval*r_hat], % 1 [velval*r_vector], % 2 [rigid body rot about rthphi-through-oxyz] % 3 [Keplerian v(r) about rthphi-through-oxyz] % X-Ray application % % (v3d_spect_... % (v3d_make_events... % (v3d_marx_model... % Ancillary: % % v3d_soft_cutoff : a soft step-down function % v3d_merge_max : a la the IDL ">" operator: c = ( a > b ) % v3d_smooth_box : return the box-car smoothed (averaged) array values, 1D. % v3d_max_box : return the box-car maximum array values, 1D. % v3d_smooth_2d : return the box-car smoothed (averaged) array values, 2D. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Public Variables in this file % v3d_na : the number of cells _per_ _radius_ along each axis. % The _total_ cells in, e.g., the x direction is: % 2*v3d_na+1 % and so to loop over all cell x-indices we'd do: % _for ix (0, 2*v3d_na, 1) { ... } % v3d_va : the spatial coordinate value along an axis, % specified at the _center_ of the cell. % v3d_cube_radius : the value in user units of the center of the % extreme cells. % v3d_user_unit : a string with the units, default "user_unit" % v3d_version_str : a string giving the v3d version. % v3d_version_date : a string giving the date of the v3d version. % V3D_HARD_CUTOFF : set to 0(soft) or 1(hard) to control operation of v3d_soft_cutoff. % save_image_file : default is NULL for screen viewing; set to a filename % string to save the v3d_project or v3d_view to a file. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Other functions defined for use in this file. % urand, seed_random % Created from gsl if urand does not already exist. % % clip_max - clip an array elements at a given value. % % Other variables defined for use in this file. variable DTOR = (PI/180.0); % Physical constants % % If the nominal ISIS-provided physical constants don't % exist, then provide them: #ifnexists Const_c % This is the contents of $ISIS_SRCDIR/share/physconst.sl % These values are in cgs system % 1998 CODATA recommended values % The following are exact public variable Const_c = 2.99792458e10; % [cm/s] speed of light % public variable Const_h = 6.62606876e-27; % [erg s] Planck const public variable Const_G = 6.673e-8; % [cm^3/g/s^2] Newtonian const public variable Const_m_e = 9.10938188e-28; % [g] electron mass public variable Const_m_p = 1.67262158e-24; % [g] proton mass public variable Const_k = 1.3806503e-16; % [erg/K] Boltzmann const % Conversion factors public variable Const_eV = 1.602176462e-12; % [erg] public variable Const_u = 1.66053873e-24; % [g] atomic mass unit % hc in [keV Angstrom] public variable Const_keV_A = Const_h * Const_c * 1.e8 / (1.e3 * Const_eV); public variable Const_hc = Const_keV_A; #endif % Additional constants: #ifnexists Const_pc % - - - Additional physical constants public variable Const_pc = 3.0857e18; % cm public variable Const_yr = 365.24 * 24.0 * 3600.0; % s / year public variable Const_SM = 1.98892e33; % g / solar mass public variable Const_amu = 1.66054e-24; % g / amu public variable Const_c_kms = Const_c * 1.e-5; #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % OK, let's do it... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public variable v3d_version_str = "1.2.1"; public variable v3d_version_date = "2011-05-03"; public variable v3d_na = 37; public variable v3d_va = Float_Type [2*v3d_na+1]; public variable v3d_cube_radius = 100.0; public variable v3d_user_unit = "user_unit"; public variable save_image_file = NULL; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_help () % {{{ %!%+{{{ %\function{v3d_help} %\synopsis{List a summary of the v3d routines} %\usage{ v3d_help; } %\description % List a summary of the v3d routines. %\notes % %\seealso{v3d_setup} %!%- }}} { message(""); message(" v3d.sl version: "+v3d_version_str+" , "+v3d_version_date); message(""); message(" Routines available:"); message(" General: v3d_ 'setup, 'help, 'list "); message(" Shapes: v3d_ 'cube, 'cone, 'cylinder, 'roche, 'sphere, 'spheroid, 'torus"); message(" ancillary: v3d_ 'soft_cutoff, 'merge_max, 'smooth_box, 'max_box, 'smooth_2d"); message(" Composite: v3d_sphere_line, v3d_sphere_ring"); message(" Advanced: v3d_sphere_rlup, '_cyl_azlup, '_2dto3d, '_histxy, 'reproject"); message(" Functions: v3d_r3dsq, v3d_f_generic"); message(" Display: v3d_view, v3d_project, save_image_file = NULL or "); message(" x,y[,e] points: v3d_evtimg"); message(" Monte-Carlo: v3d_mc_cells, v3d_mc_points,"); message(" v3d_doppler_points"); message(""); message(" see also: http://space.mit.edu/home/dd/v3d.html"); message(""); } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%v3d_help; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_list () % {{{ %!%+{{{ %\function{v3d_list} %\synopsis{List the current values of v3d cube radius and resolution} %\usage{ v3d_list; } %\description % List the values of v3d_cube_radius, v3d_user_unit and % the resolution, v3d_na, the number of cells per radius. %\notes % %\seealso{v3d_setup} %!%- }}} { message(""); message(" v3d_cube_radius = "+string(v3d_cube_radius)+ " ["+v3d_user_unit+"]"); message(" v3d_na = "+string(v3d_na)+ " cells per radius."); message(" cell-cell spacing = "+string(v3d_cube_radius/v3d_na)+ " ["+v3d_user_unit+"]"); message(" size of 3D cubes = "+string(2*v3d_na+1)+"^3 = " +string((2*v3d_na+1)^3)+" cells." ); message(""); } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_setup () % {{{ %!%+{{{ %\function{v3d_setup} %\synopsis{Initialize the v3d size and resolution values} %\usage{ v3d_setup (cube_radius, n_cells_per_radius); } %\description % The variables v3d_cube_radius and v3d_na are set and used % by other v3d routines. % Also, the variable v3d_va is filled with the 2N+1 values % of the centers of the cells along the (symmetric) axes. %\notes % The variable v3d_user_unit is available to store a desired % string describing the units of v3d_cube_radius, but it is not % required to be used. %\seealso{v3d_list} %!%- }}} { variable radius, numper; if (_NARGS != 2) usage("................. v3d_setup (radius, numper);"); numper = (); radius = (); % put a default for any NULL cases if (numper == NULL) numper = 37; if (radius == NULL) radius = 100.0; v3d_na = numper; v3d_cube_radius = radius; v3d_va = (v3d_cube_radius/double(numper))*[-1*v3d_na : v3d_na : 1]; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v3d_setup(,); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public variable V3D_HARD_CUTOFF = 0; public define v3d_soft_cutoff () % {{{ %!%+{{{ %\function{v3d_soft_cutoff} %\synopsis{Provides a continuous (smoothed) step function} %\usage{ where_negative = soft_cutoff (values - threshold); } %\description % A function that returns: % 1.0 for values << 0 % 0.0 for values >> 0 % and transitions from 1.0 to 0.0 in the % range of values -5 to +5 with % -0.5 to 0.5 going from 0.7 to 0.3 ish % -1 to +1 going from 90% to 0.1 roughly % and % -2 to +2 going from 99% to <0.01 % % Use v3d_soft_cutoff() to define slightly smooth edges: % % ... v3d_soft_cutoff( (location-boundary) / % (v3d_cube_radius/(2.0*v3d_na)) ) %\notes % The public variable V3D_HARD_CUTOFF can be set to 1 to % cause 'soft_cutoff to return only 0.0 and 1.0 values. % % It's possible to see the details of this function % by showing how much it differs from 0 or +1 : % v3d_setup(10.0,37); % plot(v3d_va, 0.5 - abs( v3d_soft_cutoff(v3d_va) - 0.5) ); %\seealso{ } %!%- }}} { if (_NARGS != 1) usage("................. v3d_soft_cutoff (values)"); variable values; values = (); variable fracs, large, trans; % Default 0.0 , for values >> 0 fracs = values*0.0; % Do it "hard" or "soft" ? if ( V3D_HARD_CUTOFF == 1) { % Negative values have fracs set to 1.0 large = where(values < 0.0); if ( length(large) > 0 ) fracs[large] = 1.0; } else { % Large negative values have fracs set to 1.0 large = where(values < -5.0); if ( length(large) > 0 ) fracs[large] = 1.0; % Transition values in -5 to 5 range trans = where(abs(values) < 5.0); if (length(trans) > 0) fracs[trans] = exp(-((values[trans]+9.55)/10.0)^12); } return fracs; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_merge_max () % {{{ %!%+{{{ %\function{v3d_merge_max} %\synopsis{Combine two 3D arrays by taking the cell-by-cell maximum} %\usage{ arr3d_max = v3d_merge_max (arr3d_a, arr3d_b); } %\description % Combines two 3D arrays by returning an array containing the % cell-by-cell maximum values. This is useful to produce the % union of two geomtric shapes. %\notes % %\seealso{ } %!%- }}} { % returns an array that contains the max of either input % arrays. if (_NARGS != 2) usage("................. v3d_merge_max (a, b);"); variable a, b; b = (); a = (); variable addins; addins = where(b > a); if ( length(addins) > 0 ) a[addins] = b[addins]; return a; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_smooth_box() % {{{ %!%+{{{ %\function{v3d_smooth_box} %\synopsis{Do box-car smoothing of a 1D array} %\usage{ smarr1d = v3d_smooth_box (arr1d, nsmooth); } %\description % Returns a smoothed version of the 1d array % smoothing is in integer units nsmooth = (1,) 2, 3,... % - for odd integers the average of nsmooth values is used % - for even integers, the next largest odd value is used % and the two extreme bins are weighted by 1/2. % %\notes % The ISIS routine "rebin" does a similar (better) % thing when the output bin scheme is well defined. % %\seealso{v3d_smooth_2d, v3d_2dto3d} %!%- }}} { if (_NARGS != 2) usage("................. v3d_smooth_box (arr1d, nsmooth);"); variable arr1d, nsmooth; nsmooth = (); arr1d = (); variable arrout, halfn, evenn, kernel, ia; % make sure the input is an integer nsmooth = int(nsmooth); % and catch the simple case... if ( nsmooth < 2 ) return arr1d; % define arrout arrout = Float_Type [length(arr1d)]; % decide on the size of the box = 2 * halfn + 1 % and set up the kernel halfn = int(nsmooth)/2; kernel = 1.0 + Float_Type [2*halfn+1]; % the extereme values are 1/2 if nsmooth is even: if ( 2*halfn == nsmooth ) { kernel[0] = 0.5; kernel[2*halfn] = 0.5; } % pad arr1d for the kernel length: arr1d = [arr1d[0]+Float_Type[halfn], arr1d, arr1d[length(arr1d)-1]+Float_Type[halfn]]; % and do it... _for ia (0, length(arrout)-1, 1) { arrout[ia] = sum( kernel * arr1d[[ia:ia+2*halfn]] ); } return arrout/sum(kernel); } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_max_box() % {{{ %!%+{{{ %\function{v3d_max_box} %\synopsis{Do box-car maximum of a 1D array} %\usage{ maxarr1d = v3d_max_box (arr1d, nwidth); } %\description % Returns the running maximum of the 1d array % - for odd integers the max of nsmooth values is used % - for even integers, the next largest odd value is used % and the two extreme bins are weighted by 1/2. % %\notes % %\seealso{} %!%- }}} { if (_NARGS != 2) usage("................. s3d_max_box (arr1d, nsmooth);"); variable arr1d, nsmooth; nsmooth = (); arr1d = (); variable arrout, halfn, evenn, kernel, ia; % make sure the input is an integer nsmooth = int(nsmooth); % and catch the simple case... if ( nsmooth < 2 ) return arr1d; % define arrout arrout = Float_Type [length(arr1d)]; % decide on the size of the box = 2 * halfn + 1 % and set up the kernel halfn = int(nsmooth)/2; kernel = 1.0 + Float_Type [2*halfn+1]; % the extereme values are 1/2 if nsmooth is even: if ( 2*halfn == nsmooth ) { kernel[0] = 0.5; kernel[2*halfn] = 0.5; } % pad arr1d for the kernel length: arr1d = [arr1d[0]+Float_Type[halfn], arr1d, arr1d[length(arr1d)-1]+Float_Type[halfn]]; % and do it... _for ia (0, length(arrout)-1, 1) { arrout[ia] = max( kernel * arr1d[[ia:ia+2*halfn]] ); } return arrout; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_smooth_2d () % {{{ %!%+{{{ %\function{v3d_smooth_2d} %\synopsis{Box-car smooth a 2D array} %\usage{ smarr2d = v3d_smooth_2d (arr2d, nsmooth); } %\description % Simple box-car smoothing in 2D. %\notes % Used by the simple v3d_evtimg routine. %\seealso{v3d_smooth_box} %!%- }}} { if (_NARGS != 2) usage("................. v3d_smooth_2d (arr2d, nsmooth);"); variable arr2d, nsmooth; nsmooth = (); arr2d = (); % catch the simple case if (nsmooth < 2) return arr2d; variable ia, arrout=0.0*arr2d; % Smooth in the two dimensions... _for ia (0,length(arr2d[*,0])-1,1) { arrout[ia,*] = v3d_smooth_box(arr2d[ia,*],nsmooth); } _for ia (0,length(arr2d[0,*])-1,1) { arrout[*,ia] = v3d_smooth_box(arrout[*,ia],nsmooth); } return arrout; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_r3dsq () % {{{ %!%+{{{ %\function{v3d_r3dsq} %\synopsis{Returns a 3d array with r-squared in each cell} %\usage{ v3darr = v3d_r3dsq (oxyz=[x,y,z]); } %\description % Given a center location, oxyz=[x,y,z], this % returns a v3d array filled with each cell's % distance squared, r3dsq, from the center location. %\notes % %\seealso{ } %!%- }}} { if (_NARGS != 1) usage("................. v3d_r3dsq (oxyz=[x,y,z]);"); variable oxyz = (); % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable ia; % Fill the array with the value of r^2 by % summing (x-xoff)^2 (y-yoff)^2 (z-zoff)^2 in turn: _for ia (0, 2*v3d_na, 1) { val[ia,*,*] = (v3d_va[ia]-oxyz[0])^2 ; } _for ia (0, 2*v3d_na, 1) { val[*,ia,*] += (v3d_va[ia]-oxyz[1])^2 ; } _for ia (0, 2*v3d_na, 1) { val[*,*,ia] += (v3d_va[ia]-oxyz[2])^2 ; } return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_cube () % {{{ %!%+{{{ %\function{v3d_cube} %\synopsis{Fill a v3darr with 1.0 within a cube} %\usage{ v3darr = v3d_cube (radii [,oxyz]); } %\description % Fill a 3d array with values of 1.0 inside a cube % of radius radii centered at OXYZ. % % Radius can be a single value (the same for all axes) % or an array of 3 values, one for each axis: % radii = radius or [radx, rady, radz] % %\notes % Note that this cube is always aligned with the axes % and hence takes no optional rthphi parameter. % %\seealso{ } %!%- }}} { % % RFE - v3d_cube - Add rthphi parameter to orient the cube. % variable radii, oxyz; switch (_NARGS) { case 1: oxyz = [0.,0.,0.]; radii = (); } { case 2: oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; radii = (); } { % otherwise: % give extra info before the usage: message(" radii = radius or [radx, rady, radz]"); usage("................. v3d_cube (radii, [,oxyz]);"); } % Defaults for NULL argument % and make radii an array of length 3 if not already % % don't want to encourage this, but could be useful: if (radii == NULL) radii = v3d_cube_radius*[1.0,1.0,1.0]; % is it length 1 ? if ( length(radii) == 1 ) radii = radii * [1.0,1.0,1.0]; % The 3-D volume array: % set to all 1's variable val = 1.0 + Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % nominal variables... variable afact, ia; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); % go through each axis and multiply by the factor for that axis... % x-axis afact = v3d_soft_cutoff( ( abs(v3d_va - oxyz[0]) - radii[0] ) / cutoff_scale ) ; _for ia ( 0, 2*v3d_na, 1) { val[ia,*,*] = val[ia,*,*] * afact[ia] ; } % y-axis afact = v3d_soft_cutoff( ( abs(v3d_va - oxyz[1]) - radii[1] ) / cutoff_scale ) ; _for ia ( 0, 2*v3d_na, 1) { val[*,ia,*] = val[*,ia,*] * afact[ia] ; } % z-axis afact = v3d_soft_cutoff( ( abs(v3d_va - oxyz[2]) - radii[2] ) / cutoff_scale ) ; _for ia ( 0, 2*v3d_na, 1) { val[*,*,ia] = val[*,*,ia] * afact[ia] ; } return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_cone () % {{{ %!%+{{{ %\function{v3d_cone} %\synopsis{Fill a v3darr with 1.0 within a cone}} %\usage{ v3darr = v3d_cone (angin, angout, axmin, axmax [,oxyz [, rthphi]]); } %\description % Creates a cone based on the parameters: % angin, angout: opening half-angles of inner and outer boundaries % in degrees % axmin, axmax: axial limits (w.r.t. center) of the emission % Centered (vertex) at OXYZ % The axis of the cone is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Y. % %\notes % %\seealso{ } %!%- }}} { variable angin, angout, axmin, axmax, oxyz, rthphi; variable def_rthphi = [0.0, 90.0]; switch (_NARGS) { case 4: rthphi = def_rthphi; oxyz = [0.,0.,0.]; axmax = (); axmin = (); angout = (); angin = (); } { case 5: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); angout = (); angin = (); } { case 6: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); angout = (); angin = (); } { % otherwise: message(" angin, angout are in degrees"); usage("................. v3d_cone (angin, angout, axmin, axmax [,oxyz [, rthphi]]);"); } % Defaults for NULL arguments: if (angin == NULL) angin = 0.0; if (angout == NULL) angout = 89.999; if (axmin == NULL) axmin = -4.0*v3d_cube_radius; if (axmax == NULL) axmax = 4.0*v3d_cube_radius; % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Define some useful constant values % compute the tangents of the angles: variable tanangin = tan(angin*DTOR); variable tanangout = tan(angout*DTOR); % nominal variables... variable axial, radial; variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; variable maxaxial = max([abs(axmin),abs(axmax)]); variable radial2limit = ((tanangout*maxaxial)+5.0*cutoff_scale)^2; variable rad2limit = radial2limit + (maxaxial+5.0*cutoff_scale)^2; % Loop over the 'toward-away' 'front-back' (z) value % Do a simple check to see if this whole plane is worth processing _for iz ( 0, 2*v3d_na, 1) { if ( (v3d_va[iz]-oxyz[2])^2 < rad2limit ) { % Loop over the 'vertical (up, North)' (y) value % and another quick check... _for iy ( 0, 2*v3d_na, 1) { if ( ( (v3d_va[iz]-oxyz[2])^2 + (v3d_va[iy]-oxyz[1])^2 ) < rad2limit) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin: % roughly: % axial = (points - origin) dotted with axisunitvect % radial = SQRT( dist_squared(points - origin) - axial^2 ) % more specifically: % axial = (px-ox)*axuvx + (py-oy)*axuvy + (pz-oz)*axuvz % radial = SQRT( (px-ox)^2 + (py-oy)^2 + (pz-oz)^2 - axial^2 ) % finally in code with "px" as an array: axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; % Only do this if some of the radial values are less than radial2limit if ( min(radial) < radial2limit ) { radial = sqrt(radial) ; % Now use the axial and radial values to define the desired shape. % Cut it off in axial dimensions axfact = v3d_soft_cutoff( (axial - axmax) / cutoff_scale ) * v3d_soft_cutoff( (axmin - axial) / cutoff_scale ) ; % Fill val with uniform value (1.0) inside the cone % val[*,iy,iz] = axfact * v3d_soft_cutoff( (radial - tanangout*abs(axial)) / cutoff_scale ) * v3d_soft_cutoff( (tanangin*abs(axial) - radial) / cutoff_scale ) ; } % if (min()) } % if ... } % y loop } % if ... } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_cylinder () % {{{ %!%+{{{ %\function{v3d_cylinder} %\synopsis{Fill a v3darr with 1.0 within a cylinder} %\usage{ v3darr = v3d_cylinder (rin, rout, axmin, axmax [,oxyz [, rthphi]]); } %\description % % Creates a cylinder using the parameters: % rin, rout: radii of inner and outer boundaries % axmin, axmax: axial limits (w.r.t. center) of the emission % Centered at OXYZ % Axis of the cylinder is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Z. % %\notes % %\seealso{v3d_cyl_azlup} %!%- }}} { variable rin, rout, axmin, axmax, oxyz, rthphi; variable def_rthphi = [270.0, 0.0]; switch (_NARGS) { case 4: rthphi = def_rthphi; oxyz = [0.,0.,0.]; axmax = (); axmin = (); rout = (); rin = (); } { case 5: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); rout = (); rin = (); } { case 6: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); rout = (); rin = (); } { % otherwise: usage("................. v3d_cylinder (rin, rout, axmin, axmax [,oxyz [, rthphi]]);"); } % NO defaults for NULL primary arguments: if ((rin == NULL) or (rout == NULL) or (axmin == NULL) or (axmax == NULL)) { usage("................. v3d_cylinder (rin, rout, axmin, axmax [,oxyz [, rthphi]]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % nominal variables... variable axial, radial; variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; variable radial2limit = (rout+5.0*cutoff_scale)^2; variable maxaxial = max([abs(axmin),abs(axmax)]); variable rad2limit = radial2limit + (maxaxial+5.0*cutoff_scale)^2; % Loop over the 'toward-away' 'front-back' (z) value % and only do relevant planes _for iz ( 0, 2*v3d_na, 1) { if ( (v3d_va[iz]-oxyz[2])^2 < rad2limit ) { % Loop over the 'vertical (up, North)' (y) value % and skip out of range ones... _for iy ( 0, 2*v3d_na, 1) { if ( ( (v3d_va[iz]-oxyz[2])^2 + (v3d_va[iy]-oxyz[1])^2 ) < rad2limit ) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; % Only do this if some of the radial values are less than radial2limit if ( min(radial) < radial2limit ) { radial = sqrt(radial) ; % Now use the axial and radial values to define the desired shape. % Cut it off in axial dimensions axfact = v3d_soft_cutoff( (axial - axmax) / cutoff_scale ) * v3d_soft_cutoff( (axmin - axial) / cutoff_scale ) ; % Fill val with uniform value (1.0) inside the cylinder % % If rin is very small or zero then only need the rout cutoff: if ( rin < cutoff_scale ) { val[*,iy,iz] = axfact * v3d_soft_cutoff( (radial - rout) / cutoff_scale ); } else { % cutoff in and out: val[*,iy,iz] = axfact * v3d_soft_cutoff( (radial - rout) / cutoff_scale ) * v3d_soft_cutoff( ( rin - radial) / cutoff_scale ); } } % if (min()...) } % if ... } % y loop } % if ... } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_cyl_azlup () % {{{ %!%+{{{ %\function{v3d_cyl_azlup} %\synopsis{Fill a cylinder using an azimuth look table} %\usage{ v3darr = v3d_cyl_azlup (rin, rout, axmin, axmax, oxyz, rthphi, azlups [, vlups]); } %\description % % Create a cylinder with aximuthal lookup values or limits. % % rin, rout: radii of inner and outer boundaries % axmin, axmax: axial limits (w.r.t. center) of the emission % Centered at OXYZ % Axis of the cylinder is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Z. % % The arguments azlups and vlups give the azimuth (in degrees, % from -180 to 180) and scalar values look-up arrays used % to assign values vs azimuth. % % If vlups is not present then the first and last elements of % azlups are used to define the azimuthal range in which the % cylinder has value 1.0 and 0.0 otherwise. % %\notes % Uses the gsl interp_linear function to do the lookup of vlups. % If interp_linear is not available then the routine will still % compile but functions as if vlups is not provided. %\seealso{ } %!%- }}} { variable rin, rout, axmin, axmax, oxyz, rthphi, azlups, vlups; variable def_rthphi = [270.0, 0.0]; switch (_NARGS) { case 7: vlups = NULL; azlups = (); rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); rout = (); rin = (); if (rin == NULL) rin = 0.0; } { case 8: vlups = (); azlups = (); rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; axmax = (); axmin = (); rout = (); rin = (); if (rin == NULL) rin = 0.0; } { % otherwise: usage("................. v3d_cyl_azlup (rin, rout, axmin, axmax, oxyz, rthphi, azlups [, vlups]);"); } % NO defaults for NULL of these primary arguments: if ( (rout == NULL) or (axmin == NULL) or (axmax == NULL) or (azlups == NULL) ) { usage("................. v3d_cyl_azlup (rin, rout, axmin, axmax, oxyz, rthphi, azlups [, vlups]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Unit vector in the y-prime (azimuth=0) direction variable ypuvx, ypuvy, ypuvz; ypuvx = sin(rthphi[1]*DTOR)*(-1.0)*cos(rthphi[0]*DTOR); ypuvy = cos(rthphi[1]*DTOR); ypuvz = sin(rthphi[1]*DTOR)*sin(rthphi[0]*DTOR); % Unit vector in the z-prime (az = 90 degrees) direction variable zpuvx, zpuvy, zpuvz; zpuvx = sin(rthphi[0]*DTOR); zpuvy = 0.0; zpuvz = cos(rthphi[0]*DTOR); % nominal variables... variable axial, radial, yprime, zprime; variable azimuth; variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; variable radial2limit = (rout+5.0*cutoff_scale)^2; variable maxaxial = max([abs(axmin),abs(axmax)]); variable rad2limit = radial2limit + (maxaxial+5.0*cutoff_scale)^2; % cutoff scale and values for the azimuths variable az_cutoff_scale = 2.0*cutoff_scale / (rin+rout); variable azstart, azstop, azsel; % Azimuths: move them into the -180 to 180 range if needed azsel = where(azlups > 180.0); if (length(azsel) > 0) azlups[azsel] = azlups[azsel] - 360.0; azsel = where(azlups < -180.0); if (length(azsel) > 0) azlups[azsel] = azlups[azsel] + 360.0; % azimuthal limits in radians azstart = DTOR * azlups[0]; azstop = DTOR * azlups[length(azlups)-1]; % want to allow a full 360 degree cylinder without a little gap % so if the start-stop range is very close (0.1 deg.) to 0 or 360 then % skip the azimuth cutoff... variable SKIP_AZ = 0; if ( abs(azstop-azstart) < DTOR*0.1 ) SKIP_AZ = 1; if ( abs(abs(azstop-azstart) - 2.0*PI) < DTOR*0.1 ) SKIP_AZ = 1; % If vlups are provided then prepare the lookups... if ( vlups != NULL ) { % azimuth lookup values: % % add 0.0's on the ends - don't need to do this... %% vlups = [0.0,vlups,0.0]; % and add ~ duplicate az s on the ends too %% azlups = [azlups[0]-0.001, azlups, azlups[length(azlups)-1]+0.001]; % % put the azlups and vlups into azlups increasing order: azsel = array_sort(azlups); azlups = azlups[azsel]; vlups = vlups[azsel]; % and include values outside the ranges: azlups = [ -200.0, azlups, 200.0 ]; vlups = [ vlups[0], vlups, vlups[length(vlups)-1] ]; % and convert to radians: azlups = DTOR * azlups; } % Loop over the 'toward-away' 'front-back' (z) value % and only do relevant planes _for iz ( 0, 2*v3d_na, 1) { if ( (v3d_va[iz]-oxyz[2])^2 < rad2limit ) { % Loop over the 'vertical (up, North)' (y) value % and skip out of range ones... _for iy ( 0, 2*v3d_na, 1) { if ( ( (v3d_va[iz]-oxyz[2])^2 + (v3d_va[iy]-oxyz[1])^2 ) < rad2limit ) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; % radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; % Only do this if some of the radial values are less than radial2limit if ( min(radial) < radial2limit ) { radial = sqrt(radial) ; % Now use the axial and radial values to define the desired shape. % Cut it off in axial dimensions axfact = v3d_soft_cutoff( (axial - axmax) / cutoff_scale ) * v3d_soft_cutoff( (axmin - axial) / cutoff_scale ) ; % Fill val with uniform value (1.0) inside the cylinder % % If rin is very small or zero then only need the rout cutoff: if ( rin < cutoff_scale ) { val[*,iy,iz] = axfact * v3d_soft_cutoff( (radial - rout) / cutoff_scale ); } else { % cutoff in and out: val[*,iy,iz] = axfact * v3d_soft_cutoff( (radial - rout) / cutoff_scale ) * v3d_soft_cutoff( ( rin - radial) / cutoff_scale ); } % Calculate the azimuths and do lookup or limits % y_prime values: yprime = (v3d_va - oxyz[0])*ypuvx + (v3d_va[iy]-oxyz[1])*ypuvy + (v3d_va[iz]-oxyz[2])*ypuvz ; % z_prime values: zprime = (v3d_va - oxyz[0])*zpuvx + (v3d_va[iy]-oxyz[1])*zpuvy + (v3d_va[iz]-oxyz[2])*zpuvz ; % % Azimuth from y-prime toward z-prime: azimuth = atan2(zprime, yprime); % special case: if azstart = azstop then no azimuth cutoff if ( SKIP_AZ == 0 ) { % if azstart < azstop then use the range between them: if ( azstop > azstart ) { % normal selection not including the 180 to -180 region val[*,iy,iz] = val[*,iy,iz] * v3d_soft_cutoff( (azimuth - azstop)/az_cutoff_scale) * v3d_soft_cutoff( (azstart - azimuth)/az_cutoff_scale) ; } else { % this selection includes the 180 to -180 region % so divide it into two pieces that are "OR"ed together val[*,iy,iz] = val[*,iy,iz] * v3d_merge_max( v3d_soft_cutoff( (azimuth - azstop)/az_cutoff_scale) , v3d_soft_cutoff( (azstart - azimuth)/az_cutoff_scale) ); } } % Finally, do the lookup if vlups is given: % (provided interp_linear is available...) #ifexists interp_linear if ( vlups != NULL ) { val[*,iy,iz] = val[*,iy,iz] * interp_linear(azimuth, azlups, vlups); } #endif } % if (min()...) } % if ... } % y loop } % if ... } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_roche () % {{{ %!%+{{{ %\function{v3d_list} %\synopsis{Fill a v3darr with 1.0 within a Roche lobe} %\usage{ v3darr = v3d_roche (r_sep, q_ratio [, thresh_adj, [,oxyz [, rthphi]]]); } %\description % Creates a roche lobe filled with 1.0 inside. % The gravitational potential is used % to determine the boundary of the lobe. % Parameters are: % % r_sep = separation between objects % q_ratio = m_lobe_object / M_companion % % thresh_adj gives the percentage increase in potential applied to set % the roche volume. A value of 2 gives a nice looking shape; % a value of 0 or minus-a-few might better represent where wind % (higher KE, loosely bound) material may be flowing in the system... % % Centered at OXYZ % Direction to companion is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along -X. %\notes % %\seealso{ } %!%- }}} { variable r_sep, q_ratio, thresh_adj, oxyz, rthphi; variable def_rthphi = [180.0, 0.0]; switch (_NARGS) { case 2: rthphi = def_rthphi; oxyz = [0.,0.,0.]; thresh_adj = 2.0; q_ratio = (); r_sep = (); } { case 3: rthphi = def_rthphi; oxyz = [0.,0.,0.]; thresh_adj = (); q_ratio = (); r_sep = (); } { case 4: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; thresh_adj = (); q_ratio = (); r_sep = (); } { case 5: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; thresh_adj = (); q_ratio = (); r_sep = (); } { % otherwise: message(" r_sep = separation between objects"); message(" q_ratio = m_lobe_object / M_companion"); message(" thresh_adj adjusts the roche threshold, 2.0 default,"); message(" smaller values (0, or -3.0) include more of the lobe."); usage("................. v3d_roche (r_sep, q_ratio [, thresh_adj, [,oxyz [, rthphi]]]);"); } % Defaults for NULL primary arguments: if ( thresh_adj == NULL ) thresh_adj = 2.0; % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Some useful variables % Find the location and value of the grav potential at the % inner L1 point - use these to determine inside/outside % of the lobe in the code further below. % Set m1 = 1.0 and use r in units of r_sep. variable m1 = 1.0; variable M2 = 1.0/q_ratio; variable r1 = sqrt(q_ratio)/(1.0+sqrt(q_ratio)); variable r2 = 1.0 - r1; % useful pre-compute: variable r_sepsq = r_sep*r_sep; % negative of the grav potential: variable npot_lobe = m1/r1 + M2/r2; % the threshold to define the lobe is thresh_adj percent % above this level: variable npot_thresh = (1.0+thresh_adj/100.0)*npot_lobe; % set the scale for the potential "cell" size: variable pot_cell_size = 0.05*npot_lobe; % nominal variables... variable axial, radial; variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; % other variables variable radialsep2, axialsep, r1_2d, r2_2d; variable pot_vals, pot_cutoff_scale; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; % * keep radial as squared by not doing: radial = sqrt(radial) ; % Now use the axial and radial values to define the desired shape. % Cut off the roche lobe of the other object (at large axial value) axfact = v3d_soft_cutoff( (axial - r1*r_sep) / cutoff_scale ) ; % Fill val with uniform value (1.0) inside the lobe % normalized coord.s (units of r_sep) % radial is really radial^2 radialsep2 = radial/r_sepsq; axialsep = axial/r_sep; % the 2d distances; keep them a bit above 0.0... r1_2d = sqrt(radialsep2 + axialsep^2 + 1.e-7); r2_2d = sqrt(radialsep2 + (1.0-axialsep)^2 + 1.e-7); % the potential values for these cells: pot_vals = m1/r1_2d + M2/r2_2d; % since we're comparing potentials, the cutoff scale is % in potential units... % to look good: vary the potential cutoff size with axial distance... pot_cutoff_scale = pot_cell_size/(1.0 + abs(axialsep/r1 + 1.0)); val[*,iy,iz] = axfact * v3d_soft_cutoff( (npot_thresh - pot_vals) / pot_cutoff_scale ); } % y loop } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_sphere () % {{{ %!%+{{{ %\function{v3d_sphere} %\synopsis{Fill a v3darray with 1.0 inside a spherical shell} %\usage{ v3darr = v3d_sphere (rin, rout [,oxyz ]); } %\description % Creates a 3D datacube describing a sphere % with radii rin and rout, centered at oxyz=[x,y,z] %\notes % %\seealso{v3d_sphere_rlup} %!%- }}} { variable rin, rout, oxyz; switch (_NARGS) { case 2: oxyz = [0.,0.,0.]; rout = (); rin = (); } { case 3: oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; rout = (); rin = (); } { % otherwise: usage("................. v3d_sphere (rin, rout, [,oxyz ]);"); } % defaults for NULL primary arguments: if ( rin == NULL ) rin = 0.0; % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the system is at: % oxyz % Variable for the r^2 values variable rsquared, rzsquared; % nominal variables... variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable rad2limit = (rout+5.0*cutoff_scale)^2; variable iz, iy; variable these_rs; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { % start rsquared: rzsquared = (v3d_va[iz]-oxyz[2])^2; % only continue this if rzsquared is small enough if ( rzsquared < rad2limit ) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { % add in the y component to rsquared rsquared = rzsquared + (v3d_va[iy]-oxyz[1])^2; % only do this if rsquared is small enough if ( rsquared < rad2limit ) { % OK, evaluate these... % the radii-squared values of these cells: these_rs = sqrt( rsquared + (v3d_va-oxyz[0])^2 ); % If rin is very small or zero then only need the rout cutoff: if ( rin < cutoff_scale ) { val[*,iy,iz] = v3d_soft_cutoff( (these_rs - rout) / cutoff_scale ) ; } else { % cutoff in and out: val[*,iy,iz] = v3d_soft_cutoff( (these_rs - rout) / cutoff_scale ) * v3d_soft_cutoff( (rin - these_rs) / cutoff_scale ) ; } } % if rsquared } % y loop } % if rzsquared } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #ifexists interp_linear public define v3d_sphere_rlup () % {{{ %!%+{{{ %\function{v3d_sphere_rlup} %\synopsis{Fill a v3darr within a spherical shell from lookup vs r.} %\usage{ v3darr = v3d_sphere_rlup (rin, rout, rlups, vlups [,oxyz ]); } %\description % Creates a spherical shell with values as a function of radius as % given by a lookup table. The parameters are: % rlups and vlups are 1d arrays of radii and values % % When the lookup tables are on a much finer scale than the v3d resolution % it is useful to use v3d_smooth_box to average the vlups array to the % appropriate scale. %\notes % Uses the gsl routine interp_linear to do the lookup and is not % compiled unless inter_linear is present. %\seealso{ } %!%- }}} { variable rin, rout, rlups, vlups, oxyz; switch (_NARGS) { case 4: oxyz = [0.,0.,0.]; vlups = (); rlups = (); rout = (); rin = (); } { case 5: oxyz = (); vlups = (); rlups = (); rout = (); rin = (); } { % otherwise: message(" rlups and vlups are 1d arrays of radii and values."); message(" uses the gsl routine interp_linear for the lookup"); message(" (use v3d_smooth_box to average the vlups if desired)"); usage("................. v3d_sphere_rlup (rin, rout, rlups, vlups [,oxyz ]);"); } % defaults for NULL primary arguments: if ( rin == NULL ) rin = 0.0; % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the system is at: % oxyz % nominal variables... variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable rad2limit = (rout+5.0*cutoff_scale)^2; variable iz, iy, rzsquared, rysquared; variable these_rs, these_lups; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { rzsquared = (v3d_va[iz]-oxyz[2])^2; if (rzsquared < rad2limit) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { rysquared = (v3d_va[iy]-oxyz[1])^2; if ( (rzsquared+rysquared) < rad2limit) { these_rs = (v3d_va - oxyz[0])^2 + (rzsquared+rysquared) ; % only do this if some are less than rad2limit if ( min(these_rs) < rad2limit ) { these_rs = sqrt(these_rs); % Now apply the lookup function to these values: these_lups = interp_linear(these_rs, rlups, vlups); % If rin is very small or zero then only need the rout cutoff: if ( rin < cutoff_scale ) { val[*,iy,iz] = these_lups * v3d_soft_cutoff( (these_rs - rout) / cutoff_scale ) ; } else { % cutoff in and out: val[*,iy,iz] = these_lups * v3d_soft_cutoff( (these_rs - rout) / cutoff_scale ) * v3d_soft_cutoff( (rin - these_rs) / cutoff_scale ) ; } } % if min(these_rs ...) } % if ... } % y loop } % if ... } % z loop return val; } % }}} #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_sphere_line () % {{{ %!%+{{{ %\function{v3d_sphere_line} %\synopsis{Fill a v3darr with a line of spheres} %\usage{ v3darr = v3d_sphere_line (p1xyz, p2xyz [, rout [, nspheres]]); } %\description % A line of nspheres is produced in the v3darr. The first and % last sphere are centered at the points p1xyz and p2xyz, % specified as a simple array: pNxyz = [px,py,pz] . % %\notes % %\seealso{ v3d_sphere_ring} %!%- }}} { variable p1xyz, p2xyz, rout, nspheres; variable ptopdist, fracs, oxs, oys, ozs, ip; switch (_NARGS) { case 2: rout = NULL; nspheres = NULL; p2xyz = (); p1xyz = (); } { case 3: nspheres = NULL; rout = (); p2xyz = (); p1xyz = (); } { case 4: nspheres = (); rout = (); p2xyz = (); p1xyz = (); } { % otherwise: message(" pNxyz = [px,py,pz]"); usage("................. v3d_sphere_line (p1xyz, p2xyz [, rout [, nspheres]]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % default sphere radius: if ( rout == NULL ) rout = v3d_cube_radius/v3d_na; % distance between the points ptopdist = sqrt(sum((p1xyz-p2xyz)^2)); if ( nspheres == NULL ) nspheres = max( [1+int(ptopdist/(4.0*rout)), 3]); % fractions along the line, nspheres from 0 to 1.0: fracs = [0:nspheres-1:1]/(nspheres - 1.0); % locations of the spheres oxs = fracs * (p2xyz[0]-p1xyz[0]) + p1xyz[0]; oys = fracs * (p2xyz[1]-p1xyz[1]) + p1xyz[1]; ozs = fracs * (p2xyz[2]-p1xyz[2]) + p1xyz[2]; % and create an array with them... ip = 0; val = v3d_sphere(0.0, rout, [oxs[ip],oys[ip],ozs[ip]]); _for ip (1, nspheres-1, 1) { val = v3d_merge_max(val, v3d_sphere(0.0, rout, [oxs[ip],oys[ip],ozs[ip]]) ); } return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_sphere_ring () % {{{ %!%+{{{ %\function{v3d_sphere_ring} %\synopsis{Fill a v3darr with a ring of spheres} %\usage{ v3darr = v3d_sphere_ring (r_ring, r_sph, nspheres [,oxyz [, rthphi]]); } %\description % A ring of nspheres is produced in the v3darr. % The circualr ring has center at oxyz and normal % in the direction rthphi. %\notes % %\seealso{ v3d_sphere_line} %!%- }}} { variable r_ring, r_sph, nspheres, oxyz, rthphi; variable def_rthphi = [270.0, 0.0]; variable angs, oxs, oys, ozs, ip; switch (_NARGS) { case 1: rthphi = def_rthphi; oxyz = [0.,0.,0.]; nspheres = NULL; r_sph = NULL; r_ring = (); } { case 2: rthphi = def_rthphi; oxyz = [0.,0.,0.]; nspheres = NULL; r_sph = (); r_ring = (); } { case 3: rthphi = def_rthphi; oxyz = [0.,0.,0.]; nspheres = (); r_sph = (); r_ring = (); } { case 4: rthphi = def_rthphi; oxyz = (); nspheres = (); r_sph = (); r_ring = (); } { case 5: rthphi = (); oxyz = (); nspheres = (); r_sph = (); r_ring = (); } { % otherwise: usage("................. v3d_sphere_ring (r_ring, r_sph, nspheres [,oxyz [, rthphi]]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable rotang, rotx, roty, rotz; % default sphere radius: if ( r_sph == NULL ) r_sph = 2.0*v3d_cube_radius/v3d_na; % default number if ( nspheres == NULL ) nspheres = 6; if ( oxyz == NULL ) oxyz = [0.,0.,0.]; if ( rthphi == NULL ) rthphi = def_rthphi;; % angles around the ring, nspheres from 0 to 1.0-1/N: angs = (2.0*PI*[0:nspheres-1:1])/(1.0*nspheres); % Locate them in the X=0 plane with first one at [y,z] = [0,-1.0*r_ring] ozs = -1.0*r_ring*cos(angs); oys = r_ring*sin(angs); oxs = 0.0 * ozs; % Now rotate about Z by phi in the X-Y plane... rotang = rthphi[1]*DTOR; rotx = oxs*cos(rotang) - oys*sin(rotang); roty = oxs*sin(rotang) + oys*cos(rotang); oxs = rotx; oys = roty; % and then about Y by theta in the X-Z plane... rotang = rthphi[0]*DTOR; rotx = oxs*cos(rotang) + ozs*sin(rotang); rotz = -1.0*oxs*sin(rotang) + ozs*cos(rotang); oxs = rotx; ozs = rotz; % and include the offsets oxs = oxs + oxyz[0]; oys = oys + oxyz[1]; ozs = ozs + oxyz[2]; % and create an array with them... ip = 0; val = v3d_sphere(0.0, r_sph, [oxs[ip],oys[ip],ozs[ip]]); _for ip (1, nspheres-1, 1) { val = v3d_merge_max(val, v3d_sphere(0.0, r_sph, [oxs[ip],oys[ip],ozs[ip]]) ); } return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_spheroid () % {{{ %!%+{{{ %\function{v3d_spheroid} %\synopsis{Fill a v3darr with 1.0 within a spheroid} %\usage{ v3darr = v3d_spheroid (rin, rout, polar_ratio [,oxyz [, rthphi]]); } %\description % Creates a spheroid in the v3darr based on the parameters: % % rin and rout are the inner and outer radii of the % spheroid in its equatorial plane. % polar_ratio is the ratio of the polar outer radius % to the equatorial outer radius. % % Centered at OXYZ % The polar axis of the spheroid is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Y. %\notes % How to transform a non-zero rin as the shape moves out of the equatorial % plane, in particular rin at the pole, doesn't have a uniquely obvious default. % In the implementation here, when polar_ratio is greater than 1.0, % the polar "thickness" of the shell is the same as the equatorial thickness, % i.e., when the poles are stretched out they keep a constant thickness. % When polar_ratio is less than 1.0 (the poles are squashed) then the polar % thickness is reduced by the polar_ratio factor. % % Thus, "nesting" of such spheroids will have gaps between them; to create % gapless nested spheroids create them as solid spheroids (rin=0.0) and % then subtract the inner one from the outer. % %\seealso{ } %!%- }}} { variable rin, rout, polar_ratio, oxyz, rthphi; variable def_rthphi = [0.0, 90.0]; switch (_NARGS) { case 3: rthphi = def_rthphi; oxyz = [0.,0.,0.]; polar_ratio = (); rout = (); rin = (); } { case 4: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; polar_ratio = (); rout = (); rin = (); } { case 5: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; polar_ratio = (); rout = (); rin = (); } { % otherwise: message(" rin, rout are in equatorial plane of spheroid,"); message(" polar_ratio is r_polar/r_equatorial"); usage("................. v3d_spheroid (rin, rout, polar_ratio [,oxyz [, rthphi]]);"); } % catch some values if (rin == NULL) rin = 0.0; % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Define some useful variables and values variable requiv, pratin; % It's tricky to decide how to define the inside spheroid limit % in the polar direction: just a polar_ratio-scaled rin ? % That would make the thickness vary from equator to pole... % Instead let's try making the polar thickness equal to % the equatorial thickness (rout-rin) by defining an % appropriate polar_ratio for the inner boundary. % Use different defaults for oblate and prolate cases... if (rin > 0.0) { if (polar_ratio < 1.0) { % have the rin in the polar direction scale as rout pratin = polar_ratio; } else { % keep roughly constant thickness from equator to pole: pratin = ((rout*polar_ratio)-(rout-rin)) / rin; } } % nominal variables... variable axial, radial; %%variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy, rad2limit; rad2limit = (rout*max([1.0,polar_ratio])+5.0*cutoff_scale)^2; % Loop over the 'toward-away' 'front-back' (z) value % and only do relevant planes _for iz ( 0, 2*v3d_na, 1) { if ( (v3d_va[iz]-oxyz[2])^2 < rad2limit ) { % Loop over the 'vertical (up, North)' (y) value % and skip out of range ones... _for iy ( 0, 2*v3d_na, 1) { if ( ( (v3d_va[iz]-oxyz[2])^2 + (v3d_va[iy]-oxyz[1])^2 ) < rad2limit ) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; radial = sqrt(radial) ; % Fill val with uniform value (1.0) inside the spheroid % The equivalent radius uses scaled axial values: requiv = sqrt(radial^2+(axial/polar_ratio)^2); val[*,iy,iz] = v3d_soft_cutoff( ( requiv - rout) / cutoff_scale ) ; % remove the inside if rin > 0.0 if (rin > 0.0) { requiv = sqrt(radial^2+(axial/pratin)^2); val[*,iy,iz] = val[*,iy,iz] * v3d_soft_cutoff( (rin - requiv) / cutoff_scale ) ; } } % end of y if } % y loop } % end of y if } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_torus () % {{{ %!%+{{{ %\function{v3d_torus} %\synopsis{Fill a v3darr with 1.0 within a torus} %\usage{ v3darr = v3d_torus (rin, rout [,oxyz [, rthphi]]); } %\description % Creates a torus in the v3darr based on the parameters: % % rin and rout are the inner and outer radii of the solid % torus in its equatorial plane. % % Centered at OXYZ % Axis of the torus is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Z. %\notes % %\seealso{ } %!%- }}} { variable rin, rout, oxyz, rthphi; variable def_rthphi = [270.0, 0.0]; switch (_NARGS) { case 2: rthphi = def_rthphi; oxyz = [0.,0.,0.]; rout = (); rin = (); } { case 3: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; rout = (); rin = (); } { case 4: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; rout = (); rin = (); } { % otherwise: message(" rin, rout are in equatorial plane of torus"); usage("................. v3d_torus (rin, rout [,oxyz [, rthphi]]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Define some useful constant values variable rave = (rout+rin)/2.0; % square of the "little" radius: variable rlittle2 = (rout - rave)^2; variable deltar2,absdr2, deltars, limin, limout; % nominal variables... variable axial, radial; %%variable axfact; variable negvals; variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; radial = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 - axial^2 ); negvals = where( radial < 0.0 ); if (length(negvals) > 0) radial[negvals] = 0.0; radial = sqrt(radial) ; % Now use the axial and radial values to define the desired shape. % Fill val with uniform value (1.0) inside the torus % inner and outer radial limits for these points are % +/- a delta r from the rave value, delta r(s) % depends on the axial value(s): deltar2 = rlittle2 - axial^2; absdr2 = abs(deltar2); deltars = (deltar2/absdr2) * sqrt(absdr2); limin = rave - deltars; limout = rave + deltars; val[*,iy,iz] = v3d_soft_cutoff( (radial - limout) / cutoff_scale ) * v3d_soft_cutoff( (limin - radial) / cutoff_scale ) ; } % y loop } % z loop return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_2dto3d () % {{{ %!%+{{{ %\function{v3d_list} %\synopsis{Fill a spherical shell with values from a 2D array} %\usage{ v3darr = v3d_2dto3d (arr2din, rlo, rhi [,oxyz [, rthphi]]); } %\description % Fill a spherical shell with values based on a 2D array % which gives the values in a grid of lattitude-radii: % arr2din [ ilattitude_0_to_90, iradius_rlo_rhi ] % Assumes the lattitude is complete from 0 to 90 degrees % and is reflected for negative lattitudes. % The user provides rlo and rhi to set the radius scale. % % Centered at OXYZ % Axis of the object (lattitude = 90 degrees) is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Y. %\notes % %\seealso{ } %!%- }}} { variable arr2din, rlo, rhi, oxyz, rthphi; variable def_rthphi = [0.0, 90.0]; switch (_NARGS) { case 3: rthphi = def_rthphi; oxyz = [0.,0.,0.]; rhi = (); rlo = (); arr2din = (); } { case 4: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; rhi = (); rlo = (); arr2din = (); } { case 5: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; rhi = (); rlo = (); arr2din = (); } { % otherwise: message(" 2D array: arr2din [ i_lattitude, i_radius ] "); message(" lat. range of i_lat. is assumed 0 to PI/2."); message(" radius range of i_radius is rlo to rhi."); usage("................. v3d_2dto3d (arr2din, rlo, rhi [,oxyz [, rthphi]]);"); } % Defaults for NULL arguments: if (rlo == NULL) rlo = 0.0; % Get info on the 2D array: variable nrads; variable latlo, lathi, nlats; % number of radii nrads = length(arr2din[0,*]); % Set the array lattitude range from 0 to PI/2. latlo = 0.0; lathi = PI/2.0; % number of lat.s nlats = length(arr2din[*,0]); % Each 3D cell is filled with the rough average of the % 2D array values that are included in the 3D cell... % Do this by pre-smoothing (averaging) the 2D array on % the appropriate scales based on the size of the 3D cells. % Note that the smoothing varies % with r because the lattitude-->distance depends on r. % will use these variables later as well. variable radius, irad, lat, ilat; variable nsmooth, cell_size, rad_spacing, lat_spacing; cell_size = v3d_cube_radius / v3d_na; rad_spacing = (rhi-rlo) / (nrads-1.0); % lat spacing is in radians lat_spacing = (lathi-latlo) / (nlats-1.0); variable arr2d = 0.0*arr2din; % indices are: arr2din [ i_lattitude, i_radius ] "); % use v3d_smooth_box (arr1d, nsmooth) % loop over the radial dimension, smoothing in lattitude radius = rlo; _for irad (0, nrads-1, 1) { % nsmooth = cell size / size betweem lat cells at this radius % ignore radius < = 0 (allows rlo to be < 0 ) if (radius > 0.0 ) { nsmooth = cell_size/(lat_spacing*radius) + 0.3; arr2d[*,irad] = v3d_smooth_box(arr2din[*,irad], nsmooth) ; } radius = radius + rad_spacing; } % now, loop over the lattitudes, smoothing along radius _for ilat (0, nlats-1, 1) { % nsmooth = cell size / radial cell size nsmooth = cell_size/rad_spacing + 0.3; arr2d[ilat,*] = v3d_smooth_box(arr2d[ilat,*], nsmooth) ; } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; % Origin of the cylinderical system is at: % oxyz % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % nominal variables... variable axial, rad3dsq; variable axfact; %%variable negvals; %%variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable rlosq = rlo^2; % catch negative rlos: if ( rlo < 0.0 ) rlosq = 0.0; variable rhisq = rhi^2; variable iz, iy, ix, sel, is; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { % Now calculate axial and radial coordinates of the % data cube center from the cylinder system origin. axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; % * Use the 3D radius rather than the radial value... rad3dsq = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 ); % Select the cells that are in the array radii limits sel = where( (rad3dsq > rlosq) and (rad3dsq < rhisq) ); if (length(sel) > 1 ) { % go through each of the selected cells ... _for is (0, length(sel)-1, 1) { % this ix value ix = sel[is]; % calculate the radius index: radius = sqrt(rad3dsq[ix]); irad = int( (radius-rlo)/rad_spacing ); % calculate the lattitude index: lat = asin( abs( axial[ix]/radius ) ); ilat = int( (lat-latlo)/lat_spacing ); val[ix,iy,iz] = arr2d[ilat,irad]; } % for is ... } % if (length(sel)) } % y loop } % z loop return val; %, arr2d; % can return the smoothed 2d array too } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #ifexists hist2d public define v3d_histxy () % {{{ %!%+{{{ %\function{v3d_histxy} %\synopsis{Returns a 3d array with a 2D histogram in central z plane.} %\usage{ v3darr = v3d_histxy (xs, ys); } %\description % The input arrays of x and y coordinates are binned into a % 2D histogram at the central-z plane of the returned 3D array. % This is a useful way of using real data to generate a % complex "shape". %\notes % The resulting array can be mapped to a 3D shape/volume using % v3d_reproject. %\seealso{v3d_reproject} %!%- }}} { if (_NARGS != 2) usage("................. v3darr = v3d_histxy (xs, ys);"); variable ys = (); variable xs = (); % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable dgrid, xygrid, maxgrid, sel; % Form the histogram grids from the v3d value: % v3d_va : the spatial coordinate value along an axis, % specified at the _center_ of the cell. dgrid=v3d_va[1]-v3d_va[0]; xygrid = v3d_va - dgrid/2.0; maxgrid = v3d_va[-1]+dgrid; % Do not let the max bins contain the excess: sel = where((xs < maxgrid) and (ys < maxgrid)); % the middle xy plane is with iz = v3d_na : val[*,*,v3d_na] = hist2d (xs[sel], ys[sel], xygrid, xygrid); return val; } % }}} #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_reproject () % {{{ %!%+{{{ %\function{v3d_reproject} %\synopsis{Map one v3d array to another retaining x,y projection.} %\usage{ arr_out = v3d_reproject (arr_in, arr_shape); } %\description % Create an output v3d array which has the same x,y projection as % the input array but has the 3D shape of the arr_shape array. % This is useful to map the (essentially 2D) output of v3d_histxy % onto a 3D shape. %\notes %\seealso{v3d_histxy} %!%- }}} { if (_NARGS != 2) usage("................. arr_out = v3d_reproject (arr_in, arr_shape);"); variable arr_shape = (); variable arr_in = (); % scale arr_shape to have max value of 1.0: arr_shape = arr_shape/max(arr_shape); variable arr_out = 0.0*arr_in; variable ix, iy, sumz, zvector; _for ix (0, 2*v3d_na, 1) { _for iy (0, 2*v3d_na, 1) { sumz = sum(arr_shape[ix,iy,*]); if (sumz < 1.e-3) { % very small shape along this ix,iy - use 0: zvector = 0.0*arr_shape[ix,iy,*]; } else { % create the normalized z-vector from the shape array zvector = arr_shape[ix,iy,*]/sumz; } % Now implement the desired output array values: % arr_out at x,y along z = % arr_in x,y summed over z * arr_shape x,y normalized along z (or 0) arr_out[ix,iy,*] = sum(arr_in[ix,iy,*]) * zvector; } } return, arr_out; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_f_generic () % {{{ %!%+{{{ %\function{v3d_f_generic} %\synopsis{Calculate a general function of coordinates} %\usage{ v3darr = v3d_f_generic ("dummy_f_string" [,oxyz [, rthphi]]); } %\description % Evaluate a user-provided function of Cartesian, cylinderical, and/or % spherical coordinates. % Available variables in the rotated system are: % (all of these are arrays along the ix cube dimension) % Cartesian: % axial [=xprime], yprime, zprime % spherical: % rad3d and rad3dsq[uared], lattitude, azimuth % cylindrical % radcyl, axial, azimuth % % At present the supplied string is a dummy argumment and the user % must use this routine as a template and modify the line: % val[*,iy,iz] = ... ; % to compute the desired function of the coordinates. % The computation of unneeded variables can be commented out or removed % as well to speed evaluation. % % The coordinates are centered at OXYZ % The axis of the object is given by: % rthphi = [theta, phi] where: % theta, phi are the polar coordinates of the axis: % theta is angle in degrees from X toward -Z axis. % phi is angle in +/- degrees from X-Z plane toward +/- Y axis % Default is Axis along +Y. %\notes % The goal was to be able to define the general function by passing % in a function string or a reference to a user-created function: % f_string = "(cos(3.0*azimuth))^2 * (sin(5.0*lattitude))^2 *" + % or "exp(-1.0*rad3dsq/(40.0)^2)"; % f_string = &f_to_use % and the routine would determine which type of f_string is given. % This may be a bit tricky to do... Hence the suggestion to use this as % a template. %\seealso{ } %!%- }}} { variable f_string, oxyz, rthphi; variable def_rthphi = [0.0, 90.0]; switch (_NARGS) { case 1: rthphi = def_rthphi; oxyz = [0.,0.,0.]; f_string = (); } { case 2: rthphi = def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; f_string = (); } { case 3: rthphi = (); if (rthphi == NULL) rthphi=def_rthphi; oxyz = (); if (oxyz == NULL) oxyz=[0.,0.,0.]; f_string = (); } { % otherwise: message(" generic function with rect, cyl, and spherical coord.s"); message(" variables available are: rad3dsq, radius, axial, radial, lattitude"); usage("................. v3d_f_generic (f_string [,oxyz [, rthphi]]);"); } % The 3-D volume array: variable val = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; message(" Function string = "+f_string); % Origin of the cylinderical system is at: % oxyz % Unit vector in axial ( = x-prime) direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); % Unit vector in the y-prime (azimuth=0) direction variable ypuvx, ypuvy, ypuvz; ypuvx = sin(rthphi[1]*DTOR)*(-1.0)*cos(rthphi[0]*DTOR); ypuvy = cos(rthphi[1]*DTOR); ypuvz = sin(rthphi[1]*DTOR)*sin(rthphi[0]*DTOR); % Unit vector in the z-prime (az = 90 degrees) direction variable zpuvx, zpuvy, zpuvz; zpuvx = sin(rthphi[0]*DTOR); zpuvy = 0.0; zpuvz = cos(rthphi[0]*DTOR); % nominal variables... variable axial, yprime, zprime; variable rad3dsq, rad3d, radcyl, lattitude, azimuth; variable axfact; variable negvals; %%variable cutoff_scale = v3d_cube_radius/(2.0*v3d_na); variable iz, iy; % Loop over the 'toward-away' 'front-back' (z) value _for iz ( 0, 2*v3d_na, 1) { % Loop over the 'vertical (up, North)' (y) value _for iy ( 0, 2*v3d_na, 1) { % Calculate the cartesian, cylinderical, and spherical coord.s % for the whole array of cells: [*,iy,iz] % % axial ( = x_prime) values: axial = (v3d_va - oxyz[0])*axuvx + (v3d_va[iy]-oxyz[1])*axuvy + (v3d_va[iz]-oxyz[2])*axuvz ; % y_prime values: yprime = (v3d_va - oxyz[0])*ypuvx + (v3d_va[iy]-oxyz[1])*ypuvy + (v3d_va[iz]-oxyz[2])*ypuvz ; % z_prime values: zprime = (v3d_va - oxyz[0])*zpuvx + (v3d_va[iy]-oxyz[1])*zpuvy + (v3d_va[iz]-oxyz[2])*zpuvz ; % % 3d radius (and squared) values rad3dsq = ( (v3d_va-oxyz[0])^2 + (v3d_va[iy]-oxyz[1])^2 + (v3d_va[iz]-oxyz[2])^2 ); % and 3d radius values rad3d = sqrt(rad3dsq); % % cylinderical coor 2d radial values: radcyl = ( rad3dsq - axial^2 ); negvals = where( radcyl < 0.0 ); if (length(negvals) > 0) radcyl[negvals] = 0.0; radcyl = sqrt(radcyl) ; % % cylinderical "lattitude" is 0 in the equator of the shape % and 90 degrees along the axi[al]s . lattitude = asin( axial/rad3d ); % % Azimuth from y-prime toward z-prime: azimuth = atan2(zprime, yprime); % - - - - - - - - - - % Generic function computation goes here. % val[*,iy,iz] = ... ; % % Available variables in the rotated system are: % (all of these are arrays along the ix cube dimension) % Cartesian: % axial [=xprime], yprime, zprime % Spherical: % rad3d and rad3dsq[uared], lattitude, azimuth % Cylindrical % radcyl, axial, azimuth val[*,iy,iz] = (cos(3.0*azimuth))^2 * (sin(5.0*lattitude))^2 * exp(-1.0*rad3dsq/(40.0)^2); % - - - - - - - - - - } % y loop } % z loop % Remove any NaNs... ( use the negvals variable ) negvals = where( isnan(val) ); if ( length(negvals) > 0 ) { % tell people about it FYI message(" v3d_f_generic: Found NaN "+string(length(negvals))+" times; replacing with 0.0s"); message(" First NaN is at index = "+string(negvals[0])); val[ negvals ] = 0.0; } return val; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % volview routine is needed here #ifexists volview public define v3d_view () % {{{ %!%+{{{ %\function{v3d_view} %\synopsis{Display a 3D volume with the volview guilet} %\usage{ v3d_view (v3darr [, file_name_h5]); } %\description % The 3D volume array is conditioned and passed to volview for display. % Specifically the "correct" xyz orientation is achieved: % +X is to the right, % +Y is "up" and % +Z is toward the observer. % Negative values are clipped to 0.0 . This allows a simple thresholding % view of the volume: % v3d_view ( v3darr - threshold_value ); % Also, a sqrt() function is applied to reduce the dynamic range of % the array. % % A useful way to save and share the viewed 3D model % is through the optional filename: if present then the array passed to % volview is written to the filename by h[df]5_write. % This file can then be re-viewed with the command: % volview(h5_read(file_name_h5)); %\notes % This routine is only compiled if volview is available. % The ability to write h5 files requires h5_write - but the routine % will compile with out it. %\seealso{ v3d_project} %!%- }}} { variable v3darr, fileout=NULL; switch (_NARGS) { case 1: v3darr = (); } { case 2: fileout=(); v3darr = (); } { % otherwise message(" Write the 3D volview'ed array to a file, if provided."); message(" ( Can re-view the saved model volume with:"); message(" vol_view(h5_read(''filename.h5'')); ) "); usage("................. v3d_view(v3darr [, ''filename.h5''])"); } variable newvol = Float_Type[2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable ix, iz, temp_vals; _for ix (0, 2*v3d_na, 1) { _for iz (0, 2*v3d_na, 1) { % switch x,y,z order to z,y,x % and reverse the y axis % and the new z-axis % *** include a sqrt() here to help better show % lower-level values *** % This also "removes" negative values, becoming NaNs... temp_vals = sqrt( reverse( v3darr[ix,*,iz] ) ); temp_vals[where( isnan(temp_vals) )] = 0.0; newvol[2*v3d_na-iz,*,ix] = temp_vals; } } % et voila volview(newvol); % The hdf routine h5_write is used here #ifexists h5_write % If we have the h5 routines available, % save it for future viewing if (fileout != NULL) { h5_write(fileout,newvol); } #endif } % }}} #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The routine imdisplay, based on gtk, is used here. #ifexists imdisplay public define v3d_project () % {{{ %!%+{{{ %\function{v3d_project} %\synopsis{Show a projected [color-coded] image of the 3D array(s)} %\usage{ v3d_project ( [rarr], [garr], [barr]); } %\description % This routine shows a "thin" projection of a v3d volume % as seen from the nominal viewing location: % the observer is at a large +Z location looking along -Z. % The resulting image has +X to the right and +Y "up". % % A single v3d array can be passed in for monochrome display, % or two or three v3d arrays for RGB colors. % % Array values that are negative are clipped to 0.0 . %\notes % This routine requires imdisplay from slgtk. % % The nominal data sent to imdisplay is a 3d array % of the form UChar_Type[nypix,nxpix,3] where the % three "planes" are R,G,B and values are from 0 to 255. %\seealso{ v3d_view} %!%- }}} { variable rarr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable garr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable barr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; switch (_NARGS) { case 1: % default to green for monochrome garr = (); } { case 2: garr = (); rarr = (); } { case 3: barr = (); garr = (); rarr = (); } { % otherwise: usage("................. v3d_project (rarr [,garr [, barr]]);"); } % allow some to be missing... if ( rarr == NULL ) rarr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; if ( garr == NULL ) garr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; if ( barr == NULL ) barr = Float_Type [2*v3d_na+1, 2*v3d_na+1, 2*v3d_na+1]; variable projr = Float_Type [2*v3d_na+1, 2*v3d_na+1]; variable projg = Float_Type [2*v3d_na+1, 2*v3d_na+1]; variable projb = Float_Type [2*v3d_na+1, 2*v3d_na+1]; variable imgclr = UChar_Type [2*v3d_na+1, 2*v3d_na+1, 3]; variable ix, iy; variable imgscale; % Remove any negative values % = clip at zero % May not always want this, but for now... rarr[where(rarr < 0.0)] = 0.0; garr[where(garr < 0.0)] = 0.0; barr[where(barr < 0.0)] = 0.0; % For each color plane, % loop over x and y and % make the projection as a sum over z _for ix (0, 2*v3d_na, 1) { _for iy (0, 2*v3d_na, 1) { % switch x and y and change direction of y: projr[2*v3d_na-iy,ix] = sum(rarr[ix,iy,*]); projg[2*v3d_na-iy,ix] = sum(garr[ix,iy,*]); projb[2*v3d_na-iy,ix] = sum(barr[ix,iy,*]); } } % convert the 2d floats to a scaled image: imgscale = 250.0/max([projr,projg,projb]); imgclr[*,*,0] = typecast(imgscale*projr, UChar_Type); imgclr[*,*,1] = typecast(imgscale*projg, UChar_Type); imgclr[*,*,2] = typecast(imgscale*projb, UChar_Type); % and show it; use the same pixel size as volview if (save_image_file == NULL) { % show to the screen: imdisplay(imgclr,"panes=one,fill=scale,size=256x256"); } else { % write the image to a file, and reset: imdisplay(imgclr,"panes=one,fill=scale,size=256x256,save="+save_image_file); save_image_file = NULL; } } % }}} #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% define clip_max (array, clip_level) { variable above = where(array > clip_level); array[above] = clip_level; return array; } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The routine imdisplay, based on gtk, is used here % as well as the hist2d routine from histogram module. #ifexists imdisplay #ifexists hist2d public define v3d_evtimg () % {{{ %!%+{{{ %\function{v3d_evtimg} %\synopsis{Bin and display a set of events with color coding} %\usage{ v3d_evtimg (xs, ys [, es [, Rrange, Grange, Brange]]); } %\description % The arrays of points xs,ys are binned in 2D and the histogram % is displayed. If the optional es array (e.g., "energy") is % provided, then the R-G-B-ranges are used to assign events to % R-G-B image planes for color coding, e.g., % Rrange = [rmin,rmax] selects the es for the R plane % % A small amount of smoothing is applied to the image planes to % reduce graininess in the image. % % The dynamic range and color balancing (plane-to-plane ratio) % can be choosen from a few built-in options for disp_mode: % % disp_mode > 0 - Use the image planes' histogram as created. % disp_mode < 0 - Apply a sqrt() to the planes' counts before further % procssing. % for different values of abs(disp_mode) do: % 1 - nominal linear scaling to a common max. % 2 - scaling to each plane's max. % 3 - Dan's special: retain color range but compress intensity range. % %\notes % This routine requires imdisplay from slgtk and also % hist2d from the histogram module. % % This is a very-lite version of the evt2img.sl tool which % displays fits files with interactive RGB selection, etc. %\seealso{ v3d_mc_points, v3d_doppler_points} %!%- }}} { variable xs, ys, es=NULL, rrange=NULL, grange=NULL, brange=NULL, disp_mode=NULL; % Do a little smoothing of the image planes... variable nsmooth=2; % if just the x,y [,e] are given then black and white % otherwise something for each color is given (can be NULL) switch (_NARGS) { case 2: ys = (); xs = (); } { case 3: es = (); ys = (); xs = (); } { case 6: brange = (); grange = (); rrange = (); es = (); ys = (); xs = (); } { case 7: disp_mode = (); brange = (); grange = (); rrange = (); es = (); ys = (); xs = (); } { % otherwise message(" Rrange = [rmin,rmax] selects on es for the R plane, etc."); usage("................. v3d_evtimg (xs, ys [, es [, Rrange, Grange, Brange [, disp_mode]]]);"); } % If no es are given, set them to 1.0 if (es == NULL) es = 1.0 + 0.0*xs; % If all three ranges are NULL, then set them all % to a large range - this will make a white/greyscale image. if (rrange==NULL and grange==NULL and brange==NULL) { % These range values refer to the es values... rrange = [0.0,1.e18]; grange = [0.0,1.e18]; brange = [0.0,1.e18]; } else { % If some plane has ranges provided, then any NULL ranges % should be exclusionary (min greater than max): if (rrange == NULL) rrange = [2.0,0.0]; if (grange == NULL) grange = [2.0,0.0]; if (brange == NULL) brange = [2.0,0.0]; } % Set a default disp_mode if not provided: if (disp_mode == NULL) { disp_mode = 1; } % Now go through the three color planes and make that plane's image % based on the selected events. % The hist2d grids are "bin lo" values: generate them from the v3d_va % bin center values and also eliminate the bin overflow of last bin. variable sel, imgr, imgg, imgb; variable grid_spacing, grid_lows, bin_highest; grid_spacing = v3d_va[1] - v3d_va[0]; grid_lows = v3d_va - 0.5*grid_spacing; bin_highest = grid_lows[-1] + grid_spacing; % do the selections on energy and bin_highest: sel = where( (es > rrange[0]) and (es < rrange[1]) and (xs < bin_highest) and (ys < bin_highest) ); imgr = double ( hist2d(ys[sel],xs[sel],grid_lows,grid_lows) ); sel = where( (es > grange[0]) and (es < grange[1]) and (xs < bin_highest) and (ys < bin_highest) ); imgg = double ( hist2d(ys[sel],xs[sel],grid_lows,grid_lows) ); sel = where( (es > brange[0]) and (es < brange[1]) and (xs < bin_highest) and (ys < bin_highest) ); imgb = double ( hist2d(ys[sel],xs[sel],grid_lows,grid_lows) ); % Smooth them a bit... imgr = v3d_smooth_2d(imgr,nsmooth); imgg = v3d_smooth_2d(imgg,nsmooth); imgb = v3d_smooth_2d(imgb,nsmooth); % sqrt scale the planes if negative type if (disp_mode < 0) { % and scale within the planes imgr = sqrt(imgr); imgg = sqrt(imgg); imgb = sqrt(imgb); } % Convert the planes to image values depending on the % selected disp_mode % Generally useful variables; compo is the composite RGB image array. variable maxall, compo = UChar_Type [2*v3d_na+1,2*v3d_na+1,3]; % if (abs(disp_mode) == 1) { % normalize to the overall max maxall = max([imgr,imgg,imgb]); compo[*,*,0] = typecast(255.0*imgr/maxall, UChar_Type); compo[*,*,1] = typecast(255.0*imgg/maxall, UChar_Type); compo[*,*,2] = typecast(255.0*imgb/maxall, UChar_Type); } if (abs(disp_mode) == 2) { % normalize each plane independently to their max maxall = max(imgr); compo[*,*,0] = typecast(255.0*imgr/maxall, UChar_Type); maxall = max(imgg); compo[*,*,1] = typecast(255.0*imgg/maxall, UChar_Type); maxall = max(imgb); compo[*,*,2] = typecast(255.0*imgb/maxall, UChar_Type); } if (abs(disp_mode) == 3) { % Retain the original color range % while compressing the intensity range with a double sqrt. % If we're not doing the sqrt on the planes % then do something intermediate between overall max % and per-plane max: bring the total counts in the planes % closer together in a sqrt sense. % (If a sqrt has been applied to the data already then this % adjustement becomes too much, hence only for mode +3.) if ( disp_mode == 3) { % bring the total counts in each plane closer together variable totr, totg, totb, totmax; totr = sum(imgr); totg = sum(imgg); totb = sum(imgb); totmax = max([totr,totg,totb]); if (totr > 0) imgr = imgr / sqrt(totr/totmax); if (totg > 0) imgg = imgg / sqrt(totg/totmax); if (totb > 0) imgb = imgb / sqrt(totb/totmax); } % normalize to total counts variable tot_intens = imgr + imgg + imgb; % set the fraction of each color in each bin: variable ifrac, fracr, fracg, fracb, fracmax; fracr = imgr / (tot_intens+1.0); fracg = imgg / (tot_intens+1.0); fracb = imgb / (tot_intens+1.0); % apply a sqrt to this intensity tot_intens = sqrt(tot_intens); % again for the 3 case: %if (disp_mode == 3) tot_intens = sqrt(tot_intens); tot_intens = sqrt(tot_intens); % load in the values retaining the pre-sqrt color fractions: % Apply a brightening factor and clip values above 255 (1.0). maxall = max(tot_intens); variable brtn = 1.6; compo[*,*,0] = typecast(255.0*clip_max(brtn*fracr*tot_intens/maxall, 1.0), UChar_Type); compo[*,*,1] = typecast(255.0*clip_max(brtn*fracg*tot_intens/maxall, 1.0), UChar_Type); compo[*,*,2] = typecast(255.0*clip_max(brtn*fracb*tot_intens/maxall, 1.0), UChar_Type); } % Need to change the direction of the compo[site] y index: % imgX(*,ix) = reverse( imgX(*,ix) ); variable ic, ix; _for ic (0,2,1) { _for ix (0,2*v3d_na,1) { compo[*,ix,ic] = reverse( compo[*,ix,ic] ); } } % and show it if (save_image_file == NULL) { % show to the screen: imdisplay(compo,"panes=one,fill=scale"); } else { % write the image to a file, and reset: imdisplay(compo, "panes=one,fill=scale,save="+save_image_file); save_image_file = NULL; } } % }}} #endif #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % - - - - - - - - Monte Carlo routines require urand - - - - - % If urand (and seed_random) doesn't exist then make a version from gsl ran_flat: % #ifnexists urand #ifexists ran_flat % Make urand and seed_random using GSL routines: % select the default generator method variable ran_gen; ran_gen = rng_alloc(); define urand () % {{{ %!%+{{{ %\function{urand} %\synopsis{Generate uniform random numbers from 0.0 to 1.0} %\usage{ rand_values = urand (n_samples); } %\description % Returns n_samples from a uniform random distribution between % 0.0 and 1.0. %\notes % urand() is the random generator available in ISIS and used here if % it is available (along with seed_random().) % If it is not available from ISIS then it and seed_random() are created % from the gsl ran_flat, rng_set, and rng_alloc routines if gsl is available. %\seealso{ } %!%- }}} { variable nrand = NULL; if (_NARGS == 1) nrand = (); if (nrand == NULL) nrand=1; return ran_flat(ran_gen,0.0,1.0,nrand); } define seed_random (seed) { rng_set (ran_gen, seed) ; } % }}} #endif #endif #ifexists urand %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_mc_cells () % {{{ %!%+{{{ %\function{v3d__mc_cells} %\synopsis{Generate Monte-Carlo random cells from a v3d array} %\usage{ (ixs, iys, izs) = v3d_mc_cells (v3darr, nevts); } %\description % Selects nevts number of random cells from the v3d array % and returns arrays of their indices: ix, iy, iz . % The probability of a cell being selected is proportional % to the cell's scalar value in the v3d array. % The v3darr is assumed to be non-negative. %\notes % Requires the urand routine. %\seealso{ v3d_mc_points, urand} %!%- }}} { if (_NARGS != 2) { message(" returns the indices of nevts randomly selected cells"); message(" with the prob of cell selection proportional its v3darr value."); usage("................. (ixs, iys, izs) = v3d_mc_cells (v3darr, nevts);"); } variable nevts = (); variable v3darr = (); variable mcixs = Long_Type [nevts]; variable mciys = Long_Type [nevts]; variable mcizs = Long_Type [nevts]; variable ix, iy, iz, iru, rus, irsort; variable dbarr = Double_Type [2*v3d_na+1,2*v3d_na+1,2*v3d_na+1]; variable dbarr_yz = Double_Type [2*v3d_na+1,2*v3d_na+1]; variable dbcum; % set the random seed value from the time... % eval( "seed_random("+str_delete_chars("9"+substr(time(),12,8),":")+");" ); seed_random(_time); % generate nevts random values: rus = urand(nevts); % and put them in sorted order - saving this to un-sort them in the end irsort = array_sort(rus); % and include an additional one above 1.0: % (so we can point to the "next one" without exceeding the array length) rus = [rus[irsort],[2.0]]; % Use Double for the cummulative sum to accomodate large dynamic % range and large number of cells. dbarr = double(v3darr); % Scale the random numbers to the cum-sum range: rus = sum(dbarr) * rus; % For each random number find where in the cum-sum it % appears/is_exceeded and return that cell location. % Start with the first (smallest) one: iru = 0L; % To speed up, get the dbarr values for the given ix,iy: variable dbarrxy; % Go through the whole array creating % the cummulative value as we go % and assigning cell indices as appropriate. dbcum=0.0; _for ix (0, 2*v3d_na, 1) { _for iy (0, 2*v3d_na, 1) { dbarrxy = dbarr[ix,iy,*]; _for iz (0, 2*v3d_na, 1) { % updated cum-sum %%dbcum = dbcum + dbarr[ix,iy,iz]; dbcum = dbcum + dbarrxy[iz]; % Faster this way % if dbcum has just become larger than the current rus % then that rus gets assigned to this ix,iy,iz % and check the next rus(s) in case they are in there too... while ( dbcum > rus[iru] ) { % save this cell location mcixs[iru] = ix; mciys[iru] = iy; mcizs[iru] = iz; % go to next random scaled value % note that it will get stuck at the last one and % not find a dbcum that exceeds it... (the 2.0 added above) iru = iru + 1; } } % end of iz loop } } % lastly, put these values back in the original random % order - the order that re-sorts the sorted indices... irsort = array_sort(irsort); return mcixs[irsort], mciys[irsort], mcizs[irsort]; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_mc_points () % {{{ %!%+{{{ %\function{v3d_mc_points} %\synopsis{Generate Monte-Carlo random points from a v3d array} %\usage{ (xs, ys, zs) = v3d_mc_points (v3darr, nevts); } %\description % Generate nevts number of Monte-Carlo x,y,z points with the % probability of the point proportional to the local value % of the v3darr cell. % The points are uniformly spread within the cells. % The v3darr is assumed to be non-negative. %\notes % Requires the urand routine. % % ( Note that optional rotations could be added to allow "viewing" the % object at some other angle besides the nominal (observer located at large +Z), % but not implemented for the moment and for overall modeling simplicity of having % a pre-defined viewing direction. ) %\seealso{ v3d_mc_cells, urand} %!%- }}} { if (_NARGS != 2) { message(" returns coordinates of nevts randomly selected points"); message(" with prob of point proportional to the local v3darr value."); usage("................. (xs, ys, zs) = v3d_mc_points (v3darr, nevts);"); } variable nevts = (); variable v3darr = (); variable mcixs = Long_Type [nevts]; variable mciys = Long_Type [nevts]; variable mcizs = Long_Type [nevts]; variable mcxs = Float_Type [nevts]; variable mcys = Float_Type [nevts]; variable mczs = Float_Type [nevts]; % get the random cells (mcixs, mciys, mcizs) = v3d_mc_cells(v3darr, nevts); % convert indices to randomized locations mcxs = v3d_va[mcixs] + (v3d_cube_radius/v3d_na)*(urand(nevts)-0.5); mcys = v3d_va[mciys] + (v3d_cube_radius/v3d_na)*(urand(nevts)-0.5); mczs = v3d_va[mcizs] + (v3d_cube_radius/v3d_na)*(urand(nevts)-0.5); return mcxs, mcys, mczs; } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #endif % - - - - - - - - - end of MC routines using urand - - - - - - % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% public define v3d_doppler_points () % {{{ %!%+{{{ %\function{v3d_doppler_points} %\synopsis{Calculate the Doppler energy correction given velocity properties} %\usage{ en_factor = v3d_doppler_points (x,y,z, vsxyz [, velval [, vtype [, oxyz [, rthphi]]]]); } %\description % This routine takes in values of x,y,z of events % and returns values of en_factor where en_factor is the energy % Doppler factor based on the line-of-sight velocity of % the points. % The velocity at x,y,z, is a sum of terms: % % - systemic velocity: vx, vy, vz [km/s] % % - radial or rigid-body rotation: velval, vtype, oxyz, rthphi % with vtype specifying the velocity properties: % -1 no additional velocity term % 0 v = velval [km/s] * r_hat with center at oxyz % 1 v = velval [km/s per unit] * r_vector, center at oxyz % 2 v = rigid body rotation: % |v| = velval [km/s per unit] * r_cylinder % with rotation axis through oxyz in direction rthphi. % 3 v = Keplerian: v_hat as in 2 but |v| ~ 1/sqrt(r) % % and for vtype = 1 or 2: set velval = v(ro)/ro % = 3: set velval = v(ro)*sqrt(ro) %\notes % The special relativistic ENERGY correction factor is given by: % en = lambda_0/lambda_los = 1/ (1+red_shift) = % % = 1 / ( gamma * (1 + beta cos(alpha)) ) % gamma = 1/sqrt(1-beta^2) % beta = v/c % beta cos(alpha) = v dot los_hat / c %\seealso{v3d_mc_points} %!%- }}} { variable x,y,z, vsxyz=NULL, velval=NULL, vtype=NULL, oxyz=NULL, rthphi=NULL; variable def_rthphi = [0.0, 90.0]; switch (_NARGS) { case 4: vsxyz = (); z = (); y = (); x = (); } { case 5: velval = (); vsxyz = (); z = (); y = (); x = (); } { case 6: vtype = (); velval = (); vsxyz = (); z = (); y = (); x = (); } { case 7: oxyz = (); vtype = (); velval = (); vsxyz = (); z = (); y = (); x = (); } { case 8: rthphi = (); oxyz = (); vtype = (); velval = (); vsxyz = (); z = (); y = (); x = (); } { % otherwise: message(" vsxyz=[vsx,vsy,vsz] in km/s"); message(" velval in k/ms [per unit]"); message(" vtype = -1 [no additional v],"); message(" 0 [velval*r_hat], "); message(" 1 [velval*r_vector],"); message(" 2 [rigid body rot about rthphi-through-oxyz]"); message(" For 1,2 set velval = v(ro)/ro "); message(" 3 v = Keplerian: v_hat as in 2 but |v| ~ 1/sqrt(r)"); message(" set velval = v(ro)*sqrt(ro) "); usage("................. en_factor = v3d_doppler_points(x,y,z, vsxyz [, " + "velval [, vtype [, oxyz [, rthphi]]]]);"); } % allow some to be missing and provide defaults if ( vsxyz == NULL ) vsxyz = [0.,0.,0.]; % velval being present signals inclusion of additional component % otherwise set vtype to -1 if ( velval == NULL ) vtype = -1; if ( vtype == NULL ) vtype = 1; if ( oxyz == NULL ) oxyz = [0.,0.,0.]; if ( rthphi == NULL ) rthphi = def_rthphi; % If it is vtype 3 then change to 2 and remember to use % a Keplerian v ~ 1/sqrt(r) instead of v ~ r . variable isKepler = 0; if (vtype == 3) { vtype = 2; isKepler = 1; } % we'll calculate (beta,) gamma, and betacosalpha for each point variable vtotal = Float_Type [length(x),3]; variable beta, oneogamma, betacosalpha, los_hat, ip; variable r_vec, r_vec_x, r_vec_y, r_vec_z, r_norm; % calculate the total velocity vector for each point % depending on the vtype % RFE - v3d_doppler_points - Need to do relativistic % addition of the vsxyz with the field velocities... if (vtype != -1) { switch (vtype) { case 0: % 0 - [velval*r_hat] r_vec_x = x - oxyz[0]; r_vec_y = y - oxyz[1]; r_vec_z = z - oxyz[2]; r_norm = sqrt( r_vec_x^2 + r_vec_y^2 + r_vec_z^2 ); vtotal[*,0] = vsxyz[0] + velval * r_vec_x/r_norm; vtotal[*,1] = vsxyz[1] + velval * r_vec_y/r_norm; vtotal[*,2] = vsxyz[2] + velval * r_vec_z/r_norm; } { case 1: % 1 - [velval*r_vector]"); % vtotal[*,0] = vsxyz[0] + velval * (x - oxyz[0]); vtotal[*,1] = vsxyz[1] + velval * (y - oxyz[1]); vtotal[*,2] = vsxyz[2] + velval * (z - oxyz[2]); } { case 2: % 2 - [rigid body rot about rthphi-through-oxyz] % v_hat = axis_hat (x) r_hat % 3 - v = Keplerian: v_hat as in 2 but |v| ~ 1/sqrt(r) % Unit vector in axial direction is: variable axuvx, axuvy, axuvz; axuvx = cos(rthphi[1]*DTOR)*cos(rthphi[0]*DTOR); axuvy = sin(rthphi[1]*DTOR); axuvz = cos(rthphi[1]*DTOR)*sin(-1.0*rthphi[0]*DTOR); variable r_len2, axial, radial; variable r_hat_x, r_hat_y, r_hat_z; variable v_hat_x, v_hat_y, v_hat_z; r_vec_x = x - oxyz[0]; r_vec_y = y - oxyz[1]; r_vec_z = z - oxyz[2]; r_len2 = r_vec_x^2 + r_vec_y^2 + r_vec_z^2; r_norm = sqrt(r_len2); % r_hat is the unit vector from O to P r_hat_x = r_vec_x/r_norm; r_hat_y = r_vec_y/r_norm; r_hat_z = r_vec_z/r_norm; % The velocity unit vector is v_hat = axuv _cross_ r_hat % x component v_hat_x = axuvy * r_hat_z - axuvz * r_hat_y; % y component v_hat_y = axuvz * r_hat_x - axuvx * r_hat_z; % z component v_hat_z = axuvx * r_hat_y - axuvy * r_hat_x; % The speed is proportional to the distance from the axis % Calculate the axial and radial coordinates of the points axial = r_vec_x*axuvx + r_vec_y*axuvy + r_vec_z*axuvz; radial = ( r_len2 - axial^2 ); radial[where(radial < 0.0)] = 0.0; radial = sqrt(radial); if (isKepler == 0) { % |v| ~ r vtotal[*,0] = vsxyz[0] + velval * radial * v_hat_x; vtotal[*,1] = vsxyz[1] + velval * radial * v_hat_y; vtotal[*,2] = vsxyz[2] + velval * radial * v_hat_z; } else { % |v| ~ 1/sqrt(r) vtotal[*,0] = vsxyz[0] + (velval/sqrt(radial)) * v_hat_x; vtotal[*,1] = vsxyz[1] + (velval/sqrt(radial)) * v_hat_y; vtotal[*,2] = vsxyz[2] + (velval/sqrt(radial)) * v_hat_z; } } { % unrecognized vtype message(" *** v3d_doppler_points: unrecognized vtype ="+string(vtype)); } } else { % vtype = -1 - just the vsxyz vtotal[*,0] = vsxyz[0]; vtotal[*,1] = vsxyz[1]; vtotal[*,2] = vsxyz[2]; } % The line of sight is assumed to be in the -z direction los_hat = [0.,0.,-1.0]; % one-over-gamma oneogamma = sqrt(1.0 - ( ( vtotal[*,0]^2 + vtotal[*,1]^2 + vtotal[*,2]^2 ) / Const_c_kms^2 ) ); betacosalpha = ( los_hat[0]*vtotal[*,0] + los_hat[1]*vtotal[*,1] + los_hat[2]*vtotal[*,2] ) / Const_c_kms; return, oneogamma / (1 + betacosalpha); } % }}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % If we made it this far, then: #ifexists provide provide ("v3d"); #endif %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % that's all folks. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%