From 080e0a9ee08d77bfa19aa85b0f33f0ba150fca54 Mon Sep 17 00:00:00 2001 From: Andreas Irrgang Date: Mon, 25 Jan 2021 10:45:22 +0100 Subject: [PATCH 1/2] Make the function 'cel2gal' return _NaNs instead of throwing an error when negative distances are provided --- src/slang/astronomy/cel2gal.sl | 40 +++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/src/slang/astronomy/cel2gal.sl b/src/slang/astronomy/cel2gal.sl index 3f1db5f0..ff0af8e7 100644 --- a/src/slang/astronomy/cel2gal.sl +++ b/src/slang/astronomy/cel2gal.sl @@ -1,12 +1,12 @@ define cel2gal() %{{{ %!%+ %\function{cel2gal} -%\synopsis{Transform celestial to cartesian Galactic coordinates} +%\synopsis{Transform celestial to Cartesian Galactic coordinates} %\usage{cel2gal(Double_Types ah[], am[], as[], dd[], dm[], ds[], dist[], vrad[], pma_cos_d[], pmd[]; qualifiers)} %\description % Transform celestial coordinates (right ascension [h, m, s], declination [deg, arcmin, arcsec], % distance [kpc], radial velocity [km/s], proper motion in right ascension times cosine of declination -% [mas/yr], proper motion in declination [mas/yr]) to right-handed, cartesian Galactic coordinates +% [mas/yr], proper motion in declination [mas/yr]) to right-handed, Cartesian Galactic coordinates % (x [kpc], y [kpc], z [kpc], vx [km/s], vy [km/s], vz [km/s]) with the Galactic center at the origin, % the Sun on the negative x-axis and the z-axis pointing to the north Galactic pole implying clockwise % Galactic rotation when seen from the half space with positive z. See function 'RD2rad' for detailed @@ -57,8 +57,12 @@ define cel2gal() %{{{ ah = [1.*ah]; am = [1.*am]; as = [1.*as]; dd = [1.*dd]; dm = [1.*dm]; ds = [1.*ds]; dist = [1.*dist]; rv = [1.*rv]; pma = [1.*pma]; pmd = [1.*pmd]; if(qualifier_exists("parallax")) dist = 1./dist; % conversion from parallax in units of mas to distance in kpc - if(length(where(dist<0))>0) - throw UsageError, sprintf("Usage error in '%s': negative distances are not possible.", _function_name()); + variable ind = where(dist<0); + if(length(ind)>0) + { + vmessage("Warning in '%s': Replaced %d negative distance(s) by _NaN(s).", _function_name(), length(ind)); + dist[ind] = _NaN; + } variable SunGCDist = qualifier("SunGCDist", 8.40); % Model I in Irrgang et al., 2013, A&A, 549, A137 variable vxs = qualifier("vxs", 11.1); % from Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829-1833 variable vys = qualifier("vys", 12.24); % from Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829-1833 @@ -67,7 +71,7 @@ define cel2gal() %{{{ % convert right ascension and declination to radians: variable ra, decl; (ra, decl) = RD2rad(ah, am, as, dd, dm, ds); - % convert polar coordinates to cartesian celestial coordinates: + % convert polar coordinates to Cartesian celestial coordinates: variable cartCC = Array_Type[3]; % [x, y, z] cartCC[0] = dist*cos(ra)*cos(decl); cartCC[1] = dist*sin(ra)*cos(decl); @@ -111,7 +115,7 @@ define cel2gal() %{{{ rotMat[2,1] = v_NGP[1]; rotMat[2,2] = v_NGP[2]; } - % multiply the Celestial cartesian data with the rotation matrix + % multiply the Celestial Cartesian data with the rotation matrix % to get the coordinates in Galactic system and add shift for the % sun-gc distance: variable x, y, z; @@ -127,17 +131,17 @@ define cel2gal() %{{{ variable mas_yr_to_rad_s = 1.536282790e-16; % conversion factor from mas/yr to radians/s = 2.*PI/(360.*60.*60.)/1000./(365.24977652625375*24.*60.*60.); pmd *= mas_yr_to_rad_s; pma *= mas_yr_to_rad_s; - % convert radial velocity to a cartesian vector in Celestial coordinates + % convert radial velocity to a Cartesian vector in Celestial coordinates variable radvCC = Array_Type[3]; % [vx, vy, vz] radvCC[0] = rv*cos(ra)*cos(decl); radvCC[1] = rv*sin(ra)*cos(decl); radvCC[2] = rv*sin(decl); - % convert angular declination velocity to a cartesian vector in Celestial coordinates + % convert angular declination velocity to a Cartesian vector in Celestial coordinates variable dvCC = Array_Type[3]; % [vx, vy, vz] dvCC[0] = -pmd*dist*cos(ra)*sin(decl); dvCC[1] = -pmd*dist*sin(ra)*sin(decl); dvCC[2] = pmd*dist*cos(decl); - % convert angular right ascention velocity to a cartesian vector in Celestial coordinates + % convert angular right ascention velocity to a Cartesian vector in Celestial coordinates variable avCC = Array_Type[3]; %[vx, vy, vz] avCC[0] = -pma*dist*sin(ra); % factor cos(decl) already in pma avCC[1] = pma*dist*cos(ra); % factor cos(decl) already in pma @@ -147,7 +151,7 @@ define cel2gal() %{{{ totalvCC[0] = radvCC[0] + dvCC[0] + avCC[0]; totalvCC[1] = radvCC[1] + dvCC[1] + avCC[1]; totalvCC[2] = radvCC[2] + dvCC[2] + avCC[2]; - % rotate the vector by the rotation matrix, and add the motions of the Sun within the local standard of rest (vxs, vys, vzss), + % rotate the vector by the rotation matrix, and add the motions of the Sun within the local standard of rest (vxs, vys, vzs), % and the motion of the local standard of rest itself (vlsr): variable vx, vy, vz; vx = rotMat[0,0]*totalvCC[0] + rotMat[0,1]*totalvCC[1] + rotMat[0,2]*totalvCC[2] + vxs; @@ -161,10 +165,10 @@ define cel2gal() %{{{ define gal2cel() %{{{ %!%+ %\function{gal2cel} -%\synopsis{Transform cartesian Galactic to celestial coordinates} +%\synopsis{Transform Cartesian Galactic to celestial coordinates} %\usage{gal2cel(Double_Types x[], y[], z[], vx[], vy[], vz[]; qualifiers)} %\description -% Transform right-handed, cartesian Galactic coordinates (x [kpc], y [kpc], z [kpc], vx [km/s], +% Transform right-handed, Cartesian Galactic coordinates (x [kpc], y [kpc], z [kpc], vx [km/s], % vy [km/s], vz [km/s]) with the Galactic center at the origin, the Sun on the negative x-axis % and the z-axis pointing to the north Galactic pole implying clockwise Galactic rotation when % seen from the half space with positive z to celestial coordinates (right ascension [h, m, s], @@ -215,12 +219,12 @@ define gal2cel() %{{{ variable vys = qualifier("vys", 12.24); % from Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829-1833 variable vzs = qualifier("vzs", 7.25); % from Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829-1833 variable vlsr = qualifier("vlsr", 242.0); % see Model I in Irrgang et al., 2013, A&A, 549, A137 - % shift cartesian coordinate system from Galactic center to Sun: + % shift Cartesian coordinate system from Galactic center to Sun: x += SunGCDist; vx -= vxs; vy -= (vys + vlsr); vz -= vzs; - % rotation matrix from Galactically aligned to celestially aligned cartesian coordinates + % rotation matrix from Galactically aligned to celestially aligned Cartesian coordinates % see below ("if(qualifier_exists("GC_NGP")==1)") for a derivation variable rotMat = Double_Type[3,3]; % initialize rotMat with matrix derived from coordinates for Galactic center (17 45 37.224, -28 56 10.23) and north Galactic pole (12 51 26.282, 27 07 42.01) @@ -259,12 +263,12 @@ define gal2cel() %{{{ rotMat[1,2] = v_NGP[1]; rotMat[2,2] = v_NGP[2]; } - % rotate cartesian coordiante system so that it is aligned to right ascension and declination: - variable cartCC = Array_Type[3]; % cartesian celestial coordinates + % rotate Cartesian coordiante system so that it is aligned to right ascension and declination: + variable cartCC = Array_Type[3]; % Cartesian celestial coordinates cartCC[0] = rotMat[0,0]*x + rotMat[0,1]*y + rotMat[0,2]*z; cartCC[1] = rotMat[1,0]*x + rotMat[1,1]*y + rotMat[1,2]*z; cartCC[2] = rotMat[2,0]*x + rotMat[2,1]*y + rotMat[2,2]*z; - variable cartCV = Array_Type[3]; % cartesian celestial velocities + variable cartCV = Array_Type[3]; % Cartesian celestial velocities cartCV[0] = rotMat[0,0]*vx + rotMat[0,1]*vy + rotMat[0,2]*vz; cartCV[1] = rotMat[1,0]*vx + rotMat[1,1]*vy + rotMat[1,2]*vz; cartCV[2] = rotMat[2,0]*vx + rotMat[2,1]*vy + rotMat[2,2]*vz; @@ -273,7 +277,7 @@ define gal2cel() %{{{ cartCC[1] *= kpc_to_km; cartCC[2] *= kpc_to_km; variable dis, ra, theta, vrad, pma, pmd; - % transform cartesian to spherical coordinates: + % transform Cartesian to spherical coordinates: (dis, ra, theta, vrad, pma, pmd) = cart2sphere(cartCC[0], cartCC[1], cartCC[2], cartCV[0], cartCV[1], cartCV[2]); variable decl = -theta + PI/2.; % due to different definitions of azimuth angle theta and declination delta pmd *= -1.; % due to different definitions of azimuth angle theta and declination delta -- GitLab From 945169c390990db268fac3c0e832259486f86564 Mon Sep 17 00:00:00 2001 From: Andreas Irrgang Date: Mon, 25 Jan 2021 10:45:48 +0100 Subject: [PATCH 2/2] Revise some qualifiers of the function 'orbit_calculator' - The qualifiers 'r_disk' and 'N_disk_crossings' are lumped together into a qualifier called 'disk', which is a structure that allows to define more generic disk geometries. - The default value of the qualifier 'MC_runs' is changed from 10^4 to 10^5. - The default value of the qualifier 'tolerance' is changed from 10^-10 to 10^-8. --- src/misc/astronomy/orbit_calculator.sl | 132 ++++++++++++++----------- 1 file changed, 73 insertions(+), 59 deletions(-) diff --git a/src/misc/astronomy/orbit_calculator.sl b/src/misc/astronomy/orbit_calculator.sl index 34e69ce2..bc8af4af 100644 --- a/src/misc/astronomy/orbit_calculator.sl +++ b/src/misc/astronomy/orbit_calculator.sl @@ -773,9 +773,9 @@ define orbit_calculator() %{{{ % the initial and final situation is of interest, but rather the entire trajectories, use % the qualifier 'set'. In this case, the additional field "tr" (trajectory) (again a structure, % one field per orbit ("o1", "o2", ...)) is added to the returned structure containing all -% orbits. For tracing back orbits to the Galactic disk, use the qualifiers 'r_disk' and -% 'N_disk_crossings' to stop individual calculations at the moment a trajectory crosses -% the xy-plane inside a circle of radius 'r_disk' for the N_disk_crossings-th time. +% orbits. For tracing back orbits to the Galactic disk, use the 'disk' qualifier to stop +% individual calculations at the moment a trajectory has crossed the disk for a certain +% number of times. %\qualifiers{ %\qualifier{coords}{ [\code{="cyl"}]: Use cylindrical ("cyl") or cartesian ("cart") coordinates for the internal % computations. Note that circular orbits like that of the Sun are computed faster in cylindrical @@ -785,15 +785,16 @@ define orbit_calculator() %{{{ %\qualifier{dt}{: Deactivate adaptive mode and use this fixed stepsize instead.} %\qualifier{model}{ [\code{="AS"}]: Function ("AS", "MN_NFW", "MN_TF", or "plummer_MW") evaluating the equations % of motion, circular velocity, and energy of the model potential.} -%\qualifier{MC_runs}{ [\code{=10000}]: Number of Monte Carlo realizations in the case that 1-sigma errors are given.} +%\qualifier{MC_runs}{ [\code{=nint(10^5)}]: Number of Monte Carlo realizations in the case that 1-sigma errors are given.} %\qualifier{ODE_solver}{ [\code{="RKCK"}]: Choose among three different Runge-Kutta (RK) integration methods: % "RKF": RK-Fehlberg, "RKCK": RK-Cash-Karp, "RKDP": RK-Dormand-Prince.} %\qualifier{parallax_pma_corr}{ [\code{=0}]: Correlation between parallax and proper motion in right ascension.} %\qualifier{parallax_pmd_corr}{ [\code{=0}]: Correlation between parallax and proper motion in declination.} %\qualifier{pma_pmd_corr}{ [\code{=0}]: Correlation between proper motion in right ascension and in declination.} -%\qualifier{r_disk and N_disk_crossings}{ [\code{=0}, \code{=1}]: If r_disk is not zero, integration will be stopped -% at the moment when the xy-plane is crossed inside a circle of radius r_disk [kpc] for the -% N_disk_crossings-th time.} +%\qualifier{disk}{ [\code{=struct{radius = 0, height = 0.1, x = 0, y = 0, z = 0, crossings = 1}}]: If disk.radius > 0, +% orbit integration will be stopped at the moment a trajectory has crossed a horizontal plane +% located at z = disk.z + disk.height * Gaussian-random-number inside a circle of radius +% disk.radius centered at (disk.x, disk.y) for the disk.crossings-th time.} %\qualifier{seed}{ [\code{=_time()}]: Seed the random number generator via the function 'seed_random'.} %\qualifier{set}{ [\code{=0}]: If present, trajectories will be saved ("Save entire trajectories"). As % this can be very memory-consuming for a large number of orbits, one can additionally @@ -802,15 +803,15 @@ define orbit_calculator() %{{{ %\qualifier{stff}{: "Save to fits files": Prefix of fits files to which initial and final structures are written.} %\qualifier{SunGCDist}{: By default, the Sun-Galactic center distance is taken from the current model. % Use this qualifier to explicitly set a distance in kpc.} -%\qualifier{tolerance}{ [\code{=1e-10}]: Absolut error control tolerance; lower limit: 1e-15.} +%\qualifier{tolerance}{ [\code{=1e-8}]: Absolute error control tolerance; lower limit: 1e-15.} %\qualifier{verbose}{: Show intermediate times t.} %\qualifier{Any qualifiers from 'cel2gal' and the model potential function except 'vlsr' and 'eomecd'.}{ % Important note: for consistency, the local standard of rest velocity vlsr is calculated from the % circular velocity of the current model potential evaluated at (r=SunGCDist, z=0).} %} %\example -% s = orbit_calculator(12,22,29.6,40,49,36,3.078,262,-13.52,16.34,-1000; r_disk=50, N_disk_crossings=1); -% s = orbit_calculator(12,22,29.6,40,49,36,3.078,[0.6,0.3],262,5,-13.52,1.31,16.34,1.37,-1000; r_disk=50, stff="HIP60350", MC_runs=100); +% s = orbit_calculator(12,22,29.6,40,49,36,3.078,262,-13.52,16.34,-1000; disk = struct{radius=50, height=0.2, x=0, y=0, z=0, crossings=1}); +% s = orbit_calculator(12,22,29.6,40,49,36,3.078,[0.6,0.3],262,5,-13.52,1.31,16.34,1.37,-1000; disk = struct{radius=50, height=0.2, x=0, y=0, z=0, crossings=1}, stff="HIP60350", MC_runs=100); % s = orbit_calculator(-8.4,0,0,0,242,0,250; set, model="MN_TF"); % plot(s.tr.o0.x, s.tr.o0.y); % s = orbit_calculator([8:10:#500], [0:1:#500], [4:5:#500], [0:0:#500], [210:230:#500], [-10:-50:#500], 1000); @@ -839,11 +840,25 @@ define orbit_calculator() %{{{ } model = __get_reference(model); % function evaluating the equations of motions, circular velocity, energy of the model potential variable dt_fixed = qualifier_exists("dt"); % fixed stepsize for integrator - variable eps = _max( qualifier("tolerance", 1e-10), 1e-15 ); % error control tolerance - variable r_disk = qualifier("r_disk", 0); % If unequal zero, integration is stopped at the moment when the xy-plane is crossed inside this radius for the $N_disk_crossings time - variable N_disk_crossings = qualifier("N_disk_crossings", 1); % Integration is stopped when the xy-plane has been crossed inside r_disk for the $N_disk_crossings time - if(typeof(N_disk_crossings)!=Integer_Type) - throw UsageError, sprintf("Usage error in '%s': qualifier 'N_disk_crossings' has to be an Integer_Type.", _function_name()); + variable eps = _max( qualifier("tolerance", 1e-8), 1e-15 ); % error control tolerance + variable disk = qualifier("disk", struct{radius = 0, height = 0.1, x = 0, y = 0, z = 0, crossings = 1}); + % If disk.radius > 0, orbit integration will be stopped at the moment a trajectory has crossed a horizontal plane located at + % z = disk.z + disk.height * Gaussian-random-number inside a circle of radius disk.radius centered at (disk.x, disk.y) for the disk.crossings-th time. + if(typeof(disk)!=Struct_Type) + { + _pop_n(_NARGS); % remove function arguments from stack + throw UsageError, sprintf("Usage error in '%s': The qualifier 'disk' has to be a Struct_Type.", _function_name()); + } + if(struct_field_exists(disk, "radius")==0 or struct_field_exists(disk, "height")==0 or struct_field_exists(disk, "x")==0 or struct_field_exists(disk, "y")==0 or struct_field_exists(disk, "z")==0 or struct_field_exists(disk, "crossings")==0) + { + _pop_n(_NARGS); % remove function arguments from stack + throw UsageError, sprintf("Usage error in '%s': The qualifier 'disk' has to be a Struct_Type with fields 'radius', 'height', 'x', 'y', 'z', and 'crossings'.", _function_name()); + } + if(typeof(disk.crossings)!=Integer_Type || disk.crossings<1) + { + _pop_n(_NARGS); % remove function arguments from stack + throw UsageError, sprintf("Usage error in '%s': The qualifier 'disk.crossings' has to be a positive Integer_Type.", _function_name()); + } variable set = qualifier_exists("set"); % "save entire trajectory" variable set_threshold = qualifier("set", 0); if(typeof(set_threshold)==Null_Type) % to cover cases where the qualifier 'set' is present but without a specified value @@ -860,8 +875,7 @@ define orbit_calculator() %{{{ if(qualifier_exists("SunGCDist")) SunGCDist = qualifier("SunGCDist"); set_struct_field(qualies, "SunGCDist", SunGCDist); % qualifier 'SunGCDist' in function cel2gal is set according to the current potential - variable temp0 = [SunGCDist,0,0,0,0,0]; - reshape(temp0, [6,1]); + variable temp0 = 0*Double_Type[6,length(SunGCDist)]; temp0[0,*] = SunGCDist; set_struct_field(qualies, "eomecd", "circ"); % important for the following line set_struct_field(qualies, "vlsr", ((@model)(0, temp0;; qualies))[0] ); % qualifier 'vlsr' in function cel2gal is set according to the current potential set_struct_field(qualies, "eomecd", "eom"); % set back to default @@ -869,7 +883,7 @@ define orbit_calculator() %{{{ set_struct_field(qualies, "coords", coords); % ============================================================================ if(qualifier_exists("stff")) % some keywords for fits-header to know what input parameters or qualifiers were used for this computation - variable keys = struct{SunGCDist = string(SunGCDist)}; + variable keys = struct{SunGCDist = string(median(SunGCDist))}; variable x, y, z, vx, vy, vz, t_end; % t_end: maximum integration time in Myrs if(_NARGS==11) @@ -884,7 +898,7 @@ define orbit_calculator() %{{{ { variable dist_sigma, rv_sigma, pma_sigma, pmd_sigma; (ah, am, as, dd, dm, ds, dist, dist_sigma, rv, rv_sigma, pma, pma_sigma, pmd, pmd_sigma, t_end) = (); - variable MC_runs = nint(qualifier("MC_runs", 10000)); + variable MC_runs = nint(qualifier("MC_runs", nint(10^5))); variable parallax_pma_corr = qualifier("parallax_pma_corr", 0); if(abs(parallax_pma_corr)>1) throw UsageError, sprintf("Usage error in '%s': Qualifier 'parallax_pma_corr' has to be between -1 and 1.", _function_name()); @@ -896,13 +910,13 @@ define orbit_calculator() %{{{ throw UsageError, sprintf("Usage error in '%s': Qualifier 'pma_pmd_corr' has to be between -1 and 1.", _function_name()); variable seed = nint(qualifier("seed", _time())); seed_random(seed); % initialize the seed according to qualifier or time - variable rv_distr, gran_numbers; + variable rv_distr; if(length(rv_sigma)==2) { - gran_numbers = grand(MC_runs); - gran_numbers[where(gran_numbers>0)] *= rv_sigma[0]; - gran_numbers[where(gran_numbers<0)] *= rv_sigma[1]; - rv_distr = rv+gran_numbers; + rv_distr = grand(MC_runs); + rv_distr[where(rv_distr>0)] *= rv_sigma[0]; + rv_distr[where(rv_distr<0)] *= rv_sigma[1]; + rv_distr += rv; } else rv_distr = rv+grand(MC_runs)*rv_sigma[0]; @@ -916,10 +930,9 @@ define orbit_calculator() %{{{ cov_mat[0,0] = 1; cov_mat[1,0] = 0; cov_mat[0,1] = 0; cov_mat[2,0] = 0; cov_mat[0,2] = 0; if(len==2) { - gran_numbers = grand(MC_runs); - gran_numbers[where(gran_numbers>0)] *= dist_sigma[0]; - gran_numbers[where(gran_numbers<0)] *= dist_sigma[1]; - distr[0,*] = gran_numbers; + distr[0,*] = grand(MC_runs); + distr[0,where(distr[0,*]>0)] *= dist_sigma[0]; + distr[0,where(distr[0,*]<0)] *= dist_sigma[1]; } else distr[0,*] = grand(MC_runs)*dist_sigma[0]; @@ -938,10 +951,9 @@ define orbit_calculator() %{{{ cov_mat[1,1] = 1; cov_mat[1,0] = 0; cov_mat[0,1] = 0; cov_mat[1,2] = 0; cov_mat[2,1] = 0; if(len==2) { - gran_numbers = grand(MC_runs); - gran_numbers[where(gran_numbers>0)] *= pma_sigma[0]; - gran_numbers[where(gran_numbers<0)] *= pma_sigma[1]; - distr[1,*] = gran_numbers; + distr[1,*] = grand(MC_runs); + distr[1,where(distr[1,*]>0)] *= pma_sigma[0]; + distr[1,where(distr[1,*]<0)] *= pma_sigma[1]; } else distr[1,*] = grand(MC_runs)*pma_sigma[0]; @@ -960,10 +972,9 @@ define orbit_calculator() %{{{ cov_mat[2,2] = 1; cov_mat[2,0] = 0; cov_mat[0,2] = 0; cov_mat[1,2] = 0; cov_mat[2,1] = 0; if(len==2) { - gran_numbers = grand(MC_runs); - gran_numbers[where(gran_numbers>0)] *= pmd_sigma[0]; - gran_numbers[where(gran_numbers<0)] *= pmd_sigma[1]; - distr[2,*] = gran_numbers; + distr[2,*] = grand(MC_runs); + distr[2,where(distr[2,*]>0)] *= pmd_sigma[0]; + distr[2,where(distr[2,*]<0)] *= pmd_sigma[1]; } else distr[2,*] = grand(MC_runs)*pmd_sigma[0]; @@ -987,16 +998,17 @@ define orbit_calculator() %{{{ variable dist_distr = distr[0,*] + dist; variable pma_distr = distr[1,*] + pma; variable pmd_distr = distr[2,*] + pmd; - variable ind = where(dist_distr>0); % omit negative distances or parallaxes - len = length(ind); + (x, y, z, vx, vy, vz) = cel2gal(ah*ones(MC_runs), am*ones(MC_runs), as*ones(MC_runs), dd*ones(MC_runs), dm*ones(MC_runs), ds*ones(MC_runs), dist_distr, rv_distr, pma_distr, pmd_distr;; qualies); + variable ind = where(dist_distr>0); len = length(ind); % omit negative distances or parallaxes if(len!=MC_runs) + { vmessage("Warning in '%s': only %d out of %d Monte Carlo runs are performed because initial conditions with negative distances have been omitted.", _function_name(), len, MC_runs); - variable dummy = 1.+Double_Type[len]; - (x, y, z, vx, vy, vz) = cel2gal(ah*dummy, am*dummy, as*dummy, dd*dummy, dm*dummy, ds*dummy, dist_distr[ind], rv_distr[ind], pma_distr[ind], pmd_distr[ind];; qualies); + x = x[ind]; y = y[ind]; z = z[ind]; vx = vx[ind]; vy = vy[ind]; vz = vz[ind]; + } if(qualifier_exists("stff")) % some keywords for fits-header to know what input parameters or qualifiers were used for this computation { - keys = struct_combine( keys, struct{ rah=string(ah), ras=string(am), ram=string(as), decd=string(dd), decm=string(dm), decs=string(ds), dist=string(dist), - dist_sigma_plus=string(dist_sigma[0]), dist_sigma_minus=string(dist_sigma[-1]), rv=string(rv), rv_sigma_plus=string(rv_sigma[0]), rv_sigma_minus=string(rv_sigma[-1]), + keys = struct_combine( keys, struct{ maximum_integration_time = string(t_end), rah=string(ah), ras=string(am), ram=string(as), decd=string(dd), decm=string(dm), decs=string(ds), + dist=string(dist), dist_sigma_plus=string(dist_sigma[0]), dist_sigma_minus=string(dist_sigma[-1]), rv=string(rv), rv_sigma_plus=string(rv_sigma[0]), rv_sigma_minus=string(rv_sigma[-1]), pma_cos_delta=string(pma), pma_cos_delta_sigma_plus=string(pma_sigma[0]), pma_cos_delta_sigma_minus=string(pma_sigma[-1]), pmd=string(pmd), pmd_sigma_plus=string(pmd_sigma[0]), pmd_sigma_minus=string(pmd_sigma[-1]) }); } @@ -1037,6 +1049,7 @@ define orbit_calculator() %{{{ ret.i.vy = vy; ret.i.vz = vz; variable disk_crossings = 0*Integer_Type[len]; + variable disk_positions = disk.z + disk.height*grand(len); variable Rmin = [sqrt(x^2+y^2+z^2)]; % initialize Rmin as Array_Type variable Rmax = [sqrt(x^2+y^2+z^2)]; % initialize Rmax as Array_Type variable m = Double_Type[6, len]; @@ -1100,8 +1113,8 @@ define orbit_calculator() %{{{ variable t = 0; variable dt = 0.001; if(dt_fixed) dt = qualifier("dt"); % use fixed stepsize when qualifier "dt" is present - variable t_eoc = t_end + 0*Double_Type[len]; % t_eoc: time at the end or moment of crossing the xy-plane - ind = [0:len-1:1]; % array of indices, important for bookkeeping of which orbits already crossed the xy-plane and thus for qualifier 'r_disk' + variable t_eoc = t_end + 0*Double_Type[len]; % t_eoc: time at the end or moment of crossing the disk + ind = [0:len-1:1]; % array of indices, important for bookkeeping of which orbits already crossed the disk and thus for qualifier 'disk.radius' if(set) { variable t_last_save = 0; @@ -1150,7 +1163,7 @@ define orbit_calculator() %{{{ % variable scale = ( 0.5*eps*dt/max(abs(dx_c-dx_e)) )^0.25; % "derivation of this formula can be found in advanced books on numerical analysis" if(scale>=1 or dt<=1e-14*t) % scale>1: accept stepsize only if correction factor is equal/larger than 1 { % dt<=abs(1e-14*t) in order to avoid infinite loops due to Double_Type precision limit at dt=1e-16*t - if(r_disk!=0) + if(disk.radius>0) { variable m_last = @(m[*,ind]); variable t_last = @t; @@ -1169,23 +1182,25 @@ define orbit_calculator() %{{{ temp0 = where(temp_R0) { - % identify orbits that have crossed the xy-plane: + % identify orbits that have crossed the disk: variable ind_cro, ind_con; % indices of orbits which have crossed plane and of those for which calculation is continued - if(coords=="cart") - ind_cro = where(m[2,ind] * m_last[2,*] <= 0 and m[2,ind]!=m_last[2,*] and sqrt( (m[0,ind])^2 + (m[1,ind])^2) <= r_disk); % z * z_last <= 0 but z!=z_last!=0 and r<=r_disk -> orbit crossed xy-plane inside radius of r_disk during last step + if(coords=="cart") % (z-z_disk) * (z_last-z_disk) <= 0 but z!=z_last!=0 and r<=disk.radius -> orbit crossed disk inside radius of disk.radius during last step: + ind_cro = where((m[2,ind]-disk_positions[ind]) * (m_last[2,*]-disk_positions[ind]) <= 0 and m[2,ind]!=m_last[2,*] and + ( (m_last[0,*]-disk.x)^2 + (m_last[1,*]-disk.y)^2 <= disk.radius^2 or (m[0,ind]-disk.x)^2 + (m[1,ind]-disk.y)^2 <= disk.radius^2)); else % coords = "cyl" - ind_cro = where(m[2,ind] * m_last[2,*] <= 0 and m[2,ind]!=m_last[2,*] and m[0,ind] <= r_disk); % z * z_last <= 0 but z!=z_last!=0 and r<=r_disk -> orbit crossed xy-plane inside radius of r_disk during last step - disk_crossings[ind[ind_cro]] += 1; - ind_cro = where(disk_crossings[ind]==N_disk_crossings, &ind_con); + ind_cro = where((m[2,ind]-disk_positions[ind]) * (m_last[2,*]-disk_positions[ind]) <= 0 and m[2,ind]!=m_last[2,*] and ( (m_last[0,*]*cos(m_last[1,*])-disk.x)^2 + (m_last[0,*]*sin(m_last[1,*])-disk.y)^2 <= disk.radius^2 or + (m[0,ind]*cos(m[1,ind])-disk.x)^2 + (m[0,ind]*sin(m[1,ind])-disk.y)^2 <= disk.radius^2)); + disk_crossings[ind[ind_cro]] += 1; + ind_cro = where(disk_crossings[ind]==disk.crossings, &ind_con); variable len_ind_cro = length(ind_cro); if(len_ind_cro>0) { % ===================================================================== % interpolate to moment of crossing and save the result in m and t_end: - variable t_travel_last = -m_last[2,ind_cro]/m_last[5,ind_cro]; % t_travel_last := t_cross - t_last = -z_last/vz_last <=> t_travel_last*vz_last + z_last = 0 <= linear interpolation - variable t_travel = m[2,ind[ind_cro]]/m[5,ind[ind_cro]]; % t_travel = t - t_cross = z/vz <=> vz*t_travel = z <= linear interpolation + variable t_travel_last = (disk_positions[ind[ind_cro]]-m_last[2,ind_cro])/m_last[5,ind_cro]; % t_travel_last := t_cross - t_last = (disk_positions - z_last)/vz_last <=> t_travel_last*vz_last = disk_positions - z_last <= linear interpolation + variable t_travel = (m[2,ind[ind_cro]]-disk_positions[ind[ind_cro]])/m[5,ind[ind_cro]]; % t_travel = t - t_cross = (z-disk_positions)/vz <=> vz*t_travel = z - disk_positions <= linear interpolation % average with weights: the larger the estimated travel time, the worse the linear approximation variable weight = t_travel / ( t_travel + t_travel_last ); t_eoc[ind[ind_cro]] = weight*t_last + (1.-weight)*t; @@ -1237,11 +1252,11 @@ define orbit_calculator() %{{{ } } } - if(r_disk!=0) + if(disk.radius>0) { if(len_ind_cro>0) { - % continue orbit calculation only for those orbits in m that did not cross the plane yet: + % continue orbit calculation only for those orbits in m that did not cross the disk yet: ind = ind[ind_con]; } } @@ -1319,19 +1334,18 @@ define orbit_calculator() %{{{ if(qualifier_exists("stff")) % some keywords for fits-header to know what input parameters or qualifiers were used for this computation { - keys = struct_combine( keys, struct{ vxs, vys, vzs, Mb, Md, Mh, bb, ad, bd, ah, exponent, cutoff, model, method, eps, r_disk} ); + keys = struct_combine( keys, struct{ vxs, vys, vzs, Mb, Md, Mh, bb, ad, bd, ah, exponent, cutoff, model, method, eps}, disk ); variable field; foreach field (["vxs", "vys", "vzs", "Mb", "Md", "Mh", "bb", "ad", "bd", "ah", "exponent", "cutoff"]) { if(qualifier_exists(field)) - set_struct_field(keys, field, string(qualifier(field))); + set_struct_field(keys, field, string(median(qualifier(field)))); else set_struct_field(keys, field, "Default value"); } keys.model = string(model); keys.method = method; keys.eps = string(eps); - keys.r_disk = string(r_disk); if(typeof(seed)!=Undefined_Type) keys = struct_combine( keys, struct{ random_seed = string(seed) } ); if(qualifier_exists("dt")) -- GitLab