From 90ceda3fe779a27ffebb05c5d9d38551cbed27f0 Mon Sep 17 00:00:00 2001 From: Thomas Dauser Date: Thu, 17 Oct 2019 18:43:03 +0200 Subject: [PATCH] Revert "Merge branch 'master' of http://www.sternwarte.uni-erlangen.de/gitlab/remeis/isisscripts" This reverts commit 24790ebc2df35350b7dc3777b0b4d4acb4884f23 --- .gitlab-ci.yml | 29 +- src/Mike/sitar.sl | 2 +- src/data/GTI/pulsar_GTI.sl | 256 ++++++++---------- src/data/binning/optimal_binning.sl | 159 ----------- src/data/data_response.sl | 74 ----- src/data/pulsar/pulseprofile_energy_map.sl | 1 - src/data/sys_err_array2string.sl | 7 - src/fitting/fit_pars.sl | 12 +- src/histogram/hist_fwhm.sl | 151 ----------- src/misc/astronomy/hardnessratio.sl | 6 +- src/misc/atomic/atom_name.sl | 43 ++- src/misc/atomic/element_name.sl | 172 +++--------- src/misc/sltable.sl | 2 +- src/misc/tcp.sl | 4 +- .../geometry => plot}/simplify_polygon.sl | 0 src/plot/xfig_plot_conf.sl | 17 +- src/sixte/fits_write_sixte_vignetting.sl | 7 +- src/slang/astronomy/coco.sl | 52 ++-- src/slang/astronomy/deg2dms.sl | 17 +- src/slang/geometry/point_in_polygon.sl | 209 -------------- src/slang/geometry/point_is_left_of_line.sl | 37 --- src/slang/math/RD2rad.sl | 4 +- src/slang/math/angular_separation.sl | 46 ---- src/slang/math/greatcircle_distance.sl | 1 - src/slang/strings/Roman.sl | 165 +++-------- src/slang/strings/angle2string.sl | 2 +- src/slang/strings/generate_iauname.sl | 51 ---- src/slang/strings/strsplit.sl | 36 --- src/slang/structs/sort_struct_fields.sl | 29 -- 29 files changed, 292 insertions(+), 1299 deletions(-) delete mode 100644 src/data/binning/optimal_binning.sl delete mode 100644 src/data/data_response.sl delete mode 100644 src/histogram/hist_fwhm.sl rename src/{slang/geometry => plot}/simplify_polygon.sl (100%) delete mode 100644 src/slang/geometry/point_in_polygon.sl delete mode 100644 src/slang/geometry/point_is_left_of_line.sl delete mode 100644 src/slang/math/angular_separation.sl delete mode 100644 src/slang/strings/generate_iauname.sl delete mode 100644 src/slang/strings/strsplit.sl delete mode 100644 src/slang/structs/sort_struct_fields.sl diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f047025a..31546606 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -11,7 +11,7 @@ stages: - build - test - deploy -# - wiki + - wiki build: stage: build @@ -20,7 +20,6 @@ build: - make help - make wiki only: - - merge_requests - master artifacts: expire_in: 5 min @@ -34,7 +33,7 @@ test: script: - make test only: - - merge_requests + - master provide: stage: deploy @@ -53,15 +52,15 @@ provide: # clone wiki with pipeline and add tm expanded help files # the markdown translation is defined in doc/doc/ISISscripts_md.tm -# wiki: -# stage: wiki -# script: -# - git clone https://stierhof:${WIKI_ACCESS}@www.sternwarte.uni-erlangen.de/gitlab/remeis/isisscripts.wiki.git -# - ./bin/collect_functions src | ./bin/software_structure share/Function-reference 1 share/isisscripts.md -# - rsync -r share/Function-reference isisscripts.wiki -# - mv share/Function-reference.md isisscripts.wiki -# - cd isisscripts.wiki -# - if ! [ -z "$(git status --porcelain)" ]; then git add -A; git commit -m "[AUTO] wiki update"; git push origin master; fi -# - cd ../ && rm -rf isisscripts.wiki -# only: -# - master \ No newline at end of file +wiki: + stage: wiki + script: + - git clone https://stierhof:${WIKI_ACCESS}@www.sternwarte.uni-erlangen.de/gitlab/remeis/isisscripts.wiki.git + - ./bin/collect_functions src | ./bin/software_structure share/Function-reference 1 share/isisscripts.md + - rsync -r share/Function-reference isisscripts.wiki + - mv share/Function-reference.md isisscripts.wiki + - cd isisscripts.wiki + - if ! [ -z "$(git status --porcelain)" ]; then git add -A; git commit -m "[AUTO] wiki update"; git push origin master; fi + - cd ../ && rm -rf isisscripts.wiki + only: + - master \ No newline at end of file diff --git a/src/Mike/sitar.sl b/src/Mike/sitar.sl index 2b55e1d8..68c01d00 100644 --- a/src/Mike/sitar.sl +++ b/src/Mike/sitar.sl @@ -2271,7 +2271,7 @@ public define sitar_avg_cpd() %\qualifiers{ %\qualifier{dt}{ [Length of evenly spaced bins. Default == 1.]} %\qualifier{times}{ [Times of measurements (otherwise presumed to have no gaps).]} -%\qualifier{norm}{ [Determine PSD normalization convention. =1, 'Leahy normalization'. Poisson noise PSD == 2, in absence of deadtime effects. !=1, rms or Belloni-Hasinger or Miyamoto normalization, i.e., PSD == (rms)^2/Hz & Poisson noise== 2/Rate.]} +%\qualifier{norm}{ [Determine PSD normalization convention. =1, 'Leahy normalization'. Poisson noise PSD == 2, in absence of deadtime effects). !=1, rms or Belloni-Hasinger or Miyamoto normalization, i.e., PSD == (rms)^2/Hz & Poisson noise== 2/Rate.]} %} %\description % Take an evenly spaced lightcurve (presumed counts vs. time), and diff --git a/src/data/GTI/pulsar_GTI.sl b/src/data/GTI/pulsar_GTI.sl index bbdaa495..b67a096b 100644 --- a/src/data/GTI/pulsar_GTI.sl +++ b/src/data/GTI/pulsar_GTI.sl @@ -1,175 +1,151 @@ -% -*- mode:slang; mode:fold -*- - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -private define phasebin_to_gti(tstart, tstop, tref, p, pdot, pddot, phase_lo, phase_hi, satellite){ %{{{ - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - variable d2s = 86400.; - variable eph = struct{t0, p0 = p, pdot = pdot, p2dot = pddot, p3dot =0.}; % orbit structure required by pulse_time +private define phasebin_to_gti(tstart, tstop, tref, p, pdot, pddot, phase_lo, phase_hi, satellite){ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - if (qualifier_exists("MJD")){ - eph.t0 = double(tref); - tref = (tref-MJDref_satellite(satellite))*d2s; % Converting the reference time in the staellite's time format - tstart = (tstart-MJDref_satellite(satellite))*d2s; % Start and stop time are converted to seconds, because - tstop = (tstop-MJDref_satellite(satellite))*d2s; % the pulse period will always be in seconds! - - }else{ - eph.t0 = tref/d2s + MJDref_satellite(satellite); % The reference time in the ephemeris structure must always be in MJD - } - variable nxp = [ int((tstart-tref)/p)-2 : int((tstop-tref)/p) +2] ; %index array of pulses during observation - - variable gti = struct { - start = Double_Type[length(nxp)] , - stop = Double_Type[length(nxp)] - }; - - %calculate pulse arrival times - gti.start = pulse_time (nxp + phase_lo; eph = eph); % pulse_time returns arrival times in MJD - gti.stop = pulse_time (nxp + phase_hi; eph = eph); % Make sure the MJD qualifier is NOT passed to pulse time!!! - - %convert to satellite's time system - gti.start = (gti.start - MJDref_satellite(satellite))*d2s; % convert back to seconds - gti.stop = (gti.stop - MJDref_satellite(satellite))*d2s; +variable eph = struct{t0 = double(tref), p0 = p, pdot = pdot, p2dot = pddot, p3dot =0.}; % orbit structure required by pulse_time - return gti; +tref = (tref-MJDref_satellite(satellite))*86400; %converting the reference time in the staellite's time format + +% if(tref>tstart){ +% variable n= int((tref-tstart)/p +0.5); +% tref = tref - n*p - 0.5*n*n*p*pdot - 1./6.*sqr(n)*n*sqr(p)*pddot; %make sure the reference time is before the first event +% } + +variable nxp = [ int((tstart-tref)/p)-2 : int((tstop-tref)/p) +2] ; %index array of pulses during observation + +variable bargti = struct {start = Double_Type[length(nxp)] , stop = Double_Type[length(nxp)]}; +variable k; + +%calculate pulse arrival times +bargti.start = pulse_time (nxp + phase_lo; eph = eph); +bargti.stop = pulse_time (nxp + phase_hi; eph = eph); + +%convert to satellite's time system +bargti.start = (bargti.start - MJDref_satellite(satellite))*86400; +bargti.stop = (bargti.stop - MJDref_satellite(satellite))*86400; + +return bargti; } -%}}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -private define convert_GTI(nobarevt, barevt, gti){ %{{{ +private define convert_GTI(nobarevt, barevt, gti){ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Do linear interpolation to calculate the GTIs in the satellite's time format - variable nobargti = struct{ - start = interpol (gti.start, barevt.time, nobarevt.time), - stop = interpol (gti.stop, barevt.time, nobarevt.time)}; +%Check the arguments + if (length(nobarevt.time) != length(barevt.time)){ + () = printf("ERROR: Event files do not contain an equal number of events! \n"); + return -1; + }; - return nobargti; +% Do linear interpolation to calculate the GTIs in the satellite's time format +variable nobargti = struct{start = interpol (gti.start, barevt.time, nobarevt.time), stop = interpol (gti.stop, barevt.time, nobarevt.time)}; +return nobargti; }; -%}}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -define pulsar_GTI(){ %{{{ +define pulsar_GTI(){ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %!%+ %\function{pulsar_GTI} -%\synopsis{Calculates GTIs for pulse phase resolved spectroscopy} -%\usage{pulsar_GTI(Double_Type tstart, tstop, String_Type satellite, basefilename, Double_Type t0, p, phase_lo, phase_hi);} +%\synopsis{Calculates GTIs for pulse phase resolved spectroscopy in the satellite's local time system} +%\usage{pulsar_GTI(String_Type local_evts, bary_evts, basefilename, Double_Type t0, p, phase_lo, phase_hi);} %\qualifiers{ +%\qualifier{satellite}{: satellite to set internal reference time (default: Suzaku) } %\qualifier{pdot}{: first derivative of the pulse period in s/s (default: 0) } %\qualifier{pddot}{: second derivative of the pulse period in s/s^2 (default: 0) } -%\qualifier{MJD}{: If set, tstart, tstop and t0 are assumed to in MJD} -%\qualifier{local}{: create GTIs in satellite's local time system. Requires qualifiers 'nobarevt' and 'barevt' to be set.} -%\qualifier{barevt}{: event file in barycentered time frame} -%\qualifier{nobarevt}{: event file in local time frame} %\qualifier{}{all other qualifiers are passed to BinaryPos} %} %\description -% This function calculates the GTIs for pulse phase resolved spectroscopy. Input arguments are\n -% -% tstart - start time of the time interval to be covered (typically the start of the observation) -% tstart - end time of the time interval to be covered (typically the end of the observation) -% satellite - the name of the satellite (needed for setting the reference MJD correctly) -% basefilename - the file created will be named 'basefilename.gti' -% t0 - the reference time (time of phase zero) -% p - the pulse period -% phase_lo - the lower phase boundary -% phase_hi - the upper phase boundary -% \n -% Upon qualifier request, tstart, tstop, and t0 can be given in MJD instead of seconds. The pulse -% period is always in seconds.\n +% This function calculates the GTIs for pulse phase resolved spectroscopy +% in the barycentric time system and interpolates them to the satellite's +% (unless specified: Suzaku) local time system, as required for some +% extraction routines. The arguments to be given are an event file, once +% in the barycentric time system and once in the local time system. These +% event files can be obtained, e.g., by using the specific keywords in +% the Remeis extraction scripts or by applying barycentric correction to +% the raw data. Alternatively, the event files can also be passed as +% structures. The GTIs are written to a FITS file, whose name is specified +% in \code{basefilename}. The pulse period \code{p} is to be given in +% seconds, the reference time \code{t0} (i.e., when the pulse phase is +% equal to 0) in MJD. The limits of the phase interval are given as +% \code{phase_lo} and \code{phase_hi}. The total length of the GTI file is +% determined by the first and last event of the event files.\n % \n % If the pulse period and the reference time are corrected for the binary % motion, the orbital parameters should be provided by qualifiers such that % the GTIs can be transformed into the observers time system. The % qualifiers are passed and equal to the \code{BinaryPos} function. % \n -% If GTIs in the satellite's local reference time system are needed (e.g., for Suzaku-XIS), -% an event file in both local and barycenterd time can be passed via qualifiers and the GTIs -% are interpolated to the local time system. Header keywords are set accordingly. +% NOTE: This tool is still under development. Therefore a thorough +% check of the created GTIs is recommended. Please report any problems to +% ralf.ballhausen@sternwarte.uni-erlangen.de +% %\example -% pulsar_GTI (2.7082e8, 2.7097e8, "nustar", "phase_0.25-0.35", 2.708242e8, 443.07, 0.25, 0.35); +% pulsar_GTI("nobary_xis3_3x3.evt","bary_xis3_3x3.evt","GTI_file", 56618.6, 18.8, 0.25, 0.5 +% ; satellite = "Suzaku");\n % \n -% creates a GTI file named 'phase_0.25-0.35.gti', ranging from 2.7082e8--2.7097e8 seconds -% in NuSTAR's mission specific time system for the pulsar 4U 1907+09 with pulse period 443.07 seconds -% for the pulse phase interval 0.25--0.35 where phase 0.0 is set to be at 2.708242e8 seconds. +% creates a GTI file named 'GTI_file.gti', based on the two event files +% 'nobary_xis3_3x3.evt' and 'bary_xis3_3x3.evt', assuming a reference +% time of MJD 56618.6 and a pulse period of 18.8 s for the phase +% interval 0.25--0.5. % -%\seealso{MJDref_satellite, BinaryPos, pulse_time} +%\seealso{MJDref_satellite, BinaryPos} %!%- - variable nobarevt, barevt, name, t0, p, phase_lo, phase_hi, tstart, tstop, satellite; - switch(_NARGS) - {case 8: (tstart, tstop, satellite, name, t0, p, phase_lo, phase_hi) = (); } + variable nobarevt, barevt, name, t0, p, phase_lo, phase_hi; + switch(_NARGS) + {case 7: (nobarevt, barevt, name, t0, p, phase_lo, phase_hi) = (); } {help(_function_name()); return; } + variable satellite = qualifier("satellite", "Suzaku"); + variable pdot = qualifier("pdot", 0.); + variable pddot = qualifier("pddot", 0.); - variable fn_local = qualifier("nobarevt", NULL); - variable local = qualifier_exists("local"); - variable fn_bary = qualifier("barevt", NULL); - variable pdot = qualifier("pdot", 0.); - variable pddot = qualifier("pddot", 0.); - variable d2s = 86400.; % days to seconds - - if (local){ - if ( fn_local == NULL or fn_bary == NULL){ - fprintf(stderr, "%s: ERROR: Conversion to LOCAL time system requires qualifiers 'nobarevt' and 'barevt' to be set!\n", _function_name); - return; - } - %Check the arguments - if(typeof(fn_local) == String_Type){ - if(fits_read_key(fn_local, "TIMEREF") != "LOCAL"){ - ()= printf("ERROR: Event file has wrong reference time! \n"); - return -1; - }; - nobarevt = fits_read_table(fn_local); - }; - if(typeof(fn_bary) == String_Type){ - if(fits_read_key(fn_bary, "TIMEREF") != "SOLARSYSTEM"){ - ()=printf("ERROR: Event file has wrong reference time! \n"); - return -1; - }; - barevt = fits_read_table(fn_bary); - }; - %Check the arguments - if (length(nobarevt.time) != length(barevt.time)){ - () = printf("ERROR: Event files do not contain an equal number of events! \n"); - return -1; - }; - - - struct_filter(nobarevt, array_sort(barevt.time)); - struct_filter(barevt, array_sort(barevt.time)); - } - - variable bargti = phasebin_to_gti(tstart, tstop, t0, p, pdot, pddot, phase_lo, phase_hi, satellite;;__qualifiers); - - % binary correction - if (qualifier_exists("asini") and qualifier_exists("porb")) { - variable qual = reduce_struct(__qualifiers, "satellite"); - % Do the binary correction in MJD - bargti.start = bargti.start/d2s + MJDref_satellite(satellite); - bargti.stop = bargti.stop/d2s + MJDref_satellite(satellite); - - bargti.start += BinaryPos(bargti.start;; qual)/d2s; - bargti.stop += BinaryPos(bargti.stop ;; qual)/d2s; - - bargti.start = (bargti.start - MJDref_satellite(satellite))*d2s; - bargti.stop = (bargti.stop - MJDref_satellite(satellite))*d2s; - } - - if (local){ - variable nobargti = convert_GTI(nobarevt, barevt, bargti); - } - - % Write the FITS file and update some keywords - fits_write_gti(name +".gti", local ? nobargti : bargti, MJDref_satellite(satellite); date = time() ); - fits_update_key(name +".gti", "TELESCOP", satellite, "Mission name"); - fits_update_key(name +".gti", "TIMEREF", local ? "LOCAL" : "SOLARSYSTEM"); - fits_update_key(name +".gti[1]", "TIMEREF", local ? "LOCAL" : "SOLARSYSTEM"); - fits_update_key(name +".gti", "MJDREFI", int(MJDref_satellite(satellite)), "MJD reference day"); - fits_update_key(name +".gti", "MJDREFF", MJDref_satellite(satellite) - int(MJDref_satellite(satellite)), "fractional part of the MJD reference"); - fits_update_key(name +".gti[1]", "MJDREFI", int(MJDref_satellite(satellite)), "MJD reference day"); - fits_update_key(name +".gti[1]", "MJDREFF", MJDref_satellite(satellite) - int(MJDref_satellite(satellite)), "fractional part of the MJD reference"); - fits_write_comment (name +".gti", "Phase_lo = " + string(phase_lo)); - fits_write_comment (name +".gti", "Phase_hi = " + string(phase_hi)); + %Check the arguments + if(typeof(nobarevt) == String_Type){ + if(fits_read_key(nobarevt, "TIMEREF") != "LOCAL"){ + ()= printf("ERROR: Event file has wrong reference time! \n"); + return -1; + }; + nobarevt = fits_read_table(nobarevt); + }; + if(typeof(barevt) == String_Type){ + if(fits_read_key(barevt, "TIMEREF") != "SOLARSYSTEM"){ + ()=printf("ERROR: Event file has wrong reference time! \n"); + return -1; + }; + barevt = fits_read_table(barevt); + }; + + struct_filter(nobarevt, array_sort(barevt.time)); + struct_filter(barevt, array_sort(barevt.time)); + + variable bargti = phasebin_to_gti(barevt.time[0], barevt.time[-1], t0, p, pdot, pddot, phase_lo, phase_hi, satellite); + variable nobargti = convert_GTI(nobarevt, barevt, bargti); + + % binary correction + if (qualifier_exists("asini") and qualifier_exists("porb")) { + variable qual = reduce_struct(__qualifiers, "satellite"); + % Do the binary correction in MJD + nobargti.start = nobargti.start/86400. + MJDref_satellite(satellite); + nobargti.stop = nobargti.stop/86400. + MJDref_satellite(satellite); + + nobargti.start += BinaryPos(nobargti.start;; qual)/86400.; + nobargti.stop += BinaryPos(nobargti.stop ;; qual)/86400.; + + nobargti.start = (nobargti.start - MJDref_satellite(satellite))*86400.; + nobargti.stop = (nobargti.stop - MJDref_satellite(satellite))*86400; + } + + + %Write the FITS file and update some keywords + fits_write_gti(name +".gti", nobargti, MJDref_satellite(satellite); date = time() ); + fits_update_key(name +".gti", "TELESCOP", satellite, "Mission name"); + fits_update_key(name +".gti", "TIMEREF", "LOCAL"); + fits_update_key(name +".gti[1]", "TIMEREF", "LOCAL"); + fits_update_key(name +".gti", "MJDREFI", int(MJDref_satellite(satellite)), "MJD reference day"); + fits_update_key(name +".gti", "MJDREFF", MJDref_satellite(satellite) - int(MJDref_satellite(satellite)), "fractional part of the MJD reference"); + fits_update_key(name +".gti[1]", "MJDREFI", int(MJDref_satellite(satellite)), "MJD reference day"); + fits_update_key(name +".gti[1]", "MJDREFF", MJDref_satellite(satellite) - int(MJDref_satellite(satellite)), "fractional part of the MJD reference"); + fits_write_comment (name +".gti", "Phase_lo = " + string(phase_lo)); + fits_write_comment (name +".gti", "Phase_hi = " + string(phase_hi)); } -%}}} \ No newline at end of file diff --git a/src/data/binning/optimal_binning.sl b/src/data/binning/optimal_binning.sl deleted file mode 100644 index 1883ed63..00000000 --- a/src/data/binning/optimal_binning.sl +++ /dev/null @@ -1,159 +0,0 @@ -%%% functions for model and data binning following certain criteria -%%% -%%% - -%%% optimal binning after Kaastra & Bleeker 2016 - -%%% -private define binWidthKaastraBleeker (numCounts, numResolver) -%%% -%%+ -%\function{binWidthKaastraBleeker} -%\synopsis{Caluculate optimized bin width per data resolution} -%\usage{Double_Type[] binWidthKaastraBleeker (Double_Type[] numCounts, -% Double Type numResolver)} -%\description -% Kaaster & Bleeker (2016) describe an aproximation to data binning to -% minimize the number of bins by keeping the resolving power for near -% resolution limit features. -% -% The approximation is based on a Monte Carlo approach to band limited -% data. The result gives the optimized bin width per FWHM for each data -% channel. It depends on the number of counts for each channel in the -% FWHM of the instrument and the number of resolution elemets -% (roughly sum 1/c_i where c_i is the FWHM in channels for -% channel i) -%\seealso{} -%%- -{ - variable r = 1 + .2*log(numResolver); - variable t = log(numCounts*r); - variable d = Double_Type[length(numCounts)]+1.; - variable n = where(t>2.119); - d[n] = (0.08 + 7./t[n] + 1.8/(t[n]^2))/(1+5.9/t[n]); - return d; -} -%%% - -%%% -private define FWHMChannelRange (response) -%%% -{ - variable rLen = length(response); - variable i; - variable c1 = Double_Type[rLen]; - variable c2 = Double_Type[rLen]; - _for i (0, rLen-1, 1) - (c1[i], c2[i]) = hist_fwhm_index(response[i], i; interp); - - return (c1, c2); -} - -%%% -private define countsInResolver (c1, c2, response, counts) -%%% -%%+ -%\function{countsInResolver} -%\synopsis{Calculate approximated number of counts per FWHM} -%\usage{Double_Type[] countsInResolver (Double_Type[] c1, -% Double_Type[] c2, -% Double_Type[][] response, -% Double_Type[] counts)} -%\description -% -%\seealso{} -%%- -{ - variable size = length(response); - variable h = Double_Type[size], c = Double_Type[size]; - variable i1, i2; - variable i,array,len; - _for i (0, size-1, 1) { - i1 = int(c1[i]); - i2 = int(c2[i]); - array = response[i]; - len = length(array); - c[i] = sum(counts[[i1:i2]]); - h[i] = sum(array)/sum(array[[(i1-((i1 == 0)?0:1)):(i2+((i2==len-1)?0:1))]]); - } - return c*h; -} - -%%% -private define buildOptimalBins (response, counts) -%%% -%%+ -%\function{buildOptimalBins} -%\synopsis{Calculate optimal binning (index)} -%\usage{(Int_Type[], Int_Type[]) buildOptimalBins (Double_Type[][] response, -% Double_Type[] counts)} -%\seealso{} -%%- -{ - variable c1,c2; - variable resWidth; - variable numResolver; - variable cnts; - variable channelWidth; - variable i, len, mini, idx; - variable loIdx, hiIdx; - - (c1,c2) = FWHMChannelRange(response); - resWidth = c2-c1; - numResolver = sum(1./resWidth); - cnts = countsInResolver(c1, c2, response, counts); - channelWidth = int(binWidthKaastraBleeker(cnts, numResolver)*resWidth); - channelWidth[where(channelWidth<1)] = 1; - len = length(counts); - loIdx = {}; - hiIdx = {}; - - i = 0; - while (i < len) { - idx = [i:int(_min(i+channelWidth[i]-1, len-1))]; - mini = int(min(idx+channelWidth[idx])); - list_append(loIdx, i); - list_append(hiIdx, mini); - i = mini; - } - hiIdx[-1] = len-1; % last index can only be size of array - return list_to_array(loIdx), list_to_array(hiIdx); -}; - -%%% -define rebin_dataset_optimal (histIdx) -%%% -%!%+ -%\function{rebin_dataset_optimal} -%\usage{rebin_dataset_optimal (Int_Type index) -%\altusage{rebin_dataset_optimal (Int_Type[] index)}} -%\synopsis{Rebin datasets addressed by based on Kaastra & Bleeker 2016} -%\seealso{rebin_data} -%!%- -{ - variable response; - variable counts; - variable idx; - variable idxLo; - variable idxHi; - variable i; - variable flip; - variable mask; - - foreach idx (histIdx) { - rebin_data(idx, 0); % reset binning, ugly, yes, but necessary - response = get_dataset_response(idx); - counts = get_data(idx).value; - - (idxLo, idxHi) = buildOptimalBins (response, counts); %idxLo and idxHi correspond to energy grid - - mask = Int_Type[length(counts)]; - flip = -1; - _for i (0, length(idxLo)-1, 1) { - mask[[idxLo[i]:idxHi[i]]] = flip; - flip *= -1; - } - - rebin_data(idx, mask); % The binning mask has again applied on a wavelength grid - } -} diff --git a/src/data/data_response.sl b/src/data/data_response.sl deleted file mode 100644 index 8a198fea..00000000 --- a/src/data/data_response.sl +++ /dev/null @@ -1,74 +0,0 @@ - -%%%%%%%%%%%%%%%%%%%%%%%%%%% -define get_dataset_response (id) -%%%%%%%%%%%%%%%%%%%%%%%%%%% -%!%+ -%\function{get_dataset_response} -%\usage{Response = get_dataset_response(id);} -%\synopsis{Compute response of dataset associated to 'id' to delta peak} -%\qualifiers{ -%\qualifier{energy}{: If given, return response on energy grid (default: wavelength)} -%} -%\description -% This function evaluates a delta peak function set at each nominal -% grid value. The returnend array contains the response for each peak -% position. Per default the response is returned as 'seen' by ISIS, that -% means on a wavelength grid. Use the 'energy' qualifier to revers the -% response. -% -%\seealso{get_rmf_data_grid, assign_rmf, assign_rsp} -%!%- -{ - % check if dataset exists - if (all(id != all_data())) - throw UsageError, sprintf("Dataset %d not defined", id); - - variable rmfId = get_data_info(id).rmfs[0]; % get first rmf id from dataset (see meanGridValues) - variable grid; % grid from rmf - variable saveIncludeData = all_data(1); % save current excluded/included data - variable saveNoticeBins = get_data_info(id).notice_list; % save noticed bins for dataset - variable saveFitFun = get_fit_fun(); % save current fit fun to recover if necessary - variable meanGridValues; % mean value from the first rmf grid (this is the one used for the data, see assign_rsp) - variable i; % running index - variable fitObject; % the fit object - variable responseResult; % calculated responses - variable energy = qualifier_exists("energy"); % reverse returning array so it matches energy grid - - % prepare delta eval - exclude(all_data()); % exclude everything but the data of interest - include(id); - notice(id); % notice all data points - fit_fun("delta(1)"); - - % what is the (unmodified) data grid when more than one rmf is in use? - % take get_data instead? - - % collect all grid means - grid = get_rmf_data_grid(rmfId); - if (NULL == grid) % with define_counts get_rmf_data_grid returns NULL use get_data in any case? - grid = get_data(id); - meanGridValues = .5*(grid.bin_hi+grid.bin_lo); - - responseResult = Array_Type[length(meanGridValues)]; - - fitObject = open_fit(); % open function - _for i (0, length(responseResult)-1, 1) { - () = fitObject.eval_statistic([1,meanGridValues[i]]); % eval function with norm=1 and lambda=mean - if (energy) - responseResult[i] = reverse(get_model(id).value); - else - responseResult[i] = get_model(id).value; - } - fitObject.close(); - - % rest to initial settings - notice_list(id,saveNoticeBins); - exclude(all_data()); - include(saveIncludeData); - fit_fun(saveFitFun); - - if (energy) - return reverse(responseResult); - else - return responseResult; -} diff --git a/src/data/pulsar/pulseprofile_energy_map.sl b/src/data/pulsar/pulseprofile_energy_map.sl index d7fcec5e..5ea06386 100644 --- a/src/data/pulsar/pulseprofile_energy_map.sl +++ b/src/data/pulsar/pulseprofile_energy_map.sl @@ -3,7 +3,6 @@ require("gsl","gsl"); require("xfig"); require("gcontour"); -require("png"); % TODO % - detailed help is missing for all functions! diff --git a/src/data/sys_err_array2string.sl b/src/data/sys_err_array2string.sl index 7e84df81..4655cdbd 100644 --- a/src/data/sys_err_array2string.sl +++ b/src/data/sys_err_array2string.sl @@ -19,7 +19,6 @@ define sys_err_array2string(arr) variable i_end = 0; variable n = length(arr), i; - if(n) { _for i(0,n-1,1) { @@ -41,11 +40,5 @@ define sys_err_array2string(arr) } % write the last line, too str_sys += delim+string(i_beg)+"-"+string(i_end)+":"+string(arr[i_beg]); - } else { - % if length(arr) is zero, there are no systematic errors defined. This - % string should result in zero systematic errors when run through - % string2sys_err_array() and set_sys_err_frac(). - str_sys = "0-0:0"; - } return str_sys; } diff --git a/src/fitting/fit_pars.sl b/src/fitting/fit_pars.sl index 787e8f02..99b76f7f 100644 --- a/src/fitting/fit_pars.sl +++ b/src/fitting/fit_pars.sl @@ -14,9 +14,6 @@ define fit_pars() %\qualifier{level}{[=1]: specifies the confidence level. Values of 0, 1, or 2 % indicate 68%, 90%, or 99% confidence levels respectively. % By default, 90% confidence limits are computed.} -%\qualifier{chi2diff}{: if given, fit_pars invokes fconf instead of -% conf with the value given to caluculate confidence in chi2 -% range. Defaults to 1.0.} %\qualifier{tolerance}{: convergence criterion for the calculation of the confidence % limits (see help for the conf command). Default: 1e-3} %} @@ -28,7 +25,7 @@ define fit_pars() % \code{buf_below} (\code{buf_above}) is the fraction of the allowed range \code{[min:max]} % which separates the lower (upper) confidence limit from \code{min} (\code{max}). % If one of these buffers is 0, your confidence interval has bounced. -%\seealso{pvm_fit_pars, conf, fconf} +%\seealso{pvm_fit_pars, conf} %!%- { variable pars; @@ -42,7 +39,6 @@ define fit_pars() variable basefilename = qualifier("basefilename", strftime("%Y-%m-%d_%H:%M:%S", localtime(_time())))+ "_"; variable level = qualifier("level",1); variable tolerance = qualifier("tolerance",1e-3); - variable chi2diff = qualifier("chi2diff", 1.); if (not(level==0 or level==1 or level==2)) { message("requested confidence level not available. Default to level =1: 90%\n"); level=1; @@ -76,11 +72,7 @@ define fit_pars() { ok = 1; for(i=0; i maxIndex) - throw UsageError, "given index not in range of array"; - - newMax = 0; - loBounced = 0; - hiBounced = 0; - while (not loBounced || not hiBounced) { - if (not loBounced && iminus>0) - iminus--; - - if (not hiBounced && iplus < histLen-1) - iplus++; - - if (iminus == 0 || Histogram[iminus] <= .5*fullMax) - loBounced = 1; - - if ((iplus == (histLen-1)) || Histogram[iplus] <= .5*fullMax) - hiBounced = 1; - - if (Histogram[iminus] > fullMax) { - fullMax = Histogram[iminus]; - maxIndex = iminus; - newMax = 1; - } - - if (Histogram[iplus] > fullMax) { - fullMax = Histogram[iplus]; - maxIndex = iplus; - newMax = 1; - } - - if (newMax) { - hiBounced = 0; - loBounced = 0; - - while (Histogram[iplus] < .5*fullMax) - iplus--; - - while (Histogram[iminus] < .5*fullMax) - iminus++; - - newMax = 0; - } - } - - if (interp) { - loSlope = Histogram[iminus+1]-Histogram[iminus]; - hiSlope = Histogram[iplus]-Histogram[iplus-1]; - iminusinterp = 1.*iminus + ((iminus == 0) ? 0. : .5*fullMax/loSlope - Histogram[iminus]/loSlope); - iplusinterp = 1.*iplus + ((iplus == histLen-1) ? 0. : .5*fullMax/hiSlope - Histogram[iplus]/hiSlope); - iminus = (loSlope == 0) ? iminus : iminusinterp; - iplus = (hiSlope == 0) ? iplus : iplusinterp; - } - - return (iminus, iplus); -} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -define get_instrument_resolution_from_data (id) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%!%+ -%\function{get_instrument_resolution_from_data} -%\usage{R = get_instrument_resolution_from_data(id)} -%\synopsis{Compute energy resolution from data set 'id'} -%\qualifiers{ -%\qualifier{energy}{: If given will return the grid as energy grid -% instead of wavelength} -%} -%\description -% This function can be used to estimate the instrument resolution -% from the data set specified with id. The return value is a struct -% { value, bin_lo, bin_hi } where bin_lo and bin_hi is the data -% grid (usually in Angstrom) and value is the resolution in keV. -% -% Note that this function removes any binning from the data set. -% -%\seealso{hist_fwhm_index} -%!%- -{ - rebin_data(id, 0); - variable response = get_dataset_response(id); - variable cmin, cmax; - variable resolution = Double_Type[length(response)]; - variable i; - variable grid = get_rmf_data_grid(get_data_info(id).rmfs[0]); - variable fmin, fmax; - variable energy = qualifier_exists("energy"); - - _for i (0, length(response)-1, 1) { - (cmin,cmax) = hist_fwhm_index(response[i], i); - fmin = (response[i][cmin+1] == response[i][cmin]) ? 0.0 : - (max(response[i][[cmin:cmax]])*.5-response[i][cmin])/ - (response[i][cmin+1]-response[i][cmin]); - fmax = (response[i][cmax] == response[i][cmax-1]) ? 0.0 : - (max(response[i][[cmin:cmax]])*.5-response[i][cmax])/ - (response[i][cmax]-response[i][cmax-1]); - resolution[i] = _A(grid.bin_lo[cmin] + (grid.bin_lo[cmin+1]-grid.bin_lo[cmin])*fmin) - -_A(grid.bin_hi[cmax] - (grid.bin_hi[cmax]-grid.bin_hi[cmax-1])*fmax); - } - - if (energy) - return struct { value = resolution, bin_lo = _A(grid.bin_hi), bin_hi = _A(grid.bin_lo) }; - else - return struct { value = resolution, @grid }; -} diff --git a/src/misc/astronomy/hardnessratio.sl b/src/misc/astronomy/hardnessratio.sl index 70910f04..559a245e 100644 --- a/src/misc/astronomy/hardnessratio.sl +++ b/src/misc/astronomy/hardnessratio.sl @@ -10,8 +10,8 @@ define hardnessratio() % \qualifier{back_s}{: background counts in the soft band} % \qualifier{back_h}{: background counts in the hard band} % \qualifier{backscale}{: ratio between the extraction regions for the source and the background} -% \qualifier{exposure}{: Exposure per bin in seconds. If given, interpret soft_count, hard_count and t -% he backgrounds as rates. Multiply the rates with exposure to get the counts.} +% \qualifier{exposure}{: if given, interpret soft_count, hard_count and the backgrounds as +% rates. Multiply the rates with exposure to get the counts.} % \qualifier{backexposure}{: if given, the background exposure. Only taken into % account if exposure is also set. If not given, it is assumed that % the source and background exposures are identical.} @@ -71,7 +71,7 @@ define hardnessratio() s=s*expo; h=h*expo; - variable backexpo=qualifier("backexposure",expo); + variable backexpo=("backexposure",expo); bs=bs*backexpo; bh=bh*backexpo; diff --git a/src/misc/atomic/atom_name.sl b/src/misc/atomic/atom_name.sl index 8a0983c1..7a8f8d66 100644 --- a/src/misc/atomic/atom_name.sl +++ b/src/misc/atomic/atom_name.sl @@ -1,13 +1,44 @@ define atom_name() %!%+ %\function{atom_name} -%\synopsis{DEPRECATED} +%\synopsis{returns the symbol for atoms with proton number Z} %\usage{String_Type atom_name(Integer_Type Z)} -%\description -% This function returns the symbol for atoms with proton number Z. -% This function is DEPRECATED, please use element_symbol. -%\seealso{element_symbol} +%\seealso{element_name} %!%- { - return element_symbol(();;__qualifiers); + variable Z; + switch(_NARGS) + { case 1: Z = (); } + { help(_function_name()); return; } + + switch(Z) + { case 1: return "H"; } + { case 2: return "He"; } + { case 3: return "Li"; } + { case 4: return "Be"; } + { case 5: return "B"; } + { case 6: return "C"; } + { case 7: return "N"; } + { case 8: return "O"; } + { case 9: return "F"; } + { case 10: return "Ne"; } + { case 11: return "Na"; } + { case 12: return "Mg"; } + { case 13: return "Al"; } + { case 14: return "Si"; } + { case 15: return "P"; } + { case 16: return "S"; } + { case 17: return "Cl"; } + { case 18: return "Ar"; } + { case 19: return "K"; } + { case 20: return "Ca"; } + { case 21: return "Sc"; } + { case 22: return "Ti"; } + { case 23: return "V"; } + { case 24: return "Cr"; } + { case 25: return "Mn"; } + { case 26: return "Fe"; } + { case 27: return "Co"; } + { case 28: return "Ni"; } + "Don't know, please tell me."; } diff --git a/src/misc/atomic/element_name.sl b/src/misc/atomic/element_name.sl index bd06d094..3d7c1834 100644 --- a/src/misc/atomic/element_name.sl +++ b/src/misc/atomic/element_name.sl @@ -1,150 +1,44 @@ - -private variable _element_sym = ["Bare","H","He","Li","Be","B","C","N","O","F","Ne", - "Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc", - "Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge", - "As","Se","Br","Kr", "Rb","Sr","Y","Zr","Nb","Mo","Tc", - "Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe", - "Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb", - "Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os", - "Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn", - "Rr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk", - "Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs", - "Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"]; - -private variable _element_name = ["Bare","Hydrogen","Helium","Lithium","Beryllium","Boron", - "Carbon","Nitrogen","Oxygen","Fluorine","Neon", - "Sodium","Magnesium","Aluminum","Silicon","Phosphorus", - "Sulfur","Chlorine","Argon","Potassium","Calcium","Scandium", - "Titanium","Vanadium","Chromium","Manganese","Iron","Cobalt", - "Nickel","Copper","Zinc","Gallium","Germanium", - "Arsenic","Selenium","Bromine","Krypton","Rubidium","Strontium", - "Yttrium","Zirconium","Niobium","Molybdenum","Technetium", - "Ruthenium","Rhodium","Palladium","Silver","Cadmium","Indium", - "Tin","Antimony","Tellurium","Iodine","Xenon", - "Cesium","Barium","Lanthanum","Cerium","Praseodymium", - "Neodymium","Promethium","Samarium","Europium","Gadolinium", - "Terbium", - "Dysprosium","Holmium","Erbium","Thulium","Ytterbium", - "Lutetium","Hafnium","Tantalum","Tungsten","Rhenium", - "Osmium", - "Oridium","Platinum","Gold","Mercury","Thallium","Lead", - "Bismuth","Polonium","Astatine","Radon", - "Francium","Radium","Actinum","Thorium","Protactinium", - "Uranium","Neptunium","Plutonium","Americium","Curium", - "Berkelium", - "Californium","Einsteinium","Fermium","Mendelevium","Nobelium", - "Lawrencium","Rutherfordium","Dubnium","Seaborgium", - "Bohrium","Hassium", - "Meitnerium","Darmstadtium","Roentgenium","Copernicium", - "Nihonium","Flerovium","Moscovium","Livermorium", - "Tennessine","Oganesson"]; - -private variable _element2Z=Assoc_Type[Int_Type]; - -define element_name() { +define element_name() %!%+ %\function{element_name} %\synopsis{returns the name of the element with proton number Z} %\usage{String_Type element_name(Integer_Type Z)} -%\qualifiers{ -%\qualifier{lc}{return symbol or name in lower case (default: Capitalized)} -%} -%\description -% This function returns the name of the element with -% nuclear charge Z for all named elements. -% -% This function is array safe. -%\seealso{element_symbol,element2Z} -%!%- - variable Z; - switch(_NARGS) - { case 1: Z = (); } - { help(_function_name()); return; } - - if (min(Z)<0 or max(Z)>=length(_element_name)) { - throw UsageError,sprintf("%s: Element Z outside of 0=length(_element_sym)) { - throw UsageError,sprintf("%s: Element Z outside of 01) { if(length(col)==1) @@ -121,15 +116,15 @@ define xfig_plot_conf() _for j(0,length(gct[i].x_list)-1,1) { pl.plot(del_x*gct[i].x_list[j]+xmin, - del_y*gct[i].y_list[j]+ymin;; struct_combine(world, - struct { color=col[i], line=line[i], width=width[i] } - )); + del_y*gct[i].y_list[j]+ymin + ; color=col[i], + line=line[i], + width=width[i] + ); } } if(cross){ - pl.plot(best_x,best_y;; struct_combine(world, - struct { color=length(col) > n_levels ? col[-1] : "black", - width=width[-1], sym=4 })); + pl.plot(best_x,best_y;color=length(col) > n_levels ? col[-1] : "black",width=width[-1],sym=4); } #ifexists png_gray_to_rgb if (qualifier_exists("plot_map")) diff --git a/src/sixte/fits_write_sixte_vignetting.sl b/src/sixte/fits_write_sixte_vignetting.sl index 7ff50851..250e4965 100644 --- a/src/sixte/fits_write_sixte_vignetting.sl +++ b/src/sixte/fits_write_sixte_vignetting.sl @@ -37,7 +37,7 @@ private define get_vignet_data(nenergies, noffsetangles, nrot){ ENERGY = Float_Type[1, nenergies], THETA = Float_Type[1, noffsetangles], PHI = Float_Type[1, nrot], - VIGNET = Float_Type[1, nrot, noffsetangles, nenergies] %% definition is not compliant with OGIP, but that's how we read it in sixte ... + VIGNET = Float_Type[1, nenergies,noffsetangles,nrot] }; } @@ -74,14 +74,13 @@ define fits_write_sixte_vignetting(){ variable ii; - for (ii=0; iiP.y)) or - ((V.y[i]>P.y) and (V.y[i+1]<=P.y))) { - % edge-ray intersect coordinate - variable vt=(P.y-V.y[i])/(V.y[i+1]-V.y[i]); - % P.x < intersect x coordinate? - if (P.x < V.x[i]+vt*(V.x[i+1]-V.x[i]) ) { - cn++; - } - } - } - return cn; -} - -define winding_number_polygon() -%!%+ -%\function{winding_number_polygon} -%\usage{wn=winding_number_polygon(P,V) -%\synopsis{return the winding number for a point P in a polygon} -%\description -% return the winding number for a point P in a (potentially -% complex polygon V) -% -% The winding number counts the number of times the polygon V winds -% around point P. The point is outside if the winding number is 0. -% -% The algorithm used here is as efficient as the determination -% of the crossing number. -% -% The point P is a struct{x,y}, while the polygon is a -% struct{x[],y[]} where the arrays are the x- and y-coordinates. -% The polygon has to be closed, i.e. V.x[n]==V.x[0] and -% V.y[n]==V.y[0] where n is the number of polygon points. -% -% See the URL below for more explanations. -% -% Based on code by Dan Sunday, -% http://geomalgorithms.com/a03-_inclusion.html -% -%\seealso{crossing_number_polygon,point_in_polygon,simplify_polygon} -%!%- -{ - if (_NARGS != 2) { - help(_function_name()); - } - - variable P,V; - (P,V)=(); - - variable n=length(V.x); - - if (n<2) { - throw UsageError,sprintf("%s: Polygon must have at least 2 vertex points.\n",_function_name()); - } - - if (V.x[0]!=V.x[n-1] or V.y[0]!=V.y[n-1]) { - throw UsageError,sprintf("%s: Polygon is not closed.\n",_function_name()); - } - - if (length(V.x) != length(V.y) ) { - throw UsageError,sprintf("%s: Vertex arrays x and y must be of same length.\n",_function_name()); - } - - variable wn=0; - % loop over all edges - variable i; - _for i(0, n-2, 1) { - if (V.y[i]<=P.y) { - % start y<=P.y - % ...upward crossing? - if (V.y[i+1] > P.y) { - % Is point left of edge? - if (point_is_left_of_line(V.x[i],V.y[i],V.x[i+1],V.y[i+1],P.x,P.y) > 0) { - wn++; - } - } - } else { - % start y>P.y - % ...downward crossing? - if (V.y[i+1]<=P.y) { - if (point_is_left_of_line(V.x[i],V.y[i],V.x[i+1],V.y[i+1],P.x,P.y) < 0) { - wn--; - } - } - } - } - return wn; -} - -define point_in_polygon() -%!%+ -%\function{point_in_polygon} -%\usage{ret=point_in_polygon(p0,V);} -%\altusage{ret=point_in_polygon(x,y,Vx,Vy);} -%\synopsis{determine whether a point is in a polygon} -%\qualifiers{ -%\qualifier{evenodd}{: use the crossing number method} -%\qualifier{crossing}{: use the crossing number method} -%\qualifier{winding}{: use the winding number method (the default)} -%} -%\description -% The function returns 1 if the point p0 is located inside the -% polygon defined by the vertices V, and 0 if it is located -% outside of the polygon. -% -% The polygon has to be closed, i.e. V.x[n]==V.x[0] and -% V.y[n]==V.y[0] where n is the number of polygon points. -% -% Either the winding number (default) or the even-odd-rule -% define what is meant by inside. -% -% The point is either defined by a struct{x,y} and the vertices -% by a struct{x[],y[]}, or the coordinates can be directly -% given in the respective arrays. -% -% See the URL below for more explanations. -% -% Based on code by Dan Sunday, -% http://geomalgorithms.com/a03-_inclusion.html -% -%\seealso{crossing_number_polygon,point_in_polygon,simplify_polygon} -%!%- -% -{ - variable P,V; - if (_NARGS==2) { - (P,V)=(); - } else { - if (_NARGS==4) { - variable x,y,Vx,Vy; - (x,y,Vx,Vy)=(); - P=struct{x=x,y=y}; - V=struct{x=@Vx,y=@Vy}; - } else { - throw UsageError,sprintf("%s: Wrong number of arguments.\n",_function_name()); - } - } - if (qualifier_exists("evenodd") or qualifier_exists("crossing")) { - variable cn=crossing_number_polygon(P,V); - printf("CN: %i\n",cn); - return ( (cn mod 2)==1 ); - } - - % winding number method - variable wn=winding_number_polygon(P,V); - return ( wn!=0); -} diff --git a/src/slang/geometry/point_is_left_of_line.sl b/src/slang/geometry/point_is_left_of_line.sl deleted file mode 100644 index e718b285..00000000 --- a/src/slang/geometry/point_is_left_of_line.sl +++ /dev/null @@ -1,37 +0,0 @@ -define point_is_left_of_line() -%!%+ -%\function{point_is_left_of_line} -%\usage{ret=point_is_left_of_line(p0,p1,p2);} -%\altusage{ret=point_is_left_of_line(p0x,p0y,p1x,p1y,p2x,p2y);} -%\synopsis{determines whether point p2 is left of a line through p0 and p1} -%\description -% This function tests if point p2 is to the left of an infinite line defined -% by points p0 and p1. The points are defined by structs p=struct{x,y}, -% where x is the x-coordinate and y is the y-coordinate. -% -% The function returns -% +1 if P2 is left of the line -% 0 if P2 is on the line -% -1 -1 if P2 is right of the line -% -% Based on code by Dan Sunday, http://geomalgorithms.com/a03-_inclusion.html -% -%\seealso{crossing_number_polygon,winding_number_polygon,point_in_polygon,simplify_polygon} -%!%- -{ - variable p0,p1,p2; - variable cross; - if (_NARGS==6) { - variable p0x,p0y,p1x,p1y,p2x,p2y; - (p0x,p0y,p1x,p1y,p2x,p2y)=(); - cross= (p1x-p0x)*(p2y-p0y) - (p2x-p0x)*(p1y-p0y); - } else { - if (_NARGS==3) { - (p0,p1,p2)=(); - cross= (p1.x-p0.x)*(p2.y-p0.y) - (p2.x-p0.x)*(p1.y-p0.y); - } else { - throw UsageError,sprintf("%s: Wrong number of arguments.\n",_function_name()); - } - } - return sign(cross); -} diff --git a/src/slang/math/RD2rad.sl b/src/slang/math/RD2rad.sl index 909feb74..26e7bec1 100644 --- a/src/slang/math/RD2rad.sl +++ b/src/slang/math/RD2rad.sl @@ -6,7 +6,7 @@ define RD2rad() %\description % Given right ascension in hours, minutes, seconds and declination in degrees, % minutes of arc, seconds of arc, the function converts them to radians. Always -% provide six numbers, i.e., fill up with zeros where necessary. Right ascension +% provide six numbers, i.e., fill up with zeros were necessary. Right ascension % is assumed to be positive. If the degrees argument of declination is negative % the minus sign is automatically applied to the minutes and seconds parameters. % If the degrees argument of declination is zero, the sign of the minutes argument @@ -18,7 +18,7 @@ define RD2rad() % (raInRad, declInRad) = RD2rad(01, 45, 12.54, 0, 0, 45.34); % (raInRad, declInRad) = RD2rad([01, 01], [45,45.209], [12.54,0], [87,-87], [55,55], [45.34,45]); % (raInRad, declInRad) = RD2rad(rad2RD(0,-PI/4.)); -%\seealso{hms2deg,dms2deg,rad2RD} +%\seealso{rad2RD} %!%- { if(_NARGS != 6) diff --git a/src/slang/math/angular_separation.sl b/src/slang/math/angular_separation.sl deleted file mode 100644 index cdc74649..00000000 --- a/src/slang/math/angular_separation.sl +++ /dev/null @@ -1,46 +0,0 @@ -define angular_separation() -%!%+ -%\function{angular_separation} -%\usage{ sep=angular_separation(ra1,dec1,ra2,dec2);} -%\synopsis{calculates the angular distance between two points on a sphere} -%\qualifiers{ -%\qualifier{deg}{: if set, the input coordinates and output are in degrees (default: radian)} -%\qualifier{radian}{: if set, the input coordinates and output are in radian (the default))} -%} -%\description -% This routine calculates the angular separation between points on the sky -% with coordinates ra1/dec1 and ra2/dec2. -% -% The function is equivalent to greatcircle_distance but is compatible in terms -% of its qualifiers with the other coordinate system routines. It also returns -% the angular separation in the units of the input parameters rather than -% always in rad. -% -% This function is array safe (for either ra1/dec1 or ra2/dec2). -% -%\seealso{hms2deg,dms2deg,angle2string,greatcircle_distance} -%!%- -{ - variable ra1,dec1,ra2,dec2; - switch(_NARGS) - { case 4: (ra1,dec1,ra2,dec2)=(); } - { help(_function_name()); return; } - - if (qualifier_exists("deg")) { - ra1*=180./PI; - dec1*=180./PI; - ra2*=180./PI; - dec2*=180./PI; - } - - variable s,c,s1,c1,s2,c2; - (s,c)=sincos(ra2-ra1); - (s1,c1)=sincos(dec1); - (s2,c2)=sincos(dec2); - variable dist= atan2(sqrt((c2*s)^2 + (c1*s2 - s1*c2*c)^2), s1*s2 + c1*c2*c); - - if (qualifier_exists("deg")) { - return dist*180./PI; - } - return dist; -} diff --git a/src/slang/math/greatcircle_distance.sl b/src/slang/math/greatcircle_distance.sl index 42411e78..9f8c4705 100644 --- a/src/slang/math/greatcircle_distance.sl +++ b/src/slang/math/greatcircle_distance.sl @@ -51,4 +51,3 @@ define greatcircle_distance() (s2,c2)=sincos(y2); return atan2(sqrt((c2*sin(Dx))^2 + (c1*s2 - s1*c2*c)^2), s1*s2 + c1*c2*c); } - diff --git a/src/slang/strings/Roman.sl b/src/slang/strings/Roman.sl index f492a41f..45fb9231 100644 --- a/src/slang/strings/Roman.sl +++ b/src/slang/strings/Roman.sl @@ -1,62 +1,46 @@ %%%%%%%%%%%% -define Roman(n) { +define Roman(n) +%%%%%%%%%%%% %!%+ %\function{Roman} -%\synopsis{translates n to upper-case string with roman numeral} -%\usage{String_Type res = Roman(Integer_Type n)} -%\qualifiers{ -%\qualifier{latex}{: typeset minus sign ("$-$"R)} -%\qualifier{toobig [=""]}{: string that is returned of n is larger -% than the largest known Roman numeral (3999) } -%} -%\description -% Converts an integer into a uppercase roman numeral. -% Even though not known in Roman times, negative numbers are allowed. -% -% Algorithm based on -% https://www.geeksforgeeks.org/converting-decimal-number-lying-between-1-to-3999-to-roman-numerals/ -% +%\synopsis{replaces an integer in the with capital Roman numeral strings} +%\usage{String_Type rom = Roman(Integer_Type);} %\seealso{roman} %!%- - variable minus = "-"; - if(qualifier_exists("latex")) {minus = "$-$"R;} - variable oob = qualifier("toobig", ""); - variable nums,w1,w2; - variable array = 1; - - if (typeof(n)==Array_Type) { - return array_map(String_Type, &Roman, n); - } - - if (typeof(n)!=Integer_Type) { - throw UsageError,sprintf("%s: argument must be of Integer type",_function_name()); - } - - if (n==0) { - return "0"; - } - - variable num= [1,4,5,9,10,40,50,90,100,400,500,900,1000]; - variable sym = ["I","IV","V","IX","X","XL","L","XC","C","CD","D","CM","M"]; - - variable i=12; - variable res=""; - if (n<0) { - res=minus; - n=-n; - } - - while(n>0) { - variable div = n/num[i]; - n = n mod num[i]; - while(div>0) { - div--; - res+=sym[i]; - } - i--; - } - - return(res); +{ + switch(n) + { case 0: return "0"; } + { case 1: return "I"; } + { case 2: return "II"; } + { case 3: return "III"; } + { case 4: return "IV"; } + { case 5: return "V"; } + { case 6: return "VI"; } + { case 7: return "VII"; } + { case 8: return "VIII"; } + { case 9: return "IX"; } + { case 10: return "X"; } + { case 11: return "XI"; } + { case 12: return "XII"; } + { case 13: return "XIII"; } + { case 14: return "XIV"; } + { case 15: return "XV"; } + { case 16: return "XVI"; } + { case 17: return "XVII"; } + { case 18: return "XVIII"; } + { case 19: return "XIX"; } + { case 20: return "XX"; } + { case 21: return "XXI"; } + { case 22: return "XXII"; } + { case 23: return "XXIII"; } + { case 24: return "XXIV"; } + { case 25: return "XXV"; } + { case 26: return "XXVI"; } + { case 27: return "XXVII"; } + { case 28: return "XXVIII";} + { case 29: return "XXIX"; } + { case 30: return "XXX"; } + return "Don't know, please tell me."; } %%%%%%%%%%%% @@ -66,80 +50,9 @@ define roman(n) %\function{roman} %\synopsis{replaces an integer with lowercase Roman numeral strings} %\usage{String_Type rom = roman(Integer_Type);} -%\synopsis{translates n to lower-case string with roman numeral} -%\usage{String_Type roman = romann(Integer_Type n)} -%\qualifiers{ -%\qualifier{latex}{: typeset minus sign ("$-$"R)} -%\qualifier{toobig [=""]}{: string that is returned of n is larger -% than the largest known Roman numeral} -%} %\seealso{Roman} %!%- { - return strlow(Roman(n;; __qualifiers)); + return strlow(Roman(n)); } -define roman2int() -%!%+ -%\function{roman2int} -%\usage{Integer_Type=roman2int(String_Type)} -%\synopsis{translates a roman numeral into an integer} -%\description -% This function converts roman numerals into integers. -% The function is case insensitive and array safe. -%\seealso{roman,Roman} -%!%- -{ - variable r=(); - - if (typeof(r)==Array_Type) { - return array_map(Integer_Type,&roman2int,r); - } - - r=strup(r); - - % deal with negative numbers - variable sn=+1; - if (substr(r,1,1)=="-") { - sn=-1; - r=substr(r,2,strlen(r)); - } - - variable vals=Assoc_Type[Integer_Type]; - vals["0"]=0; - vals["I"]=1; - vals["V"]=5; - vals["X"]=10; - vals["L"]=50; - vals["C"]=100; - vals["D"]=500; - vals["M"]=1000; - - variable res=0; - variable i; - for (i=1; i<=strlen(r); i++) { - variable s1=vals[substr(r,i,1)]; - if (i=s2) { - res+=s1; - } else { - res+=s2-s1; - i++; - } - } else { - res+=s1; - } - } - - res=sn*res; - - % sanity check (otherwise we silently return erroneous results for - % illegal roman numerals) - if (r!=Roman(res)) { - throw UsageError,sprintf("%s: %s is not a valid roman numeral",_function_name(),r); - } - - - return sn*res; -} diff --git a/src/slang/strings/angle2string.sl b/src/slang/strings/angle2string.sl index 3420c201..828d1b27 100644 --- a/src/slang/strings/angle2string.sl +++ b/src/slang/strings/angle2string.sl @@ -35,7 +35,7 @@ define angle2string() % % This routine is array safe. % -%\seealso{hms2deg,dms2deg,generate_iauname} +%\seealso{hms2deg,dms2deg} %!%- { variable angle=(); diff --git a/src/slang/strings/generate_iauname.sl b/src/slang/strings/generate_iauname.sl deleted file mode 100644 index fe2e3da2..00000000 --- a/src/slang/strings/generate_iauname.sl +++ /dev/null @@ -1,51 +0,0 @@ -define generate_iauname() { -%!%+ -%\function{generate_iauname} -%\synopsis{for a given position, generate a coordinate string obeying the IAU convention} -%\usage{string=generate_iauname(ra,dec)} -%\qualifiers{% -% \qualifier{prefix}{: prefix for the string (e.g., "XMMU_")} -% \qualifier{radian}{: if set, the angles are in radians, not degrees} -%} -%\description -%This is a convenience routine to produce strings of the type -%XMMU Jhhmmss.s+ddmmss -%from J2000.0 positions. The routine obeys the IAU convention of -%truncating (not rounding!) the coordinate to the digits shown. -% -%The default-string is "Jhhmmss.s+ddmmss", use the prefix qualifier -%to prepend the mission name. -% -%Use angle2string to format coordinates with appropriate rounding. -% -%\seealso{dms2deg, hms2deg, angle2string, deg2dms} -%!%- - variable ra,dec; - - switch(_NARGS) - {case 2: (ra,dec)=(); } - {return help(_function_name()); - } - - variable hh,mm,ss; - variable dd,dm,ds; - - - if (qualifier_exists("radian")) { - ra*=180./PI; - dec*=180./PI; - } - - variable prefix=qualifier("prefix",""); - - variable sgn= (dec<0.) ? "-" : "+"; - - (hh,mm,ss)=deg2dms(ra;hours); - (dd,dm,ds)=deg2dms(abs(dec)); - - % chop off digits per IAU convention - ss=int(ss*10.)/10.; - - return sprintf("%sJ%02i%02i%04.1f%1s%02i%02i%02i", - prefix,int(hh),int(mm),ss,sgn,int(dd),int(dm),int(ds)); -} diff --git a/src/slang/strings/strsplit.sl b/src/slang/strings/strsplit.sl deleted file mode 100644 index e64d21b4..00000000 --- a/src/slang/strings/strsplit.sl +++ /dev/null @@ -1,36 +0,0 @@ -define strsplit() -%!%+ -%\function{strsplit} -%\synopsis{split a string at a separator and return contents as an array} -%\usage{Array_Type strsplit(String_Type s, String_Type separator);} -%\description -%This routine takes a string s and splits it into parts separated by -%a separator character. -%\example -%variable ra="12:34:56.78"; -%variable rastr,hh,mm,ss; -%rastr=strsplit(ra,":"); -%hh=atof(rastr[0]); mm=atof(rastr[1]); ss=atof(rastr[2]); -%ra=hms2deg(hh,mm,ss); -%!%- -{ - variable s,separator; - switch(_NARGS) - { case 2: (s,separator) = (); } - { return help(_function_name()); } - - if (length(separator)!=1) { - throw RunTimeError,sprintf("%s: Separator must be exactly one character\n",_function_name()); - } - variable retarr={}; - variable str=s; - variable pos=is_substr(str,separator); - while ( pos != 0 ) { - list_append(retarr,substr(str,1,pos-1)); - str=substr(str,pos+1,strlen(str)-pos); - pos=is_substr(str,separator); - } - list_append(retarr,str); - - return list_to_array(retarr); -} diff --git a/src/slang/structs/sort_struct_fields.sl b/src/slang/structs/sort_struct_fields.sl deleted file mode 100644 index 0e4d91ae..00000000 --- a/src/slang/structs/sort_struct_fields.sl +++ /dev/null @@ -1,29 +0,0 @@ -define sort_struct_fields() { -%!%+ -%\function{sort_struct_fields} -%\synopsis{sorts structure tags in a predefined order} -%\usage{ Struct_Type newstruc = empty_struct(Struct_Type s, String_Type taglist[]);} -%\description -% This function creates a new structure newstruc with the fiels given by the -% array taglist, and fills them with the values from structure s. -% -% The order in which the fields are added is given by taglist, fields not -% listed in taglist are ignored. -% -% The function is useful, e.g., when using fits_binary_table to create a FITS-file -% from a structure. In this case the order in which the FITS columns are created -% corresponds to the order in which they are contained in the structure. One can then -% use this function to create a FITS-file with the desired order of columns. -% -%!%- - variable s,fieldarr; - (s,fieldarr)=(); - - variable snew=@Struct_Type(fieldarr); - variable i; - _for i(0,length(fieldarr)-1,1) { - set_struct_field(snew,fieldarr[i],get_struct_field(s,fieldarr[i])); - } - - return snew; -} -- GitLab