1 #! /bin/csh -f 2 # 3 # Demo1 : create an exponential disk, embed it in some halo and integrate 4 # then extract a given time, freeze potential, and check orbits 5 # also project and make velocity field, and determine rotation curve 6 # and compare model with model-model 7 # 8 9 # parameters to be changed by the user 10 set run=run1 # identification, and basename for all files 11 set nbody=1000 # number of bodies 12 set Qtoomre=1.2 # 'Q' for the disk 13 set potname=plummer ## not used yet (for external halo) 14 set potpars=0,0 ## not used yet 15 set xymovie=0 # show a movie? 16 set orbit= # orbit to be selected in xymovie 17 set pa=30 # PA of disk on sky (E through N) 18 set inc=60 # INC of disk 19 set phi=30 # angle of bar with major axis of galaxy (on sky) 20 set beam=0.1 # FWHM of beam 21 set range=2 # gridding from -range:range 22 set nsize=64 # number of pixels in gridding (px=py=2*range/nx) 23 set radii=5:120:5 # set some rings for rotcur 24 set yanc=0 # use YANC or NEMO's default integrator (hackcode1) 25 26 set goto=start 27 28 # some other parameters, probably not useful for the user at an initial stage 29 set freq=32 # integration frequency (i.e. 1/timestep) 30 set freqout=2 # output frequency (i.e. 1/dumptimes) 31 set tstop=10 # total integration time 32 set tview=10 # viewing time for velocity field analysis 33 set eps=0.05 # softening length 34 set tol=0.7 # opening angle for treecode 35 set tstart=0 # starting time for analysis of pattern speed 36 set seed=123 # random seed for sprinkling particles 37 set rcut=2 # cutoff radius of exponential disk 38 39 # parse the command line arguments (they will then override the above defaults) 40 41 foreach a ($*) 42 set $a 43 end 44 45 46 # derive some parameters that appear common or logically belong together 47 48 set pot_pars=(potname=$potname potpars=$potpars) 49 set tree_pars=(eps=$eps tol=$tol) 50 set grid_pars=(xrange=-${range}:${range} yrange=-${range}:${range} nx=$nsize ny=$nsize) 51 set cell=`nemoinp "2*$range/$nsize*60"` 52 set step=`nemoinp 1/$freqout` 53 54 # keep a log, incase we call this routine multiple times 55 echo `date` :: $* >> $run.history 56 57 goto $goto 58 59 # ================================================================================ START 60 start: 61 62 rm -f $run.* >& /dev/null 63 64 65 echo Creating exponential disk with $nbody particles 66 mkexpdisk out=$run.01 nbody=$nbody Qtoomre=$Qtoomre seed=$seed rcut=$rcut tab=t > $run.tab 67 echo "Integrating to time=$tstop ..." 68 if ($yanc) then 69 echo "Using YANC with $tree_pars" 70 set y=$run.01y 71 echo "## Created for YANC by NEMO demo1 script" > $y 72 echo "# theta $tol" >> $y 73 echo "# eps $eps" >> $y 74 echo "# time 0.0" >> $y 75 echo "# kernel 30" >> $y 76 echo "# h0 5" >> $y 77 echo "# N $nbody" >> $y 78 echo "# data_format 0" >> $y 79 echo "# start_data" >> $y 80 snapprint $run.01 m,x,y,z,vx,vy,vz >> $y 81 echo YANC $y $y.out $tstop $step 0 82 time YANC $y $y.out $tstop $step 0 > $run.2.log 83 set t=0.0 84 tabcomment $y - delete=t |\ 85 tabtos - - block1=mass,x,y,z,vx,vy,vz nbody=$nbody times=$t > $run.02 86 rm -f $y 87 foreach file ($y.out.*) 88 tabcomment $file - delete=t |\ 89 tabtos - - block1=mass,x,y,z,vx,vy,vz,phi nbody=$nbody times=$t |\ 90 csf - - item=SnapShot >> $run.02 91 rm -f $file 92 set t=`nemoinp $t+$step` 93 end 94 else 95 echo "Using hackcode1 with $tree_pars" 96 time hackcode1 in=$run.01 out=$run.02 \ 97 options=mass,phi,acc \ 98 $tree_pars \ 99 freq=$freq freqout=$freqout tstop=$tstop > $run.2.log 100 endif 101 102 # 103 # =========================================================================== MOVIE 104 # 105 movie: 106 107 if ($xymovie) then 108 (snapxyz $run.02 - | xyzview - orbit=$orbit maxpoint="2*$nbody") & 109 endif 110 111 # 112 # =========================================================================== TRIM 113 # 114 trim: 115 rm -f $run.1? $run.1?? >& /dev/null 116 echo "Trimming and centering at given time $tview" 117 118 #snaptrim $run.2 $run.3 times=$tview 119 snaptrim $run.02 - times=$tview |\ 120 snapcenter - $run.12 weight="-phi*phi*phi" 121 echo :::SNAPRECT 122 snaprect $run.12 $run.13 weight="-phi*phi*phi" 123 echo :::SNAPINERT 124 snapinert $run.12 - weight="-phi*phi*phi" 125 snapinert $run.12 - weight="-phi*phi*phi" tab=t 126 echo :::SNAPKINEM 127 snapkinem $run.12 128 129 echo "Could study pattern speed by running do2" 130 echo "Followed by align.awk and tablsqfit" 131 132 snapgrid $run.13 $run.14 $grid_pars 133 ccdprint $run.14 x= y= newline=t label=x,y > $run.15 134 tabtos $run.15 $run.16 block1=x,y 135 snapmass $run.16 $run.17 mass=1 136 137 echo "Compute potentials on a regular grid" 138 139 hackforce $run.13 $run.19 $run.17 $tree_pars fcells=2 140 snapgrid $run.19 $run.110 $grid_pars evar=phi mean=t 141 ccdflip $run.110 - x | ccdflip - $run.111 y 142 ccdmath $run.110,$run.111 $run.112 "(%1+%2)/2" 143 ccdmath $run.110,$run.111 $run.113 "(%1-%2)/(%1+%2)" 144 145 echo "Study sample orbits by running 'do1'" 146 # do1 [run=run1] [iorbit=0] [range=2] [omega=0] 147 148 # 149 # =========================================================================== ROTCUR 150 # 151 rotcur: 152 153 echo "Creating a velocity field" 154 rm -f $run.2? $run.2?? $run.rotcur >& /dev/null 155 snaprotate $run.13 $run.20 "atand(tand($phi)/cosd($inc)),$inc,$pa" zyz 156 snapgrid $run.20 $run.21 $grid_pars moment=0 157 snapgrid $run.20 $run.22 $grid_pars moment=1 158 snapgrid $run.20 $run.23 $grid_pars moment=2 159 ccdsmooth $run.21 $run.21s $beam 160 ccdsmooth $run.22 $run.22s $beam 161 ccdsmooth $run.23 $run.23s $beam 162 ccdmath $run.21s,$run.22s,$run.23s $run.20d %1 163 ccdmath $run.21s,$run.22s,$run.23s $run.20v "%2/%1*100" 164 ccdmath $run.21s,$run.22s,$run.23s $run.20s "100*sqrt(%3/%1-%2*%2/(%1*%1))" 165 166 rotcur $run.20v units=arcmin \ 167 tab=$run.rotcur \ 168 radii=$radii pa=$pa inc=$inc vrot=100 fixed=vsys,pa,inc,xpos,ypos > $run.rotcur.log 169 tabcomment $run.rotcur - > $run.rotcur1 170 171 tabcomment $run.rotcur - | tabmath - - %1,%4,%5 all > $run.rottab 172 173 174 175 # care about units ... we're doing arcmin vs. arcsec here 176 echo Creating model from rotcur, with cell=$cell 177 ccdvel $run.24 rotcurfit=run1.rotcur1 size=$nsize cell=$cell 178 echo And subtracting it, look at $run.24d for the difference map 179 ccdmath $run.20v,$run.24 $run.24d fie="ifeq(%1,0.0,0.0,ifeq(%2,0.0,0.0,%1-%2))" 180 181 # 182 # =========================================================================== OMEGA 183 # 184 omega: 185 echo Determining the pattern speed over times ${tstart}:${tstop}:$step 186 ## get the pattern speed from an analysis of the direction of the moment of inertia tensor 187 do2 run=$run times=${tstart}:${tstop}:$step > $run.ome 188 # DEMO PLOT omega.gif/ps 189 align.awk $run.ome | tabplot - 1 2,3 2 10 -1080 360 line=1,1 color=2,3 point=2,0.1 \ 190 xlab=Time ylab=Phi headline="$*" 191 align.awk $run.ome | tablsqfit - 1 2 192 align.awk $run.ome | tablsqfit - 1 3 193 194 exit 0 195 196 ###### demo commands for plots: ####################################################### 197 # 198 # DEMO PLOT evolution.gif/ps 199 # overview of evolution 200 snapplot run1.02 nxy=5,5 nxticks=3 nyticks=3 yapp=/vps 201 202 # get angle and shape for a set of times 203 do2 times=0:10:0.5 > run1.ome 204 tabmath run1.ome - %1,%4,%5,%6,%5/%4,%6/%4 all |\ 205 tabplot - 1 5,6 0 10 0 1 line=1,1 color=2,3 point=2,0.2 yapp=/vps 206 207 # DEMO PLOT: shape.gif/ps 208 align.awk run1.ome |\ 209 tabplot - 1 2,3 0 10 line=1,1 color=2,3 point=2,0.2 yapp=/vps 210 align.awk run1.ome | tablsqfit - 1 2 211 align.awk run1.ome | tablsqfit - 1 3 212 # get 'a =', and multiply by pi/180 to get the pattern speed omega 213 214 ## orbit analysis via 'do1' - which is tkrun enabled... 215 216 tkrun do1