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

docs/source/fits.rst

@@ -0,0 +1,10 @@
+FITS
+----
+
+
+FITS is one of the earliest data formats used in (observational) astronomy, and
+NEMO handles some of the conversions from and to FITS. 
+
+
+
+For the SDFITS extension there is limited support in `sdinfo` and `scanfits`.

docs/source/glossary.rst

@@ -38,7 +38,7 @@ Glossary
       "Flexible Interchange Transport System", a
       standard dataformat used to interchange data between
       machines. Commonly used for images, but has since its inception
-      be modified to store tables as well.
+      be modified to store tables as well. See also :term:`SDFITS`
 
   FWHM
       (Full Width Half Max): the effective resolution of the
@@ -73,7 +73,8 @@ Glossary
       or end with *orb*.
 
   pixel
-      PIXEL/voxel: an area in 2D or 3D representing
+      PIXEL/voxel: an area in 2D or 3D representing selected
+      phase space coordinates.
 
   potential
       Dataformat in NEMO used to represent a
@@ -89,7 +90,11 @@ Glossary
       A small user interface that pops up when a
       program is interrupted. Type ``quit`` to exit it, or ``?``
       for help. This feature of the user interface may not be
-      installed in your version. 
+      installed in your version.
+
+  SDFITS
+      A FITS extension that allows efficient storage of single dish
+      spectra. See also :term:`FITS`
 
   set
       Compound hierarchical data-structure of a
@@ -120,7 +125,10 @@ Glossary
   table
       A table consists of rows and columns of values, numbers or text.
       Most commonly stored in ASCII. Less well defined, it is one of the
-      four data types in NEMO.
+      four data types commonly used in NEMO.
+
+  voxel
+      See pixel
 
   yapp
       "Yet Another Plotting Package", the library

man/man1/ccdellint.1

@@ -36,7 +36,7 @@ For metric=2 these would be called true rings, for metric=1 a square.
 input data cube. No default.   
 .TP
 \fBradii=\fP
-radii of the ring boundaries (Nring+1). See \fBmetrix=\fP below for
+radii of the ring boundaries (Nring+1). See \fBmetric=\fP below for
 different ways to change the metric of computing a radius. No default.
 .TP
 \fBpa=\fP

man/man1/sdinfo.1

@@ -1,4 +1,4 @@
-.TH SDINFO 1NEMO "13 May 2023"
+.TH SDINFO 1NEMO "18 October 2023"
 
 .SH "NAME"
 sdinfo \- sdfits (BINTABLE) info and benchmark
@@ -79,6 +79,15 @@ Another benchmark option: allowing all data to appear in one long one dimensiona
 DATA[nrows*nchan] array vs. a two dimensional DATA[nrows][nchan] array.
 [1]
 
+.TP
+\fBtab=\fP
+If given, this will be the (ascii) table to which either 1dim or 2dim data (see datadim=) is written.
+For 1dim the x,y values printed are assumed Hz converted to GHz and raw intensity. For 2dim the row,x,y are printed,
+where row=1 is the first row. Apologies to non-radio users.
+Default: none.
+
+
+
 .SH "CAVEATS"
 By default NEMO compiles in double precision. The DATA1[nrows*nchan] "ONEDIM" mode will this occupy half the
 memory of a DATA2[nrows][nchan] array. NEMO will need to be compiled in single mode to have them use the same
@@ -139,4 +148,5 @@ Peter Teuben
 18-dec-2019	V0.7 Handle CIMAFITS and MATRIX extensions	PJT
 27-sep-2021	V0.9 Improved column reading, display TDIM	PJT
 15-mar-2023	V1.0 Finalized more benchmarks		PJT
+18-oct-2023	V1.1 add tab=	PJT
 .fi

man/man1/xrandom.1

@@ -113,14 +113,16 @@ Median               : 0.000425562
 .fi
 .PP
 xrandom can also be use to estimate the wall clock time of a complex pipe. Here
-is an example:
+is an example using the \fIcentiseconds since boot\fP timer from the \fBseed=-1\fP
+mode:
 .nf
+
   t0=$(xrandom -1);mkplummer - 1000000 | snapsort - - r | hackdens - .; t1=$(xrandom -1); nemoinp "($t1-$t0)/100"
   6.09
+
 .fi
-where 6.09 was measured on the fastest  machine I had available at the time of this writing
-(Jansky, with nemobench5 ~ 1350). Turns out for this example most of the time is spent in \fPhackdens(1NEMO)\fP,
-about 90%.
+where 6.09sec was measured on the fastest machine I had available at the time of this writing
+(Jansky, with nemobench5 ~ 1350). Turns out for this example about 90% of the CPU time is spent in \fPhackdens(1NEMO)\fP.
 
 .SH "GSL"
 If the GNU Scientific Library (GSL) has been enabled during installation, the

scripts/csh/mkmh97.sh

@@ -40,26 +40,26 @@
 #          22-jun-2023   added code=3 for princeton hackathon comparisons
 
 _script=mkmh97
-_version=22-jun-2023
+_version=6-oct-2023
 _pars=nemopars.rc
 _date=$(date +%Y-%m-%dT%H:%M:%S)
 
 #            text between #--HELP and #--HELP is displayed when --help is used
 #--HELP
-#            parameters for the integration
+#            parameters for the integration   [defaults are following Makino & Hut 1997]
 run=run0        # directory and basename of the files belonging to this simulation  [ empty means self] @todo
 nbody=2048      # number of bodies in one model
 m=1             # mass of second galaxy (mass of first will always be 1)
-em=0            # equal mass particles? (em=0 means nbody same for both galaxies)
+em=0            # equal mass particles? (em=0 means nbody same for both galaxies, thus individual masses not equal)
 fixed=0         # fixed potential for galaxy-1 ?
 step=1          # step in time when to dump snapshots
 v0=1.0          # initial impact/circular speed
 rp=0.0          # impact offset radius
 r0=16.0         # initial offset position for v0 > 0   [note v0 > 2/sqrt(r0)]
 eps=0.03125     # gravitational softening
 kmax=7          # integration timestep is 1/2**kmax
-eta=0.01        # alternative to timestep (only used by gravidy)
-code=1          # 0=hackcode1 1=gyrfalcON  2=bonsai2 (GPU)  3=rungravidy
+eta=0.01        # extra timestep accuracy (only used by gravidy)
+code=1          # 0=hackcode1  0q=hackcode1_qp  1=gyrfalcON   2=bonsai2 (GPU)  3=rungravidy  4=directcodeb
 seed=0          # random seed (use seed=123 for the benchmark)
 trim=0          # trim the data after analysis so only the last snapshot is kept
 save=1          # save the "final" plots in a subdirectory "movies" labeled with tstop
@@ -74,6 +74,7 @@ bsigma=0.0001                      # asinh/log breakover point
 tplot=0,5,10,15,20,25,30,40,50     # times to plot in evolution
 hackforce=hackforce                # which hackforce to use in the force re-calculation
 yapp=png                           # pick png, or vps (for yapp_pgplot) _ps for native ps
+progress=0                         # spawn a progress bar for the nbody run
 debug=1                            # 1=set -x,-e,-u   0=nothing
 #
 #--HELP
@@ -172,10 +173,10 @@ if [ $restart = 1 ]; then
 	snapstack $run.1 $run.2 $run.3 deltar=$r0,0,0   deltav=0,$v0,0   zerocm=t
     fi
 
-    # integrator:  0:  hackcode1 is O(NlogN) code
+    # integrator:  0:  hackcode1 is O(NlogN) code 
     #              1:  gyrfalcON is O(N)
     #              2:  bonsai2 is O(N) but scales faster for "small" N
-    #              3:  rungravidy
+    #              3:  rungravidy O(N^2) code
     echo "Use:   tail -f $run/$run.4.log     to monitor progress of the integrator"
     if [ $fixed = 1 ]; then
 	echo Fixed potential now, only code=0 is supported
@@ -193,7 +194,11 @@ if [ $restart = 1 ]; then
 	if [ $code = 0 ]; then
 	    hackcode1 $run.3 $run.4 eps=$eps freq=2**$kmax freqout=1/$step fcells=2 tstop=$tstop options=mass,phase,phi,acc > $run.4.log    
 	    snapdiagplot $run.4 tab=$run.4.etot
+	elif [ $code = 0q ]; then
+	    hackcode1_qp $run.3 $run.4 eps=$eps freq=2**$kmax freqout=1/$step fcells=2 tstop=$tstop options=mass,phase,phi,acc > $run.4.log    
+	    snapdiagplot $run.4 tab=$run.4.etot
 	elif [ $code = 1 ]; then
+	    progress_bar.sh $progress $run.4.log 1 $tstop &
 	    gyrfalcON $run.3 $run.4 eps=$eps kmax=$kmax step=$step tstop=$tstop give=mxvap > $run.4.log
 	    tabcols $run.4.log 1,2 > $run.4.etot
 	elif [ $code = 2 ]; then
@@ -215,16 +220,22 @@ if [ $restart = 1 ]; then
 	    fi	    
 	    echo Wrote final combined snapshot in $run.4
 	elif [ $code = 3 ]; then
+	    progress_bar.sh $progress $run.4.d/$run.4.d.out.log 2 $tstop &	    
 	    /usr/bin/time rungravidy $run.3 $run.4.d tcrit=$tstop deltat=$step eps=$eps eta=$eta > $run.4.log
 	    ln -s $run.4.d/OUT3.snap $run.4
 	    tabcols $run.4.d/$run.4.d.out.log 2,5 > $run.4.etot
+	    # @todo ### Fatal error [snapcenter]: total weight is zero
+	elif [ $code = 4 ]; then
+	    directcode $run.3 $run.4 eps=$eps freq=2**$kmax freqout=1/$step tstop=$tstop options=mass,phase,phi,acc > $run.4.log    
+	    snapdiagplot $run.4 tab=$run.4.etot
 	else
 	    set +x
-	    echo Unknown code=$code, valid are:
-	    echo 0 = hackcode1
-	    echo 1 = gyrfalcON
-	    echo 2 = bonsai2
-	    echo 3 = rungravidy
+	    echo "Unknown code=$code, valid are:"
+	    echo "0 = hackcode1"
+	    echo "1 = gyrfalcON"
+	    echo "2 = bonsai2"
+	    echo "3 = rungravidy"
+	    echo "4 = directcode"
 	fi
     fi
 else
@@ -248,7 +259,7 @@ if [ ! -e $run.xv.tab ]; then
 
     #
     echo "# x1 y1 z1 vx1 vy1 vz1"                                                                                   >$run.4.g1.tab
-    if [ $code == 2 ]; then
+    if [ $code == 2 ] || [ $code == 3 ]; then
 	snapcopy $run.4 - "select=i<$nbody?1:0" | $hackforce - - debug=-1 | snapcenter - . "-phi*phi*phi" report=t >>$run.4.g1.tab
     else
 	snapcenter $run.4 . "weight=i<$nbody?-phi*phi*phi:0" report=t                                              >>$run.4.g1.tab

src/image/fits/sdinfo.c

@@ -41,7 +41,8 @@ string defv[] = {
     "bench=1\n           How many times to run benchmark",
     "mode=-1\n           Mode how much to process the SDFITS files (-1 means all)",
     "datadim=2\n         1: DATA[nrows*nchan]   2: DATA[nrows][nchan]",
-    "VERSION=0.9h\n      2-may-2023 PJT",
+    "tab=\n              If given, produce tabular output (datadim=1 or 2 only)",
+    "VERSION=1.1\n       18-oct-2023 PJT",
     NULL,
 };
 
@@ -125,6 +126,8 @@ void average2(int ndata, real *data1, real *data2, real *ave)
 
 /*
  *   perform a baseline subtraction
+ *
+ *   @todo the X coordinate is not w.r.t. some center (or multiple centra for each order polynomial)
  */
 
 
@@ -267,8 +270,9 @@ void nemo_main(void)
     fitsfile *fptr;       /* pointer to the FITS file; defined in fitsio.h */
     int status, fmode, ii, jj, i, j, k, data_col;
     long int nrows;
-    long int nrows0;
-    int ncols, nchan, nfiles, found, row, nhdus;
+    // long int nrows0;
+    int ncols, nfiles, found, row, nhdus;
+    int nchan = -1;
     string fname = getparam("in"), *fnames;
     string *colnames;
     int nsize, ndims, dims[MDMAXDIM];
@@ -362,6 +366,8 @@ void nemo_main(void)
 	}
 	dprintf(1,"%3d:  %-16s %s\n",ii,colname,data_fmt);
       } //for(i)
+      if (data_col < 0) warning("No DATA column found");
+      if (nchan < 0) warning("DATA with no nchan???");      
       real np = nrows * nchan * 4 / 1e6;
       dprintf(0,"%s : Nrows: %d   Ncols: %d  Nchan: %d  nP: %g Mp\n",fname,nrows,ncols,nchan,np);
       int col_tcal  = keyindex(ncols, colnames, "TCAL") + 1;
@@ -373,7 +379,6 @@ void nemo_main(void)
       dprintf(0,"tcal: %d   cal: %d sig: %d fdnum: %d ifnum: %d plnum: %d\n",
 	      col_tcal, col_cal, col_sig, col_fdnum, col_ifnum, col_plnum);
 
-
       // get optional columns
       int *fdnum_data = get_column_int(fptr, "FDNUM", nrows, ncols, colnames);
       int *ifnum_data = get_column_int(fptr, "IFNUM", nrows, ncols, colnames);
@@ -413,6 +418,17 @@ void nemo_main(void)
       printf("SIG:   %d  = %c\n", nsig, nsig==1 ? sig_data[0][0] : ' ');
       printf("CAL:   %d  = %c\n", ncal, ncal==1 ? cal_data[0][0] : ' ');
 
+
+      // the 1st coordinate for SDFITS is important
+      int col_crval = keyindex(ncols, colnames, "CRVAL1") + 1;
+      int col_crpix = keyindex(ncols, colnames, "CRPIX1") + 1;
+      int col_cdelt = keyindex(ncols, colnames, "CDELT1") + 1;
+      dprintf(0,"wcs1: %d %d %d \n", col_crval, col_crpix, col_cdelt);
+
+      double *crval_data = get_column_dbl(fptr, "CRVAL1", nrows, ncols, colnames);
+      double *crpix_data = get_column_dbl(fptr, "CRPIX1", nrows, ncols, colnames);
+      double *cdelt_data = get_column_dbl(fptr, "CDELT1", nrows, ncols, colnames);
+
       if (mode >= 0 && mode < 2) continue;            
       dprintf(1,"MODE=2\n");
 
@@ -477,6 +493,16 @@ void nemo_main(void)
 	  else
 	    dprintf(0,"DATA1 %g %g ... %g (only 1 row)\n",data1[0],data1[1],data1[nchan-1]);
 	  if (blfit >= 0) warning("ONEDIM mode cannot subtract baselines");  // also due to float/real
+	  if (hasvalue("tab")) {
+	    double xval;
+	    stream tstr = stropen(getparam("tab"),"w");	    
+	    for(int i=0; i<nchan; i++) {
+	      xval = (i+1-crpix_data[0])*cdelt_data[0] + crval_data[0];
+	      xval /= 1e9;
+	      fprintf(tstr,"%g %g\n",xval,data1[i]);
+	    }
+	    strclose(tstr);	    
+	  }
 	} else {
 	  // 2dim array, suitable for waterfall plot. Note this could be double precision in NEMO
 	  dprintf(0,"TWODIM DATA2: make Waterfall\n");
@@ -492,6 +518,18 @@ void nemo_main(void)
 	    for (ii=0; ii<nrows; ii++)
 	      baseline(nchan, NULL, data2[ii], blfit, coeffs, errors);
 	  }
+	  if (hasvalue("tab")) {
+	    double xval;
+	    stream tstr = stropen(getparam("tab"),"w");
+	    for (j=0; j<nrows; j++) {
+	      for (i=0; i<nchan; i++) {
+		xval = (i+1-crpix_data[j])*cdelt_data[j] + crval_data[j];
+		xval /= 1e9;
+		fprintf(tstr,"%d %g %g\n", j+1, xval, data2[j][i]);
+	      }
+	    }
+	    strclose(tstr);
+	  }
 	}
 
 	if (mode >= 0 && mode < 5) continue;

src/scripts/Makefile

@@ -5,7 +5,8 @@ include $(NEMOLIB)/makedefs
 
 SCRIPTS = mknemo mkpdoc manlaser nemoman nemobench ldso nemo.version \
 	fitsedit fitsmerge nds9 shtool nemo.coverage showvar mir2nemo getline pedit \
-	wgetc yapp_query nemoplot nemopars sbatch_nemo.sh nemo_functions.sh
+	wgetc yapp_query nemoplot nemopars \
+        sbatch_nemo.sh nemo_functions.sh progress_bar.sh
 MAKES = Makefile.lib 
 
 help:

src/scripts/install_anaconda3

@@ -16,6 +16,7 @@ version=2022.10    # 3.9.13 sphinx
 version=2023.03-0  # 3.10.9
 version=2023.03-1  # 3.10.9
 version=2023.07-2  # 3.11.4
+version=2023.09-0  # 3.11.5
 
 
    dir=$(pwd)/anaconda3       # where anaconda will be located

src/scripts/progress_bar.sh

@@ -0,0 +1,49 @@
+#! /usr/bin/env  bash
+#
+
+# usage:   progress table column max
+#
+#  progress=0     none
+#           1     zenity
+#           2     qt5
+
+pro=$1
+tab=$2
+col=$3
+max=$4
+
+
+echo "progress.sh: $pro $tab $col $max"
+if [ $pro = 0 ]; then
+    exit 0
+fi
+
+#  some programs take some startup time to create the progress file
+for i in $(seq 5); do
+    if [ ! -e $tab ]; then
+	echo SLEEPING for $tab
+	sleep 1
+    else
+	echo Finally saw $tab
+	break
+    fi
+done
+
+#   if progress table still doesn't exist, give up
+if [ ! -e $tab ]; then
+    echo "$tab does not exist, giving up"
+    exit 0
+fi
+
+#    another minor sleep
+sleep 2
+
+#    looping to check progress, passing on the percentage completeness to zenity
+(while [ 1 = 1 ]; do
+ now=$(tail -1 $tab | tabcols - $col);
+ p=$(nemoinp "100*${now}/${max}" format=%d);
+ printf "%d\n"  $p;
+ sleep 1;
+done)   | zenity --progress --auto-close --time-remaining --title "NEMO: $tab" --text "$tab" --percentage=0 
+
+