diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
index d78b3725d..5cf909658 100644
--- a/.github/workflows/ccpp.yml
+++ b/.github/workflows/ccpp.yml
@@ -37,7 +37,7 @@ jobs:
build2:
- runs-on: macos-12
+ runs-on: macos-13
steps:
- uses: actions/checkout@v1
diff --git a/VERSION b/VERSION
index cbe06cdbf..a84947d6f 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-4.4.4
+4.5.0
diff --git a/docs/examples/spectrum-fit.html b/docs/examples/spectrum-fit.html
index ac2e3a307..741aee28b 100755
--- a/docs/examples/spectrum-fit.html
+++ b/docs/examples/spectrum-fit.html
@@ -13,6 +13,11 @@ Spectrum Fit
tabplot spectrum1fit.tab 1 3 -50 50 -1 5 line=1,1 ycoord=0 yapp=2/xs
+ raw spectrum
+ 
+ baseline subtracted spectrum
+ 
+
This page was last modified on
diff --git a/docs/install_nemo.sh b/docs/install_nemo.sh
index 071fe3178..a48a444ba 100755
--- a/docs/install_nemo.sh
+++ b/docs/install_nemo.sh
@@ -18,7 +18,7 @@
# opt=1 python=1 14'00"
#
-echo "install_nemo.sh: Version 1.9 -- 19-feb-2024"
+echo "install_nemo.sh: Version 1.10 -- 20-oct-2024"
opt=0 # install some optional mknemos= packages
nemo=nemo # root directory where NEMO will be installed (. = here)
@@ -141,7 +141,7 @@ fi
# pick a configure
-CC=gcc$v CXX=g++$v FC=gfortran$v ./configure $opt $with_yapp
+CC=gcc$v CXX=g++$v F77=gfortran$v ./configure $opt $with_yapp
#./configure $opt --enable-debug --with-yapp=pgplot --with-pgplot-prefix=/usr/lib # ok
#./configure $opt --enable-debug --with-yapp=pgplot --with-pgplot-prefix=/usr/lib --enable-pedantic
#./configure $opt --enable-debug --with-yapp=pgplot
diff --git a/docs/requirements.txt b/docs/requirements.txt
index e69e1df3c..b7ba3a4c8 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -2,3 +2,6 @@ numpydoc>=0.9.1
sphinx>=4.1.1
sphinx-rtd-theme
sphinx_inline_tabs
+# dysh
+sphinx-rtd-theme
+sphinxcontrib-mermaid
diff --git a/docs/source/codes.rst b/docs/source/codes.rst
index 872200eae..fdc8b14b7 100644
--- a/docs/source/codes.rst
+++ b/docs/source/codes.rst
@@ -20,6 +20,8 @@ that have a tighter connection to NEMO:
.. include:: starlab.rst
+.. include:: falcon2.rst
+
.. include:: rebound.rst
.. include:: galpy.rst
diff --git a/docs/source/falcon2.rst b/docs/source/falcon2.rst
new file mode 100644
index 000000000..db21b8ebb
--- /dev/null
+++ b/docs/source/falcon2.rst
@@ -0,0 +1,150 @@
+FALCON2
+-------
+
+
+The falcon V2 package is still under development for a public release. The main
+integrator is now called ``griffin``, replacing the older ``gyrfalcON`` from falcon V1.
+Although falcon V1 is distributed with NEMO, it is not maintained and may eventually
+succumb to software rot. Users should switch to using ``griffin``.
+
+A paper describing the ``griffin`` FMM code is in Dehnen (2014)
+
+.. warning::
+ Adding FALCON2 to NEMO will result in some programs that have duplicated names, e.g. **mkplummer**.
+
+
+Installation
+~~~~~~~~~~~~
+
+Hopefully soon we can use
+
+.. code-block::
+
+ mknemo falcON2
+
+but unless you are a close collaborator, this will likely not be working.
+
+This type of install would place the package in ``$NEMO/local/falcon2``, the
+file ``INSTALL.md`` describes the installation procedure.
+
+
+
+Command Line Interface
+~~~~~~~~~~~~~~~~~~~~~~
+
+Although the CLI will look familiar to NEMO users, there are some salient differences.
+Some highlights:
+
+- the CLI is a set of ``key=val``, like in NEMO. If given in order, the ``key=`` portion can be
+ omitted, leaving just a series of values. This not recommended in scripts.
+- ``--help`` describes the keywords and some help, much like the ``help=h`` option in NEMO
+- ``help=1`` shows hidden CLI options (NEMO doesn't have hidden options)
+- ``debug=``
+
+
+
+
+Examples
+~~~~~~~~
+
+
+For details on a specific program, type
+
+.. code-block::
+
+ program --help
+
+
+- Create a Plummer sphere with 100 particlces
+
+.. code-block::
+
+ mkplummer p100.in 100
+
+- View the contents of this HDF5 dataset
+
+.. code-block::
+
+ dump p100.in
+
+- Convert falcon2 HDF5 files to NEMO snapshot, and review like the ``dump`` program
+
+.. code-block::
+
+ s2a p100.in | tabcomment - - delete=t | tabtos - p100.bsf block1=m,pos,vel,skip nbody=100
+ tsf p100.bsf
+
+- Comparing performance of griffin in parallel (oneTBB) mode:
+
+.. code-block::
+
+ rm -f p10k.*
+ mkplummer p10k.in 10000 time=0
+ /usr/bin/time griffin p10k.in p10k.out step=1 tstop=10 eps=0.05 tau=2^-4
+ threads=0 98.41user 19.02system 0:24.23elapsed 484%CPU
+ threads=1 21.19user 0.79system 0:22.03elapsed 99%CPU
+ threads=2 24.06user 2.90system 0:15.05elapsed 179%CPU
+ threads=4 33.37user 4.73system 0:13.77elapsed 276%CPU 40.70user 5.90system 0:16.47elapsed 282%CPU
+ threads=8 62.97user 11.25system 0:21.23elapsed 349%CPU
+ 12 82.25user 14.55system 0:24.61elapsed 393%CPU
+ 16 93.44user 16.78system 0:25.63elapsed 430%CPU
+ 20 97.90user 18.67system 0:23.95elapsed 486%CPU
+
+
+- Comparing falcon1 with falcon2
+
+.. code-block::
+
+ # some parameters
+ nbody=10
+ eps=0.05
+ kmax=8
+ tstop=1
+ #
+ rm -f p1.*
+ mkplummer p1.in $nbody time=0
+ griffin p1.in p1.out step=1 tstop=$tstop eps=$eps tau=2^-$kmax
+ # 8.7sec
+
+ s2a p1.in | tabcomment - - delete=t | tabtos - p1.bsf block1=m,pos,vel,skip nbody=$nbody
+ gyrfalcON p1.bsf p1.out2 step=1 tstop=$tstop eps=$eps kmax=$kmax
+ # 0.6 sec
+
+ # comparing initials (notice p1.out does not contain times=0)
+ echo "=== Initial conditions ==="
+ s2a p1.in times=0 | tabcols - 2:7
+ snapprint p1.out2 times=0 format=%11.9e
+
+
+ # comparing final
+ echo "=== Final snapshot ==="
+ s2a p1.out times=$tstop | tabcols - 2:7
+ snapprint p1.out2 times=$tstop format=%11.9e
+
+
+but griffin is more accurate and such a direct comparison may not be fair. In particular, is tau=2^-$kmax even fair?
+
+- Comparing X and Y
+
+
+Surprises
+~~~~~~~~~
+
+to be resolved
+
+- dump --help
+
+ does not show the help we usually see, but seems to think --help is a file
+
+
+- "griffin --help"
+
+ says: "please provide 'out', 'tau', 'eps'"
+ why complain, we didn't attempt to run.
+
+- mkplummer has a keyword 'q-ran', but itsn't it better to use q_ran, since that's
+ a more common one used in all falcon programs. keep it consistent.
+
+- mkplummer has a default time=1
+
+- ..
diff --git a/docs/source/iface.rst b/docs/source/iface.rst
index 65fae831e..8d54da895 100644
--- a/docs/source/iface.rst
+++ b/docs/source/iface.rst
@@ -163,13 +163,16 @@ results in something like:
fill_circle=t frame= VERSION=1.3f
-As you see, ``snapplot`` happens to be a program
+Compare the VERSION= with your version, and it will likely be different.
+
+
+As you can see, ``snapplot`` happens to be a program
with quite an extensive parameter list.
Also note that ``help`` itself is not listed in the above list of program
keywords because it is a **system keyword**
(more on these later).
-There are a few *short-cut*
+There are a few *short-cuts*
in this user interface worth mentioning
at this
stage. First of all, keywords don't have to be specified
@@ -178,7 +181,7 @@ order, they will be associated by the appropriate keyword.
The order of program keywords can be seen with
the keyword ``help=``.
The moment you deviate from
-this order, or leave gaps, all
+this order, or skip keywords, all
values must be accompanied by their keywords, *i.e.* in
the example
@@ -301,7 +304,7 @@ In summary, the system keywords are:
Deprecated
- **np=**
- Number of processors (for OpenMP) to maximally use. Default is max.
+ Number of processors (e.g. for OpenMP) to maximally use. Default is the max.
For a more detailed description of the system keywords and all their options
see :ref:`aiface`. The actual degree of implementation of the system
@@ -326,15 +329,8 @@ shells like ``tcsh`` and ``bash``
can be used very efficiently in this mode.
In batch mode shell scripts, if used properly, can provide a very
powerful method of running complex simulations.
-Other plug-compatible
-interfaces that are available are ``mirtool`` and ``miriad``,
-described in more detail in
-Appendix~\ref{s:mirtool} and \ref{s:miriad} There was also a
-Khoros (cantata, under khoros V1)
-interface (``http://www.khoral.com``) available, but this product is not
-open source anymore.
Lastly, lets not forget scripting languages like python, perl and ruby.
-Although the class UNIX (c)sh shell is very WYSIWYG, with a modest amount
+Although the classic UNIX (c)sh shell is very WYSIWYG, with a modest amount
of investment the programmability of higher level scripts can give you
a very powerful programming environment.
@@ -391,10 +387,14 @@ caveat, here are various help options:
keyword descriptions and more vertical space.
The special ``--help`` option is allowed for those with gnu fingers.
+ Using ``-h`` is a short-cut.
The special ``--man`` option delivers the unix style man
page (see next item).
+ The special ``--version`` option displays the current version.
+ Using ``-v`` is a short-cut.
+
- Unix manual pages
for programs, functions, and file
formats, all in good old UNIX tradition. All these
@@ -437,7 +437,8 @@ caveat, here are various help options:
as well as advanced users, and at is being coverted to this RST
manual, outdated.
-- This manual, in **reStructuredText** might be available in many different formats. html and pdf are the common ones.
+- The manual you are reading here,
+ in **reStructuredText** might be available in many different formats. html and pdf are the common ones.
.. _aiface:
diff --git a/docs/source/progr.rst b/docs/source/progr.rst
index f0f464a8d..ecded0afa 100644
--- a/docs/source/progr.rst
+++ b/docs/source/progr.rst
@@ -24,17 +24,17 @@ We start by looking at the *Hello Nemo* program:
.. code-block:: C
- #include /* standard (NEMO) definitions */
+ #include /* 1. standard (NEMO) definitions */
- string defv[] = { /* standard keywords and default values and help */
+ string defv[] = { /* 2. standard keywords and default values and help */
"n=10\n Number of interations", /* key1 */
"VERSION=1.2\n 25-may-1992 PJT", /* key2 */
NULL, /* standard terminator of defv[] vector */
};
- string usage = "Example NEMO program 'hello'"; /* usage help text */
+ string usage = "Example NEMO program 'hello'"; /* 3. usage help text */
- void nemo_main() /* standard start of any NEMO program */
+ void nemo_main() /* 4. standard start of any NEMO program */
{
int n = getiparam("n"); /* get n */
@@ -42,6 +42,18 @@ We start by looking at the *Hello Nemo* program:
if (n < 0) /* deal with fatal */
error("n=%d is now allowed, need >0",n); /* errors */
}
+
+
+Although this looks like C, there are four odd looking ways a typical NEMO program
+is structured.
+
+1. An include file that captures NEMO definitions (such as string)
+2. A global string array defv/[] that defines the command line user interface
+3. A global string that captures a usage
+4. A global function nemo_main() that defines the starting point of the NEMO program
+
+The standard C main() is present in the NEMO library, with which each program
+is linked, but ensures the program is set up for the NEMO environment.
@@ -50,20 +62,25 @@ The NEMO Programming Environment
The modifications necessary to your UNIX environment in order to
access NEMO are extensively described elsewhere. This not only
-applies to a user, but also to the application programmer.
+applies to a user, but also to developers.
In summary, the essential changes to your environment consist of
one simple additions to your local shell startup file
-or you can do it interactively in your shell (be it sh or csh like).
+or you can do it interactively in your shell, viz.
.. code-block:: bash
% source /opt/nemo/nemo_start.sh
- or
- % source /opt/nemo/nemo_start.csh
where the location of NEMO=/opt/nemo is something that will likely
-be different for you.
+be different for you. This example is for bash-like shells, there
+is a similar file for csh, and there is also an option so use modules,
+such that a command like this will work:
+
+.. code-block:: bash
+
+ % module load nemo
+
The NEMO Macro Packages
diff --git a/docs/source/scripting.rst b/docs/source/scripting.rst
index e12ad3737..ed727107e 100644
--- a/docs/source/scripting.rst
+++ b/docs/source/scripting.rst
@@ -1,14 +1,11 @@
Scripting
---------
-
-.. todo:: more scripting examples needed here
-
-We will focus on **bash** scripting here, basing comments
-on the **mkmh97.sh** script.
+We will focus on **bash** scripting here, for an example see the
+**mkmh97.sh** script.
A typical script will contain the following components
-1. *hash-bang* the first line of the script, and optionally include nemo_functions.sh
+1. *hash-bang* the first line of the script, and optionally include ``nemo_functions.sh``
where various convenience functions are maintained
@@ -32,8 +29,9 @@ and make the script executable
run=run1 # directory and run ID #> ENTRY
nbody=1000 # number of bodies #> RADIO 1000,10000,100000,1000000
+ eps=0.05 # softening #> SCALE 0:1:0.01
-and the parsing could be achieved as follows
+and then parsing could be achieved as follows
.. code-block:: bash
@@ -42,7 +40,7 @@ and the parsing could be achieved as follows
done
3. optionally, but highly recommended, some matching **--HELP** tags that can be retrieved using
-a --help command line option
+the --help command line option
.. code-block:: bash
@@ -77,3 +75,49 @@ a --help command line option
5. the body of the script, shell variables, nemopars variables etc. should be used
where possible.
+
+
+
+Example: Extracting results from run directories
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+A common workflow is to run a series of simulations and walking over a multi-dimensional parameter space.
+This usually results into running each simulations in its own run-directory, plus storing
+parameters and possibly resulting values in a parameter file for later extraction. The directory name
+could even contain the value of these parameters as well. Several common strategies have been seen in
+the wild, we list three:
+
+.. code-block::
+
+ # 1. simply enumerated (e.g id=001,002,003,....) parameters stored inside
+ run_${id}
+
+ # 2. linear list of directories, parameters stored inside
+ run_${a}_${b}_${c}
+
+ # 3. hierarchical directories, parameters stored inside
+ run/$a/$b/$c
+
+Although easier to visually identify the values of the parameters in 2. and 3., they don't scale very
+well if a new parameter is introduced. In the first case a simple lookup table can be created using
+``nemopars``, thus making it easier to find which parameters are used in while run directory. Here's
+an example:
+
+.. code-block::
+
+ nemopars id,a,b,c run_*/nemopars.rc > run.pars
+
+
+
+Summary
+~~~~~~~
+
+Summarizing, here are the recommended methods to maintain and extract NEMO variables.
+
+1. **nemopars**: extract parameters from a bash-style rc file (python should also be able to use it)
+
+2. **nemovar**: get and set NEMOVAR variables
+
+3. **show_vars**: alias via nemo_functions.sh
+
+
diff --git a/docs/whatsnew.html b/docs/whatsnew.html
index decf3c9ea..4714e666f 100644
--- a/docs/whatsnew.html
+++ b/docs/whatsnew.html
@@ -21,13 +21,23 @@ What's New in NEMO
-->
+
+ 4.5.0: in git
+
+
- 4.4.3/4.4.4: in git
+ 4.4.3/4.4.4/4.4.5 (a.k.a. Malta Release)
+ - (Dec 22, 2024) define TIMEFUZZ in stdinc.h, not per program
+
- (Dec 12, 2024) drafting two codes: diskspectrum and lineid
+
- (Dec 7, 2024) added falcon2 docs.
+
- (Oct 20, 2024) more code-rot fixed for gcc-14
- (Jun 21, 2024) enhanced tabsmooth with more filters
- (Apr 20, 2024) experimental falcon program
to process falcON style (HDF5) snapshots for falcON2. "mknemo falcon2" has
- been created to get the source code for falcON2.
+ been created to get the source code for falcON2 and hopefully guide the install
- (Apr 1, 2024) experimental nemovar version
to communicate variables, phase 1.
- (Mar 19, 2024) the nemobench script has an extra first argument, a label. this will
diff --git a/inc/stdinc.h b/inc/stdinc.h
index 467c5c2da..54228af67 100644
--- a/inc/stdinc.h
+++ b/inc/stdinc.h
@@ -330,6 +330,13 @@ typedef real (*float_proc)(float);
#define DR2AS 206264.8062470963551564733
#define DAS2R 4.848136811095359935899141e-6
+/*
+ * TIMEFUZZ
+ * fuzzy time for snapshots; ideally really should not be hardcoded
+ *
+ */
+#define TIMEFUZZ 0.000001
+
/*
* POS_ANGLE, SYM_ANGLE: bring angular variables into standard form.
* To map an angular variable 'phi' into the range [0, 2pi),
diff --git a/inc/version.h b/inc/version.h
index 1bf89aab6..eb2faa7a0 100644
--- a/inc/version.h
+++ b/inc/version.h
@@ -1,5 +1,5 @@
/* do not edit - created by the src/scripts/version script */
#if defined(GETPARAM_VERSION_ID)
-static char *version_h = "git_11690"; // git rev-list --count HEAD
+static char *version_h = "git_12231"; // git rev-list --count HEAD
#endif
-#define NEMO_VERSION "4.4.3"
+#define NEMO_VERSION "4.5.0"
diff --git a/man/man1/diskspectrum.1 b/man/man1/diskspectrum.1
new file mode 100644
index 000000000..524464282
--- /dev/null
+++ b/man/man1/diskspectrum.1
@@ -0,0 +1,131 @@
+.TH DISKSPECTRUM 1NEMO "28 November 2024"
+
+.SH "NAME"
+diskspectrum \- global spectrum of a rotating thin disk
+
+.SH "SYNOPSIS"
+\fBdiskspectrum out=\fPsnapshot [parameters=values ...]
+
+.SH "DESCRIPTION"
+\fIdiskspectrum\fP sets up a ``cold'' disk of test particles
+in the gravitational potential of a user-supplied potential in
+\fIpotential(5NEMO)\fP format. The disk is imagined to be inclined by
+a given inclination and a global spectral profile is computed.
+.PP
+This program was cloned from \fImkdisk(1NEMO)\fP, with which it shares
+a number of keywords. Instead of producing a \fIsnapshot(5NEMO)\fP, it
+produces a table, containing the global spectrum.
+.PP
+Warning: this program is under development, for example, currently it
+produces the unsmoothed half spectrum shifted to 0.
+
+.SH "PARAMETERS"
+The following parameters are recognized in any order if the keyword is also
+given:
+.TP 25
+\fBout=\fIout_file\fP
+The data are written to this file in \fIsnapshot(5NEMO)\fP
+format . Currently the potential and total energy are stored in the
+Potential and Aux slots of the snapshot. [no default].
+.TP
+\fBnbody=\fInum_bodies\fP
+Number of bodies particles [default: \fB2048\fP]
+.TP
+\fBpotname=\fIpotential_name\fP
+name of the potential, see \fB$NEMOOBJ/potential\fP for the current
+object repository. The user can supply his own, see \fIpotential(5NEMO)\fP.
+[no default].
+.TP
+\fBpotpars=\fIpot_pars\fP
+Paramaters to the user supplied potential. The number of parameters
+depends on the potential supplied, the first parameter is reserved
+for the pattern speed.
+[default: not supplied, parameters as defined by potential(5)].
+.TP
+\fBpotfile=\fIpot_data\fP
+Optional data file(s) to the user supplied potential.
+The number of files depends
+on the potential supplied.
+[default: not supplied, or datafile(s) defined by potential(5)].
+.TP
+\fBrmin=\fImin_disk_radius\fP
+Inner cutoff radius of test-particle disk. [Default: \fB0\fP].
+.TP
+\fBrmax=\fImax_disk_radius\fP
+Outer cutoff radius of test-particle disk. [Default: \fB1\fP].
+.TP
+\fBmass=\fItot_mass\fP
+Total mass of the disk. The default total mass is 0, since this is a
+testdisk. Since there are a few programs that will not like mass-less
+particles, this keyword can be used, or \fIsnapmass\fP
+will have to be used before further action is taken.
+[default: \fB0\fP].
+.TP
+\fBfrac=\fIfraction\fP
+The fraction that the local velocity dispersion has w.r.t.
+the local rotation velocity. If a number is given, radial and tangential velocity are
+by default taken to be the same, with no Z-velocities. If two numbers are given,
+the first one is the radial dispersion, the second one the tangential
+(currently wihtout any 2*Omega/Kappa correction), and the Z dispersion 0.
+Finally, when 3 numbers are given, random motions in Z are also added.
+[Default: \fB0\fP]
+.TP
+\fBseed=\fIrandom_seed\fP
+Use random number seed. If 0 is given, a random number is generated
+from the time of the day. [default: \fB0\fP].
+.TP
+\fBsign=-1|1\fP
+Sign of the angular momentum. 1 means counter clockwise rotation in the
+XY-plane (in case you were wondering,
+our galaxy has sign=-1). [Default: \fB1\fP].
+.TP
+\fBin=\fIinfile\fP
+Input file, in \fIsnapshot(5NEMO)\fP format, of which the positions
+are used as initial conditions. [Default: not used].
+*** \fI lost the implementation ** \fP
+.TP
+\fBangle=t|f\fP
+Normally the particles are distributed randomly in phase, but setting this
+to true, they will be regularly spread between 0 and TWO_PI.
+Default: f
+.TP
+\fBvrad=\fP
+A constant radial velocity (>0 means outward motion) added to all particles.
+Default: 0.
+.TP
+\fBenergy=t|f\fP
+If random motions are added (see frac=), the energy of each particle will
+change. By setting \fBenergy=t\fP the energy is conserved and the "circular"
+speed will be reduced. This way the guiding center for epicyclic orbits
+will also be conserved. Note that this only applies to motion in XY.
+For stars with too much random motion (if \fBfrac\fP too high), the
+random number generator is run again to ensure conserving the planar
+motion kinetic energy. Default: false.
+.TP
+\fBabs=t|f\fP
+Use absolute value of velocity dispersion instead of fractional.
+Default: f
+.TP
+\fBz0=z0_d,z0_v\fP
+Scaleheight for density dropoff, and circular velocity dropoff. Default: 0,0
+.TP
+\fBvloss=\fP
+(deprecated) Fractional loss in circular streaming as z increased. Use z0_v now.
+.TP
+\fBheadline=\fImessage\fP
+Text headline for output file [default: not used].
+
+.SH "EXAMPLES"
+Under development, but here's a quick test that should produce a plot:
+.nf
+ diskspectrum - 1000000 vrange=1 | tabplot - line=1,1
+.fi
+
+.SH "SEE ALSO"
+mkdisk(1NEMO)
+
+.SH "UPDATE HISTORY"
+.nf
+.ta +1.5i +4.5i
+28-nov-2024 V0.1: cloned of mkdisk to check performance PJT/DiWen
+.fi
diff --git a/man/man1/lineid.1 b/man/man1/lineid.1
new file mode 100644
index 000000000..2acdcf089
--- /dev/null
+++ b/man/man1/lineid.1
@@ -0,0 +1,89 @@
+.TH LINEID 1NEMO "15 December 2024"
+
+.SH "NAME"
+lineid \- Identify peaks in a spectrum and optionally find lines
+
+.SH "SYNOPSIS"
+\fBlineid\fP [parameter=value]
+
+.SH "DESCRIPTION"
+\fBlineid\fP identifies peaks in a spectrum, and
+attempts to identify spectral lines. The units of the frequency
+axis will need to be specified (the default is GHz) in order for the line identification to work.
+.PP
+Operation can occur in several modes.
+.br
+1. If the VLSR is given, a LineID is attempted. Either from a given
+restfreq, or from a linelist.
+.br
+2. If the RESTFREQ is given, a VLSR is determined from the frequency of the
+fitted peak.
+.br
+3. If neither is given, just the peak(s) are reported.
+.PP
+See also \fItablsqfit(1NEMO)\fP with the option \fBfit=peak\fP, and \fIccdmom(1NEMO)\fP has
+an algorithm to find (properties of) the Nth peak using mom=3,30,31,32,33,34
+.PP
+More sophisticated algorithms can be find in for example ADMIT's LineID task.
+
+.SH "PARAMETERS"
+.so man1/parameters
+.TP 20
+\fBin=\fP
+Input (table) file name. No default.
+.TP
+\fBxcol=\fP
+X coordinate column, representing the frequency or wavelength axis.
+Units are GHz, unless specified with the \fBxunit\fP keyword.
+[1]
+.TP
+\fBycol=\fP
+Y coordinate column, reprepenting an intensity. [2]
+.TP
+\fBxunit=\fP
+X axis unit (GHz, MHz, A, nm) -- not implemented yet.
+[GHz]
+.TP
+\fBminchan=\fP
+Minimum channels for a peak [3] - not used yet
+.TP
+\fBmaxchan=\fP
+Maximum channels for a peak [] - not used yet
+.TP
+\fBvlsr=\fP
+VLSR of object. Not set by default.
+.TP
+\fBrestfreq=\fP
+If line is known, give restfreq in xunits. Not set by default.
+.TP
+\fBlinelist=\fP
+Table with frequencies and lines. Not set by default.
+
+.SH "EXAMPLES"
+The following example smoothes a set of random values by 10 channels to create a handful
+of peaks and investigates the strongest peak:
+.EX
+
+ tabgen - 115 1 2 123 | tabsmooth - 1 12 10,5| lineid - 0 1
+ ### nemo Debug Info: Peak: 40.1222 0.169591
+
+ tabgen - 115 1 2 123 | tabsmooth - 1 12 10,5| lineid - 0 1 vlsr=10000
+ ### nemo Debug Info: Line at 40.1222, Look near freq = 41.5067 for a line
+
+ tabgen - 115 1 2 123 | tabsmooth - 1 12 10,5| lineid - 0 1 restfreq=115
+ ### nemo Debug Info: Line at 40.1222 has z=0.651112 or vlsr=...
+.EE
+
+.SH "SEE ALSO"
+tablsqfit(1NEMO), table(5NEMO)
+.PP
+admit::LineID_AT()
+
+.SH "AUTHOR"
+Peter Teuben
+
+.SH "UPDATE HISTORY"
+.nf
+.ta +1.25i +4.5i
+13-dec-2024 0.1 Draft created PJT
+.fi
diff --git a/man/man1/nemoinp.1 b/man/man1/nemoinp.1
index 1de9fd8c9..3e3620367 100644
--- a/man/man1/nemoinp.1
+++ b/man/man1/nemoinp.1
@@ -25,6 +25,8 @@ nemoinp 1:10:2
matlab 1:2:10
seq 1 2 10
numpy.arange 1,10,2 last element (10) is never used
+IDL findgen(5)*2+1 -or- indgen(5)*2+1
+
.fi
.SH "PARAMETERS"
@@ -97,7 +99,7 @@ seq(1), tabgen(1NEMO), tabmath(1NEMO), nemoinp(3NEMO), getrange(3NEMO), function
nemoinp 'iflt(1,2,3,4)' print 3, not 4, because 1 is less than 2
nemoinp 30:30:15 dms=t convert hexasegimal 30:30:15 into decimal 30.5043
nemoinp c,G,h,p,m a few good constants
- nemoinp 'p*pi/(3600*180)' 1 AU in km
+ nemoinp 'p*pi/(3600*180)' 1 AU in m
.fi
.SH "AUTHOR"
diff --git a/man/man1/perorb.1 b/man/man1/perorb.1
index f75159dae..2a177b015 100644
--- a/man/man1/perorb.1
+++ b/man/man1/perorb.1
@@ -250,7 +250,7 @@ In the 1st example we establish that energy is roughly conserved (and the value
the center. The 4th example tries to reproduce the periodic orbit with that specified energy.
.SH "SEE ALSO"
-orbstat(1NEMO), potlist(1NEMO), orbint(1NEMO), newton0(1NEMO), potential(5NEMO)
+orbstat(1NEMO), potlist(1NEMO), orbint(1NEMO), newton0(1NEMO), henyey(1NEMO), potential(5NEMO)
.nf
Particle Swarm Optimization (PSO): Skokos et al. - 2005MNRAS.359..251S
.fi
diff --git a/man/man1/quadcode.1 b/man/man1/quadcode.1
index c4cbb9669..0b1f3bef9 100644
--- a/man/man1/quadcode.1
+++ b/man/man1/quadcode.1
@@ -1,9 +1,12 @@
-.TH QUADCODE 1NEMO "6 May 1992"
-.SH NAME
+.TH QUADCODE 1NEMO "20 Oct 2024"
+
+.SH "NAME"
quadcode \- global quadrupole-order N-body code integrator
-.SH SYNOPSIS
+
+.SH "SYNOPSIS"
\fBquadcode\fP \fBin=\fP\fIsnapshot\fP [\fIparameter\fP=\fIvalue\fP] .\|.\|.
-.SH DESCRIPTION
+
+.SH "DESCRIPTION"
\fIquadcode\fP is an equal-timestep implementation of an
N-body code where the forces and potential are computed
from a potential expansion in spherical harmonics
@@ -12,8 +15,9 @@ it multipole expansion). See also a followup in
Bontekoe & van Albada (MNRAS 224, 349 (1987)) and a rebuttal
in Zaritsky & White (MNRAS 235, 289 (1988).
See also Hernquist, L. & Barnes, J. \fIAp.J.\fP \fB349\fP, 562 (1990)
-.SH PARAMETERS
-The following parameters are recognized; they may be given in any order.
+
+.SH "PARAMETERS"
+.so man1/parameters
.TP 24
\fBin\fP=\fIin-file\fP
Snapshot with initial conditions. No default.
@@ -82,25 +86,32 @@ to be recompiled, for example:
(Note that all \fIquad*\fP programs need to (should)
be recompiled if the maximum allowed number of bodies will
be changed.
-.SH SEE ALSO
-quadforce(1NEMO), quadinter(1NEMO), CSG(1NEMO), octcode(1NEMO)
-.SH ADS
+
+.SH "SEE ALSO"
+quadforce(1NEMO), quadinter(1NEMO), CSG(1NEMO), octcode(1NEMO)
+
+.SH "ADS"
@ads 1983ApJ...274...53W closely derived code? or should we link White's 1978MNRAS.184..185W paper
-.SH BUGS
+
+.SH "BUGS"
When a \fBrestart\fP is specified, values must be given for \fBALL\fP
legal parameters which do not take default values.
-.SH AUTHOR
+
+.SH "AUTHOR"
Joshua E. Barnes.
-.SH FILES
+
+.SH "FILES"
.nf
.ta +2i
src/nbody/evolve/multicode/ exported source code
usr/josh/nbody/multicode/ Josh' original source code
.fi
-.SH HISTORY
+
+.SH "HISTORY"
.nf
-.ta +1i +4i
+.ta +1.25i +4.5i
4-mar-89 V1.2 some formal NEMO version JEB
12-nov-91 V1.3 new NEMO V2. location in $NEMO/src tree PJT
6-may-92 document improved PJT
+20-oct-2024 V1.4a gcc-14 prototype fixed PJT
.fi
diff --git a/man/man1/quadforce.1 b/man/man1/quadforce.1
index 425c2984b..2f87ae3e6 100644
--- a/man/man1/quadforce.1
+++ b/man/man1/quadforce.1
@@ -1,17 +1,20 @@
-.TH QUADFORCE 1NEMO "6 May 1992"
-.SH NAME
+.TH QUADFORCE 1NEMO "20 Oct 2024"
+
+.SH "NAME"
quadforce \- quadrupole-order force calculation of an N-body system.
-.SH SYNOPSIS
+
+.SH "SYNOPSIS"
\fBquadforce\fP \fBin=\fPsnapshot [parameter=value] ...
-.SH DESCRIPTION
+
+.SH "DESCRIPTION"
\fIquadforce\fP computes a quadrupole-order force calculation
of a snapshot N-body system. (see also \fIquadcode(1NEMO)\fP). The
output can be either in the form of the snapshot, or a table, in
which case quadrupole-oder forces can be inserted in other snapshots
(see \fIquadinter(1NEMO)\fP).
-.SH PARAMETERS
-The following parameters are recognized in any order if the keyword
-is also given:
+
+.SH "PARAMETERS"
+.so man1/parameters
.TP 20
\fBin=\fP\fIsnaphot\fP
Input file with N-body snapshot. No default.
@@ -28,18 +31,23 @@ Radial softening parameter. [Default: \fB0.05\fP]
.TP 20
\fBeps_t=\fP
Tangential softening parameter. [Default: \fB0.07\fP]
-.SH LIMITATIONS
+
+.SH "LIMITATIONS"
The code has a hardcoded maximum number of particles, through the
macro \fBMBODY\fP. (currently 4096). If more are needed, it needs
to be recompiled. See comments in \fIquadcode(1NEMO)\fP
-.SH AUTHOR
-Joshue Barnes
-.SH SEE ALSO
+
+.SH "AUTHOR"
+Joshua Barnes
+
+.SH "SEE ALSO"
quadcode(1NEMO), quadinter(1NEMO)
-.SH UPDATE HISTORY
+
+.SH "HISTORY"
.nf
-.ta +1.0i +4.0i
+.ta +1.25i +4.5i
4-mar-89 V1.0 Old version JEB
12-nov-91 V1.1 FOr new nemo V2. PJT
6-may-92 V1.1a man page written, added some warnings PJT
+20-oct-2024 V1.4a fix various prototypes for modern gcc PJT
.fi
diff --git a/man/man1/quadinter.1 b/man/man1/quadinter.1
index 67a9158d2..9d8027239 100644
--- a/man/man1/quadinter.1
+++ b/man/man1/quadinter.1
@@ -1,15 +1,18 @@
-.TH QUADINTER 1NEMO "6 May 1992"
-.SH NAME
+.TH QUADINTER 1NEMO "20 Oct 2024"
+
+.SH "NAME"
quadinter \- quadrupole-order force calculation from tabulated field.
-.SH SYNOPSIS
+
+.SH "SYNOPSIS"
\fBquadinter\fP \fBquad=\fIquad_table\fP \fBin=\fIsnapshot\fP [parameter=value] ...
-.SH DESCRIPTION
+
+.SH "DESCRIPTION"
\fIquadinter\fP computes quadrupole-order forces and potentials for a set of
particules supplies in an N-body snapshot from a quad-table that must have
been prepared by \fIquadforce\fP.
-.SH PARAMETERS
-The following parameters are recognized in any order if the keyword
-is also given:
+
+.SH "PARAMETERS"
+.so man1/parameters
.TP 20
\fBquad=\fP\fIquad_table\fP
Input file with quadrupole-order field tables. No default.
@@ -20,18 +23,22 @@ Input file with an N-body snapshot. No default.
\fBout=\fP
output file with added forces and potentials. This will be a copy
of the input \fIsnapshot\fP. Default: none.
-.SH LIMITATIONS
+
+.SH "LIMITATIONS"
The code has a hardcoded maximum number of particles, through the
macro \fBMBODY\fP. (currently 4096). If more are needed, it needs
to be recompiled. See comments in \fIquadcode(1NEMO)\fP
-.SH SEE ALSO
+
+.SH "SEE ALSO"
quadcode(1NEMO), quadforce(1NEMO)
-.SH AUTHOR
+
+.SH "AUTHOR"
Joshua Barnes
-.SH UPDATE HISTORY
+.SH "HISTORY"
.nf
-.ta +1.0i +4.0i
+.ta +1.25i +4.5i
4-mar-89 V1.0 Old version JEB
12-nov-91 V1.1 For new nemo V2. PJT
6-may-92 V1.1a man page written PJT
+20-oct-2024 V1.4a gcc-14 prototype fixes PJT
.fi
diff --git a/man/man1/reb2s.1 b/man/man1/reb2s.1
index b04276dcb..78be411a1 100644
--- a/man/man1/reb2s.1
+++ b/man/man1/reb2s.1
@@ -1,4 +1,4 @@
-.TH REB2S 1NEMO "4 July 2024"
+.TH REB2S 1NEMO "8 July 2024"
.SH "NAME"
reb2s \- convert REBOUND SimulationArchive to NEMO snapshot
@@ -40,7 +40,8 @@ Since verbose output is on standard out, the following pipe will fail in unpredi
.EE
since now NEMO's snapshot as well as the verbose output are in the pipe.
.PP
-
+Cannot select by time yet. This should be implemented, see \fIsnapcopy(1NEMO)\fP or \fIsnaptrim(1NEMO)\fP
+for a NEMO filter.
.SH "SEE ALSO"
s2reb(1NEMO), reb_viewer(1NEMO), reb_integrate(1NEMO), snapshot(5NEMO), rebound(8NEMO)
diff --git a/man/man1/rotcur.1 b/man/man1/rotcur.1
index f9e98fb0e..2d8bdd811 100644
--- a/man/man1/rotcur.1
+++ b/man/man1/rotcur.1
@@ -349,7 +349,7 @@ P. J. Teuben (NEMO C version)
.SH "UPDATE HISTORY"
.nf
-.ta +1i +4i
+.ta +1.5i +7.0i
19/jul/83 original program KGB
9/mar/85 revision of program KGB
23/may/86 migrated to VAX-VMS KGB
diff --git a/man/man1/snapbench.1 b/man/man1/snapbench.1
index 6b5cadf51..caa8aebaa 100644
--- a/man/man1/snapbench.1
+++ b/man/man1/snapbench.1
@@ -1,4 +1,4 @@
-.TH SNAPBENCH 1NEMO "12 March 2024"
+.TH SNAPBENCH 1NEMO "9 July 2024"
.SH "NAME"
snapbench \- benchmark (re)assigning masses to a snapshot
@@ -14,6 +14,13 @@ It also uses two memory models of keeping a snapshot:
an array of structs (AOS) vs. a struct of arrays (SOA). snapshots'
are normally stored on disk as SOA's, but most programs actually
transpose them to an AOS, which can be very in-efficient.
+.PP
+A hidden code \fIsnaprun\fP compares defining simple variables (x,y,z) in a
+structure with an array pos[3], and also looping over pos[] insted of
+hardcoding a simpler Eulerian type of update. Counter intuitive, looping
+made the code about 4-8 times slower. Simple variables are slightly faster
+than using directly addressing the array. Results vary a bit on different
+processors.
.SH "PARAMETERS"
.so man1/parameters
@@ -63,11 +70,13 @@ with the following examples
Only reading and assignment is measured. No options for output (if out= is given)
.SH "SEE ALSO"
-snapshot(5NEMO)
+snaprun, snapshot(5NEMO)
.SH "FILES"
+.nf
src/nbody/trans/snapbench.c
-
+src/tutor/nbody/snaprun.c
+.fi
.SH "AUTHOR"
Peter Teuben
diff --git a/man/man1/tabcols.1 b/man/man1/tabcols.1
index 07c6946d0..f27cfca0f 100644
--- a/man/man1/tabcols.1
+++ b/man/man1/tabcols.1
@@ -1,4 +1,4 @@
-.TH TABCOLS 1NEMO "20 November 2022"
+.TH TABCOLS 1NEMO "10 September 2024"
.SH "NAME"
tabcols \- Select columns, and optionally rearrange, from a table
@@ -7,13 +7,13 @@ tabcols \- Select columns, and optionally rearrange, from a table
\fPtabcols\fP [parameter=value] ...
.SH "DESCRIPTION"
-\fBtabcols\fP selects columns from a file, and can also re-arrange the
+\fBtabcols\fP selects columns from a file, and can thus also re-arrange the
output column order. New column separators can be used on output.
This way a newline on input can be converted into a space on
output and therefore creating a file with one long column from
-a table with mulitple columns.
+a table with multiple columns.
.PP
-To select specific lines, use \fItabrows(1NEMO)\fP or \fItabmath(1NEMO)\fP.
+To select specific rows/lines, use \fItabrows(1NEMO)\fP or \fItabmath(1NEMO)\fP.
.PP
To align a table by column, use the \fIcolumn(1)\fP program with the \fB-t\fP
option.
@@ -22,13 +22,16 @@ option.
.so man1/parameters
.TP 20
\fBin=\fP
-input file name(s). No default.
+input file name(s).
+.br
+No default.
.TP
\fBselect=\fP
list of columns to select. Any \fInemoinp(1NEMO)\fP expression can be used.
Columns can be repeated, and in any order. The output will be in exactly
the same order as given with this keyword. First column is 1. Column 0 can also
be added to the output, which is the row number (first row is 1).
+.br
[all]
.TP
\fBcolsep=\fP
@@ -36,14 +39,24 @@ Column separator. Valid are: \fBs\fP (space), \fBn\fP (newline),
\fBt\fP (tab),
\fBr\fP (carriage return), as well as the literal symbols '\fB,\fP'
and '\fB:\fP'.
+.br
[Default: \fBs\fP]
+
+.TP
+\fBcolsepin=\fP
+Column separator for input. By default an internal set (space, comma, tab, vertical bar)
+is used. A string can be specified here consistent of multiple characters.
+
+
.TP
\fBout=-\fP
output file name. By default standard output it used.
+.br
+[-]
.SH "CAVEATS"
Although this program uses the new table interface, there are
-still some internal variables that limits the column count.
+still some internal variables that limit the column count.
.SH "EXAMPLE"
Here is an example of how to reverse the columns from a table with 4 columns:
@@ -55,9 +68,13 @@ re-arranges all the numbers in a table in one long column
.nf
tabcols in.tab all n > out.tab
.fi
+or in one long row:
+.nf
+ tabcols in.tab all n | tabtranspose - - > out.tab
+.fi
.SH "SEE ALSO"
-tabrows(1NEMO), tabmath(1NEMO), tabs(1NEMO), tabtab(1NEMO), awk(1), column(1)
+tabrows(1NEMO), tabtranspos(1NEMO), tabmath(1NEMO), tabs(1NEMO), tabtab(1NEMO), awk(1), column(1)
.SH "FILES"
src/kernel/tab/tabcols.c source code
@@ -71,4 +88,5 @@ Peter Teuben
26-Jan-00 V1.0 Created PJT
9-apr-09 V2.0 out=- now default and last argument PJT
5-may-2022 V2.1 converted to table V2 interface PJT
+10-sep-2024 V2.4 added colsepin= PJT
.fi
diff --git a/man/man1/tabnllsqfit.1 b/man/man1/tabnllsqfit.1
index 2da8e515a..d97b59f41 100644
--- a/man/man1/tabnllsqfit.1
+++ b/man/man1/tabnllsqfit.1
@@ -1,10 +1,13 @@
-.TH TABNLLSQFIT 1NEMO "9 October 2013"
-.SH NAME
+.TH TABNLLSQFIT 1NEMO "3 January 2024"
+
+.SH "NAME"
tabnllsqfit \- general purpose non-linear least squares fitting program
-.SH SYNOPSIS
+
+.SH "SYNOPSIS"
.PP
\fBtabnllsqfit in=\fPtable-file [parameter=value]
-.SH DESCRIPTION
+
+.SH "DESCRIPTION"
\fBtabnllsqfit\fP takes a non-linear least squares fit (minimized ChiSquared)
of a set of datapoints (x(i),y(i)) where a functional
form y=f(x;p0,p1,...,pN) is assumed. The function is not required to
@@ -31,9 +34,8 @@ the residuals are still computed correctly.
This routine can also be used to plot the function, although only in 1D cases,
see \fBx=\fP.
-.SH PARAMETERS
-The following parameters are recognized in any order if the keyword is also
-given:
+.SH "PARAMETERS"
+.so man1/parameters
.TP 20
\fBin=\fIin-file\fP
(ascii) input file, a table of values from which data is taken. No default.
@@ -155,7 +157,7 @@ Fitting method. Gipsy uses their \fInllsqfit(3NEMO)\fP, NumRec uses their
\fImrqfit\fP and MINPACK uses \fImpfit(3NEMO)\fP. Only first character
is used, case insensitive. [Default: \fBg\fP]
-.SH FIT PARAMETERS
+.SH "FIT PARAMETERS"
Parameters are referred to as p0,p1,p2,p3,.....
.nf
.ta +1.5i
@@ -171,7 +173,8 @@ arm3 y=p0+p1*cos(x)+p2*sin(x)+p3*cos(3*x)+p4*sin(3*x)
loren y=p1/PI( (x-p0)^2 + p1^2 )
psf y=p1*x^p2*sin(y)^p3+p4
.fi
-.SH EXAMPLE
+
+.SH "EXAMPLES"
Here is an example of a linear fit: a straight
line, with some added noise and random weights between 1 and 2:
.nf
@@ -259,7 +262,25 @@ bootstrap= 3.11603 0.164922 1.96464 0.086429 0.00293632 0.00820007
P0 dP0 P1 dP2 P3 dP3
.fi
-.SH LOAD FUNCTIONS
+.PP
+Here is an example of calling tabnllsqfit twice, once to fit a baseline, extract the
+fit parameters, and a second
+time to subtract the baseline using the parameters without a fit and using out=.
+The input table has two columns, velocity
+and intensity, but a bad baseline between -400 and +500. There is a spectral line
+between 0 and 25 km/s, so we fit a local high order baseline between -100 and +100 but
+excluding the line:
+.EX
+
+tabplot MonR2_110433__0.txt line=1,1
+tabnllsqfit MonR2_110433__0.txt xrange=-100:-20,45:150 fit=poly order=3 > pars0
+p0=$(txtpar pars0 p0=p0=,1,2 p1=p1=,1,2 p2=p2=,1,2 p3=p3=,1,2)
+tabnllsqfit MonR2_110433__0.txt fit=poly order=3 out=- free=0,0,0,0 par="$p0" | tabcomment - > tab0
+tabplot tab0 1 3 -50 50 -1 15 line=1,1 ycoord=0
+
+.EE
+
+.SH "LOAD FUNCTIONS"
With the \fBload=\fP keyword dynamic object files can be loaded using the
\fIloadobj(3NEMO)\fP mechanism. The convention is that two functions
must be externally visible, and named \fIfunc_\fP\fImethod\fP and
@@ -295,7 +316,8 @@ One word of caution: if you find the program having a hard time finding
a solution in complex cases, it is quite possible that this is not due to
the fact that the function is complex, but due to noise or bad initial
conditions.
-.SH DERIVATIVES CHECK
+
+.SH "DERIVATIVES CHECK"
A common debug problem is to check the analytical derivatives in
your external loaded function. The \fBx=\fP keyword can be used
to check the function values (useful to plot up the function
@@ -336,10 +358,12 @@ is incremented by a decreasing factor of 0.1
You can see that the offset (par 0) and amplitude (par 1) of the gauss are well behaved (as they are linear),
but the centroid (par 2) and width (par 3) of the gauss only slowly converge linearly with the step. The last
number in square brackets is the error between numerical and analytical derivative.
-.SH CAVEATS
+
+.SH "CAVEATS"
It will not recognize linear fits if the non-linear parameters are kept fixed,
e.g. the offset 0 in \fBfit=gauss1d\fP.
-.SH SEE ALSO
+
+.SH "SEE ALSO"
tablsqfit(1NEMO), tabhist(1NEMO), tabmath(1NEMO),
gaussfit(1NEMO), linreg(1NEMO), nllsqfit(3NEMO), fit.dc1(GIPSY)
.PP
@@ -357,17 +381,19 @@ A new scheme for calculating weights and describing
correlations in nonlinear least squares fits.
(Hessler, Curent & Ogren, C.I.P. 10, 186, 1996):
http://dx.doi.org/10.1063/1.168569
-.SH AUTHOR
+
+.SH "AUTHOR"
Peter Teuben
-.SH FILES
+
+.SH "FILES"
.nf
.ta +2.5i
~/src/kernel/tab tabnllsqfit.c
~/src/kernel/tab/fit example fitting functions
.fi
-.SH "UPDATE HISTORY"
+.SH "HISTORY"
.nf
-.ta +1.0i +4.0i
+.ta +1.5i +5.5i
12-jul-02 V1.0 cloned off tablsqfit PJT
17-jul-02 V1.1 added load=, x=, numrec= PJT
11-sep-02 V1.1e changes error/warning to accomodate residual writen PJT
diff --git a/man/man3/timers.3 b/man/man3/timers.3
index 4da6afba1..6e1d0a461 100644
--- a/man/man3/timers.3
+++ b/man/man3/timers.3
@@ -38,12 +38,12 @@ could be useful.
stamp_timers(i);
stamp_timers(n);
for (i=0; i SCA
r1=0.1 # central unresolved bulge, bar or black hole
v1=0 # representative rotation speed at r1
re=20 # exponential scalelength of disk (arcsec) #> SCALE 1:100:1
+n=-1 # n<1 exp disk; n>=0 PLEC disk #> SCALE -1:10:0.1
+rm=-1 # if > 0: use it has R_mol^2 - AMT model #> SCALE -1:50:0.1
rmax=60 # edge of disk (arcsec) #> SCALE 1:80:1
inc=60 # INC of disk #> SCALE 0:90:1
z0=0 # scaleheight [@todo buggy] #> SCALE 0:10:0.5
@@ -34,6 +42,8 @@ sigma=0 # random motion in plane (km/s) #> SCA
noise=0 # add optional noise to cube #> ENTRY
vlsr=0 # optional VLSR if non-zero #> ENTRY
+debug=-1 # add debugging #> SCALE -1:9:1
+
plot=profile,rotcur # show which plots can be made #> CHECK rotcur,profile,density
#--HELP
@@ -51,11 +61,18 @@ for arg in "$@"; do
export "$arg"
done
+# handle debugging
+debug=$(nemoinp $debug format=%d)
+if [ $debug -gt 0 ]; then
+ set -x
+fi
+export DEBUG=$debug
+echo DEBUG=$DEBUG
+
# fixed now compared to edge_aca.sh
range=$rmax # gridding size
nsize=1 # number of spatial pixels (px=py=2*range/nx)
-pa=0 # PA of receding side disk (E through N)
# derive some parameters that appear common or logically belong together
@@ -64,13 +81,22 @@ grid_pars="xrange=-${range}:${range} yrange=-${range}:${range} nx=$nsize ny=$nsi
cell=`nemoinp "2*$range/$nsize*60"`
cen=`nemoinp $nsize/2-0.5`
restfreq=230.53800 # CO(2-1) in GHz
-nbody=`nemoinp 10**$logn format=%d`
+nbody=`nemoinp "10**$logn" format=%d`
+sininc=`nemoinp "sind($inc)"`
+noise=`nemoinp $noise/$nbody` # kludge
+
+function nemo_stamp {
+ if [ $debug -ge 0 ]; then
+ echo "$(date +%H:%M:%S.%N) $*"
+ fi
+}
# ================================================================================ START
# Announce:
echo "$0 version $_version"
+nemo_stamp start
# Clear old model
rm -f $run.* >& /dev/null
@@ -78,74 +104,99 @@ rm -f $run.* >& /dev/null
# keep a log
echo "`date` :: $*" > $run.history
+
+mmode=$(nemoinp "ifge($n,0,1,0)")
+if [ $mmode = 0 ]; then
+ mass="exp(-r/$re)"
+else
+ mass="pow(r/$re*exp(1-r/$re),$n)"
+fi
+mmode=$(nemoinp "ifge($rm,0,2,$mmode)")
+if [ $mmode = 2 ]; then
+ mass="exp(-r/$re)/(1+$rm*exp(-1.6*r/$re))"
+fi
+echo MMODE=$mmode MASS=$mass
+echo NBODY=$nbody
+nemo_stamp begin
if [ $m0 != 0 ]; then
- echo "Creating brandt disk with $nbody particles times"
+ echo "Creating brandt disk with $nbody particles."
mkdisk out=- nbody=$nbody seed=$seed z0=$z0,$z0 \
potname=brandt potpars=0,$v0,$r0,$m0 mass=1 sign=-1 frac=$sigma abs=t rmax=$rmax |\
- snapmass - - "mass=exp(-r/$re)" norm=1 |\
- snaprotate - $run.20 "$inc,$pa" yz
+ snapmass - $run.20 "mass=$mass" norm=1
+ nemo_stamp mkdisk
rotcurves name1=brandt pars1=0,$v0,$r0,$m0 tab=t radii=0:${rmax}:0.01 plot=f |\
tabmath - - "%1,1,%2,$sigma,exp(-%1/$re)" all > $run.tab
+ nemo_stamp rotcurves
vmax=$(tabstat $run.tab 3 | txtpar - p0=max:,1,2)
echo "Max rotcur: vmax=$vmax m0=$m0"
+ nemo_stamp tabstat
elif [ $v1 = 0 ]; then
- echo "Creating homogeneous disk with $nbody particles times"
+ echo "Creating homogeneous disk with $nbody particles. (need optimizing test)"
mkdisk out=- nbody=$nbody seed=$seed z0=$z0,$z0 \
potname=rotcur0 potpars=0,$v0,$r0 mass=1 sign=-1 frac=$sigma abs=t rmax=$rmax |\
- snapmass - - "mass=exp(-r/$re)" |\
- snapscale - - mscale=$nbody |\
- snaprotate - $run.20 "$inc,$pa" yz
+ snapmass - - "mass=$mass" |\
+ snapscale - $run.20 mscale=$nbody
rotcurves name1=rotcur0 pars1=0,$v0,$r0 tab=t radii=0:${rmax}:0.01 plot=f |\
tabmath - - %1,1,%2,$sigma all > $run.tab
else
m1=`nemoinp "$r1*($v1/0.62)**2"`
- echo "Creating homogeneous disk with $nbody particles times and a nuclear component m1=$m1"
+ echo "Creating homogeneous disk with $nbody particles and a nuclear component m1=$m1 (needs optimizing test)"
rotcurves name1=rotcur0 pars1=0,$v0,$r0 name2=plummer pars2=0,$m1,$r1 tab=t radii=0:${rmax}:0.01 |\
tabmath - - %1,1,%2,$sigma all > $run.tab
mktabdisk $run.tab - nbody=$nbody seed=$seed rmax=$rmax sign=-1 |\
- snapmass - - "mass=exp(-r/$re)" |\
- snapscale - - mscale=$nbody |\
- snaprotate - $run.20 "$inc,$pa" yz
+ snapmass - - "mass=$mass" |\
+ snapscale - $run.20 mscale=$nbody
fi
-echo "Creating velocity moments"
-snapgrid $run.20 $run.21 $grid_pars moment=0 evar=m
-snapgrid $run.20 $run.22 $grid_pars moment=1 evar=m
-snapgrid $run.20 $run.23 $grid_pars moment=2 evar=m
-# skip smoothing
-ccdmath $run.21,$run.22,$run.23 $run.20d %1
-ccdmath $run.21,$run.22,$run.23 $run.20v "%2/%1"
-ccdmath $run.21,$run.22,$run.23 $run.20s "sqrt(%3/%1-%2*%2/(%1*%1))"
-# clipping
-ccdmath $run.20d,$run.21 $run.21d "ifgt(%2,0,%1,0)"
-ccdmath $run.20v,$run.21 $run.21v "ifgt(%2,0,%1,0)"
-ccdmath $run.20s,$run.21 $run.21s "ifgt(%2,0,%1,0)"
-
+#snapsort $run.20 - r help=c | snapshell - radii=0:${rmax}:0.01 pvar=m help=c > $run.shell
+snapshell $run.20 radii=0:${rmax}:0.01 pvar=m help=c > $run.shell
+nemo_stamp snapshell
echo "Creating a velocity field - method 2"
-snapgrid $run.20 - $grid_pars \
- zrange=-${vrange}:${vrange} nz=$nvel mean=f evar=m |\
- ccdflip - $run.30 flip=$flip wcs=t
+snapgrid $run.20 $run.30 $grid_pars \
+ zrange=-${vrange}:${vrange} nz=$nvel mean=f evar=m zvar=vy*$sininc
+nemo_stamp snapgrid
ccdstat $run.30
+nemo_stamp ccdstat
+
if [ $vbeam = 0 ]; then
ccdmath $run.30 $run.31 "%1+rang(0,$noise)"
else
ccdmath $run.30 - "%1+rang(0,$noise)" | ccdsmooth - $run.31 $vbeam dir=z
fi
+nemo_stamp ccdmath
# single dish profile
+if [ 1 = 1 ]; then
ccdmom $run.31 - axis=1 mom=0 |\
ccdmom - - axis=2 mom=0 |\
ccdprint - x= y= z= label=z newline=t |\
tabcomment - - punct=f delete=t > $run.spec
+else
+ # faster, but normalization is off now
+ ccdprint $run.31 x= y= z= label=z newline=t |\
+ tabcomment - - punct=f delete=t > $run.spec
+fi
+nemo_stamp ccdprint
echo -n "Total integral over spectrum: "
tabint $run.spec
+nemo_stamp tabint
+echo -n "Estimate of signal/noise in spectrum: "
+peak=$(tabstat $run.spec 2 qac=t | txtpar - p0=QAC,1,6)
+rms=$(tabtrend $run.spec 2 | tabstat - qac=t robust=t | txtpar - p0=QAC,1,4)
+if [ $rms = 0 ]; then
+ echo "$peak / 0 = 999999"
+else
+ echo "$peak / $rms = $(nemoinp $peak/$rms)"
+fi
+
-# export for other programs, in decent units
+# export to FITS, in decent units
# this way the input spatial scale is in arcsec and km/s
-ccdfits $run.31 $run.fits radecvel=t scale=1/3600.0,1/3600.0,1.0 crval=$crval restfreq=$restfreq
+# SKIP for production
+# ccdfits $run.31 $run.fits radecvel=t scale=1/3600.0,1/3600.0,1.0 crval=$crval restfreq=$restfreq
echo "PLOT: $plot"
@@ -153,9 +204,10 @@ if [[ "$plot" == *"rotcur"* ]]; then
tabplot $run.tab 1 3 headline="Rotation Curve" yapp=1/xs
fi
if [[ "$plot" == *"profile"* ]]; then
- tabplot $run.spec line=1,1 headline="Spectrum around VLSR=$vlsr" yapp=2/xs
+ tabplot $run.spec line=1,1 ycoord=0 headline="Spectrum around VLSR=$vlsr" yapp=2/xs
fi
if [[ "$plot" == *"density"* ]]; then
- tabplot $run.tab 1 5 line=1,1 ymin=0 headline="Surface Density" yapp=3/xs
+ tabplot $run.shell 1 4 line=1,1 ymin=0 headline="Surface Density" yapp=3/xs
fi
+nemo_stamp end
diff --git a/src/image/misc/ccdbench.c b/src/image/misc/ccdbench.c
index f9fcd53b1..1637a5841 100644
--- a/src/image/misc/ccdbench.c
+++ b/src/image/misc/ccdbench.c
@@ -12,13 +12,15 @@
#include
#include
#include
+#include
string defv[] = {
"in=???\n Input image file",
"order=xyz\n Order to loop over cube, where the last axis is innermost loop",
"iter=1\n How often to loop",
"mode=0\n 0=CubeValue(i,j,k) 1=cube[i][k][k]",
- "VERSION=0.2\n 19-mar-2024 PJT",
+ "test=f\n Show memory addresses for a small 2x3 array",
+ "VERSION=0.3\n 13-apr-2024 PJT",
NULL,
};
@@ -27,6 +29,8 @@ string usage = "benchmark taking moments along an axis of an image";
#define LOOP(i,n) for(i=0;i>2) or pointy (<<1)",
"VERSION=0.6\n 6-dec-2022 PJT",
@@ -168,6 +168,7 @@ void nemo_main(void)
xpos = (Nx(velptr)-1.0)/2.0;
ypos = (Ny(velptr)-1.0)/2.0;
}
+ dprintf(0,"Center pixel: %g %g\n",xpos,ypos);
pa = getdparam("pa");
inc = getdparam("inc");
vsys = getdparam("vsys");
diff --git a/src/image/trans/ccdfill.c b/src/image/trans/ccdfill.c
index c95227a2a..b114c3870 100644
--- a/src/image/trans/ccdfill.c
+++ b/src/image/trans/ccdfill.c
@@ -9,6 +9,7 @@
* 2-jun-03 1.4 trying out xmirror/ymirror pjt
* 1-aug-05 1.5 make it work on cubes too pjt
* 2-may-2017 2.0 also allow deviate pixels values pjt
+ * 7-dec-2023 2.1 allow mask= to designate bad areas pjt
*
*/
@@ -23,14 +24,15 @@ string defv[] = {
"in=???\n Input image file",
"out=???\n Output image file",
"n=1\n Number of neighbor cells on all sides to use",
- "bad=0.0\n Value of a bad pixel to be patched",
+ "bad=0.0\n Value of a bad pixel to be patched (see also mask=)",
"all=f\n Force all points to be refitted?",
"m=3\n Minimum number of neighbor pixels needed",
"iter=1\n Number of iterations",
"xmirror=f\n Use points mirrored in X to get near a border",
"ymirror=f\n Use points mirrored in Y to get near a border",
"spike=\n Spike value above the neighbors to be filled",
- "VERSION=2.0\n 2-may-2017 PJT",
+ "mask=\n If given, this file designates the bad pixels", // @todo like maskfill?
+ "VERSION=2.1\n 7-dec-2023 PJT",
NULL,
};
diff --git a/src/image/trans/ccdpixelclip.c b/src/image/trans/ccdpixelclip.c
index 12d160278..18786ee0e 100644
--- a/src/image/trans/ccdpixelclip.c
+++ b/src/image/trans/ccdpixelclip.c
@@ -13,7 +13,7 @@
* state. It restores dynamic range to that of the camera and recovers
* the blown-out star cores.
*
- * 24-aug-2022 drafted - Peter Teuben
+ * 24-aug-2022 drafted for NEMO - Peter Teuben
*
*/
@@ -31,7 +31,7 @@ string defv[] = {
"size=1\n size around clip pixel to average around (-size..+size)",
"matchclip=t\n Also fix data where it matches the clip value?",
"jwst=t\n JWST mode? (only true is accepted now)",
- "VERSION=0.2\n 24-aug-2022 PJT",
+ "VERSION=0.3\n 25-aug-2022 PJT",
NULL,
};
@@ -90,7 +90,7 @@ void nemo_main()
}
}
dprintf(0,"Clipped %d/%d %g%% values %s %g\n",
- countclip,nx*ny*nz, (1.0*countclip)/(nx*ny*nz),
+ countclip,nx*ny*nz, (100.0*countclip)/(nx*ny*nz),
Qjwst ? (Qmatch ? "equal or below" : "below") :
(Qmatch ? "equal or above" : "above"),
clip);
diff --git a/src/kernel/cores/allocate.c b/src/kernel/cores/allocate.c
index 2476cd9d3..60d6e0146 100644
--- a/src/kernel/cores/allocate.c
+++ b/src/kernel/cores/allocate.c
@@ -23,6 +23,9 @@
#include
#include
+/* _FL are versions that report File and Line ; see WD notes in stdinc.h */
+/* @todo consider mm_malloc() for efficiency ? */
+
void *allocate_FL(size_t nb, const_string file, int line)
{
void *mem;
diff --git a/src/kernel/io/nemovar.c b/src/kernel/io/nemovar.c
index 134eec114..85ef2b7fd 100644
--- a/src/kernel/io/nemovar.c
+++ b/src/kernel/io/nemovar.c
@@ -11,7 +11,8 @@
string defv[] = {
"var=\n Name of variable(s)",
"val=\n Optional new value",
- "VERSION=0.4\n 25-jun-2024 PJT",
+ "comment=\n Optional comment when new value was set",
+ "VERSION=0.5\n 26-jun-2024 PJT",
NULL,
};
@@ -25,10 +26,14 @@ void nemo_main()
bool hasVar = hasvalue("var");
string value = getparam("val");
bool hasVal = hasvalue("val");
+ string comment = getparam("comment");
+ bool hasComment = hasvalue("comment");
string nemovar = getenv("NEMOVAR");
if (nemovar == NULL) error("$NEMOVAR was not set");
dprintf(1,"nemovar: %s\n", nemovar);
+ if (hasComment)
+ warning("comment %s not yet used", comment);
// bug or feature, nemovar currently needs to exist
if (!fexist(nemovar)) {
diff --git a/src/kernel/misc/xrandom.c b/src/kernel/misc/xrandom.c
index 3e69aa3e3..0a1ecaef7 100644
--- a/src/kernel/misc/xrandom.c
+++ b/src/kernel/misc/xrandom.c
@@ -190,6 +190,8 @@ double xrandom(double xl, double xh)
{
double retval;
+ if (xl==xh) return xl;
+
for(;;) {
#if defined(HAVE_GSL)
if (my_r == NULL) error("GSL init_xrandom was never called");
@@ -237,6 +239,8 @@ double grandom(double mean, double sdev)
static double v1, v2, s;
static int gcount = 0;
+ if (sdev <= 0) return mean;
+
if (gcount) {
gcount = 0;
return mean + sdev * v2 * s;
diff --git a/src/kernel/tab/lineid.c b/src/kernel/tab/lineid.c
new file mode 100644
index 000000000..224c5ec53
--- /dev/null
+++ b/src/kernel/tab/lineid.c
@@ -0,0 +1,378 @@
+/*
+ * LINEID: draft
+ *
+ * 13-dec-2024 0.2 simple brightest peak finder.
+ */
+
+/**************** INCLUDE FILES ********************************/
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#define MAXCOL 2
+
+/**************** COMMAND LINE PARAMETERS **********************/
+
+string defv[] = { /* DEFAULT INPUT PARAMETERS */
+ "in=???\n Input (table) file name",
+ "xcol=1\n x coordinate column",
+ "ycol=2\n y coordinate column",
+ "xunit=GHz\n X axis unit (GHz, MHz, A, nm)",
+ "minchan=3\n Minimum channels for a peak",
+ "maxchan=\n Maximum channels for a peak",
+ "vlsr=\n VLSR of object, if known (km/s)",
+ "restfreq=\n RESTFREQ in xunits, if known",
+ "linelist=\n ASCII linelist (freq,label)",
+ "rms=\n Do not fit peaks when this RMS is reached",
+ "VERSION=0.2\n 13-dec-2024 PJT",
+ NULL
+};
+
+string usage = "lineid draft";
+
+
+/**************** GLOBAL VARIABLES ************************/
+
+local string input; /* filename */
+local stream instr; /* input file */
+local table *tptr; /* table */
+local mdarray2 d2; /* data[col][row] */
+
+local int xcol, ycol;
+
+local real *x, *y; /* data from file */
+local int npt; /* actual number of data points */
+
+local bool Qvlsr, Qrest;
+local real vlsr;
+local real restfreq;
+local string linelist = NULL;
+
+#define UNIT_MHZ 1
+#define UNIT_GHZ 2
+#define UNIT_A 3
+#define UNIT_nm 4
+
+void setparams(void);
+void read_data(void);
+void peak_data(void);
+void do_peak(real *xpeak, real *xerr, real *ypeak);
+
+void nemo_main()
+{
+ warning("Program in draft stage, don't expect it to be correct");
+ setparams();
+ read_data();
+ peak_data();
+}
+
+void setparams()
+{
+ input = getparam("in"); /* input table file */
+ instr = stropen (input,"r");
+ tptr = table_open(instr,0);
+
+ xcol = getiparam("xcol");
+ ycol = getiparam("ycol");
+
+ Qvlsr = hasvalue("vlsr");
+ if (Qvlsr) vlsr = getdparam("vlsr");
+
+ Qrest = hasvalue("restfreq");
+ if (Qrest) restfreq = getdparam("restfreq");
+
+ if (hasvalue("linelist"))
+ linelist = getparam("linelist");
+}
+
+#define MVAL 64
+#define MLINELEN 512
+
+void read_data(void)
+{
+ int colnr[1+MAXCOL];
+
+ dprintf(2,"Reading datafile, xcol,ycol=%d..,%d,...\n",xcol,ycol);
+
+ colnr[0] = xcol;
+ colnr[1] = ycol;
+ d2 = table_md2cr(tptr, 2, colnr,0,0);
+ npt = table_nrows(tptr);
+ dprintf(0,"Found %d channels\n", npt);
+
+ x = &d2[0][0];
+ y = &d2[1][0];
+
+
+
+ // get minmax in X and Y
+}
+
+void peak_data(void)
+{
+ real xpeak, xerr, ypeak;
+ real z;
+
+ dprintf(1,"C(mks)=%g\n",c_MKS);
+
+ // simple brightest peak finder
+ do_peak(&xpeak, &xerr, &ypeak);
+
+ if (Qrest) {
+ z = 1 - xpeak/restfreq; // @todo fix
+ dprintf(0,"Line at %g has z=%g or vlsr=...\n",xpeak,z);
+ } else if (Qvlsr) {
+ z = vlsr/c_MKS*1000.0;
+ restfreq = xpeak / (1-z); // @todo fix
+ dprintf(0,"Line at %g, Look near freq = %g for a line\n", xpeak, restfreq);
+ } else {
+ dprintf(0,"Peak: %g %g\n", xpeak, ypeak);
+ }
+
+}
+
+
+/* from: tablsqfit */
+void do_peak(real *xpeak, real *xerr, real *ypeak)
+{
+ real mat[(MAXCOL+1)*(MAXCOL+1)], vec[MAXCOL+1], sol[MAXCOL+1], a[MAXCOL+2];
+ int i, j, k, range=1;
+ int order = 2;
+
+ for (i=1, j=0; i npeak+3) break;
+
+ for(i=0; i peakvalue) { /* find largest values and remember where it was */
+ peakvalue = data[i];
+ apeak = i;
+ }
+ }
+ }
+ dprintf(1,"apeak=%d\n",apeak);
+ if (apeak==0) { /* never allow first point to be the peak ... */
+ mask[apeak] = 0;
+ apeak = -1;
+ for (i=1, oldvalue=peakvalue; mask[i]==0 && i oldvalue) break; /* util data increase again */
+ oldvalue = data[i];
+ mask[i] = npeak;
+ }
+ continue;
+ }
+ if (apeak==(n-1)) { /* ...or last point, since they have no neighbors */
+ mask[apeak]= 0;
+ apeak = -1;
+ // compiler: peakvalue may be used uninitialized
+ for (i=n-2; oldvalue=peakvalue, mask[i]==0 && i>0; i--) { /* walk down, and mask out */
+ if (data[i] > oldvalue) break; /* until data increase again */
+ oldvalue = data[i];
+ mask[i] = npeak;
+ }
+ // break;
+ continue;
+ }
+ if (apeak > 0) { /* done, found a peak, but mask down the peak */
+ mask[apeak] = npeak;
+ for (i=apeak+1, oldvalue=peakvalue; mask[i]==0 && i oldvalue) break;
+ oldvalue = data[i];
+ mask[i] = npeak;
+ }
+ for (i=apeak-1, oldvalue=peakvalue; mask[i]==0 && i>=0; i-- ) {
+ if (data[i] > oldvalue) break;
+ oldvalue = data[i];
+ mask[i] = npeak;
+ }
+ break;
+ }
+
+ } /* while() */
+ return apeak;
+}
+
+static inline int outside(real x,int i0,int i1) {
+ if(xi1) return 0;
+ return 1;
+}
+
+local void peak_assign(int n, real *d, int *s)
+{
+ int i;
+ real p;
+
+ for (i=1; i0) {
+ dprintf(1,"re-assign-1 @%d from %d to %d\n",i,s[i],s[i-1]);
+ s[i] = s[i-1]; /* re-assign */
+ }
+ } else {
+ p = peak_spectrum(n,d,i-1);
+ //dprintf(0,"assign-2: %d %d %d %g %g %g\n",i,s[i],s[i-1],p,d[i-1],d[i]);
+ if (p<0) {
+ dprintf(1,"re-assign-2 @%d from %d to %d\n",i,s[i-1],s[i]);
+ s[i-1] = s[i]; /* re-assign */
+ }
+ }
+ }
+ }
+}
+
+
diff --git a/src/kernel/tab/tabcols.c b/src/kernel/tab/tabcols.c
index 0fdc95934..d8a2179dd 100644
--- a/src/kernel/tab/tabcols.c
+++ b/src/kernel/tab/tabcols.c
@@ -13,9 +13,10 @@
string defv[] = { /* DEFAULT INPUT PARAMETERS */
"in=???\n input file name(s)",
"select=all\n columns to select",
- "colsep=SP\n Column separator (SP,TAB,NL)",
+ "colsep=SP\n Column separator for output (SP,TAB,NL)",
+ "colsepin=\n Optional enforced character to separate columns in input",
"out=-\n output file name",
- "VERSION=2.3\n 30-mar-2024 PJT",
+ "VERSION=2.5\n 10-sep-2024 PJT",
NULL
};
@@ -35,7 +36,8 @@ int ninput; /* number of input files */
int keep[MAX_COL+1]; /* column numbers to keep */
int nkeep; /* actual number of skip columns */
int maxcol = 0; /* largest column number */
-char colsep; /* character to separate columns */
+string colsepin; /* character(s) to separate columns on input */
+char colsep; /* character to separate columns on output */
bool Qall;
local void setparams(void);
@@ -104,6 +106,10 @@ local void setparams(void)
case '/': colsep = separator[0]; break;
default: error("Illegal separator: s t n r : ,");
}
+ dprintf(1,"colsep:'%c'\n", colsep);
+
+ colsepin = getparam("colsepin");
+ dprintf(1,"colsepin:'%s'\n", colsepin);
}
@@ -112,8 +118,13 @@ local void convert(stream instr, stream outstr)
char *line;
int i, nlines, noutv;
string *outv; /* pointer to vector of strings to write */
- char *seps=", |\t"; /* column separators */
-
+ char *seps;
+
+ if (strlen(colsepin) == 0)
+ seps = strdup(", |\t");
+ else
+ seps = strdup(colsepin);
+
nlines=0; /* count lines read so far */
for (;;) {
diff --git a/src/kernel/tab/tabpeak.c b/src/kernel/tab/tabpeak.c
index 4ff977323..ecc3e4832 100644
--- a/src/kernel/tab/tabpeak.c
+++ b/src/kernel/tab/tabpeak.c
@@ -38,7 +38,7 @@ string defv[] = {
"npeak=0\n extract the Nth peak (N>0)",
"epeak=1\n expand around the peak by this factor",
"nmax=100000\n max size if a pipe",
- "VERSION=0.8\n 25-may-2023 PJT",
+ "VERSION=0.8a\n 25-dec-2024 PJT",
NULL
};
@@ -163,12 +163,20 @@ local void singles(void)
n=0;
for (i=edge; i delta && dcol[i+1] > delta) {
+#if 0
+ // try to reject astrophysical signals
+ // also needs: edge = edge + pmin/2 ?
+ if (pmin > 1) { /* peak needs to be a single channel peak */
+ if (dcol[i-1] > delta) continue;
+ if (dcol[i+2] > delta) continue;
+ }
+#endif
printf("%f %f %f %f %f\n",xcol[i], ycol[i],
ycol[i] - 0.5*(ycol[i-1]+ycol[i+1]),dcol[i],dcol[i+1]);
n++;
}
}
- dprintf(1,"Found %d singles\n",n);
+ dprintf(1,"Found %d singles; pmin=%d\n",n,pmin);
}
local void peak_data(void)
diff --git a/src/nbody/evolve/flowcode/Makefile b/src/nbody/evolve/flowcode/Makefile
index 4cc8bff77..53496a526 100644
--- a/src/nbody/evolve/flowcode/Makefile
+++ b/src/nbody/evolve/flowcode/Makefile
@@ -72,6 +72,8 @@ clean:
tar:
tar cvf flowcode.tar $(TARFILES)
+all: $(BINFILES)
+
# $(CC) $(CFLAGS) $(LOCAL_INC) -o $* $*.o $(NEMO_LIBS) $(LOCAL_LIB) $(FORLIBS) $(EL)
flowcode: $(OBJFILES)
diff --git a/src/nbody/evolve/multicode/quadcode.c b/src/nbody/evolve/multicode/quadcode.c
index bfc632a77..26a5651d7 100644
--- a/src/nbody/evolve/multicode/quadcode.c
+++ b/src/nbody/evolve/multicode/quadcode.c
@@ -24,13 +24,12 @@ string defv[] = { /* DEFAULT INPUT PARAMETERS */
"minor_freqout=32.0\n minor data-output frequency",
"options=\n misc options",
"headline=\n random mumble for humans",
- "VERSION=1.4\n 6-jan-2022 PJT",
+ "VERSION=1.4a\n 20-oct-2024 PJT",
NULL,
};
string usage = "Global quadrupole N-body integrator";
-string cvsid="$Id$";
local void force(Body *, int , real);
diff --git a/src/nbody/evolve/multicode/quaddefs.h b/src/nbody/evolve/multicode/quaddefs.h
index aa31b8a42..5cd62195a 100644
--- a/src/nbody/evolve/multicode/quaddefs.h
+++ b/src/nbody/evolve/multicode/quaddefs.h
@@ -67,7 +67,7 @@ extern void pcstep(body *btab, int nb, real *tptr, force_proc force, real dt);
extern void moveaccel(body *btab, int nb);
/* quadforce.c */
-extern int quadforce(body *btab, int nb, real eps1, real eps2);
+extern void quadforce(body *btab, int nb, real eps1, real eps2);
/* quadinter.c */
-extern int quadinter(Body *btab, int nb, real eps1, real eps2);
+extern void quadinter(Body *btab, int nb, real eps1, real eps2);
diff --git a/src/nbody/evolve/multicode/quadforce.c b/src/nbody/evolve/multicode/quadforce.c
index f092060a9..e9fc76277 100644
--- a/src/nbody/evolve/multicode/quadforce.c
+++ b/src/nbody/evolve/multicode/quadforce.c
@@ -6,6 +6,7 @@
* 16-mar-90 PJT made GCC happy
* 20-may-94 pjt remove allocate() decl. into headers
* 25-dec-02 pjt better ANSI C
+ * 20-oct-24 pjt gcc-14 function pointer syntax fix
*/
#include "quaddefs.h"
@@ -20,13 +21,13 @@ typedef struct {
local void init_quad_field(shadow *, int, real);
local void int_quad_force(shadow *, int);
local void ext_quad_force(shadow *, int);
-local int rankrad(shadowptr, shadowptr);
+local int rankrad(const void *, const void *);
/*
* QUADFORCE: compute the force-field of a spheroidal distribution.
*/
-quadforce(Body *btab, int nb, real eps1, real eps2)
+void quadforce(Body *btab, int nb, real eps1, real eps2)
{
shadowptr shad;
Body *b;
@@ -52,8 +53,11 @@ quadforce(Body *btab, int nb, real eps1, real eps2)
* RANKRAD: ranking function for quicksort.
*/
-local int rankrad(shadowptr sp1, shadowptr sp2)
+//local int rankrad(shadowptr sp1, shadowptr sp2)
+local int rankrad(const void *p1, const void *p2)
{
+ shadowptr sp1 = (shadowptr) p1;
+ shadowptr sp2 = (shadowptr) p2;
return (sp1->rad0 < sp2->rad0 ? -1 : 1); /* compare body radii */
}
diff --git a/src/nbody/evolve/multicode/quadforce_main.c b/src/nbody/evolve/multicode/quadforce_main.c
index e451ba9c5..c88c01b27 100644
--- a/src/nbody/evolve/multicode/quadforce_main.c
+++ b/src/nbody/evolve/multicode/quadforce_main.c
@@ -22,13 +22,15 @@ string defv[] = {
"quad=\n output file with field tables",
"eps_r=0.05\n radial softening parameter",
"eps_t=0.07\n tangential softening parameter",
- "VERSION=1.4\n 6-jan-2022 PJT",
+ "VERSION=1.4a\n 20-oct-2024 PJT",
NULL,
};
string usage="quadrupole-order force calculation of an N-body system.";
-nemo_main()
+void put_quadfield(stream quadstr, real tsnap, real eps1, real eps2);
+
+void nemo_main()
{
stream instr, outstr, quadstr;
Body *btab = NULL, *bp;
@@ -64,9 +66,7 @@ nemo_main()
}
}
-put_quadfield(quadstr, tsnap, eps1, eps2)
-stream quadstr;
-real tsnap, eps1, eps2;
+void put_quadfield(stream quadstr, real tsnap, real eps1, real eps2)
{
put_set(quadstr, "QuadField");
put_data(quadstr, "Time", RealType, &tsnap, 0);
diff --git a/src/nbody/evolve/multicode/quadinter.c b/src/nbody/evolve/multicode/quadinter.c
index 7213f86f1..f578bbdc1 100644
--- a/src/nbody/evolve/multicode/quadinter.c
+++ b/src/nbody/evolve/multicode/quadinter.c
@@ -9,14 +9,14 @@
#include "quaddefs.h"
-local inter_field(real r);
-local force_eval(Body *b, real eps1, real eps2);
+local void inter_field(real r);
+local void force_eval(Body *b, real eps1, real eps2);
/*
* QUADINTER: interpolate field and evaluate force on array of particles.
*/
-quadinter(Body *btab, int nb, real eps1, real eps2)
+void quadinter(Body *btab, int nb, real eps1, real eps2)
{
Body *b;
@@ -59,7 +59,7 @@ local matrix Q22, P22;
* INTER_FIELD: interpolate field tables to given radius.
*/
-local inter_field(real r)
+local void inter_field(real r)
{
int i, j, k;
real f;
@@ -94,7 +94,7 @@ local inter_field(real r)
* FORCE_EVAL: compute force and potential on particle.
*/
-local force_eval(Body *b, real eps1, real eps2)
+local void force_eval(Body *b, real eps1, real eps2)
{
real rsq, r1i, r2i, r2is, r2iq, q11r, rq22r, tmp;
vector q22r, p22r, tmpv;
diff --git a/src/nbody/evolve/multicode/quadinter_main.c b/src/nbody/evolve/multicode/quadinter_main.c
index b39a446b3..8caa7e5c6 100644
--- a/src/nbody/evolve/multicode/quadinter_main.c
+++ b/src/nbody/evolve/multicode/quadinter_main.c
@@ -20,7 +20,7 @@ string defv[] = {
"quad=???\n input file with field tables",
"in=???\n input file with N-body snapshot",
"out=\n output file with forces and potentials",
- "VERSION=1.4\n 12-nov-91 PJT",
+ "VERSION=1.4a\n 20-oct-2024 PJT",
NULL,
};
diff --git a/src/nbody/image/snapccd.c b/src/nbody/image/snapccd.c
index 4946c836a..c2a006648 100644
--- a/src/nbody/image/snapccd.c
+++ b/src/nbody/image/snapccd.c
@@ -36,11 +36,11 @@ string defv[] = {
"cell=0.0625\n cellsize : 64 pixels if size=4",
"vrange=-infinity:infinity\n range in velocity space",
"moment=0\n velocity moment to weigh with",
- "VERSION=4.7\n 27-jan-2021 PJT",
+ "VERSION=4.7a\n 9-oct-2023 PJT",
NULL,
};
-string usage = "simple conversion of snapshot to image";
+string usage = "simple mass-based conversion of snapshot to image";
#define EPS 0.00001
#ifndef HUGE
@@ -55,10 +55,9 @@ local int nobj;
local real tnow;
local real *mass=NULL;
local real *phase=NULL;
-local string headline;
/* IMAGE INTERFACE */
local imageptr iptr=NULL; /* will be allocated dynamically */
-local int nx,ny,nsize; /* map-size */
+local int nsize; /* map-size */
local real origin[2]; /* of map in sky coordinates */
local real size; /* size of frame (square) */
local real cell; /* cell or pixel size (square) */
@@ -92,7 +91,7 @@ void setparams()
size = getdparam("size");
cell = getdparam("cell");
- nsize = size/cell + EPS; /* fudge factor to prevent unfavourable roundoff */
+ nsize = rint(size/cell); /* use nearest integer */
dprintf (1,"size of CCD frame = %d square\n",nsize);
moment=getiparam("moment");
@@ -204,7 +203,7 @@ void bin_data()
{
real xsky, ysky, vrad;
real m_min, m_max, brightness, inv_surden, total;
- int i, k, ix, iy, nx, ny, cnt, noutside, noutvel, ndata;
+ int i, k, ix, iy, nx, ny, noutside, noutvel, ndata;
real *pptr;
/* (re)initialize CCD and some other local variables */
@@ -289,4 +288,4 @@ void nemo_main ()
strclose(instr);
strclose(outstr);
-}
\ No newline at end of file
+}
diff --git a/src/nbody/image/snapgrid.c b/src/nbody/image/snapgrid.c
index 25bbb59bf..b2d471f9f 100644
--- a/src/nbody/image/snapgrid.c
+++ b/src/nbody/image/snapgrid.c
@@ -77,7 +77,7 @@ string defv[] = { /* keywords/default values/help */
"stack=f\n Stack all selected snapshots?",
"integrate=f\n Sum or Integrate along 'dvar'?",
"proj=\n Sky projection (SIN, TAN, ARC, NCP, GLS, CAR, MER, AIT)",
- "VERSION=6.1a\n 9-feb-2024 PJT",
+ "VERSION=6.2\n 22-dec-2024 PJT",
NULL,
};
@@ -87,7 +87,6 @@ string usage="grid a snapshot into a 2/3D image-cube";
#ifndef HUGE
#define HUGE 1.0e20 /* don't use INF, ccdfits writes bad headers */
#endif
-#define TIMEFUZZ 0.000001
#define CUTOFF 4.0 /* cutoff of gaussian in terms of sigma */
#define MAXVAR 16 /* max evar's */
@@ -461,8 +460,7 @@ typedef struct point {
local Point **map =NULL;
-local int pcomp(Point **a, Point **b);
-//local int pcomp(void *, void *);
+local int pcomp(const void *, const void *);
void bin_data(int ivar)
{
@@ -724,9 +722,11 @@ void los_data(void)
}
-int pcomp(Point **a, Point **b)
+int pcomp(const void *app, const void *bpp)
{
- return (*a)->depth > (*b)->depth;
+ Point *a = (Point *) app;
+ Point *b = (Point *) bpp;
+ return (a)->depth > (b)->depth;
}
void free_snap()
diff --git a/src/nbody/image/snapgridsmooth.c b/src/nbody/image/snapgridsmooth.c
index 6d4dca0dc..e5cdde287 100644
--- a/src/nbody/image/snapgridsmooth.c
+++ b/src/nbody/image/snapgridsmooth.c
@@ -38,7 +38,7 @@ string defv[] = {
"stack=f\n Stack all selected snapshots?",
"periodic=f\n Periodic boundary conditions for smoothing?",
"normalize=t\n Normalize smoothing to conserve mass (evar)",
- "VERSION=0.5\n 19-mar-2022 PJT",
+ "VERSION=0.6\n 22-dec-2024 PJT",
NULL,
};
@@ -46,7 +46,6 @@ string usage="grid a snapshot into a 3D image-cube with adaptive smoothing";
#define HUGE 1.0e20 /* don't use INF, ccdfits writes bad headers */
-#define TIMEFUZZ 0.000001
#define CUTOFF 4.0 /* cutoff of gaussian in terms of sigma */
#define MAXVAR 16 /* max evar's */
diff --git a/src/nbody/image/snapifu.c b/src/nbody/image/snapifu.c
index 97d67fb8b..cfc67160d 100644
--- a/src/nbody/image/snapifu.c
+++ b/src/nbody/image/snapifu.c
@@ -38,16 +38,14 @@ string defv[] = { /* keywords/default values/help */
"moment=0\n moment in zvar (-2,-1,0,1,2...)",
"mean=f\n mean (moment=0) or sum per cell",
"stack=f\n Stack all selected snapshots?",
- "VERSION=1.0\n 8-apr-09 PJT",
+ "VERSION=1.1\n 22-dec-2024 PJT",
NULL,
};
string usage="take spectra from a snapshot at a set of specified grid points";
-string cvsid="$Id$";
#define HUGE 1.0e20 /* don't use INF, ccdfits writes bad headers */
-#define TIMEFUZZ 0.000001
#define CUTOFF 4.0 /* cutoff of gaussian in terms of sigma */
#define MAXVAR 16 /* max evar's */
#define MAXPOINTS 1024 /* max # gridpoints */
diff --git a/src/nbody/image/snapmap.c b/src/nbody/image/snapmap.c
index 8f3cda114..aee7f54dc 100644
--- a/src/nbody/image/snapmap.c
+++ b/src/nbody/image/snapmap.c
@@ -46,15 +46,13 @@ string defv[] = { /* keywords/default values/help */
"stack=f\n Stack all selected snapshots?",
"proj=\n Sky projection (SIN,TAN,ARC,NCP,GLS,MER,AIT)",
"emax=10\n Maximum exp smoothing parameter",
- "VERSION=2.1\n 23-sep-2021 PJT",
+ "VERSION=2.2\n 22-dec-2024 PJT",
NULL,
};
string usage="grid a snapshot into a 2D map";
-string cvsid="$Id$";
-#define TIMEFUZZ 0.000001
#define MAXVAR 16 /* max evar's */
#define SMODE_GAUSS 1
diff --git a/src/nbody/init/diskspectrum.c b/src/nbody/init/diskspectrum.c
new file mode 100644
index 000000000..8543775c9
--- /dev/null
+++ b/src/nbody/init/diskspectrum.c
@@ -0,0 +1,351 @@
+/*
+ * DISKSPECTRUM.C: combining mkdisk + snapmass + snaprotate + snapgriof
+ * to create a faster method for edge_gbt.sh
+ *
+ * 24-nov-2024 cloned off mkdisk
+ *
+ */
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+string defv[] = {
+ "out=???\n Output file name (spectrum)",
+ "nbody=2048\n Number of disk particles",
+
+ "potname=plummer\n Name of potential(5)",
+ "potpars=\n Parameters to potential(5); omega needed but not used",
+ "potfile=\n Optional data file with potential(5)",
+
+ "rmin=0\n Inner disk radius",
+ "rmax=1\n Outer cutoff radius",
+ "model=0\n Mass model (0=CONST 1=EXP, 2=PLEC... 3=AMT)",
+ "mpars=\n Parameters for the mass model",
+
+ "frac=0\n Relative vel.disp w.r.t. local rotation speed",
+ "seed=0\n Usual random number seed",
+
+ "angle=f\n Regular angular distribution?",
+ "vrad=0\n radial velocity",
+ "energy=f\n preserve energy if random motions added?",
+ "abs=f\n Use absolute vel.disp instead of fractional?",
+ "z0=0,0\n ** Vertical scaleheight for density; use 2nd one for velocity dropoff if needed",
+ "vloss=-1\n ** Fractional loss of orbital speed at the scaleheight (<1 => Burkert)",
+
+ "inc=30\n Inclination to observe the disk at",
+
+ "vbeam=5\n FWHM of smoothing beam in velocity",
+ "vrange=400\n velocity gridding will be -vrange:vrange",
+ "nvel=200\n number of spectral pixels",
+
+ "headline=\n Text headline for output",
+ "VERSION=0.2\n 27-nov-2024 PJT",
+ NULL,
+};
+
+string usage="global spectrum of a rotating thin disk";
+
+local real rmin, rmax, mass;
+local bool Qangle;
+local bool Qenergy;
+local bool Qabs;
+
+local int ndisk;
+local real frac[NDIM], vrad;
+local Body *disk;
+
+local proc potential;
+
+local real z0_d; /* the old z0 */
+local real z0_v; /* dropoff in velocity */
+local real vloss;
+local bool Qrandom;
+
+local real inc, sininc;
+local real vbeam;
+local real vrange;
+local int nvel;
+
+local int mmodel;
+local real mpars[16];
+// 1: EXP: re
+// 2: PLEC: rd,n
+// 3: AMT: rd,Rm
+
+/* old style */
+// #define OLD_BURKERT 1
+
+
+extern double xrandom(double,double), grandom(double,double);
+
+local real took(real);
+local void testdisk(void);
+local void spectrum(void);
+
+local real mysech2(real z)
+{
+ real y = exp(z);
+ return sqr(2.0/(y+1.0/y));
+}
+
+local real pick_z(real z0)
+{
+ real z = frandom(-6.0,6.0,mysech2) * z0;
+ return z;
+}
+
+local real pick_dv(real r, real z, real z0)
+{
+#ifdef OLD_BURKERT
+ real dv = 1 - (1 + (z/z0) * tanh(z/z0))*(z0/r);
+#else
+ //real dv = tanh(z/z0);
+ //dv = sqrt(1 - ABS(dv));
+ real dv = ABS(z/z0);
+ //dv = sqrt(exp(-dv));
+ dv = exp(-dv);
+#endif
+ dprintf(1,"PJT %g %g %g\n",r,z,dv);
+ return dv;
+}
+
+void nemo_main()
+{
+ int nfrac, seed, nz;
+ real z0[2];
+
+ warning("This program is under development");
+
+ potential = get_potential(getparam("potname"),
+ getparam("potpars"), getparam("potfile"));
+ rmin = getdparam("rmin");
+ rmax = getdparam("rmax");
+ vrad = getdparam("vrad");
+ ndisk = getiparam("nbody");
+
+ /* z0= is now split in a density and velocity */
+ nz = nemoinpr(getparam("z0"),z0,2);
+ if (nz == 0)
+ z0[0] = z0[1] = 0.0;
+ else if (nz == 1)
+ z0[1] = 0.0;
+ z0_d = z0[0];
+ z0_v = z0[1];
+
+
+ vloss = getrparam("vloss");
+ if (z0_d > 0) {
+ if (vloss < 0)
+#ifdef OLD_BURKERT
+ dprintf(0,"Burkert et al 2010 streaming(z) model\n");
+#else
+ dprintf(0,"new style tanh(z) streaming model\n");
+#endif
+ else
+ dprintf(0,"Toy streaming(z) model with vloss=%g\n",vloss);
+ }
+ nfrac = nemoinpr(getparam("frac"),frac,NDIM);
+ switch (nfrac) {
+ case 1:
+ frac[1] = frac[0];
+ frac[2] = 0.0;
+ break;
+ case 2:
+ frac[2] = 0.0;
+ break;
+ case 3:
+ break;
+ default:
+ error("%d: bad parsing frac=%s",nfrac,getparam("frac"));
+ }
+ dprintf(1,"frac: %g %g %g\n",frac[0],frac[1],frac[2]);
+ Qrandom = (frac[0]>0 || frac[1]>0 || frac[2]>0);
+ if (!Qrandom)
+ dprintf(0,"No random motions\n");
+
+ mass = 1.0 / ndisk;
+ seed = init_xrandom(getparam("seed"));
+ dprintf(1,"Seed=%d\n",seed);
+ Qangle = getbparam("angle");
+ Qenergy = getbparam("energy");
+ Qabs = getbparam("abs");
+ vbeam = getdparam("vbeam");
+ vrange = getdparam("vrange");
+ nvel = getiparam("nvel");
+ inc = getdparam("inc");
+ sininc = sin(inc/DR2D);
+
+ testdisk();
+ spectrum();
+}
+
+/*
+ * TESTDISK: use forces due to a potential to make a uniform
+ * density test disk.
+ */
+
+void testdisk(void)
+{
+ Body *dp;
+ real rmin2, rmax2, r_i, theta_i, vcir_i;
+ real dv_r, dv_t, sint, cost, vrandom;
+ real sigma_r, sigma_t, sigma_z;
+ vector acc_i;
+ int i, ncirc, ndim=NDIM;
+ double pos_d[NDIM], acc_d[NDIM], pot_d, time_d = 0.0;
+
+ disk = (Body *) allocate(ndisk * sizeof(Body));
+ rmin2 = rmin * rmin;
+ rmax2 = rmax * rmax;
+ theta_i = xrandom(0.0, TWO_PI);
+ for (dp=disk, i = 0, ncirc=0; i < ndisk; dp++, i++) { /* loop over all stars */
+ Mass(dp) = mass;
+ if (ndisk == 1)
+ r_i = rmin;
+ else
+ r_i = sqrt(rmin2 + i * (rmax2 - rmin2) / (ndisk - 1.0));
+ if (Qangle) {
+ theta_i += TWO_PI/ndisk;
+ } else {
+ theta_i = xrandom(0.0, TWO_PI);
+ }
+ cost = cos(theta_i);
+ sint = sin(theta_i);
+ Pos(dp)[0] = pos_d[0] = r_i * cost; /* set positions */
+ Pos(dp)[1] = pos_d[1] = r_i * sint;
+ Pos(dp)[2] = pos_d[2] = 0.0; /* it's a DISK ! */
+ (*potential)(&ndim,pos_d,acc_d,&pot_d,&time_d); /* get forces */
+ SETV(acc_i,acc_d);
+ vcir_i = sqrt(r_i * absv(acc_i)); /* v^2 / r = force */
+ /* now cheat and rotate slower away from the plane */
+ if (z0_d > 0) {
+ if (vloss >= 0.0) {
+ // toy model (early sep 2017)
+ Pos(dp)[2] = grandom(0.0, 0.5 * z0_d);
+ vcir_i *= (1-vloss*ABS(Pos(dp)[2]/z0_d));
+ if (vcir_i < 0) vcir_i = 0.0; /* added dec 2017 */
+ } else {
+ // Burkert et al. 2010 model, doesn't need vloss= anymore, triggered when vloss < 0
+ Pos(dp)[2] = pick_z(z0_d);
+ vcir_i = vcir_i*vcir_i;
+ vcir_i *= pick_dv(r_i,Pos(dp)[2],z0_v);
+ if (vcir_i > 0)
+ vcir_i = sqrt(vcir_i);
+ else
+ vcir_i = 0.0;
+ }
+ }
+
+#if 1
+ if (Qabs) {
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]);
+ sigma_t = grandom(0.0,frac[1]);
+ sigma_z = grandom(0.0,frac[2]);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
+ //Vel(dp)[0] = -vcir_i * sint ;
+ Vel(dp)[1] = vcir_i * cost ;
+ //Vel(dp)[0] += cost*sigma_r - sint*sigma_t; /* add dispersions */
+ Vel(dp)[1] += sint*sigma_r + cost*sigma_t;
+ } else {
+ do { /* iterate, if needed, to get vrandom */
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]*vcir_i);
+ sigma_t = grandom(0.0,frac[1]*vcir_i);
+ sigma_z = grandom(0.0,frac[2]*vcir_i);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
+ dv_t = sigma_t;
+ dv_r = sigma_r * took(r_i) ;
+ vrandom = sqrt(dv_t*dv_t + dv_r*dv_r);
+ if (vrandom > vcir_i) ncirc++;
+ } while (Qenergy && vrandom > vcir_i);
+ vcir_i = sqrt((vcir_i-vrandom)*(vcir_i+vrandom));
+ dv_r += vrad;
+ //Vel(dp)[0] = -vcir_i * sint ;
+ Vel(dp)[1] = vcir_i * cost ;
+ //Vel(dp)[0] += cost*dv_r - sint*dv_t; /* add dispersions */
+ Vel(dp)[1] += sint*dv_r + cost*dv_t;
+ }
+#else
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]*vcir_i);
+ sigma_t = grandom(0.0,frac[1]*vcir_i);
+ sigma_z = grandom(0.0,frac[2]*vcir_i);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
+ dv_t = sigma_t;
+ dv_r = sigma_r * took(r_i) ;
+
+ /* Qenergy only uses radial motion: thus preserving the
+ * guiding center for epicycles ?? (Olling 2003)
+ */
+
+ if (Qenergy)
+ vcir_i = sqrt((vcir_i-dv_r)*(vcir_i+dv_r));
+ //Vel(dp)[0] = -vcir_i * sint ;
+ Vel(dp)[1] = vcir_i * cost ;
+ //Vel(dp)[0] += cost*dv_r;
+ Vel(dp)[1] += sint*dv_r;
+ if (!Qenergy) {
+ //Vel(dp)[0] += -sint*dv_t;
+ Vel(dp)[1] += cost*dv_t;
+ }
+
+#endif
+ Vel(dp)[2] = sigma_z;
+ Vel(dp)[1] *= sininc;
+ }
+ if (ncirc) dprintf(0,"Had to respin random %d times\n",ncirc);
+}
+
+
+/*
+ * 2*omega/kappa; from spline interpolation.....
+ */
+
+real took(real r)
+{
+ return 1.0;
+}
+
+
+void spectrum(void)
+{
+ real mtot = 0.0;
+ real vrad, vrad_max = -1.0;
+ Body *dp;
+ int i, iv;
+ Grid g;
+ real *spectrum = (real *) allocate(nvel*sizeof(real));
+ // nvel needs to be even for symmetric profiles
+
+ inil_grid(&g, nvel, 0, vrange);
+
+ for (dp=disk, i = 0; i < ndisk; dp++, i++) { /* loop over all stars */
+ vrad = ABS(Vel(dp)[1]);
+ if (vrad > vrange) continue;
+ if (vrad > vrad_max) vrad_max = vrad;
+ mtot += Mass(dp);
+ iv = index_grid(&g, vrad);
+ spectrum[iv] += Mass(dp);
+ }
+ printf("# Total mass: %g ndisk=%d vrad_max: %g\n", mtot, ndisk, vrad_max);
+ for (i=0; i 0) warning("no smoothing done yet");
+}
diff --git a/src/nbody/init/mkdisk.c b/src/nbody/init/mkdisk.c
index 02c795a48..ce9be7aaa 100644
--- a/src/nbody/init/mkdisk.c
+++ b/src/nbody/init/mkdisk.c
@@ -56,7 +56,7 @@ string defv[] = {
"z0=0,0\n Vertical scaleheight for density; use 2nd one for velocity dropoff if needed",
"vloss=-1\n Fractional loss of orbital speed at the scaleheight (<1 => Burkert)",
"headline=\n Text headline for output",
- "VERSION=4.9i\n 19-jan-2018 PJT",
+ "VERSION=4.9j\n 23-nov-2024 PJT",
NULL,
};
@@ -77,6 +77,7 @@ local proc potential;
local real z0_d; /* the old z0 */
local real z0_v; /* dropoff in velocity */
local real vloss;
+local bool Qrandom;
/* old style */
// #define OLD_BURKERT 1
@@ -166,6 +167,9 @@ void nemo_main()
error("%d: bad parsing frac=%s",nfrac,getparam("frac"));
}
dprintf(1,"frac: %g %g %g\n",frac[0],frac[1],frac[2]);
+ Qrandom = (frac[0]>0 || frac[1]>0 || frac[2]>0);
+ if (!Qrandom)
+ dprintf(0,"No random motions\n");
mass = getdparam("mass") / ndisk;
if (mass==0.0) {
@@ -266,18 +270,24 @@ void testdisk(void)
#if 1
if (Qabs) {
- sigma_r = grandom(0.0,frac[0]);
- sigma_t = grandom(0.0,frac[1]);
- sigma_z = grandom(0.0,frac[2]);
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]);
+ sigma_t = grandom(0.0,frac[1]);
+ sigma_z = grandom(0.0,frac[2]);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
Vel(dp)[0] = -vcir_i * sint * jz_sign;
Vel(dp)[1] = vcir_i * cost * jz_sign;
Vel(dp)[0] += cost*sigma_r - sint*sigma_t; /* add dispersions */
Vel(dp)[1] += sint*sigma_r + cost*sigma_t;
} else {
do { /* iterate, if needed, to get vrandom */
- sigma_r = grandom(0.0,frac[0]*vcir_i);
- sigma_t = grandom(0.0,frac[1]*vcir_i);
- sigma_z = grandom(0.0,frac[2]*vcir_i);
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]*vcir_i);
+ sigma_t = grandom(0.0,frac[1]*vcir_i);
+ sigma_z = grandom(0.0,frac[2]*vcir_i);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
dv_t = sigma_t;
dv_r = sigma_r * took(r_i) ;
vrandom = sqrt(dv_t*dv_t + dv_r*dv_r);
@@ -291,9 +301,12 @@ void testdisk(void)
Vel(dp)[1] += sint*dv_r + cost*dv_t;
}
#else
- sigma_r = grandom(0.0,frac[0]*vcir_i);
- sigma_t = grandom(0.0,frac[1]*vcir_i);
- sigma_z = grandom(0.0,frac[2]*vcir_i);
+ if (Qrandom) {
+ sigma_r = grandom(0.0,frac[0]*vcir_i);
+ sigma_t = grandom(0.0,frac[1]*vcir_i);
+ sigma_z = grandom(0.0,frac[2]*vcir_i);
+ } else
+ sigma_r = sigma_t = sigma_z = 0.0;
dv_t = sigma_t;
dv_r = sigma_r * took(r_i) ;
diff --git a/src/nbody/init/mkflowdisk.c b/src/nbody/init/mkflowdisk.c
index 77dc58468..3785070a3 100644
--- a/src/nbody/init/mkflowdisk.c
+++ b/src/nbody/init/mkflowdisk.c
@@ -46,7 +46,7 @@ string defv[] = {
"test=f\n test shape of spiral (all particles at 0 phase offset)",
"constant=f\n force vt constant in rotating frame ?",
"headline=\n text headline for output ",
- "VERSION=1.5\n 1-jan-04 PJT",
+ "VERSION=1.5a\n 20-oct-2024 PJT",
NULL,
};
@@ -71,6 +71,11 @@ local real theta[361], dens[361];
local proc potential;
+void writegalaxy(stream outstr);
+void setdensity(void);
+double density(double t);
+void testdisk(int n);
+
void nemo_main()
{
@@ -123,7 +128,7 @@ void nemo_main()
* WRITEGALAXY: write galaxy model to output.
*/
-writegalaxy(stream outstr)
+void writegalaxy(stream outstr)
{
real tsnap = 0.0;
int bits = MassBit | PhaseSpaceBit | TimeBit;
@@ -132,7 +137,7 @@ writegalaxy(stream outstr)
put_snap(outstr, &disk, &ndisk, &tsnap, &bits);
}
-setdensity(void)
+void setdensity(void)
{
int i, ndim=NDIM;
double pos_d[NDIM], vel_d[NDIM], den_d, time_d = 0.0;
@@ -197,10 +202,10 @@ double density(double t)
* density test disk.
*/
-testdisk(int n)
+void testdisk(int n)
{
Body *dp;
- real rmin2, rmax2, r_i, theta_i, mass_i, rring;
+ real rmin2, rmax2, r_i, theta_i, mass_i;
real cost, sint;
int i, ndim=NDIM;
double pos_d[NDIM], vel_d[NDIM], pot_d, time_d = 0.0;
@@ -209,7 +214,6 @@ testdisk(int n)
if (disk == NULL) disk = (Body *) allocate(ndisk * sizeof(Body));
rmin2 = rmin * rmin;
rmax2 = rmax * rmax;
- rring = 0.5*(rmin+rmax);
mass_i = 1.0 / (real)ndisk;
for (dp=disk, i = 0; i < ndisk; dp++, i++) {
diff --git a/src/nbody/init/mkplummer.c b/src/nbody/init/mkplummer.c
index 1487e9de9..6bc773668 100644
--- a/src/nbody/init/mkplummer.c
+++ b/src/nbody/init/mkplummer.c
@@ -66,7 +66,7 @@ string defv[] = { /* DEFAULT INPUT PARAMETERS */
"headline=\n Verbiage for output",
"nmodel=1\n number of models to produce",
"mode=1\n 0=no data, 1=data, no analysis 2=data, analysis",
- "VERSION=3.0b\n 2-sep-2021 PJT",
+ "VERSION=3.0c\n 6-jun-2024 PJT",
NULL,
};
@@ -74,7 +74,7 @@ string usage="construct a Plummer model";
/*---------------------------------------------------------------------------
* main -- a tool to make a Plummer model, by invoking mkplummer()
- * note: only for 3 dimensions and Carthesian coordinates
+ * note: only for 3 dimensions and Cartesian coordinates
*---------------------------------------------------------------------------
*/
void nemo_main(void)
@@ -136,8 +136,9 @@ void nemo_main(void)
init_xrandom(sseed);
}
btab[i] = mkplummer(nbody, mlow, mfrac, rfrac, seed, snap_time, zerocm, scale,
- quiet,mrange,mfunc);
+ quiet,mrange,mfunc);
}
+
if (mode > 0) {
bits = (MassBit | PhaseSpaceBit | TimeBit);
sprintf(hisline,"init_xrandom: seed used %d",seed);
diff --git a/src/nbody/io/testfio.c b/src/nbody/io/testfio.c
index a150fad20..a42cd6929 100644
--- a/src/nbody/io/testfio.c
+++ b/src/nbody/io/testfio.c
@@ -21,11 +21,10 @@
string defv[] = {
"in=???\n input file name",
"out=???\n output file name",
- "VERSION=1.0", 15-sep-90 PJT",
+ "VERSION=1.1\n 22-dec-2024 PJT",
NULL,
};
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
main(argc, argv)
int argc;
diff --git a/src/nbody/io/testio.c b/src/nbody/io/testio.c
index dbd1d7e34..a7f366ccc 100644
--- a/src/nbody/io/testio.c
+++ b/src/nbody/io/testio.c
@@ -27,11 +27,10 @@ string defv[] = {
"out=\n output file name (optional)",
"io=get_snapshot\n mode of I/O {get_snapshot, get_snap}",
"work=bodywork\n which work routine?",
- "VERSION=1.1\n 16-sep-90 PJT",
+ "VERSION=1.2\n 22-dec-2024 PJT",
NULL,
};
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
main(argc, argv)
int argc;
diff --git a/src/nbody/io_nemo/src/io_get_put.c b/src/nbody/io_nemo/src/io_get_put.c
index 0e5ff2e0c..0315d7f7c 100644
--- a/src/nbody/io_nemo/src/io_get_put.c
+++ b/src/nbody/io_nemo/src/io_get_put.c
@@ -42,7 +42,6 @@
#include "parameters.h"
-#define TIMEFUZZ 0.0000001
extern int maxbodies[];
extern int CURRENT_IO;
/* -----------------------------------------------------------------
diff --git a/src/nbody/io_nemo/src/io_get_put_f.c b/src/nbody/io_nemo/src/io_get_put_f.c
index 2851cb2eb..3e0ae1113 100644
--- a/src/nbody/io_nemo/src/io_get_put_f.c
+++ b/src/nbody/io_nemo/src/io_get_put_f.c
@@ -37,7 +37,6 @@
#include "parameters.h"
-#define TIMEFUZZ 0.0000001
extern int maxbodies[];
extern int CURRENT_IO;
/* ----------------------------------------------------------------
diff --git a/src/nbody/io_nemo/test_src/snapmask2.c b/src/nbody/io_nemo/test_src/snapmask2.c
index 2cfac42f8..c8a94090e 100644
--- a/src/nbody/io_nemo/test_src/snapmask2.c
+++ b/src/nbody/io_nemo/test_src/snapmask2.c
@@ -14,6 +14,7 @@
|* V2.0 : imported in teuben's nemo CVS tree
|* V2.1 : PosTag VelTag and KeyTag added
|* V2.1a : note that within() uses real, not double
+|* V2.2 : TIMEFUZZ now from stdinc.h
\* -------------------------------------------------------------- */
@@ -40,13 +41,12 @@ string defv[] = { /* DEFAULT INPUT PARAMETERS */
"select=all\n selected particles",
"times=all\n selected times",
"step=1\n copy only step by step snapshots",
- "VERSION=2.1\n 09-Jul-02 JCL",
+ "VERSION=2.2\n 22-dec-2024 PJT",
NULL,
};
string usage = "mask out particles while copying snapshot data";
-#define TIMEFUZZ 0.00001
/* -------------------------------------------------------------- *\
|* save_parameter :
|* Save ParametersTag
diff --git a/src/nbody/reduc/Testfile b/src/nbody/reduc/Testfile
index 619d99f95..c5d6d9dc9 100644
--- a/src/nbody/reduc/Testfile
+++ b/src/nbody/reduc/Testfile
@@ -1,5 +1,5 @@
DIR = src/nbody/reduc
-BIN = snapplot snapplot3 snapdiagplot snapplotv snapmradii radprof real snapfit snapprint
+BIN = snapplot snapplot3 snapdiagplot snapplotv snapmradii snapmstat radprof real snapfit snapprint snapcmp
NEED = $(BIN) hackcode1 mkplummer tabplot snapfour snapgrid snaprotate
help:
@@ -10,7 +10,7 @@ need:
clean:
@echo Cleaning $(DIR)
- @rm -f snap.in cube.in hack.out hack2.out
+ @rm -f snap.in snap2.in cube.in hack.out hack2.out
NBODY = 10
@@ -21,6 +21,11 @@ snap.in:
$(EXEC) mkplummer snap.in $(NBODY) seed=123
@bsf snap.in '0.0140845 0.896875 -4.6523 4.80925 71'
+snap2.in:
+ @echo Creating $@
+ $(EXEC) mkplummer snap2.in $(NBODY) seed=1234
+ @bsf snap2.in '0.0140845 0.536539 -1.63453 1.89357 71'
+
cube.in: snap.in
@echo Creating $@
$(EXEC) snaprotate snap.in - 30,45,20 zyz | $(EXEC) snapgrid - cube.in nx=32 ny=32 nz=32 zrange=-2:2
@@ -54,6 +59,18 @@ snapdiagplot: hack.out
snapmradii: hack2.out
@echo Running $@
$(EXEC) snapmradii hack2.out | tabplot - 1 2:9 line=1,1 ; nemo.coverage snapmradii.c
+ $(EXEC) snapmradii snap.in
+ @echo '0 0.276882 0.43351 0.5547 0.715477 0.789928 0.790888 1.05668 1.35295 4.76802'
+
+snapmstat: snap.in
+ @echo Running $@
+ $(EXEC) snapmstat snap.in ; nemo.coverage snapmradii.c
+
+
+snapcmp: snap.in snap2.in
+ @echo Running $@
+ $(EXEC) snapcmp snap.in snap2.in ; nemo.coverage snapcmp.c
+ @echo "0 0.548814 0.824836 1.96007 2.14353 5.33402 2.21031 1.5212 expected"
radprof: snap.in
@echo Running $@
diff --git a/src/nbody/reduc/snapcmp.c b/src/nbody/reduc/snapcmp.c
index ac576e8d3..e6ce383f1 100644
--- a/src/nbody/reduc/snapcmp.c
+++ b/src/nbody/reduc/snapcmp.c
@@ -55,7 +55,7 @@ string defv[] = { /* DEFAULT INPUT PARAMETERS */
"headline=\n header",
#endif
"time0=t\n Also print time=0 ?",
- "VERSION=1.7a\n 20-jan-2024 PJT",
+ "VERSION=1.7b\n 20-oct-2024 PJT",
NULL,
};
@@ -81,7 +81,7 @@ local bool Qlog, Qtime0;
local real *snapcmp(Body*, Body*, int, real);
-local int cmpreal(real*, real*);
+local int cmpreal(const void*, const void*);
local void printquart(real*, int, real);
extern rproc btrtrans(string);
@@ -136,6 +136,8 @@ void nemo_main()
if (dt > small_dt)
warning("times = %f, %f are different (%g)", tsnap1, tsnap2, dt);
}
+ dprintf(0,"time1,2: %g %g %g %g nout %d\n",time1,time2,tsnap1,tsnap2,nout);
+
if (!Qtime0 && time1==0 && time2==0) continue;
if (!Qtime0 && nout==0) continue;
if (bits1 != bits2)
@@ -151,7 +153,7 @@ void nemo_main()
snapcmpplot(result, btab1, nbody1, tsnap1);
#endif
nout++;
- }
+ } // for(;;)
}
local real *snapcmp(Body *btab1, Body *btab2, int nbody, real tsnap)
@@ -217,17 +219,16 @@ local void printquart(real result[], int nbody, real tsnap)
}
-/* should have prototype :: int (*compar)(const void *, const void *)) */
-
-local int cmpreal(real *ap, real *bp)
+local int cmpreal(const void *ap, const void *bp)
{
- return (*ap < *bp ? -1 : *ap > *bp ? 1 : 0);
+ real *a = (real *) ap;
+ real *b = (real *) bp;
+ return (*a < *b ? -1 : *a > *b ? 1 : 0);
}
-
#endif
#if HISTOGRAM
diff --git a/src/nbody/reduc/snapdiagplot.c b/src/nbody/reduc/snapdiagplot.c
index 9f52cc2b8..745c69035 100644
--- a/src/nbody/reduc/snapdiagplot.c
+++ b/src/nbody/reduc/snapdiagplot.c
@@ -27,7 +27,7 @@ string defv[] = { /* DEFAULT INPUT PARAMETERS */
"eps=0.05\n using this softening",
"formal=false\n publication-style plot",
"tab=\n If given, tabulate Time,Etot",
- "VERSION=1.5b\n 12-feb-2024 PJT",
+ "VERSION=1.6\n 22-dec-2024 PJT",
NULL,
};
@@ -38,8 +38,6 @@ local string times;
local bool exactpot;
local real eps;
local bool formal;
-
-#define TIMEFUZZ 0.00001
/*
* Arrays used to store diagnostic info.
diff --git a/src/nbody/reduc/snapfit.c b/src/nbody/reduc/snapfit.c
index 8e168085b..19892235b 100644
--- a/src/nbody/reduc/snapfit.c
+++ b/src/nbody/reduc/snapfit.c
@@ -58,7 +58,7 @@ string defv[] = {
"contour=\n Optional Output image file with chi-squared",
"iter=0\n number of more iterations after best on matrix",
"times=all\n Times selected from snapshot models",
- "VERSION=0.7d\n 15-jul-04 PJT",
+ "VERSION=0.7e\n 20-oct-2024 PJT",
NULL,
};
@@ -102,8 +102,12 @@ void snap_fit(void);
int write_snapshot(stream outstr, int nbody, body *btab, real t1, real t2, real rscale, real vscale);
void yrotate(matrix mat, real theta);
void zrotate(matrix mat, real theta);
-int eigenframe(vector frame[], matrix mat);
-int invert(vector frame[]);
+void eigenframe(vector frame[], matrix mat);
+void invert(vector frame[]);
+
+extern void eigsrt(); // (d,v,n)
+extern void matinv(); // real *, int, int,real *);
+extern void jacobi(); // (a,n,d,v,nrot)
/* -------------------------------------------------------------------------- */
@@ -381,9 +385,7 @@ real start, incr;
return a;
}
-void printvec(name, vec)
-string name;
-vector vec;
+void printvec(string name, vector vec)
{
dprintf(my_debug,"%s %10.5f %10.5f %10.5f %10.5f\n",
name, absv(vec), vec[0], vec[1], vec[2]);
@@ -545,12 +547,8 @@ void snap_fit()
rscale/w_rscale, vscale/w_vscale);
}
-
-write_snapshot( outstr, nbody, btab, t1, t2, rscale, vscale)
-stream outstr;
-int nbody;
-Body *btab;
-real t1, t2, rscale,vscale;
+//outstr, nbody, btab, t1, t2, rscale, vscale)
+int write_snapshot(stream outstr, int nbody, Body *btab, real t1, real t2, real rscale, real vscale)
{
warning("output snapshot not supported yet");
}
@@ -564,9 +562,7 @@ real t1, t2, rscale,vscale;
* around the y-axis
*/
-void yrotate(mat,theta)
-matrix mat;
-real theta;
+void yrotate(matrix mat,real theta)
{
SETMI(mat); /* unit matrix */
mat[2][2] = mat[0][0] = cos(DEG2RAD * theta);
@@ -580,9 +576,7 @@ real theta;
* around the y-axis
*/
-void zrotate(mat,theta)
-matrix mat;
-real theta;
+void zrotate(matrix mat, real theta)
{
SETMI(mat); /* unit matrix */
mat[1][1] = mat[0][0] = cos(DEG2RAD * theta);
@@ -595,9 +589,7 @@ real theta;
#if 1
-eigenframe(frame, mat) /* float version */
-vector frame[];
-matrix mat;
+void eigenframe(vector frame[], matrix mat)
{
float **q, *d, **v;
int i, j, nrot;
@@ -617,9 +609,7 @@ matrix mat;
#else
-eigenframe(frame, mat) /* double version */
-vector frame[];
-matrix mat;
+void eigenframe(vector frame[], matrix mat)
{
double **q, *d, **v;
int i, j, nrot;
@@ -638,8 +628,7 @@ matrix mat;
}
#endif
-invert(frame)
-vector frame[];
+void invert(vector frame[])
{
real mat[NDIM*NDIM], det;
diff --git a/src/nbody/reduc/snapfour.c b/src/nbody/reduc/snapfour.c
index e00fac367..3f9007369 100644
--- a/src/nbody/reduc/snapfour.c
+++ b/src/nbody/reduc/snapfour.c
@@ -31,13 +31,12 @@ string defv[] = {
"weight=1\n Weight applied to observable",
"amode=t\n Display sin/cos amps or amp/phase if possible?",
"times=all\n Snapshots to select",
- "VERSION=1.2d\n 3-apr-2020 PJT",
+ "VERSION=1.3\n 22-dec-2024 PJT",
NULL,
};
string usage = "Fourier coefficients of an N-body distribution";
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
#define MAXRAD 513
#define MAXORDER 8
diff --git a/src/nbody/reduc/snapmradii.c b/src/nbody/reduc/snapmradii.c
index 470214c42..b261a8e5c 100644
--- a/src/nbody/reduc/snapmradii.c
+++ b/src/nbody/reduc/snapmradii.c
@@ -32,19 +32,17 @@ string defv[] = {
"tab=f\n Full table of r,m(r) ? ",
"log=f\n Print radii in log10() ? ",
"sort=r\n Observerble to sort masses by",
- "VERSION=1.6\n 1-apr-2021 PJT",
+ "VERSION=1.6a\n 20-oct-2024 PJT",
NULL,
};
string usage="Langrangian mass-fraction radii of a snapshot";
-string cvsid="$Id$";
-
#define MFRACT 256
local void snapsort(Body *, int , real, rproc_body);
-local int rank_aux(Body *, Body *);
+local int rank_aux(const void *, const void *);
void nemo_main()
@@ -126,8 +124,10 @@ void snapsort(Body *btab, int nbody, real tsnap, rproc_body sortptr)
qsort(btab, nbody, sizeof(Body), rank_aux);
}
-int rank_aux(Body *a, Body *b)
+int rank_aux(const void *ap, const void *bp)
{
+ Body *a = (Body *) ap;
+ Body *b = (Body *) bp;
return (Aux(a) < Aux(b) ? -1 : Aux(a) > Aux(b) ? 1 : 0);
}
diff --git a/src/nbody/reduc/snapmstat.c b/src/nbody/reduc/snapmstat.c
index 1df871379..99af7938c 100644
--- a/src/nbody/reduc/snapmstat.c
+++ b/src/nbody/reduc/snapmstat.c
@@ -28,7 +28,7 @@ string defv[] = {
"sort=f\n Sort masses before processing?",
"species=100\n Maximum number of species",
"show=f\n Show only the select= for glnemo2",
- "VERSION=2.1\n 28-nov-2023 PJT",
+ "VERSION=2.1a\n 20-oct-24 PJT",
NULL,
};
@@ -36,7 +36,7 @@ string usage="report some statistics of the masses in a snapshot";
local void snapsort(Body *, int);
-local int rank_mass(Body *, Body *);
+local int rank_mass(const void *, const void *);
void nemo_main(void)
{
@@ -113,8 +113,11 @@ local void snapsort(Body *btab, int nbody)
qsort(btab, nbody, sizeof(Body), rank_mass);
}
-local int rank_mass(Body *a, Body *b)
+//local int rank_mass(Body *a, Body *b)
+local int rank_mass(const void *ap, const void *bp)
{
+ Body *a = (Body *) ap;
+ Body *b = (Body *) bp;
return (Mass(a) < Mass(b) ? -1 : Mass(a) > Mass(b) ? 1 : 0);
}
diff --git a/src/nbody/reduc/snapstab.c b/src/nbody/reduc/snapstab.c
index a8ee30242..f8af21cba 100644
--- a/src/nbody/reduc/snapstab.c
+++ b/src/nbody/reduc/snapstab.c
@@ -23,15 +23,14 @@
string defv[] = {
"in=???\n Input file name (snapshot)",
"times=all\n Times to select snapshots from",
- "VERSION=1.2\n 16-feb-97 PJT",
+ "VERSION=1.3\n 22-dec-2024 PJT",
NULL,
};
string usage = "report disk stability/virial";
-#define TIMEFUZZ 0.00001 /* tolerance in time comparisons */
-nemo_main()
+void nemo_main()
{
stream instr;
real inv_rscale, e_pot, e_kin, tsnap;
diff --git a/src/nbody/trans/snapcopy.c b/src/nbody/trans/snapcopy.c
index ee573b0b7..463fb1a5d 100644
--- a/src/nbody/trans/snapcopy.c
+++ b/src/nbody/trans/snapcopy.c
@@ -35,15 +35,13 @@ string defv[] = {
"precision=double\n Precision of results to store (double/single) [unused]",
"keep=all\n Items to copy in snapshot",
"i=-1\n Select one body to select (overrides select=)",
- "VERSION=2.0a\n 8-aug-2022 PJT",
+ "VERSION=2.1\n 22-dec-2024 PJT",
NULL,
};
string usage="copy an N-body snapshot";
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
-
void nemo_main(void)
{
stream instr, outstr;
diff --git a/src/nbody/trans/snapdens.c b/src/nbody/trans/snapdens.c
index ca6f70db5..6498425dd 100644
--- a/src/nbody/trans/snapdens.c
+++ b/src/nbody/trans/snapdens.c
@@ -39,14 +39,12 @@ string defv[] = {
"tfactor=-1.0\n conversion factor v->r [virial=sqrt(2)]",
"nn=f\n add NN index to the Key field?",
"ndim=3\n 3dim or 2dim densities?",
- "VERSION=1.5c\n 5-apr-2006 PJT",
+ "VERSION=1.5d\n 20-oct-2024 PJT",
NULL,
};
string usage="density estimator using Kth-nearest neighbor";
-string cvsid="$Id$";
-
#define FAC1 4.188790203 /* 4.pi/3 */
#define FAC2 15.74960994 /* (2.pi)^(3/2) */
@@ -73,7 +71,7 @@ local real raddif(Body *, Body *);
local void stat_nn(Body *);
-nemo_main()
+void nemo_main()
{
stream instr, outstr;
real tsnap, dm;
diff --git a/src/nbody/trans/snapdisk.c b/src/nbody/trans/snapdisk.c
index 6974b7f60..157370f51 100644
--- a/src/nbody/trans/snapdisk.c
+++ b/src/nbody/trans/snapdisk.c
@@ -24,20 +24,19 @@ string defv[] = {
"sign=1\n Sign of rotcur",
"outflow=f\n Outflow or Rotation/Spin",
"times=all\n times of snapshots to copy",
- "VERSION=1.0\n 3-aug-06 PJT",
+ "VERSION=1.1\n 22-dec-2024 PJT",
NULL,
};
string usage="assign rotation to a disk";
-string cvsid="$Id$";
-
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
void get_rotcur(string);
real rotcur(real);
-nemo_main()
+#define MAXLINES 1000000
+
+void nemo_main()
{
stream instr, outstr;
real tsnap, omega, omega2, sign, r,v,f,v2, cost,sint;
@@ -126,7 +125,7 @@ void get_rotcur(string name)
int i,colnr[2], nmax;
real *coldat[2];
- nmax = nemo_file_lines(name);
+ nmax = nemo_file_lines(name, MAXLINES);
if (nmax<1) error("No data to read from %s",name);
rad = (real *) allocate(nmax*sizeof(real));
vel = (real *) allocate(nmax*sizeof(real));
diff --git a/src/nbody/trans/snapgalview.c b/src/nbody/trans/snapgalview.c
index 208d961bd..d7993ca78 100644
--- a/src/nbody/trans/snapgalview.c
+++ b/src/nbody/trans/snapgalview.c
@@ -22,13 +22,12 @@ string defv[] = {
"vel=\n Velocity vector of LSR",
"center=\n Star # to represent LSR (0=first)",
"times=all\n Times to select snapshots from",
- "VERSION=1.0\n 16-feb-03 PJT",
+ "VERSION=1.1\n 22-dec-2024 PJT",
NULL,
};
string usage="coordinate transformations for galactic viewing";
-#define TIMEFUZZ 0.000001 /* tolerance in time comparisons */
void nemo_main()
{
diff --git a/src/nbody/trans/snapinert.c b/src/nbody/trans/snapinert.c
index f667f9985..11675a9be 100644
--- a/src/nbody/trans/snapinert.c
+++ b/src/nbody/trans/snapinert.c
@@ -36,14 +36,12 @@ string defv[] = {
"weight=m\n factor to use in computing center/inertia ",
"per_weight=t\n flag to give inertia per weight basis ",
"tab=f\n flag to produce one-line output",
- "VERSION=2.1a\n 9-may-2023 PJT",
+ "VERSION=2.2\n 22-dec-2024 PJT",
NULL,
};
string usage="get inertia tensor & its eigenvectors, eigenvalues";
-#define TIMEFUZZ 0.001 /* slop tolerated in time comparisons */
-
void snapinert(Body *, int, real, rproc_body, float i[3][3]);
// nr.h
diff --git a/src/nbody/trans/snapkey.c b/src/nbody/trans/snapkey.c
index 9b2a24397..933084400 100644
--- a/src/nbody/trans/snapkey.c
+++ b/src/nbody/trans/snapkey.c
@@ -26,16 +26,13 @@ string defv[] = {
"keys=\n Manually supply all keys",
"keyfile=\n ascii table with keys",
"key=\n bodytrans expression",
- "VERSION=0.3\n 17-feb-2020 PJT",
+ "VERSION=0.4\n 22-dec-2024 PJT",
NULL,
};
string usage="(re)assign keys to a snapshot";
-string cvsid="$Id$";
-
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
-#define MAXNBODY 100000
+#define MAXNBODY 1000000
extern btiproc btitrans(string);
diff --git a/src/nbody/trans/snapmask.c b/src/nbody/trans/snapmask.c
index 3f99f4f6c..ef7850f9f 100644
--- a/src/nbody/trans/snapmask.c
+++ b/src/nbody/trans/snapmask.c
@@ -36,7 +36,7 @@ string defv[] = {
"times=all\n times select string ",
"keyfile=\n file with Key field to select from ",
"keyoffset=0\n offsets to be applied extra to outkey ",
- "VERSION=1.9a\n 9-feb-2024 PJT",
+ "VERSION=1.10\n 22-dec-2024 PJT",
NULL,
};
@@ -55,7 +55,6 @@ local bool *Qsel=NULL; /* pointer to select */
local bool *Qself=NULL; /* pointer to select from file */
local int *sel=NULL;
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
void nemo_main()
diff --git a/src/nbody/trans/snapmass.c b/src/nbody/trans/snapmass.c
index 92d2fccf8..4e11d15a9 100644
--- a/src/nbody/trans/snapmass.c
+++ b/src/nbody/trans/snapmass.c
@@ -15,6 +15,7 @@
* 29-aug-06 c prototypes for frandom() and getrfunc() properly used
* 12-jul-2019 2.2 added ccd=
* 14-aug-2020 2.3 allow Aux to be used
+ * 5-nov-2023 2.4 add timers
*
*/
#include
@@ -46,15 +47,12 @@ string defv[] = {
"ccd=\n Input CCD with mapvalues that represent mass",
"norm=\n Normalization value for the total mass (if used)",
"aux=f\n Store in Aux instead?",
- "VERSION=2.3\n 14-aug-2020 PJT",
+ "VERSION=2.5\n 22-dec-2024 PJT",
NULL,
};
string usage="(re)assign masses to a snapshot";
-string cvsid="$Id$";
-
-#define TIMEFUZZ 0.0001 /* tolerance in time comparisons */
extern real_proc getrfunc(string , string , string , int *);
@@ -72,6 +70,10 @@ void nemo_main(void)
bool Qnorm, first = TRUE;
bool Qaux = getbparam("aux");
string what = (Qaux ? strdup("aux") : strdup("mass"));
+ int itimers=0, ntimers = 100;
+
+ init_timers(ntimers+1);
+ stamp_timers(itimers++); // TIMERS
instr = stropen(getparam("in"), "r");
if (hasvalue("mass")) {
@@ -110,7 +112,6 @@ void nemo_main(void)
Qnorm = hasvalue("norm");
if (Qnorm) norm = getdparam("norm");
-
for(;;) {
/* input main data */
get_history(instr); /* skip history & comments */
@@ -119,6 +120,7 @@ void nemo_main(void)
error("Snapmass (in): Need a snapshot");
}
get_snap(instr, &btab, &nbody, &tsnap, &bits);
+ stamp_timers(itimers++); // TIMERS
if (bits&MassBit) {
dprintf(0,"Warning: existing masses overwritten\n");
if (!Qnorm) {
@@ -132,7 +134,7 @@ void nemo_main(void)
tsnap = 0.0;
}
- if (inmassstr) { /* input mass data from a file */
+ if (inmassstr) { /* input mass data from a different file */
get_history(inmassstr);
if (!get_tag_ok(inmassstr, SnapShotTag))
error("Snapmass (inmass): Need a snapshot");
@@ -148,6 +150,8 @@ void nemo_main(void)
error("essential mass data missing\tbits = %x\n", bits);
}
}
+
+ stamp_timers(itimers++); // TIMERS
/* output */
if (first) {
put_history(outstr);
@@ -212,10 +216,17 @@ void nemo_main(void)
bits |= MassBit; /* turn mass bit on anyways */
if (Qaux)
bits |= AuxBit;
+ stamp_timers(itimers++); // TIMERS
put_snap(outstr, &btab, &nbody, &tsnap, &bits);
}
-
+ stamp_timers(itimers++); // TIMERS
strclose(instr);
if (inmassstr) strclose(inmassstr);
strclose(outstr);
+ stamp_timers(itimers++); // TIMERS
+ for (i=0; i $dir/python_start.csh
+echo "export PATH=${dir}/bin:"'$PATH' > $dir/python_start.sh
+
+if [ $quick = 1 ]; then
+ echo Done.
+fi
+
+# ensure we have a bash kernel in jupyter notebooks
+pip install bash_kernel
+python -m bash_kernel.install
+
+pip freeze > $dir/freeze.log
+echo "Created python_start files for $version in $dir ; no modifications were made to your HOME startup files!"
+python --version
+echo "typically you would need to "
+echo " source $dir/python_start.sh"
+echo "to add this python to your shell environment."
+
+# conda install -y ipython numpy scipy matplotlib notebook jupyter astropy
+
+# pip install astroquery
+# conda -c astropy astroquery
+
+# pip install image_registration radio_beam reproject spectral_cube
diff --git a/src/scripts/mknemo.d/astromatic b/src/scripts/mknemo.d/astromatic
index d8b6c8857..dcfdce89b 100755
--- a/src/scripts/mknemo.d/astromatic
+++ b/src/scripts/mknemo.d/astromatic
@@ -6,7 +6,10 @@
# Packages covered here: skymaker swarp
# Todo: sextractor needs ATLAS or --enable-openblas
# --disable-model-fitting bypasses them
-# ubuntu pacakes: fftw, libopenblas-dev
+# ubuntu packages: fftw, libopenblas-dev
+#
+# failing with the method here: stiff, stuff, eye, psfex
+#
#
version=git
url=https://github.com/astromatic
diff --git a/src/scripts/mknemo.d/boost b/src/scripts/mknemo.d/boost
new file mode 100755
index 000000000..e1ce1680e
--- /dev/null
+++ b/src/scripts/mknemo.d/boost
@@ -0,0 +1,43 @@
+#! /usr/bin/env bash
+#
+# How to Install the Boost Library in C++ on Ubuntu or any other Linux Distribution
+# https://trendoceans.com/install-boost-library-for-c/
+
+# sudo apt install libboost-all-dev
+
+# but if you MUST do it from source, here's the recipe
+# on a fast machine this takes about 2 mins
+
+v=1_81_0 # dec-2022
+v=1_86_0 # aug-2024
+w=https://boostorg.jfrog.io/artifactory/main/release/%s/source/boost_%s.tar.gz
+
+
+#https://github.com/boostorg/boost/tags
+
+cd $NEMO/local
+
+
+vv=$(echo $v | sed s/_/./g)
+url=$(printf $w $vv $v)
+echo url=$url
+tgz=$(basename $url)
+
+if [ ! -e $tgz ]; then
+ wget $url
+else
+ echo "Reusing existing $tgz"
+fi
+
+
+
+tar xf $tgz
+cd boost_$v
+
+# --with-toolset
+./bootstrap.sh --prefix=$NEMO/opt/
+
+./b2
+./b2 install
+
+echo boost $v `date` >> $NEMO/opt/mknemo.log
diff --git a/src/scripts/mknemo.d/falcON2 b/src/scripts/mknemo.d/falcON2
index 73c1dcc09..541c8bbbd 100755
--- a/src/scripts/mknemo.d/falcON2
+++ b/src/scripts/mknemo.d/falcON2
@@ -1,18 +1,55 @@
-#! /bin/csh -f
+#! /bin/bash
#
-# private falcON2
+# this is still the private falcON2
+#
+
+url=git@bitbucket.org:WalterDehnen/falcon.git
+url=git@bitbucket.org:astroteuben/falcon2.git
+url=git@bitbucket.org:peter_teuben/falcon.git # during development
+dir=falcon2
+branch=master
+branch=autoconf # during development
+cgal=0
-set url=git@bitbucket.org:WalterDehnen/falcon.git
-set url=git@bitbucket.org:astroteuben/falcon2.git
-set dir=falcon2
+
+echo "falcON2: using url=$url and git branch=$branch"
cd $NEMO/local
-if (! -e $dir) then
- git clone $url
+if [ ! -e $dir ]; then
+ git clone $url $dir
cd $dir
+ git remote add upstream git@bitbucket.org:WalterDehnen/falcon.git
else
cd $dir
git pull
-endif
+fi
+git checkout $branch
+
+# only for linux for now
+if [ ! -e tbb ]; then
+ echo "TBB not present, grabbing the linux one"
+ wget https://github.com/uxlfoundation/oneTBB/releases/download/v2020.3/tbb-2020.3-lin.tgz
+ tar zxf tbb-2020.3-lin.tgz
+fi
+
+./configure
+source falcon_start.sh
+
+# if a local version is needed
+if [ $cgal = 1 ]; then
+ v=5.6.2
+ wget https://github.com/CGAL/cgal/releases/download/v${v}/CGAL-${v}.tar.xz
+ tar xf CGAL-${v}.tar.xz
+ cd CGAL-${v}
+ mkdir build
+ cd build
+ cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$FALCON ..
+ make install
+fi
+
-echo "Installation is still manual. See INSTALL.md"
+make -C utils2 inc/cachesize.h
+make -C h5pp
+make inc/utils inc/kepler inc/h5pp
+make clean
+make -j 4
diff --git a/src/scripts/mknemo.d/gcc b/src/scripts/mknemo.d/gcc
new file mode 100755
index 000000000..03c2d6b2f
--- /dev/null
+++ b/src/scripts/mknemo.d/gcc
@@ -0,0 +1,60 @@
+#! /usr/bin/env bash
+#
+# sudo apt install build-essential
+# sudo apt install libmpfr-dev libgmp3-dev libmpc-dev -y
+
+# --enable-checking=release (esp. useful for compile-time tests
+# as it disables some consistency checks in the compile),
+# --disable-bootstrap (builds faster, but disables some consistency checks).
+#
+# ~2h dart76 just l=c and j=1
+
+V=10.5.0 # July 2023
+v=11.5.0
+v=12.4.0
+v=13.3.0
+v=14.2.0 # version (mid 2024)
+a=x86_64-linux-gnu # architecture
+b=86 # lib name
+j=4 # number of processors for makea
+l=c,c++,fortran # languages to compile
+
+for arg in $*; do
+ export $arg
+done
+
+cd $NEMO/local
+
+root=$NEMO/opt/gcc-${v}
+lib=lib$b
+tgz=gcc-${v}.tar.gz
+if [ ! -e $tgz ]; then
+ wget http://ftp.gnu.org/gnu/gcc/gcc-${v}/$tgz
+fi
+
+tar -xf $tgz
+cd gcc-${v}
+./configure -v \
+ --build=${a} \
+ --host=${a} \
+ --target=${a} \
+ --prefix=$NEMO/opt/gcc-${v} \
+ --enable-checking=release \
+ --enable-languages=$l \
+ --disable-multilib \
+ --program-suffix=-${v}
+make clean
+make -j $j
+make install
+
+
+
+rc=$NEMO/opt/gcc-${v}/gcc_start.sh
+
+echo "_r=$root" > $rc
+echo 'export PATH=$_r/bin:$PATH' >> $rc
+echo 'export LD_LIBRARY_PATH=$_r/$lib:$LD_LIBRARY_PATH' >> $rc
+
+echo "Wrote $rc"
+
+echo gcc '"'$v'"' `date` >> $NEMO/opt/mknemo.log
diff --git a/src/scripts/mknemo.d/hdf5 b/src/scripts/mknemo.d/hdf5
index 6d0526fee..2cdf2928a 100755
--- a/src/scripts/mknemo.d/hdf5
+++ b/src/scripts/mknemo.d/hdf5
@@ -7,12 +7,10 @@
# 17-mar-2020 wgetc
#
-# if ($?NEMO == 0) set NEMO=`pwd`
-
-
-version=1.10.7 # 11-sep-2020
+# highest versions in respective 1.xx versions
+version=1.10.11 # oct-2021
version=1.12.1 # 06-jul-2021
-version=1.13.0 # 01-dec-2021
+version=1.13.0 # 01-dec-2021 falcon2 cannot yet use 1.13 and up
version=1.13.3 # 27-oct-2022
version=1.14.3 # 30-oct-2023
url=https://github.com/HDFGroup/hdf5
@@ -31,9 +29,12 @@ $wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-$mversion/hdf5-$versio
tar zxf hdf5-$version.tar.gz
cd hdf5-$version
-./configure --prefix=$NEMO/opt --enable-build-mode=production
+./configure --prefix=$NEMO/opt --enable-build-mode=production --enable-cxx
make -j
make install
echo hdf5 $version `date` >> $NEMO/opt/mknemo.log
echo Example to test linking: mknemo
+
+# future 2.0 release will use cmake, e.g.
+# cmake .. -DCMAKE_INSTALL_PREFIX=${NEMO}/opt
diff --git a/src/scripts/nemo.bench b/src/scripts/nemo.bench
index c6521ee61..ea08ff6f7 100755
--- a/src/scripts/nemo.bench
+++ b/src/scripts/nemo.bench
@@ -21,6 +21,7 @@
ften=1
f5=1
np=1
+ sleep=0
for arg in $*; do\
export $arg
@@ -29,7 +30,7 @@ done
mkdir $tmp
cd $tmp
-echo "NEMOBENCH: (2022-09-14) $tmp : nbody0=$nbody0 nbody1=$nbody1 nbody3=$nbody3 size=$size ften=$ften clean=$clean bsf=$bsf tmp=$tmp"
+echo "NEMOBENCH: (2024-10-18) $tmp : nbody0=$nbody0 nbody1=$nbody1 nbody3=$nbody3 size=$size ften=$ften clean=$clean bsf=$bsf tmp=$tmp"
echo `hostname`
echo `uname -a`
echo `date`
@@ -90,21 +91,26 @@ if [ $mode == 5 ]; then
echo "Tstop: $t1 $t2 $t3 $t4 $t5 $t6 [5sec CPU for each on a 11th Gen Intel(R) Core(TM) i5-1135G7 up to 4.2GHz using gcc-11]"
echo -n .;
echo -n 1;$time directcode nbody=$nbody0 out=. tstop=$t1 seed=123 $help > bench2.directcode.log 2>&1
+ sleep $sleep
echo -n .;$time mkplummer p0 $nbody2 seed=123 $help > bench0.mkplummer.log 2>&1
echo -n 2;$time gyrfalcON p0 . kmax=6 tstop=$t2 eps=0.05 $help > bench2.gyrfalcon.log 2>&1
+ sleep $sleep
echo -n 3;$time hackcode1 nbody=$nbody1 out=. seed=123 tstop=$t3 $help > bench2.hackcode1.log 2>&1
+ sleep $sleep
echo -n .;$time mkorbit o0 x=1 e=1 lz=1 potname=log $help > bench0.mkorbit.log 2>&1
echo -n 4;$time orbint o0 . tstop=$t4 dt=0.01 nsave=\$tstop ndiag=\$tstop/100 $help > bench2.orbint.log 2>&1
+ sleep $sleep
echo -n 5;$time potcode p0 . freqout=10 freq=1000 tstop=$t5 potname=plummer $help > bench2.potcode.log 2>&1
+ sleep $sleep
echo -n 6;$time treecode1 nbody=$nbody1 out=. seed=123 tstop=$t6 $help > bench2.treecode1.log 2>&1
echo -n .;
echo ""
cat bench2.*.log | perl -pe 's/(CPU_USAGE)/\n$1/' | grep CPU_USAGE > cpu_usage.log
column -t cpu_usage.log
- cpu_mean=$(awk '{print $4}' cpu_usage.log | tabstat - | grep mean: | awk '{print $2}')
- score=$(nemoinp 5000/$cpu_mean)
- ntest=$(cat cpu_usage.log | wc -l)
+ _out=($(tabstat cpu_usage.log 4 bad=0 | txtpar - 5000/%1,%2 p0=mean,1,2 p1=npt,1,2))
+ score=${_out[0]}
+ ntest=${_out[1]}
echo "NEMOBENCH5 score: $score"
echo "$score $ntest $(hostname) $(date)" >> $NEMO/nemobench5.log
@@ -113,7 +119,7 @@ if [ $mode == 5 ]; then
cd ..
rm -rf $tmp
else
- echo
+ echo "More detailed results in $tmp"
fi
exit 0
diff --git a/src/scripts/python/fitsplot.py b/src/scripts/python/fitsplot.py
index c2e45b72c..d471e4e66 100755
--- a/src/scripts/python/fitsplot.py
+++ b/src/scripts/python/fitsplot.py
@@ -4,11 +4,16 @@
import os
import sys
-import aplpy
import argparse
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
+try:
+ import aplpy
+except:
+ print("sorry, aplpy not working")
+ sys.exit(1)
+
help_main = ["Simple color plot of a FITS image",
diff --git a/src/tutor/mp/Makefile b/src/tutor/mp/Makefile
index 1c6864096..14ce94f3d 100644
--- a/src/tutor/mp/Makefile
+++ b/src/tutor/mp/Makefile
@@ -81,3 +81,4 @@ bench9: scaling2
@echo Long integration, see laptop performance drop, see core swaps
$(TIME) ./scaling2 iter=500 > bench9.log 2>&1
grep cputime bench9.log | sed s/###// | tabplot - 6 8 line=1,1 ycoord=0 ymin=0
+ @echo "Jansky 0.077 k2 0.090 0.098"