From d311b01119f9a1af1a1ecaf439e24cb62b705568 Mon Sep 17 00:00:00 2001 From: Matthias Kuehnel Date: Thu, 3 May 2018 21:05:08 +0200 Subject: [PATCH 1/3] SPLIT pulsarorbit INTO dopplerorbit AND pulsartorque Split the fit-function 'pulsarorbit' into the separate fit-functions 'dopplerorbit' and 'pulsartorque'. Move 'pulsarorbit_fluxerror_mc' into separate file. Require more testing and backward compatibility! --- src/fitting/fit-functions/pulsar.sl | 420 ------------------ .../fit-functions/pulsar/dopplerorbit.sl | 104 +++++ .../fit-functions/pulsar/pulsarorbit.sl | 73 +++ .../pulsar/pulsarorbit_fluxerror_mc.sl | 144 ++++++ .../fit-functions/pulsar/pulsartorque.sl | 216 +++++++++ 5 files changed, 537 insertions(+), 420 deletions(-) delete mode 100644 src/fitting/fit-functions/pulsar.sl create mode 100644 src/fitting/fit-functions/pulsar/dopplerorbit.sl create mode 100644 src/fitting/fit-functions/pulsar/pulsarorbit.sl create mode 100644 src/fitting/fit-functions/pulsar/pulsarorbit_fluxerror_mc.sl create mode 100644 src/fitting/fit-functions/pulsar/pulsartorque.sl diff --git a/src/fitting/fit-functions/pulsar.sl b/src/fitting/fit-functions/pulsar.sl deleted file mode 100644 index a7650828..00000000 --- a/src/fitting/fit-functions/pulsar.sl +++ /dev/null @@ -1,420 +0,0 @@ - -%%%%%%%%%%%%%%%%%%%%% -private define pulsarorbit_fit(lo, hi, par) -%%%%%%%%%%%%%%%%%%%%% -%!%+ -%\function{pulsarorbit} -%\synopsis{fit-function modelling a neutron star pulse period evolution} -%\usage{fit_fun("pulsarorbit");} -%\description -% The spin-up of an accreting neutron star is connected -% with its luminosity via -% \\dot{P} = - b * P^2 L^alpha -% as found by Ghosh & Lamb (1979). This fit-function -% solves this differential equation to calculate the -% spin-evolution of the neutron star. Finally, the -% Doppler shift of orbital motion is applied to the -% spin-periods to get the observed pulse-periods. -% -% The parameters of the fit-function are: -% p - spin-period at t0 (s) -% t0 - reference time t0 (MJD), has to be fixed -% a - constant spin-up or -down (s/s) -% b - torque strength (s/s @ L0 and p) -% alpha - exponent of b (try to freeze it) -% L0 - reference luminosity, has to be fixed -% asini - semi major axis (lt-s) -% porb - orbital period (days) -% tau - time of periastron passage (MJD) -% ecc - eccentricity -% omega - longitude of periastron (degrees) -% -% In order to express the torque strength, b, in the -% usual unit of a spin-change (s/s), the equation is -% modfied by two normalization constants: -% \\dot{P} = - b * (P/P0)^2 (L/L0)^alpha -% where P0 is the spin-period p at t0 and L0 a -% reference luminosity, for instance at t0 as well. -% -% The exponent, alpha, of the torque strength, b, -% depends on the accretion mechanism. After Ghosh & -% Lamb, alpha=6/7 for disk-, and alpha=1 for wind- -% accretion. -% -% In order to calculate the spin-evolution, the flux -% over time has to be assigned to the dataset using -% 'set_dataset_metadata' (see example below). The period -% measurements have to be within the time range of the -% flux evolution. The structure assigned to the dataset -% must contain the following fields: -% time - time of flux measurements (MJD, ordered) -% flux -% eps - numerical precision for solving the differential -% equation: any period at any time is more precise -% than eps (optional, default = DOUBLE_EPSILON) -% maxiter - maximum number of iterations solving the -% equation (optional, default = 10) -% -% Usually, the differential equation is found after a -% few (< 10) iterations with the default configuration. -% If the orbital motion is not needed, one might use the -% reduced fit-function 'pulsar'. -% -% To include the uncertainties of the flux measurements, -% one might use 'pulsarorbit_fluxerror_mc' to estimate -% the induced uncertainties. -%\example -% % define measurements given as time vs. period -% % as ISIS dataset -% id = define_counts(time, make_hi_grid(time), period, period_err); -% -% % assign flux measurements, e.g., a lightcurve lc, -% % to that dataset -% set_dataset_metadata(id, struct { -% time = lc.time, -% flux = lc.rate, -% eps = 1e-6, % overwrite numerical precision, e.g. the -% % uncertainties of your data is >1e-6 -% maxiter = 20 % overwrite iteration limit -% }); -% -% % set fit-function -% fit_fun("pulsarorbit(1)"); -%\seealso{pulsar, pulsarorbit_fluxerror_mc, radial_velocity_binary, -% BinaryPos, set_dataset_metadata} -%!%- -{ - % process parameters - variable p0 = par[0], t0 = par[1], a = par[2], b = par[3], - alpha = par[4], Lnorm = par[5]; - variable asini = par[6], porb = par[7], tau = par[8], - ecc = par[9], omega = par[10]; - - % get and check meta-structure - variable meta = get_dataset_metadata(Isis_Active_Dataset); - variable L = NULL, cache = 0, errorprop = 0; - if (meta != NULL) { - if (struct_field_exists(meta, "time") - + struct_field_exists(meta, "flux") < 2) { - message("warning: metadata-structure does not contain required fields, setting b = 0"); - b = 0; - } else { - L = struct { time = @(meta.time), flux = @(meta.flux) }; - cache = struct_field_exists(meta, "cache"); - errorprop = struct_field_exists(meta, "errorprop"); - if (errorprop) { - if (struct_field_exists(meta, "dflux") == 0) { - message("warning: metadata-structure does not have a 'dflux'-field required to propagate the errors"); - errorprop = 0; - } else { - L = struct_combine(L, struct { dflux = @(meta.dflux) }); - } - } - variable eps = struct_field_exists(meta, "eps") ? meta.eps : DOUBLE_EPSILON; - variable maxiter = struct_field_exists(meta, "maxiter") ? meta.maxiter : 10; - variable Lint = Double_Type[length(L.time)], dLint = @Lint; - } - } - if (meta == NULL and b != 0.) { - message("warning: metadata-structure is not set, setting b = 0"); - b = 0; - } - % initialize cache - if (cache && meta.cache == NULL) { - meta.cache = struct { pars = _Inf + Double_Type[6], pit }; - } - - % transform times into binary's barycenter - variable lot = COPY(lo); - if (porb > 0 and asini > 0) { - lot -= BinaryPos(lot; porb = porb, eccentricity = ecc, omega = omega, - asini = asini, t0 = tau)/86400; - L.time -= BinaryPos(L.time; porb = porb, eccentricity = ecc, omega = omega, - asini = asini, t0 = tau)/86400; - } - - % evaluate flux dependance if necessary - if (b != 0.) { - variable pit; % will be the intrinsic spin period - % use cache - if (cache && sum(abs(meta.cache.pars - par[[0:5]])) == 0.) { - pit = meta.cache.pit; - } - % calculate intrinsic spin period - else { - ifnot (L.time[0] <= lot[0] && lot[-1] <= L.time[-1]) { - vmessage("error (%s): flux time range does not include data", _function_name); - } else { - % index of reference time t0 - variable l, n = int(round(where_min(abs(L.time - t0))[0])); - % iterative loop because P(t) ~ int L(t)^L0 * P(t)^2 dt - pit = a*(L.time - t0)*86400 + p0; - variable tmp = Double_Type[length(pit)]; - variable pi = 0; - while (pi < maxiter and max(abs(tmp-pit)) > eps) { - % loop over flux and calculate Pdot - _for l (0, length(L.time)-1, 1) { - % indices of times before the actual time - if (l == n) { Lint[l] = 0.; } - else { - % integrate over flux (pulse period normalized to reference one) - variable i = (n < l ? [n:l] : [l:n]); - Lint[l] = (n < l ? 1 : -1) * integrate_trapez(L.time[i], - (pit[i]/p0)^2 * sign(L.flux[i])*abs(L.flux[i])^alpha - ); - } - } - % divide by normalization flux - Lint /= Lnorm^alpha; - % update pulse period - tmp = COPY(pit); - pit = (a*(L.time - t0) - b*Lint)*86400 + p0; - pi++; - } - if (pi == maxiter) { - vmessage("warning: period didn't converge after %d iterations", maxiter); - } - } - % error propagation - if (errorprop) { - variable dt = diff(L.time); dt = [dt, mean(dt)]; - _for l (0, length(L.time)-1, 1) { - if (l == n) { dLint[l] = 0.; } - else { - i = (n < l ? [n:l] : [l:n]); - dLint[l] = sqrt(sum(sqr( - .5 * (pit[i]/p0)^2 * sign(L.flux[i])*abs(L.flux[i])^(alpha-1) * L.dflux[i] * dt[i] - ))); - } - } - dLint = abs(dLint) / Lnorm^(alpha-1); - % period uncertainties - meta.errorprop = [interpol(lot, L.time, b*alpha/Lnorm * dLint*86400)]; - } - % interpolate flux times to period times - pit = [interpol(lot, L.time, pit)]; - % cache - if (cache) { meta.cache.pit = pit; } - } - } - % b = 0 - else { - pit = a*(lot - t0)*86400 + p0; - if (errorprop) { meta.errorprop = 0; } - } - - % apply orbital motion - if (porb > 0 and asini > 0) { - variable vrad = radial_velocity_binary(lot; - asini = asini, P = porb, T0 = tau, e = ecc, omega = omega, degrees); - pit *= (1. + vrad); - if (errorprop) { meta.errorprop *= (1. + vrad); } - } - - % return model - return pit; -} -add_slang_function("pulsarorbit", [ - "p [s]", "t0 [MJD]", "a [s/s]", "b [s/s @ (L0,p)]", "alpha", "L0", - "asini [lt-s]", "porb [days]", "tau [MJD]", "ecc", "omega [degrees]" -]); -set_function_category("pulsarorbit", ISIS_FUN_ADDMUL); -variable __pulsarorbit_Defaults = [ - struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 300, step = 1e-6, relstep = 1e-8 }, - struct { value = 56000.0, freeze = 1, hard_min = 15000, hard_max = int(ceil(UNIXtime2MJD(_time))), min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, - struct { value = 0.0, freeze = 1, hard_min = -1e-7, hard_max = 1e-7, min = 0, max = 1e-7, step = 1e-10, relstep = 1e-12 }, - struct { value = 0.0, freeze = 1, hard_min = 0, hard_max = 1e-7, min = 0, max = 1e-7, step = 1e-10, relstep = 1e-12 }, - struct { value = 6./7., freeze = 1, hard_min = 0, hard_max = 10, min = 6./7., max = 1., step = 1e-2, relstep = 1e-3 }, - struct { value = 1., freeze = 1, hard_min = 0, hard_max = DOUBLE_MAX, min = 0, max = 0, step = 0, relstep = 0 }, - struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, - struct { value = 10.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, - struct { value = 56000.0, freeze = 0, hard_min = 15000, hard_max = int(ceil(UNIXtime2MJD(_time))), min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, - struct { value = 0.0, freeze = 0, hard_min = 0.0, hard_max = 1.0, min = 0, max = 1.0, step = 1e-2, relstep = 1e-4 }, - struct { value = 0.0, freeze = 0, hard_min = -360, hard_max = 360, min = -180, max = 180, step = 1e-1, relstep = 1e-3 } -]; -private define pulsarorbit_default_hook(i) { return __pulsarorbit_Defaults[i]; } -set_param_default_hook("pulsarorbit", &pulsarorbit_default_hook); - -% chi-square statistics taking the model uncertainties into account -define pulsarorbit_fit_constraint(stat, pars) { - variable dmdl = array_flatten(array_struct_field(array_map( - Struct_Type, &get_dataset_metadata, all_data - ), "errorprop")); - variable data = array_map( - Struct_Type, &get_data_counts, all_data - ); - variable err = array_flatten(array_struct_field(data, "err")); - data = array_flatten(array_struct_field(data, "value")); - variable mdl = array_flatten(array_struct_field(array_map( - Struct_Type, &get_model_counts, all_data - ), "value")); - return sum(sqr(data - mdl) / (sqr(err) + sqr(dmdl))); -} - -%%%%%%%%%%%%%%%%%%%%% -private define pulsar_fit(lo, hi, par) -%%%%%%%%%%%%%%%%%%%%% -%!%+ -%\function{pulsar} -%\synopsis{fit-function modelling a neutron star pulse period evolution} -%\usage{fit_fun("pulsar");} -%\description -% Same as 'pulsarorbit' except that the Doppler shift -% of orbital motion is neglected. Thus, the corresponding -% parameters do not exist here. -%\seealso{pulsarorbit} -%!%- -{ - return pulsarorbit_fit(lo, hi, [par, 0, 0, 0, 0, 0]); -} -add_slang_function("pulsar", [ - "p [s]", "t0 [MJD]", "a [s/s]", "b [s/s @ (L0,p)]", "alpha", "L0"] -); -set_function_category("pulsar", ISIS_FUN_ADDMUL); -set_param_default_hook("pulsar", &pulsarorbit_default_hook); - - -%%%%%%%%%%%%%%%%%%%%% -define pulsarorbit_fluxerror_mc() -%%%%%%%%%%%%%%%%%%%%% -%!%+ -%\function{pulsarorbit_fluxerror_mc} -%\synopsis{estimates additional uncertainties caused by the flux measurements} -%\usage{Double_Type[] pulsarorbit_fluxerror_mc(Integer_Type dataset);} -%\qualifiers{ -% runs - number of MC loops (default: 100) -% save - save the uncertainty estimation to the -% assigned FITS-filename -% collect - file-pattern used by 'glob' to read and -% merge all FITS-files from previous runs -% modify - add the estimated uncertainties to the -% dataset(s) assigned to this qualifier -% (see below for a detailed description) -% chatty - be chatty if > 0 (default: 0) -%} -%\description -% The flux measurements used to calculate the evolution -% of the spin-period in the fit-function 'pulsarorbit' -% usually have uncertainties. These uncertainties are -% not taken into account during an ordinary fit. This -% function performs Monte Carlo simulations to estimate -% the uncertainties of the modelled period induced by -% the flux uncertainties. -% -% During each run the flux evolution associated to the -% dataset 'id' is varied using -% flux = flux + grand * flux_err -% such that synthetic flux evolutions are created. Then -% a fit is performed using only the given dataset 'id'. -% This results in many modelled period evolutions. Their -% standard deviation at each time is finally considered -% as an additional uncertainty in the period space. -% These uncertainties are returned by the function. -% -% If the 'modify'-qualifier is set, the dataset(s) -% assigned to that qualifier or, in case of NULL, the -% dataset 'id' is modified as follows: -% new_error = sqrt(sqr(data_error) + sqr(mc_error)) -% where data_error is the current 'error'-field of the -% dataset and mc_error is the estimated additional -% uncertainty as calculated by this function. In case -% of multiple given datasets, the time range of 'id' -% should include all of these datasets to get a proper -% period evolution. -%\seealso{pulsarorbit} -%!%- -{ - variable id; - if (_NARGS != 1) { help(_function_name); return; } - id = (); - - variable chatty = qualifier("chatty", 0); - - % either perform the simulation or load from previous runs - variable sim, n; - if (qualifier_exists("collect")) { - if (chatty) { message("collecting previous runs"); } - sim = array_map(Array_Type, &get_struct_field, - array_map(Struct_Type, &fits_read_table, glob(qualifier("collect"))), "periods" - ); - variable dim = [0,array_shape(sim[0])[1]]; - _for n (0, length(sim)-1, 1) { - dim[0] += array_shape(sim[n])[0]; - }; - sim = _reshape(array_flatten(sim), dim); - } else { - % get flux evolution - variable L = get_dataset_metadata(id); - if (L == NULL) { - vmessage("error (%s): metadata-structure is not set", _function_name); return; - } - if (not struct_field_exists(L, "time") or not struct_field_exists(L, "flux") or not struct_field_exists(L, "flux_err")) { - vmessage("error (%s): metadata-structure does not have the required fields", _function_name); return; - } - - variable pars = get_params; - variable incl = all_data[where(merge_struct_arrays(get_data_info(all_data)).exclude == 0)]; - exclude(all_data); - include(id); - - % Monte Carlo uncertainty estimation - seed_random(_time); - if (chatty) { message("starting MC runs"); } - variable runs = qualifier("runs", 100); - sim = Array_Type[runs]; - _for n (1, runs, 1) { - variable mcflux = L.flux + grand(length(L.flux))*L.flux_err; - mcflux[where(mcflux < 0)] = 0.; - - set_dataset_metadata(id, struct_combine(L, struct { flux = mcflux } )); - - ()=fit_counts(; fit_verbose = -1); - - sim[n-1] = get_model_counts(id).value; - } - - % restore original flux data and parameter - set_dataset_metadata(id, L); - set_params(pars); - exclude(id); - include(incl); - - % reshape - sim = _reshape(array_flatten(sim), [length(sim), length(sim[0])]); - } - - if (chatty) { - vmessage("time length = %d, runs = %d", length(sim[0,*]), length(sim[*,0])); - } - % eventually save - if (qualifier_exists("save")) { - if (chatty) { message("saving runs"); } - fits_write_binary_table(qualifier("save"), "mcflux", struct { periods = sim }); - } - - % calculate standard deviation - variable err = Double_Type[length(sim[0,*])]; - _for n (0, length(err)-1, 1) { - err[n] = moment(sim[*,n]).sdev; - } - - % modify dataset - if (qualifier_exists("modify")) { - variable mid = qualifier("modify"); - if (mid == NULL) { mid = [id]; } - if (typeof(mid) != Array_Type or (_typeof(mid) != Integer_Type and _typeof(mid) != UInteger_Type)) { - vmessage("error (%s): 'modify' must be an array of integers", _function_name); return; - } - - % loop and modify - variable t = get_data_counts(id).bin_lo; - foreach n (mid) { - if (chatty) { vmessage("modifying dataset %d", n); } - variable cts = get_data_counts(n); - cts.err = sqrt(sqr(cts.err) + sqr(interpol(cts.bin_lo, t, err))); - put_data(n, cts); - } - } - - return err; -} diff --git a/src/fitting/fit-functions/pulsar/dopplerorbit.sl b/src/fitting/fit-functions/pulsar/dopplerorbit.sl new file mode 100644 index 00000000..9109c4df --- /dev/null +++ b/src/fitting/fit-functions/pulsar/dopplerorbit.sl @@ -0,0 +1,104 @@ +%%%%%%%%%%%%%%%%%%%%% +private define dopplerorbit_fit(lo, hi, par) +%%%%%%%%%%%%%%%%%%%%% +%!%+ +%\function{dopplerorbit} +%\synopsis{fit-function for the orbital Doppler shift factor} +%\usage{fit_fun("dopplerorbit(1; qualifiers) * ...");} +%\qualifiers{ +% t90: switch to the time of mean longitude 90 degrees (see below) +% K : switch to velocity semi amplitude (see below) +%} +%\description +% Calculates the radial velocoity, v_rad, of a star in a binary and +% returns the corresponding Doppler factor, beta, given by +% beta = (1 + V_rad/Const_c) +% assuming that V_rad/Const_c << 1. The grid of the data (lo,hi) +% the model is fitted to needs has to be the time in MJD. +% +% The parameters of the fit-function are: +% asini - semi major axis (lt-s) +% or velocity semi amplitude (km/s) +% porb - orbital period (days) +% tau - time of periastron passage (MJD) +% or time of mean longitude 90 degrees (MJD) +% ecc - eccentricity +% omega - longitude of periastron (degrees) +% v0 - systemic velocity (km/s) +% +% In case the 't90'-qualifier is set, the parameter 'tau' has the +% meaning of the time of mean longitude of 90 degrees. In case the +% 'K'-qualifiers is set, the parameter 'asini' has the meaning of +% the velocity semi amplitude. +%\seealso{radial_velocity_binary, BinaryPos, KeplerEquation} +%!%- +{ + % input parameters + variable asini = par[0], porb = par[1], tau = par[2], + ecc = par[3], omega = par[4], v0 = par[5]; + + % sanity checks + if (porb <= 0 || asini <= 0) { + throw InvalidParmError, "porb and asini have to be greater zero"$; + } + + % build qualifier struct for 'radial_velocity_binary' + % and 'BinaryPos' + variable rvb = struct { + P = porb, omega = omega, e = ecc, degrees, v0 = v0 + }; + variable bp = struct { + porb = porb, eccentricity = ecc, omega = omega + }; + if (qualifier_exists("t90")) { + rvb = struct_combine(rvb, struct { T90 = tau }); + bp = struct_combine(bp, struct { t90 = tau }); + } else { + rvb = struct_combine(rvb, struct { T0 = tau }); + bp = struct_combine(bp, struct { t0 = tau }); + } + if (qualifier_exists("K")) { + rvb = struct_combine(rvb, struct { K = asini }); + bp = struct_combine(bp, struct { % need to convert K to asini + asini = asini / (2*PI) * (porb*86400. * sqrt((1.-ecc)*(1.+ecc))) / (Const_c*1e-5) + }); + } else { + rvb = struct_combine(rvb, struct { asini = asini*Const_c*1e-5 }); + bp = struct_combine(bp, struct { asini = asini }); + } + + % transform times into binary's barycenter + variable cache = fitfun_get_cache(); + if (cache.binarytime == NULL || any(cache.lastpar[[:4]] != par[[:4]])) { + cache.binarytime = lo - BinaryPos(lo;; bp)/86400.; + cache.lastpar = par; + } + + % compute radial velocity and corresponding Doppler factor + return (1. + radial_velocity_binary(cache.binarytime;; rvb) / (Const_c*1e-5)); +} + + +% add fit-function +add_slang_function("dopplerorbit", [ + "asini [lt-s]", "porb [days]", "tau [MJD]", "ecc", "omega [degrees]", "v0 [km/s]" +]); +set_function_category("dopplerorbit", ISIS_FUN_ADDMUL); + + +% define parameter defaults +private define dopplerorbit_defaults(i) { + return [ + struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, + struct { value = 10.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, + struct { value = 56000.0, freeze = 0, hard_min = 15000, hard_max = 100000, min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, + struct { value = 0.0, freeze = 0, hard_min = 0.0, hard_max = 1.0, min = 0, max = 1.0, step = 1e-2, relstep = 1e-4 }, + struct { value = 0.0, freeze = 0, hard_min = -360, hard_max = 360, min = -180, max = 180, step = 1e-1, relstep = 1e-3 }, + struct { value = 0.0, freeze = 1, hard_min = -1000.0, hard_max = 1000.0, min = -200, max = 200, step = 1e-2, relstep = 1e-4 } + ][i]; +} +set_param_default_hook("dopplerorbit", &dopplerorbit_defaults); + + +% setup the fit-function's cache +fitfun_init_cache("dopplerorbit", &dopplerorbit_fit, struct { binarytime = NULL, lastpar }); diff --git a/src/fitting/fit-functions/pulsar/pulsarorbit.sl b/src/fitting/fit-functions/pulsar/pulsarorbit.sl new file mode 100644 index 00000000..9398fc14 --- /dev/null +++ b/src/fitting/fit-functions/pulsar/pulsarorbit.sl @@ -0,0 +1,73 @@ +%%%%%%%%%%%%%%%%%%%%% +private define pulsarorbit_fit(lo, hi, par) +%%%%%%%%%%%%%%%%%%%%% +%!%+ +%\function{pulsarorbit} +%\synopsis{fit-function modelling a neutron star pulse period evolution} +%\usage{fit_fun("pulsarorbit");} +%\description +% Deprecated, use +% fit_fun("dopplerorbit*pulsartorque"); +% from now on. +%!%- +{ + message("*** DEPRECATED ***"); + return eval_fun2("dopplerorbit", lo, hi, [par[[6:]], 0]) + * eval_fun2("pulsartorque", lo, hi, par[[:5]]); +} + + +% add fit-function +add_slang_function("pulsarorbit", [ + "p [s]", "t0 [MJD]", "a [s/s]", "b [s/s @ (L0,p)]", "alpha", "L0", + "asini [lt-s]", "porb [days]", "tau [MJD]", "ecc", "omega [degrees]" +]); +set_function_category("pulsarorbit", ISIS_FUN_ADDMUL); + + +% define parameter defaults +private define pulsarorbit_defaults(i) { + return [ + struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 300, step = 1e-6, relstep = 1e-8 }, + struct { value = 56000.0, freeze = 1, hard_min = 15000, hard_max = int(ceil(UNIXtime2MJD(_time))), min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, + struct { value = 0.0, freeze = 1, hard_min = -1e-7, hard_max = 1e-7, min = 0, max = 1e-7, step = 1e-10, relstep = 1e-12 }, + struct { value = 0.0, freeze = 1, hard_min = 0, hard_max = 1e-7, min = 0, max = 1e-7, step = 1e-10, relstep = 1e-12 }, + struct { value = 6./7., freeze = 1, hard_min = 0, hard_max = 10, min = 6./7., max = 1., step = 1e-2, relstep = 1e-3 }, + struct { value = 1., freeze = 1, hard_min = 0, hard_max = DOUBLE_MAX, min = 0, max = 0, step = 0, relstep = 0 }, + struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, + struct { value = 10.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 1000, step = 1e-2, relstep = 1e-4 }, + struct { value = 56000.0, freeze = 0, hard_min = 15000, hard_max = int(ceil(UNIXtime2MJD(_time))), min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, + struct { value = 0.0, freeze = 0, hard_min = 0.0, hard_max = 1.0, min = 0, max = 1.0, step = 1e-2, relstep = 1e-4 }, + struct { value = 0.0, freeze = 0, hard_min = -360, hard_max = 360, min = -180, max = 180, step = 1e-1, relstep = 1e-3 } + ][i]; +} +set_param_default_hook("pulsarorbit", &pulsarorbit_defaults); + + +%%%%%%%%%%%%%%%%%%%%% +private define pulsar_fit(lo, hi, par) +%%%%%%%%%%%%%%%%%%%%% +%!%+ +%\function{pulsar} +%\synopsis{fit-function modelling a neutron star pulse period evolution} +%\usage{fit_fun("pulsar");} +%\description +% Deprecated, use +% fit_fun("pulsartorque"); +% from now on. +%!%- +{ + message("*** DEPRECATED ***"); + return eval_fun2("pulsartorque", lo, hi, par); +} + + +% add fit-function +add_slang_function("pulsar", [ + "p [s]", "t0 [MJD]", "a [s/s]", "b [s/s @ (L0,p)]", "alpha", "L0"] +); +set_function_category("pulsar", ISIS_FUN_ADDMUL); + + +% define parameter defaults +set_param_default_hook("pulsar", &pulsarorbit_defaults); diff --git a/src/fitting/fit-functions/pulsar/pulsarorbit_fluxerror_mc.sl b/src/fitting/fit-functions/pulsar/pulsarorbit_fluxerror_mc.sl new file mode 100644 index 00000000..53a69449 --- /dev/null +++ b/src/fitting/fit-functions/pulsar/pulsarorbit_fluxerror_mc.sl @@ -0,0 +1,144 @@ +%%%%%%%%%%%%%%%%%%%%% +define pulsarorbit_fluxerror_mc() +%%%%%%%%%%%%%%%%%%%%% +%!%+ +%\function{pulsarorbit_fluxerror_mc} +%\synopsis{estimates additional uncertainties caused by the flux measurements} +%\usage{Double_Type[] pulsarorbit_fluxerror_mc(Integer_Type dataset);} +%\qualifiers{ +% runs - number of MC loops (default: 100) +% save - save the uncertainty estimation to the +% assigned FITS-filename +% collect - file-pattern used by 'glob' to read and +% merge all FITS-files from previous runs +% modify - add the estimated uncertainties to the +% dataset(s) assigned to this qualifier +% (see below for a detailed description) +% chatty - be chatty if > 0 (default: 0) +%} +%\description +% The flux measurements used to calculate the evolution +% of the spin-period in the fit-function 'pulsarorbit' +% usually have uncertainties. These uncertainties are +% not taken into account during an ordinary fit. This +% function performs Monte Carlo simulations to estimate +% the uncertainties of the modelled period induced by +% the flux uncertainties. +% +% During each run the flux evolution associated to the +% dataset 'id' is varied using +% flux = flux + grand * flux_err +% such that synthetic flux evolutions are created. Then +% a fit is performed using only the given dataset 'id'. +% This results in many modelled period evolutions. Their +% standard deviation at each time is finally considered +% as an additional uncertainty in the period space. +% These uncertainties are returned by the function. +% +% If the 'modify'-qualifier is set, the dataset(s) +% assigned to that qualifier or, in case of NULL, the +% dataset 'id' is modified as follows: +% new_error = sqrt(sqr(data_error) + sqr(mc_error)) +% where data_error is the current 'error'-field of the +% dataset and mc_error is the estimated additional +% uncertainty as calculated by this function. In case +% of multiple given datasets, the time range of 'id' +% should include all of these datasets to get a proper +% period evolution. +%\seealso{pulsarorbit} +%!%- +{ + variable id; + if (_NARGS != 1) { help(_function_name); return; } + id = (); + + variable chatty = qualifier("chatty", 0); + + % either perform the simulation or load from previous runs + variable sim, n; + if (qualifier_exists("collect")) { + if (chatty) { message("collecting previous runs"); } + sim = array_map(Array_Type, &get_struct_field, + array_map(Struct_Type, &fits_read_table, glob(qualifier("collect"))), "periods" + ); + variable dim = [0,array_shape(sim[0])[1]]; + _for n (0, length(sim)-1, 1) { + dim[0] += array_shape(sim[n])[0]; + }; + sim = _reshape(array_flatten(sim), dim); + } else { + % get flux evolution + variable L = get_dataset_metadata(id); + if (L == NULL) { + vmessage("error (%s): metadata-structure is not set", _function_name); return; + } + if (not struct_field_exists(L, "time") or not struct_field_exists(L, "flux") or not struct_field_exists(L, "flux_err")) { + vmessage("error (%s): metadata-structure does not have the required fields", _function_name); return; + } + + variable pars = get_params; + variable incl = all_data[where(merge_struct_arrays(get_data_info(all_data)).exclude == 0)]; + exclude(all_data); + include(id); + + % Monte Carlo uncertainty estimation + seed_random(_time); + if (chatty) { message("starting MC runs"); } + variable runs = qualifier("runs", 100); + sim = Array_Type[runs]; + _for n (1, runs, 1) { + variable mcflux = L.flux + grand(length(L.flux))*L.flux_err; + mcflux[where(mcflux < 0)] = 0.; + + set_dataset_metadata(id, struct_combine(L, struct { flux = mcflux } )); + + ()=fit_counts(; fit_verbose = -1); + + sim[n-1] = get_model_counts(id).value; + } + + % restore original flux data and parameter + set_dataset_metadata(id, L); + set_params(pars); + exclude(id); + include(incl); + + % reshape + sim = _reshape(array_flatten(sim), [length(sim), length(sim[0])]); + } + + if (chatty) { + vmessage("time length = %d, runs = %d", length(sim[0,*]), length(sim[*,0])); + } + % eventually save + if (qualifier_exists("save")) { + if (chatty) { message("saving runs"); } + fits_write_binary_table(qualifier("save"), "mcflux", struct { periods = sim }); + } + + % calculate standard deviation + variable err = Double_Type[length(sim[0,*])]; + _for n (0, length(err)-1, 1) { + err[n] = moment(sim[*,n]).sdev; + } + + % modify dataset + if (qualifier_exists("modify")) { + variable mid = qualifier("modify"); + if (mid == NULL) { mid = [id]; } + if (typeof(mid) != Array_Type or (_typeof(mid) != Integer_Type and _typeof(mid) != UInteger_Type)) { + vmessage("error (%s): 'modify' must be an array of integers", _function_name); return; + } + + % loop and modify + variable t = get_data_counts(id).bin_lo; + foreach n (mid) { + if (chatty) { vmessage("modifying dataset %d", n); } + variable cts = get_data_counts(n); + cts.err = sqrt(sqr(cts.err) + sqr(interpol(cts.bin_lo, t, err))); + put_data(n, cts); + } + } + + return err; +} diff --git a/src/fitting/fit-functions/pulsar/pulsartorque.sl b/src/fitting/fit-functions/pulsar/pulsartorque.sl new file mode 100644 index 00000000..69e0e91f --- /dev/null +++ b/src/fitting/fit-functions/pulsar/pulsartorque.sl @@ -0,0 +1,216 @@ +% define the simplified ordinary differential equation +% \dot{P} = a - b * P^2 L^alpha +% after Ghosh & Lamb (1979) +private define _ode_torque(p, t) { + variable p0 = qualifier("p0"); + variable a = qualifier("a"); + variable b = qualifier("b"); + variable alpha = qualifier("alpha"); + variable lc = qualifier("lc"); + + variable rate = interpol(t, lc.time, lc.rate); + + return a - b * (p/p0)^2 * sign(rate)*abs(rate)^alpha; +} + +%%%%%%%%%%%%%%%%%%%%% +private define pulsartorque_fit(lo, hi, par) +%%%%%%%%%%%%%%%%%%%%% +%!%+ +%\function{pulsartorque} +%\synopsis{fit-function modelling the accretion torque of a neutron star (simple way)} +%\usage{fit_fun("pulsartorque"); +% or fit_fun("dopplerorbit*pulsartorque);} +%\qualifiers{ +% eps - numerical precision for solving the differential +% equation: any period at any time is more precise than +% eps (default: DOUBLE_EPSILON) +% maxiter - maximum number of iterations solving the equation +% (default: 10) +%} +%\description +% The spin-up of an accreting neutron star is connected with its +% luminosity via +% \\dot{P} = - b * P^2 L^alpha +% as found by Ghosh & Lamb (1979). This fit-function solves this +% differential equation to calculate the spin-evolution of the +% neutron star. +% +% The parameters of the fit-function are: +% p - spin-period at t0 (s) +% t0 - reference time t0 (MJD), has to be fixed +% a - constant spin-up or -down (s/s) +% b - torque strength (s/s @ L0 and p) +% alpha - exponent of b (try to freeze it) +% L0 - reference luminosity, has to be fixed +% +% In order to express the torque strength, b, in the usual unit of +% a spin-change (s/s), the equation is modfied by two normalization +% constants: +% \\dot{P} = - b * (P/P0)^2 (L/L0)^alpha +% where P0 is the spin-period p at t0 and L0 a reference +% luminosity, for instance at t0 as well. +% +% The exponent, alpha, of the torque strength, b, depends on the +% accretion mechanism. After Ghosh & Lamb, alpha=6/7 for disk-, and +% alpha=1 for wind- accretion. +% +% In order to calculate the spin-evolution, the count rate over +% time, i.e., a light curve has to be assigned to the dataset using +% 'set_dataset_metadata' (see example below). The period +% measurements have to be within the time range than of the light +% curve. The structure assigned to the dataset must contain the +% following fields: +% time - time (MJD, ordered) +% rate - count rate at 'time' +% +% Usually, the differential equation is found after a +% few (< 10) iterations with the default configuration. +% +% NOTE: the full Ghosh & Lamb accretion torque theory provides +% equations for the torque strength, b, as a function of the +% neutron star parameters. This is not implemented here for +% simplicity. The full theory is implemented by the 'pulsarGL79' +% fit-function. +% +% NOTE: in case you add the Doppler shift of orbital motion to the +% model using the 'dopplerorbit' fit-function, make sure that this +% fit-function is calculated *first* in order to transform the time +% ('lo' parameter of the fit-function) into the barycenter of the +% binary! This corrected time is shared among these fit-functions +% via the ISISscripts caching extension. +%\example +% % define the measured period evolution as ISIS dataset +% id = define_counts(time, make_hi_grid(time), period, period_err); +% +% % assign the light curve to this dataset +% set_dataset_metadata(id, struct { +% time = lc.time, +% rate = lc.rate +% }); +% +% % examples for setting the fit-function +% fit_fun("pulsarorbit(1)"); +% fit_fun("pulsarorbit(1; eps = 1e-4)"); % decrease the precision +% fit_fun("dopplerorbit(1)*pulsarorbit(1)"); % add orbital Doppler shift +%\seealso{pulsarGL79, dopplerorbit, set_dataset_metadata, fitfun_cache} +%!%- +{ + % input parameters + variable p0 = par[0], t0 = par[1], a = par[2]*86400, b = par[3]*86400, + alpha = par[4], Lnorm = par[5]; + + % get and check meta-structure + variable meta = get_dataset_metadata(Isis_Active_Dataset); + if (typeof(meta) != Struct_Type) { + throw InvalidParmError, "metadata-structure is not defined!"$; + } + if (sum(array_map(Integer_Type, &struct_field_exists, meta,["time", "rate", "flux"])) < 2) { + throw InvalidParmError, "metadata-structure does not contain required fields!"$; + } + + % define variables + variable lcrate; + variable lctime = @(meta.time); + if (struct_field_exists(meta, "flux")) { % backward compatibility + lcrate = @(meta.flux) / Lnorm; + } else { lcrate = @(meta.rate) / Lnorm; } + if (lo[0] < meta.time[0] || lo[-1] > meta.time[-1]) { + throw InvalidParmError, "light curve does not include pulse period data"$; + } + + % define cache + variable cache = fitfun_get_cache(); + + % get binary corrected times + variable binarytime = fitfun_get_cache("dopplerorbit"); + if (binarytime != NULL && binarytime.binarytime != NULL) { + binarytime = binarytime.binarytime; +% lc.time = interpol_points(lc.time, lo, binarytime); +% print(length(lc.time)); + } + else { binarytime = lo; } + + % calculate intrinsic spin period evolution by solving the + % differential equation + if (b != 0. && (cache.spin == NULL || any(cache.lastpar != par))) { + cache.spin = solveODEbyIntegrate( + &_ode_torque, binarytime; t0 = t0, x0 = p0, qualifiers = struct { + p0 = p0, a = a, b = b, alpha = alpha, lc = struct { time = lctime, rate = lcrate } + } + ); + } + % b = 0 + else if (cache.spin == NULL || any(cache.lastpar[[:2]] != par[[:2]])) { + cache.spin = a*(binarytime - t0)*86400 + p0; +% if (errorprop) { meta.errorprop = 0; } + } + + + % % error propagation + % if (errorprop) { + % variable dt = diff(L.time); dt = [dt, mean(dt)]; + % _for l (0, length(L.time)-1, 1) { + % if (l == n) { dLint[l] = 0.; } + % else { + % i = (n < l ? [n:l] : [l:n]); + % dLint[l] = sqrt(sum(sqr( + % .5 * (pit[i]/p0)^2 * sign(L.flux[i])*abs(L.flux[i])^(alpha-1) * L.dflux[i] * dt[i] + % ))); + % } + % } + % dLint = abs(dLint) / Lnorm^(alpha-1); + % % period uncertainties + % meta.errorprop = [interpol(lot, L.time, b*alpha/Lnorm * dLint*86400)]; + % } + % cache + % if (cache) { meta.cache.pit = pit; } + + % update cached parameters + cache.lastpar = par; + + % return model + return cache.spin * (hi-lo); +} + + +% add fit-function +add_slang_function("pulsartorque", [ + "p [s]", "t0 [MJD]", "a [s/s]", "b [s/s @ (L0,p)]", "alpha", "L0", +]); +set_function_category("pulsarorbit", ISIS_FUN_ADDMUL); + + +% define parameter defaults +private define pulsartorque_defaults(i) { + return [ + struct { value = 100.0, freeze = 0, hard_min = 0.0, hard_max = 10000, min = 0, max = 300, step = 1e-6, relstep = 1e-8 }, + struct { value = 56000.0, freeze = 1, hard_min = 15000, hard_max = 100000, min = 15000, max = int(ceil(UNIXtime2MJD(_time))), step = 1e-1, relstep = 1e-3 }, + struct { value = 0.0, freeze = 1, hard_min = -1e-7, hard_max = 1e-7, min = 0, max = 1e-7, step = 1e-10, relstep = 1e-12 }, + struct { value = 0.0, freeze = 1, hard_min = 0, hard_max = 1e-5, min = 0, max = 1e-5, step = 1e-10, relstep = 1e-12 }, + struct { value = 6./7., freeze = 1, hard_min = 0, hard_max = 10, min = 6./7., max = 1., step = 1e-2, relstep = 1e-3 }, + struct { value = 1., freeze = 1, hard_min = 0, hard_max = DOUBLE_MAX, min = 0, max = 0, step = 0, relstep = 0 } + ][i]; +} +set_param_default_hook("pulsartorque", &pulsartorque_defaults); + + +% setup the fit-function's cache +fitfun_init_cache("pulsartorque", &pulsartorque_fit, struct { spin = NULL, lctime, lastpar }); + + +% chi-square statistics taking the model uncertainties into account +define pulsarorbit_fit_constraint(stat, pars) { + variable dmdl = array_flatten(array_struct_field(array_map( + Struct_Type, &get_dataset_metadata, all_data + ), "errorprop")); + variable data = array_map( + Struct_Type, &get_data_counts, all_data + ); + variable err = array_flatten(array_struct_field(data, "err")); + data = array_flatten(array_struct_field(data, "value")); + variable mdl = array_flatten(array_struct_field(array_map( + Struct_Type, &get_model_counts, all_data + ), "value")); + return sum(sqr(data - mdl) / (sqr(err) + sqr(dmdl))); +} -- GitLab From 9a96cc62fe84ea6e2b2b26eadd47e00cbc4e53bc Mon Sep 17 00:00:00 2001 From: Matthias Kuehnel Date: Thu, 3 May 2018 23:42:44 +0200 Subject: [PATCH 2/3] ADD BINARYCOR-CORRECTION OF METADATA TO 'binaryorbit' The fit-function 'binaryorbit' is able to correct specific "time"-fields fo the datasets metadata for binary motion --- .../fit-functions/pulsar/dopplerorbit.sl | 84 +++++++++++++++++-- 1 file changed, 75 insertions(+), 9 deletions(-) diff --git a/src/fitting/fit-functions/pulsar/dopplerorbit.sl b/src/fitting/fit-functions/pulsar/dopplerorbit.sl index 9109c4df..6f4ff435 100644 --- a/src/fitting/fit-functions/pulsar/dopplerorbit.sl +++ b/src/fitting/fit-functions/pulsar/dopplerorbit.sl @@ -8,6 +8,10 @@ private define dopplerorbit_fit(lo, hi, par) %\qualifiers{ % t90: switch to the time of mean longitude 90 degrees (see below) % K : switch to velocity semi amplitude (see below) +% metacor: string array of structure-field-paths inside the meta- +% data of the current dataset, which will be assumed to +% be time in MJD and corrected for binary motion (default: +% ["time"]). Set to NULL to disable any correction. %} %\description % Calculates the radial velocoity, v_rad, of a star in a binary and @@ -30,7 +34,20 @@ private define dopplerorbit_fit(lo, hi, par) % meaning of the time of mean longitude of 90 degrees. In case the % 'K'-qualifiers is set, the parameter 'asini' has the meaning of % the velocity semi amplitude. -%\seealso{radial_velocity_binary, BinaryPos, KeplerEquation} +% +% The influence of the binary motion is removed from the input +% time-array 'lo' using 'BinaryCor'. The corrected times are saved +% as 'binarytime' into the fitfun_cache of the fit-function. This +% can be used by other fit-functions, like 'pulsartaylor'. +% +% The fit-function also corrects time arrays inside the metadata of +% the current dataset (Isis_Active_Dataset). The name of the fields +% are set by the 'metacor'-qualifier, which is an array of strings +% giving the "path" to the field. By default, the 'time'-field +% inside the metadata is corrected in order to simplify the usage +% of, e.g., 'pulsartorque'. The corrected times are save into the +% 'metacor'-struct inside the cache using the same "path". +%\seealso{radial_velocity_binary, BinaryCor, KeplerEquation} %!%- { % input parameters @@ -43,35 +60,82 @@ private define dopplerorbit_fit(lo, hi, par) } % build qualifier struct for 'radial_velocity_binary' - % and 'BinaryPos' + % and 'BinaryCor' variable rvb = struct { P = porb, omega = omega, e = ecc, degrees, v0 = v0 }; - variable bp = struct { + variable bc = struct { porb = porb, eccentricity = ecc, omega = omega }; if (qualifier_exists("t90")) { rvb = struct_combine(rvb, struct { T90 = tau }); - bp = struct_combine(bp, struct { t90 = tau }); + bc = struct_combine(bc, struct { t90 = tau }); } else { rvb = struct_combine(rvb, struct { T0 = tau }); - bp = struct_combine(bp, struct { t0 = tau }); + bc = struct_combine(bc, struct { t0 = tau }); } if (qualifier_exists("K")) { rvb = struct_combine(rvb, struct { K = asini }); - bp = struct_combine(bp, struct { % need to convert K to asini + bc = struct_combine(bc, struct { % need to convert K to asini asini = asini / (2*PI) * (porb*86400. * sqrt((1.-ecc)*(1.+ecc))) / (Const_c*1e-5) }); } else { rvb = struct_combine(rvb, struct { asini = asini*Const_c*1e-5 }); - bp = struct_combine(bp, struct { asini = asini }); + bc = struct_combine(bc, struct { asini = asini }); } % transform times into binary's barycenter variable cache = fitfun_get_cache(); + variable dometacor = 0; if (cache.binarytime == NULL || any(cache.lastpar[[:4]] != par[[:4]])) { - cache.binarytime = lo - BinaryPos(lo;; bp)/86400.; + cache.binarytime = BinaryCor(lo;; bc); cache.lastpar = par; + dometacor = 1; + } + + % purge binary corrected metadata + variable metafields = qualifier("metacor", ["time"]); + variable metacor = length(get_struct_field_names(cache.metacor)); + if (metafields == NULL) { + set_struct_field(cache, "metacor", @Struct_Type(String_Type[0])); + } + else if (length(metafields) != metacor) { dometacor = 1; } + % binary correct metadata + if (dometacor) { + % transform metadata fields + if (metafields != NULL) { + variable f, tfield, cfield, cname; + % loop over fields to correct + _for f (0, length(metafields)-1, 1) { + tfield = get_dataset_metadata(Isis_Active_Dataset); + cname = "metacor"; cfield = cache; + if (tfield == NULL) { continue; } + % walk along the path and recursively set 'tfield' and 'cfield' + variable path; + foreach path (strchop(metafields[f], '.', 0)) { + ifnot (struct_field_exists(tfield, path)) { tfield = NULL; break; } + % create path inside the cache if necessary + ifnot (struct_field_exists(get_struct_field(cfield, cname), path)) { + set_struct_field(cfield, cname, struct_combine( + get_struct_field(cfield, cname), path + )); + } + % set 'tfield' and 'cfield' + tfield = get_struct_field(tfield, path); + cfield = get_struct_field(cfield, cname); + cname = path; + % create empty structure in 'cfield' + if (typeof(tfield) == Struct_Type) { + if (get_struct_field(cfield, path) == NULL) { + set_struct_field(cfield, path, @Struct_Type(String_Type[0])); + } + } + } + if (tfield == NULL) { continue; } + % correct the times + set_struct_field(cfield, cname, BinaryCor(tfield;; bc)); + } + } } % compute radial velocity and corresponding Doppler factor @@ -101,4 +165,6 @@ set_param_default_hook("dopplerorbit", &dopplerorbit_defaults); % setup the fit-function's cache -fitfun_init_cache("dopplerorbit", &dopplerorbit_fit, struct { binarytime = NULL, lastpar }); +fitfun_init_cache("dopplerorbit", &dopplerorbit_fit, struct { + binarytime = NULL, lastpar, metacor = @Struct_Type(String_Type[0]) +}); -- GitLab From 3ea156fb04d87fe1f0d0f28be9450917678ddc00 Mon Sep 17 00:00:00 2001 From: Matthias Kuehnel Date: Fri, 4 May 2018 12:50:54 +0200 Subject: [PATCH 3/3] FIX the binning and binary-correction in 'pulsartorque'-fit-function --- .../fit-functions/pulsar/pulsartorque.sl | 107 ++++++++++-------- 1 file changed, 58 insertions(+), 49 deletions(-) diff --git a/src/fitting/fit-functions/pulsar/pulsartorque.sl b/src/fitting/fit-functions/pulsar/pulsartorque.sl index 69e0e91f..f400258b 100644 --- a/src/fitting/fit-functions/pulsar/pulsartorque.sl +++ b/src/fitting/fit-functions/pulsar/pulsartorque.sl @@ -1,6 +1,7 @@ % define the simplified ordinary differential equation % \dot{P} = a - b * P^2 L^alpha -% after Ghosh & Lamb (1979) +% after Ghosh & Lamb (1979), which will be solved by the +% fit-function below. private define _ode_torque(p, t) { variable p0 = qualifier("p0"); variable a = qualifier("a"); @@ -8,9 +9,11 @@ private define _ode_torque(p, t) { variable alpha = qualifier("alpha"); variable lc = qualifier("lc"); - variable rate = interpol(t, lc.time, lc.rate); + % not needed because the input time, t, is equal to the light-curve + % time grid + % variable rate = interpol(t, lc.time, lc.rate); - return a - b * (p/p0)^2 * sign(rate)*abs(rate)^alpha; + return a - b * (p/p0)^2 * sign(lc.rate)*abs(lc.rate)^alpha; } %%%%%%%%%%%%%%%%%%%%% @@ -21,13 +24,6 @@ private define pulsartorque_fit(lo, hi, par) %\synopsis{fit-function modelling the accretion torque of a neutron star (simple way)} %\usage{fit_fun("pulsartorque"); % or fit_fun("dopplerorbit*pulsartorque);} -%\qualifiers{ -% eps - numerical precision for solving the differential -% equation: any period at any time is more precise than -% eps (default: DOUBLE_EPSILON) -% maxiter - maximum number of iterations solving the equation -% (default: 10) -%} %\description % The spin-up of an accreting neutron star is connected with its % luminosity via @@ -57,28 +53,32 @@ private define pulsartorque_fit(lo, hi, par) % % In order to calculate the spin-evolution, the count rate over % time, i.e., a light curve has to be assigned to the dataset using -% 'set_dataset_metadata' (see example below). The period +% 'set_dataset_metadata' (see example below). The period % measurements have to be within the time range than of the light % curve. The structure assigned to the dataset must contain the % following fields: % time - time (MJD, ordered) % rate - count rate at 'time' % -% Usually, the differential equation is found after a -% few (< 10) iterations with the default configuration. +% NOTE: in case you add the Doppler shift of orbital motion to the +% model using the 'dopplerorbit' fit-function, make sure that +% this fit-function is calculated *first* in order to +% transform the times ('lo' parameter of the fit-function +% /and/ light curve 'time' grid) into the barycenter of the +% binary! These corrected times are shared among these +% fit-functions via the ISISscripts caching extension. Then +% the reference time, t0, is interpreted as binary corrected! % % NOTE: the full Ghosh & Lamb accretion torque theory provides -% equations for the torque strength, b, as a function of the -% neutron star parameters. This is not implemented here for -% simplicity. The full theory is implemented by the 'pulsarGL79' -% fit-function. +% equations for the torque strength, b, as a function of the +% neutron star parameters. This is not implemented here for +% simplicity. The full theory is implemented in the +% 'pulsarGL79' fit-function. % -% NOTE: in case you add the Doppler shift of orbital motion to the -% model using the 'dopplerorbit' fit-function, make sure that this -% fit-function is calculated *first* in order to transform the time -% ('lo' parameter of the fit-function) into the barycenter of the -% binary! This corrected time is shared among these fit-functions -% via the ISISscripts caching extension. +% NOTE: the model is calculated on the light curve time grid in +% order to ensure that the integration includes the reference +% time, t0. This model is interpolated onto the period time +% grid in the end.% %\example % % define the measured period evolution as ISIS dataset % id = define_counts(time, make_hi_grid(time), period, period_err); @@ -90,9 +90,8 @@ private define pulsartorque_fit(lo, hi, par) % }); % % % examples for setting the fit-function -% fit_fun("pulsarorbit(1)"); -% fit_fun("pulsarorbit(1; eps = 1e-4)"); % decrease the precision -% fit_fun("dopplerorbit(1)*pulsarorbit(1)"); % add orbital Doppler shift +% fit_fun("pulsartorque(1)"); +% fit_fun("dopplerorbit(1)*pulsartorque(1)"); % add orbital Doppler shift %\seealso{pulsarGL79, dopplerorbit, set_dataset_metadata, fitfun_cache} %!%- { @@ -123,26 +122,34 @@ private define pulsartorque_fit(lo, hi, par) variable cache = fitfun_get_cache(); % get binary corrected times - variable binarytime = fitfun_get_cache("dopplerorbit"); - if (binarytime != NULL && binarytime.binarytime != NULL) { - binarytime = binarytime.binarytime; -% lc.time = interpol_points(lc.time, lo, binarytime); -% print(length(lc.time)); + variable orbitcache = fitfun_get_cache("dopplerorbit"); + variable binarytime = lo; + if (orbitcache != NULL) { + if (orbitcache.binarytime != NULL) { + binarytime = orbitcache.binarytime; + } + if (struct_field_exists(orbitcache.metacor, "time") && orbitcache.metacor.time != NULL) { + lctime = orbitcache.metacor.time; + } } - else { binarytime = lo; } % calculate intrinsic spin period evolution by solving the - % differential equation + % differential equation on the light curve time grid. this is + % necessary because the reference time might be outside of the + % period time range -> using the light curve we still can predict + % the period at these times. a user-defined time grid using + % set_eval_grid_method does not work because ISIS rebins (=sums up!) + % the resulting model onto the data time grid. if (b != 0. && (cache.spin == NULL || any(cache.lastpar != par))) { cache.spin = solveODEbyIntegrate( - &_ode_torque, binarytime; t0 = t0, x0 = p0, qualifiers = struct { + &_ode_torque, lctime; t0 = t0, x0 = p0, qualifiers = struct { p0 = p0, a = a, b = b, alpha = alpha, lc = struct { time = lctime, rate = lcrate } } ); } % b = 0 else if (cache.spin == NULL || any(cache.lastpar[[:2]] != par[[:2]])) { - cache.spin = a*(binarytime - t0)*86400 + p0; + cache.spin = a*(lctime - t0) + p0; % if (errorprop) { meta.errorprop = 0; } } @@ -169,8 +176,8 @@ private define pulsartorque_fit(lo, hi, par) % update cached parameters cache.lastpar = par; - % return model - return cache.spin * (hi-lo); + % return model on the data time grid (as an array always!) + return [ interpol(binarytime, lctime, cache.spin) ]; } @@ -201,16 +208,18 @@ fitfun_init_cache("pulsartorque", &pulsartorque_fit, struct { spin = NULL, lctim % chi-square statistics taking the model uncertainties into account define pulsarorbit_fit_constraint(stat, pars) { - variable dmdl = array_flatten(array_struct_field(array_map( - Struct_Type, &get_dataset_metadata, all_data - ), "errorprop")); - variable data = array_map( - Struct_Type, &get_data_counts, all_data - ); - variable err = array_flatten(array_struct_field(data, "err")); - data = array_flatten(array_struct_field(data, "value")); - variable mdl = array_flatten(array_struct_field(array_map( - Struct_Type, &get_model_counts, all_data - ), "value")); - return sum(sqr(data - mdl) / (sqr(err) + sqr(dmdl))); + vmessage("error (%s): is being re-implemented", _function_name); + return; + % variable dmdl = array_flatten(array_struct_field(array_map( + % Struct_Type, &get_dataset_metadata, all_data + % ), "errorprop")); + % variable data = array_map( + % Struct_Type, &get_data_counts, all_data + % ); + % variable err = array_flatten(array_struct_field(data, "err")); + % data = array_flatten(array_struct_field(data, "value")); + % variable mdl = array_flatten(array_struct_field(array_map( + % Struct_Type, &get_model_counts, all_data + % ), "value")); + % return sum(sqr(data - mdl) / (sqr(err) + sqr(dmdl))); } -- GitLab