#! /bin/csh -f
#
#  Demo2 : create an exponential disk, let a bar form and compute how fast
#          that bar is spinning. Should make a polynomial fit how much the
#          slowdown is... Based on demo1 for the GH2001 workshop
#  Needs:    align.awk
#  NEMOBIN:  mkexpdisk snapscale snaprotate snapadd snapprint hackcode1 YancNemo
#            tabcomment tablsqfit tabmath tabplot tabhist tabtos nemoinp
#            snaprect snapinert snapsort 
#
#  Updates:
#    10-jun-2001   added yappext=, store keywords in a $run.keys
#    11-jun-2001   bit more documentation
#    13-jun-2001   even more documentation for shell :-)
#

#================================================================================
#  parameters to be changed by the user  (if you add a new parameter, also add
#  it to keywords() below
set run=run0               # identification, and basename for all files
set file=""                # alternate way of running script where simulation had been done
                           # e.g. file=XXX would set run=XXX also
set nbody=1000             # number of bodies
set Qtoomre=1.2            # initial 'Q' for the disk 
set xymovie=0              # show a movie? (needs vogl based xyzview)
set yanc=2                 # use YANC or NEMO's default integrator (hackcode1)
set gv=0                   # view resulting plots on screen?

#   alternate entry point in script (use with care)
#   movie and omega are really the only other entry points left here ....
set goto=start

#  some other parameters, probably not useful for the user at an initial stage
set freq=32                # integration frequency (i.e. 1/timestep)
set hmin=6                 # integration step for YANC
set freqout=2              # output frequency (i.e. 1/dumptimes)
set tstop=10               # total integration time
set tview=10               # viewing time for velocity field analysis
set eps=0.05               # softening length
set kernel=20              # softening kernel (only used in YANC)
set tol=0.7                # opening angle for treecode
set tstart=0               # starting time for analysis of pattern speed
set seed=123               # random seed for sprinkling particles
set rcut=2                 # cutoff radius of exponential disk
set hackdens=0             # override potentials with hackdens based densities
set pfract=0.5             # fraction of bound particles to use for shape
set symm=0                 # only seed 1/2 particles, set other half from symmetry
set log=RESULTS.log        # logfile of overall results
set yappext=/vcps          # any optional directons for the YAPP device (PGPLOT, ..)

#================================================================================
#  parse the command line arguments (they will then override the above defaults)
#  don't add any user modifyable keywords below this line, and if you add new
#  ones, also add it to keyword() below

foreach a ($*)
  set $a
end

#  each keyword from the list above should be in this list !!!
set keywords=(run file nbody Qtoomre xymovie yanc gv goto freq hmin freqout tstop tview eps kernel tol tstart seed rcut hackdens pfract symm log yappext)

#  before we do anything, if "file"  is used, "run" is overriden, and will not be used
#  and also the simulation is now assumed to exist in "file"
if (X$file != X) then
  set run=$file
  if (! -e $run) then
    echo File $file does not exist
    exit 1
  endif
  if ($goto == start) set goto=omega
endif

#  Save all the current settings in $run.keys
set keycsh=$run.keycsh
rm -rf $keycsh
foreach k ($keywords)
  echo echo KEYWORD: $k='$'$k >> $keycsh
end
(source $keycsh) > $run.keys

#  derive some parameters that appear common or logically belong together

set tree_pars=(eps=$eps tol=$tol)
set step=`nemoinp 1/$freqout`
set nfract=`nemoinp "$nbody*$pfract" format=%d`

#  keep a log, in case we call this routine multiple times
if (-e $run.history) then
  echo `date` :: $* > $run.history
else
  echo `date` :: $* >> $run.history
endif

goto $goto

# ================================================================================ START
start:

rm -f $run.* >& /dev/null


echo Creating exponential disk with $nbody particles
if ($symm == 1) then
  echo Creating a mirror image and then adding it
  mkexpdisk out=$run.01 nbody=$nbody/2 Qtoomre=$Qtoomre seed=$seed rcut=$rcut \
	    tab=t time=0 headline="$*" > $run.tab
  snapscale $run.01  $run.01a mscale=0.5
  snapscale $run.01a $run.01b rscale=-1 vscale=-1 
  rm $run.01
  snapadd $run.01a,$run.01b $run.01
else if ($symm == 2) then
  echo Creating 3 images by rotation and then adding them
  mkexpdisk out=$run.01 nbody=$nbody/4 Qtoomre=$Qtoomre seed=$seed rcut=$rcut \
	    tab=t time=0 headline="$*" > $run.tab
  snapscale  $run.01  $run.01a mscale=0.25
  snaprotate $run.01a $run.01b  90 z
  snaprotate $run.01a $run.01c 180 z
  snaprotate $run.01a $run.01d 270 z
  rm $run.01
  snapadd $run.01a,$run.01b,$run.01c,$run.01d $run.01
else if ($symm == 3) then
  echo Creating 2 mirror images and rotating two others, then adding them
  mkexpdisk out=$run.01 nbody=$nbody/4 Qtoomre=$Qtoomre seed=$seed rcut=$rcut \
	    tab=t time=0 headline="$*" > $run.tab
  snapscale  $run.01  $run.01a mscale=0.25
  snaprotate $run.01a $run.01b  90 z
  snapscale  $run.01a $run.01c rscale=-1 vscale=-1
  snaprotate $run.01c $run.01d  90 z
  rm $run.01
  snapadd $run.01a,$run.01b,$run.01c,$run.01d $run.01
else
  mkexpdisk out=$run.01 nbody=$nbody Qtoomre=$Qtoomre seed=$seed rcut=$rcut tab=t \
	    headline="$*" time=0 > $run.tab
endif


echo "Integrating to time=$tstop ..."
if ($yanc == 1) then
  echo "Using YANC with $tree_pars kernel=$kernel hmin=$hmin"
  set y=$run.01y
  echo "## Created for YANC by NEMO $0 script"  > $y
  echo "#  theta $tol"    >> $y
  echo "#  eps $eps"      >> $y
  echo "#  time 0.0"      >> $y
  echo "#  kernel $kernel">> $y
  echo "#  hmin $hmin"    >> $y
  echo "#  Nbref 8"        >> $y
  echo "#  N  $nbody"     >> $y
  echo "#  give_pot 1"    >> $y
  echo "#  give_rho 1"    >> $y
  echo "#  data_format 0" >> $y
  echo "#  start_data"    >> $y
  snapprint $run.01 m,x,y,z,vx,vy,vz >> $y
  echo YANC $y $y.out $tstop $step 0
  time YANC $y $y.out $tstop $step 0 > $run.2.log
  set t=0.0
  if (0) then                          # we fixed YANC locally (27may)
    tabcomment $y - delete=t |\
	tabtos - - block1=mass,x,y,z,vx,vy,vz nbody=$nbody times=$t > $run.02
  endif
  rm -f $y
  foreach file ($y.out.*)
    tabcomment $file - delete=t |\
	tabtos - - block1=mass,x,y,z,vx,vy,vz,phi,aux nbody=$nbody times=$t |\
	csf - - item=SnapShot >> $run.02
    rm -f $file
    set t=`nemoinp $t+$step`
  end
else if ($yanc == 2) then
  echo "Using YancNewmo with $tree_pars kernel=$kernel hmin=$hmin"
  time YancNemo in=$run.01 out=$run.02 \
    eps=$eps theta=$tol kernel=$kernel \
    tstop=$tstop step=$step hmin=$hmin give_pot=1 give_rho=1 > $run.2.log
else
  echo "Using hackcode1 with $tree_pars"
  time hackcode1 in=$run.01 out=$run.02 \
    options=mass,phi,acc \
    $tree_pars \
    freq=$freq freqout=$freqout tstop=$tstop > $run.2.log
endif

#  
# =========================================================================== MOVIE
# show a movie, you need libvogl to compile xyzview
movie:

if ($xymovie) then
  (snapxyz $run.02 - | xyzview - orbit=$orbit maxpoint="2*$nbody") &
endif

#  
# =========================================================================== OMEGA
#
omega:
echo Determining the pattern speed over times ${tstart}:${tstop}:$step 
echo Using $nfract / $nbody of the particles with lowest potential

set times=${tstart}:${tstop}:$step 

#    we allow a few different weighting schemes, -phi**3 seems to be the best
if ($hackdens == 1) then
  set weight=("phi")
else if ($hackdens == 2) then
  set weight=("phi*phi")
else
  set weight=("-phi*phi*phi")
endif

#                      loop over all times requested
rm -f $run.ome
foreach t (`nemoinp $times`)
 rm -f $run.snaptmp
 if ($hackdens) then
    snaptrim $run.02 - times=$t debug=-1 |\
        hackdens - $run.snaptmp 
    snaprect $run.snaptmp . weight="$weight" > $run.tmp
 else
    #  only if hackdens=0 will we sort the particles by potential
    #  and cut out the deepest fraction
    snaptrim $run.02 $run.snaptmp times=$t debug=-1 
    snaprect $run.snaptmp . weight="$weight" > $run.tmp
    snapsort $run.snaptmp - phi |\
	snapmask - - 0:$nfract |\
	snaprect - . weight="$weight" > $run.tmpf
 endif
 # get the angles of the bar in X and Y for all particles
 set ex=(`grep e_x $run.tmp | awk -F: '{print $2}'`) 
 set ey=(`grep e_y $run.tmp | awk -F: '{print $2}'`) 
 if ($#ex != 6) continue
 # and also for the fraction of particles most bound
 set exf=(`grep e_x $run.tmpf | awk -F: '{print $2}'`) 
 set eyf=(`grep e_y $run.tmpf | awk -F: '{print $2}'`) 
 if ($#exf != 6) continue
 # and run this to get the axis ratios of the moments of inertia
 # BUG?! we're still using all particles here, not the 'pfract' fraction!!!
 snapinert $run.snaptmp - weight="$weight" tab=t > $run.tmp 
 set si=(`cat $run.tmp`)

 #   time PHI   PHIf     Ixx     Iyy     Izz
 echo $t $ex[6] $exf[6] $si[11] $si[12] $si[13] >> $run.ome
end

#                      clean up the mess
rm -f $run.tmp $run.snaptmp


#  align the angles in columns 2 and 3 (clockwise bars only)
#  and plot PHI(t)
align.awk $run.ome | tabplot - 1 2,3 0 10 -1080 360 line=1,1 color=2,3 point=2,0.1 \
    xlab=Time ylab=Phi headline="$*"  yapp=$run.plot1.ps$yappext
if ($gv) gv $run.plot1.ps &

#   fit polynomials of degrees 1,2,3 and 4 to PHI(t)
#   often 1 is not  good enough, 2 is mostly ok.
foreach o (1 2 3 4)
  rm -f $run.ome_$o
  align.awk $run.ome |\
    tablsqfit - 1 2 fit=poly order=$o out=$run.ome_$o tab=t  > $run.fit.ome_$o
end
#   paste the different fits and residuals back together for a plot
paste $run.ome_1 $run.ome_2 $run.ome_3 $run.ome_4 |\
  tabmath - - %1,%3-%2,%6-%5,%9-%8,%12-%11 all > $run.ome_0
tabplot $run.ome_0 1 2,3,4,5 0 10 -40 40 line=1,1 color=2,3,4,5 point=2,0.1 ycoord=0 \
    headline="$*" yapp=$run.plot2.ps$yappext
if ($gv) gv $run.plot2.ps &

#   compute Iyy/Ixx and Izz/Ixx and plot the shapes as function of time
tabmath $run.ome - %1,%4,%5,%6,%5/%4,%6/%4 all |\
   tabplot - 1 5,6 0 10 0 1 line=1,1 color=2,3 point=2,0.2 \
	headline="$*" yapp=$run.plot3.ps$yappext
set shape=(`tabmath $run.ome - %5/%4 all`)
echo $shape >> shape.$run.tab

if ($gv) gv $run.plot3.ps &
# plot of energy conservation, different programs store this differently
# so currenlky only do this for YANC.
if ($yanc) then
    tabplot $run.2.log 1 3 line=1,1 ycoord=0.5 yapp=$run.plot4.ps$yappext
    tabhist $run.2.log 3   yapp=$run.plot5.ps$yappext >& $run.e
    tabplot $run.2.log 1 2 line=1,1 yapp=$run.plot6.ps$yappext
else
    echo no non-yanc plots: need to figure out how to make hackcode1 plots
    echo or replace with new treecode....
endif

#  Make summary table (some kind of RESULTS.log file)
#  this table is appended only, so useful if you run the script
#  many times and remove the data.

if ($yanc) then
  echo -n "$yanc $eps $kernel $tol $seed $nbody $symm" >> $log
  foreach o (1 2 3 4)
    echo -n " `cat $run.fit.ome_$o | tail -1` " >> $log
  end
  echo -n `grep Mean $run.e|awk -F: '{print $2}'` >> $log
  echo "" >> $log
endif