Merge branch 'master' of https://github.com/teuben/nemo

inc/table.h

@@ -58,6 +58,7 @@ typedef struct {
 
   string *comments; // pointer to 'nh' (header/comment) lines
   string *lines;    // pointer to 'nr' lines (depends on mode)   lines[0], lines[1], ....
+  char   separ;     // optional enforced column separator
 
   size_t linelen;   // see Posix getline(3)
   char  *line;      // see Posix getline(3)
@@ -69,6 +70,7 @@ table  *table_open(stream instr, int mode);
 table  *table_open1(stream instr, int mode, int nlines);
 table  *table_cat(table *tp1, table *tp2, int mode);
 void    table_close(tableptr tptr);
+void    table_set(int mode, void *value);
 string  table_line(tableptr tptr);
 ssize_t table_line1(tableptr tptr, char **line, size_t *linelen, int newline);
 size_t  table_nrows(tableptr tprt);

man/man1/tabint.1

@@ -1,4 +1,4 @@
-.TH TABINT 1NEMO "6 April 2024"
+.TH TABINT 1NEMO "17 January 2025"
 
 .SH "NAME"
 tabint \- integrate a sorted table, optionally higher order moments
@@ -10,18 +10,21 @@ tabint \- integrate a sorted table, optionally higher order moments
 \fBtabint\fP integrates a function Y=F(X) that is defined by a set of (X,Y) points
 selected from two columns from an ASCII table.
 The function can either be resampled using a spline interpolation,
-or directly summed using the trapezoid rule.
+or directly summed using the trapezoid rule. X values need to sorted
+increasing or decreasing. 
 .PP
 See also \fItabhist(1NEMO)\fP with \fBcumul=t\fP for a visual display, though the
 cumulative option here can also produce a table for display to check the
 convergence of the integral.
 .PP
 An optional scale factor can be applied, which folds in any X or Y scaling that might
-want to be applied. This is particularly useful if you integrate the points themselves
+need to be applied. This is particularly useful if you integrate the points themselves
 (i.e. without giving a \fBstep=\fP)  but still want to scale the X coordinate.
 .PP
-Optionally the 1st or 2nd moment can be computed (step=, normalize= and cumulative=
+Optionally higher order moments (1 through 4) can be also computed (step=, normalize= and cumulative=
 are then ignored).
+.br
+0: integrated sum  1: mean X value  2: dispersion in X  3: skewness (0=symmetric)  4: kurtosis (0 means normal)
 
 .SH "PARAMETERS"
 .so man1/parameters
@@ -32,23 +35,23 @@ Input table file.
 No default.
 .TP
 \fBxcol=\fP
-Column with X coordinate. If the X coordinates are not increasing, the
+Column with sorted X coordinates. If the X coordinates are not increasing, the
 array will be reversed.
 .br
 [Default: 1]
 .TP
 \fBycol=\fP
-Column with Y coordinate of function
+Column with Y coordinate of the function
 .br
 [Default: 2]
 .TP
 \fBxmin=\fP
-if given, value below which data to be discarded
+If given, value below which data to be discarded
 .br
 By default it uses the smallest datapoint from the table.
 .TP
 \fBxmax=\fP
-if given, value above which data to be discarded
+If given, value above which data to be discarded
 .br
 By default it uses the largest datapoint from the table.
 .TP
@@ -78,7 +81,7 @@ if integrating the individual points.
 .br
 [Default: 1]
 .TP
-\fBmom=0|1|2\fP
+\fBmom=0|1|2|3|4\fP
 The moment along the profile. mom=0 is the flux, mom=1 the first moment along the X axis,
 and mom=2 the dispersion.
 .br
@@ -106,11 +109,12 @@ and intensity in Kelvin. We want to get an answer in K.km/s instead. We know the
    tabplot - 1 2 line=1,1 point=2,0.1 ycoord=0
 
 .EE
-recalling that \fBc\fP is in m/s and the requested line integral was requested in km/s.
+Recall that \fBc\fP is in m/s and the requested line integral was requested in km/s.
 
 .SH "CAVEATS"
 The table needs to be (reverse) sorted in X. A future option might sort, though
-\fIsort(1)\fP might already work. Here's an example where the 3rd column is an unsorted X coordinate.
+\fIsort(1)\fP might already work.
+Here's an example where the 3rd column is an unsorted X coordinate.
 The 6th column is the data value:
 .EX
 
@@ -120,6 +124,8 @@ perhaps easier
 
    % tabcols my.tab 3,6 | sort -n | tabint -
 .EE
+All moments, expecially the higher order ones, are susceptible to the correct window if the shape
+of a line feature is to be determined accurately.
 
 .SH "SEE ALSO"
 tabfilter(1NEMO), tabhist(1NEMO), tabmath(1NEMO), tabnllsqfit(1NEMO)¸ tabpeak(1NEMO), table(5NEMO)
@@ -138,4 +144,6 @@ Peter Teuben
 22-feb-2022	V0.6 added xmin=, xmax=, scale=	PJT
 31-jul-2023	V0.7 poor man's implementation of mom=	PJT
 6-apr-2024	V0.9 clarifications	PJT
+17-jan-2025	V1.0 added mom=3,4 options	PJT
+
 .fi

src/kernel/tab/tabint.c

@@ -7,6 +7,7 @@
  *      22-feb-23    add xmin,xmax,scale keywords
  *      31-jul-23    add mom=
  *       6-apr-24    cleanup
+ *      17-jan-25    add higher moments 3,4
  */
 
 #include <stdinc.h> 
@@ -26,13 +27,13 @@ string defv[] = {
   "normalize=f\n    Normalize integral",
   "cumulative=f\n   Show accumulation of integral",
   "scale=1\n        Scale factor to apply to integral",
-  "mom=0\n          0=flux 1=weighted mean 2=dispersion",
-  "VERSION=0.9\n    6-apr-2024 PJT",
+  "mom=0\n          0=flux 1=weighted mean 2=dispersion 3=skewness 4=kurtosis",
+  "VERSION=1.0\n    17-jan-2025 PJT",
   NULL,
 
 };
 
-string usage="integrate a (sorted) table";
+string usage="integrate or compute moments of a (sorted) table";
 
 
 extern int minmax(int, real *, real *, real *);
@@ -155,7 +156,8 @@ void nemo_main()
       }
       dprintf(1,"dx min/max: %g %g\n",dxmin,dxmax);
     } else {
-      double sum0 = 0.0, sum1 = 0.0, sum2 = 0.0, retval=0.0;
+      // this is for mom > 0
+      double sum0 = 0.0, sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4 = 0.0, retval=0.0;
 
       for (i=0; i<n; i++) {
 	sum0 +=  ydat[i];
@@ -164,12 +166,29 @@ void nemo_main()
       }
       sum1 /= sum0;
       sum2 /= sum0;
+      dprintf(1,"dx=%g\n", xdat[1]-xdat[0]);
+      dprintf(1,"sum0=%g\n",      sum0);
       dprintf(1,"sum1/sum0=%g\n", sum1);
       dprintf(1,"sum2/sum0=%g\n", sum2);
       sum2 = sqrt(sum2 - sum1*sum1);
       dprintf(1,"dispersion=%g\n", sum2);
       if (mom==1) retval=sum1;
       if (mom==2) retval=sum2;
+      if (mom>2) warning("new unchecked feature");
+      if (mom==3) {
+	// skewness is 0 for symmetric
+	sum3 = 0.0;
+	for (i=0; i<n; i++)
+	  sum3 += ydat[i]*qbe(xdat[i]-sum1);
+	retval = sum3/qbe(sum2)/sum0;
+      }
+      if (mom==4) {
+	// kurtosis is 0 for gaussian,   < 0 
+	sum4 = 0.0;
+	for (i=0; i<n; i++)
+	  sum4 += ydat[i]*sqr(sqr(xdat[i]-sum1));
+	retval = -3.0 + sum4/sqr(sqr(sum2))/sum0;
+      }
       printf("%g\n",retval);
       return;
     }

src/scripts/install_miniconda3

@@ -8,13 +8,14 @@
     url=https://repo.anaconda.com/miniconda/Miniconda3-py%s-%s-%s.sh
    url1=https://repo.anaconda.com/miniconda/Miniconda3-latest-%s-%s.sh
     # various miniconda versions available:
-   version=latest                 # !! special. watch out for caching this file !!
+version=latest                 # !! special. watch out for caching this file !!
 version=37_23.1.0-1            # 3.7.16
 version=38_23.11.0-2           # 3.8.18
 version=39_24.7.1-0            # 3.9.19
 version=310_24.7.1-0           # 3.10.15
 version=311_24.7.1-0           # 3.11.9
 version=312_24.7.1-0           # 3.12.4
+version=312_24.11.1-0          # 3.12.8
 
     dir=$(pwd)/miniconda3      # where miniconda will be located
    wget=wget                   # use wgetc is you have my cashing version (wget=curl is also allowed)