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