From a970dfa462e94c1014349bd30f5aae2d22ddf4e5 Mon Sep 17 00:00:00 2001
From: Ole Koenig <ole.koenig@fau.de>
Date: Tue, 10 Nov 2020 14:35:47 +0100
Subject: [PATCH] Add new function segment_lc_for_psd

Function to segment a lightcurve in preparation for a PSD, filtering
out the segments which do not fit into multiples of the segmentation
length.
---
 src/data/timing/segment_lc_for_psd.sl | 99 +++++++++++++++++++++++++++
 1 file changed, 99 insertions(+)
 create mode 100644 src/data/timing/segment_lc_for_psd.sl

diff --git a/src/data/timing/segment_lc_for_psd.sl b/src/data/timing/segment_lc_for_psd.sl
new file mode 100644
index 00000000..bc51462c
--- /dev/null
+++ b/src/data/timing/segment_lc_for_psd.sl
@@ -0,0 +1,99 @@
+define segment_lc_for_psd(lc,dt,dimseg) { %{{{
+%!%+
+%\function{segment_lc_for_psd}
+%\synopsis{Function to segment a lightcurve in preparation for a PSD,
+%    filtering out the segments which do not fit into multiples of
+%    the segmentation length.}
+%\usage{segment_lc_for_psd(lc,dt,dimseg)}
+%\qualifiers{
+%    gapfactor - relative factor, telling how large the gap is relative
+%                to dt [default: 1.1]
+%    ratefield - Name of rate field in lightcurve [default: rate]
+%    verbose   - Increase verbosity to display which parts of the
+%                lightcurve and what proportion was rejected.
+%}
+%\description
+%    Given a lightcurve with gaps, this function segments it such that
+%    only intervals of length "dimseg" (in bins) are present. Data
+%    that does not fit into integer multiples of the segmentation
+%    length are cut off.
+%    
+%    ATTENTION: Never use a discrete Fourier transform (e.g., 
+%    foucalc) with a segmentation length ("dimseg") larger than chosen
+%    in this function!
+%\notes
+%    The time array is not altered, so the lightcurve will still
+%    contain gaps, however none that have different length than
+%    "dimseg". In other words: The output data arrays can begin at
+%    arbitrary times, so are *not* sorted to always start at integer
+%    multiples of the segmentation length.
+%\example
+%    dt=1.; % (s)
+%    len=64.; % (s)
+%    offset=10; % background level
+%    T=5.; % (s), sinusoid periodidity
+%    omega=2*PI/T; % Angular frequency (rad/s)
+%    time_arr=[0:len:dt];
+%    lc = struct{ time = time_arr,
+%                 rate = offset+sin(omega*time_arr)+0.1*grand(int(len/dt)) };
+%    lc_gaps=struct{
+%             time=[lc.time[[0:20]],lc.time[[30:55]]],
+%             rate=[lc.rate[[0:20]],lc.rate[[30:55]]] };
+%    dimseg=16; % (bins)
+%    lc_split = segment_lc_for_psd(lc_gaps, dt, dimseg);
+%    res=foucalc(struct{time=lc_split.time,rate1=lc_split.rate},dimseg);
+%\seealso{split_lc_at_gaps, get_smallest_segment_in_lc, foucalc}
+%!%-
+
+  % define relative factor, telling how much the gap is relative to dt; default = 1.1;
+  variable gapfactor = qualifier("gapfactor",1.1); 
+
+  %allow qualifier that defines the name of the "rate" field; if not set, default is "rate"
+  variable ratefield = qualifier("ratefield","rate");
+    
+  % split lightcurve at gaps
+  variable split_lc = split_lc_at_gaps(lc,dt*gapfactor);
+
+  variable segmented_lc = struct{
+    time=Double_Type[0],
+    rate=Double_Type[0],
+    error=Double_Type[0]
+  };
+  
+  %% keep a record of the total time and the time thrown away
+  variable total_time=0.; % (s)
+  variable rejected_time=0.; 
+
+  variable ii;
+  _for ii (0,length(split_lc)-1,1) {
+    
+    total_time += length(split_lc[ii].time);
+    
+    % get how often the segment fits into a given part of the split lightcurve
+    variable nn = int(length(split_lc[ii].time)/dimseg);
+
+    % only do this if there is at least one segment
+    if (nn > 0) {
+      % note that each segment goes from 0 to nn*dimseg-1
+      % (the -1 is the important part!)	  
+      segmented_lc.time=[segmented_lc.time,split_lc[ii].time[[0:nn*dimseg-1]]];      
+      segmented_lc.rate=[segmented_lc.rate,  
+			 get_struct_field(split_lc[ii],ratefield)[[0:nn*dimseg-1]]];
+      if (struct_field_exists(split_lc[ii],"error")){
+	segmented_lc.error=[segmented_lc.error, split_lc[ii].error[[0:nn*dimseg-1]]];
+      }
+    } else {
+      rejected_time += length(split_lc[ii].time);
+      if (qualifier_exists("verbose")) {
+	vmessage("  ***Skipping lightcurve segment with index %i, length %.2fs",ii,length(split_lc[ii].time));
+      }
+    }
+  }
+
+  if (qualifier_exists("verbose")) {
+    vmessage("\nRejected time: %.2f s (%.3f percent of total data time=%.2fs)\n",
+	     rejected_time, rejected_time/total_time, total_time);
+  }
+  
+  return segmented_lc; 
+} %}}}
\ No newline at end of file
-- 
GitLab